Login

News

GxE with ordinal manifest in RAM

Posted on
No user picture. Leo Joined: 01/09/2020
Forums
Behavioral Genetics Models
Hi,

I'm trying to build a Purcell style GxE model in RAM. I have not found any script for this so I created it from scratch. Right now I have a two ordinal manifests (twin 1, twin 2, same phenotype, 4 cats, 3threshs). The moderator will be continuous in the end. However, right now for experimental reasons I have created a binary moderator (sex) for me to easily check if my model does produce accurate results.

I have build two models: 1) All thresholds free, two mxconstraints, means=0, variances=1 and 2) 2 thresholds fixed (0,1), third one free, mean free, variance free, no constraint.

1) is working almost perfectly, at least regarding the results. However, I don't really want to rely on this method as it is kind of messy. And I'm really unsure if I set all the ubounds, lbounds and constraints right. All I can really tell is the results are right.

2) produces sometimes similar results as 1) but does differ wildly often. Particularly the standard errors go crazy, sometimes as much as 500x estimate.

I was wondering if I could get some feedback what would have to be changed.

Here are both scripts:

1):


ordDVs = c("pop0100_w1_T1", "pop0100_w1_T2")

cgr2[, ordDVs] = mxFactor(cgr2[, ordDVs], 1:4)

# subset MZ DZ data
mzData = subset(cgr2, zyg0102w1 == 1)
dzData = subset(cgr2, zyg0102w1 == 2)

latents = c("A1", "C1", "E1", "A2", "C2", "E2")
manifests = ordDVs

latVariances <- mxPath( from=latents, arrows=2,
free=F, values=1)

# means of latent variables
latMeans <- mxPath( from="one", to=latents, arrows=1,
free=F, values=0 )
# means of observed variables
obsMeans <- mxPath( from="one", to=manifests, arrows=1,
free=F, values=0, labels="mean[1,1]")
# path coefficients for twin 1
pathAceT1 <- mxPath( from=c("A1","C1","E1"), to="pop0100_w1_T1", arrows=1,
free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))
# path coefficients for twin 2
pathAceT2 <- mxPath( from=c("A2","C2","E2"), to="pop0100_w1_T2", arrows=1,
free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))
# covariance between C1 & C2
covC1C2 <- mxPath( from="C1", to="C2", arrows=2,
free=FALSE, values=1 )
# covariance between A1 & A2 in MZ twins
covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2,
free=FALSE, values=1 )
# covariance between A1 & A2 in DZ twins
covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2,
free=FALSE, values=.5 )

beta1 <- mxMatrix(type = "Full", 1,1, free = F, labels = "beta1", values = 0)
mean1 <- mxMatrix(type = "Full", 1,1, free = F, labels = "mean1", values = 0)
blub <- mxAlgebra(name="mean", expression = mean1+beta1*data.sex_w1_T1)

a <- mxMatrix(type = "Full", 1,1, free = T, labels = "a", lbound = 0)
beta2 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta2", ubound = 1, lbound = -1)
patha <- mxAlgebra(name="patha", expression = a+beta2*data.sex_w1_T1)

c <- mxMatrix(type = "Full", 1,1, free = T, labels = "c", lbound = 0)
beta3 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta3", ubound = 1, lbound = -1)
pathc <- mxAlgebra(name="pathc", expression = c+beta3*data.sex_w1_T1)

e <- mxMatrix(type = "Full", 1,1, free = T, labels = "e", lbound = 0)
beta4 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta4", ubound = 1, lbound = -1)
pathe <- mxAlgebra(name="pathe", expression = e+beta4*data.sex_w1_T1)

paths <- list( latVariances, latMeans, obsMeans,
pathAceT1, pathAceT2, covC1C2, blub, beta1, mean1,
a, beta2, patha, c, beta3, pathc, e, beta4, pathe)

dataMZ <- mxData(observed=mzData, type = "raw")
dataDZ <- mxData(observed=dzData, type = "raw")

threshold <- mxThreshold(vars=c("pop0100_w1_T1", "pop0100_w1_T2"), nThresh=c(3, 3), free = c(T, T, T),
labels = c("true1", "true2", "true3"), values = c(0,1,2))

A <- mxAlgebra( expression=patha*patha, name="A2" )
C <- mxAlgebra( expression=pathc*pathc, name="C2" )
E <- mxAlgebra( expression=pathe*pathe, name="E2" )
V <- mxAlgebra( expression=A2+C2+E2, name="V" )
a2 <- mxAlgebra( expression=A2/V, name="a2")
c2 <- mxAlgebra( expression=C2/V, name="c2")
e2 <- mxAlgebra( expression=E2/V, name="e2")

constraint <- mxConstraint(a*a+c*c+e*e == 1, name = "constraintordinal")
constraint2 <- mxConstraint((beta2+a)*(beta2+a)+(beta3+c)*(beta3+c)+(beta4+e)*(beta4+e) == 1, name = "constraintordinal2")

modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=manifests,
latentVars=latents, paths, covA1A2_MZ, dataMZ, threshold,
A, C, E, V, a2, c2, e2, constraint, constraint2)
modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=manifests,
latentVars=latents, paths, covA1A2_DZ, dataDZ, threshold,
A, C, E, V, a2, c2, e2)

obj <- mxFitFunctionMultigroup(c("MZ", "DZ"))

modelACE <- mxModel("ACE", modelMZ, modelDZ, obj)
fitACE <- mxRun(modelACE)

summary(fitACE)

script 2:


ordDVs = c("pop0100_w1_T1", "pop0100_w1_T2")

cgr2[, ordDVs] = mxFactor(cgr2[, ordDVs], 1:4)

