##### Latent Trajectory Modeling in lavaan

library(lavaan)
library(AICcmodavg)

setwd("E:\ppt_files\_education\SEM.11.7-LatentTrajectoryModels")
dat <- read.csv("SEM_11_7_Latent_Trajectory_Modeling_in_Lavaan_data.csv")
names(dat)

dat2 <- with(dat, data.frame(age))
dat2$r1 <- dat$sr94adj
dat2$r2 <- dat$sr95adj
dat2$r3 <- dat$sr96adj
dat2$r4 <- dat$sr97adj
dat2$r5 <- dat$sr98adj
dat2$c1 <- dat$covtot94
dat2$c2 <- dat$covtot95
dat2$c3 <- dat$covtot96
dat2$c4 <- dat$covtot97
dat2$c5 <- dat$covtot98
dat2$abio <- dat$abio
dat2$fire <- dat$fireint1
summary(dat2)


### Model 101: simple latent curve
mod.101 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5'

fit.101 <- growth(mod.101, data=dat2)
print(fit.101)
summary(fit.101)
subset(modindices(fit.101), mi > 3.8)

### diagnostic plot for mod.101
# observed means for richness
r1.mean<-mean(dat2$r1); r2.mean<-mean(dat2$r2); r3.mean<-mean(dat2$r3); 
r4.mean<-mean(dat2$r4); r5.mean<-mean(dat2$r5)
rich.mean <- c(r1.mean, r2.mean, r3.mean, r4.mean, r5.mean)
# model predictions for richness
i.101 = 43.953
s.101 = -3.991
r1.pred = i.101 +s.101*0
r2.pred = i.101 +s.101*1
r3.pred = i.101 +s.101*2
r4.pred = i.101 +s.101*3
r5.pred = i.101 +s.101*4
rich.pred.101 <- c(r1.pred, r2.pred, r3.pred, r4.pred, r5.pred)
plot(rich.mean ~ rich.pred.101, pch=19); abline(lm(rich.mean ~ rich.pred.101))
plot(rich.mean, xlab="time", pch=19)
# several possibilities: (1) time-varying covariate, (2) autoregressive effect
# plot rich.mean.101 ~ cover.mean.101
c1.mean<-mean(dat2$c1); c2.mean<-mean(dat2$c2); c3.mean<-mean(dat2$c3); 
c4.mean<-mean(dat2$c4); c5.mean<-mean(dat2$c5)
cover.mean <- c(c1.mean, c2.mean, c3.mean, c4.mean, c5.mean); 
plot(cover.mean, xlab="time", pch=19)
plot(rich.mean ~ cover.mean, pch=19); abline(lm(rich.mean ~ cover.mean))


### Model 102: Include cover as time-varying covariate 
mod.102 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5
# time-varying regressions
r1 ~ c1
r2 ~ c2
r3 ~ c3
r4 ~ c4
r5 ~ c5'

fit.102 <- growth(mod.102, data=dat2)
print(fit.102)
summary(fit.102)
subset(modindices(fit.102), mi > 3.8)


### Model 103: Include autoregressive effect from r1 to r2
mod.103 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5
# time-varying regressions
r1 ~ c1
r2 ~ c2
r3 ~ c3
r4 ~ c4
r5 ~ c5
# autoregressive effects
r2 ~ r1'

fit.103 <- growth(mod.103, data=dat2)
print(fit.103)
summary(fit.103)
subset(modindices(fit.103), mi > 3.8)
fitMeasures(fit.108, "gfi")


### Model 104: Include lag effect of c1 on r2 - doesn't help
mod.104 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5
# time-varying regressions
r1 ~ c1
r2 ~ c2 +c1
r3 ~ c3
r4 ~ c4
r5 ~ c5
# autoregressive effect
r2 ~ r1'

fit.104 <- growth(mod.104, data=dat2)
print(fit.104)
summary(fit.104)
subset(modindices(fit.104), mi > 3.8)


### Model 105: Allow nonlinear effect between times 2 and 3 - doesn't help
mod.105 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +lambda3*r3 +3*r4 +4*r5
# time-varying regressions
r1 ~ c1
r2 ~ c2
r3 ~ c3
r4 ~ c4
r5 ~ c5
# autoregressive effect
r2 ~ r1'

fit.105 <- growth(mod.105, data=dat2)
print(fit.105)
summary(fit.105)
subset(modindices(fit.105), mi > 3.8)


### Model 106: Add time-invariant covariate to Model 103
mod.106 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5
# time-varying regressions
r1 ~ c1
r2 ~ c2
r3 ~ c3
r4 ~ c4
r5 ~ c5
# autoregressive effects
r2 ~ r1
# time-invariant effects
i ~ abio
s ~ abio'

fit.106 <- growth(mod.106, data=dat2)
print(fit.106)
summary(fit.106)
subset(modindices(fit.106), mi > 3.8)


### Model 107: Add fire severity and age effects
mod.107 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5
# time-varying regressions
r1 ~ c1
r2 ~ c2
r3 ~ c3
r4 ~ c4
r5 ~ c5
# autoregressive effects
r2 ~ r1
# time-invariant effects
i ~ abio
s ~ abio
# fire severity effects
r1 ~ fire +abio
c1 ~ fire +abio
fire ~ age +abio'

fit.107 <- growth(mod.107, data=dat2)
print(fit.107)
summary(fit.107, rsq=TRUE)
subset(modindices(fit.107), mi > 3.8)
fitMeasures(fit.107, "gfi")

### Model 108: Prune Model 107, remove abio effect on slope
mod.108 <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5
# time-varying regressions
r1 ~ c1
r2 ~ c2
r3 ~ c3
r4 ~ c4
r5 ~ c5
# autoregressive effects
r2 ~ r1
# time-invariant effects
i ~ abio
#s ~ abio
# fire severity effects
r1 ~ fire +abio
c1 ~ fire +abio
fire ~ age +abio'

fit.108 <- growth(mod.108, data=dat2)
print(fit.108)
summary(fit.108, rsq=TRUE)
subset(modindices(fit.108), mi > 3.8)
fitMeasures(fit.108, "gfi")

# model predictions for richness from model 108
i.108 = 3.271
s.108 = -1.356
r1.pred = i.108 +s.108*0 +0.129*mean(dat2$c1) -1.883*mean(dat2$fire) +0.227*mean(dat2$abio)
r2.pred = i.108 +s.108*1 +0.028*mean(dat2$c2)
r3.pred = i.108 +s.108*2 +0.084*mean(dat2$c3)
r4.pred = i.108 +s.108*3 +0.018*mean(dat2$c4)
r5.pred = i.108 +s.108*4 +0.017*mean(dat2$c5)
rich.pred.108 <- c(r1.pred, r2.pred, r3.pred, r4.pred, r5.pred)
plot(rich.mean ~ rich.pred.108, pch=19); abline(lm(rich.mean ~ rich.pred.108))


### Model 108alt: Prune Model 108 of cover effects
mod.108alt <- ' 
# intercept and slope with fixed coefficients
i =~ 1*r1 +1*r2 +1*r3 +1*r4 +1*r5
s =~ 0*r1 +1*r2 +2*r3 +3*r4 +4*r5
# autoregressive effects
r3 ~ r2
# time-invariant effects
i ~ abio
#s ~ abio
# fire severity effects
r1 ~ fire +abio
fire ~ age +abio'

fit.108alt <- growth(mod.108alt, data=dat2)
print(fit.108alt)
summary(fit.108alt, rsq=TRUE)
subset(modindices(fit.108alt), mi > 3.8)
fitMeasures(fit.108alt, "gfi")