Question about covariates/definition variables
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
Complete exogeneity
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.
Log in or register to post comments
In reply to Complete exogeneity by mhunter
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!
Log in or register to post comments
In reply to Hi Mike, by JuanJMV
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!
Log in or register to post comments
In reply to Hi again, by JuanJMV
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.
Log in or register to post comments
In reply to Complete exogeneity by mhunter
Question on manifest variable as a covariate
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.
Log in or register to post comments
In reply to Complete exogeneity by mhunter
My model does look like this
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
Log in or register to post comments
In reply to My model does look like this by marijal
script?
Log in or register to post comments
In reply to script? by AdminRobK
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
Log in or register to post comments
In reply to Sure, I can do this by marijal
Last code lines..
Log in or register to post comments
Unclear
Does this model run? The way you're handling the age covariate in this model is unfamiliar to me.
Log in or register to post comments
In reply to Unclear by mhunter
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!
Log in or register to post comments
In reply to Dear Mike, by JuanJMV
I believe the script below
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:
Log in or register to post comments
In reply to I believe the script below by Hermine
Thank you so much, Hermine.
This is really helpful! I really appreciate it.
With all good wishes!
Log in or register to post comments