Question about covariates/definition variables

14 posts / 0 new

Last post

Sat, 08/12/2023 - 10:52

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Question about covariates/definition variables

Hi all,

I have a question about covariates/definition variables. I have noticed that in most of the scripts (Boulder workshop) covariates are included as definition variables (also in GxE models). I just wondered if covariates could be included as manifest(measured) variables and a regression path with the studied phenotype. I know both methods should be equivalent but this would allow to include missing values. Is there any reason/advantage for including covariates always as definition variables?

I have read the OpenMx user guide and the examples about how to fit a regression but I am not sure if this is related about the way to fit the models (path vs matrix). I guess the matrix specification is easier as compared to the path specification.

Thank you so much in advance!

With all good wishes

Fri, 08/18/2023 - 14:10

mhunter

Offline

Joined: 07/31/2009 - 15:26

Complete exogeneity

When you use a definition variable as a covariate you are assuming complete exogeneity for that variable. That is, you make no distributional assumptions about the definition variable; the definition variable is measured without error; missing data is not allowed for that variable.

When you use a manifest variable as a covariate you are assuming endogeneity. That is, you assume the covariate is part of multivariate normal distribution with all the other manifest variables; the covariate might be measured with error; missing data is handled via full information maximum likelihood.

If you use the same variable as both a definition variable and a manifest variable, you are apply the exogeneity assumptions for the definition variable uses and the endogeneity assumptions for the manifest variable uses.

As a side note, many of the regression examples for OpenMx make the endogeneity assumption; whereas classic ordinary least squares regression makes the exogeneity assumption for all predictors. When the data are multivariate normal, this doesn't matter much, but there is a technical difference between these.

Sat, 08/19/2023 - 06:40

(Reply to #2) #3

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Hi Mike,

This is really helpful! Thank you so much for your great response.

I am trying to add a covariate (manifest) to an ACE model, but I have not been able to do it. Are you aware of any ACE script with a covariate as a manifest variable?

I have tried by adding this part from the OpenMx User guide:

matrA <- mxMatrix( type="Full", nrow=2, ncol=2,
                   free=c(F,F,T,F), values=c(0,0,1,0),
                   labels=c(NA,NA,"beta1",NA), byrow=TRUE, name="A" )
matrS <- mxMatrix( type="Symm", nrow=2, ncol=2,
                   free=c(T,F,F,T), values=c(1,0,0,1),
                   labels=c("varx",NA,NA,"residual"), byrow=TRUE, name="S" )
matrF <- mxMatrix( type="Iden", nrow=2, ncol=2, name="F" )
matrM <- mxMatrix( type="Full", nrow=1, ncol=2,
                   free=c(T,T), values=c(0,0),
                   labels=c("meanx","beta0"), name="M")
expRAM <- mxExpectationRAM("A","S","F","M", dimnames=c("X","Y"))

and adding b1 to the ExpMean

but I am not able to make it work. Should I remove the residual variance? I am also not sure if I should use mxExpectationRAM(). Any guidance would be greatly appreciated.

Thank you so much in advance.
With all good wishes!

Mon, 08/21/2023 - 06:11

(Reply to #3) #4

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Hi again,

I have been working on this and I have fitted these two models:

First Model

# ----------------------------------------------------------------------------------------------------------------------
# Select Variables for Analysis
vars <- c('PHQ8_181','PHQ8_182','Age') # list of variables names
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- vars
# Select Data for Analysis
mzData    <- subset(twinData, Zyg==1|Zyg==3, c(selVars))
dzData    <- subset(twinData, Zyg==2|Zyg==4|Zyg==5, c(selVars))
# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
# Set Starting Values
svMe <- 20 # start value for means
svPa <- .2 # start value for path coefficient
svPe <- .5 # start value for path coefficient for e
# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=3, free=TRUE, values=svMe, labels=c("meanPHQ", "meanPHQ", "MeanAge"), name="meanG" )
# Create Matrices for Variance Components
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" )
covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VE11", name="VE" )
 
covAge <- mxMatrix( type="Symm", nrow=1, ncol=1, free=TRUE, values=5, label="VAge11", name="VAge" )
covB1 <- mxMatrix( type="Symm", nrow=1, ncol=1, free=TRUE, values=5, label="VB11", name="VB" )
 
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= VA+VC+VE, name="V" )
covP2 <- mxAlgebra( expression= V+(VAge*(VB^2)), name="V2" )
 
covMZ <- mxAlgebra( expression= VA+VC, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%VA+ VC, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V2, cMZ,0), cbind(t(cMZ), V2,0),cbind(0,0,VAge)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V2, cDZ,0), cbind(t(cDZ), V2,0),cbind(0,0,VAge)), name="expCovDZ" )
# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list( meanG, covA, covC, covE, covP,covP2,covAge,covB1 )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Unstandardized and Standardized Variance Components
rowUS <- rep('US',nv)
colUS <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list(rowUS,colUS) )
# Create Confidence Interval Objects
ciACE <- mxCI( "US[1,1:3]" )
# Build Model with Confidence Intervals
modelACE <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE )
# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
sumACE <- summary( fitACE )

