why-natural-selection.Rmd

---
title: "Human capital mediates natural selection in contemporary humans"
author: "David Hugh-Jones \\thanks{Corresponding author. School of Economics, University of East Anglia, Norwich, UK. Email: D.Hugh-Jones@uea.ac.uk}, Abdel Abdellaoui \\thanks{Department of Psychiatry, Amsterdam UMC, University of Amsterdam, Amsterdam, The Netherlands. Email: a.abdellaoui@amsterdamumc.nl}"
date: "`r format(Sys.Date(), '%d %B %Y')`"
abstract: "\\noindent Natural selection has been documented in contemporary humans, but  little is known about the mechanisms behind it. We test for natural selection through the association between 33 polygenic scores and fertility, across two generations, using data from UK Biobank (N = 409,629 British subjects with European ancestry). Consistently over time, polygenic scores that predict higher earnings, education and health also predict lower fertility. Selection effects are concentrated among lower SES groups, younger parents, people with more lifetime sexual partners, and people not living with a partner. The direction of natural selection is reversed among older parents, or after controlling for age at first live birth. These patterns are in line with the economic theory of fertility, in which earnings-increasing human capital may either increase or decrease fertility via income and substitution effects in the labour market. Studying natural selection can help us understand the genetic architecture of health outcomes: we find evidence in modern day Great Britain for multiple natural selection pressures that vary between subgroups in the direction and strength of their effects, that are strongly related to the socio-economic system, and that may contribute to health inequalities across income groups."
bibliography: "why-natural-selection.bib"
output: 
  bookdown::pdf_document2:
    toc: false
    latex_engine: xelatex
    number_sections: true
    keep_tex: true
header-includes:
  - \usepackage{subfig}
  - \usepackage{setspace}\doublespacing
  - \usepackage{placeins}
  - \usepackage[format=plain, labelfont={bf,it}, textfont=it]{caption}
  - \usepackage{titlesec}
  - \titleformat*{\section}{\sffamily\LARGE}
  - \titleformat*{\subsection}{\sffamily\itshape\Large}
  - \hypersetup{colorlinks = true, linkcolor = {blue}, linkbordercolor = {white}}
  - \usepackage{amsthm}
  - \theoremstyle{plain}
  - \newtheorem{lem}{\protect\lemmaname}
  - \providecommand{\lemmaname}{Lemma}
  - \usepackage{etoc}
  - \usepackage{lineno}
  - \linenumbers
editor_options: 
  markdown: 
    wrap: 72
  chunk_output_type: console
mainfont: Baskerville
mathfont: Baskerville
sansfont: "Gill Sans"
---


```{r setup, include = FALSE}

library(drake)
library(magrittr)
library(dplyr)
library(forcats)
library(ggplot2)
library(tidyr)
library(santoku) # keep this after tidyr as it masks chop
library(purrr)
library(glue)
library(broom)
library(huxtable)
loadNamespace("scales")
loadNamespace("shades")

# nasty hack to make add_ashe_income work later:
if (! exists("add_ashe_income")) source("~/import-ukbb-data/import-ukbb-data.R")

drake::loadd(rgs)
drake::loadd(score_names)
drake::loadd(famhist)

options(digits = 2)
options(dplyr.summarise.inform = FALSE) 

knitr::opts_chunk$set(
        echo       = FALSE, 
        warning    = FALSE, 
        cache      = FALSE, 
        error      = FALSE,
        fig.height = 3.5,
        dev        = "cairo_pdf"
      )

knitr::knit_hooks$set(
  inline = function (x) {
    if (is.numeric(x)) x <- as.character(round(x, getOption("digits")))
    x <- gsub("-", "\u2212", x)
    paste(as.character(x), collapse = ", ")
  }
)
                      
options(huxtable.long_minus = TRUE)
huxtable::set_default_properties(latex_float = "h!")

theme_set(theme_minimal())
theme_update(
  text = element_text(family = "Abadi MT Condensed Light")
)

my_hline <- geom_hline(yintercept = 0, colour = "grey20", linetype = "dotted")
my_vline <- geom_vline(xintercept = 0, colour = "grey20", linetype = "dotted")


standard_ggplot <- function (
  dfr, 
  fill_col, 
  n_regs, 
  ...,
  score_col = quo(term), 
  order_idx = 1, 
  fill_direction = 1,
  n_cats    = NULL,
  conf_int  = NULL
) {
  if (! missing(score_col)) score_col <- enquo(score_col)
  
  mfc <- missing(fill_col)
  fill_col <- if (mfc) quo(NULL) else enquo(fill_col)
  if (missing(n_cats)) {
    n_cats <- if (mfc) 1 else length(unique(pull(dfr, {{fill_col}})))
  }
  
  if (is.null(conf_int)) conf_int <- "conf.low" %in% names(dfr)
  
  if (conf_int) {
    confint_scale <- scale_linetype_manual(
                       name   = "", 
                       values = c("95% c.i. uncorrected" = 1),
                       guide  = guide_legend(order = 3)
                     )
    confint_geom_segment <- geom_segment(
                              aes(
                                x        = conf.low, 
                                xend     = conf.high,
                                linetype = "95% c.i. uncorrected"
                              ), 
                              color   = "grey45", 
                              alpha   = 0.2, 
                              size    = 1.1, 
                              lineend = "round"
                            )
  } else {
    confint_scale <- NULL
    confint_geom_segment <- NULL
  }
  
  my_shape_fill_guide <- guide_legend(order = 1)
  
  n_regs <- as.double(n_regs)
  shape_values <- paste(c("circle", "square", "triangle", "diamond", 
                            "triangle down"), "filled")
  shape_values <- shape_values[1:n_cats]
  if (fill_direction == -1) shape_values <- rev(shape_values)
  dfr %>% mutate(
    !! score_col := pretty_names(!! score_col),
    !! score_col := fct_reorder(!! score_col, estimate, order_abs(order_idx))
  ) %>% 
    rename(p = p.value) %>% 
    ggplot(aes(estimate, {{score_col}}, yend = {{score_col}}, 
               shape = {{fill_col}}, ...)) +
    my_vline + 
    confint_geom_segment +
    geom_point(size = 1.8, alpha = 0.8, aes(color = {{fill_col}}, 
                 fill = stage(p < 0.05/{{n_regs}}, after_scale = ifelse(
                 fill == "white", fill, color)))) +
    my_fill_scale(aesthetics = "color", n = n_cats, direction = fill_direction, 
                    guide = my_shape_fill_guide) +
    scale_fill_manual(values = c("TRUE" = "black", "FALSE" = "white"), 
                        guide = guide_legend(
                                  title = sprintf("p < 0.05/%s", n_regs),
                                  override.aes = list(shape = "circle filled"),
                                  order = 2
                                )) + 
    scale_shape_manual(values = shape_values, guide = my_shape_fill_guide) +
    my_labs() +
    my_theme(
      legend.justification = "top", 
      panel.grid.major.y   = element_blank()
    ) +
    confint_scale
}


my_theme <- function (...) {
  theme(
          axis.text.y      = element_text(size = 7),
          panel.grid.major = element_line(size = rel(0.5)),
          ...
        )
}

my_fill_scale <- function (n, aesthetics = "fill", direction = 1, ...) {
  
  fill_colors <- c("steelblue4", "darkred")
  
  scale <- if (n <= 2) {
    dots <- list(...)
    dots$direction <- NULL
    do.call(scale_colour_manual, c(list(values = fill_colors[1:n],
                                        aesthetics = aesthetics), 
                                   dots))
  } else {
    # scale <- scale_color_viridis_d(aesthetics = aesthetics, option = "D",
    #                                  direction = direction, end = 0.6, ...)
    fill_colors <- c(fill_colors[1], "#008b8b", fill_colors[2])
    if (direction != -1) fill_colors <- rev(fill_colors)
    scale_color_manual(aesthetics = aesthetics, 
                         values = shades::gradient(fill_colors,
                                                     steps = n, space = "Lab"), 
                         ...) 
  }
  
  scale
}


my_labs <- function (x = "Effect size", y = "", ...) labs(x = x, y = y, ...)


# for use when reordering factors. Order by the nth score in a group
order_abs <- function (n = 1) {
  function (x) x[n] 
}


comma_num <- function (x) prettyNum(x, big.mark = ",")


pretty_names <- function (names) {
  pretty <- c(
    ADHD_2017               = "ADHD",
    age_at_menarche         = "Age at menarche",
    age_at_menopauze        = "Age at menopause",
    agreeableness           = "Agreeableness",
    ai_substance_use        = "Age at smoking initiation",
    alcohol_schumann        = "Alcohol use",
    alzheimer               = "Alzheimer",
    autism_2017             = "Autism",
    bipolar                 = "Bipolar",
    bmi_combined            = "BMI",
    body_fat                = "Body Fat",
    caffeine                = "Caffeine",
    cannabis                = "Cannabis (ever vs. never)",
    cognitive_ability       = "Cognitive Ability",
    conscientiousness       = "Conscientiousness",
    coronary_artery_disease = "Coronary Artery Disease",
    cpd_substance_use       = "Cigarettes per day",
    diagram_T2D             = "Type 2 Diabetes",
    dpw_substance_use       = "Drinks per week",
    EA2_noUKB               = "Educ. attainment 2 (no UKBB)",
    EA3_excl_23andMe_UK     = "Educ. attainment 3 (no UK)",
    eating_disorder         = "Eating disorder",
    extraversion            = "Extraversion",
    height_combined         = "Height",
    hip_combined            = "Hip circumference",
    MDD_PGC2_noUKB          = "Major Depressive Disorder",
    neuroticism             = "Neuroticism",
    openness                = "Openness",
    sc_substance_use        = "Smoking cessation",
    SCZ2                    = "Schizophrenia",
    si_substance_use        = "Smoking initiation",
    wc_combined             = "Waist circumference",
    whr_combined            = "Waist-hip ratio"
  )
  
  pretty[names]
}


n_in_regs <- function (var, dv = "RLRS", data = famhist) {
  fh_subset <- data
  if (dv == "RLRS") fh_subset %<>% filter(kids_ss)
  # all PGS have the same NA pattern, so we just use one
  tbl <- table(fh_subset[
                  ! is.na(fh_subset[dv]) & ! is.na(fh_subset$whr_combined), 
                  var
                ])
  
  c(tbl)
}


add_n <- function (var, famhist_var = var, dv = "RLRS", data = famhist, 
                     reverse = FALSE) {
  var <- as.factor(var)
  n <- n_in_regs(famhist_var, dv, data)
  # this is a horrible hack which just happens to work
  # check all Ns manually!
  if (reverse) n <- rev(n)
  levels(var) <- sprintf("%s (N = %s)", levels(var), 
                         prettyNum(n[seq_len(nlevels(var))], big.mark = ","))
  # var_n <- sprintf("%s (N = %s)", var, prettyNum(n[var], big.mark = ","))
  # var_n[is.na(var)] <- NA
  var
}


# from https://github.com/dgrtwo/drlib/blob/master/R/reorder_within.R
# for ordering within facets
reorder_within <- function(x, by, within, fun = mean, sep = "___", ...) {
  new_x <- paste(x, within, sep = sep)
  stats::reorder(new_x, by, FUN = fun)
}


scale_x_reordered <- function(..., sep = "___") {
  reg <- paste0(sep, ".+$")
  ggplot2::scale_x_discrete(labels = function(x) gsub(reg, "", x), ...)
}

scale_y_reordered <- function(..., sep = "___") {
  reg <- paste0(sep, ".+$")
  ggplot2::scale_y_discrete(labels = function(x) gsub(reg, "", x), ...)
}

```

\normalem


# Introduction

