# R script of Structural equation modelling as used in the manuscript "Prokaryotic community composition and extracellular polymeric substances affect soil microaggregation in carbonate containing semiarid grasslands"
# Details on data is provided in file "Data_1st_field_paper"
# Last modified by J.H.T. Zethof, 07/01/2020
# Load necessary package
install.packages("lavaan")
library(lavaan)
# Set working direction (where data file is stored)
Workdir <- "G:/TUDresden/Project/Processing_R"
# Load data
data1<- read.table(paste(Workdir,"/Dataset_SEM.txt", sep = ""), header=TRUE,sep="",dec=",")
# Normalize data with calculating Z-score, per site
data_Z_score_RH <-data.frame(data1[data1$Site==1,1:6],(scale(data1[data1$Site==1,7:(ncol(data1))],center=TRUE,scale=TRUE)))
data_Z_score_AB <-data.frame(data1[data1$Site==2,1:6],(scale(data1[data1$Site==2,7:(ncol(data1))],center=TRUE,scale=TRUE)))
# Create Model for Structural Equation Modelling, by exchanging MWD_0Jml for MWD_200Jml the SEM of the supplementary figure 13 can be obtained
MyModel <- '
PC1_site ~ a*OC + b*TN + c*IC
X16S_rRNA ~ d*OC + e*TN + f*IC
PC1_site ~ g*X16S_rRNA
OC ~ h*TN
Clay ~ i*MWD_text
EPS_Sacc ~ j*PC1_site + k*OC + l*TN + m*X16S_rRNA
Con_ang ~ o*EPS_Sacc + p*OC
MWD_0Jml ~ q*MWD_text + r*Clay + s*EPS_Sacc + t*OC + u*IC + v*Con_ang
'
# Structural equation modelling, note that by changing data input, sites can be switched
fit.MM <- lavaan::sem(MyModel, estimator = "MLR", missing = "ML", data = data_Z_score_AB, fixed.x=FALSE, cluster="Subplot") #error is neglectable as it approximates 0 very well
summary(fit.MM,fit.measures=TRUE, standardized = TRUE)