#MetaSEM twin modelling

#Main lavaan model and basic set up borrowed from...

#-------------------------------------Twin modelling script-------------------------------------#
#-------------------Chapter 4.6, Beaujean 2014 Latent Variable Modeling Using R-----------------#


# Environment set up
rm(list = ls()) # clear R
options(max.print=1000000) # set max print

# Load lavaan
library(lavaan)
library(metaSEM)

# MZ twins
MZ <- lav_matrix_lower2full(c(1, 0.7735, 1)) 
rownames(MZ) <- colnames(MZ) <- c("P1", "P2") 
# DZ twins
DZ <- lav_matrix_lower2full(c(1, 0.3871, 1))
rownames(DZ) <- colnames(DZ) <- c("P1", "P2")

# MZ twins asymp cov
MZacov <- 0.00625681


# DZ twins asymp cov
DZacov <- 0.00311364


# combine the covariances, sample sizes, asympt covs
twin.cov <- list(MZ = MZ, DZ = DZ) 
twin.acov <- list(MZ = MZacov, DZ = DZacov)
twin.n <- list(MZ =  8308, DZ = 11983) 


ace.model<-"
A1=~ NA*P1 + c(a,a)*P1 
A2=~ NA*P2 + c(a,a)*P2 
C1 =~ NA*P1 + c(c,c)*P1
C2 =~ NA*P2 + c(c,c)*P2
# variances
A1 ~~ 1*A1
A2 ~~ 1*A2
C1 ~~ 1*C1
C2 ~~ 1*C2 
P1~~c(e2,e2)*P1 
P2~~c(e2,e2)*P2
# covariances
A1 ~~ c(1,.5)*A2 
A1 ~~ 0*C1 + 0*C2 
A2 ~~ 0*C1 + 0*C2 
C1 ~~ c(1,1)*C2"

#standard ace model ignoring heterogeneity
ace.fit<-cfa(ace.model, sample.cov=twin.cov, sample.nobs=twin.n)
summary(ace.fit, standardized = TRUE, ci = TRUE)
standardizedsolution(ace.fit, ci = TRUE)

#metaSEM model
RAM1 <- lavaan2RAM(ace.model, obs.variables=c('P1', 'P2'), ngroups=2)

wls(Cov = twin.cov, aCov = twin.acov, 
    n = twin.n, RAM = RAM1,  cor.analysis = TRUE)