Living organisms evolve through natural selection, in which allele frequencies
change in the population through differential reproduction rates. Studying the
mechanisms behind natural selection can help us better understand how individual
differences in complex traits and disease risk arise [@benton2021influence].
Recent work confirms that natural selection is taking place in modern human
populations, using genome-wide analysis
[@Barban_2016;@beauchamp2016genetic;@conley2016assortative;@kong2017selection;@sanjak2018evidence;@FIEDER202216]. In
particular, genetic variants associated with higher educational attainment are
being selected against, although effect sizes appear small.

As yet we know little about the social mechanisms behind natural selection. The
economic theory of fertility [@becker1960economic] offers a potential
explanation. Higher potential earnings have two opposite effects on fertility: a
fertility-increasing *income effect* (higher income makes children more
affordable), and a fertility-lowering *substitution effect* (time spent on
childrearing has a higher cost in foregone earnings). Thus, an individual's
*human capital* -- skills and personality traits which are valuable in labour
markets -- can increase or decrease their fertility. Genetic variants  which are 
linked to human capital will then be selected for or against. Also, the economic
theory predicts that the relative strength of income and substitution effects
will vary systematically across different social groups.

This study uses data from UK Biobank [@bycroft2018uk] to learn more about
contemporary natural selection. We test for natural selection on 33 different
polygenic scores by estimating their correlation with fertility. We extend the
analysis over two generations, using data on respondents' number of siblings as
well as their number of children. This is interesting because consistent natural
selection over multiple generations could lead to substantive effects in the
long run. Next, we examine correlations with fertility in different subgroups.
Across the board, selection effects are stronger in groups with lower income and
less education, among younger parents, people not living with a partner, and
people with more lifetime sexual partners. Outside these groups, effects are
weaker and often statistically insignificant. In some subgroups, the direction
of selection is even reversed.
 
We then show that a simple model of human capital, education and fertility
choices can give rise to these empirical results. At higher incomes, the income
and substitution effects are balanced, while among lower-income people, or
single parents who face a bigger time burden from childcare, the substitution
effect dominates. The theory predicts that polygenic scores' correlation
with fertility is associated with their correlation with education and earnings,
and we confirm this. We then run a mediation analysis, which shows that part of the
correlation with fertility is indeed mediated by educational attainment. 
Thus, contemporary natural selection on polygenic scores can be explained 
by scores' correlation with earnings-increasing human capital.

Lastly, we discuss the effects of natural selection. While our estimated effects 
on measured polygenic scores are small, natural selection substantially 
increases the correlation between polygenic scores and income, increasing 
genetic differences between different social groups, and thus making the 
"genetic lottery" [@harden2021genetic] more unfair.


# Results


We created polygenic scores for 33 traits in `r comma_num(nrow(famhist))`
individuals of European descent, corrected for ancestry using 100 genetic
principal components (see Materials and Methods). Figure
\@ref(fig:plot-means-over-time) plots mean polygenic scores in the sample by
5-year birth intervals. Several scores show consistent increases or declines
over this 30-year period, of the order of 5% of a standard deviation. These
changes could reflect natural selection within the UK population, but also
emigration, or ascertainment bias in the sample [@10.1093/aje/kwx246]. 


```{r plot-means-over-time, fig.height = 7, fig.cap = "Mean polygenic scores (PGS) by birth year in UK Biobank. Symbols show means for 5-year intervals. Bars are 95% confidence intervals. Triangles denote a significant linear increase or decrease over time (p < 0.05/33)."}

drake::loadd(pgs_over_time)

pgs_over_time %>% 
      mutate(score_name = pretty_names(score_name)) %>% 
      ggplot(aes(YOB, group = score_name, colour = Change, shape = Change)) + 
      geom_linerange(aes(ymin = score_lo, ymax = score_hi), color = "grey70", 
                       size = 3) +
      geom_point(aes(y = score), size = 1.5) +
      facet_wrap(~score_name, ncol = 5) +
      my_hline +
      scale_color_manual(values = c("-" = "darkred", "+" = "darkgreen", "o" = "black")) +
      scale_shape_manual(values = c(
                                      "-" = "triangle down open", 
                                      "+" = "triangle open", 
                                      "o" = "circle open"
                                    )) +
      theme(
        line        = element_line(size = 0.5, lineend = "round"),
        strip.text  = element_text(size = 8), 
        axis.text.x = element_text(size = 7),
        legend.position = "none"
      ) +
      labs(y = "PGS", x = "Birth year (5 year intervals)")

```


To test for natural selection more directly, we regress respondents' relative 
lifetime reproductive success (RLRS) on each polygenic score (PGS):


```{=tex}
\begin{equation}
\mathrm{RLRS}_i = \alpha + \beta\mathrm{PGS}_i + \varepsilon_i
\label{eq:regression}
\end{equation}
```

RLRS is defined as respondent $i$'s number of children, divided by the mean number
of children of people born in the same year. The "selection effect", $\beta$,
reflects the strength of natural selection within the sample. In fact, since
polygenic scores are normalized, $\beta$ is the expected polygenic score among
children of the sample [@beauchamp2016genetic].[^normalization] Note that
equation \eqref{eq:regression} does not control for many environmental and genetic
factors that could affect fertility, and as a result, $\beta$ is not an estimate
of the causal effect of a polygenic score on fertility. However, natural
selection is a matter of correlation not causation: polygenic scores which 
correlate with high fertility are being selected for, whatever the underlying 
causal mechanism.

[^normalization]:The selection effect $\beta$ equals $Cov(RLRS, PGS)/Var(PGS)$.
Since PGS are normalized to variance 1 and mean 0, this reduces to $Cov(RLRS,
PGS) = E(RLRS\times PGS) - E(RLRS) E(PGS) = E(RLRS \times PGS)$. This is the
polygenic score weighted by relative lifetime reproductive success, which is the
average polygenic score in the next generation [@robertson1966mathematical].


```{r calc-res-weighted}

drake::loadd(res_wt_flb_weights)
drake::loadd(res_wt_msoa_weights)
drake::loadd(res_wt_age_qual_weights)
drake::loadd(res_wt_van_alten_weights)
drake::loadd(res_sibs_parent_weights)
drake::loadd(res_unweighted)
drake::loadd(res_sex)

res_wt_combined <- bind_rows(
  None                = res_unweighted %>% 
                          select(-score_name) %>% 
                          filter(term != "(Intercept)"),
  Main                = res_wt_van_alten_weights,
  Geographical        = res_wt_msoa_weights,
  `Age/Qualification` = res_wt_age_qual_weights,
  `Age/Qual/AFLB`     = res_wt_flb_weights,
  .id = "Weights"
)

n_sig <- sum(res_wt_van_alten_weights$p.value < 0.05/33)

res_wt_cb_wide <- res_wt_combined %>% 
                    tidyr::pivot_wider(
                      id_cols     = term,
                      names_from  = Weights,
                      values_from = estimate
                    ) %>% 
                    rowwise() %>% 
                    mutate(
                      consistent = all(c_across(-term) > 0) ||
                        all(c_across(-term) < 0)
                    ) %>% 
                    ungroup() %>% 
                    left_join(
                      res_wt_combined %>% 
                        filter(Weights == "None") %>% 
                        select(term, p.value),
                      by = "term"
                    ) %>% 
                    left_join(
                      res_sex %>% 
                        filter(sex == "Female") %>% 
                        select(term, estimate_females = estimate),
                      by = "term"
                    )

rel_sizes <- res_wt_cb_wide %>% 
                filter(
                  consistent, 
                  p.value < 0.05/length(score_names)
                ) %>% 
                mutate(
                  Geographical        = Geographical/None,
                  `Age/Qualification` = `Age/Qualification`/None,
                  Main                = Main/None,
                  `Age/Qual/AFLB`     = `Age/Qual/AFLB`/estimate_females
                ) %>% 
                select(-consistent, -None, -p.value, -estimate_females) %>% 
                mutate(
                  term = pretty_names(term)
                )

rel_sizes_summary <- rel_sizes %>% 
                summarize(
                  across(-term, ~ c(mean(.x), median(.x)))
                ) %>% 
                mutate(
                  term = c("Mean", "Median")
                )

mean_rel <- c(rel_sizes_summary[1, 1:3])
```

 
Figure \@ref(fig:plot-res-wt) plots selection effects in the whole
sample.[^balance-diverse-def] To correct for ascertainment bias, we use
participant weights from @vanalten2022reweighting, which match the UK Biobank eligible
population on sex, birth year, location, education, employment, health,
household size and tenure, number of cars and age at death. Weighting makes a
large difference: effect sizes go up by a mean of 
`r scales::percent(mean_rel$Main - 1)`.[^weighting] `r n_sig` out of 33 weighted
selection effects are significant at $p$ < 0.05/33.

[^balance-diverse-def]: We also check for stabilizing and disruptive selection by
estimating \@ref(eq:regression) with a quadratic term. Stabilizing selection
selects for intermediate values, while disruptive selection selects for extreme 
values. In particular, we find disruptive selection for educational attainment
polygenic scores: at higher values of these scores, the negative effect on
fertility is smaller (Appendix Figure \@ref(fig:plot-purifying)).

[^weighting]: We use these weights throughout. All our qualitative results are
robust if we run unweighted regressions. Appendix Table \@ref(tab:tbl-res-wt)
shows results from alternative weighting schemes. 


```{r plot-res-wt, fig.cap = "Selection effects: weighted and unweighted regressions.  Each point represents a single bivariate regression of RLRS on a polygenic score. P value threshold is 0.05, Bonferroni-corrected for multiple comparisons. Confidence intervals are uncorrected."}

n_regs <- as.double(length(score_names))

res_wt_combined %>% 
      filter(Weights %in% c("Main", "None")) %>% 
      mutate(
        Weights  = recode(Weights, "Main" = "Weighted", "None" = "Unweighted"),
      ) %>% 
      standard_ggplot(fill_col = Weights, n_regs = n_regs, fill_direction = -1,
                        order_idx = 2)

```


We now show the empirical puzzles which motivate our economic model. Each
concerns differences in the strength of natural selection across different
subgroups in the sample. We re-estimate \@ref(eq:regression) splitting the
sample by demographic and social variables, including income and education, and
family structure variables including age at first live birth, presence of a
partner, and lifetime number of sexual partners.


```{r calc-pct-diffs}
drake::loadd(res_edu)
drake::loadd(res_income)
drake::loadd(res_partners_joint)
drake::loadd(res_with_partner)
res_partners_joint %<>% filter(grepl(":", term))
res_with_partner %<>% filter(grepl(":", term))

res_edu_wide <- full_join(
                  res_edu %>% filter(age_fte_cat == "< 16"),
                  res_edu %>% filter(age_fte_cat == "> 18"),
                  by = "score_name"
                )

res_income_wide <- full_join(
                     res_income %>% filter(income_cat == 1),
                     res_income %>% filter(income_cat == 5),
                     by = "score_name"
                   )

res_partners_wide <- full_join(
                       res_partners_joint %>% filter(grepl("lo_partnersFALSE", term)),
                       res_partners_joint %>% filter(grepl("lo_partnersTRUE", term)),
                       by = "score_name"
                     )

res_with_partner_wide <- full_join(
                           res_with_partner %>% 
                             filter(grepl("with_partnerFALSE", term)),
                           res_with_partner %>% 
                             filter(grepl("with_partnerTRUE", term)),
                           by = "score_name"
                         )

median_pct_diff <- function (dfr) {
   prop <- dfr %>% filter(
                    p.value.x < 0.05/33,
                    sign(estimate.x) == sign(estimate.y)
                  ) %>% 
                  mutate(
                    pct_diff = estimate.x/estimate.y
                  ) %>% 
                  pull(pct_diff) %>% 
                  median()
   # percentage difference:
   prop - 1
}

pct_diff_edu <- median_pct_diff(res_edu_wide)
pct_diff_income <- median_pct_diff(res_income_wide)
pct_diff_n_partners <- median_pct_diff(res_partners_wide)
pct_diff_with_partner <- median_pct_diff(res_with_partner_wide)
```

