--- title: "aceb" author: "Mirko Ruks" date: "2022-09-20" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # ACE-beta models First, we load the packages using the `library()` function. In the first 6 lines of the code chunk above, you can find the code to install these packages (via `install.packages()` for packages stored in CRAN and `devtools:install_github()` for packages stored in GitHub). If you don't have all packages installed, please uncomment the corresponding lines in the code chunk. We also load some helper functions stored in the file `functions.R`. It is saved in the main directory of the folder you downloaded from sciebo. You can run a code chunk by clicking on the red triangle in the upper right of the code chunk. Alternatively, you can click into the code chunk and press `Ctrl+Alt+C` to run the selected code chunk. Or you click in one line and press `Ctrl+Enter` to run the selected line. Before we run the script: - You have to tell RStudio the location of the default working directory when working in RMarkdown. We want wo use the project directory as the default working directory. You can do so by clicking on "Tools" -> "Global Options..." -> "R Markdown" -> Change "Evaluate chunks in directory" to "Project". - Rename the name of the TwinLife data set in recode.R if necessary (default is `ZA6701_person_wid1_v6-0-0.dta`). Please, run the following code chunk and install the packages if necessary. ```{r, message = FALSE} # install.packages("tidyverse") # install.packages("haven") # install.packages("OpenMx") # install.packages("lavaan") # install.packages("psych") # install.packages("devtools") # devtools::install_github("mirkoruks/twinflex") rm(list = ls()) library(twinflex) library(tidyverse) library(haven) library(lavaan) library(psych) library(OpenMx) source("functions.R") set.seed(1) ``` 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 2/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 ACE-beta Model We want to know whether there is an effect of cognitive ability on math grades when controlling for genetic and shared environmental confounding. Therefore, the ACE-beta model is a good choice. The `twinflex()` function has several arguments: - `acevars`: String vector where you specify the variables you want to decompose into ACE components - `sep`: String vector of length 1 where you specify the separator between the variable name and the suffix reffering two the twins - `zyg`: String vector of length 1. At the moment it has to be `zyg = "zyg"`. So the zygosity must be stored in a column named "zyg" in your data frame. And: The levels need to be 1 = MZ twins and 2 = DZ twins. - `type`: Here, we choose `type = "aceb"` because we want to run a ACE-beta model. Alternatively, we can choose `type = "Chol"` to run a Cholesky-type model. - `data`: Specify the name of your data frame here - `TryHard`: If `TryHard = TRUE`, different starting values are used to estimate the model in order to prevent ending in a local minimum instead of a global minimum. Good choice for more complex models. ```{r, echo = TRUE, message=FALSE} model1 <- twinflex(acevars = c("CA","math"), sep = "_", zyg = "zyg", type = "aceb", data = df_final_wide, TryHard = TRUE) df_final_wide_test <- df_final_wide %>% mutate(CA_test_1 = CA_1) %>% mutate(CA_test_2 = CA_2) model1_test <- twinflex(acevars = c("CA_test","math"), sep = "_", zyg = "zyg", type = "aceb", data = df_final_wide_test, TryHard = TRUE) ``` Now, we want to inspect the model results. There are two functions: - `summary()`: Uses the model fit object as input (i.e. the output of `twinflex()`) providing a good overview, including model fit indices. Drawback: Only Standard Errors, no stars! - `makeOutput()`: Uses the model fit object as input (i.e. the output of `twinflex()`) providing a good overview, including easy to assess inference stats. Drawback: At the moment no model fit indices (e.g. -2LL)! ```{r} summary(model1) makeOutput(summary(model1), rounded = TRUE, n = 2) ``` We see that there are genetic, shared and non-shared environmental effects on both, cognitive ability and math grades. Cognitive ability has the index 1 (as it comes first in the `acevars()` argument) and math grades the index 2. The notation of the effects is always the same: `a/c/e/b+ Index_DepVar + Index_IndepVar`. `a` = Genetic effects, `c` = Shared environmental effects, `e` = Non-shared environmental effects. `mu` is the mean/intercept followed by the index of the variable. Great, so some questions for you. We will discuss them together: 1. What does the parameter `a21` tells us? 2. Which parameter informs us about the effect of the shared environment unique to math grades on math grades? 3. Does cognitive ability affect math grades? In other words, is there a significant phenotypic effect? 4. Is there genetic and shared environmental confounding? If so: What are the parameters? Is there a upward or downward bias of the cross-sectional OLS estimate (to answer this question you just have to look at the effects, no formulas are needed!)? 5. Why is the parameter `e21` missing? ## Calculating the bias Cool. Going back to question 4. Let's calculate the cross-sectional effect based on the model estimates. You can look up the formula in the slides of the presentation. Keep in mind that we assume `e21=0` `` ```{r} a11 <- 0.84 a21 <- 0.32 c11 <- 0.65 c21 <- 0.12 e11 <- 0.81 b21 <- 0.09 varx <- (a11^2 + c11^2 + e11^2) covyx <- (b21*varx)+(a21*a11)+(c21*c11) b_ols <- covyx/varx b_ols ``` ## Model reduction Let's see if we can drop the parameter `c21` from the model. Before we start: What does it mean *theoretically* to drop `c21` from the model? For the model reduction, we need some functions from the `OpenMx` package. We will use the function `mxModel()` to define a new model, the function `omxSetParameters()` to modify the model parameters, the function `mxTryHard()` to fit the new modified model and the function `mxCompare()` to compare the models. The most important function here is `omxSetParameters()` as here we specify what we want to do with the model parameters. These are the arguments: - `model`: Here, you specify the mxModel object (mxFit object works as well) of the model you want to modify - `labels`: Here, you specify the parameters you want to modify - `free`: Specify whether a parameter should be estimated freely (TRUE) or not (FALSE) - `values`: Specify the starting value (if `free = TRUE`) or the value wo which you want to fix the parameter (if `free = FALSE`) Print the object `model1_red1_cf` to get the results of the Likelihood-Ratio test of the full and constrained model. What does it tell? ```{r} model1_red1_model <- mxModel(model1, name="Drop c21") model1_red1_model <- omxSetParameters(model = model1_red1_model, labels=c("c21"), free=FALSE, values=0) model1_red1_fit <- mxTryHard(model1_red1_model, exhaustive = TRUE) model1_red1_cf <- mxCompare(model1, model1_red1_fit) model1_red1_cf # p > 0.05 -> Ok ! makeOutput(summary(model1_red1_fit), rounded = T) summary(model1_red1_fit) ``` ## Calculating the heritability In the bivariate model, we have to variables, X (Cognitive Ability) and Y (Math Grades). In the presentation slides, there is a formula that shows how you can calculate the heritability of Y. Let's go more into detail to understand the formula. First, where does the formula "come from"? You can derive it from you model parameters using covariance algebra or the path tracing rules. So, using e.g. the path tracing rules, it is possible to show that in the bivariate ACE-beta model the variance of X (Var(X)) and the variance of Y (Var(Y)) are defined as: - `Var(X)=a11^2+c11^2+e11^2` - `Var(Y)=(b21^2*Var(X))+(2*b21*a11*a21+2*b21*c11*c21)+(a21^2+a22^2+c21^2+c22^2+e22^2)` The heritability of Y is simply the proportion of the variance of Y that is due to genetic factors. So, selecting all the terms containing a genetic effects, leads us to: - `GenVar(Y)=(b21^2*a11^2)+(2*b21*a11*a21)+(a21^2+a22^2)` So, we get: - `h2(Y)=((b21^2*a11^2)+(2*b21*a11*a21)+(a21^2+a22^2))/(var(Y))` Now, that we have the formula, let's calculate the heritability of Y. Use the formula derived above and the model parameters of model `model1_red1_fit` to derive the heritability of Y. ```{r} #Var_X <- #Var_Y <- #GenVar_Y <- #h2_Y <- ``` So, the heritability is the proportion of a variance that is explained by genetic factors. So by dividing the genetic variance by the total variance, we get the heritability. What is the genetic variance of X? What is the heritability of X? Uncomment the following code chunk and calculate both statistics. ```{r} #Var_X <- #GenVar_X <- #h2_X <- ``` ## Adding covariates Normally, one controls for age and sex when estimating twin models. If not the shared environmental effects may be too big. So, let's control for age and sex in the ACE-beta model. To do this, we added two new arguments: - `covvars`: String vector with the covariates - `covariance`: Logical to define whether covariates should be part of covariance matrix (random) or of mean vector (fixed). If you choose `covariance = TRUE` you can use FIML for the covariates as well but the estimation is more complicated because now the (co-)variances of the covariates are model parameters as well. ```{r, echo = TRUE, message=FALSE} model2a <- twinflex(acevars = c("CA","math"), sep = "_", zyg = "zyg", covvars = c("age","female"), covariance = FALSE, type = "aceb", data = df_final_wide, TryHard = TRUE) makeOutput(summary(model2a), rounded = TRUE) ``` We can also model the covariates as part of the covariance matrix. (This is also a good option if you want to control for variables that vary within twin pairs (because of endogeneity issues, remember, you assume that `e21=0`) but without decomposing them into ACE components.) ```{r} model2b <- twinflex(acevars = c("CA","math"), sep = "_", zyg = "zyg", covvars = c("age","female"), covariance = TRUE, type = "aceb", data = df_final_wide, TryHard = TRUE) makeOutput(summary(model2b), rounded = TRUE) ``` ## A stepwise modeling approach In the presentations, you heard that a stepwise modeling approach may be appropriate to tackle a question of bias and causality: We start with a "naive" model that assumes that there is no genetic nor shared environmental confounding. Statistically, this means that this model assumes `a21` and `c21` to be zero. In a second step, we estimate a more realistic model that allows for genetic and shared environmental confounding. Statistically, this means that we estimate `a21` and `c21` as free parameters. Create the full model: ```{r, echo = TRUE, message=FALSE} model3 <- twinflex(acevars = c("CA","math"), sep = "_", zyg = "zyg", covvars = c("age","female"), covariance = TRUE, type = "aceb", data = df_final_wide, TryHard = TRUE) ``` Create the naive model by fixing `a21` and `c21` to zero: ```{r, echo = TRUE, message=FALSE} model3_red1_model <- mxModel(model3, name="Drop c21 and a21") model3_red1_model <- omxSetParameters(model = model3_red1_model, labels=c("c21","a21"), free=FALSE, values=0) model3_red1_fit <- mxTryHard(model3_red1_model, exhaustive = TRUE) model3_red1_cf <- mxCompare(model3, model3_red1_fit) model3_red1_cf # p > 0.05 -> Ok ! makeOutput(summary(model3_red1_fit), rounded = T) summary(model3_red1_fit) ``` Based on the naive model, what is our naive estimate of the effect of cognitive abilities on math grades? Based on the model comparison, what is the better model? ## Tasks 1. You want to estimate the effect of externalizing problem behavior (hyperactivity) (`exthyp`) on german grades (`german`) by calculating an ACE-beta model. You want to control for age and sex in the mean vector. Please uncomment and complete the following code chunk and run the model. ```{r, echo = TRUE, message=FALSE} model4 <- twinflex(acevars = c("exthyp","german"), sep = "_", zyg = "zyg", covvars = c("age","female"), covariance = FALSE, type = "aceb", data = df_final_wide, TryHard = TRUE) makeOutput(summary(model4)) ``` Reduce the model by dropping insignificant paths from the model. Do not drop paths of the control variables/covariates age and sex. Focus on the paths of the ACE decomposed variables specified in the `acevars` argument. Uncomment, complete and run the following code chunk. ```{r, echo = TRUE, message=FALSE} model4_red1_model <- mxModel(model4, name="Drop ") model4_red1_model <- omxSetParameters(model = model4_red1_model, labels=c(), free=, values=) model4_red1_fit <- mxTryHard(model4_red1_model, exhaustive = TRUE) model4_red1_cf <- mxCompare(model4, model4_red1_fit) model4_red1_cf # p > 0.05 -> Ok ! makeOutput(summary(model4_red1_fit), rounded = T) summary(model4_red1_fit) ``` 2. You want to know whether self perceived math ability (`spabmath`) mediates the effect of cognitive ability (`CA`) on math grades (`math`). You want to be sure that you get valid estimates of the phenotypic effects so that you choose to use an ACE-beta model to control for genetic and shared environmental confounding. First, estimate a bivariate ACE-beta model to test whether there is an effect of self perceived math ability on math grades. To do so, uncomment, complete and run the following code chunk: ```{r, echo = TRUE, message=FALSE} model5 <- twinflex(acevars = c("spabmath","math"), sep = "_", zyg = "zyg", covvars = c("age","female"), covariance = TRUE, type = "aceb", data = df_final_wide, TryHard = TRUE) makeOutput(summary(model5), rounded = TRUE) ``` Reduce the model by dropping insignificant parameters. Uncomment, complete and run the following code chunk. ```{r, echo = TRUE, message=FALSE} model5_red1_model <- mxModel(model5, name="Drop ") model5_red1_model <- omxSetParameters(model = model5_red1_model, labels=c("c11","c21"), free=FALSE, values=0) model5_red1_fit <- mxTryHard(model5_red1_model, exhaustive = TRUE) model5_red1_cf <- mxCompare(model5, model5_red1_fit) model5_red1_cf # p > 0.05 -> Ok ! makeOutput(summary(model5_red1_fit), rounded = T) summary(model5_red1_fit) ``` Does self perceived math ability affect math grades when controlling for genetic and shared environmental confounding? Is there genetic or shared environmental confounding? If so, can you quantify the bias? ```{r, echo = TRUE, message=FALSE} model5 <- twinflex(acevars = c("CA","spabmath","math"), sep = "_", zyg = "zyg", covvars = c("age","female"), covariance = TRUE, type = "aceb", data = df_final_wide, TryHard = TRUE) makeOutput(summary(model5)) ``` ```{r, echo = TRUE, message=FALSE} model5_red1_model <- mxModel(model5, name="Drop ") model5_red1_model <- omxSetParameters(model = model5_red1_model, labels=c("a32","c32","b21","c22","c21","c31"), free=FALSE, values=0) model5_red1_fit <- mxTryHard(model5_red1_model, exhaustive = TRUE) model5_red1_cf <- mxCompare(model5, model5_red1_fit) model5_red1_cf # p > 0.05 -> Ok ! makeOutput(summary(model5_red1_fit), rounded = T) summary(model5_red1_fit) ```