Hi.
I've been trying everything to get an ordinal manifest working in RAM models but I can't get it done. Also, didn't find anything online about this, as most examples have matrix specification. Would be very appreciated to get some help. Besides from the solution below, I have also tried every possible way with mxthreshold and umxthresholdmatrix.
obesityLevels = c('normal', 'obese') cutPoints = quantile(twinData[, "bmi1"], probs = .2, na.rm = TRUE) twinData$obese1 = cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) twinData$obese2 = cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) # Step 2: Make the ordinal variables into umxFactors (ordered, with the levels found in the data) selVars = c("obese1", "obese2") twinData[, selVars] = umxFactor(twinData[, selVars]) #tmp = umxThresholdMatrix(twinData, selDVs = selVars, sep = "", verbose = TRUE) tmp = umx_make_twin_data_nice(data=twinData, sep="", zygosity="zygosity", numbering=1:2) tmp = umx_scale_wide_twin_data(varsToScale= c("wt", "ht"), sep= "_T", data= tmp) mzData = subset(tmp, zygosity %in% c("MZFF", "MZMM")) dzData = subset(tmp, zygosity %in% c("DZFF", "DZMM")) latents = c("A1", "C1", "E1", "A2", "C2", "E2") manifests = c("obese_T1", "obese_T2") latVariances <- mxPath( from=latents, arrows=2, free=FALSE, values=1 ) manifestsVariances <- mxPath( from=manifests, arrows=2, free=FALSE, values=1 ) # means of latent variables latMeans <- mxPath( from="one", to=latents, arrows=1, free=FALSE, values=0 ) # means of observed variables obsMeans <- mxPath( from="one", to=manifests, arrows=1, free=F, values=0, labels="mean" ) # path coefficients for twin 1 pathAceT1 <- mxPath( from=c("A1","C1","E1"), to="obese_T1", arrows=1, free=TRUE, values = 1, label=c("a","c","e") ) # path coefficients for twin 2 pathAceT2 <- mxPath( from=c("A2","C2","E2"), to="obese_T2", arrows=1, free=TRUE, values = 1, label=c("a","c","e") ) # 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 ) paths <- list( latVariances, latMeans, obsMeans, pathAceT1, pathAceT2, covC1C2, manifestsVariances ) dataMZ <- mxData(observed=mzData, type = "raw") dataDZ <- mxData(observed=dzData, type = "raw") #threshold <- mxThreshold(vars=c("obese_T1", "obese_T2"), nThresh=c(1), free = T) threshold <- mxMatrix("Full", nrow=1, ncol=2, byrow=T, name="thresh", dimnames = list(c(), manifests), free = c(T, T), values = c(-.8, -.8), ) expect <- mxExpectationRAM(A="A", S="S", F="F", M="M", thresholds = "thresh") fitfunction = mxFitFunctionML() modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=manifests, latentVars=latents, paths, covA1A2_MZ, dataMZ, threshold, expect, fitfunction) modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=manifests, latentVars=latents, paths, covA1A2_DZ, dataDZ, threshold, expect, fitfunction) obj <- mxFitFunctionMultigroup(c("MZ", "DZ")) modelACE <- mxModel("ACE", modelMZ, modelDZ, obj) fitACE <- mxRun(modelACE)
I have tried fixing the threshold across twins and zygosity, but to no avail. The results either are extremely high estimates for the paths with high SEs, or they are just plainly wrong (if compared to umx results). Help is great appreciated,
thank you
I think the issue is that your variances don't add up. You probably need to free the variances and add a constraint to target the ordinal variances at 1.0. If I fit these data with umx,
you can inspect the content of the model and find such a constraint,
Thank you for your reply. I've added this now:
This addition is successful in constraining the variance to one. However, results are still vastly wrong and the standard errors can't be estimated. You mentioned freeing variances. I have tried that to no avail. Is there anything else to do? Thank you
P.S. the model is easily replicable as it's simply the twinData dataset in umx/OpenMx
I got it working now. In addition to the Mxconstraint, I had to delete the "manifestsVariances" path. What is the explanation for this? I just want to be sure to understand. Thank you
This is the working script:
> What is the explanation for this?
You are accounting for the error variance in the E component. There is no additional error variance to account for or estimate.
Thank you for your helpful reply. When switching to the bivariate case (umxACEv style), how do I put in the constraint? Particularly if I add another continuous manifest, everything blows up.
This is my script: I have followed the documentation on Joint Ordinal continuous manifests.
I have been continuously working on it throughout the weekend but haven't been able to get it to work. So my aim is to build a bivariate ACE Ram model using the variance components approach. For simplicity I selected two ordinal variables now (2 for each twin).
In a first step I just want to replicate this:
which works excellently.
What I think I can understand is the following:
- mean is not fixed
- variance is not fixed
- first two thresholds are fixed but third one is estimated
- paths from latents to manifests are fixed at 1
So I have tried integrating this, as well as trying different combinations, but they all results either in bogus results or errors. I have tried estimating paths, fixing means, fixing variance, different thresholds and combinations thereof but nothing is working. Any ideas?
Here is the code: Thank you very much
I don't see anyplace in your syntax where you put the MxConstraints into an
mxModel()
statement. They won't do anything if they're sitting in R's workspace, and aren't part ofmodelACE
's namespace. I'd suggest creating them before creatingmodelMZ
andmodelDZ
, and putting them into only one of those two MxModels.Edit: well, I guess that doesn't apply if you're trying to identify the model by fixing thresholds. If you're going to fix two of the thresholds, you can freely estimate both the means and the variances of the ordinal phenotypes. Also, when people fix two thresholds, they usually fix them to 0 and 1 (but that's arbitrary).
Specifically what "bogus results or errors" are you getting?
Hi Rob,
thanks for your reply.
I'm aware that I have to put the mxconstraints into the actual mxmodel, they were just leftovers from a previous version. I have tried two versions:
Method 2 produces the error message "In model 'ACE' Optimizer returned a non-zero status code 5. The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED). "
mxcheckidentification:
If I read correctly,
1. You're trying to implement an ACE model in path specification.
2. You're struggling with thresholds.
In answer to 1, The
umxTwinMaker
function is really helpful for this; You just describe the model for one twin, and the relationship between the twins (defaults mostly work for this) and you get back a two-group twin model implementing the proper paths connecting twins in a pair and differences in relationships between pairs. It's great!ordinal variables and thresholds introduce a handful of things to handle, with different decisions to be made for binary, ordinal, continuous, the thresholds matrix, offsets matrix, algebras, and implications to keep straight for constrained latent variables in the model that will map out to your thresholds. Additional concerns about modelling covariates on the means etc. etc. Can give head aches.
umxACE and umxACEv take care of all of this for you.
Hi Tim,
thank you and umx is very much appreciated. Need to understand the backend, though, before publishing anything.
gotcha: I've not followed this long thread closely. But might be useful to look inside the xmuRAM2Ordinal and umxThresholdMatrix functions to see how factor variables with 2, 3, and more levels are handled, creating lowerOnes_for_thresh, deviations_for_thresh, and thresholdsAlgebra. Also see the implications for the model implementation that depend on these - in terms of fixed variance and means in dummy latent variables. You could adapt that approach to your long-form script.