Second model

# ----------------------------------------------------------------------------------------------------------------------
# Select Variables for Analysis
vars <- c('PHQ8_181','PHQ8_182','Age') # list of variables names
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- vars
# Select Data for Analysis
mzData    <- subset(twinData, Zyg==1|Zyg==3, c(selVars))
dzData    <- subset(twinData, Zyg==2|Zyg==4|Zyg==5, c(selVars))
# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
# Set Starting Values
svMe <- 20 # start value for means
svPa <- .2 # start value for path coefficient
svPe <- .5 # start value for path coefficient for e
# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=3, free=TRUE, values=svMe, labels=c("meanPHQ", "meanPHQ", "MeanAge"), name="meanG" )
expMean <- mxAlgebra( expression= meanG +(VAge*(VB^2)) , name="expMeanG" )
# Create Matrices for Variance Components
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" )
covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VE11", name="VE" )
 
covAge <- mxMatrix( type="Symm", nrow=1, ncol=1, free=TRUE, values=5, label="VAge11", name="VAge" )
covB1 <- mxMatrix( type="Symm", nrow=1, ncol=1, free=TRUE, values=5, label="VB11", name="VB" )
 
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= VA+VC+VE, name="V" )
covP2 <- mxAlgebra( expression= V+(VAge*(VB^2)), name="V2" )
 
covMZ <- mxAlgebra( expression= VA+VC, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%VA+ VC, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ,0), cbind(t(cMZ), V,0),cbind(0,0,VAge)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ,0), cbind(t(cDZ), V,0),cbind(0,0,VAge)), name="expCovDZ" )
# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list( meanG,expMean, covA, covC, covE, covP,covP2,covAge,covB1 )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Unstandardized and Standardized Variance Components
rowUS <- rep('US',nv)
colUS <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list(rowUS,colUS) )
# Create Confidence Interval Objects
ciACE <- mxCI( "US[1,1:3]" )
# Build Model with Confidence Intervals
modelACE <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE )
# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
sumACE <- summary( fitACE )

I guess the second one is the correct one, but I am not sure. Is the second one first subtracting all the variance due to Age? And then partitioning the remaining variance into A, C and E?

Thank you so much!

Fri, 09/08/2023 - 11:39

(Reply to #4) #5

AdminNeale

Offline

Joined: 03/01/2013 - 14:09

Model specification issue

There are more direct and easier ways to do what you do - see response from Prof. Maes. I note some issues with the specification that you have been trying to use.

expCovMZ <- mxAlgebra( expression= rbind( cbind(V2, cMZ,0), cbind(t(cMZ), V2,0),cbind(0,0,VAge)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V2, cDZ,0), cbind(t(cDZ), V2,0),cbind(0,0,VAge)), name="expCovDZ" )

This seems to assume that the covariance between age and the phenotypes of the twins should be zero. But if Age is causing some variation in the phenotype, their covariance will typically not be zero. The covariance between age and phenotype should be betaAge * var(Age) - from rules of path analysis, where age has variance parameter var(Age) and the regression of the phenotype on Age is betaAge. The added variation of betaAge^2 times var(Age) affects both twin 1 and twin 2's phenotypic variances and their covariance. This is the case for both MZ and DZ pairs - which is why variance due to Age mimics that due to C. Indeed, adding the variance due to age to C would be one way to implement the model.

