Age as covariate or moderator
I'm running a set of univariate ACE models, on a range of outcomes. I have two age groups (3-4 and 7-8), that I would like to analyse separately, by also combined.
So far, I ran a model that includes age as a continuous covariate, the script for which (a slightly edited version of one that was provided at last year's IBG) is pasted below. My question regarding this script is just: how can I see the influence of age on my model, can I get something like a regression coefficient from this model?
And my second question, if I'm interested in the question of whether my ACE model estimates are different across my age groups, would it make more sense to run a moderated model instead, with age group as a binary rather than continuous variable?
Script:
# ------------------------------------------------------------------------------
# Program: oneACEca.R
# Author: Hermine Maes
# Date: 02 25 2016
#
# Twin Univariate ACE model to estimate means and (co)variances across multiple groups
# Matrix style model - Raw data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|
setwd("H:...")
# Load Libraries & Options
library(OpenMx)
library(psych)
library(foreign)
source("miFunctions.R")
twinDataAll<- read.spss("....sav",use.value.labels=F, max.value.labels=Inf,to.data.frame=TRUE)
twinDataAll<-data.frame(twinDataAll)
head(twinDataAll)
describe(twinDataAll)
twinData<-twinDataAll[c("Zygosity", "AgeTwinYears", "TMCQ_C1", "TMCQ_C2")]
colnames(twinData)<- c("Zygosity", "Age", "TMCQ_C1", "TMCQ_C2")
Selvars = c("TMCQ_C1", "TMCQ_C2")
# Select Variables for Analysis
vars <- 'TMCQ_C' # list of variables names
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
# Select Covariates for Analysis
twinData[,'Age'] <- twinData[,'Age']/1
twinData <- twinData[-which(is.na(twinData$Age)),]
covVars <- 'Age'
head(twinData)
describe(twinData)
## Split file in mono en dizy
mzData <- subset(twinData, Zygosity==1, c(selVars, covVars))
dzData <- subset(twinData, Zygosity==2, c(selVars, covVars))
describe(mzData)
describe(dzData)
# ------------------------------------------------------------------------------
# Set Starting Values
svMe <- 2 # start value for means
svPa <- .4 # start value for path coefficient
svPe <- .8 # start value for path coefficient for e
lbPa <- .0001 # start value for lower bounds
# PREPARE MODEL
# ACE Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE,
labels=c("data.Age"), name="Age" )
pathB <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values=.01, label="b11", name="b" )
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels="xbmi", name="meanG" )
expMean <- mxAlgebra( expression= meanG + cbind(b%*%Age,b%*%Age),
name="expMean" )
# Create Matrices for Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="a11", name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="c11", name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="e11", name="e" )
# Create Algebra for 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" )
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= A+C+E, name="V" )
covMZ <- mxAlgebra( expression= A+C, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(A+C+E, A+C),
cbind(A+C, A+C+E)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(A+C+E, 0.5%x%A+C),
cbind(0.5%x%A+C, A+C+E)), 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="expMean",
dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean",
dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list(pathB, meanG, pathA, pathC, pathE, covA, covC, covE, covP)
defs <- list(defAge)
modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ,
dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ,
dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Variance Components
rowVC <- rep('VC',nv)
colVC <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V),
name="VC", dimnames=list(rowVC,colVC))
mxOption(NULL,"Default optimizer","SLSQP")
# Build Model with Confidence Intervals
modelACE <- mxModel( "oneACEca", pars, modelMZ, modelDZ, multi, estVC )
# ------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxRun( modelACE )
sumACE <- summary( fitACE )
# Print Goodness-of-fit Statistics & Parameter Estimates
fitEsts(fitACE)
# ------------------------------------------------------------------------------
# RUN SUBMODELS
# Run AE model
modelAE <- mxModel( fitACE, name="oneAEca" )
modelAE <- omxSetParameters( modelAE, labels="c11", free=FALSE, values=0 )
fitAE <- mxRun( modelAE, intervals=T )
mxCompare( fitACE, fitAE )
fitEsts(fitAE)
# Run CE model
modelCE <- mxModel( fitACE, name="oneCEca" )
modelCE <- omxSetParameters( modelCE, labels="a11", free=FALSE, values=0 )
fitCE <- mxRun( modelCE, intervals=T )
mxCompare( fitACE, fitCE )
fitEsts(fitCE)
# Run E model
modelE <- mxModel( fitAE, name="oneEca" )
modelE <- omxSetParameters( modelE, labels="a11", free=FALSE, values=0 )
fitE <- mxRun( modelE, intervals=T )
mxCompare( fitAE, fitE )
fitEsts(fitE)
# ------------------------------------------------------------------------------
# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list( fitAE, fitCE, fitE) )
round(rbind(fitACE$VC$result,fitAE$VC$result,fitCE$VC$result,fitE$VC$result),4)