Hi everyone,
I am trying to fit a moderation model. I have modified one script but in this script the moderation variable is the same for twin1 and twin 2 (age). I want to fit a moderation model with 1 variable with different values for twin1 and twin2 (like IQ). I have done all the changes but I am not sure if it is correct. Specially this line:
meanG <- mxAlgebra( expression= cbind((int + DR1),(int + DR2)), name="expMeanG" )
Also, I don’t know if the plot is correct since I don’t know how to add the two variables in this line (I think the plot doesn’t change if I use SL1 or SL2 but I am not sure if it is the correct way to do the plot):
SL <- sort(unique(c(mzData$SL2)))
Thank you so much in advance I am stuck with that a little bit of time. Any suggestion and guidance are very welcome
Here the full script:
# Matrices to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Full", nrow=nv, ncol=na, free=TRUE, values=.6, label="a11", name="a" ) pathC <- mxMatrix( type="Full", nrow=nv, ncol=nc, free=TRUE, values=.6, label="c11", name="c" ) pathE <- mxMatrix( type="Full", nrow=nv, ncol=ne, free=TRUE, values=.6, label="e11", name="e" ) # Matrices to store the moderating a, c, and e Path Coefficients modPathA <- mxMatrix( type='Full', nrow=nv, ncol=na, free=TRUE, values=.1, label="aM11", name="aM" ) modPathC <- mxMatrix( type='Full', nrow=nv, ncol=nc, free=TRUE, values=.1, label="cM11", name="cM" ) modPathE <- mxMatrix( type='Full', nrow=nv, ncol=ne, free=TRUE, values=.1, label="eM11", name="eM" ) # Matrices for the expected means vector intercept <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, label="interc", name="int" ) PathM <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0, label=c("m11"), name="m" ) # Matrices to compute the non-moderated A, C, and E Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="nA" ) covC <- mxAlgebra( expression=c %*% t(c), name="nC" ) covE <- mxAlgebra( expression=e %*% t(e), name="nE" ) # Algebra to compute the total variance per twin covP <- mxAlgebra( expression=nA+nC+nE, name="nV" ) # Matrix for the moderator variable mod1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.SL1"), name="Mod1") mod2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.SL2"), name="Mod2") # Matrices to compute the moderated A, C, and E variance components covAmod1 <- mxAlgebra( expression=(a+ Mod1%*%aM) %*% t(a+ Mod1%*%aM), name="A" ) covCmod1 <- mxAlgebra( expression=(c+ Mod1%*%cM) %*% t(c+ Mod1%*%cM), name="C" ) covEmod1 <- mxAlgebra( expression=(e+ Mod1%*%eM) %*% t(e+ Mod1%*%eM), name="E" ) covAmod2 <- mxAlgebra( expression=(a+ Mod2%*%aM) %*% t(a+ Mod2%*%aM), name="AI2" ) covCmod2 <- mxAlgebra( expression=(c+ Mod2%*%cM) %*% t(c+ Mod2%*%cM), name="CI2" ) covEmod2 <- mxAlgebra( expression=(e+ Mod2%*%eM) %*% t(e+ Mod2%*%eM), name="EI2" ) covAmod12 <- mxAlgebra( expression=(a+ Mod1%*%aM) %*% t(a+ Mod2%*%aM), name="AI12" ) covCmod12 <- mxAlgebra( expression=(c+ Mod1%*%cM) %*% t(c+ Mod2%*%cM), name="CI12" ) covAmod21 <- mxAlgebra( expression=(a+ Mod2%*%aM) %*% t(a+ Mod1%*%aM), name="AI21" ) covCmod21 <- mxAlgebra( expression=(c+ Mod2%*%cM) %*% t(c+ Mod1%*%cM), name="CI21" ) # Algebra to compute the total variance per twin covPmod <- mxAlgebra( expression=A+C+E, name="V" ) # Algebra for the expected mean vector and expected variance/covariance matrices 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" ) wmod1 <- mxAlgebra( expression= m %*% Mod1, name="DR1" ) wmod2 <- mxAlgebra( expression= m %*% Mod2, name="DR2" ) meanG <- mxAlgebra( expression= cbind((int + DR1),(int + DR2)), name="expMeanG" ) # 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="expMeanG", dimnames=selVars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanG", dimnames=selVars ) fitMZ <- mxFitFunctionML() fitDZ <- mxFitFunctionML() # Combine Groups parsZ <- list( pathA, pathC, pathE, modPathA, modPathC, modPathE, PathM, intercept) modelMZ <- mxModel( parsZ, meanG, covAmod1, covCmod1, covEmod1,covAmod2,covCmod2,covEmod2,covAmod12,covCmod12,covAmod21,covCmod21, covPmod, mod1, mod2, wmod1,wmod2, covMZ, dataMZ, objMZ, fitMZ, name="MZ" ) modelDZ <- mxModel( parsZ, meanG, covAmod1, covCmod1, covEmod1,covAmod2,covCmod2,covEmod2,covAmod12,covCmod12,covAmod21,covCmod21, covPmod, mod1,mod2, wmod1,wmod2, covDZ, dataDZ, objDZ, fitDZ, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective+ DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) ModAceModel <- mxModel( "ModACE", parsZ, modelMZ, modelDZ, minus2ll, obj ) #=======================================================================# # RUN MODEL # #=======================================================================# # Run Moderation ACE model ModAceFit <- mxRun(ModAceModel) ModAceSumm <- summary(ModAceFit) ModAceSumm round(ModAceFit@output$estimate,4) # Generate Output using Helper Functions tableFitStatistics(ModAceFit) # Generate Table of Parameter Estimates using mxEval pathEstimatesACE <- round(mxEval(cbind(a,c,e,aM,cM,eM), ModAceFit),4) rownames(pathEstimatesACE) <- 'pathEstimates' colnames(pathEstimatesACE) <- c('a','c','e','aM','cM','eM') pathEstimatesACE varComponentsACE <- round(mxEval(cbind(A/V,C/V,E/V), ModAceFit@submodels$MZ),4) #varComponentsACE <- round(mxEval(cbind(A,C,E), ModAceFit@submodels$MZ),4) rownames(varComponentsACE) <- 'varComponents' colnames(varComponentsACE) <- c('%a^2','%c^2','%e^2') varComponentsACE fitGofs(ModAceFit) fitEstS(ModAceFit) #=======================================================================# # PLOT MODERATION # #=======================================================================# a <- pathEstimatesACE[,'a'] c <- pathEstimatesACE[,'c'] e <- pathEstimatesACE[,'e'] aM <- pathEstimatesACE[,'aM'] cM <- pathEstimatesACE[,'cM'] eM <- pathEstimatesACE[,'eM'] SL <- sort(unique(c(mzData$SL2))) Va <- (a + aM*SL)^2 Vc <- (c + cM*SL)^2 Ve <- (e + eM*SL)^2 Vt <- Va + Vc + Ve SVa <- Va/Vt SVc <- Vc/Vt SVe <- Ve/Vt out <- as.matrix(cbind(Va,Vc,Ve,Vt,SVa,SVc,SVe)) rownames(out)=paste('SL=',1:2,sep='') out # par(mfrow=c(2,1)) matplot(SL, out[,1:4], type="l", lty=4:1, col=4:1, xlab="SL", ylab="Variance Component", main="SL Moderation - unstandardized variance components") legend(2.7, 1.2, c("Va","Vc","Ve","Vt"), lty=4:1, col=4:1)