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.
Your model in script 2 is unidentified, because of this:
You want to fix those variance to 1.0, right? Change
free=T
tofree=F
.I thought I could estimate variances if the first two thresholds are fixed?
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.
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.
See my comment here.
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.
it does differ between twins, but only minimally. yes, maybe i will go full bivariate, but first I have to get this running