# Clear Workspace
rm(list = ls())

# Switch off scientific notation 
options(scipen=999)

# Install packages
# install.packages("tidyverse") # for recoding variables
# install.packages("haven") # for load dta files
# install.packages("lme4") # for mixed effects modeling
# install.packages("lmerTest") # for sig. stars in model results
# install.packages("performance") # to compute ICC statistics
# install.packages("plm") # for FE modeling
# install.packages("sjPlot") # for plots
# install.packages("effects") # for plots
# install.packages("Hmisc") # for correlations
# install.packages("lavaan") # for SEM

# Load packages - if not already installed, use install.packages("packagename")
library(tidyverse) 
library(haven) 
library(lme4) 
library(lmerTest) 
library(performance) 
library(plm) 
library(sjPlot) 
library(effects) 
library(Hmisc) 
library(lavaan)

# Load useful functions
source("Session 2/Scripts/recode.R")

########################
# Create a balanced data frame

df_final <- df_final %>% 
  group_by(fid) %>% 
  mutate(ntwins = sum(n())) %>% 
  ungroup()

df_final <- df_final %>% 
  filter(ntwins == 2)

df_final <- df_final %>% 
  select(-ntwins)

# check if data frame is balanced
is.pbalanced(df_final, index = c("fid","ptyp"))

########################
# Inspect the data frame
########################
summary(df_final)
# Note: Variable zyg -> 0 = MZ; 1 = DZ

########################
# within-transformation of predictors
########################

df_final <- df_final %>% 
  group_by(fid) %>% 
  mutate(CA_m = mean(CA, na.rm = TRUE)) %>% 
  mutate(CA_dev = CA-CA_m) %>% 
  mutate(notlazy_m = mean(notlazy, na.rm = TRUE)) %>% 
  mutate(notlazy_dev = notlazy-notlazy_m) %>% 
  mutate(nervous_m = mean(nervous, na.rm = TRUE)) %>% 
  mutate(nervous_dev = nervous-nervous_m) %>% 
  mutate(notquiet_m = mean(notquiet, na.rm = TRUE)) %>% 
  mutate(notquiet_dev = notquiet-notquiet_m) %>% 
  mutate(knowledge_m = mean(knowledge, na.rm = TRUE)) %>% 
  mutate(knowledge_dev = knowledge-knowledge_m) %>% 
  mutate(selfeff_m = mean(selfeff, na.rm = TRUE)) %>% 
  mutate(selfeff_dev = selfeff-selfeff_m) %>% 
  mutate(worries_m = mean(worries, na.rm = TRUE)) %>% 
  mutate(worries_dev = worries-worries_m) %>% 
  mutate(help_m = mean(help, na.rm = TRUE)) %>% 
  mutate(help_dev = help-help_m) %>% 
  ungroup()

summary(df_final)

########################
# Analyze the data
########################

###########
# Model 1: Random-Intercept-Only
###########
mem1 <- lmer(formula = math ~ 1 + (1|fid),
             data    = df_final) 
# Model summary
summary(mem1)
  # Variance on the family level = 0.4342
  # Variance on the individual level = 0.3785

# ICC
icc(mem1)
# ICC by hand
0.4342/(0.4342+0.3785)


###########
# Model 2: Add cognitive ability as a predictor
###########
mem2 <- lmer(formula = math ~ 1 + CA + (1|fid),
             data    = df_final) 
summary(mem2) # cross-sectional effect = 0.25427

# compare model1 and model2 with LR test
anova(mem1, mem2) # significant better fit to the data


###########
# Model 3: Estimating the within-pair effect
###########
# Model 3a: Mixed effects parametrization
mem3 <- lmer(formula = math ~ 1 + CA_dev + CA_m +  (1|fid),
             data    = df_final) 

# Model 3b: Fixed effects parametrization
fem3 <- plm(formula = math ~ CA , data = df_final, # fem3 <- grun.fe <- plm
                       index = c("fid","ptyp"),
                       model = "within")

summary(mem3) # within-pair effect = 0.17948
summary(fem3) # within-pair effect = 0.179479

###########
# Model 4: Add zygosity covariate
###########
mem4 <- lmer(formula = math ~ 1 + CA_dev + CA_m + zyg + (1|fid),
             data    = df_final) 
summary(mem4) # As expected: there are no zygosity differences in math grades

###########
# Model 5: Decompose the confounders via an interaction between zygosity and within-pair effect
###########
# Model 5a: Mixed effects parametrization
mem5 <- lmer(formula = math ~ 1 + CA_dev + CA_m + zyg + CA_dev*zyg + (1|fid),
             data    = df_final) 

# Model 5b: Fixed effects parametrization
fem5 <- plm(formula = math ~ CA*zyg , data = df_final, 
            index = c("fid","ptyp"),
            model = "within") 
summary(mem5) 
summary(fem5) 

# MZ within-pair effect                           = Main Effect                    = IQ_dev              = 0.08531*
# Difference between MZ and DZ within-pair effect = Interaction term               = IQ_dev:zygDZ        = 0.13813**
# DZ within-pair effect                           = Main Effect + Interaction term = IQ_dev+IQ_dev:zygDZ = +  = 0.08531 + 0.13813 = 0.22344

###########
# Plots
###########

# Plot of all effect sizes
sjPlot::plot_model(mem5)
# with labels
sjPlot::plot_model(mem5, 
                   axis.labels=c("Diff Within-pair DZ-MZ", "Zygosity (DZ)","CA (Between)","CA (Within-pair) MZ"),
                   show.values=TRUE, show.p=TRUE,
                   title="Effects on Math Grades")

# Prediction plots of interaction effects
eff.p1 <- effect("CA_dev*zyg", mem5)
plot(eff.p1)
plot_model(mem5, type = "pred", terms = c("CA_dev","zyg"),
           title = "Differences in Within-effect by zygosity",
           axis.title = c("CA (demeaned)","Predicted math grade"))

# for more information, see ?sjPlot or https://strengejacke.github.io/sjPlot/articles/plot_interactions.html
# alternatively, you could access the model object and use ggplot

########################################################################################################################

## Practice

# 1. The effect of cognitive ability may be still biased by uncontrolled factors that vary within twin pairs.
  # Control for additional covariates. Interact the covariate with zygosity to see whether the effect varies by zygosity

# 2. Use a different predictor and a different outcome. What do you find? Is the effect biased? If so, by genetic or shared environmental factors?
    # For example: Use a measure of parental investment (help) as the outcome in order to see whether parental investment depends on some child characteristics.
# 3. Do the effects vary by sex or SES? (you could try to answer this question with a three-way interaction)