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)
I refer you to my comments in another thread.
Thank you so much Rob! The document and the script are extremely helpful. Now I am a little bit less confused.
I have been working on it and I have fitted the full bivariate model with 17parameters. But I still have some doubts.
1)I don’t know if I have to tried to fit the extended univariate moderation model. I think in my data set M and T are correlated but I don’t know how I can check it.
In your script you checked it with the true values for "a0m","c0m","e0m…." but how can I get the true values with my data?
2) Also, I don’t know how can I get the results from my model, I mean I want to know a c e for SL=1 and for SL=2. Have I to create something like in that? But what values should I use?:
3)I am trying to learn it fitting several models but now I am focused in one of them where the moderation variable is dichotomous. This model (the script below) is ok for that? If the trait is dichotomous or ordinal I have to use a threshold model?
4) the order of the variables in the data set should be moderation variable 1, trait 1, moderation variable 2, trait 2, right?
Thank you so much Rob and sorry for the inconveniences I am not very good with that but I need it and I like it!
Here the full script adapted for my data:
First off, I can't take credit (or blame) for the example script, since I didn't write it--Conor V. Dolan did. Secondly, you may find Conor's slides from Boulder 2016 useful. An important take-home message is that you should always use the "full bivariate" model (in which the moderator is an endogenous random regressor), unless (1) the moderator correlates 1.0 between twins, (2) the moderator correlates 0 between twins, or (3) the covariance between the moderator and trait is equally due to A, C, and E. In practice, cases (2) and (3) are exceedingly rare, and case (3) can really only be checked by model-fitting anyhow.
Once you have parameter estimates from your selected moderation model, you can plot the values of standardized and raw A, C, and E variance components at different values of the moderator. Take a look at Figure 2 in my relevant publication. The exact code used to create those plots is:
The scale of your moderator and trait will be different, the symbol for your MxModel will be different, the labels for your free parameters will be different, and the function of the free parameters that maps x to y will be different, but the essential idea is there.
I haven't looked over your script yet., but if you're fitting a moderation model in which the moderator itself is biometrically decomposed, then treating an ordinal moderator as such (i.e., using thresholds) is a good idea. However, you might also be able to side-step all this Purcell-style moderation modeling with savvy use of multiple groups. You'd have the following groups: MZ twins, concordant for score of 1 on the moderator; MZ twins, concordant 0; MZ twins, discordant; DZ twins, concordant 1; DZ twins, concordant 0; DZ twins, discordant. I'm not entirely sure of the details offhand, but in principle this should work.
Yes, that works.
Thank you so much!!! I am working on it.
This reminds me...two miscellaneous points.
First, replace
with
. They are nearly equivalent, but the important thing is that an MxFitFunctionMultigroup is compatible with certain utility functions in OpenMx, whereas an algebra fitfunction isn't.
Second, does your script actually run? I had thought that it would be necessary to enter the moderator, 'SL', into the dataset twice, with one copy used as an endogenous variable, and the other as a "definition" variable. But I could be wrong about that. It looks like you're using the same column as an endogenous variable AND as a definition variable, right?