Figure \@ref(fig:plot-income-educ-level) plots selection effects for each
polygenic score, grouping respondents by age of completing full-time education,
and by household income. Effects are larger and more significant for the lowest
education category, and for the lowest income category. The median percentage
difference between the lowest and highest education categories, among scores
which are significant for the lowest category and have the same sign across
categories, is `r scales::percent(pct_diff_edu)`. Between the lowest and
highest income categories, it is `r scales::percent(pct_diff_income)`. These
results are robust to controlling for respondents' age (Appendix section
\@ref(sec:age-control)). Turning to family structure, we split respondents by
lifetime number of sexual partners, at the median value of 3 (Figure
\@ref(fig:plot-n-partners)a). Now, selection effects are larger and more
significant among those with more than 3 lifetime partners, with a median percentage
difference of `r scales::percent(pct_diff_n_partners)`. Next we split
respondents by whether they were living with a spouse or partner at the time of
interview (Figure \@ref(fig:plot-n-partners)b).  Effects are larger among those not living with a spouse or partner. The median percentage difference is 
`r scales::percent(pct_diff_with_partner)`.[^sex]

[^sex]: The same pattern holds if we analyse men and women separately (Appendix Figure
\@ref(fig:plot-n-partners-sex)). We also directly compared selection effects 
between men and women (Appendix Figure \@ref(fig:plot-sexes)). 


```{r plot-income-educ-level, fig.cap = "Selection effects by education and income.", fig.subcap = c("Age left full-time education", "Household income"), fig.ncol = 1}

drake::loadd(res_edu)
n_regs <- as.double(nrow(res_edu))

res_edu %>% 
      mutate(
        "Age left FTE" = fct_relevel(age_fte_cat, "< 16", "16-18", "> 18"),
        "Age left FTE" = add_n(`Age left FTE`, "age_fte_cat")
      ) %>% 
      standard_ggplot(n_regs = n_regs, fill_col = `Age left FTE`)
   
drake::loadd(res_income)

n_regs <- as.double(nrow(res_income))
res_income %>% 
      mutate(
        Income = factor(income_cat, 
          labels = c("< £18K", "£18-30K", "£31-51K", "£52-100K", "> £100K")
        ),
        Income = add_n(Income, "income_cat")
      ) %>% 
      standard_ggplot(n_regs = n_regs, fill_col = `Income`)
```


```{r plot-n-partners, fig.cap = "Selection effects by number of sexual partners and presence of a partner.", fig.subcap = c("Lifetime number of sexual partners", "Presence of a partner"), fig.ncol = 1, fig.align = "center", fig.height = 4}

drake::loadd(res_partners_joint)
res_partners_joint %<>%       
      filter(grepl(":", term)) %>% 
      mutate(`N partners` = ifelse(grepl("TRUE:", term), "3 or less", "4 or more")) 


n_regs <- as.double(nrow(res_partners_joint))
res_partners_joint %>% 
      mutate(
        `N partners` = add_n(`N partners`, "n_partners_split", reverse = TRUE,
                               data = famhist %>%
                               mutate(n_partners_split = n_partners <= 3))
      ) %>% 
      standard_ggplot(fill_col = `N partners`, n_regs = n_regs, 
                        score_col = score_name, fill_direction = -1)


drake::loadd(res_with_partner)
res_with_partner %<>% filter(grepl(":", term)) 

n_regs <- as.double(nrow(res_with_partner))

res_with_partner %>% 
      mutate(
        Household = ifelse(grepl("TRUE", term), "With partner", 
                      "Without partner"),
        Household = add_n(Household, "with_partner", reverse = TRUE)
      ) %>% 
      standard_ggplot(fill_col = Household, n_regs = n_regs, 
                        score_col = score_name, fill_direction = -1)
```


```{r calc-age-flb-cross}

drake::loadd(res_age_flb_cross)

cor_age_flb <- res_age_flb_cross %>% 
      mutate(category = gsub(":.*", "", term)) %>% 
      pivot_wider(score_name, names_from = category, values_from = estimate) %>%              {cor(.$`age_flb_cat10-22`, .$`age_flb_cat28-52`)}

```


```{r calc-corr-age-flb}

drake::loadd(res_age_flb)
drake::loadd(res_sex)

res_sex_comparison <- res_sex %>% 
      filter(sex == "Female") %>% 
      mutate(score_name = term) %>% 
      select(score_name, term, estimate:conf.high)

res_age_flb %<>% filter(term != "age_flb") 

res_combined <- bind_rows(
        raw      = res_sex_comparison, 
        with_age = res_age_flb %>% select(score_name, term, estimate:conf.high), 
        .id      = "reg.type"
      ) %>% 
      arrange(score_name)

raw_flb_comparison <- res_combined %>% 
      pivot_wider(
        id_cols     = score_name, 
        names_from  = reg.type, 
        values_from = estimate
      )

cor_raw_flb <- cor(raw_flb_comparison$raw, raw_flb_comparison$with_age)
opp_signed <- sum(
        sign(raw_flb_comparison$raw) != sign(raw_flb_comparison$with_age)
      )

```


Lastly, we split female respondents by age at first live birth
(AFLB).[^aflb-no-men] There is evidence for genetic effects on AFLB
[@Barban_2016], and there is a close link between this variable and number of
children born. Figure \@ref(fig:plot-age-flb) shows effect sizes estimated
separately for each tercile of AFLB. Effects are strikingly different
across terciles. Educational attainment, ADHD and MDD are selected for amongst
the youngest third of mothers, but selected against among the oldest two-thirds.
Similarly, several polygenic scores for body measurements are selected against
only among older mothers. The correlation between effect sizes for the youngest
and oldest terciles is `r cor_age_flb`. To investigate this further, we estimate
equation \@ref(eq:regression) among females, *controlling* for AFLB. In 
`r opp_signed` out of `r length(score_names)` cases, effects change sign when
controls are added. The correlation between effect sizes controlling for AFLB,
and raw effect sizes, is `r cor_raw_flb`. Thus, selection effects seem to come
through two opposing channels: a correlation with AFLB, and an opposite-signed
correlation with number of children after AFLB is controlled for.

[^aflb-no-men]: AFLB is unavailable for men.


```{r plot-age-flb, fig.cap = "Selection effects by age at first live birth terciles (women only).", fig.align = "center"}

drake::loadd(res_age_flb_cross)

# this counts each cross term as a separate test:
n_regs <- as.double(nrow(res_age_flb_cross)) 

res_age_flb_cross %>% 
      mutate(
        `Age at first live birth` = gsub("age_flb_cat(.*):.*", "\\1", term),
        `Age at first live birth` = add_n(`Age at first live birth`, "age_flb_cat")
      ) %>% 
      standard_ggplot(score_col = score_name, 
                      fill_col = `Age at first live birth`, n_regs = n_regs)

```


```{r calc-siblings-children}
sib_chn_cor <- cor.test(famhist$RLRS_parents, famhist$RLRS, use = "complete")
stopifnot(sib_chn_cor$p.value < 2e-16)
```


We emphasize that these categories are not exogenous to polygenic scores. For
example -- both in the data (Appendix Figure \@ref(fig:plot-age-flb-dv)) and in
our theoretical model -- education and age at first live birth are choice
variables, which are endogenous to a person's human capital and to relevant
polygenic scores. Nevertheless, differences in selection effects across
subgroups constrain the set of possible explanations. A good theory of contemporary
natural selection needs to show how these differences come about. As we describe
below, a model based on the economic theory of fertility can do just that.

We also examine selection effects among respondents' parents, using information
on respondents' number of siblings to calculate parents' RLRS. Effect sizes of
polygenic scores are highly correlated across the two generations (Appendix
Figure \@ref(fig:plot-res-sibs)). Median-splitting respondents by year of birth,
we find little evidence of change in effect sizes among the parents' generation.
There is some evidence that selection effect sizes are increasing in the 
respondents' generation, with 8 polygenic scores showing a significant increase. 
We also check whether selection effects vary by
AFLB and socio-economic status in the parents' generation, using the 1971 
Townsend deprivation score of respondents' birthplace as a proxy for income
[@townsend1987deprivation]. Results show the same pattern as for the
respondents' generation. Effect sizes are larger and more often significant in
the most deprived areas (Appendix Figure \@ref(fig:plot-siblings-townsend)).
Effects are larger among younger fathers and mothers, and change sign when
controlling for AFLB (Appendix Figures \@ref(fig:plot-age-birth-parents-cross),
\@ref(fig:plot-age-birth-parents)). Lastly, we check for a "quantity-quality
tradeoff" between parents' number of children and number of grandchildren. We
don't find any: in fact, the correlation between respondents' and parents' RLRS
is positive ($\rho$ = `r sib_chn_cor$estimate`, $p < 2 \times 10^{-16}$).


\FloatBarrier

# Human capital and natural selection

These results show that selection effects are weaker, absent, or even reversed
among some subgroups of the population. A possible explanation for this comes 
from the economic theory of fertility 
[@becker1960economic; @willis1973new;@becker1976child]. According to this theory,
increases in a person's wage affect their fertility via two opposing channels.
There is an *income effect* by which children become more affordable, like any
other good. There is also a *substitution effect*: since childrearing has a cost
in time, the opportunity cost of childrearing increases if one's market wage is
higher. The income effect leads higher earners to have more children. The
substitution effect leads them to have fewer. 

Suppose that certain genetic variants correlate with *human capital*: skills or
other characteristics that affect an individual's earnings in the labour market
[@mincer1958investment; @becker1964human]. These variants may then be associated
with opposing effects on fertility. The income effect will lead to natural
selection in favour of earnings-increasing variants (or variants that are merely
associated with higher earnings). The substitution effect will do the reverse.

