Linear Growth Model in One Twin Sample
Jump to navigation
Jump to search
#this is an example R script to run a basic linear growth model in one twin sample
#using time-varying definition variables for slope factor loadings
library(psych)
library(OpenMx)
library(plyr)
###########################################################
#begin with a dataframe in wide format
mtfs.w <- reshape(mtfsT, v.names=c("ADHDres", "ADres", "CDres", "mj_12m_freq_14", "age14", "mj_12m_freq_17", "age17", "mj_12m_freq_20", "age20", "mj_12m_freq_24", "age24", "mj_12m_freq_29", "age29"), timevar="twin", idvar="famid", direction="wide", sep="")
#format age and definition variables
#get an average age within each twin pair
mtfs.w$age14avg <- (mtfs.w$age140 + mtfs.w$age141)/2
mtfs.w$age17avg <- (mtfs.w$age170 + mtfs.w$age171)/2
mtfs.w$age20avg <- (mtfs.w$age200 + mtfs.w$age201)/2
mtfs.w$age24avg <- (mtfs.w$age240 + mtfs.w$age241)/2
mtfs.w$age29avg <- (mtfs.w$age290 + mtfs.w$age291)/2
#if one twin in the twin pair is missing on age, use the available twin’s age
mtfs.w$age14twinpair <- ifelse(is.na(mtfs.w$age140), mtfs.w$age141, mtfs.w$age14avg)
mtfs.w$age14twinpair <- ifelse(is.na(mtfs.w$age141), mtfs.w$age140, mtfs.w$age14avg)
summary(mtfs.w$age140)
summary(mtfs.w$age141)
summary(mtfs.w$age14avg)
summary(mtfs.w$age14twinpair)
mtfs.w$age17twinpair <- ifelse(is.na(mtfs.w$age170), mtfs.w$age171, mtfs.w$age17avg)
mtfs.w$age17twinpair <- ifelse(is.na(mtfs.w$age171), mtfs.w$age170, mtfs.w$age17avg)
summary(mtfs.w$age170)
summary(mtfs.w$age171)
summary(mtfs.w$age17avg)
summary(mtfs.w$age17twinpair)
mtfs.w$age20twinpair <- ifelse(is.na(mtfs.w$age200), mtfs.w$age201, mtfs.w$age20avg)
mtfs.w$age20twinpair <- ifelse(is.na(mtfs.w$age201), mtfs.w$age200, mtfs.w$age20avg)
summary(mtfs.w$age200)
summary(mtfs.w$age201)
summary(mtfs.w$age20avg)
summary(mtfs.w$age20twinpair)
mtfs.w$age24twinpair <- ifelse(is.na(mtfs.w$age240), mtfs.w$age241, mtfs.w$age24avg)
mtfs.w$age24twinpair <- ifelse(is.na(mtfs.w$age241), mtfs.w$age240, mtfs.w$age24avg)
summary(mtfs.w$age240)
summary(mtfs.w$age241)
summary(mtfs.w$age24avg)
summary(mtfs.w$age24twinpair)
mtfs.w$age29twinpair <- ifelse(is.na(mtfs.w$age290), mtfs.w$age291, mtfs.w$age29avg)
mtfs.w$age29twinpair <- ifelse(is.na(mtfs.w$age291), mtfs.w$age290, mtfs.w$age29avg)
summary(mtfs.w$age290)
summary(mtfs.w$age291)
summary(mtfs.w$age29avg)
summary(mtfs.w$age29twinpair)
#if both twins are missing age for a wave, use the average age at that wave
mtfs.w$age14twinpair[is.na(mtfs.w$age14twinpair)] <- mean(mtfs.w$age14twinpair, na.rm=T)
mtfs.w$age17twinpair[is.na(mtfs.w$age17twinpair)] <- mean(mtfs.w$age17twinpair, na.rm=T)
mtfs.w$age20twinpair[is.na(mtfs.w$age20twinpair)] <- mean(mtfs.w$age20twinpair, na.rm=T)
mtfs.w$age24twinpair[is.na(mtfs.w$age24twinpair)] <- mean(mtfs.w$age24twinpair, na.rm=T)
mtfs.w$age29twinpair[is.na(mtfs.w$age29twinpair)] <- mean(mtfs.w$age29twinpair, na.rm=T)
#center each age at 14
mtfs.w$age14twinpairC <- mtfs.w$age14twinpair - 14
mtfs.w$age17twinpairC <- mtfs.w$age17twinpair - 14
mtfs.w$age20twinpairC <- mtfs.w$age20twinpair - 14
mtfs.w$age24twinpairC <- mtfs.w$age24twinpair - 14
mtfs.w$age29twinpairC <- mtfs.w$age29twinpair - 14
### Create dummy intercept definition variable and dummy 0 variable
mtfs.w$intercept <- 1
mtfs.w$zero <- 0
########################
### Biometric Models ###
########################
### Select observed variables
VarsMN <- c("mj_12m_freq_14", "mj_12m_freq_17", "mj_12m_freq_20", "mj_12m_freq_24", "mj_12m_freq_29")
nvMN <- length(VarsMN)
ntvMN <- nvMN*2
selVarsMN <- paste(VarsMN, c(rep(0,nvMN), rep(1,nvMN)), sep="")
############
###Select definition variables
DefVarsMN <- c("age14twinpairC", "age17twinpairC", "age20twinpairC", "age24twinpairC", "age29twinpairC", "intercept", "sex")
### Select Data for Analysis
mtfsMzData <- subset(mtfs.w, zygosity==1, select=c(selVarsMN, DefVarsMN))
mtfsDzData <- subset(mtfs.w, zygosity==2, select=c(selVarsMN, DefVarsMN))
mtfsDataMZ <- mxData( observed=mtfsMzData, type="raw" )
mtfsDataDZ <- mxData( observed=mtfsDzData, type="raw" )
###################################################################
#base models
###################################################################
nl <- 2 # number of latent factors
# Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s)
XMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE,
values=c(.6,
.1, .1),
labels=c("x11MN",
"x21MN", "x22MN"),
lbound=.00001, name="XMN" )
YMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE,
values=c(.6,
.1, .1),
labels=c("y11MN",
"y21MN", "y22MN"),
lbound=.00001, name="YMN" )
ZMN <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE,
values=c(.6,
.1, .1),
labels=c("z11MN",
"z21MN", "z22MN"),
lbound=.00001, name="ZMN" )
AlMN <- mxAlgebra(XMN %*% t(XMN), name="AlMN")
ClMN <- mxAlgebra(YMN %*% t(YMN), name="ClMN")
ElMN <- mxAlgebra(ZMN %*% t(ZMN), name="ElMN")
# Matrices as, cs, and es to store a, c, and e path coefficients for specific factors
pathAsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3,
labels=c("a11MN", "a22MN", "a33MN", "a44MN", "a55MN"),
lbound=.00001, name="asMN" )
pathCsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3,
labels=c("c11MN","c22MN", "c33MN", "c44MN", "c55MN"),
lbound=.00001, name="csMN" )
pathEsMN <- mxMatrix(type="Diag", nrow=nvMN, ncol=nvMN, free=TRUE, values=.3,
labels=c("e11MN", "e22MN", "e33MN", "e44MN", "e55MN"),
lbound=.00001, name="esMN" )
pathFlMN <- mxMatrix(type="Full", nrow=5, ncol=nl, free=FALSE,
labels=c("data.intercept", "data.age14twinpair",
"data.intercept", "data.age17twinpair",
"data.intercept", "data.age20twinpair",
"data.intercept", "data.age24twinpair",
"data.intercept", "data.age29twinpair"), name="flMN", byrow=T)
covAMN <- mxAlgebra( expression=flMN %*% AlMN %*% t(flMN) + asMN %*% t(asMN), name="AMN" )
covCMN <- mxAlgebra( expression=flMN %*% ClMN %*% t(flMN) + csMN %*% t(csMN), name="CMN" )
covEMN <- mxAlgebra( expression=flMN %*% ElMN %*% t(flMN) + esMN %*% t(esMN), name="EMN" )
#estimate means for latent factors
FacMeansMN <- mxMatrix("Full",nrow=nl,ncol=1, free=T, values=1, labels=c("mean1MN", "mean2MN", "mean3MN"), name="FacMeansMN")
meansMN <- mxAlgebra(t(flMN %*% FacMeansMN), name="meansMN")
regCoefMN <- mxMatrix(type="Full", nrow=1, ncol=nvMN, free=TRUE, values=0,
labels=c("beta1MN", "beta2MN", "beta3MN", "beta4MN", "beta5MN"), name="betaMN")
sexCovMN <- mxMatrix(type="Full", nrow=1, ncol=nvMN, free=FALSE, labels="data.sex", name="sexMN")
MuMN <- mxAlgebra(expression= meansMN+(betaMN*sexMN), name="MuMN")
expMeanMN <- mxAlgebra(cbind(MuMN, MuMN), name="expMeanMN")
covMZMN <- mxAlgebra( expression= rbind( cbind(AMN+CMN+EMN, AMN+CMN),
cbind( AMN+CMN, AMN+CMN+EMN)), name="expCovMZMN" )
covDZMN <- mxAlgebra( expression= rbind( cbind(AMN+CMN+EMN, 0.5%x%AMN+CMN),
cbind(0.5%x%AMN+CMN, AMN+CMN+EMN)), name="expCovDZMN" )
### Combine Groups
objMZMN <- mxExpectationNormal( covariance="expCovMZMN", means="expMeanMN", dimnames=selVarsMN )
objDZMN <- mxExpectationNormal( covariance="expCovDZMN", means="expMeanMN", dimnames=selVarsMN )
parsMN <- list( XMN, YMN, ZMN, AlMN, ClMN, ElMN, pathAsMN, pathCsMN, pathEsMN, pathFlMN, covAMN, covCMN,
covEMN, expMeanMN, regCoefMN, sexCovMN, MuMN,
FacMeansMN, meansMN)
funML <- mxFitFunctionML()
mtfsModelMZ <- mxModel( covMZMN, mtfsDataMZ, objMZMN, funML, name="MZMN", parsMN)
mtfsModelDZ <- mxModel( covDZMN, mtfsDataDZ, objDZMN, funML, name="DZMN", parsMN)
#run model
fitML <- mxFitFunctionMultigroup(c("MZMN.fitfunction","DZMN.fitfunction"))
LG <- mxModel("LinearGrowthACE", mtfsModelMZ, mtfsModelDZ, fitML)
mxOption(NULL,"Default optimizer","SLSQP")
LGFit <- mxRun(LG) ### Run model
RefModLG <- mxRefModels(LGFit, run=TRUE)
summary(LGFit, refModels=RefModLG)
LGFit$output$status$code
#examine Output of interest
round(LGFit$output$algebras$MZMN.AlMN, 3)
round(LGFit$output$algebras$MZMN.ClMN, 3)
round(LGFit$output$algebras$MZMN.ElMN, 3)
round(LGFit$output$algebras$MZMN.MuMN, 3)
round(LGFit$output$matrices$MZMN.FacMeansMN, 3)
round(LGFit$output$matrices$MZMN.flMN, 3)
#compute variance of each latent factor and correlate latent factors
VarMN <- LGFit$output$algebras$MZMN.AlMN + LGFit$output$algebras$MZMN.ClMN + LGFit$output$algebras$MZMN.ElMN
round(VarMN, 3)
round(cov2cor(VarMN), 2)
round(cov2cor(LGFit$output$algebras$MZMN.AlMN), 2)
round(cov2cor(LGFit$output$algebras$MZMN.ClMN), 2)
round(cov2cor(LGFit$output$algebras$MZMN.ElMN), 2)
#variance of each factor as a proportion
#intercept
LGFit$output$algebras$MZMN.AlMN[1,1]/(LGFit$output$algebras$MZMN.AlMN[1,1] + LGFit$output$algebras$MZMN.ClMN[1,1] + LGFit$output$algebras$MZMN.ElMN[1,1])
LGFit$output$algebras$MZMN.ClMN[1,1]/(LGFit$output$algebras$MZMN.AlMN[1,1] + LGFit$output$algebras$MZMN.ClMN[1,1] + LGFit$output$algebras$MZMN.ElMN[1,1])
LGFit$output$algebras$MZMN.ElMN[1,1]/(LGFit$output$algebras$MZMN.AlMN[1,1] + LGFit$output$algebras$MZMN.ClMN[1,1] + LGFit$output$algebras$MZMN.ElMN[1,1])
#slope
LGFit$output$algebras$MZMN.AlMN[2,2]/(LGFit$output$algebras$MZMN.AlMN[2,2] + LGFit$output$algebras$MZMN.ClMN[2,2] + LGFit$output$algebras$MZMN.ElMN[2,2])
LGFit$output$algebras$MZMN.ClMN[2,2]/(LGFit$output$algebras$MZMN.AlMN[2,2] + LGFit$output$algebras$MZMN.ClMN[2,2] + LGFit$output$algebras$MZMN.ElMN[2,2])
LGFit$output$algebras$MZMN.ElMN[2,2]/(LGFit$output$algebras$MZMN.AlMN[2,2] + LGFit$output$algebras$MZMN.ClMN[2,2] + LGFit$output$algebras$MZMN.ElMN[2,2])