# 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)