To show this, consider a simple model of fertility choices. $h$ is
an individual's level of human capital. For now, we simply identify this with
his or her wage $W$. Raising a child takes time $b$. People maximize utility $U$
from the number of children $N$ and from income $Y\equiv(1-bN)W$:
\[
U = u(Y)+aN.
\]
Here $a$ captures the strength of preference for children. $u(\cdot)$ captures
the taste for income, and is increasing and concave. We treat $N$ as continuous,
in line with the literature: this can be thought of as the expected number of
children among people with a given $a$, $b$ and $W$. The marginal
benefit of an extra child is $\frac{dU}{dN} = -bWu'(Y)+a$. The effect of an
increase in human capital on this marginal benefit is
\[
\frac{d^{2}U}{dNdW}=\underbrace{-bu'(Y)}_{\textrm{Substitution effect}}\underbrace{-bYu''(Y)}_{\textrm{Income effect}}.
\]
The *substitution effect* is negative and reflects that when wages increase,
time devoted to childcare costs more in foregone income. The positive *income
effect* depends on the curvature of the utility function, and reflects that when
income is higher, the marginal loss of income from children is less painful.

To examine education and fertility timing, we extend the model to
two periods. For convenience we ignore time discounting, and assume
that credit markets are imperfect so that agents cannot borrow. Write
\begin{equation}
U(N_{1},N_{2}) = u(Y_{1}) + u(Y_{2}) + aN_{1} + aN_{2}\label{eq:U}
\end{equation}
Instead of identifying human capital with wages, we now allow individuals to
spend time $s \in [0,1]$ on education in period 1.
Education is complementary to human capital $h > 0$, and increases period 2
wages, which take the simple functional form $w(s,h) = sh$. We normalize period
1 wages to 1, and let $u(\cdot)$ take the constant relative risk aversion form
$u(y)=\frac{y^{1 - \sigma} - 1}{1 - \sigma}$. $\sigma > 0$ measures the
curvature of the utility function, i.e. the decline in marginal utility of
income as income increases. We examine total fertility 
$N^{*} = N_{1}^{*} + N_{2}^{*}$ and the *fertility-human capital relationship*,
$\frac{dN^{*}}{dh}$. For $\sigma < 1$ and close enough to 1, Table
\@ref(tab:theory) shows five theoretical predictions, along with our
corresponding empirical results for the correlation between polygenic scores and
RLRS.[^one-period-results] The key insight of the model is that for middling
levels of $\sigma$, the substitution effect dominates at low income levels, but
as income increases, the income and substitution effect balance out.

[^one-period-results]: Predictions 1-3 also hold in the one-period model with
constant relative risk aversion. Our empirical results are actually stronger
than prediction 5, in that correlations with fertility are *reversed* at higher
AFLB. This prediction can be accommodated in the model if children have a money
cost as well as a time cost (Appendix Figure \@ref(fig:N-plot-with-m)).


------------------------------------------------------------------------------------
    **Theory: the fertility-human capital         **Empirical results**
    relationship is...**                                      
--- --------------------------------------------- ----------------------------------
1.  Negative: $\frac{dN^{*}}{dh} < 0$.            Figures \@ref(fig:plot-means-over-time) 
                                                  and \@ref(fig:plot-res-wt).
            
2.  Weaker (closer to zero) at higher wages       Figure 
    and/or levels of human capital.               \@ref(fig:plot-income-educ-level)a.
                                                  Selection effects are also weaker at
                                                  higher polygenic scores for 
                                                  educational attainment 
                                                  (Appendix Figure
                                                  \@ref(fig:plot-purifying)).
            
3.  More negative when the time burden            Stronger effects for single parents
    of children $b$ is larger.                    (Figure \@ref(fig:plot-n-partners)).
            
4.  Weaker at higher levels of education $s$.     Figure \@ref(fig:plot-income-educ-level)b.
            
5.  Weaker among those who start fertility        Effects weaker among those starting 
    in period 2 ($N_{1}^{*} = 0$) than among      fertility later 
    those who start fertility in period 1         (Figure \@ref(fig:plot-age-flb)).
    ($N_{1}^{*} > 0$). 
    
------------------------------------------------------------------------------------

Table: (\#tab:theory) Predictions from the theoretical model and corresponding 
empirical results.


Thus, a simple economic model can explain many of our results. Other empirical
work in economics also supports the link from human capital to fertility.
@caucutt2002women and @monstad2008education show that education and skills
affect age at first birth and fertility. Income decreases fertility at low
income levels, but increases it at higher income levels [@cohen2013financial].
US fertility decreases faster with education among single mothers than married
mothers [@baudin2015fertility], in line with our prediction 3 and as predicted
by @becker1981treatise. A related literature shows negative correlations between 
IQ and fertility [e.g. @lynn2004new; @reeve2018systematic].


```{r plot-cors-earnings-educ, fig.cap="Selection effects by correlations with earnings and educational attainment. Each point represents one polygenic score. Selected scores are annotated."}

drake::loadd(res_all)
drake::loadd(res_age_flb)
drake::loadd(res_cor_educ)
drake::loadd(res_cor_income)

effect_size <-  res_all %>% 
                  filter(dep.var == "RLRS", reg.type == "controlled") %>% 
                  pull(estimate)

effect_size_flb <-  res_age_flb %>% 
                      filter(term != "age_flb") %>% 
                      pull(estimate)

dfr <- data.frame(
        PGS          = rep(score_names, 4),
        Correlation  = rep(c("Earnings", "Education"), each = 33, 2),
        cor          = rep(c(res_cor_income[,1], res_cor_educ[,1]), 2),
        effect_size  = c(rep(effect_size, 2), rep(effect_size_flb, 2)),
        Controls     = rep(c("None", 
                        "Age at first live birth"), each = 66)
      )

dfr %>% 
      mutate(
        PGS      = ifelse(abs(effect_size) > 0.027, PGS, ""),
        PGS      = sub("(^\\w*?)_.*", "\\1", PGS)
      ) %>% 
      filter(
        Controls == "None"
      ) %>% 
      ggplot(aes(effect_size, cor)) + 
        geom_point(color = "steelblue4") +
        geom_text(aes(label = PGS), size = 3, colour = "black", 
                    family = "Abadi MT Condensed Light",
                    check_overlap = TRUE, nudge_x = 0.004) +
        facet_wrap( ~ Correlation, scales = "free") + 
        my_vline + 
        my_hline +
        labs(x = "Selection effect", y = "Correlation") +
        coord_cartesian(clip = "off") +
        theme(
          legend.position = "none",
          panel.spacing = unit(20, "pt")
        )

```


# Testing the theory

We test the economic theory in two ways. First, it predicts that
genetic variants will be selected for (or against) in proportion to their
correlation with human capital. Figure \@ref(fig:plot-cors-earnings-educ) plots
selection effects on each polygenic score against that score's correlation with
two measures of human capital: earnings in a respondent's first job, and
educational attainment. The relationships are strongly negative. Thus, human
capital appears to be relevant to natural selection. The negative relationship
suggests that substitution effects dominate income effects, which fits the known
negative association between income and fertility
[@becker1960economic;@jones2006economic]. The correlations reverse when we
control for age at first live birth, suggesting that within AFLB categories, the
income effect dominates. 

```{r calc-mediation}
drake::loadd(res_mediation)

res_mediation %<>% mutate(
                     # 0.025 not 0.05 because two-sided
                     sig     = abs(statistic_ind) > qnorm(1 - 0.025/
                                                            nrow(res_mediation)),
                     pos_sig = sig & prop_ind > 0,
                     bounded_ci = sign(estimate_total_conf_low) ==
                                  sign(estimate_total_conf_high)
                   )


n_pos_sig <- sum(res_mediation$pos_sig)
median_prop <- median(res_mediation$prop_ind[res_mediation$pos_sig])

```


Second, we run a mediation analysis to directly test whether the correlation
between each polygenic score and fertility is mediated by educational attainment
(Appendix Table \@ref(tab:tbl-res-mediation)). We use the 
`r nrow(res_mediation)` scores where the selection effect is significant at $p$ <
0.05/33. Figure \@ref(fig:plot-mediation) shows estimated
proportions explained by educational attainment, along with bootstrap 95%
confidence intervals (uncorrected; 100 bootstraps). For `r n_pos_sig` scores,
the indirect effect of the score on fertility *via* educational attainment takes
the same sign as the overall effect, and is significantly different from zero
($p$ < 0.05/`r nrow(res_mediation)`). Among these scores, the median proportion of
the total effect explained by the indirect effect is 
`r scales::percent(median_prop)`. The educational attainment variable is a
relatively crude measure of human capital: more accurate measures would likely
explain more of the total effect.


(ref:franzcite) [@franz2007ratios]

```{r plot-mediation, fig.align = "center", fig.cap = "Proportion of selection effect mediated by educational attainment, among polygenic scores with significant selection effects. Bootstrap confidence intervals for the proportion are shown only where the interval is bounded (ref:franzcite)."}

res_mediation <- res_mediation %>% 
      mutate(
        term = pretty_names(term),
        term = fct_reorder(term, prop_ind, order_abs(1)),
        `Indirect effect p < 0.05/23` = sig 
      )

res_mediation %>% 
      ggplot(aes(y = term)) + 
        geom_col(aes(x = prop_ind, fill = `Indirect effect p < 0.05/23`), 
                   width = .5) + 
        geom_errorbar(aes(
                          xmin = prop_ind_conf_low, 
                          xmax = prop_ind_conf_high,
                          linetype = "Bootstrap 95% c.i. uncorrected",
                        ),
                        data = res_mediation %>% filter(bounded_ci),
                        width = 0.2,
                        color = "grey20"
                      ) +
        scale_fill_manual(values = c("FALSE" = "grey65", "TRUE" = "steelblue4")) +
        my_vline +
        scale_x_continuous(labels = scales::percent, n.breaks = 7) +
        coord_cartesian(xlim = c(-0.5, 1)) +
        labs(x = "Proportion mediated", y = "") +
        scale_linetype_manual(name = "", 
          values = c("Bootstrap 95% c.i. uncorrected" = 1)) +
        my_theme(
          panel.grid.major.y = element_blank(),
          panel.grid.minor.x = element_blank(),
          legend.justification = "top"
        )


```


```{r calc-risk-attitude}

drake::loadd(res_risk_control)
drake::loadd(res_all)

n_risk_sig <- res_risk_control %>% 
                filter(term == 'f.2040.0.0') %>% 
                {sum(.$p.value <= .Machine$double.eps)}

n_pgs_sig <- res_risk_control %>% 
               filter(term == score_name) %>% 
               {sum(.$p.value <= 0.05/33)}

res_compare <- res_all %>% 
      filter(dep.var == "RLRS", reg.type == "controlled") %>% 
      left_join(res_risk_control %>% filter(term == score_name), 
                  by = "score_name")

ratio <- res_compare$estimate.y/res_compare$estimate.x
n_less_sig <- sum(res_compare$p.value.y <= 0.05/33 & res_compare$p.value.x > 0.05/33)
stopifnot(n_less_sig == 0)
```

We consider three alternative theories that might explain our results. First,
welfare benefits which incentivize child-bearing might be taken up more among
low-income people. However, the majority of effect sizes appear unchanged
over a large span of twentieth-century history (Appendix Table
\@ref(tab:tbl-change-by-period)), during
which government spending on child-related benefits varied considerably
[@hoc1999socsec]. In general, there is only weak evidence that welfare benefits
affect fertility [@gauthier2007impact; see also @bergsvik2021can]. Future
work could test this theory more explicitly. A second alternative theory is that
polygenic scores correlate with the motivation to have children, i.e. parameter
$a$ in the model [cf. @jones2008fertility]. This theory would not explain why
selection effects are smaller at higher incomes and education levels. In fact,
in the model, $a$'s effect on fertility gets stronger at higher levels
of human capital. A third alternative is that traits under selection are linked
to externalizing behaviour and risk-seeking. This might be partially captured by
our parameter $\sigma$, which can be interpreted as a measure of risk aversion
over income; a more direct channel is risky sexual behaviour
[@mills2021identification]. The data here provide some support for this story: scores
which might plausibly be linked to externalizing behaviour, like ADHD and
younger age at smoking initiation, are selected for. However, risk-seeking 
seems unlikely to explain variation in fertility across the full range of scores
under selection, including physical measures like waist-hip ratio and BMI. We
test this theory directly by re-estimating equation \@ref(eq:regression)
controlling for a measure of risk attitude (UK Biobank field 2040). The median
ratio of effect sizes between regressions with and without controls is 
`r median(ratio)`; all scores which are significant at $p < 0.05/33$ in
uncontrolled regressions remain so when controlling for risk attitude. This
non-result could simply reflect the imprecision of the risk attitude measure,
which is a single yes/no question. But this measure *does* predict the overall
number of children, highly significantly ($p < 2 \times 10^{-16}$ in 
`r n_risk_sig` out of 33 regressions). Given that, and the statistical power we
get from our sample size, we believe that the non-result is real: while risk
attitude does predict fertility in the sample, it is not an important channel
for natural selection.


# Discussion

Previous work has documented natural selection in modern populations on variants
underlying polygenic traits [@beauchamp2016genetic; @kong2017selection;
@sanjak2018evidence]. We show that correlations between polygenic scores and
fertility are highly concentrated among specific subgroups of the population,
including people with lower income, lower education, younger first parenthood,
and more lifetime sexual partners. Among mothers aged 22+, selection effects
are reversed. Furthermore, the size of selection effects on a polygenic score
correlates with that score's association with labour market earnings. Strikingly,
some of these results were predicted by @fisher1930genetical, pp. 253-254. The
economic theory of fertility gives a parsimonious explanation for these
findings. Because of the substitution effect of earnings on fertility, scores
are selected for when they correlate with low human capital, and this effect is
stronger at lower levels of income and education.

Polygenic scores which correlate with lower earnings and less education are
being selected for. In addition, many of the phenotypes under positive selection
are linked to disease risk. Many people would probably prefer to have high
educational attainment, a low risk of ADHD and major depressive disorder, and a
low risk of coronary artery disease, but natural selection is pushing against
genes associated with these traits. Potentially, this could increase the health
burden on modern populations, but that depends on effect sizes. Our results show
that naïve estimates can be affected by sample ascertainment bias. There may be
remaining sources of ascertainment bias after our weighting; if so, we expect
that, like the sources of ascertainment we have controlled for, they probably
bias our results towards zero. Researchers should be aware of the risks of
ascertainment when studying modern natural selection.

We also do not know how estimated effect sizes of natural selection will change
as more accurate polygenic scores are produced, or whether genetic
variants underlying other phenotypes will show a similar pattern to those
studied here. Also, effects of polygenic scores may be inflated in
population-based samples, because of indirect genetic effects, gene-environment
correlations, and/or assortative mating [@lee2018gene;
@selzam2019comparing; @kong2018nature; @Howe2021.03.05.433935], although we do
not expect that this should change their association with number of offspring,
or the resulting changes in allele frequencies. Although effects on
our measured polygenic scores are small even after weighting, individually small
disadvantages can cumulate to create larger effects. Lastly, note that our
data comes from people born before 1970. Recent evidence suggests that fertility
patterns may be changing [@doepke2022economics]. Overall, it is probably 
too early to tell whether modern natural selection has a substantively important 
effect on population averages of phenotypes under selection.


```{r calc-inequality}
drake::loadd(res_ineq)

n_increases <- sum(res_ineq$ratio > 1)
median_pct_change <- (median(res_ineq$ratio) - 1) * 100

```


Because selection effects are concentrated in lower-income groups, they may also
increase inequality with respect to polygenic scores. For example, Figure
\@ref(fig:plot-mean-EA3-by-income) plots mean polygenic scores for educational
attainment (EA3) among children from households of different income groups. The
blue bars show the actual means, i.e. parents' mean polygenic score weighted by
number of children. The grey bars show the hypothetical means if all households
had equal numbers of children. Natural selection against genes associated with
educational attainment is stronger at the bottom of the income distribution, and
this increases the differences between groups. Overall, natural selection
increases the correlation of polygenic scores with income for `r n_increases`
out of 33 polygenic scores, with a median percentage increase of 
`r median_pct_change`% in the respondents' generation (Appendix Table
\@ref(tab:tbl-inequality-effects)).  If inequalities in polygenic scores are
important for understanding social structure and mobility [@belsky2018genetic;
@Rimfeld_2018; @harden2021genetic], then these increases are substantive.
Similarly, since many polygenic scores are predictive of disease risk, they
could potentially increase health inequalities. In general, the evolutionary
history of anatomically modern humans is related to disease risk
[@benton2021influence]; understanding the role of contemporary natural selection
may help researchers to map the genetic architecture of health disparities.


```{r plot-mean-EA3-by-income, fig.align = "center", fig.cap = "Mean polygenic score for educational attainment (EA3) of children by household income group. Blue is actual. Grey is hypothetical in the absence of selection effects.", fig.width = 4}
income_EA3 <- famhist %>% 
      filter(! is.na(n_children), ! is.na(income_cat)) %>% 
      mutate(
        `Household income` = factor(income_cat, 
                          labels = c("< £18K", "£18-30K", "£31-51K", "£52-100K", 
                                      "> £100K"))
      ) %>% 
      group_by(`Household income`) %>% 
      summarize(
        `Without selection` = mean(EA3_excl_23andMe_UK, na.rm = TRUE),
        Actual              = weighted.mean(EA3_excl_23andMe_UK, n_children, 
                                              na.rm = TRUE)
      ) %>% 
      tidyr::pivot_longer(-`Household income`) %>% 
      rename(
        `Children's mean EA3` = value
      )