Tue, 09/05/2023 - 12:40

(Reply to #5) #6

marijal

Offline

Joined: 05/12/2023 - 06:59

Question on manifest variable as a covariate

Hello,

this comment makes some points to me more clear and I have a question on that. I do specify a single-factor model with 5 indicators and a definition variable. I would like to regress the latent factor on a manifest variable (which is no definition variable). My model is build by specifying different matrices and matrix algebra equations (mxAlgebra) before I do run the model with ML (mxRun). My definition variable is "Age" and I do specify the matrix this way:

matV1 <- mxMatrix(type="Full", nrow=1, ncol=1,
free=FALSE,
labels="data.Age",
name = "Age")

However, I am not sure, how I could specify a simple manifest variable as a covariate (e.g., well-being) which predicts the latent factor in a separate matrix.

Any professional hint would be appreciated.

Wed, 09/06/2023 - 10:02

(Reply to #6) #7

marijal

Offline

Joined: 05/12/2023 - 06:59

My model does look like this

My model does look like this (attached file) where Age is a definition variable with effects on factor variances and means, respectivly. I have defined them this way:

matV1 <- mxMatrix(type="Full", nrow=1, ncol=1,     # Modeling Age as a definition variable
                  free=TRUE, 
                  labels="data.Age", 
                  name = "Age")
 
matH1 <- mxMatrix(type ="Symm", nrow = 1, ncol = 1,  # matH1: contain the direct effects of Age on the common-factor variances and correlation (denoted by h)
                  free = T,  
                  values =0,
                  name ="matH1")
 
matG1 <- mxMatrix(type="Full", nrow=1, ncol=1,     # matG1: contain the direct effects of Age on the common-factor means (represented by g); these direct effects are fixed to zero in the configural model
                  free=T, 
                  values=0,
                  name="matG1")
 
matA <- mxAlgebra(expression=matA0+matG1*Age+matG2*Age2,      # matA: linear moderation function of the common factor means
                  name="matA")

I would like to integrate positive affect as a covariate and I am doint it this way:

matPA <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = FALSE, values = 0, name = "pa")
matH2 <- mxMatrix(type ="Symm", nrow = 1, ncol = 1,  
                  free = T,  
                  values =0,
                  name ="matH2")
 
matG2 <- mxMatrix(type="Full", nrow=1, ncol=1,     
                  free=T, 
                  values=0,
                  name="matG2")
 
matY_PA <- mxAlgebra(expression = matG3*PA , name = "matY_PA")

However, when I run the model with further matrices (mxRun()), there is no SE for G2 (NA) and I get a value for Age2. I do not understand it and I am not sure, if I am doing something wrong. Can anyone explain me, what I am doing wrong?

Best

File attachments:

model.jpg

Wed, 09/06/2023 - 11:32

(Reply to #7) #8

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

script?

I think you'll need to provide your full script before anyone can answer your questions.

Wed, 09/06/2023 - 12:41

(Reply to #8) #9

marijal

Offline

Joined: 05/12/2023 - 06:59

Sure, I can do this

Sure, I can do this:

#Step 1: packages
library(dplyr)
library(OpenMx)
 
#Step 2: Prepare data
df <- readRDS("data.RDS")
 
## Standardize
df$Age <- as.numeric(scale(df$age))
df$panas <- as.numeric(scale(df$panas))
 
## data needed
df_ER <- df[, c("Age", "ER_ANGM", "ER_IMPU", "ER_CONF", "ER_OPTI", "ER_STRE",  "sex", "panas")]
 
## Create data object
mxdataER <- mxData(observed = df_ER, type = "raw")
 
##Indicate names and number of indicators
manVarsER <- colnames(df_ER[, c("ER_ANGM", "ER_IMPU", "ER_CONF", "ER_OPTI", "ER_STRE")])
nv <- length(manVarsER)
 
#Step 3: Full scalar invariance model (MNLFA)
### MxMatrix objects for indicator intercepts of the configural model
matT0 <- mxMatrix(type="Full", nrow=1, ncol=nv,   # matT0: baseline intercepts tau-0
                  free=TRUE,
                  values=1,
                  name="matT0")
matB1 <- mxMatrix(type="Full", nrow=1, ncol=nv,    # matB1: a full matrix containing direct effects of the background variable Age on the intercepts fixed to zero
                  free=FALSE,
                  values=0,
                  name="matB1")
 
## MxMatrix objects for loadings of the configural model
matL0 <- mxMatrix(type="Full", nrow=nv, ncol=1,    # matL0 full matrix containing the baseline factor loadings Lambda0
                  free=c(c(T, T, T, T, T)),
                  values=c(1, 1, 1, 1, 1),
                  byrow=TRUE,
                  name="matL0")
 
matC1 <- mxMatrix(type="Full", nrow=nv, ncol=1,   # matC1: full matrix  containing direct effects of age on the factor loadings fixed to zero
                  free=FALSE
                  byrow=TRUE,
                  values=0,
                  name="matC1")
 
## Matrices of residual variances
matE0 <- mxMatrix(type="Diag", nrow=nv, ncol=nv,   # matE0: diagonal matrix containing the baseline residual variances theta0
                  free=TRUE,
                  values=1,
                  name="matE0")
matD1 <- mxMatrix(type="Diag", nrow=nv, ncol=nv,   # matD1: diagonal matrix containing the effects of age on residual variances (represented as d)
                  free=TRUE,
                  values=0,
                  name="matD1")
 
# Matrices for common-factor variances (and correlations)
matP0 <- mxMatrix(type ="Symm", nrow=1, ncol=1,      
                  free = FALSE, 
                  values = 1,
                  name = "matP0")
 
matH1 <- mxMatrix(type ="Symm", nrow = 1, ncol = 1,  # matH1: contain the direct effects of Age on the common-factor variances 
                  free = TRUE,  
                  values =0,
                  name ="matH1")
 
matH2 <- mxMatrix(type ="Symm", nrow = 1, ncol = 1,  # matH2: contain the direct effects of panas
                  free = TRUE,  
                  values =0,
                  name ="matH2")
 
# Matrices for common factor means
matA0 <- mxMatrix(type="Full", nrow=1, ncol=1
                  free=FALSE,
                  values=0,
                  name="matA0") ,      # matA0: Matrix containing the baseline common-factor means which are fixed to zero to identify the model
 
matG1 <- mxMatrix(type="Full", nrow=1, ncol=1,    
                  free=TRUE, 
                  values=0,
                  name="matG1")  # matG1: contain the direct effects of Age on the common-factor means (represented by g); these direct effects are fixed to zero in the configural model
 
matG2 <- mxMatrix(type="Full", nrow=1, ncol=1,     
                  free=TRUE, 
                  values=0,
                  name="matG2") # matG2: contain the direct effects of panas on the common-factor means (represented by g); these direct effects are fixed to zero in the configural model
 
# Specification of matrices for the matrix algebra
matV1 <- mxMatrix(type="Full", nrow=1, ncol=1,     # Modeling Age as a definition variable
                  free=FALSE, 
                  labels="data.Age", 
                  name = "Age")
 
matPA <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = FALSE, name = "panas") #panas as covariate
 
matT <- mxAlgebra(expression=matT0+matB1*Age,     
                  name="matT") # matT: linear moderation function for the indicator intercepts of a person i
 
matL <- mxAlgebra(expression=matL0+matC1*Age,
                  name="matL")  # matL: linear moderation function for the factor loadings of a person i
 
 
matE <- mxAlgebra(expression=matE0*exp(matD1*Age), 
                  name="matE") # matE: linear moderation function for the residual variances of a person i
 
matA <- mxAlgebra(expression=matA0+matG1*Age+ matG2*panas,  name="matA") # matA: linear moderation function of the common factor means
 
#(Co)Variance
matVar <- mxAlgebra(expression=(matP0*exp(matH1*Age+matH2*panas)),        
                    name="matVar")
 
matR <- mxAlgebra(expression=(exp(2*(matP0+matH1*Age+matH2*panas))-1)/                        (exp(2*(matP0+matH1*Age+matH2*panas))+1), 
                  name="matR")
 
matIa <- mxMatrix(type="Diag", nrow=1, ncol=1,      
                  free=FALSE,
                  values=1, 
                  name="matIa")
matIb <- mxMatrix(type="Full", nrow=1, ncol=1,      
                  free=FALSE, 
                  values=0,
                  name="matIb")
 
matCov <- mxAlgebra(expression=(matIa*sqrt(matVar))%*%matR%*%(matIa*sqrt(matVar)), 
                    name="matCov")
 
matP <- mxAlgebra(expression=matIa*matVar+matIb*matCov, 
                  name="matP")
matM <- mxAlgebra(expression=matT+t(matL%*%matA),           
                  name="matM") # matM: model implied means
 
matC <- mxAlgebra(expression=matL%*%matP%*%t(matL)+matE,    
                  name="matC") # matC: model implied variances and covariances
 
## Specify expectation and fit function
expF <- mxExpectationNormal(covariance="matC", # expF: definition of the way how the model expectations are calculated
                            means ="matM",
                            dimnames=manVarsER)
 
fitF <- mxFitFunctionML()  
 
modConfig <- mxModel(model="Scalar", 
                     matT, matT0, matB1, 
                     matL, matL0, matC1,  
                     matE, matE0, matD1, 
                     matP, matP0, matH1, matH2, 
                     matA, matA0, matG1, matG2, 
                     matIa, matIb, matV1, matPA,
 
                     matVar, matR, matCov, matM, matC, 
                     expF, fitF, mxdataER)
fitScalar <- mxTryHard(modScalar)
summary(fitScalar)

Everything works, except the panas variable

Wed, 09/06/2023 - 12:44

(Reply to #9) #10

marijal

Offline

Joined: 05/12/2023 - 06:59

Last code lines..

It is modScalar <- ... instead of modConfig <- ...

Tue, 08/22/2023 - 01:55

#11

mhunter

Offline

Joined: 07/31/2009 - 15:26

Unclear

Based on your code it's hard for me to tell what you are trying to do. Consequently, I'm not sure if you're succeeding in your intention or not.

Does this model run? The way you're handling the age covariate in this model is unfamiliar to me.

Tue, 08/22/2023 - 10:45

(Reply to #11) #12

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Dear Mike,

Dear Mike,
Thank you so much for your help with this. I just want to add age as a covariate (manifest variable) but I am struggling to do it.
Here you can see my code with some simulated data. I hope this makes things much easier.

library(OpenMx)
library(psych); library(polycor)
source("miFunctions.R")
mxOption(NULL,"Default optimizer","CSOLNP")
set.seed(1)
cov.matone <- matrix(c(1, .7,.5,
                       .7, 1,.5,
                       .5,.5,1), nrow = 3)
 
dfMZ <- data.frame(MASS::mvrnorm(n = 1e4, 
                                  mu = c(0, 0,0), 
                                  Sigma = cov.matone))
 
cov.mattwo <- matrix(c(1, .35,.5,
                       .35, 1,.5,
                       .5,.5,1), nrow = 3)
 
dfDZ <- data.frame(MASS::mvrnorm(n = 1e4, 
                                 mu = c(0, 0,0), 
                                Sigma = cov.mattwo))
 
cov(dfMZ)
cov(dfDZ)
# ----------------------------------------------------------------------------------------------------------------------
# Select Variables for Analysis
selVars <- c('X1','X2','X3') # list of variables names
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
# Select Data for Analysis
mzData    <- dfMZ
dzData    <- dfDZ
# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
# Set Starting Values
svMe <- 0 # start value for means
svPa <- .3 # start value for path coefficient
svPe <- .5 # start value for path coefficient for e
# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=3, free=TRUE, values=svMe, labels=c("meanPHQ", "meanPHQ", "MeanAge"), name="meanG" )
# Create Matrices for Variance Components
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" )
covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VE11", name="VE" )
 
covAge <- mxMatrix( type="Symm", nrow=1, ncol=1, free=TRUE, values=1, label="VAge11", name="VAge" )
covB1 <- mxMatrix( type="Symm", nrow=1, ncol=1, free=T, values=0.3, label="VB11", name="VB" )
 
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= VA+VC+VE, name="V" )
covP2 <- mxAlgebra( expression= V+(VAge*(VB^2)), name="V2" )
 
covMZ <- mxAlgebra( expression= VA+VC+VAge*(VB^2), name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%VA+ VC+VAge*(VB^2), name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V2, cMZ,VAge*(VB)), cbind(cMZ, V2,VAge*(VB)),cbind(VAge*(VB),VAge*(VB),VAge)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V2, cDZ,VAge*(VB)), cbind(cDZ, V2,VAge*(VB)),cbind(VAge*(VB),VAge*(VB),VAge)), name="expCovDZ" )
# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list( meanG, covA, covC, covE, covP,covP2,covAge,covB1 )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Unstandardized and Standardized Variance Components
rowUS <- rep('US',nv)
colUS <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list(rowUS,colUS) )
# Create Confidence Interval Objects
ciACE <- mxCI( "US[1,1:3]" )
# Build Model with Confidence Intervals
modelACE <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE )
# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
sumACE <- summary( fitACE )

