# Adding covariates to Cholesky script

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!

## please provide a reproducable example

Log in or register to post comments

In reply to please provide a reproducable example by jpritikin

## In this case, I think the

Log in or register to post comments

## some advice

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 put `meanG`, `Male`, `Age`, `Edu`, `bMale`, `bAge`, and `bEdu` into the MZ and DZ MxModels somehow. As you discovered, you can't put them into the `pars` list if you're going to put `pars` 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`, and `expMean`

directlyinto the MZ and DZ MxModels, or put them into `pars` and don't put `pars` 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:

$labels

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]

[1,] "data.Age1" "data.Age2" "data.Age1" "data.Age2" "data.Age1" "data.Age2" "data.Age1" "data.Age2" "data.Age1" "data.Age2"

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

# create matrices to store observed values for each covariate

Male <- mxMatrix(type="Full", nrow=1, ncol=2, 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=2, free=FALSE,

labels=c("data.Age1","data.Age2"), name="Age")

Edu <- mxMatrix(type="Full", nrow=1, ncol=2, free=FALSE,

labels=c("data.Education1","data.Education2"), name="Edu")

# create matrices to store coefficients

bMale <- mxMatrix(type="Full", nrow=1, ncol=5, free=TRUE,

values= .01, label=paste("betaMale",1:5,sep=""), name="bMale")

bAge <- mxMatrix(type="Full", nrow=1, ncol=5, free=TRUE,

values= .01, label=paste("betaAge",1:5,sep=""), name="bAge")

bEdu <- mxMatrix(type="Full", nrow=1, ncol=5, free=TRUE,

values= .01, label=paste("betaEdu",1:5,sep=""), name="bEdu")

`## Add Covariates to the mean`

## Note the use of the Kronecker-product operator, %x% :

expMean <- mxAlgebra(meanG + (Male%x%bMale) + (Age%x%bAge) + (Edu%x%bEdu), name="expMean")

Log in or register to post comments

In reply to some advice by AdminRobK

## Thank you so much for your

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

`> Age$labels`

[,1] [,2]

[1,] "data.Age1" "data.Age2"

But I assume this pattern is just repeating x5 for the phenotypes. Could I just read in my variables using:

`> Vars <- c('Z_Wellbeing','Alpha','Beta','Theta','Delta') %>%`

+ rep(rep(1), each = 2)

> selVars <- paste(Vars,c(rep(c(1,2),nv)),sep="")

> selVars

[1] "Z_Wellbeing1" "Z_Wellbeing2" "Alpha1" "Alpha2" "Beta1" "Beta2" "Theta1"

[8] "Theta2" "Delta1" "Delta2"

So that each variable is repeated for each twin? I've implemented this for now.

Thanks again,

Miranda

Log in or register to post comments

In reply to Thank you so much for your by Mirch

## I should have run this all

`> # Print the rG, rC, rE matrices (above, in the summary output you will find them with 95% CI)`

> CholAceFit$Acor

mxAlgebra 'Acor'

$formula: iSDa %&% A

$result:

[,1] [,2] [,3] [,4] [,5]

[1,] 1 1 1 1 1

[2,] 1 1 1 1 1

[3,] 1 1 1 1 1

[4,] 1 1 1 1 1

[5,] 1 1 1 1 1

dimnames: NULL

The matrix I get with the original code (errors fixed) is:

`mxAlgebra 'Acor'`

$formula: iSDa %&% A

$result:

[,1] [,2] [,3] [,4] [,5]

[1,] 1.00000000 -0.02157329 -0.1882037 -0.2731677 -0.6819023

[2,] -0.02157329 1.00000000 0.6261026 0.6672235 0.7141160

[3,] -0.18820375 0.62610265 1.0000000 0.5182056 0.6364864

[4,] -0.27316768 0.66722346 0.5182056 1.0000000 0.7815057

[5,] -0.68190235 0.71411600 0.6364864 0.7815057 1.0000000

dimnames: NULL

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

Log in or register to post comments

In reply to I should have run this all by Mirch

## I think you're confused

If you're using the syntax I provided you, then `Male`, `Age`, and `Edu` 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,

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

only works if the phenotypes are ordered as they were in your opening post. Trust me, the additional effort entailed isn't worth it.

Log in or register to post comments

In reply to I think you're confused by AdminRobK

## Okay, great! I've gone back

Thanks for the clarification!

Log in or register to post comments

## Thanks!

Log in or register to post comments

In reply to Thanks! by DryLabRebel

## good to hear

Log in or register to post comments