income_EA3 %>% 
              ggplot(aes(`Household income`, `Children's mean EA3`, fill = name)) +
                geom_col(position = "dodge") +
                scale_fill_manual(values = c("steelblue4", "grey65")) +
                theme(
                  legend.background = element_rect(
                                        fill = "white", 
                                        colour = "black",
                                        size = rel(0.5)
                                      ),
                  legend.margin = margin(t = 0, b = 5.5, l = 5.5, 
                                           r = 5.5),
                  legend.title = element_blank(), 
                  legend.position = c(.25, .85)
                )
```


Existing evidence on human natural selection has led some to "biocosmic
pessimism" [@sarraf2019modernity]. Others are more sanguine, and argue that
natural selection's effects are outweighed by environmental improvements, like
those underlying the Flynn effect [@flynn1987massive]. The evidence
here may add some nuance to this debate. Patterns of natural selection have been
relatively consistent across the past two generations, but they are not the
outcome of a single, society-wide phenomenon. Instead they result from 
opposing forces, operating in different parts of society and pulling in
different directions.

Any model of fertility is implicitly a model of natural selection, but so far,
the economic and human genetics literatures have developed in parallel.
Integrating the two could deepen our understanding of natural selection in
modern societies. Economics possesses a range of theoretical models on the
effects of skills, education and income [see @hotz1997economics;
@lundberg2007american]. One perennial problem is how to test these theories in a
world where education, labour and marriage markets all interact. Genetic data,
such as polygenic scores, could help to pin down the direction of causality,
for example via Mendelian randomization [@davey2003mendelian].
Conversely, economic theories and empirical results can shine a light on the
mechanisms behind natural selection, and thereby on the nature of individual
differences in complex traits and disease risk.


# Materials and methods

We use participant data from UK Biobank [@bycroft2018uk], which has received
ethical approval from the National Health Service North West Centre for Research
Ethics Committee (reference: 11/NW/0382). We limit the sample to white British
participants of European descent, as defined by genetic estimated ancestry and
self-identified ethnic group, giving a sample size of
`r comma_num(nrow(famhist))`. For regressions on number of children we use
participants over 50 (males)/45 (females), since most fertility is completed by this age. This gives
a sample size of `r comma_num(nrow(famhist %>% filter(kids_ss)))`.

Polygenic scores were chosen so as to cover a reasonably broad range of traits,
and based on the availability of a large and powerful GWAS which did not include
UK Biobank. Scores were computed by summing the alleles across ~1.3 million
genetic variants weighted by their effect sizes as estimated in 33 genome-wide
association studies (GWASs) that excluded UK Biobank. To control for population
stratification, we corrected the polygenic scores for 100 principal components
(PCs). To compute polygenic scores and PCs, the same procedures were followed as 
described in @abdellaoui2019genetic.

Earnings in first job are estimated from mean earnings in the 2007 Annual Survey
of Hours and Earnings, using the SOC 2000 job code (Biobank field 22617).

Weighting data was kindly provided by @vanalten2022reweighting.

Code for this paper is available at 
https://github.com/hughjonesd/why-natural-selection.

## Acknowledgements 

AA is supported by the Foundation Volksbond Rotterdam and by ZonMw grant 849200011 from The Netherlands Organisation for Health Research and Development. This study was conducted using UK Biobank resources under application numbers 40310 and 19127.


\FloatBarrier
\clearpage
# Appendix

\localtableofcontents
\clearpage

## Selection effects by sex

Figure \@ref(fig:plot-sexes) plots selection effects by sex. Differences are
particularly large for educational attainment, height, ADHD and MDD. Several
polygenic scores for mental illness and personality traits are more selected for
(or less against) among women, including major depressive disorder (MDD),
schizophrenia and neuroticism, while extraversion is more selected for among
men.

```{r plot-sexes, fig.cap = "Selection effects by sex. Solid lines are significant differences at p < 0.05/33. Solid points are significantly different from 0 at p < 0.05/66."}

drake::loadd(res_sex)
n_regs <- as.double(nrow(res_sex))


res_sex_wide <- res_sex %>% 
      select(term, estimate, sex, diff.p.value) %>% 
      pivot_wider(names_from = "sex", values_from = "estimate") %>% 
      mutate(
        min_est = pmin(Female, Male),
        max_est = pmax(Female, Male),
        term = pretty_names(term),
        term = fct_reorder(term, Female),
        p = diff.p.value
      ) 
n_comps <-as.double(nrow(res_sex_wide)) 

res_sex %>% 
      standard_ggplot(fill_col = sex, n_regs = n_regs, score_col = term,
                        conf_int = FALSE, fill_direction = -1) +
      geom_linerange(data = res_sex_wide, inherit.aes = FALSE, 
                     color = "steelblue4", alpha = 1,
              aes(
                y      = term,
                x      = NULL,
                xmin   = min_est + 0.0005, 
                xmax   = max_est - 0.0005,
                linetype = p < 0.05/{{n_comps}}
              )) +
      scale_linetype_manual(values = c("TRUE" = 1, "FALSE" = 3))
      

      

```


\FloatBarrier
\clearpage
## Alternative weighting schemes

We compare results for our main weights to 3 alternative weighting schemes:
weighting by age/qualification; geographical (weighting by Middle Super Output
Area); and for women only, age, qualification and age at first live birth.
Population data for weighting is taken from the 2011 UK Census and the 2006
General Household Survey (GHS). Weighting for Age/Qualification and
Age/Qualification/AFLB weights was done using marginal totals from a linear
model, using the `calibrate()` function in the R "survey" package
[@lumley2020survey]. Geographical weighting was done with iterative
post-stratification using the `rake()` function, on Census Middle Layer Super
Output Areas, sex and presence/absence of a partner.

Table \@ref(tab:tbl-res-wt) gives effect sizes as a proportion of the unweighted
effect size, for all polygenic scores which are consistently signed and which
are significantly different from zero in unweighted regressions.


```{r tbl-res-wt}

rel_sizes %>% 
      arrange(desc(Main)) %>% 
      bind_rows(rel_sizes_summary) %>% 
      rename(PGS = term) %>% 
      as_hux() %>% 
      insert_row("", "Weighting", fill = "") %>% 
      set_caption(
        "Weighted effect sizes as a proportion of unweighted effect sizes."
      ) %>% 
      merge_cells(1, 2:4) %>% 
      set_align(1, 2, "centre") %>% 
      set_header_rows(1, TRUE) %>% 
      theme_compact() %>% 
      set_number_format(-1, -1, 2) %>% 
      style_headers(bold = TRUE) %>% 
      set_bottom_border(1, -1) %>% 
      set_bottom_border(2, everywhere) %>% 
      set_italic(final(2), 1) %>% 
      set_font_size(8) %>% 
      set_width(0.8) %>% # ensures line wrapping below
      set_col_width(c(.28, .18, .18, .18, .18)) %>% 
      add_footnote(paste0(
        "Only consistently-signed and significant (when unweighted) ",
        "estimates are shown. ",
        "Age/Qual/AFLB as a proportion of unweighted regressions ",
        "including females only."),
        top_padding  = 0,
        left_padding = 0
      )
      