# subset MZ DZ data
mzData = subset(cgr2, zyg0102w1 == 1)
dzData = subset(cgr2, zyg0102w1 == 2)

latents = c("A1", "C1", "E1", "A2", "C2", "E2")
manifests = ordDVs

latVariances <- mxPath( from=latents, arrows=2,
free=T, values=1, labels = c("AV", "CV", "EV", "AV", "CV", "EV"))

# means of latent variables
latMeans <- mxPath( from="one", to=latents, arrows=1,
free=F, values=0 )
# means of observed variables
obsMeans <- mxPath( from="one", to=manifests, arrows=1,
free=F, values=0, labels="mean[1,1]")
# path coefficients for twin 1
pathAceT1 <- mxPath( from=c("A1","C1","E1"), to="pop0100_w1_T1", arrows=1,
free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))
# path coefficients for twin 2
pathAceT2 <- mxPath( from=c("A2","C2","E2"), to="pop0100_w1_T2", arrows=1,
free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))
# covariance between C1 & C2
covC1C2 <- mxPath( from="C1", to="C2", arrows=2,
free=FALSE, values=1 )
# covariance between A1 & A2 in MZ twins
covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2,
free=FALSE, values=1 )
# covariance between A1 & A2 in DZ twins
covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2,
free=FALSE, values=.5 )

beta1 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta1", values = 0)
mean1 <- mxMatrix(type = "Full", 1,1, free = T, labels = "mean1", values = 0)
blub <- mxAlgebra(name="mean", expression = mean1+beta1*data.sex_w1_T1)

a <- mxMatrix(type = "Full", 1,1, free = T, labels = "a", lbound = 0)
beta2 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta2", ubound = 1, lbound = -1)
patha <- mxAlgebra(name="patha", expression = a+beta2*data.sex_w1_T1)

c <- mxMatrix(type = "Full", 1,1, free = T, labels = "c", lbound = 0)
beta3 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta3", ubound = 1, lbound = -1)
pathc <- mxAlgebra(name="pathc", expression = c+beta3*data.sex_w1_T1)

e <- mxMatrix(type = "Full", 1,1, free = T, labels = "e", lbound = 0)
beta4 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta4", ubound = 1, lbound = -1)
pathe <- mxAlgebra(name="pathe", expression = e+beta4*data.sex_w1_T1)

paths <- list( latVariances, latMeans, obsMeans,
pathAceT1, pathAceT2, covC1C2, blub, beta1, mean1,
a, beta2, patha, c, beta3, pathc, e, beta4, pathe)

dataMZ <- mxData(observed=mzData, type = "raw")
dataDZ <- mxData(observed=dzData, type = "raw")

threshold <- mxThreshold(vars=c("pop0100_w1_T1", "pop0100_w1_T2"), nThresh=c(3, 3), free = c(F, F, T),
labels = c("true1", "true2", "true3"), values = c(0,1,2))

A <- mxAlgebra( expression=patha*patha, name="A2" )
C <- mxAlgebra( expression=pathc*pathc, name="C2" )
E <- mxAlgebra( expression=pathe*pathe, name="E2" )
V <- mxAlgebra( expression=A2+C2+E2, name="V" )
a2 <- mxAlgebra( expression=A2/V, name="a2")
c2 <- mxAlgebra( expression=C2/V, name="c2")
e2 <- mxAlgebra( expression=E2/V, name="e2")

modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=manifests,
latentVars=latents, paths, covA1A2_MZ, dataMZ, threshold,
A, C, E, V, a2, c2, e2)
modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=manifests,
latentVars=latents, paths, covA1A2_DZ, dataDZ, threshold,
A, C, E, V, a2, c2, e2)

obj <- mxFitFunctionMultigroup(c("MZ", "DZ"))

modelACE <- mxModel("ACE", modelMZ, modelDZ, obj)
fitACE <- mxRun(modelACE)

summary(fitACE)

tryhard = mxTryHardOrdinal(modelACE, exhaustive = T, extraTries = 50)

summary(tryhard)

any help would be much appreciated.

Replied on Tue, 07/28/2020 - 11:00
Picture of user. AdminRobK Joined: 01/24/2014

Your model in script 2 is unidentified, because of this:

latVariances <- mxPath( from=latents, arrows=2,
free=T, values=1, labels = c("AV", "CV", "EV", "AV", "CV", "EV"))

You want to fix those variance to 1.0, right? Change `free=T` to `free=F`.
Log in or register to post comments
Replied on Thu, 07/30/2020 - 10:59
Picture of user. AdminRobK Joined: 01/24/2014

You can freely estimate the variances of the *manifest* variables if you fix 2 thresholds. You're parameterizing the variances of the manifest variables in terms of one-headed paths going from latent variables to manifest variables. The model won't be identified unless the variances of the latent variables are fixed. Alternately, you could fix those one-headed paths, and free the latent-variable variances. But, you can't free both sets of parameters.
Log in or register to post comments
Replied on Thu, 07/30/2020 - 18:13
No user picture. Leo Joined: 01/09/2020

ah okay, thank you, now I understand. And what is the difference to the free-all-thresholds + mxconstraint method? The results are the same, but I can estimate all thresholds but in turn have to put in a constraint? I have read about it often on this forum, but I am not sure when to choose one over the other.
Log in or register to post comments
Replied on Tue, 07/28/2020 - 11:11
Picture of user. AdminRobK Joined: 01/24/2014

BTW, what is the moderator variable? There's a good chance you ought to be using the so-called "full bivariate" form of the Purcell model. See my previous comments here and here.
Log in or register to post comments
  Login

News