You are here

GxE with ordinal manifest in RAM

8 posts / 0 new
Last post
Leo's picture
Leo
Offline
Joined: 01/09/2020 - 14:36
GxE with ordinal manifest in RAM

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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
model identification

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.

Leo's picture
Leo
Offline
Joined: 01/09/2020 - 14:36
I thought I could estimate

I thought I could estimate variances if the first two thresholds are fixed?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
manifest variables

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.

Leo's picture
Leo
Offline
Joined: 01/09/2020 - 14:36
ah okay, thank you, now I

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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
See my comment here.

See my comment here.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
moderator?

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.

Leo's picture
Leo
Offline
Joined: 01/09/2020 - 14:36
it does differ between twins,

it does differ between twins, but only minimally. yes, maybe i will go full bivariate, but first I have to get this running