```


\FloatBarrier
\clearpage
## Stabillizing and disruptive selection

Stabilizing selection reduces variance in the trait under selection, while
disruptive selection increases variance. To check for these, we rerun equation
\@ref(eq:regression), adding a quadratic term in $PGS_i$. Scores for hip
circumference show significant stabilizing selection ($p < 0.05/33$, negative
coefficient on quadratic term). The EA2 score for educational attainment shows
significant disruptive selection ($p < 0.05/33$, positive coefficient), which
reduces the strength of selection against educational attainment at very high
levels of the PGS. (The quadratic on the EA3 score has a similar coefficient but
is not significant at $p < 0.05/33$.) Figure \@ref(fig:plot-purifying) plots
predicted number of children against polygenic score from these regressions.

We also checked for stabilizing selection in the parents' generation, using
weights multiplied by the inverse of number of siblings. Scores for EA2 and EA3
show significant disruptive selection ($p < 0.05/33$, positive coefficient on
quadratic). Other scores including hip circumference were not significant.


```{r plot-purifying, fig.align = "center", fig.cap = "Stabilizing/disruptive selection: predicted number of children by polygenic score."}

drake::loadd(res_quadratic)
sig <- res_quadratic %>% 
         filter(grepl("^2", term, fixed = TRUE)) %>% 
         filter(p.value < 0.05/33) %>% 
         pull(score_name)

stopifnot(sig == c("ADHD_2017", "bmi_combined", "EA2_noUKB", 
                     "hip_combined", "wc_combined"))

sig_all_terms <- res_quadratic %>% filter(score_name %in% sig)


drake::loadd(res_sibs_quadratic)
sig_sibs <- res_sibs_quadratic %>% 
         filter(grepl("^2", term, fixed = TRUE)) %>% 
         filter(p.value < 0.05/33) %>% 
         pull(score_name)
stopifnot(sig_sibs == c("EA2_noUKB", "EA3_excl_23andMe_UK"))

op <- par(mfrow = c(2, 3))

for (score in sig) {
  coef <- sig_all_terms %>% filter(score_name == score) %>% pull(estimate)
  curve(
         coef[1] + coef[2] * x + coef[3] * x^2,
         xlim = c(-4, 4), 
         xlab = pretty_names(score),
         ylab = "Predicted RLRS"
       )
}
par(op)

```


\FloatBarrier
\clearpage
## Controlling for age {#sec:age-control}

```{r calc-income-edu-controlled}

drake::loadd(res_income_controlled, res_edu_controlled)
n_regs_ic <- nrow(res_income_controlled)
n_sig_ic <- sum(res_income_controlled$p.value < 0.05/n_regs_ic)

n_regs_educ <- nrow(res_edu_controlled)
n_sig_educ <- sum(res_edu_controlled$p.value < 0.05/n_regs_educ)
```

Results in Figure \@ref(fig:plot-income-educ-level) could be
explained by age, if older respondents have lower income and are less educated,
and also show more natural selection on polygenic scores. However, when we rerun
the regressions, interacting the polygenic score with income category and also
with a quadratic in age, the interaction with income remains significant at
$p$ < 0.05/`r n_regs_ic` for `r n_sig_ic` out of `r n_regs_ic` regressions. Similarly
if we interact the PGS with age of leaving full time education and a quadratic
in age, the interaction with age leaving full time education remains significant
at $p$ < 0.05/`r n_regs_educ` for `r n_sig_educ` out of `r n_regs_educ` regressions.


\FloatBarrier
\clearpage
## Number of partners and presence of partner by sex

Figure \@ref(fig:plot-n-partners-sex) splits up Figure \@ref(fig:plot-n-partners)
by sex. The pattern of results is the same in both sexes: selection effects are stronger among those with more lifetime sexual partners, and among those not
currently living with a partner.


```{r plot-n-partners-sex, fig.cap = "Selection effects by number of sexual partners and presence of a partner, for men and women separately.", fig.subcap = c("Lifetime number of sexual partners", "Presence of a partner"), fig.ncol = 1, fig.align = "center", fig.height = 4}

sexes_theme <-  my_theme(
                  plot.margin          = unit(c(1, 1, 0, 1), "lines"),
                  legend.position      = "bottom", 
                  #legend.box.spacing   = unit(0, "lines"),
                  legend.box.margin    = margin(-5, 0, 0, -85),
                  legend.margin        = margin(0, 0, 0, 0),
                  legend.box.just      = "left",
                  legend.justification = "left",
                  legend.spacing.y     = unit(0, "lines"),
                  #legend.spacing       = unit(0.5, "lines"),
                  legend.title         = element_text(size = 10)
                  #legend.text          = element_text(size = 8)
                )
sexes_guides <- guides(
                  fill  = guide_legend(ncol = 1),
                  color = guide_legend(ncol = 1)
                )

drake::loadd(res_partners)
res_partners %<>%       
      filter(grepl(":", term)) %>% 
      mutate(`N partners` = ifelse(grepl("TRUE:", term), "3 or less", "4 or more")) 


n_regs <- as.double(nrow(res_partners))
res_partners %>% 
      standard_ggplot(fill_col = `N partners`, score_col = score_name, 
                        n_regs = n_regs, fill_direction = -1) +
        facet_grid(cols = vars(sex))

drake::loadd(res_with_partner_sex)
res_with_partner_sex %<>% filter(grepl(":", term)) 

n_regs <- as.double(nrow(res_with_partner_sex))

res_with_partner_sex %>% 
      mutate(
        Household = ifelse(grepl("TRUE", term), "With partner", 
                      "Without partner")
      ) %>% 
      standard_ggplot(fill_col = Household, score_col = score_name, 
                        n_regs = n_regs, fill_direction = -1) +
        facet_grid(cols = vars(sex)) 
  
```


\FloatBarrier
\clearpage
## Parents' generation {#sec-parents}
### Selection effects and change over time

```{r calc-res-sibs}
drake::loadd(res_wt_van_alten_weights)
drake::loadd(res_children_comparison)

cmp <- left_join(res_wt_van_alten_weights, res_children_comparison, by = "term")

cor_inc_vs_excl_childless <- cor(cmp$estimate.x, cmp$estimate.y)
n_smaller <- sum(abs(cmp$estimate.y) < abs(cmp$estimate.x))

pct_diff <- abs(cmp$estimate.y)/abs(cmp$estimate.x) - 1
median_pct_diff <- median(pct_diff)

```


The UK Biobank data contains information on respondents' number of siblings
(including them), i.e. their parents' number of children. Since respondents'
polygenic scores are equal in expectation to the mean scores of their parents,
we can use this to look at selection effects in the parents' generation. We
estimate equation (\ref{eq:regression}) using parents' RLRS as the dependent
variable.[^parent-RLRS] The parents' generation has an additional source of ascertainment
bias: sampling parents of respondents overweights parents who have many
children. For instance, parents of three children will have, on average, three
times more children represented in UK Biobank than parents of one child. Parents
of no children will by definition not be represented. To compensate, we multiply
our weights by the inverse of *number of siblings*.

[^parent-RLRS]: We don't have data on parents' year of birth for most 
respondents. To create parents' RLRS, we divide respondents' number of siblings 
by the average number of siblings of all respondents born in the same year,
weighting the average by respondents' inverse of number of siblings to 
compensate for ascertainment bias.

Figure \@ref(fig:plot-res-sibs) shows regressions of parents' RLRS on
polygenic scores. For a clean comparison with the respondents' generation, we
rerun regressions on respondents' RLRS excluding those with no children, and
show results in the figure. Selection effects are highly correlated across the
two generations, and most share the same sign. Absolute effect size estimates
are larger for the parents' generation. We treat this result cautiously, because
effect sizes in both generations may depend on polygenic scores' correlation with
childlessness, and we cannot estimate this for the parents' generation.

To learn more about this, we compare effect sizes excluding and including
childless people in the *current* generation. The correlation between the two
sets of effect sizes is `r cor_inc_vs_excl_childless`. So, patterns across
different scores are broadly similar whether the childless are counted or not.
However, absolute effect sizes are smaller when the childless are excluded, for
`r n_smaller` out of 33 scores; the median percentage change is 
`r scales::percent(median_pct_diff)`.


```{r plot-res-sibs, fig.cap = "Selection effects, respondents' parents vs. respondents. Parental generation weights multiplied by 1/number of siblings. Respondents' regression excludes childless respondents."}
drake::loadd(res_sibs_parent_weights)
drake::loadd(res_children_comparison)

res_sibs_kids <-  bind_rows(
                    Parents = res_sibs_parent_weights, 
                    Respondents = res_children_comparison,
                    .id = "Generation"
                  )

standard_ggplot(res_sibs_kids, Generation, n_regs = length(score_names),
                  fill_direction = -1)
```



```{r calc-period-regs}

drake::loadd(res_period_parents)
drake::loadd(res_period_children)

res_period <- bind_rows(res_period_parents, res_period_children)

res_period %<>% 
      tidyr::extract(term, c("period", "score_name"),
        regex="year_split(.*):(.*)")
n_tests <- nrow(res_period)/2


summary_period <- res_period %>% 
      select(-std.error, -statistic) %>% 
      pivot_wider(
        names_from  = period, 
        values_from = c(estimate, conf.low, conf.high, p.value)
      ) %>% 
      group_by(score_name, children) %>% 
      mutate(
        significant = diff.p.value < 0.05/n_tests,
        sign_early  = sign(estimate_early),
        sign_late   = sign(estimate_late),
        size_inc    = abs(estimate_early) < abs(estimate_late),
        change_sign = sign_early != sign_late,
        Change      = case_when(
                        ! significant ~ "Insignificant",
                        change_sign   ~ "Change sign",
                        size_inc      ~ "Size increasing",
                        ! size_inc    ~ "Size decreasing"
                      )
      ) %>% 
      select(score_name, children, Change, significant, change_sign, 
        sign_early, sign_late) %>% 
      mutate(children = ifelse(children, "Respondents", "Parents")) %>% 
      pivot_wider(names_from = children, 
        values_from = c(Change, significant, change_sign, sign_early, 
        sign_late)) %>% 
      mutate(overall = case_when(
              ! change_sign_Parents & ! change_sign_Respondents & 
                sign_late_Parents == sign_early_Parents ~ 
                "Consistent",
              TRUE ~ "Changes direction"
            ))


stopifnot(sum(summary_period$Change_Parents == "Insignificant") == 32)
stopifnot(sum(summary_period$Change_Respondents == "Size increasing") == 8)
stopifnot(sum(summary_period$Change_Respondents == "Insignificant") == 25)
```

The fact that childless people have such a strong effect on estimates makes it
hard to compare total effect sizes across generations. In particular, since the
parents' generation has a different distribution of numbers of children,
childless people may have had more or less effect in that generation. Another
issue is that we are estimating parents' polygenic scores by the
scores of their children. This introduces noise into our independent
variable, which might lead to errors-in-variables and bias coefficients towards
zero.

As an alternative approach, we run regressions interacting polygenic scores with
birth year, median split at 1950 ("early born" versus "late born"). We use both
respondents' RLRS and parents' RLRS as a dependent variable. We use our standard
weights, and further adjust for selection in the parents' generation (see above).

Table \@ref(tab:tbl-change-by-period) summarizes the results. We report the
number of scores showing significant changes over time (i.e. a significant
interaction between polygenic score and the "late born" dummy): either a
significant change in sign, a significant increase in effect size, or a
significant decrease in size. There is little evidence for changes in selection
effects within the parents' generation, with just one score showing a
significant decrease in size. In the respondents' generation, effect sizes were
significantly larger in absolute size among the later-born for eight polygenic
scores: ADHD, age at menopause, cognitive ability, Coronary Artery Disease, EA2,
EA3, extraversion and Major Depressive Disorder. These changes are inconsistent
with the intergenerational change, where estimated effect sizes were larger
among the earlier, parents' generation.

Overall, while there is some suggestive evidence for an increase in the strength
of selection in recent history, the clearest result is that the pattern of
relative effect sizes across scores is broadly consistent over time.


```{r tbl-change-by-period}

tbl_resp <- summary_period %>% 
              group_by(Change_Respondents) %>% 
              count()

