I looked around the forum, and found referrals to a post that might address this issue, but that page wouldn't load anymore, so I thought I'll post my question here.
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)
Printer-friendly version