| Attachment | Size |
|---|---|
 Postscript.txt | 4.57 KB |
I still consider myself a BG novice, and just when I think I've gone as deep as I'll ever go, something else comes along! Awhile ago, the folks on this board graciously helped me do a bivariate Cholesky decomposition of two variables (caster and ogo). The two variables were assessed at different time points, so one variable needed to be adjusted for age1 and sex, the other variable for age2 and sex. We came up with a successful script (attached).
I am now interested in adapting this work so that the "ogo" variable is an ordinal variable, while caster is still continuous. I've read some about thresholds and mxFactor and such, and I've found a number of example scripts that seem to be doing what I'm trying to do (minus the covariates)....still, I know that something is clearly wrong in what I've put together, as I keep getting the error message "Unknown thresholds name 'expThre' detected in the expectation function of model 'MZ'."
mzData <- subset(mydata, ZYGOSITY=="0") dzData <- subset(mydata, ZYGOSITY=="1") # Select Variables for Analysis Vars <- c('caster','ogo') nv <- 2 # number of variables ntv <- nv*2 # number of total variables selVars <- c("caster1", "ogo1", "caster2", "ogo2") bivACEModel <- mxModel("bivACE", mxModel("ACE", # Matrices a, c, and e to store a, c, and e path coefficients mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("a11","a21","a22"), name="a" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("c11","c21","c22"),name="c" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("e11","e21","e22"),name="e" ), # Matrices a, c, and e to store a, c, and e path coefficients mxAlgebra( expression=a %*% t(a), name="A" ), mxAlgebra( expression=c %*% t(c), name="C" ), mxAlgebra( expression=e %*% t(e), name="E" ), # Declare a matrix for the definition variable regression parameters, called beta mxMatrix( type="Full", nrow=1, ncol=2, free=c(TRUE,FALSE), values= c(.4,0), label=c("bage_caster", "bage_ogo"), name="betaage"), mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .1, label=c("bsex_caster", "bsex_ogo"), name="betasex"), mxMatrix( type="Full", nrow=1, ncol=2, free=c(FALSE,TRUE), values= c(0,.1), label=c("bage2_caster", "bage2_ogo"), name="betaage2"), # Matrices A, C, and E compute variance components mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), mxConstraint(V[1,2]==1, name="Vconst"), mxAlgebra(iSD %*% a, name="sta"), mxAlgebra(iSD %*% c, name="stc"), mxAlgebra(iSD %*% e, name="ste"), mxAlgebra(A/V, name="StandardizedA"), mxAlgebra(C/V, name="StandardizedC"), mxAlgebra(E/V, name="StandardizedE"), mxAlgebra(solve(sqrt(I*A)) %&% A, name="CorA"), mxAlgebra(solve(sqrt(I*C)) %&% C, name="CorC"), mxAlgebra(solve(sqrt(I*E)) %&% E, name="CorE"), mxAlgebra(solve(sqrt(I*V)) %&% V, name="CorP"), # Matrix & Algebra for expected means vector mxMatrix( type="Full", nrow=1, ncol=ntv, free=c(T,F,T,F), values=c(.6,0,.6,0), labels=c("caster1","ogo1","caster2","ogo2"), name="expMean"), mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=0.5, label="thre", name="expThre"), # Algebra for expected variance/covariance matrix in MZ mxAlgebra( expression= rbind ( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ), # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) ), mxModel("MZ", mxData( observed=mzData, type="raw" ), # Algebra for making the means a function of the definition variables age and sex mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="obs_age2" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ), mxAlgebra( expression= ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2 , name="expMean"), mxFIMLObjective( covariance="ACE.expCovMZ", means="expMean", dimnames=selVars, thresholds="expThre", threshnames=c('ogo1','ogo2') ) ), mxModel("DZ", mxData( observed=dzData, type="raw" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="obs_age2" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ), mxAlgebra( expression=ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2 , name="expMean"), mxFIMLObjective( covariance="ACE.expCovDZ", means="expMean", dimnames=selVars, thresholds="expThre", threshnames=c('ogo1','ogo2') ) ), mxAlgebra( expression=MZ.objective + DZ.objective, name="minus2sumloglikelihood" ), mxAlgebraObjective("minus2sumloglikelihood"), mxCI(c('ACE.sta', 'ACE.stc', 'ACE.ste')), mxCI(c('ACE.StandardizedA', 'ACE.StandardizedC', 'ACE.StandardizedE')), mxCI(c('ACE.CorA', 'ACE.CorC', 'ACE.CorE', 'ACE.CorP')) )
Printer-friendly version
Hi
You put the line
in the ACE 'container' model, not in the MZ and DZ data models. You could refer to it as
or you could put it in BOTH the MZ and DZ groups without the group identifier prefix ACE. in the mxFIMLObjective() call. Incidentally, it would be better to switch to mxFitFunctionML() and mxExpectationNormal() calls. This would make it easier to switch to WLS, for example, should you wish to try that out.
Thank you! I'm making progress, I hope. I fixed this issue. However, now I am having a problem with specifying the number of thresholds. As I mentioned, my ordinal variable is "ogo." In twin1, values range from 0 to 6 (although only one person endorsed "6"), in twin2, values range from 0 to 5. I'm making notes about places that are different from my original script where I analyzed two continuous variables.
So, to set this up, I now have:
So, above, I tried to tell openMX how many levels my ordinal variables have.
Below, I tried to add the constraint on the variance that needs to happen for the ordinal variable "ogo."
Then onward.
Below, the expected means vector is different than my original script, as I'm fixing the mean of the ordinal variable to some value...then I have the threshold matrix.
And so on...
Threshnames = ogo1 and ogo2
mxAlgebra( expression= ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2 , name="expMean"), mxFIMLObjective( covariance="ACE.expCovMZ", means="expMean", dimnames=selVars, thresholds="ACE.expThre", threshnames=c('ogo1', 'ogo2') ) ), mxModel("DZ", mxData( observed=dzData, type="raw" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="obs_age2" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ), mxAlgebra( expression=ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2 , name="expMean"), mxFIMLObjective( covariance="ACE.expCovDZ", means="expMean", dimnames=selVars, thresholds="ACE.expThre", threshnames=c('ogo1', 'ogo2') ) ), mxAlgebra( expression=MZ.objective + DZ.objective, name="minus2sumloglikelihood" ), mxAlgebraObjective("minus2sumloglikelihood"), mxCI(c('ACE.sta', 'ACE.stc', 'ACE.ste')), mxCI(c('ACE.StandardizedA', 'ACE.StandardizedC', 'ACE.StandardizedE')), mxCI(c('ACE.CorA', 'ACE.CorC', 'ACE.CorE', 'ACE.CorP')) )When I try to run, I get an error message that says "Error: The job for model 'bivACE' exited abnormally with the error message: fit is not finite (Found 4 thresholds too close together in column 2.)"
Hi Julia
I strongly recommend making the categories equal for the two twins. You could make both 6 or both 5, but different numbers of categories for twin 1 and 2 does something that you probably don't want. Technically, OpenMx will calculate the integral from the last threshold to infinity for anyone who is in the topmost category. This means that twin 2 in category 5 will have a likelihood based on the integral from threshold 5 to infinity, but twin 1 in the same category (5) would have a likelihood calculated from threshold 5 to 6. This is bad. Either make both 0-5 or both 0-6. It should not hurt ML estimation to use 0-6 but it could adversely affect WLS estimation, so I tend to favor 0-5 recoding for both.
Second, we recommend trying to keep the thresholds separate and in their correct order with a little matrix algebra trick. If you change
to
then you should be good to go. Note how originally all the thresholds were set to have values of .1 (hence the error message) but the modification makes the expThreDev parameters the first threshold, the deviation from threshold 1 to 2, deviation from 2 to 3 etc. The deviations are bounded to keep them positive, so the thresholds stay in strictly increasing order after the algebra forms the expThre result.
Thank you so much for both the explanation and the modifications!
Now, I am starting out my code with:
with the understanding that the "cut" function will create 6 breaks in each variable.
I included this:
and am now getting this error message:
"Error: 'lbound' argument to mxMatrix function must be of numeric type in mxMatrix(type = "Full", nrow = 7, ncol = 2, free = T, values = 0.1, lbound = c(NA, rep(0.01, 6), label = "thre", name = "expThreDev"))"
Try changing
to
I feel like every new error message is some progress, at least! :) Now I'm getting:
Error: The free parameter 'thre' has been assigned multiple lower bounds! See matrix 'ACE.expThreDev' at location (2, 1) and matrix 'ACE.expThreDev' at location (1, 1)
There's a helper function for that issue. Run your MxModel through
omxAssignFirstParameters(), e.g.OK - so after running the syntax for bivACEModel, I did
and I am still getting the error message - "Error: The free parameter 'thre' has been assigned multiple lower bounds! See matrix 'ACE.expThreDev' at location (2, 1) and matrix 'ACE.expThreDev' at location (1, 1)"
Could this have something to do with my version of openmx?
My bad, perhaps. Note how there's a single labels='thre' part. If a list isn't long enough then R re-uses the list from the beginning, so all your parameters ended up constrained to be equal. Since that is not what you want, just get rid of that argument from the function call.
It ran!
Does this output look right, though? It looks like I am estimating a different mean for "caster1" than for "caster2", and I have two sets of thresholds....."caster" is my continuous variable and "ogo" is my ordinal variable.
Looking back at this - did I somehow estimate a different threshold for each twin...?