Common Path Model in One Twin Sample

###############
## This script merges runs a common path model in one twin sample


library(psych)
library(OpenMx)
mxOption(NULL, "Default optimizer", "SLSQP")



#start with a data frame in wide format where one row is a twin pair


########################
### Biometric Models ###
########################
###
### Select variables
Vars <- c("CDres", "ADres", "ADHDres", "MJres")
nv <- length(Vars)
ntv <- nv*2
selVars <- paste(Vars, c(rep(0,nv), rep(1,nv)), sep="")

### Select Data for Analysis
mtfsMzData    <- subset(mtfs.w, zygosity==1, selVars)
mtfsDzData    <- subset(mtfs.w, zygosity==2, selVars)


# Generate Descriptive Statistics
colMeans(mtfsMzData,na.rm=TRUE)
colMeans(mtfsDzData,na.rm=TRUE)
round(cor(mtfsMzData,use="pairwise.complete.obs"),2)
round(cor(mtfsDzData,use="pairwise.complete.obs"),2)

mtfsDataMZ    <- mxData( observed=mtfsMzData, type="raw" )
mtfsDataDZ    <- mxData( observed=mtfsDzData, type="raw" )




### ------------------------------------------------------------------------------
### Common Pathway ACE model
### ------------------------------------------------------------------------------
nl <- 1 # number of latent factors
# Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s)
Xm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels="xm", lbound=.00001, name="Xm" )
Ym <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels="ym", lbound=.00001, name="Ym" )
Zm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels="zm", lbound=.00001, name="Zm" )


Alm <- mxAlgebra(Xm %*% t(Xm), name="Alm")
Clm <- mxAlgebra(Ym %*% t(Ym), name="Clm")
Elm <- mxAlgebra(Zm %*% t(Zm), name="Elm")

# Matrices as, cs, and es to store a, c, and e path coefficients for specific factors
pathAsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("asm1", "asm2", "asm3", "asm4"), name="asm" )
pathCsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("csm1", "csm2", "csm3", "csm4"), name="csm" )
pathEsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("esm1", "esm2", "esm3", "esm4"), name="esm" )

# Matrix and Algebra for constraint on variance of latent phenotype
covarLPm <- mxAlgebra( expression= Alm + Clm + Elm, name="CovarLPm" )
varLPm <- mxAlgebra( expression= diag2vec(CovarLPm), name="VarLPm" )
unitm <- mxMatrix( type="Unit", nrow=nl, ncol=1, name="Unitm")
varLP1m <- mxConstraint( expression=VarLPm == Unitm, name="varLP1m")


### Matrix f for factor loadings on latent phenotype
pathFlm <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=TRUE, values=.2, lbound=.000001, labels=c("f1m", "f2m", "f3m", "f4m"), name="flm")

# Matrices A, C, and E compute variance components
covAm <- mxAlgebra( expression=flm %*% Alm %*% t(flm) + asm %*% t(asm), name="Am" )
covCm <- mxAlgebra( expression=flm %*% Clm %*% t(flm) + csm %*% t(csm), name="Cm" )
covEm <- mxAlgebra( expression=flm %*% Elm %*% t(flm) + esm %*% t(esm), name="Em" )


covPm <- mxAlgebra( expression=Am+Cm+Em, name="Vm" )
matIm <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="Im")
invSDm <- mxAlgebra( expression=solve(sqrt(Im*Vm)), name="iSDm")
meanGm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=1, labels=c("m1m", "m2m", "m3m", "m4m","m1m", "m2m", "m3m", "m4m"), name="expMeanm" )


covMZm <- mxAlgebra( expression= rbind( cbind(Am+Cm+Em, Am+Cm),
                                        cbind(Am+Cm,  Am+Cm+Em)), name="expCovMZm" )
covDZm <- mxAlgebra( expression= rbind( cbind(Am+Cm+Em, 0.5%x%Am+Cm),
                                        cbind(0.5%x%Am+Cm, Am+Cm+Em)), name="expCovDZm" )

   
### Combine Groups
objMZm <- mxExpectationNormal( covariance="expCovMZm", means="expMeanm", dimnames=selVars )
objDZm <- mxExpectationNormal( covariance="expCovDZm", means="expMeanm", dimnames=selVars )


mCI <- mxCI("f4m")

parsm <- list( Xm, Ym, Zm, Alm, Clm, Elm, covarLPm, varLPm, unitm, pathFlm, pathAsm, pathCsm, pathEsm, covAm, covCm, covEm, covPm, matIm, invSDm, meanGm, mCI )

funML <- mxFitFunctionML()

mtfsModelMZ <- mxModel( parsm, covMZm, mtfsDataMZ, objMZm, funML, name="MZm" )
mtfsModelDZ <- mxModel( parsm, covDZm, mtfsDataDZ, objDZm, funML, name="DZm" )



fitML <- mxFitFunctionMultigroup(c("MZm.fitfunction","DZm.fitfunction"))
Cf <- mxModel( "ComACE", parsm, varLP1m, mtfsModelMZ, mtfsModelDZ,  fitML, mxCI(c('Alm', 'Clm', 'Elm', 'Alc','Clc','Elc')))

CfFit <- mxRun(Cf, intervals=TRUE) ### Run ComACE model
summary(CfFit, verbose=T)

RefModCom <- mxRefModels(CfFit, run=TRUE)
summary(CfFit, refModels=RefModCom)
summary(RefModCom$Saturated)
### Get standardized loadings and specific factors
round(CfFit$iSDm$result %*% CfFit$flm$values, 2)