The model runs OK and the estimates for the variance due to age match (25%). However, the estimates for A, C and E are not reasonable (C should be 0?). And in the AE model there is no VE and all is VA.

Thank you so much!

Fri, 09/08/2023 - 21:17

(Reply to #12) #13

Hermine

Offline

Joined: 07/31/2009 - 14:06

I believe the script below

I believe the script below does what you're attempting:

library(OpenMx)
library(psych); library(polycor)
source("miFunctions.R")
mxOption(NULL,"Default optimizer","CSOLNP")
set.seed(1)
cov.matone <- matrix(c(1, .7,.5,
                       .7, 1,.5,
                       .5,.5,1), nrow = 3)
 
dfMZ <- data.frame(MASS::mvrnorm(n = 1e4, 
                                  mu = c(0, 0,0), 
                                  Sigma = cov.matone))
 
cov.mattwo <- matrix(c(1, .35,.5,
                       .35, 1,.5,
                       .5,.5,1), nrow = 3)
 
dfDZ <- data.frame(MASS::mvrnorm(n = 1e4, 
                                 mu = c(0, 0,0), 
                                Sigma = cov.mattwo))
 
cov(dfMZ)
cov(dfDZ)
# ----------------------------------------------------------------------------------------------------------------------
# Select Variables for Analysis
selVars <- c('X1','X2','X3') # list of variables names
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
nc <- 1 # number of covariates (assuming covariate is shared across twins)
ntvc <- ntv+nc
# Select Data for Analysis
mzData    <- dfMZ
dzData    <- dfDZ
# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
# Set Starting Values
svMe <- 0 # start value for means
svPa <- .3 # start value for path coefficient
svPe <- .5 # start value for path coefficient for e
# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntvc, free=TRUE, values=svMe, labels=c("meanPHQ", "meanPHQ", "MeanAge"), name="meanG" )
# Create Matrices for Variance Components
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" )
covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VE11", name="VE" )
 
bAge <- mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=1, label="bAge11", name="bAge" )
covB <- mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=0.3, label="VB11", name="VB" )
covAB <- mxAlgebra( expression= bAge*VB, name="covAB" )
 
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= VA+VC+VE +bAge^2, name="V" )
covMZ <- mxAlgebra( expression= VA+VC +bAge^2, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%VA+VC +bAge^2, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ, covAB), cbind(cMZ, V, covAB), cbind(covAB, covAB, VB)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ, covAB), cbind(cDZ, V, covAB), cbind(covAB, covAB, VB)), name="expCovDZ" )
# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list( meanG, covA, covC, covE, covP, bAge, covB, covAB )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Unstandardized and Standardized Variance Components
rowUS <- rep('US',nv)
colUS <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list(rowUS,colUS) )
# Create Confidence Interval Objects
ciACE <- mxCI( "US[1,1:3]" )
# Build Model with Confidence Intervals
modelACE <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE )
# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
sumACE <- summary( fitACE )

I also attach the corresponding path diagram:

File attachments:

oneACE+B.png

Sun, 09/10/2023 - 07:31

(Reply to #13) #14

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Thank you so much, Hermine.

This is really helpful! I really appreciate it.

With all good wishes!

Main menu

Navigation

You are here

Question about covariates/definition variables

Main menu

User login

Navigation

You are here

Question about covariates/definition variables

Search form