Continuous Effects of Early Adversities

Library and Data

library(survival)
library(kableExtra)
library(tidyverse)
library(GGally)
library(penalized)

raw_data <- read_csv("data/Aim3-early-adversity-subjects-F-survival-2021-02-17.csv")
raw_data_sub <- raw_data %>%
  dplyr::select(dead, age,
                rain, density, mom_rank, mom_sci,
                adv_mom, adv_sib)

covariate_names <- c("adv_mom", "adv_sib",
                     "rain", "density", "mom_rank", "mom_sci")
covariate_types <- c(rep("Logical", 2), rep("Numeric", 4))
covariate_meanings <- c("Mom dies before age 4",
                       "A close sibling is born",
                       "Group size at birth",
                       "Rainfall in the first year of life",
                       "Mom rank",
                       "Mom Social connection")
covariate_summary <- c("139/562",
                       "149/562",
                       "340.5(134.7)",
                       "26.6(8.6)",
                       "0.5(0.3)",
                       "0.1(0.75)")

data.frame(
  list(types = covariate_types,
       meanings = covariate_meanings,
       summary = covariate_summary)
) %>%
  `row.names<-`(covariate_names) %>%
  kable(col.names = c("Type", "Description", "Summary")) %>%
  kable_styling(full_width = FALSE, position = "center", font_size = 18)

	Type	Description	Summary
adv_mom	Logical	Mom dies before age 4	139/562
adv_sib	Logical	A close sibling is born	149/562
rain	Numeric	Group size at birth	340.5(134.7)
density	Numeric	Rainfall in the first year of life	26.6(8.6)
mom_rank	Numeric	Mom rank	0.5(0.3)
mom_sci	Numeric	Mom Social connection	0.1(0.75)

ggpairs(data = dplyr::select(raw_data_sub, - c(dead, age)),
        lower = list(continuous = wrap("points", size = 0.5, alpha = 0.5),
                     combo = wrap("facethist", bins = 20))) +
  theme_classic()

Process

age_below_1 <- raw_data_sub %>%
  mutate(dead = if_else(age > 1, FALSE, dead)) %>%
  mutate(age = if_else(age > 1, 1, age))

age_above_1_below_5 <- raw_data %>%
  filter(age >= 1) %>%
  mutate(dead = if_else(age > 5, FALSE, dead)) %>%
  mutate(age = if_else(age > 5, 5, age))

age_above_5 <- raw_data %>%
  filter(age >= 5)

As before, the dataset is divided into infant, juvenile, and adults phase based on the cut of age at 1 and 5.
For each subgroup
- First run a simple regression with one variable at a time. (e.g. survival ~ rain + (maybe rain^2)).
- Then put all variables in one model.

Age below 1

rain

g1_rain_coxph <- coxph(Surv(age, dead) ~ pspline(rain, df = 4), data = age_below_1)
g1_rain_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(rain, df = 4), data = age_below_1)
## 
##                               coef se(coef)      se2    Chisq   DF     p
## pspline(rain, df = 4), li 0.001393 0.000655 0.000655 4.516923 1.00 0.034
## pspline(rain, df = 4), no                            1.299434 3.05 0.737
## 
## Iterations: 4 outer, 12 Newton-Raphson
##      Theta= 0.751 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=6.12  on 4.05 df, p=0.2
## n= 562, number of events= 120

density

g1_density_coxph <- coxph(Surv(age, dead) ~ pspline(density, df = 4), data = age_below_1)
g1_density_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(density, df = 4), data = age_below_1)
## 
##                              coef se(coef)     se2   Chisq   DF    p
## pspline(density, df = 4), -0.0109   0.0108  0.0108  1.0054 1.00 0.32
## pspline(density, df = 4),                           3.5655 3.07 0.32
## 
## Iterations: 4 outer, 11 Newton-Raphson
##      Theta= 0.777 
## Degrees of freedom for terms= 4.1 
## Likelihood ratio test=6.67  on 4.07 df, p=0.2
## n= 562, number of events= 120

mom_rank

g1_momrank_coxph <- coxph(Surv(age, dead) ~ pspline(mom_rank, df = 4), data = age_below_1)
g1_momrank_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(mom_rank, df = 4), 
##     data = age_below_1)
## 
##                             coef se(coef)    se2  Chisq   DF    p
## pspline(mom_rank, df = 4) -0.206    0.284  0.283  0.523 1.00 0.47
## pspline(mom_rank, df = 4)                         4.231 3.06 0.25
## 
## Iterations: 5 outer, 11 Newton-Raphson
##      Theta= 0.857 
## Degrees of freedom for terms= 4.1 
## Likelihood ratio test=5.56  on 4.06 df, p=0.2
## n= 562, number of events= 120

