Login

News

Latent Growth Curves and Definition Variables in a Twin Sample

Posted on
No user picture. szellers Joined: 04/04/2018
Forums
Behavioral Genetics Models
I am trying to fit a latent growth model in twin data and I am having trouble writing code for definition variables. I want to fix the loading from the intercept to the observed vars to 1 and I want to use age at assessment as the loading for the slope. This model runs when I use just a number corresponding to wave (1, 2, 3, 4, 5) for the slope loadings, but I am having errors when trying to use the definition variable of age at each assessment as the factor loading for the slope. My code here fails at the flValues line.


########################
### Biometric Models ###
########################
###
### Select observed variables
Vars <- c("mj_12m_freq_14", "mj_12m_freq_17", "mj_12m_freq_20", "mj_12m_freq_24", "mj_12m_freq_29")
nv <- length(Vars)
ntv <- nv*2
selVars <- paste(Vars, c(rep(0,nv), rep(1,nv)), sep="")

############
#definition variables
DefVars <- c("age14twinpair", "age17twinpair", "age20twinpair", "age24twinpair", "age29twinpair")

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

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

### ------------------------------------------------------------------------------
### Common Pathway ACE model
### ------------------------------------------------------------------------------
nl <- 2 # 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=c("xm11", "xm21", "xm22"), lbound=.00001, name="Xm" )
Ym <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6,
labels=c("ym11", "ym21", "ym22"), lbound=.00001, name="Ym" )
Zm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6,
labels=c("zm11", "zm21", "zm22"), 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", "asm5"), name="asm" )
pathCsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("csm1", "csm2", "csm3", "csm4", "csm5"), name="csm" )
pathEsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("esm1", "esm2", "esm3", "esm4", "esm5"), 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
flFree <- c(FALSE, FALSE,
FALSE, FALSE,
FALSE, FALSE,
FALSE, FALSE,
FALSE, FALSE)

flValues <- c(1, data.age14twinpair,
1, data.age17twinpair,
1, data.age20twinpair,
1, data.age24twinpair,
1, data.age29twinpair)
###
pathFlm <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=flFree, values=flValues,
lbound=.000001, labels=c("Int1", "Slope1",
"Int2", "Slope2",
"Int3", "Slope3",
"Int4", "Slope4",
"Int5", "Slope5"), name="flm", byrow=T)

# 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","m5m", "m1m", "m2m", "m3m", "m4m", "m5m"), 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 )

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

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')))

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

RefModCom <- mxRefModels(CfFit, run=TRUE)
summary(CfFit, refModels=RefModCom)
summary(RefModCom$Saturated)

Replied on Thu, 04/12/2018 - 11:21
Picture of user. AdminRobK Joined: 01/24/2014

Try changing these lines,

flValues <- c(1, data.age14twinpair,
1, data.age17twinpair,
1, data.age20twinpair,
1, data.age24twinpair,
1, data.age29twinpair)
###
pathFlm <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=flFree, values=flValues,
lbound=.000001, labels=c("Int1", "Slope1",
"Int2", "Slope2",
"Int3", "Slope3",
"Int4", "Slope4",
"Int5", "Slope5"), name="flm", byrow=T)

, to this:

flValues <- c(1, 0,
1, 0,
1, 0,
1, 0,
1, 0)
###
pathFlm <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=flFree, values=flValues,
lbound=.000001, labels=c("Int", "data.age14twinpair",
"Int", "data.age17twinpair",
"Int", "data.age20twinpair",
"Int", "data.age24twinpair",
"Int", "data.age29twinpair"), name="flm", byrow=T)
Log in or register to post comments
Replied on Thu, 04/12/2018 - 12:12
No user picture. szellers Joined: 04/04/2018

Thank you! I am no longer getting an error from the pathFlm matrix line, but now when I try to run the model I get the error "Error: A definition variable 'data.age14twinpair' has been declared in model 'ComACE' that does not contain a data set". The age variables are in both my data objects mtfsDataMz and mtfsDataDz, which are included in my model, so I am uncertain why I am getting this error. Could it be that I have two submodels (MZ and DZ) with their own data objects, and if so, how would I need to specify which data object that each model should use?

Updated code below

########################
### Biometric Models ###
########################
###
### Select observed variables
Vars <- c("mj_12m_freq_14", "mj_12m_freq_17", "mj_12m_freq_20", "mj_12m_freq_24", "mj_12m_freq_29")
nv <- length(Vars)
ntv <- nv*2
selVars <- paste(Vars, c(rep(0,nv), rep(1,nv)), sep="")

############
#definition variables
DefVars <- c("age14twinpair", "age17twinpair", "age20twinpair", "age24twinpair", "age29twinpair")

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

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

### ------------------------------------------------------------------------------
### Common Pathway ACE model
### ------------------------------------------------------------------------------
nl <- 2 # 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=c("xm11", "xm21", "xm22"), lbound=.00001, name="Xm" )
Ym <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6,
labels=c("ym11", "ym21", "ym22"), lbound=.00001, name="Ym" )
Zm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6,
labels=c("zm11", "zm21", "zm22"), 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", "asm5"), name="asm" )
pathCsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("csm1", "csm2", "csm3", "csm4", "csm5"), name="csm" )
pathEsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("esm1", "esm2", "esm3", "esm4", "esm5"), 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
flFree <- c(FALSE, FALSE,
FALSE, FALSE,
FALSE, FALSE,
FALSE, FALSE,
FALSE, FALSE)

flValues <- c(1, 0,
1, 0,
1, 0,
1, 0,
1, 0)
###
pathFlm <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=flFree, values=flValues,
lbound=.000001, labels=c("Int", "data.age14twinpair",
"Int", "data.age17twinpair",
"Int", "data.age20twinpair",
"Int", "data.age24twinpair",
"Int", "data.age29twinpair"), name="flm", byrow=T)

# 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","m5m", "m1m", "m2m", "m3m", "m4m", "m5m"), 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 )

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

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')))

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

RefModCom <- mxRefModels(CfFit, run=TRUE)

Log in or register to post comments
Replied on Thu, 04/12/2018 - 13:35
Picture of user. AdminRobK Joined: 01/24/2014

Thank you! I am no longer getting an error from the pathFlm matrix line, but now when I try to run the model I get the error "Error: A definition variable 'data.age14twinpair' has been declared in model 'ComACE' that does not contain a data set". The age variables are in both my data objects mtfsDataMz and mtfsDataDz, which are included in my model, so I am uncertain why I am getting this error. Could it be that I have two submodels (MZ and DZ) with their own data objects, and if so, how would I need to specify which data object that each model should use?

You're seeing that error because the object with R symbol pathFlm is being put directly into the "container" model (it's an element of list parsm), which has R symbol Cf and does not directly contain an MxData object. Modify your script so that pathFlm only ends up in the submodels (which are named 'MZm' and 'DZm'), and not directly inside the container model.

Log in or register to post comments
  Login

News