R version 3.5.1 (2018-07-02) -- "Feather Spray" Copyright (C) 2018 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # Load required libraries > library(lavaan) This is lavaan 0.6-3 lavaan is BETA software! Please report any bugs. > > ########################################################################################################## > > > # Read in data > bacs = read.csv(file.choose(), row.names=1) > bacs <- na.omit(bacs) # remove NAs: two samples don't have soil data > > head(bacs) Site Invasion Season pH P C NH4 NO3 N FV_11_In_AUT FlowerValley Invaded Autumn 4.4 8 2.05 0.097 20.3 0.06 FV_12_In_AUT FlowerValley Invaded Autumn 4.4 8 1.79 0.094 15.1 0.05 FV_8_In_AUT FlowerValley Invaded Autumn 5.1 17 1.33 0.168 18.8 0.08 FV_9_In_AUT FlowerValley Invaded Autumn 4.8 7 2.26 0.104 35.1 0.06 FV_10_In_SPR FlowerValley Invaded Spring 5.5 7 3.90 0.124 1.0 0.11 FV_7_In_SPR FlowerValley Invaded Spring 5.0 11 1.60 0.045 2.6 0.06 NMDS1 NMDS2 S H FV_11_In_AUT -0.113286880 -0.125007956 3027 1009.1662 FV_12_In_AUT -0.108981745 -0.006953161 3369 1132.3253 FV_8_In_AUT -0.053937509 -0.053159631 3150 927.5607 FV_9_In_AUT -0.056755292 -0.019473018 3539 1285.1441 FV_10_In_SPR -0.011440393 -0.066204960 3379 1291.5111 FV_7_In_SPR -0.000830224 -0.095768133 2889 863.1870 > str(bacs) # so P is picked up as an integer --> fix! 'data.frame': 78 obs. of 13 variables: $ Site : Factor w/ 5 levels "FlowerValley",..: 1 1 1 1 1 1 1 1 1 1 ... $ Invasion: Factor w/ 2 levels "Invaded","Uninvaded": 1 1 1 1 1 1 1 1 2 2 ... $ Season : Factor w/ 2 levels "Autumn","Spring": 1 1 1 1 2 2 2 2 1 1 ... $ pH : num 4.4 4.4 5.1 4.8 5.5 5 4.5 4.7 4.5 5.6 ... $ P : int 8 8 17 7 7 11 9 5 8 9 ... $ C : num 2.05 1.79 1.33 2.26 3.9 1.6 3.94 2.03 2.4 3.06 ... $ NH4 : num 0.097 0.094 0.168 0.104 0.124 0.045 0.158 0.066 0.126 0.171 ... $ NO3 : num 20.3 15.1 18.8 35.1 1 2.6 35.1 0.3 5.2 31.8 ... $ N : num 0.06 0.05 0.08 0.06 0.11 0.06 0.1 0.07 0.07 0.1 ... $ NMDS1 : num -0.1133 -0.109 -0.0539 -0.0568 -0.0114 ... $ NMDS2 : num -0.12501 -0.00695 -0.05316 -0.01947 -0.0662 ... $ S : int 3027 3369 3150 3539 3379 2889 3275 3278 3453 2933 ... $ H : num 1009 1132 928 1285 1292 ... - attr(*, "na.action")= 'omit' Named int 56 77 ..- attr(*, "names")= chr "VM_5_In_SPR" "WA_25_Un_SPR" > > # Correct value types > bacs$P <- as.numeric(bacs$P) > str(bacs) 'data.frame': 78 obs. of 13 variables: $ Site : Factor w/ 5 levels "FlowerValley",..: 1 1 1 1 1 1 1 1 1 1 ... $ Invasion: Factor w/ 2 levels "Invaded","Uninvaded": 1 1 1 1 1 1 1 1 2 2 ... $ Season : Factor w/ 2 levels "Autumn","Spring": 1 1 1 1 2 2 2 2 1 1 ... $ pH : num 4.4 4.4 5.1 4.8 5.5 5 4.5 4.7 4.5 5.6 ... $ P : num 8 8 17 7 7 11 9 5 8 9 ... $ C : num 2.05 1.79 1.33 2.26 3.9 1.6 3.94 2.03 2.4 3.06 ... $ NH4 : num 0.097 0.094 0.168 0.104 0.124 0.045 0.158 0.066 0.126 0.171 ... $ NO3 : num 20.3 15.1 18.8 35.1 1 2.6 35.1 0.3 5.2 31.8 ... $ N : num 0.06 0.05 0.08 0.06 0.11 0.06 0.1 0.07 0.07 0.1 ... $ NMDS1 : num -0.1133 -0.109 -0.0539 -0.0568 -0.0114 ... $ NMDS2 : num -0.12501 -0.00695 -0.05316 -0.01947 -0.0662 ... $ S : int 3027 3369 3150 3539 3379 2889 3275 3278 3453 2933 ... $ H : num 1009 1132 928 1285 1292 ... - attr(*, "na.action")= 'omit' Named int 56 77 ..- attr(*, "names")= chr "VM_5_In_SPR" "WA_25_Un_SPR" > > ########################################################################################################## > > > > > ########################################################################################################## > > scaledBacs <- scale(bacs[,4:13]) #this is because the variances of some variables much higher than others, and gives error when fitting > scaledBacs <- cbind(bacs[,1:3], scaledBacs) # add the other columns back > > > ID <- bacs$Season=="Spring" # to subset either of the two sampling seasons > > > # MODEL 1: > > model1 <- ' + # latent variables + Soil =~ pH + P + NO3 + N + C + NH4 + + # regressions + NMDS1 ~ Soil + Invasion + Soil ~ Invasion + ' > > fit <- sem(model1, data=scaledBacs[ID,], cluster="Site") Warning message: In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING: The variance-covariance matrix of the estimated parameters (vcov) does not appear to be positive definite! The smallest eigenvalue (= -5.902741e-14) is smaller than zero. This may be a symptom that the model is not identified. > summary(fit, standardized=TRUE) lavaan 0.6-3 ended normally after 66 iterations Optimization method NLMINB Number of free parameters 23 Number of observations 38 Number of clusters [Site] 5 Estimator ML Robust Model Fit Test Statistic 201.106 NA Degrees of freedom 19 19 P-value (Chi-square) 0.000 NA Scaling correction factor NA for the Yuan-Bentler correction (Mplus variant) Parameter Estimates: Information Observed Observed information based on Hessian Standard Errors Robust.cluster Latent Variables: Estimate Std.Err z-value P(>|z|) Std.lv Std.all Soil =~ pH 1.000 0.648 0.675 P 0.598 0.075 8.005 0.000 0.388 0.439 NO3 0.072 0.037 1.971 0.049 0.047 0.295 N 1.479 0.393 3.762 0.000 0.959 0.990 C 1.068 0.288 3.711 0.000 0.692 0.882 NH4 1.315 0.371 3.549 0.000 0.853 0.982 Regressions: Estimate Std.Err z-value P(>|z|) Std.lv Std.all NMDS1 ~ Soil 0.949 0.071 13.298 0.000 0.616 0.642 Invasion 0.102 0.115 0.887 0.375 0.102 0.053 Soil ~ Invasion -0.218 0.142 -1.535 0.125 -0.337 -0.168 Intercepts: Estimate Std.Err z-value P(>|z|) Std.lv Std.all .pH 0.315 0.491 0.641 0.522 0.315 0.327 .P 0.158 0.427 0.370 0.711 0.158 0.179 .NO3 -0.303 0.025 -11.927 0.000 -0.303 -1.916 .N 0.684 0.347 1.974 0.048 0.684 0.706 .C 0.334 0.235 1.420 0.156 0.334 0.426 .NH4 0.443 0.306 1.446 0.148 0.443 0.510 .NMDS1 0.211 0.586 0.360 0.719 0.211 0.220 .Soil 0.000 0.000 0.000 Variances: Estimate Std.Err z-value P(>|z|) Std.lv Std.all .pH 0.504 0.246 2.048 0.041 0.504 0.545 .P 0.628 0.391 1.608 0.108 0.628 0.807 .NO3 0.023 0.020 1.128 0.259 0.023 0.913 .N 0.018 0.013 1.457 0.145 0.018 0.019 .C 0.137 0.098 1.396 0.163 0.137 0.222 .NH4 0.027 0.029 0.949 0.343 0.027 0.036 .NMDS1 0.548 0.274 1.997 0.046 0.548 0.596 .Soil 0.409 0.385 1.062 0.288 0.972 0.972 > round(fitMeasures(fit)[c("chisq", "pvalue", "rmsea")],4) chisq pvalue rmsea 201.1061 0.0000 0.5022 Warning message: In pchisq(XX2, df = df2, ncp = ncp) : NaNs produced > > > > > # MODEL 2: > > model2 <- ' + # regressions + NMDS1 ~ pH + NH4 + NO3 + N + P + C + Invasion + pH ~ Invasion + P ~ Invasion + NH4 ~ Invasion + NO3 ~ Invasion + N ~ Invasion + C ~ Invasion + ' > > fit <- sem(model2, data=scaledBacs[ID,], cluster="Site") Warning message: In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING: The variance-covariance matrix of the estimated parameters (vcov) does not appear to be positive definite! The smallest eigenvalue (= -2.549433e-14) is smaller than zero. This may be a symptom that the model is not identified. > summary(fit, standardized=TRUE) lavaan 0.6-3 ended normally after 54 iterations Optimization method NLMINB Number of free parameters 27 Number of observations 38 Number of clusters [Site] 5 Estimator ML Robust Model Fit Test Statistic 267.422 14231.309 Degrees of freedom 15 15 P-value (Chi-square) 0.000 0.000 Scaling correction factor 0.019 for the Yuan-Bentler correction (Mplus variant) Parameter Estimates: Information Observed Observed information based on Hessian Standard Errors Robust.cluster Regressions: Estimate Std.Err z-value P(>|z|) Std.lv Std.all NMDS1 ~ pH 1.087 0.048 22.676 0.000 1.087 0.971 NH4 -0.223 0.083 -2.678 0.007 -0.223 -0.180 NO3 0.401 0.077 5.220 0.000 0.401 0.059 N 0.010 0.055 0.180 0.857 0.010 0.009 P -0.036 0.039 -0.939 0.348 -0.036 -0.030 C 0.106 0.048 2.217 0.027 0.106 0.077 Invasion -0.052 0.100 -0.521 0.602 -0.052 -0.024 pH ~ Invasion -0.087 0.275 -0.317 0.751 -0.087 -0.045 P ~ Invasion -0.130 0.307 -0.423 0.673 -0.130 -0.073 NH4 ~ Invasion -0.341 0.268 -1.273 0.203 -0.341 -0.196 NO3 ~ Invasion -0.084 0.047 -1.779 0.075 -0.084 -0.267 N ~ Invasion -0.312 0.280 -1.114 0.265 -0.312 -0.161 C ~ Invasion -0.022 0.312 -0.069 0.945 -0.022 -0.014 Intercepts: Estimate Std.Err z-value P(>|z|) Std.lv Std.all .NMDS1 0.277 0.123 2.248 0.025 0.277 0.257 .pH 0.118 0.703 0.168 0.867 0.118 0.123 .P 0.157 0.742 0.211 0.833 0.157 0.177 .NH4 0.524 0.249 2.107 0.035 0.524 0.603 .NO3 -0.200 0.089 -2.251 0.024 -0.200 -1.265 .N 0.667 0.344 1.942 0.052 0.667 0.689 .C 0.017 0.328 0.051 0.959 0.017 0.021 Variances: Estimate Std.Err z-value P(>|z|) Std.lv Std.all .NMDS1 0.018 0.003 6.325 0.000 0.018 0.016 .pH 0.922 0.289 3.191 0.001 0.922 0.998 .P 0.774 0.276 2.803 0.005 0.774 0.995 .NH4 0.725 0.494 1.469 0.142 0.725 0.961 .NO3 0.023 0.019 1.223 0.221 0.023 0.929 .N 0.914 0.596 1.533 0.125 0.914 0.974 .C 0.616 0.398 1.549 0.121 0.616 1.000 > round(fitMeasures(fit)[c("chisq", "pvalue", "rmsea")],4) chisq pvalue rmsea 267.4222 0.0000 0.6655