# 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 on 14/04/2020, with R version 3.6.2 (2019-12-12)

# Load necessary package
if(!require(lavaan)){install.packages("lavaan")}
library(lavaan)

# Set working direction (where data file is stored)
Workdir <- "C:/Users/Gebruiker/Downloads"
# 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
# For Rambla Honda (Figure 6a):
MyModel_RH <- ' 
    PC1_site ~ a*OC + c*IC + w*Graphitic_C 
    X16S_rRNA ~ d*OC + f*IC 
    PC1_site ~ g*X16S_rRNA 
    Graphitic_C ~ x*IC 
    Clay ~ i*MWD_text
    EPS_Sacc  ~  j*PC1_site + k*OC + 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 for the site Alboloduy, as presented in Figure 6b
fit.MM.RH <- lavaan::sem(MyModel_RH, estimator = "MLR", missing = "ML", data = data_Z_score_RH, fixed.x=FALSE, cluster="Subplot") 
#error is neglectable as is a roundoff error resulting from the fact that there are less datapoints than parameters in the model.
summary(fit.MM.RH,fit.measures=TRUE, standardized = TRUE)
# the proportion of variances 
inspect(fit.MM.RH,"r2")

# For Alboloduy (Figure 6b):
MyModel_AB <- ' 
    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 for the site Alboloduy, as presented in Figure 6b
fit.MM.AB <- lavaan::sem(MyModel_AB, estimator = "MLR", missing = "ML", data = data_Z_score_AB, fixed.x=FALSE, cluster="Subplot") 
#error is neglectable as is a roundoff error resulting from the fact that there are less datapoints than parameters in the model.
summary(fit.MM.AB,fit.measures=TRUE, standardized = TRUE)
# the proportion of variances 
inspect(fit.MM.AB,"r2")