GxE Interaction affects heritability of IQ depending on socio-economic status

Introduction

A large genotyped sample of the population of Wisconsin, USA (The Wisconsin Longitudinal Study, N=8,509) are examined for evidence of the Scarr-Rowe effect, an adverse gene x environment (GxE) interaction that reduces the heritability of IQ among those with low childhood socioeconomic status (SES), using the Continuous Parameter Estimation Method (CPEM):

\[\begin{equation} SES \sim PGS \cdot IQ \end{equation}\]

A large correlation indicates increasing genetic expressivity with increasing SES

Imports and settings

library(pacman)
p_load(knitr, tidyverse, magrittr, moments, psych, corrr, pander, broom, glue, ggpubr)
PanderOpts(digits = 3, table.split.table = Inf, missing = "-", table.alignment.default = "left")
theme_set(theme_bw(base_size=14))

Read data

A data set with PGS is prepared here. See link for details about data preparation.

df <- read_delim("pgs_data.csv", delim=";") %>% 
  select(idpub, pgs_std, iq_std = z_iq_std, iq100 = z_iq100, 
         ses57, sexrsp=z_sexrsp, rtype) %>% 
  na.omit()

tribble(~"Sample size", df %>% nrow)

Find CPEM, the inner product og pgs and iq

With standardized versions of all variables

df %<>% mutate(
  ses57_std = as.numeric(scale(ses57)),
  pgs_iq = as.numeric(scale(pgs_std * iq_std)),
  logses57_std = as.numeric(scale(log(ses57))))

Regression of pgs_iq with ses

m <- df %$% lm(pgs_iq ~ 0 + ses57_std)
m

Scatter plot of ses vs pgs_iq

ggplot(df) +
  geom_point(aes(y=pgs_iq, x=ses57_std), color="grey40", size=1, alpha=0.3) + 
  geom_abline(aes(intercept=0, slope=coef(m)[["ses57_std"]])) +
  ylab("CPE (IQxEA3)") + xlab("Parental SES")

Check normality of ses variable

Since the R square is low, these are essentially also the results for the residuals

plot_normality <- function(df, var){
  p1 <-ggplot(df, aes_string(x=var)) + 
    geom_histogram(fill="turquoise4", col="black", binwidth=0.4)
  p2 <- ggplot(df, aes_string(sample=var)) + stat_qq()
  
  ggarrange(p1, p2, ncol=2)
}

plot_normality(df, "ses57_std")

Analysis with log(ses)

df %$% lm(pgs_iq ~ logses57_std)

Check normality of log ses variable

plot_normality(df, "logses57_std")

Regressions for all combinations of sex, intercept, log transformation and outlier removal

rm_outliers <- function(df, sd){
  df %>% filter_at(vars(iq_std, ses57_std, pgs_iq), all_vars(between(., -sd, sd)))
}

getregr <- function(ses, sex, outlier_rm, intercept, df){
  intercept_s <- if (intercept == "0") "0 + " else ""

  df <- if(sex == "both") df else df %>% filter(sexrsp == sex)
  t_index <- if (intercept == "normal") 2 else 1

  if (outlier_rm == "SD < 3") df %<>% rm_outliers(3)
  else if (outlier_rm == "SD < 4") df %<>% rm_outliers(4)

  m <- lm(as.formula(glue("pgs_iq ~ {intercept_s} {ses}")), df)
  est <- tidy(m) %>% filter(term == ses)
  tibble(ses, outlier_rm, sex, intercept,
         estimate = est$estimate,
         t.value = tidy(m)$statistic[t_index],
         std.error = est$std.error,
         fstat = glance(m)$statistic,
         adjr2 = glance(m)$adj.r.squared,
         skewness = skewness(resid(m)),
         p.value = est$p.value)
}

crossing(ses = c("ses57_std", "logses57_std"), 
         sex = c("male", "female", "both"),
         outlier_rm = c("SD < 3", "SD < 4", "none"),
         intercept = c("normal", "0")) %>% 
  pmap(getregr, df) %>% bind_rows %>%
  arrange(desc(ses), outlier_rm)

Decile means (>4SD outliers removed)

df %>% rm_outliers(4) %>% mutate(decile = ntile(ses57_std, 10)) %>% 
  group_by(decile) %>%
  summarize(
    decile_sd_iq = sd(iq_std),
    decile_sd_pgs = sd(pgs_std),
    decile_mean_pgs_iq = mean(pgs_iq),
    cor = cor(pgs_std, iq_std))

Correlation between variables

show_cor <- function(df, title) {
  df %>% select(iq_std, ses57_std, pgs_std) %>% 
    correlate(quiet=T) %>% 
    pander(caption=title)
}

df %>% show_cor("Full sample")
df %>% filter(sexrsp == "male") %>% show_cor("Males")
df %>% filter(sexrsp == "female") %>% show_cor("Females")

Descriptions of variables

tribble(~variable, ~code_name, ~WLS_name,
       "IQ", "iq_std", "gwiiq_bm for graduates and swiiq_t for siblings",
       "Polygenic score", "pgs_std", "pgs_ea3_mtag",
       "Socio-economic status of family", "ses57", "ses57")

df %>% select(iq_std, pgs_std, ses57) %>% describe()

CPEM vs linear regression

Under which conditions is CPEM more powerful than using linear regression and looking at the interaction term? I perform some simulations to investigate this question here