--- title: "aceb" author: "Mirko Ruks" date: "2022-09-20" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # Gene-environment interactions First, we load the packages using the `library()` function. These are the same as in the session about the ACE-beta model + the package `umx` which is also a wrapper for `OpenMx`. Please, run the following code chunk and install the packages if necessary. ```{r, message = FALSE} rm(list = ls()) # install.packages("tidyverse") # install.packages("haven") # install.packages("OpenMx") # install.packages("umx") # install.packages("lavaan") # install.packages("psych") # install.packages("devtools") # devtools::install_github("mirkoruks/twinflex") set.seed(1) library(twinflex) library(tidyverse) library(haven) library(lavaan) library(psych) library(OpenMx) library(umx) source("functions.R") ``` Now, we have installed and loaded the packages. Please run the next code chunk. It takes some time as all the variables we need are recoded. We get a data frame back that is called `df_final_wide`. It is a wide-formatted (i.e. one row = one *twin pair*) data set ready to use for SEM-based twin modeling. ```{r, cache = TRUE} source("Session 3/Scripts/recode.R") ``` Run the following code chunk to get some summary statistics of your data. ```{r, message = FALSE} summary(df_final_wide) head(df_final_wide) describe(df_final_wide) ``` ## The first Gene-Environment Interaction Model We want to know whether the effects of the ACE components on cognitive ability vary by age. So, we want to do a GxE analysis. We start with a simple univariate Purcell model. We can do this with `umx` and with `twinflex`. Let's start with `umx` which provides the function `umxGxE()`. As in every interaction model, it is nice if there is a meaningful interpretation for the 0 value of the moderator. So, before we run the model, we standardize the age variable. Here, we have to do some special steps in the data preparation and create a MZ and a DZ data frame. The moderator has two be splitted into twin specific variables. So we split `age` into `age_1` and `age_2`. The `umxGxE()` function has the following arguments: - `selDVs`: the ACE decomposed variable - `sep`: the separator (as in `twinflex()`) - `dzData`: the data set containing the DZ twins - `mzData`: the data set containing the MZ twins - `tryHard`: try out different starting values So, run the code in the following chunk and have a look at the model output. ```{r, echo = TRUE, message=FALSE} df_final_wide <- df_final_wide %>% mutate(age = (age - mean(age, na.rm = TRUE)) / sd(age, na.rm = TRUE)) umx_data <- df_final_wide %>% mutate(age_1 = age) %>% mutate(age_2 = age) df_mz <- umx_data %>% filter(zyg == 1) df_dz <- umx_data %>% filter(zyg == 2) umxGxE_GxE_CA <- umxGxE(selDVs = "math", sep = "_", selDefs = "age", dzData = df_dz, mzData = df_mz, tryHard = "yes", lboundACE = 0.001) summary(umxGxE_GxE_CA) ``` Let's see the model output Here, the parameters with the m before the separator ("_") refer to the interaction effects. The other parameters are the main effects, i.e. they give you the conditional effect for M=0 ```{r} makeOutput(summary(umxGxE_GxE_CA)) ``` Now, let's run the same model with `twinflex`. You can control the moderation via the `moderation` argument. It is comparable to the syntax of `Mplus`. `CA BY age` means that all effects on the variable `CA` are moderated by `age`. In theory, you can include up to five moderators with `twinflex`: `Outcome BY M1 + M2 + M3 + M4`. In practice, these models may fail to converge. Importantly, you can imagine the interaction as a regression of the effect on the moderator. E.g. if you moderate the effect `a11`. Then, you define `a11` as a function of the moderator: `a11 = am11 + bam11*Mod` ```{r} tf_GxE_CA <- twinflex(acevars = "math", sep = "_", zyg = "zyg", data = df_final_wide, moderation = "math BY age") ``` Compare the models estimated with `umx` and `twinflex`. What do you think? ```{r} summary(umxGxE_GxE_CA) summary(tf_GxE_CA) ``` Now, we drop the quadrativ main effect from the model estimated with `umx`: ```{r} umxGxE_GxE_CA_red1_model <- mxModel(umxGxE_GxE_CA, name="Drop quad11") umxGxE_GxE_CA_red1_model <- omxSetParameters(model = umxGxE_GxE_CA_red1_model, labels=c("quad11"), free=FALSE, values=0) umxGxE_GxE_CA_red1_fit <- mxTryHard(umxGxE_GxE_CA_red1_model, exhaustive = TRUE) umxGxE_GxE_CA_red1_cf <- mxCompare(umxGxE_GxE_CA, umxGxE_GxE_CA_red1_fit) umxGxE_GxE_CA_red1_cf ``` and compare the models again: ```{r} summary(tf_GxE_CA) summary(umxGxE_GxE_CA_red1_fit) ``` ## Bivariate GxE Now we want to moderate different ACE paths by different moderators in a bivariate setting. We can use `twinflex()` to do this. We can specify in the `moderation` argument which paths should be moderated by which moderator. There are two "styles": - `Var BY Mod`: Here the effects of the ACE components unique to `Var` are allowed to vary by `Mod` - `Var1 -> Var2 BY Mod`: Here the effect of the ACE components associated with `Var1` on `Var2` are allowed to vary by `Mod` ```{r} tfBivGxE <- twinflex(acevars = c("CA","math"), sep = "_", zyg = "zyg", data = df_final_wide, ModCov = "DV", moderation = c("CA BY female","CA -> math BY age","math BY female"), TryHard = TRUE, lboundACE = FALSE, exh = FALSE) ``` See the model results: ```{r} summary(tfBivGxE) makeOutput(summary(tfBivGxE), rounded = TRUE) ``` Clearly, there are many noisy parameters we want to drop from the model. We start by dropping the parameters of the moderation of all shared environmental effects: ```{r} tfBivGxE_red1_model <- mxModel(tfBivGxE, name="Drop c moderation") tfBivGxE_red1_model <- omxSetParameters(model = tfBivGxE_red1_model, labels=c("bm1c11","bm1c22","bm2c21"), free=FALSE, values=0) tfBivGxE_red1_fit <- mxTryHard(tfBivGxE_red1_model) tfBivGxE_red1_cf <- mxCompare(tfBivGxE, tfBivGxE_red1_fit) tfBivGxE_red1_cf makeOutput(summary(tfBivGxE_red1_fit), rounded = TRUE) ``` Still, we can reduce the model by dropping all moderations of the genetic effects: ```{r} tfBivGxE_red2_model <- mxModel(tfBivGxE_red1_fit, name="Drop a moderation") tfBivGxE_red2_model <- omxSetParameters(model = tfBivGxE_red2_model, labels=c("bm1a11","bm1a22","bm2a21"), free=FALSE, values=0) tfBivGxE_red2_fit <- mxTryHard(tfBivGxE_red2_model) tfBivGxE_red2_cf <- mxCompare(tfBivGxE_red1_fit, tfBivGxE_red2_fit) tfBivGxE_red2_cf makeOutput(summary(tfBivGxE_red2_fit), rounded = TRUE) ``` Finally, we drop the moderation of the `e22` effect: ```{r} tfBivGxE_red3_model <- mxModel(tfBivGxE_red2_fit, name="Drop e22 moderation") tfBivGxE_red3_model <- omxSetParameters(model = tfBivGxE_red3_model, labels=c("bm1e22"), free=FALSE, values=0) tfBivGxE_red3_fit <- mxTryHard(tfBivGxE_red3_model) tfBivGxE_red3_cf <- mxCompare(tfBivGxE_red2_fit, tfBivGxE_red3_fit) tfBivGxE_red3_cf # p > 0.05 -> Ok ! makeOutput(summary(tfBivGxE_red3_fit), rounded = TRUE) ``` ## Moderated ACE-beta model We can also fit a moderated ACE-beta models. We want to see whether the effects linking `CA` and `math` vary by sex: ```{r} tfAceBGxE <- twinflex(acevars = c("CA","math"), sep = "_", zyg = "zyg", type = "aceb", data = df_final_wide, moderation = c("CA -> math BY female"), TryHard = TRUE, exh = FALSE) summary(tfAceBGxE) makeOutput(summary(tfAceBGxE), rounded = TRUE) ``` Now we drop the moderation effects: ```{r} tfAceBGxE_red1_model <- mxModel(tfAceBGxE, name="Drop all moderations") tfAceBGxE_red1_model <- omxSetParameters(model = tfAceBGxE_red1_model, labels=c("bm1b21","bm1a21","bm1c21"), free=FALSE, values=0) tfAceBGxE_red1_fit <- mxTryHard(tfAceBGxE_red1_model) tfAceBGxE_red1_cf <- mxCompare(tfAceBGxE, tfAceBGxE_red1_fit) tfAceBGxE_red1_cf # p > 0.05 -> Ok ! makeOutput(summary(tfAceBGxE_red1_fit), rounded = TRUE) ``` Neither the beta effect nor the genetic/shared environmental confounding varies by sex. ## Two Moderators We also can use multiple moderators. Here, we just add another "covariate" to the regression of the moderated effect on the moderator. See again the case of the effect `a11`: `a11 = ba11 + b1am11*Mod1 + b2am11*Mod2` ```{r} tfAceBGxE2 <- twinflex(acevars = c("CA","math"), sep = "_", zyg = "zyg", type = "aceb", data = df_final_wide, moderation = c("CA -> math BY female + age"), TryHard = TRUE, exh = FALSE) summary(tfAceBGxE2) makeOutput(summary(tfAceBGxE2), rounded = TRUE) ``` ```{r} tfAceBGxE2_red1_model <- mxModel(tfAceBGxE2, name="Drop all sex moderations") tfAceBGxE2_red1_model <- omxSetParameters(model = tfAceBGxE2_red1_model, labels=c("bm1b21","bm1a21","bm1c21"), free=FALSE, values=0) tfAceBGxE2_red1_fit <- mxTryHard(tfAceBGxE2_red1_model) tfAceBGxE2_red1_cf <- mxCompare(tfAceBGxE2, tfAceBGxE2_red1_fit) tfAceBGxE2_red1_cf # p > 0.05 -> Ok ! makeOutput(summary(tfAceBGxE2_red1_fit), rounded = TRUE) ``` ```{r} tfAceBGxE2_red2_model <- mxModel(tfAceBGxE2_red1_fit, name="Drop all age moderations") tfAceBGxE2_red2_model <- omxSetParameters(model = tfAceBGxE2_red2_model, labels=c("bm2b21","bm2a21","bm2c21"), free=FALSE, values=0) tfAceBGxE2_red2_fit <- mxTryHard(tfAceBGxE2_red2_model) tfAceBGxE2_red2_cf <- mxCompare(tfAceBGxE2_red1_fit, tfAceBGxE2_red2_fit) tfAceBGxE2_red2_cf # p > 0.05 -> Ok ! makeOutput(summary(tfAceBGxE2_red2_fit), rounded = TRUE) ``` ## Tasks: Testing the Scarr-Rowe hypothesis One famous hypothesis in behavioral genetics is that the genetic influence on cognitive ability increases with social origin. So, a positive GxSES is assumed. We can test this hypothesis for three outcomes: cognitive ability, math grades and german grades. Run a univariate ACE model for the three outcomes. In a final step, run one or two bivariate models (CA -> math or CA -> german) 1. ... for cognitive ability 2. ... for math grades 3. ... for german grades 4. in a bivariate model