Login

News

Bivariate ACE model with covariates for ordinal variables

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

Replied on Fri, 01/31/2020 - 10:08
Picture of user. AdminNeale Joined: 03/01/2013

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

Log in or register to post comments
Replied on Fri, 01/31/2020 - 12:01
Picture of user. AdminRobK Joined: 01/24/2014

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.
Log in or register to post comments
Replied on Wed, 02/05/2020 - 12:43
Picture of user. AdminRobK Joined: 01/24/2014

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?

Log in or register to post comments
Replied on Mon, 02/10/2020 - 08:03
No user picture. YiTan Joined: 10/08/2014

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

Log in or register to post comments
Replied on Mon, 02/10/2020 - 10:06
Picture of user. AdminRobK Joined: 01/24/2014

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()`.

Log in or register to post comments
Replied on Mon, 02/10/2020 - 10:23
Picture of user. AdminRobK Joined: 01/24/2014

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.
Log in or register to post comments
Replied on Mon, 02/10/2020 - 10:26
Picture of user. AdminRobK Joined: 01/24/2014

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?
Log in or register to post comments
Replied on Tue, 02/11/2020 - 00:12
No user picture. YiTan Joined: 10/08/2014

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?

Log in or register to post comments
Replied on Tue, 02/11/2020 - 14:55
Picture of user. AdminRobK Joined: 01/24/2014

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.

Log in or register to post comments
Replied on Thu, 02/20/2020 - 10:31
Picture of user. AdminRobK Joined: 01/24/2014

For the record, I should mention that future versions of OpenMx will give a more-helpful error message in cases like this. See here.
Log in or register to post comments
  Login

News