You are here

Bivariate ACE model with covariates for ordinal variables

13 posts / 0 new
Last post
YiTan's picture
Offline
Joined: 10/08/2014 - 05:28
Bivariate ACE model with covariates for ordinal variables

Hi,

I couldn't find any existing OpenMx codes to conduct bivariate genetic modelling for ordinal variables with covariates, so I've adapted Hermine Mae's twoACEvo.R code (bivariate ACE model for ordinal variables) by adding covariates to the code. I thought it was quite straightforward, but when I ran the code I got the following error:

Error in as.vector(data) : 
  no method for coercing this S4 class to a vector

If it helps, I have included a segment of my code with how I incorporated the covariates below. Any advice or comments will be greatly appreciated (or perhaps point me to some other relevant codes to adapt)! Thanks!

# Matrix for moderating/interacting variable
defSex    <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
                       labels=c("data.Sex1","data.Sex2"), name="Sex")
 
# Matrices declared to store linear Coefficients for covariate
B_Sex     <- mxMatrix( type="Full", nrow=nth, ncol=2, free=c(F,T), 
                       values= c(0,.1), label=c('bSexV1','bSexV2'), name="bSex", byrow = TRUE)
 
meanSex   <- mxAlgebra(  bSex%x%Sex, name="SexR")
 
#age
defAge   <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
                      labels=c("data.Age1","data.Age2"), name="Age")
 
# Matrices declared to store linear Coefficients for covariate
B_Age     <- mxMatrix( type="Full", nrow=nth, ncol=2, free=c(T,F), 
                       values= c(.01,0), label=c('bAgeV1','bAgeV2'), name="bAge", byrow = TRUE)
 
meanAge   <- mxAlgebra(  bAge%x%Age, name="AgeR")
 
 
#YrsEd
defYEd   <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
                      labels=c("data.yrsEd1","data.yrsEd2"), name="YEd")
 
# Matrices declared to store linear Coefficients for covariate
B_YEd     <- mxMatrix( type="Full", nrow=nth, ncol=2, free=c(T,F),
                       values= c(.01,0), labels=c('bYEdV1','bYEdV2'), name="bYEd", byrow = TRUE)
 
meanYEd   <- mxAlgebra(  bYEd%x%YEd, name="YEdR")
 
#Age-related hearing condition
defHearing   <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
                          labels=c("data.AHearing1","data.AHearing2"), name="Hearing")
 
# Matrices declared to store linear Coefficients for covariate
B_Hearing     <- mxMatrix( type="Full", nrow=nth, ncol=2, free=TRUE, 
                           values= .01, labels=c('bAHearV1','bAHearV2'), name="bHearing", byrow = TRUE)
 
meanHearing   <- mxAlgebra(  bHearing%x%Hearing, name="HearingR")
#***************************************************************
 
defs      <- list( defSex, B_Sex, meanSex, defAge, B_Age, meanAge, defYEd, B_YEd, meanYEd, defHearing, B_Hearing, meanHearing)
 
#setting up the regression
# Matrix & Algebra for expected means vector and expected thresholds
meanG    <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" )
threG    <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=TRUE, values=svTh, lbound=lbTh, labels=labThZ, name="Thre" )
inc      <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" )
threT    <- mxAlgebra( expression = Inc %*% Thre, name="expThre" )
threC <- mxAlgebra( expression = expThre + AgeR + SexR + YEdR + HearingR, name = "expThreC")  #with covariates
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=c('Vars1','PVars1','Vars2','PVars2'), thresholds="expThreC" )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=c('Vars1','PVars1','Vars2','PVars2'), thresholds="expThreC" )
funML     <- mxFitFunctionML()
AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
data reserved word

Hi

data() is a function in R - I suspect that changing the name of your datafile from data to, .e.g., mydata may solve the issue. There doesn't seem to be an mxData() function call in your code. I am also unsure as to why you want your data to be a vector - mxData takes matrices and data frames.

HTH
Mike

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
full script?

YiTan, there's not enough information in your OP to diagnose what's wrong. I doubt you're getting that error message because of how you're including the covariates. We'd need to see your full script.

YiTan's picture
Offline
Joined: 10/08/2014 - 05:28
Full script attached

Hi Rob,

Please find attached my R script.

Thanks for your help,
Yi Ting

File attachments: 
YiTan's picture
Offline
Joined: 10/08/2014 - 05:28
R script attached

Hi Rob,

Please find attached my R script.

Thanks for your help,
Yi Ting

File attachments: 
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
Kronecker products; traceback?

You need to reverse the order of the factors in your Kronecker products. For example, this algebra,

meanSex   <- mxAlgebra(  bSex%x%Sex, name="SexR")

is equivalent to this,

