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

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

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.

Hi Rob,

Please find attached my R script.

Thanks for your help,

Yi Ting

Hi Rob,

Please find attached my R script.

Thanks for your help,

Yi Ting

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

is equivalent to this,

. 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?

Hi Rob,

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

Here's what I get from traceback():

Yes, correct. Thank you for the

`backtrace()`

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

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.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?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:

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

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.For the record, I should mention that future versions of OpenMx will give a more-helpful error message in cases like this. See here.