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.