##### CODE FOR ILLUSTRATIONS IN PAPER
# version 2

### load libraries
library(lavaan)

########## Simulation Study #1 ##########
# Covariance matrix for input
sim.cov <- '
 1.2472 
-0.1492  1.019 
 0.8442 -0.178  1.6417  
 0.0358  0.704  0.0357  1.5512  
-0.2922 -0.233 -0.1217 -0.5276   1.488 
 0.5235  0.149  0.5335  0.3277  -0.655  1.029'

sim.cov.dat <- getCov(sim.cov, names = c("BioA1", "BioB1", "BioA2", "BioB2", "StemChg12", "Stems1"))


##### Scenario 1 - Set N=150 #####
### Model 1 - initial model
mod1 <- '
# regressions
BioA1 ~ Stems1
BioB1 ~ Stems1
BioA2 ~ BioA1 +Stems1 +BioB1
BioB2 ~ BioB1 +Stems1 +BioA1
StemChg12 ~ Stems1 +BioA2 +BioB2'

mod1.fit <- sem(mod1, sample.cov = sim.cov.dat, sample.nobs = 150)
# Global Fit Measures
summary(mod1.fit, fit.measures=T)
# Residuals
resid(mod1.fit, type="cor")
# Modification Indices
subset(modindices(mod1.fit), mi > 2, c("lhs", "op", "rhs", "mi", "epc"))
# For examination of individual parameter support
parameterEstimates(mod1.fit)


### Model 2 (parameter numbers correspond to those in Fig.2)
mod2 <- '
# regressions
BioA1 ~ b2*Stems1
BioB1 ~ b7*Stems1
BioA2 ~ b3*BioA1 +b5*Stems1 +b11*BioB1
BioB2 ~ b8*BioB1 +b10*Stems1 +b6*BioA1
StemChg12 ~ b1*Stems1 +b4*BioA2 +b9*BioB2
# error covariance
BioA1 ~~ b12*BioB1'

mod2.fit <- sem(mod2, sample.cov = sim.cov.dat, sample.nobs = 150)
# Global Fit Measures
summary(mod2.fit, fit.measures=T)
# Residuals
resid(mod2.fit, type="cor")
# Modification Indices
subset(modindices(mod2.fit), mi > 2, c("lhs", "op", "rhs", "mi", "epc"))
# For examination of individual parameter support
parameterEstimates(mod2.fit)


### Model 3 (parameter b4 now estimates an error correlation)
mod3 <- '
# regressions
BioA1 ~ b2*Stems1
BioB1 ~ b7*Stems1
BioA2 ~ b3*BioA1 +b5*Stems1 +b11*BioB1
BioB2 ~ b8*BioB1 +b10*Stems1 +b6*BioA1
StemChg12 ~ b1*Stems1 +b9*BioB2
# error covariance
BioA1 ~~ b12*BioB1
StemChg12 ~~ b4*BioA2'

mod3.fit <- sem(mod3, sample.cov = sim.cov.dat, sample.nobs = 150)
# Global Fit Measures
summary(mod3.fit, fit.measures=T)
# Residuals
resid(mod3.fit, type="cor")
# Modification Indices
subset(modindices(mod3.fit), mi > 2, c("lhs", "op", "rhs", "mi", "epc"))
# For examination of individual parameter support
parameterEstimates(mod3.fit)


### Model 3B (parameter b6 now set to zero)
mod3B <- '
# regressions
BioA1 ~ b2*Stems1
BioB1 ~ b7*Stems1
BioA2 ~ b3*BioA1 +b5*Stems1 +b11*BioB1
BioB2 ~ b8*BioB1 +b10*Stems1 +b6*BioA1
StemChg12 ~ b1*Stems1 +b9*BioB2
# error covariance
BioA1 ~~ b12*BioB1
StemChg12 ~~ b4*BioA2
# set parameters to zero
b6 == 0'

mod3B.fit <- sem(mod3B, sample.cov = sim.cov.dat, sample.nobs = 150)
# Global Fit Measures
summary(mod3B.fit, fit.measures=T)
# Residuals
resid(mod3B.fit, type="cor")
# Modification Indices
subset(modindices(mod3B.fit), mi > 2, c("lhs", "op", "rhs", "mi", "epc"))
# For examination of individual parameter support
parameterEstimates(mod3B.fit)


### Model 3C (parameter b11 now set to zero)
mod3C <- '
# regressions
BioA1 ~ b2*Stems1
BioB1 ~ b7*Stems1
BioA2 ~ b3*BioA1 +b5*Stems1 +b11*BioB1
BioB2 ~ b8*BioB1 +b10*Stems1 +b6*BioA1
StemChg12 ~ b1*Stems1 +b9*BioB2
# error covariance
BioA1 ~~ b12*BioB1
StemChg12 ~~ b4*BioA2
# set parameters to zero
b6 == 0
b11 == 0'

mod3C.fit <- sem(mod3C, sample.cov = sim.cov.dat, sample.nobs = 150)
# Global Fit Measures
summary(mod3C.fit, fit.measures=T)
# Residuals
resid(mod3C.fit, type="cor")
# Modification Indices
subset(modindices(mod3C.fit), mi > 2, c("lhs", "op", "rhs", "mi", "epc"))
# For examination of individual parameter support
parameterEstimates(mod3C.fit)


### Model 3D (parameter b4 now set to zero)
mod3D <- '
# regressions
BioA1 ~ b2*Stems1
BioB1 ~ b7*Stems1
BioA2 ~ b3*BioA1 +b5*Stems1 +b11*BioB1
BioB2 ~ b8*BioB1 +b10*Stems1 +b6*BioA1
StemChg12 ~ b1*Stems1 +b9*BioB2
# error covariance
BioA1 ~~ b12*BioB1
StemChg12 ~~ b4*BioA2
# set parameters to zero
b6 == 0
b11 == 0
b4 == 0'

mod3D.fit <- sem(mod3D, sample.cov = sim.cov.dat, sample.nobs = 150)
# Global Fit Measures
summary(mod3D.fit, fit.measures=T)
# Residuals
resid(mod3D.fit, type="cor")
# Modification Indices
subset(modindices(mod3D.fit), mi > 2, c("lhs", "op", "rhs", "mi", "epc"))
# For examination of individual parameter support
parameterEstimates(mod3D.fit)


##### Multimodel Comparisons using AICc #####
library (AICcmodavg)

aictab(list(mod3.fit, mod3B.fit, mod3C.fit, mod3D.fit), 
       c("MOD3", "MOD3B","MOD3C","MOD3D"))


### Full Results for Model 3B
summary(mod3B.fit, fit.measures=T, rsq=T, standardized=T)
parameterEstimates(mod3B.fit, standardized=T)
standardizedsolution(mod3B.fit)


##### Additional Illustration: Setting b7 to zero for Model 3B #####
### Model 3B.alt (parameter b7 now set to zero)
mod3Balt <- '
# regressions
BioA1 ~ b2*Stems1
BioB1 ~ b7*Stems1
BioA2 ~ b3*BioA1 +b5*Stems1 +b11*BioB1
BioB2 ~ b8*BioB1 +b10*Stems1 +b6*BioA1
StemChg12 ~ b1*Stems1 +b9*BioB2
# error covariance
BioA1 ~~ b12*BioB1
StemChg12 ~~ b4*BioA2
# set parameters to zero
b6 == 0
b7 == 0'

mod3Balt.fit <- sem(mod3Balt, sample.cov = sim.cov.dat, sample.nobs = 150)
summary(mod3Balt.fit, fit.measures=T, rsq=T, standardized=T)