require(OpenMx) require(umx) require(MuMIn) twinData <- read.csv("Moderating_covariances_data.dat.txt", sep="") #=======================================================================# # PREPARE DATA % #=======================================================================# twinData$ses1 = twinData$ses2 = twinData$ses selDVs <- c('iq1','iq2') # phenotype variables to be analyzed selDefs <- c('ses1','ses2') # moderator variable selVars <- c(selDVs, selDefs) ### Move the data into an object for monzygotic and dizygotic twins. mzData <- subset(twinData, zyg == 1, selVars) dzData <- subset(twinData, zyg == 2, selVars) ### Run the umx model for gene by environment interaction. full_moderation <- umxGxE(selDVs = selDVs, selDefs = selDefs, dzData = dzData, mzData = mzData, dropMissingDef = TRUE) #=======================================================================# # * OMG * we are basically done!!! # #=======================================================================# ### Set some parameters to generate figures umxSummaryGxE(full_moderation) umxSummary(full_moderation, location = "topright") umxSummary(full_moderation, separateGraphs = TRUE) # procude a figure plot(full_moderation) # obtain a summary of the results summary(full_moderation) #=======================================================================# # If you publish your GxE twin model, always include the output for # # the full model we ran above! This helpd mets-analyseds, and # # inclusion in meta-analysis means citations for you! # # Citations -> ??? -> profit # <- if this bothers you, you may be conscientious... #=======================================================================# ### We can proceed to drop the moderator of variance components CE_moderation <- umxModify(full_moderation, "am_.*", regex=TRUE, comparison = FALSE) AE_moderation <- umxModify(full_moderation,"cm_.*", regex=TRUE, comparison = FALSE) AC_moderation <- umxModify(full_moderation, "em_.*", regex=TRUE, comparison = FALSE) ### We can drop a second moderator: A_moderation <- umxModify(AE_moderation, "em_.*", regex=TRUE, comparison = FALSE) C_moderation <- umxModify(CE_moderation, "em_.*", regex=TRUE, comparison = FALSE) E_moderation <- umxModify(AE_moderation, "am_.*", regex=TRUE, comparison = FALSE) ### We can drop all variance moderators: means_moderation <- umxModify(E_moderation, "em_.*", regex=TRUE, comparison = FALSE) ### We could also plot the model where only the mean value is moderated umxSummaryGxE(means_moderation) umxSummary(means_moderation, location = "topright") umxSummary(means_moderation, separateGraphs = TRUE) ### We the use umxCompare to sumirize a comparison between two models where A is dropped umxCompare(full_moderation,CE_moderation) umxCompare(AE_moderation,E_moderation) umxCompare(AC_moderation,C_moderation) umxCompare(A_moderation,means_moderation) #=======================================================================# # Inference based on the models above would not be equal # # and chosing the "right" order in whihc to drop parameters, # # or even whewther to drop parameters at al, is a complex question # # the awnser to which is dependent on power, your hypothesis, # # your philospohy about inference... % <- if your bothered again you might want to come work for me, I am not very conscientious, so I could use you... #=======================================================================# #=======================================================================# # One option is to not pick a best model at al, but to use AIC weights # # AIC weights represent the relative probability a model is closest to # # the "true" model, gien all the models you consider # # see: Wagemakers et al https://www.ncbi.nlm.nih.gov/pubmed/15117008 # #=======================================================================# # AIC your.AIC <- AIC(full_moderation,AC_moderation,AE_moderation,CE_moderation,A_moderation,C_moderation,E_moderation,means_moderation) # AICc (AICc is better suited for smaller samples than the AIC) your.AICc <- AICc(full_moderation,AC_moderation,AE_moderation,CE_moderation,A_moderation,C_moderation,E_moderation,means_moderation) # Convert the "raw" AIC into AIC weights aic.weights <- Weights(your.AIC) aiccweights <- Weights(your.AICc) # posterior probability for all models where A is moderated: sum(aic.weights[c(1,2,3,5)]) #=======================================================================# # Try this for yourself! # # GO to this article on GxE: # # journals.plos.org/plosone/article?id=10.1371/journal.pone.0030320 # # Download Table S1, and for 1 of the models, compute the AIC weights # #=======================================================================# # I already computed the AIC weights for the model at age 2 Age2 <- Weights(c(2667.258,2666.556,2675.752,2668.588,2674.027,2666.637,2675.250,2675.001)) Age2 #=======================================================================# # The model the authors identify as the best model at age2, has a 32% # # conditioonal probabillity.. # # two other models have conditional prob of 23% an 31% # # # # Would you want to base public policy on the best model, or would # # you want to consider all three models? # #=======================================================================# #=======================================================================# # Just to be completely fair, I have myself made judgements based on # # model comparison that could be improved if I weree to use the AIC # # Here are the AICs'sfrom 4 models from # # Nivard et al.(2017) Schizophrenia Bulletin # #=======================================================================# NivardEtAl <- Weights(c(-1200.476,-1216.885, -1230.258,-1231.304 )) NivardEtAl # Model 3 and 4 are hard to distinguish, I could have more strongly emphasized this in the paper.... sessionInfo()