Hi!
I have a script from SGDP summer school, I can analyze using 2 continuous variables with this script. However, when I try using 3 variables (2 continuous 1 dummy coded), I get stuck while setting up the matrix. Before starting the analysis, I standardized my variables so that the SDs were 1 and their mean was 0.
Part of my code is as follows:
nv <- 3 # number of variables for a twin ntv <- 2*nv # number of variables for a pair Vars <- c('bmi', 'age','dummy_x') selVars <- c('bmi1','age1', 'dummy_x1', 'bmi2','age2', 'dummy_x2') # Select Data for Analysis mzData <- subset(TwinData, zygosity=='MZ', selVars) dzData <- subset(TwinData, zygosity=='DZ', selVars) # Descriptive Statistics of data by zygosity group describe(mzData) describe(dzData) colMeans(mzData,na.rm=TRUE) cov(mzData,use="complete") cor(mzData,use="complete") colMeans(dzData,na.rm=TRUE) cov(dzData,use="complete") cor(dzData,use="complete") # ----------------------------------------------------------------------------------------------------------------------------------- # Now we estimate the MZ and DZ covariances and Means by fitting a 'model' to the data # 1a) Specify Bivariate Saturated Model (Cholesky Decomposition) # ----------------------------------------------------------------------------------------------------------------------------------- # Matrix & Algebra for expected means and covariances MZlow <-mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=T, values=5, name="LowMZ" ) CovMZ <-mxAlgebra( expression=LowMZ %*% t(LowMZ), name="ExpCovMZ" ) MeanMZ <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=c(0,0,0,0,0,0), name="ExpMeanMZ" ) DZlow <-mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=T, values=5, name="LowDZ" ) CovDZ <-mxAlgebra( expression=LowDZ %*% t(LowDZ), name="ExpCovDZ" ) MeanDZ <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=c(0,0,0,0,0,0), name="ExpMeanDZ" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="ExpCovMZ", means="ExpMeanMZ", dimnames=selVars) objDZ <- mxExpectationNormal( covariance="ExpCovDZ", means="ExpMeanDZ", dimnames=selVars) fitFunction <- mxFitFunctionML() # Combine Groups modelMZ <- mxModel( MZlow, CovMZ, MeanMZ, dataMZ, objMZ, fitFunction, name="MZ") modelDZ <- mxModel( DZlow, CovDZ, MeanDZ, dataDZ, objDZ, fitFunction, name="DZ") minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) SatModel <- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj) # ------------------------------------------------------------------------------------------------------------------------------- # 1a) RUN Saturated Model (Cholesky Decomposition) SatFit <- mxRun(SatModel) (SatSum <- summary(SatFit)) # Generate SatModel Output round(SatFit@output$estimate,4) mxEval(MZ.ExpMeanMZ, SatFit) mxEval(MZ.ExpCovMZ, SatFit) mxEval(DZ.ExpMeanDZ, SatFit) mxEval(DZ.ExpCovDZ, SatFit) # ----------------------------------------------------------------------------------------------------------------------------------- # 1b) Specify Bivariate Saturated Model (Gaussian Decomposition) # We will use this specification to fit a constrained model # ----------------------------------------------------------------------------------------------------------------------------------- svM <-c(0,0,0,0,0,0) # these are start values svSD <-c(1,1,1,1,1,1) svMZ<-c(-.04,.2,-.4,-.1,.2,-.4) svDZ <-c(-.03,-.1,-.1,-.1,.5,-.3) # Matrix & Algebra for expected means and covariances MeanMZ <-mxMatrix( "Full", 1, ntv, free=T, values=svM, labels=c("MZm11", "MZm21","MZm31","MZm12", "MZm22",'Mzm32'), name="ExpMeanMZ" ) MZsd <-mxMatrix( "Diag", ntv, ntv, free=T, values=svSD,labels=c("MZs11", "MZs21","MZs31","MZs12", "MZs22",'Mzs32'), name="sdMZ" ) Cormz <-mxMatrix( "Stand", ntv, ntv, free=T, values=svMZ, labels=c("rPhMZ1","rMZ1","MZxtxt1","MZxtxt2","rMZ2","rPhMZ2"),name="MZCor") CovMZ <-mxAlgebra( expression=sdMZ %*% MZCor %*% t(sdMZ), name="ExpCovMZ" ) MeanDZ <-mxMatrix( "Full", 1, ntv, free=T, values=svM, labels=c("DZm11", "DZm21", "DZm12", "DZm22"), name="ExpMeanDZ" ) DZsd <-mxMatrix( "Diag", ntv, ntv, free=T, values=svSD, labels=c("DZs11", "DZs21", "DZs12", "DZs22"), name="sdDZ" ) Cordz <-mxMatrix( "Stand", ntv, ntv, free=T, values=svDZ, labels=c("rPhDZ1","rDZ1","DZxtxt1","DZxtxt2","rDZ2","rPhDZ2"),name="DZCor") CovDZ <-mxAlgebra( expression=sdDZ %*% DZCor %*% t(sdDZ), name="ExpCovDZ" )
First thing I don't understand is how do we determine the svMZ and svDZ values? How many values do we need to enter? Originally there were 6 for 2 variables, in this case do I have to enter 9 values for 3 variables?
Secondly, how do we write the elements of the matrix in the corMZ and corDZ part?
Thanks!
Printer-friendly version
First of all, I notice that you're treating your dummy-coded variable as though it's continuous, and not as a "threshold trait". Are you sure that's what you want to do?
That's not necessary, and I'm not sure it's even a good idea.
They're start values, so they should just be reasonable guesses at what the corresponding statistics are in your datasets. For instance, you could use the output of
cor(mzData,use="complete")andcor(dzData,use="complete")to guide your choice of start values. Remember that, by default, R puts values into matrices in column-major order.No. The number of unique elements in a kxk correlation matrix is k(k-1)/2. So, for a 6x6 correlation matrix (for 3 traits and 2 twins), you would need 6(5)/2 = 15 values.
Sorry, I don't understand this question...?
Finally, consider replacing these two lines,
, with some syntax that uses
mxFitFunctionMultigroup.Hi Rob!
Thank you for your answers. Well, how can I add it (dummy-coded variable) with the threshold trait. Frankly, I couldn't find how to include a dummy coded variable in the analysis.
Secondly, how do we write the elements of the matrix in the corMZ and corDZ part?
What I meant by this question was how I should write the corDZ and corMZ matrices using the svMZ and svDZ values. I wrote 11 21 31 12 22 32 (for 3 traits and 2 twins) in the MeanMZ line. But I am confused how to write this in CorMZ line. Should it be written as "rPhMZ1","rMZ1","MZxtxt1","MZxtxt2","rMZ2","rPhMZ2","rPhMZ3","rMZ3","MZxtxt3".
Thank you for your time and interest!
Again, the two correlation matrices are 6x6, so they have (at most) 15 unique elements. However, due to the twin-related, intraclass structure of your data, there are only nine unique elements in each, since some elements are constrained equal to each other. You need to define the MZ correlation matrix as something like this,
, and correspondingly for the DZ correlation matrix. I didn't write the script you're modifying, so I don't know for sure, but I can guess at what the labels are supposed to mean. "rMZ" probably means "MZ correlation for a given phenotype", "rPhMZ" probably means "correlation between two phenotypes as estimated in MZ twins", and "xtxt" almost certainly means "cross-twin, cross-trait".
This script provides an example of a joint analysis of one continuous and one threshold trait.
Thanks! I tried the script you mentioned in the comment, and was able to get all the parameters. However, when looking at the confidence interval, some values look NA. How can I eliminate this?
lbound estimate ubound
twoACEvj.US[1,10] -1.0535 -0.1931 0.6014
twoACEvj.US[1,13] -0.2559 0.3224 0.8897
twoACEvj.US[1,16] 0.5613 0.8707 1.2305
twoACEvj.US[2,10] NA -3.2858 NA
twoACEvj.US[2,13] NA 0.4030 NA
twoACEvj.US[2,16] NA 3.8829 NA
twoACEvj.US[2,11] -0.0382 0.3953 0.8303
twoACEvj.US[2,14] -0.0805 0.3184 0.6339
twoACEvj.US[2,17] 0.1786 0.2863 0.4704
Also, I want to report it like as in the image [Ace.h2 [1,1]....], so I tried the script I posted as the first comment. How can I do this with the script you shared?
Thanks!
OpenMx is reporting some of those confidence limits as
NAbecause it thinks they might be wrong. You can see the masked confidence limits, and why OpenMx thinks they might be wrong, if you pass argumentverbose=TRUEtosummary(), e.g.summary(myFittedModel,verbose=TRUE). For a 95% confidence interval, the fit at a correct confidence limit should be 3.841 higher than the fit at the maximum-likelihood solution; with that in mind, you can judge how wrong some of the masked confidence limits are.Sorry, could you rephrase the question?
Also, I want to report it like as in the image [Ace.h2 [1,1]....], so I tried the script I posted as the first comment. How can I do this with the script you shared?
I don't understand why the confidence intervals are not reported sequentially across the 3 variables. I think it stems from this part:
And it is reported as:
It confused me to interpret this. Instead I want to report one by one such as twoACEvj.US[1,1], twoACEvj.US[1,2]....Also, what can I do to examine the confidence intervals of h2, c2,e2,Rph,Ra,Re.
Thank you for your answers and your time!
Some advice: one of the guiding principles in OpenMx's design is WYSIWID--"What You Say Is What It Does". If you want confidence intervals for every element of the MxAlgebra named "US", then just reference it in your call to
mxCI(), e.g.ciACE <- mxCI("US"). Of course, if you want confidence intervals for other matrices or algebras, you'll need to include them in the reference as well, e.g.,. However, you first need to put MxAlgebras or MxMatrices with those additional names into your MxModel (which your current script does not do).
In addition to my previous post, I have also attached my script: