Hi there, I am a new PhD student working on a project with twins for the first time. I've been provided with some older scripts from past students and some example code from a twin analysis workshop in Queensland that I haven't been able to attend yet (which I believe was run by Sarah Medland).
The code from the workshop is for a Cholesky model without covariates. I would like to regress out sex, age, and education from the mean, but otherwise use the same script. I am having trouble incorporating this and can't see what I'm doing wrong based on the other scripts I have that do have covariates (but are not in the same format as what I would like).
This is my edited code from the workshop:
require(OpenMx) require(psych) library(tidyverse) # Read in data file wellbeing <- read.csv("CleanData/FINAL_EC_Cent_twins.csv", header = TRUE) glimpse(wellbeing) #summary(wellbeing) #describe(wellbeing) nv <- 5 # number of variables ntv <- 2*nv # number of variables*max family size # identify variables Vars <- c('Z_Wellbeing','Alpha','Beta','Theta','Delta') selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") # identify covariates # Covs <- c('male' 'Age','Education','Beta','Theta','Delta') # use if you want to use same method as for Vars selCovs <- c('male', 'Age1', 'Education1', 'Age2', 'Education2') #allVars <- c(selVars, selCovs) # list with variabels and covariates # This yields a list of the 8 variables for analysis # Satisfaction1 Positivity1 Wellbeing1 Mastery1 Satisfaction2 Positivity2 Wellbeing2 Mastery2 mzData <- subset(wellbeing, mz==1, c(selVars, selCovs)) dzData <- subset(wellbeing, mz==0, c(selVars, selCovs)) # Create start values Stmean <-colMeans(mzData[,1:nv],na.rm=TRUE) # Create Labels for Lower Triangular Matrices aLabs <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") cLabs <- paste("c", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") eLabs <- paste("e", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") # Create Labels for Column and Diagonal Matrices mLabs <- paste("m",1:nv,sep="") # Specify Cholesky ACE Decomposition # ------------------------------------------------------------------------------------------------- # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.4, labels=aLabs, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.1, labels=cLabs, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.2, labels=eLabs, name="e" ) # Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+C+E, name="V" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebra to compute Correlations invSDa<- mxAlgebra( expression=solve(sqrt(I*A)), name="iSDa") invSDc<- mxAlgebra( expression=solve(sqrt(I*C)), name="iSDc") invSDe<- mxAlgebra( expression=solve(sqrt(I*E)), name="iSDe") Rph <- mxAlgebra( expression=iSD %&% V , name="Phcor") Rg <- mxAlgebra( expression=iSDa %&% A , name="Acor") Rc <- mxAlgebra( expression=iSDc %&% C , name="Ccor") Re <- mxAlgebra( expression=iSDe %&% E , name="Ecor") # Algebras to compute Standardized Variance Components rowVars <- rep('v',nv) colVars <- rep(c('A','C','E','h2','c2','e2'),each=nv) estVars <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="est", dimnames=list(rowVars,colVars)) # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins ## ---------------------------------------------------------------------- # MeanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=Stmean, labels=c(mLabs,mLabs), name="expMean" ) # Regress covariates out of the mean # find intercept meanG <- mxMatrix(type="Full", nrow=1, ncol=ntv, free=TRUE, values=Stmean, labels=c(mLabs,mLabs), name="meanG" ) # create matrices to store observed values for each covariate Male <- mxMatrix(type="Full", nrow=1, ncol=ntv, free=FALSE, labels=c("data.male","data.male"), name="Male") # note I only have 1 column for male, so no need to specify 'male1' and 'male2' Age <- mxMatrix(type="Full", nrow=1, ncol=ntv, free=FALSE, labels=c("data.Age1","data.Age2"), name="Age") Edu <- mxMatrix(type="Full", nrow=1, ncol=ntv, free=FALSE, labels=c("data.Education1","data.Education1"), name="Edu") # create matrices to store coefficients bMale <- mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="betaMale", name="bMale") bAge <- mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="betaAge", name="bAge") bEdu <- mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="betaEdu", name="bEdu") # regression algebra for each covariate # regMale <- mxAlgebra(bMale%*%Male, name="MaleR") # regAge <- mxAlgebra(bAge%*%Age, name="AgeR") # regEdu <- mxAlgebra(bEdu%*%Edu, name="EduR") ## Add Covariates to the mean expMean <- mxAlgebra(meanG + (bMale%*%Male) + (bAge%*%Age) + (bEdu%*%Edu), name="expMean") covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) # 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="expMean", dimnames=selVars) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars) fit <- mxFitFunctionML() # Combine Groups pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, invSDa, invSDc, invSDe, Rph, Rg, Rc, Re, estVars ) modelMZ <- mxModel( pars, covMZ, expMean, dataMZ, objMZ, fit, name="MZ" ) modelDZ <- mxModel( pars, covDZ, expMean, dataDZ, objDZ, fit, name="DZ" ) fitmg <- mxFitFunctionMultigroup(c('MZ', 'DZ')) ## the below lines were updated using mxFitFunctionMultigroup() ## minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) ## obj <- mxFitFunctionAlgebra( "m2LL" ) conf1 <- mxCI (c ('MZ.est[1,13]', 'MZ.est[2,14]', 'MZ.est[3,15]', 'MZ.est[4,16]') ) # h2 conf2 <- mxCI (c ('MZ.est[1,17]', 'MZ.est[2,18]', 'MZ.est[3,19]', 'MZ.est[4,20]') ) # c2 conf3 <- mxCI (c ('MZ.est[1,21]', 'MZ.est[2,22]', 'MZ.est[3,23]', 'MZ.est[4,24]') ) # e2 conf4 <- mxCI (c ('MZ.Acor[2,1]', 'MZ.Acor[3,1]', 'MZ.Acor[4,1]', 'MZ.Acor[3,2]', 'MZ.Acor[4,2]', 'MZ.Acor[4,3]')) #Rg conf5 <- mxCI (c ('MZ.Ccor[2,1]', 'MZ.Ccor[3,1]', 'MZ.Ccor[4,1]', 'MZ.Ccor[3,2]', 'MZ.Ccor[4,2]', 'MZ.Ccor[4,3]')) #Rc conf6 <- mxCI (c ('MZ.Ecor[2,1]', 'MZ.Ecor[3,1]', 'MZ.Ecor[4,1]', 'MZ.Ecor[3,2]', 'MZ.Ecor[4,2]', 'MZ.Ecor[4,3]')) #Re CholAceModel<- mxModel( "CholACE", pars, modelMZ, modelDZ, fitmg, conf1, conf2, conf3, conf4, conf5, conf6) # ---------------------------------------------------------------------------------------------------------------------------------- # (model 2) RUN Cholesky Decomposition ACE Model CholAceFit <- mxRun(CholAceModel, intervals=F) CholAceSumm <- summary(CholAceFit) CholAceSumm
When I run this, the error I get is:
Error: Unknown expected means name 'expMean' detected in the expectation function of model 'MZ'
I've also tried this:
# Combine Groups pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, invSDa, invSDc, invSDe, Rph, Rg, Rc, Re, estVars, bMale, bEdu, bAge, Male, Age, Edu ) modelMZ <- mxModel( pars, covMZ, meanG, expMean, dataMZ, objMZ, fit, name="MZ" ) modelDZ <- mxModel( pars, covDZ, meanG, expMean, dataDZ, objDZ, fit, name="DZ" ) fitmg <- mxFitFunctionMultigroup(c('MZ', 'DZ'))
Which results in the error code:
Error: A definition variable 'data.male' has been declared in model 'CholACE' that does not contain a data set
When I check the names of my data I get:
> names(mzData) [1] "Z_Wellbeing1" "Alpha1" "Beta1" "Theta1" "Delta1" "Z_Wellbeing2" "Alpha2" [8] "Beta2" "Theta2" "Delta2" "male" "Age1" "Education1" "Age2" [15] "Education2" > names(dzData) [1] "Z_Wellbeing1" "Alpha1" "Beta1" "Theta1" "Delta1" "Z_Wellbeing2" "Alpha2" [8] "Beta2" "Theta2" "Delta2" "male" "Age1" "Education1" "Age2" [15] "Education2"
I spent hours playing around with it trying to get it to run on Friday and everything I've tried has had an error so any help with the code would be appreciated. This is my first time posting on any code-related forum so please let me know if I need to supply any more details or data.
Thank you!
You'll need to attach a complete example that we can run. Include the R script and a few representative rows of data that will allow us to reproduce the error.
In this case, I think the issue is diagnosable solely from the information in the OP.
Hmm. That's not a very helpful error message, is it? Anyhow, it looks like the problem is that none of
expMean
's dependencies is going into the MxModel object anyplace. You need to putmeanG
,Male
,Age
,Edu
,bMale
,bAge
, andbEdu
into the MZ and DZ MxModels somehow. As you discovered, you can't put them into thepars
list if you're going to putpars
into the container MxModel "CholACE", because the container MxModel contains no dataset of its own.So, either put
meanG
,Male
,Age
,Edu
,bMale
,bAge
,bEdu
, andexpMean
directly into the MZ and DZ MxModels, or put them intopars
and don't putpars
into the container MxModel.You will probably also want to re-think how you're setting up the regression onto age, sex, and education level. For instance, you have the phenotypes ordered so that twin #1's block of 5 is followed by twin #2's block of 5, but the labels in your
Age
MxMatrix aren't consistent with that:Furthermore, you only use "data.Education1", but not "data.Education2", as a label for
Edu
(is that deliberate?). Finally, each of your matrices of regression slopes is 1x1, and is being post-multiplied by a 1x10 matrix of definition variables. Do you really want the regression slopes w/r/t a given covariate to be constrained equal across all 5 phenotypes? Perhaps you want to try something like this instead?:Thank you so much for your help, Rob.
You are correct, the double "data.Education1" was a mistake (I've swapped one for "data.Education2"), and I do want my covariates to vary across phenotypes. I've implemented those changes and moved the covariates matrices directly into the MZ and DZ models (and am now using the Kronecker-product as you pointed out).
It all runs now, which is excellent. The only thing I'm not sure about is how to fix the labels? I can't seem to replicate your list of age variables (with 10 values) as I only get
But I assume this pattern is just repeating x5 for the phenotypes. Could I just read in my variables using:
So that each variable is repeated for each twin? I've implemented this for now.
Thanks again,
Miranda
I should have run this all the way through before posting. Changing this is definitely doing something different but I don't think it's correct because I got the following rA matrix:
The matrix I get with the original code (errors fixed) is:
The examples I have of multivariate models with covariates (from previous lab members) have done this both ways - one person had variables as c("A1", "A2", "B1", "B2") and another had the variables as c("A1", "B1", "A2", "B2"). Both use the same format when setting up the regression variables (same as what I have used). I didn't use their code because both used a 'classic' Mx style and written in a way that made it quite difficult to add in my additional variables (and it was hard to tell what was going on!)
If you're using the syntax I provided you, then
Male
,Age
, andEdu
are supposed to only have a two-element vector of labels, as you describe.You don't want to do that, because if you do, then you'll need to completely respecify your model's covariance structure. For instance, this,
only works if the phenotypes are ordered as they were in your opening post. Trust me, the additional effort entailed isn't worth it.
Okay, great! I've gone back to my original running code and will keep the variables in that order.
Thanks for the clarification!
Thanks, this thread was exactly what I needed. I had the same problem!
Glad to hear it. That's what these forums are for!