You are here

Bivariate Cholesky - Continuous & Ordinal, with Covariates

12 posts / 0 new
Last post
JuliaJoplin's picture
Offline
Joined: 10/12/2014 - 16:44
Bivariate Cholesky - Continuous & Ordinal, with Covariates
AttachmentSize
Plain text icon Postscript.txt4.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'))
)
AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
ACE.expThre

Hi

You put the line

mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=0.5, label="thre", name="expThre"),
 

in the ACE 'container' model, not in the MZ and DZ data models. You could refer to it as

mxFIMLObjective( covariance="ACE.expCovDZ", means="expMean", dimnames=selVars, thresholds="ACE.expThre", threshnames=c('ogo1','ogo2') )    ),
 

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.

JuliaJoplin's picture
Offline
Joined: 10/12/2014 - 16:44
Thanks!

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:

mydata$ogo1 <- mxFactor(mydata$ogo1, levels=c(0:6), ordered=TRUE)
mydata$ogo2 <- mxFactor(mydata$ogo2, levels=c(0:5), ordered=TRUE)

So, above, I tried to tell openMX how many levels my ordinal variables have.

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"),

Below, I tried to add the constraint on the variance that needs to happen for the ordinal variable "ogo."

 mxConstraint(V[1,2]==1, name="Vconst"),
 

Then onward.

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"),    
 

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.

   # 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=7, ncol=2, free=T, values=0.1, label="thre", name="expThre"),
 

And so on...

    # 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" ),

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.)"

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Recommend collapsing categories 5&6

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

mxMatrix(type="Full", nrow=7, ncol=2, free=T, values=0.1, label="thre", name="expThre"),
 

to

mxMatrix(type="Full", nrow=7, ncol=2, free=T, values=0.1, lbound=c(NA,rep(.01,6), label="thre", name="expThreDev"),
mxMatrix(type="Lower", nrow=7, ncol=7, free=F, values=1, name="unitLower"),
mxAlgebra(unitLower %*% expThreDev, name="expThre"),
 

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.

JuliaJoplin's picture
Offline
Joined: 10/12/2014 - 16:44
thanks again!

Thank you so much for both the explanation and the modifications!

Now, I am starting out my code with:

mydata$ogo1 <- cut(mydata$ogo1, breaks=6, ordered_result=T, labels=c(0,1,2,3,4,5))
 
mydata$ogo2 <- cut(mydata$ogo2, breaks=6, ordered_result=T, labels=c(0,1,2,3,4,5))
 
mydata$ogo1 <- mxFactor(mydata$ogo1, levels=c(0:5), ordered=TRUE)
mydata$ogo2 <- mxFactor(mydata$ogo2, levels=c(0:5), ordered=TRUE)

with the understanding that the "cut" function will create 6 breaks in each variable.

I included this:

  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=7, ncol=2, free=T, values=0.1, lbound=c(NA,rep(.01,6), label="thre", name="expThreDev")),
        mxMatrix(type="Lower", nrow=7, ncol=7, free=F, values=1, name="unitLower"),
        mxAlgebra(unitLower %*% expThreDev, name="expThre"),

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"))"

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
parentheses out-of-place

Try changing

mxMatrix(type="Full", nrow=7, ncol=2, free=T, values=0.1, lbound=c(NA,rep(.01,6), label="thre", name="expThreDev")),

to

mxMatrix(type="Full", nrow=7, ncol=2, free=T, values=0.1, lbound=c(NA,rep(.01,6)), label="thre", name="expThreDev"),
JuliaJoplin's picture
Offline
Joined: 10/12/2014 - 16:44
Progress

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)

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
omxAssignFirstParameters()

There's a helper function for that issue. Run your MxModel through omxAssignFirstParameters(), e.g.

bivACEModel <- omxAssignFirstParameters(bivACEModel)
JuliaJoplin's picture
Offline
Joined: 10/12/2014 - 16:44
hmmm

OK - so after running the syntax for bivACEModel, I did

bivACEModel <- omxAssignFirstParameters(bivACEModel)
bivACEFit <- mxRun(bivACEModel)

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?

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Don't label

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.

JuliaJoplin's picture
Offline
Joined: 10/12/2014 - 16:44
Ah ha!!

It ran!

free parameters:
                  name         matrix row col     Estimate lbound ubound
1                  a11          ACE.a   1   1  0.573918262              
2                  a21          ACE.a   2   1  0.859416207              
3                  a22          ACE.a   2   2  1.333587765              
4                  c11          ACE.c   1   1  0.580043867              
5                  c21          ACE.c   2   1  0.654597474              
6                  c22          ACE.c   2   2 -1.149262954              
7                  e11          ACE.e   1   1  0.510430764              
8                  e21          ACE.e   2   1  0.248946777              
9                  e22          ACE.e   2   2  1.610223048              
10         bage_caster    ACE.betaage   1   1  0.360890048              
11         bsex_caster    ACE.betasex   1   1 -0.238669021              
12            bsex_ogo    ACE.betasex   1   2 -0.940551484              
13           bage2_ogo   ACE.betaage2   1   2 -0.004847995              
14             caster1    ACE.expMean   1   1 -6.306266346              
15             caster2    ACE.expMean   1   3 -6.279644224              
16 ACE.expThreDev[1,1] ACE.expThreDev   1   1  1.314117229              
17 ACE.expThreDev[2,1] ACE.expThreDev   2   1  1.003624657   0.01       
18 ACE.expThreDev[3,1] ACE.expThreDev   3   1  0.982194334   0.01       
19 ACE.expThreDev[4,1] ACE.expThreDev   4   1  1.170978906   0.01       
20 ACE.expThreDev[5,1] ACE.expThreDev   5   1  2.695439602   0.01       
21 ACE.expThreDev[6,1] ACE.expThreDev   6   1  0.100000000   0.01       
22 ACE.expThreDev[7,1] ACE.expThreDev   7   1  0.100000000   0.01       
23 ACE.expThreDev[1,2] ACE.expThreDev   1   2  0.556288850              
24 ACE.expThreDev[2,2] ACE.expThreDev   2   2  0.551343467   0.01       
25 ACE.expThreDev[3,2] ACE.expThreDev   3   2  0.820331867   0.01       
26 ACE.expThreDev[4,2] ACE.expThreDev   4   2  1.031453384   0.01       
27 ACE.expThreDev[5,2] ACE.expThreDev   5   2  1.539475496   0.01       
28 ACE.expThreDev[6,2] ACE.expThreDev   6   2  0.100000000   0.01       
29 ACE.expThreDev[7,2] ACE.expThreDev   7   2  0.100000000   0.01  

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.

JuliaJoplin's picture
Offline
Joined: 10/12/2014 - 16:44
hmmm

Looking back at this - did I somehow estimate a different threshold for each twin...?