tbl_par <- summary_period %>% 
             group_by(Change_Parents) %>% 
             count()

tbl_both <- full_join(tbl_par, tbl_resp, 
                        by = c("Change_Parents" = "Change_Respondents"))

tbl_both %>% 
            rename(
              Change              = Change_Parents, 
              `Parents' RLRS`     = n.x, 
              `Respondents' RLRS` = n.y
            ) %>% 
            as_hux() %>%
            set_bold(1, everywhere) %>% 
            set_bottom_border(1, everywhere) %>% 
            add_footnote(glue("Significance is measured at p < 0.05/{n_tests}."), 
              border = 0.4, font_size = 8) %>% 
            set_top_padding(1) %>% 
            set_bottom_padding(1) %>%
            set_left_padding(everywhere, 1, 0) %>% 
            set_right_padding(everywhere, 3, 0) %>% 
            set_caption("Numbers of polygenic scores showing changes in selection effects between early and late born. Parental generation weights multiplied by 1/number of siblings.") %>% 
            set_caption_width(1)
        
  
      
```


\clearpage
\FloatBarrier
### Area deprivation


Figure \@ref(fig:plot-siblings-townsend) plots effects on parents' RLRS by
Townsend deprivation quintile of birth area.

```{r plot-siblings-townsend, fig.cap = "Selection effects (parents' RLRS) by Townsend deprivation quintile of birth area. Higher = more deprived. Weights multiplied by 1/number of siblings.", fig.align = "center"}

drake::loadd(res_townsend_parents)
n_regs <- nrow(res_townsend_parents)

res_townsend_parents %>% 
    mutate(
        n = prettyNum(n, big.mark = ","),
        `Townsend quintile` = paste0(quintile, " (N = ", n, ")")
      ) %>% 
  standard_ggplot(fill_col    = `Townsend quintile`,
                    n_regs    = n_regs, 
                    score_col = score_name, 
                    order_idx = 5,
                    fill_direction = -1
                  )

```

For comparison, Figure \@ref(fig:plot-townsend) plots effects on 
respondents' RLRS by Townsend deprivation quintile of birth area.

```{r plot-townsend, fig.cap = "Selection effects in the respondents' generation by Townsend deprivation quintile of birth area. Higher = more deprived.", fig.align = "center"}

drake::loadd(res_townsend)
n_regs <- nrow(res_townsend)

res_townsend %>% 
    mutate(
        n = prettyNum(n, big.mark = ","),
        `Townsend quintile` = paste0(quintile, " (N = ", n, ")")
      ) %>% 
  standard_ggplot(fill_col    = `Townsend quintile`,
                    n_regs    = n_regs, 
                    score_col = score_name, 
                    order_idx = 5,
                    fill_direction = -1
                  )

```


\clearpage
\FloatBarrier
### Age at first live birth


```{r calc-age-birth-parents}

drake::loadd(res_age_birth_parents)
drake::loadd(res_all)

res_sibs_unc <- filter(res_all, dep.var == "RLRS_parents", reg.type == "controlled")
compare_abp_unc <- res_age_birth_parents %>% 
      filter(! term %in% c("fath_age_birth", "moth_age_birth")) %>% 
      left_join(res_sibs_unc, by = "term")

cor_unc_fath <- with(
  compare_abp_unc %>% filter(control == "fath_age_birth"), 
  cor(estimate.x, estimate.y))
cor_unc_moth <- with(
  compare_abp_unc %>% filter(control == "moth_age_birth"), 
  cor(estimate.x, estimate.y))
```

Among the parents' generation, we can control for age at first live birth using
the subsets of respondents who reported their mother's or father's age, and who
had no elder siblings. We run regressions on parents' RLRS on these
subsets. Figure \@ref(fig:plot-age-birth-parents-cross) shows selection effects
by terciles of age at first live birth, for mothers and fathers. As in the
respondents' generation, effect sizes are smaller, or even oppositely signed, for older parents. Importantly, this holds for both sexes.

```{r plot-age-birth-parents-cross, fig.cap = "Selection effects (parents' RLRS) among eldest siblings, by parents' age at first live birth terciles.  Weights multiplied by 1/number of siblings.", fig.subcap = c("Mothers", "Fathers"), fig.ncol = 1}

drake::loadd(res_age_flb_mothers_cross)
drake::loadd(res_age_flb_fathers_cross)

famhist_bo1 <- famhist %>% 
                 filter(birth_order == 1) %>%
                 select(RLRS_parents, whr_combined, moth_age_birth_cat, 
                          fath_age_birth_cat)

# this counts each cross term as a separate test:
n_regs <- as.double(nrow(res_age_flb_mothers_cross)) 

res_age_flb_mothers_cross %>% 
      mutate(
        `Age at first live birth` = gsub("moth_age_birth_cat(.*):.*", "\\1", term),
        `Age at first live birth` = add_n(`Age at first live birth`, 
                                            "moth_age_birth_cat", 
                                            dv   = "RLRS_parents",
                                            data = famhist_bo1
                                          )
      ) %>% 
      standard_ggplot(score_col = score_name, 
                      fill_col = `Age at first live birth`, n_regs = n_regs)



# this counts each cross term as a separate test:
n_regs <- as.double(nrow(res_age_flb_fathers_cross)) 

res_age_flb_fathers_cross %>% 
      mutate(
        `Age at first live birth` = gsub("fath_age_birth_cat(.*):.*", "\\1", term),
        `Age at first live birth` = add_n(`Age at first live birth`, 
                                            "fath_age_birth_cat", 
                                            dv   = "RLRS_parents",
                                            data = famhist_bo1
                                          )
      ) %>% 
      standard_ggplot(score_col = score_name, 
                      fill_col = `Age at first live birth`, n_regs = n_regs)

rm(famhist_bo1)

```

Figure \@ref(fig:plot-age-birth-parents) shows the regressions controlling for
either parent's age at their birth. Effect sizes are very similar, whether
controlling for father's or mother's age. As in the respondents' generation,
effect sizes are negatively correlated with the effect sizes from bivariate
regressions without the control for age at birth (father's age at birth: $\rho$ 
`r cor_unc_fath`; mother's age at birth: $\rho$ `r cor_unc_moth`).

```{r plot-age-birth-parents, fig.cap = "Selection effects (parents' RLRS) among eldest siblings, controlling for parents' age at birth. Weights multiplied by 1/number of siblings."}


n_regs <- as.double(nrow(res_age_birth_parents)/2)

res_age_birth_parents %>% 
    filter(! term %in% c("fath_age_birth", "moth_age_birth")) %>% 
    mutate(
      Control = ifelse(control == "fath_age_birth", 
                  "Father's age at birth", 
                  "Mother's age at birth"
                )
    ) %>% 
    standard_ggplot(fill_col = Control, n_regs = n_regs, order_idx = 2,
                      fill_direction = -1)
    

```


\clearpage
\FloatBarrier
## Effects of polygenic scores on age at first live birth


```{r calc-age-flb-dv}

drake::loadd(res_age_flb_dv)
drake::loadd(res_age_birth_parents_dv)

compare_flb_dv <- left_join(
        res_age_flb_dv           %>% select(term, estimate),
        res_age_birth_parents_dv %>% select(term, estimate, dep.var),
        by = "term"
      )

cor_flb_fath <- with(compare_flb_dv %>% filter(dep.var == "fath_age_birth"),
        cor(estimate.x, estimate.y)
      )
cor_flb_moth <- with(compare_flb_dv %>% filter(dep.var == "moth_age_birth"),
        cor(estimate.x, estimate.y)
      )
```

Our results suggest that polygenic scores may directly correlate with age at
first live birth. Figure \@ref(fig:plot-age-flb-dv) plots estimated effect sizes
from bivariate regressions for respondents. Figure
\@ref(fig:plot-age-birth-parents-dv) does the same for their parents, using only
eldest siblings.[^only-eldest] Effect sizes are reasonably large. They are also
highly correlated across generations. Effect sizes of polygenic scores on
father's age at own birth, and on own age at first live birth, have a
correlation of `r cor_flb_fath`; for mother's age and own age it is 
`r cor_flb_moth`.

[^only-eldest]: Parental AFLB can only be calculated for this group.


```{r plot-age-flb-dv, fig.cap = "Effects of polygenic scores on age at first live birth.", fig.align = "center"}

drake::loadd(res_age_flb_dv)

n_regs <- as.double(nrow(res_age_flb_dv))

res_age_flb_dv %>% standard_ggplot(n_regs = n_regs)

```


```{r plot-age-birth-parents-dv, fig.cap = "Effects of polygenic scores on parents' age at respondent's birth, eldest siblings. Weights multiplied by 1/number of siblings.", fig.align = "center"}

drake::loadd(res_age_birth_parents_dv)

n_regs <- as.double(nrow(res_age_birth_parents_dv))

res_age_birth_parents_dv %>% 
      mutate(
        "Dependent variable" = dplyr::recode(dep.var, 
          "fath_age_birth" = "Father's AFLB", 
          "moth_age_birth" = "Mother's AFLB"
        )
      ) %>% 
      standard_ggplot(fill_col = `Dependent variable`, n_regs = n_regs, 
                        fill_direction = -1)

```


\FloatBarrier
\clearpage
## Mediation analysis

We run a standard mediation analysis in the framework of @baron1986moderator.
For each polygenic score where the bivariate correlation with RLRS is 
significant at $p$ < 0.05/33, we estimate 
\begin{align}
RLRS_i & = \alpha + \beta PGS_i + \gamma EA_i + X_i\mu + \varepsilon_i \\
EA_i & = \delta + \zeta PGS_i + X_i\mu + \eta_i
\end{align}
where $RLRS_i$ is relative lifetime reproductive success, $PGS_i$ is the
polygenic score, $EA_i$ is educational attainment (age of leaving fulltime
education), and $X_i$ is a vector of controls. The total effect of $PGS$ on
$RLRS$ is $\beta + \gamma \zeta$. The "indirect effect" mediated by
$EA$ is $\gamma \zeta$. The standard error of the indirect effect can be
calculated as

\[
\sqrt{\hat\gamma^2 \hat\sigma_\zeta^2 + \hat\zeta^2\hat\sigma_\gamma^2}
\]

where $\hat\sigma_\zeta$ is the standard error of $\hat\zeta$, etc. We include 
controls for age and sex in $X$. 

```{r tbl-res-mediation}

drake::loadd(res_mediation)

n_regs <- nrow(res_mediation)
res_mediation %<>% mutate(
                     # 0.025 not 0.05 because two-sided
                     sig  = abs(statistic_ind) > qnorm(1 - 0.025/n_regs),
                     star = ifelse(sig, "*", "")
                   )


n_pos_sig <- sum(res_mediation$sig & res_mediation$prop_ind > 0)

footnote <- sprintf("* p < 0.05/%s. Analysis run on %s PGS which correlated significantly with fertility.", n_regs, n_regs)

res_mediation %>%
      arrange(desc(estimate_total)) %>%
      transmute(
        PGS = pretty_names(term),
        `Total effect`    = estimate_total,
        `Indirect effect` = paste(estimate_ind, star),
        `Proportion (%)`  = prop_ind * 100,
        bounded_ci        = sign(estimate_total_conf_low) ==
                              sign(estimate_total_conf_high),
        `Proportion 95% c.i. (%)` = glue::glue(
                     "[{prop_ind_conf_low * 100}, {prop_ind_conf_high * 100}]"),
        `Proportion 95% c.i. (%)` = ifelse(bounded_ci, 
                                             `Proportion 95% c.i. (%)`,
                                             "Unbounded")
      ) %>% 
      select(-bounded_ci) %>%
      as_hux() %>%
      set_align(everywhere, -1, "right") %>% 
      set_number_format(-1, 2:3, 4) %>% 
      set_number_format(-1, 4:5, 1) %>% 
      set_caption("Mediation analysis") %>% 
      theme_compact() %>% 
      set_width(0.95) %>%
      set_col_width(c(.25, .15, .15, .15, .25)) %>% 
      set_bottom_border(final(), everywhere) %>%
      set_font_size(9) %>% 
      add_footnote(footnote, font_size = 8)
