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!