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 ofSex
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 yourmxVersion()
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.