```

Table \@ref(tab:tbl-res-mediation) shows results. For `r n_pos_sig` out of 
`r n_regs` scores, the indirect effect on fertility via human capital is
significantly different from 0 at $p$ = 0.05/`r n_regs` and has the same sign
as the total effect. We also calculate the proportion of the total effect that
is mediated via the indirect effect, along with uncorrected 95% confidence
intervals (100 bootstraps). Note that if the confidence interval for the total
effect contains zero, the confidence interval for the proportion may be unbounded
[@franz2007ratios].

\FloatBarrier
\clearpage
## Within-siblings regressions


```{r calc-fe-fertility}
drake::loadd(res_fe_fertility)
drake::loadd(res_all)

educ_reduction <- with(res_fe_fertility, 
       estimate[Regression == "Controlled" & term == "EA3_excl_23andMe_UK"] / 
       estimate[Regression == "Raw" & term == "EA3_excl_23andMe_UK"]
     )
educ_reduction <- 1 - educ_reduction

res_all %<>% filter(dep.var == "RLRS", reg.type == "controlled")
res_fe_fertility %<>% filter(Regression == "Raw")
res_all %<>% left_join(res_fe_fertility, by = "term")

lm_pool_w <- tidy(lm(estimate.y ~ estimate.x, res_all))
coef_pool_w <- lm_pool_w[[2, "estimate"]]
se_pool_w <- lm_pool_w[[2, "std.error"]]

```



```{r plot-res-fe-fertility, fig.cap = "Selection effects controlling for sibling-group fixed effects, with and without a control for education (left education before 16, 16-18, or after 18). Each set of 29 results is from a single regression of RLRS on 29 polygenic scores. Standard errors clustered by sibling group."}

drake::loadd(res_fe_fertility)
n_regs <- nrow(res_fe_fertility)/2 # Raw and Controlled treated separately

res_fe_fertility %>% 
  mutate(
    Controls = ifelse(Regression == "Controlled", "FE + education", "FE only")
  ) %>% 
  standard_ggplot(fill_col = Controls, n_regs = n_regs, fill_direction = -1) +
    guides(alpha = guide_legend(
                   override.aes = list(
                             shape = 1,
                             color = "black",
                             alpha = 1,
                             size = 1.8
                   )
                 )
    )

n_groups <- attr(res_fe_fertility, "n_groups")
n <- attr(res_fe_fertility, "n")

```


Results in the main text support our theory that natural selection on polygenic scores is
driven by their *correlation* with human capital. Here, we test whether polygenic
scores *cause* fertility by running within-siblings regressions. We run a single regression on 29 polygenic scores within `r n_groups` sibling
groups (N = `r n`). Thus, we control both for environmental confounds (since
scores are randomly allocated within sib-groups by meiosis), and for genetic
confounds captured by our polygenic scores. We remove four scores which
correlate highly with other scores (educational attainment 2, hip circumference,
waist circumference and waist-hip ratio). Figure \@ref(fig:plot-res-fe-fertility)
shows the results.

With a reduced sample size, all within-sibling effects are insignificant after
Bonferroni correction. However, effect sizes are positively correlated with
effect sizes from the pooled model, and about 70% smaller (regressing
within-sibling on pooled effect sizes, $b$ = `r format(coef_pool_w, digits = 3)`).
This attenuation is broadly consistent with the decrease in heritability
in within-sibling GWASs on age at first birth and educational attainment
[@Howe2021.03.05.433935]. We see these results as providing tentative evidence
that polygenic scores cause fertility, with effects being partly driven by
correlations with environmental variation in human capital. We also reran
within-siblings regressions adding a control for education. Most effect sizes
barely change, suggesting that our measure of education does not in general
mediate differences in fertility among siblings.


\clearpage
\FloatBarrier
## Effects on inequality

Table \@ref(tab:tbl-inequality-effects) shows correlations between children's
polygenic scores and household income (UKB data field 738). Column "With
selection" uses respondents' scores, multiplying weights by number
of children. Column "Without selection" uses our standard weights, i.e. 
it estimates the counterfactual correlation if all respondents
had the same number of children.


```{r tbl-inequality-effects}


res_ineq %>% 
      select(
        PGS                      = score, 
        `Cor. with selection`    = actual, 
        `Cor. without selection` = cf,
        `Ratio`                  = ratio
      ) %>% 
      mutate(
        PGS = pretty_names(PGS)
      ) %>% 
      arrange(desc(`Cor. with selection`)) %>% 
      as_hux() %>% 
      set_align(everywhere, -1, "right") %>% 
      set_number_format(everywhere, 2:3, 3) %>% 
      set_number_format(everywhere, 4, 2) %>% 
      set_caption("Correlations of polygenic scores with income group.") %>% 
      theme_compact() %>% 
      set_width(0.75) %>% 
      set_col_width(c(.4, .2, .2, .2)) %>% 
      set_bottom_border(final(1), everywhere) %>% 
      set_position("center")
```


\clearpage
\FloatBarrier
## Further results
### Selection effects on raw polygenic scores


```{r calc-regs-controlled}

drake::loadd(res_all)

res_all_controlled <- res_all %>% 
  select(reg.type, dep.var, term, estimate) %>% 
  pivot_wider(names_from = c(reg.type, dep.var), values_from = estimate) %>% 
  mutate(
    consistent = sign(`raw_RLRS_parents`) == sign(`raw_RLRS`) &
                 sign(`raw_RLRS_parents`) == sign(`controlled_RLRS_parents`) &
                 sign(`raw_RLRS_parents`) == sign(`controlled_RLRS`)
  )

n_pgs  <- as.double(nrow(res_all_controlled))
prop_consistent_controlled <- mean(res_all_controlled$consistent)
controlled_smaller_sibs <- abs(res_all_controlled[["controlled_RLRS_parents"]]) < 
  abs(res_all_controlled[["raw_RLRS"]])
controlled_smaller_children <- abs(res_all_controlled[["controlled_RLRS"]]) < 
  abs(res_all_controlled[["raw_RLRS"]])
median_prop_sibs <- median(abs(res_all_controlled[["controlled_RLRS_parents"]]) /
    abs(res_all_controlled[["raw_RLRS_parents"]]))
median_prop_children <- median(abs(res_all_controlled[["controlled_RLRS"]]) /
    abs(res_all_controlled[["raw_RLRS"]]))

```

Figure \@ref(fig:plot-regs-controlled-pcs) compares selection effects on
polygenic scores residualized for the top 100 principal components of the
genetic data, to selection effects on raw, unresidualized polygenic scores. In
siblings regressions, effect sizes are larger for raw scores  -- sometimes much
larger, as in the case of height. `r sum(controlled_smaller_sibs)` out of 
`r n_pgs` "raw" effect sizes have a larger absolute value than the corresponding
"residualized" effect size. The median proportion between raw and controlled
effect sizes is `r median_prop_sibs`. Among the children regressions, this no
longer holds. Effect sizes are barely affected by controlling for principal
components.

Overall, `r prop_consistent_controlled * 100` per cent of effect sizes
are consistently signed across all four regressions (on children and
siblings, and with and without residualization).

```{r plot-regs-controlled-pcs, fig.cap = "Selection effects using unresidualized polygenic scores. Parental generation weights multiplied by 1/number of siblings.", fig.subcap = c("Respondents", "Parents"), fig.ncol = 1}

drake::loadd(res_all)

n_regs <- as.double(nrow(res_all))
res_all %<>% mutate(
               `PGS` = ifelse(reg.type == "controlled", 
                                "Residualized for 100 PCs", "Unresidualized")
             )

res_all %>% 
      filter(dep.var == "RLRS") %>% 
      standard_ggplot(fill_col = PGS, n_regs = n_regs, fill_direction = -1)

res_all %>% 
      filter(dep.var == "RLRS_parents") %>% 
      standard_ggplot(fill_col = PGS, n_regs = n_regs, fill_direction = -1)

```


To get a further insight into this we regress respondents' and parents' RLRS
on individual principal components. Figure \@ref(fig:plot-pc-regs)
shows the results. Labels show the top principal components. These
have larger effect sizes in siblings regressions. One possibility is that 
the parents' generation was less geographically mobile, and so geographic
patterns of childrearing were more correlated with principal components, which
partly capture the location of people's ancestors.


```{r plot-pc-regs, fig.cap = "Selection effects of 100 principal components of genetic data. Each dot represents one bivariate regression. Parental generation weights multiplied by 1/number of siblings. Absolute effect sizes are plotted. Points are jittered on the Y axis. Top principal components are labelled."}

drake::loadd(res_pcs)
n_regs <- as.double(nrow(res_pcs))

pj <- position_jitter(height = 0.15, width = 0, seed = 12345)

res_pcs %>% 
      mutate(
        estimate = abs(estimate),
        dep.var  = fct_relevel(dep.var, "RLRS"),
        dep.var  = fct_recode(dep.var,
          "Parents"     = "RLRS_parents",
          "Respondents" = "RLRS"
        ),
        top_pcs  = ifelse(as.numeric(pc_name) <= 5, term, "")
      ) %>% 
      rename(p = p.value) %>% 
      ggplot(aes(estimate, dep.var, shape = p < 0.05/{{n_regs}})) + 
        geom_point(alpha = 0.8, color = "steelblue4", position = pj) +
        geom_text(aes(estimate, dep.var, label = top_pcs), position = pj, 
                    hjust = "right", color = "black", size = 3,
                    family = "Abadi MT Condensed Light") +
        scale_x_log10(labels = scales::label_number(accuracy = 0.0001)) +
        labs(x = "Effect size (log scale)", y = "") + 
        scale_shape_manual(values = c("FALSE" = "circle open", "TRUE" = "circle"))
        

```


\FloatBarrier
\clearpage

### Genetic correlations with EA3


```{r plot-rgs-by-effect-size, fig.cap = "Selection effects plotted against genetic correlation with EA3."}

# TODO: ref source of data.

drake::loadd(rgs)
drake::loadd(res_all)

rgs <- left_join(rgs, 
        res_all %>% filter(reg.type == "controlled", dep.var == "RLRS"), 
        by = c("p2" = "term"))

rgs_cor <- cor.test(~ rg + estimate, 
                      data = rgs %>%
                               filter(
                                 p2 != "EA3_excl_23andMe_UK", 
                                 p2 != "EA2_noUKB"
                               )
                    )
rgs %>% 
      filter(p2 != "EA3_excl_23andMe_UK", p2 != "EA2_noUKB") %>% 
      ggplot(aes(rg, estimate)) + 
        geom_point() +
        my_vline +
        my_hline +
        labs(x = "Genetic correlation with EA3", y = "Effect on RLRS")
          
```

Another way to examine the "earnings" theory of natural selection is to compare
selection effects of polygenic scores with their genetic correlation with educational
attainment (EA3). Since EA3 strongly predicts earnings, if earnings drives
differences in fertility, we'd expect a correlation between the two sets of
results. Figure \@ref(fig:plot-rgs-by-effect-size) shows this is so: the 
correlation, after excluding EA2, is `r rgs_cor$estimate`. Genetic
correlations were calculated using LD score regression from GWAS summary
statistics.


\clearpage
\FloatBarrier

## Model proofs

```{r include-model-appendix, child="model-appendix.Rmd"}

```

\clearpage
\FloatBarrier

# References