meanSex   <- mxAlgebra(
  cbind( bSex[1,1]*Sex[1,1], bSex[1,1]*Sex[1,2], bSex[2,1]*Sex[1,1], bSex[2,1]*Sex[1,2] ),
  name="SexR")

. That's backwards, assuming that the elements of bSex are respectively the sex effects on trait 1 and on trait 2, and that the elements of Sex are respectively the sex of twin 1 and twin 2 (is that the case?).

Concerning the error message in your OP: what do you get from traceback() right after you see the message?

YiTan's picture
Offline
Joined: 10/08/2014 - 05:28
Hi Rob,

Hi Rob,

Thanks for your reply. Just to confirm, do you mean I should do the following for my Kronecker products instead:

meanSex   <- mxAlgebra(  Sex%x%bSex, name="SexR")

Here's what I get from traceback():

> traceback()
17: as.vector(data)
16: array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), 
        NULL) else NULL)
15: as.matrix.default(result)
14: as.matrix(result)
13: EvalInternal(expression, model, modelvariable, compute, show, 
        defvar.row, cache, cacheBack, 2L)
12: mxEval(SD, flatModel, compute = TRUE)
11: eval(substitute(mxEval(x, flatModel, compute = TRUE), list(x = as.symbol(entityName))))
10: eval(substitute(mxEval(x, flatModel, compute = TRUE), list(x = as.symbol(entityName))))
9: FUN(X[[i]], ...)
8: lapply(flatModel@intervals, generateIntervalListHelper, flatModel, 
       modelname, parameters, labelsData)
7: generateIntervalList(flatModel, model@name, parameters, labelsData)
6: runHelper(model, frontendStart, intervals, silent, suppressWarnings, 
       unsafe, checkpoint, useSocket, onlyFrontend, useOptimizer, 
       beginMessage)
5: mxRun(BivAceModel, intervals = T) at #3
4: doTryCatch(return(expr), name, parentenv, handler)
3: tryCatchOne(expr, names, parentenv, handlers[[1L]])
2: tryCatchList(expr, classes, parentenv, handlers)
1: tryCatch(expr = {
       mxRun(BivAceModel, intervals = T)
   }, warning = function(w) {
       mxTryHard(BivAceModel, extraTries = 20, bestInitsOutput = T, 
           intervals = T)
   })
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
yes; thanks
Thanks for your reply. Just to confirm, do you mean I should do the following for my Kronecker products instead:

meanSex   <- mxAlgebra(  Sex%x%bSex, name="SexR")

Yes, correct. Thank you for the backtrace().

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
I think I know what the

I think I know what the problem is now. In your script, you're requesting confidence intervals for elements of 'SD':

ci <- mxCI (c("SA[1,1]","SA[2,2]","SD[1,1]","SD[2,2]","SC[1,1]","SC[2,2]","SE[1,1]","SE[2,2]","rA[2,1]","rC[2,1]","rE[2,1]","corrA","corrC","corrE"))

But, you never put an algebra named 'SD' into your MxModel. You define it in your script, but in terms of the dominance-genetic covariance matrix 'D', which is never created in the first place. Try deleting the references to 'SD' from your mxCI() statement and see if that makes a difference.

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

Admittedly, the error message you're seeing isn't very informative. I would have thought OpenMx would have checked for the existence of what's in the mxCI() reference. What is your mxVersion() output?

YiTan's picture
Offline
Joined: 10/08/2014 - 05:28
Hi Rob,

Hi Rob,

Thanks for your help, the code finally ran properly without error after I removed the 'SD' from mxCI.
My mxVersion output is as follows:

OpenMx version: 2.15.4 [GIT v2.15.4]
R version: R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 
MacOS: 10.14.6
Default optimizer: NPSOL
NPSOL-enabled?: Yes
OpenMP-enabled?: Yes

Btw, it took 40 minutes just to run the entire code. I noticed that the codes (whether univariate or bivariate) involving ordinal variables tend to take longer to run compared to continuous variables. Is there any way I can speed that up?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
increase number of threads
Btw, it took 40 minutes just to run the entire code. I noticed that the codes (whether univariate or bivariate) involving ordinal variables tend to take longer to run compared to continuous variables. Is there any way I can speed that up?

Really, the only thing you can do is increase the number of threads OpenMx uses. If you want to maximize the number of threads, put

mxOption(key='Number of Threads', value=parallel::detectCores())

in your script, after you load OpenMx, but before any calls to mxRun().

Edit: Actually, come to think of it, there are two other things you could do, but they might not be good ideas. One would be to increase mxOption "mvnRelEps", and the other would be to decrease any of the "mvnMaxPoints*" mxOptions. The reason why they might be bad ideas is that they will reduce the numerical accuracy of likelihood evaluation for ordinal data.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
fixed for future versions

For the record, I should mention that future versions of OpenMx will give a more-helpful error message in cases like this. See here.