Login

News

Bivariate Cholesky - Continuous & Ordinal, with Covariates

Posted on
No user picture. JuliaJoplin Joined: 10/12/2014
Forums
Behavioral Genetics Models
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'))
)

Replied on Tue, 03/07/2017 - 09:35
Picture of user. AdminNeale Joined: 03/01/2013

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.

Log in or register to post comments
Replied on Tue, 03/07/2017 - 12:00
No user picture. JuliaJoplin Joined: 10/12/2014

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

Log in or register to post comments
Replied on Tue, 03/07/2017 - 14:20
Picture of user. AdminNeale Joined: 03/01/2013

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.

Log in or register to post comments
Replied on Tue, 03/07/2017 - 17:08
No user picture. JuliaJoplin Joined: 10/12/2014

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

Log in or register to post comments
Replied on Tue, 03/07/2017 - 18:03
Picture of user. AdminRobK Joined: 01/24/2014

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"),
Log in or register to post comments
Replied on Tue, 03/07/2017 - 18:12
No user picture. JuliaJoplin Joined: 10/12/2014

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)

Log in or register to post comments
Replied on Wed, 03/08/2017 - 13:34
No user picture. JuliaJoplin Joined: 10/12/2014

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?

Log in or register to post comments
Replied on Wed, 03/08/2017 - 16:15
Picture of user. AdminNeale Joined: 03/01/2013

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.
Log in or register to post comments
Replied on Wed, 03/08/2017 - 17:49
No user picture. JuliaJoplin Joined: 10/12/2014

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.

Log in or register to post comments
  Login

News