# Bivariate ACE model with covariates for ordinal variables

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

## data reserved word

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

## full script?

Log in or register to post comments

In reply to full script? by AdminRobK

## Full script attached

Please find attached my R script.

Thanks for your help,

Yi Ting

Log in or register to post comments

## R script attached

Please find attached my R script.

Thanks for your help,

Yi Ting

Log in or register to post comments

## Kronecker products; traceback?

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

In reply to Kronecker products; traceback? by AdminRobK

## 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

In reply to Hi Rob, by YiTan

## yes; thanks

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

Log in or register to post comments

## I think I know what the

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

## mxVersion()?

Log in or register to post comments

In reply to mxVersion()? by AdminRobK

## 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

In reply to Hi Rob, by YiTan

## increase number of threads

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

## fixed for future versions

Log in or register to post comments