mom_sci

g1_momsci_coxph <- coxph(Surv(age, dead) ~ pspline(mom_sci, df = 4), data = age_below_1)
g1_momsci_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(mom_sci, df = 4), data = age_below_1)
## 
##                            coef se(coef)   se2 Chisq DF     p
## pspline(mom_sci, df = 4), 0.262    0.141 0.141 3.443  1 0.064
## pspline(mom_sci, df = 4),                      5.254  3 0.154
## 
## Iterations: 5 outer, 17 Newton-Raphson
##      Theta= 0.52 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=12.4  on 4 df, p=0.01
## n= 562, number of events= 120

All in one model

g1_all <- coxph(
  Surv(age, dead) ~ adv_mom + rain + density + mom_rank + mom_sci,
  data = age_below_1
)
print(g1_all)

## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom + rain + density + 
##     mom_rank + mom_sci, data = age_below_1)
## 
##                   coef  exp(coef)   se(coef)      z      p
## adv_momTRUE  0.4688391  1.5981378  0.1949505  2.405 0.0162
## rain         0.0014374  1.0014384  0.0006481  2.218 0.0266
## density     -0.0068457  0.9931777  0.0108498 -0.631 0.5281
## mom_rank    -0.2823644  0.7539988  0.3030402 -0.932 0.3515
## mom_sci      0.2791812  1.3220469  0.1253020  2.228 0.0259
## 
## Likelihood ratio test=15.83  on 5 df, p=0.007349
## n= 562, number of events= 120

Age above 1 and below 5

rain

g2_rain_coxph <- coxph(Surv(age, dead) ~ pspline(rain, df = 4), data = age_above_1_below_5)
g2_rain_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(rain, df = 4), data = age_above_1_below_5)
## 
##                                coef  se(coef)       se2     Chisq   DF    p
## pspline(rain, df = 4), li -0.001057  0.000927  0.000927  1.301235 1.00 0.25
## pspline(rain, df = 4), no                                3.123290 3.04 0.38
## 
## Iterations: 2 outer, 6 Newton-Raphson
##      Theta= 0.616 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=6.07  on 4.04 df, p=0.2
## n= 440, number of events= 92

density

g2_density_coxph <- coxph(Surv(age, dead) ~ pspline(density, df = 4), data = age_above_1_below_5)
g2_density_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(density, df = 4), data = age_above_1_below_5)
## 
##                               coef se(coef)      se2    Chisq   DF     p
## pspline(density, df = 4), -0.00255  0.01226  0.01225  0.04322 1.00 0.835
## pspline(density, df = 4),                             9.36752 3.08 0.026
## 
## Iterations: 3 outer, 11 Newton-Raphson
##      Theta= 0.694 
## Degrees of freedom for terms= 4.1 
## Likelihood ratio test=12.8  on 4.08 df, p=0.01
## n= 440, number of events= 92

mom_rank

g2_momrank_coxph <- coxph(Surv(age, dead) ~ pspline(mom_rank, df = 4), data = age_above_1_below_5)
g2_momrank_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(mom_rank, df = 4), 
##     data = age_above_1_below_5)
## 
##                             coef se(coef)    se2  Chisq   DF    p
## pspline(mom_rank, df = 4) -0.285    0.342  0.341  0.692 1.00 0.41
## pspline(mom_rank, df = 4)                         1.668 3.04 0.65
## 
## Iterations: 5 outer, 12 Newton-Raphson
##      Theta= 0.823 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=3.33  on 4.04 df, p=0.5
## n= 440, number of events= 92

mom_sci

g2_momsci_coxph <- coxph(Surv(age, dead) ~ pspline(mom_sci, df = 4), data = age_above_1_below_5)
g2_momsci_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(mom_sci, df = 4), data = age_above_1_below_5)
## 
##                             coef se(coef)    se2  Chisq   DF    p
## pspline(mom_sci, df = 4), 0.0852   0.1348 0.1348 0.3996 1.00 0.53
## pspline(mom_sci, df = 4),                        1.7927 3.05 0.62
## 
## Iterations: 3 outer, 7 Newton-Raphson
##      Theta= 0.67 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=3.02  on 4.05 df, p=0.6
## n= 440, number of events= 92

All in one model

g2_all <- coxph(
  Surv(age, dead) ~ adv_mom + adv_sib + rain + density + mom_rank + mom_sci,
  data = age_above_1_below_5
)
print(g2_all)

## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom + adv_sib + rain + 
##     density + mom_rank + mom_sci, data = age_above_1_below_5)
## 
##                   coef  exp(coef)   se(coef)      z        p
## adv_momTRUE  0.7685575  2.1566531  0.2175745  3.532 0.000412
## adv_sibTRUE -0.0428078  0.9580956  0.2873361 -0.149 0.881568
## rain        -0.0008726  0.9991278  0.0008628 -1.011 0.311831
## density     -0.0054389  0.9945758  0.0122728 -0.443 0.657645
## mom_rank    -0.3355305  0.7149587  0.3470404 -0.967 0.333627
## mom_sci      0.1455543  1.1566805  0.1389974  1.047 0.295020
## 
## Likelihood ratio test=14.52  on 6 df, p=0.02432
## n= 440, number of events= 92

Age above 5

rain

g3_rain_coxph <- coxph(Surv(age, dead) ~ pspline(rain, df = 4), data = age_above_5)
g3_rain_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(rain, df = 4), data = age_above_5)
## 
##                                coef  se(coef)       se2     Chisq   DF    p
## pspline(rain, df = 4), li -0.000418  0.000728  0.000728  0.330270 1.00 0.57
## pspline(rain, df = 4), no                                4.648255 3.03 0.20
## 
## Iterations: 4 outer, 11 Newton-Raphson
##      Theta= 0.726 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=6.71  on 4.03 df, p=0.2
## n= 277, number of events= 110

density

g3_density_coxph <- coxph(Surv(age, dead) ~ pspline(density, df = 4), data = age_above_5)
g3_density_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(density, df = 4), data = age_above_5)
## 
##                              coef se(coef)     se2   Chisq   DF     p
## pspline(density, df = 4), -0.0022   0.0130  0.0130  0.0289 1.00 0.865
## pspline(density, df = 4),                           8.9465 3.09 0.032
## 
## Iterations: 2 outer, 7 Newton-Raphson
##      Theta= 0.637 
## Degrees of freedom for terms= 4.1 
## Likelihood ratio test=12.7  on 4.09 df, p=0.01
## n= 277, number of events= 110

mom_rank

g3_momrank_coxph <- coxph(Surv(age, dead) ~ pspline(mom_rank, df = 4), data = age_above_5)
g3_momrank_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(mom_rank, df = 4), 
##     data = age_above_5)
## 
##                              coef se(coef)     se2   Chisq   DF    p
## pspline(mom_rank, df = 4) -0.0919   0.3241  0.3230  0.0804 1.00 0.78
## pspline(mom_rank, df = 4)                           2.8263 3.05 0.43
## 
## Iterations: 5 outer, 11 Newton-Raphson
##      Theta= 0.838 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=3.5  on 4.05 df, p=0.5
## n= 277, number of events= 110

mom_sci

g3_momsci_coxph <- coxph(Surv(age, dead) ~ pspline(mom_sci, df = 4), data = age_above_5)
g3_momsci_coxph

## Call:
## coxph(formula = Surv(age, dead) ~ pspline(mom_sci, df = 4), data = age_above_5)
## 
##                             coef se(coef)    se2  Chisq   DF    p
## pspline(mom_sci, df = 4), -0.154    0.130  0.130  1.395 1.00 0.24
## pspline(mom_sci, df = 4),                         2.306 3.04 0.52
## 
## Iterations: 2 outer, 6 Newton-Raphson
##      Theta= 0.617 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=4.86  on 4.04 df, p=0.3
## n= 277, number of events= 110

All in one model

g3_all <- coxph(
  Surv(age, dead) ~ adv_mom + adv_sib + rain + density + mom_rank + mom_sci,
  data = age_above_5
)
print(g3_all)

## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom + adv_sib + rain + 
##     density + mom_rank + mom_sci, data = age_above_5)
## 
##                   coef  exp(coef)   se(coef)      z       p
## adv_momTRUE  0.6955337  2.0047787  0.2222889  3.129 0.00175
## adv_sibTRUE  0.4574650  1.5800635  0.2548730  1.795 0.07267
## rain        -0.0006612  0.9993390  0.0007203 -0.918 0.35862
## density     -0.0014830  0.9985181  0.0131040 -0.113 0.90989
## mom_rank    -0.1304891  0.8776660  0.3307282 -0.395 0.69317
## mom_sci     -0.2204143  0.8021864  0.1228480 -1.794 0.07278
## 
## Likelihood ratio test=14.36  on 6 df, p=0.02582
## n= 277, number of events= 110

Discussion

Few of rain, density, mom_rank, and mom_sci have shown significance in the cox model on the three datasets, either as a single variable or put together.
For the infant group(age < 1), rain is significant with p-value 0.03. However, no significance is shown for the nonlinear part.
Nonlinear part of density is significant for the juvenile(1 < age < 5) and adult(age > 5) group, with p-value 0.023 and 0.03 respectively.