Library and Data

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

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) %>%
  mutate(adv_mom = as.numeric(adv_mom),
         adv_sib = as.numeric(adv_sib)) %>%
  filter(age > 0)
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 = raw_data_sub %>%
          select(- c(dead, age)) %>%
          mutate(adv_mom = as.factor(adv_mom),
                 adv_sib = as.factor(adv_sib)),
        lower = list(continuous = wrap("points", size = 0.5, alpha = 0.5),
                     combo = wrap("facethist", bins = 20))) +
  theme_classic()


Process

  • 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.
    • Lastly, include the interactions which is significant when modeled individually, apply variable selection using lasso, and run the model on picked variables.


Age below 1

summary

covariate_names <- c("adv_mom", "rain", "density", "mom_rank", "mom_sci")

covariate_lin_pval <- c(0.0198, 0.041, 0.36, 0.51, 0.048)
covariate_non_pval <- c(NA, 0.736, 0.29, 0.25, 0.157)

options(knitr.kable.NA = '')
data.frame(
  list(linear_pval = covariate_lin_pval,
       nonlinear_pval = covariate_non_pval)
) %>%
   `row.names<-`(covariate_names) %>%
  kable(col.names = c("linear p-val", "nonlinear p-val")) %>%
  kable_styling(full_width = FALSE, position = "center", font_size = 18)
linear p-val nonlinear p-val
adv_mom 0.0198
rain 0.0410 0.736
density 0.3600 0.290
mom_rank 0.5100 0.250
mom_sci 0.0480 0.157

adv_mom

g1_adv_mom_coxph <- coxph(Surv(age, dead) ~ adv_mom, data = age_below_1)
g1_adv_mom_coxph
## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom, data = age_below_1)
## 
##           coef exp(coef) se(coef)     z      p
## adv_mom 0.4524    1.5721   0.1941 2.331 0.0198
## 
## Likelihood ratio test=5.14  on 1 df, p=0.0234
## n= 561, number of events= 119

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.001345 0.000657 0.000657 4.190073 1.00 0.041
## pspline(rain, df = 4), no                            1.303057 3.05 0.736
## 
## Iterations: 4 outer, 12 Newton-Raphson
##      Theta= 0.75 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=5.85  on 4.05 df, p=0.2
## n= 561, number of events= 119

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.00988  0.01088  0.01087  0.82551 1.00 0.36
## pspline(density, df = 4),                             3.81747 3.06 0.29
## 
## Iterations: 4 outer, 11 Newton-Raphson
##      Theta= 0.776 
## Degrees of freedom for terms= 4.1 
## Likelihood ratio test=6.79  on 4.06 df, p=0.2
## n= 561, number of events= 119

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.188    0.285  0.284  0.434 1.00 0.51
## pspline(mom_rank, df = 4)                         4.227 3.06 0.25
## 
## Iterations: 5 outer, 12 Newton-Raphson
##      Theta= 0.856 
## Degrees of freedom for terms= 4.1 
## Likelihood ratio test=5.46  on 4.06 df, p=0.3
## n= 561, number of events= 119

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.285    0.144 0.144 3.923  1 0.048
## pspline(mom_sci, df = 4),                      5.211  3 0.157
## 
## Iterations: 5 outer, 17 Newton-Raphson
##      Theta= 0.51 
## Degrees of freedom for terms= 4 
## Likelihood ratio test=12.8  on 4 df, p=0.01
## n= 561, number of events= 119

Age below 1: All in one model

Main effects

g1_all <- coxph(
  Surv(age, dead) ~ adv_mom + rain + mom_sci,
  data = age_below_1
)
print(g1_all)
## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom + rain + mom_sci, data = age_below_1)
## 
##              coef exp(coef)  se(coef)     z     p
## adv_mom 0.4948759 1.6402947 0.1947243 2.541 0.011
## rain    0.0014728 1.0014739 0.0006479 2.273 0.023
## mom_sci 0.2928659 1.3402630 0.1248956 2.345 0.019
## 
## Likelihood ratio test=15.21  on 3 df, p=0.001648
## n= 561, number of events= 119

Lasso path

penalized_covariates <- model.matrix(
  ~ adv_mom + rain + mom_sci +
    rain:mom_sci + adv_mom:rain + adv_mom:mom_sci,
  data = age_below_1
) %>% .[, colnames(.) != "(Intercept)"]

opt <- optL1(
  Surv(age_below_1$age, age_below_1$dead),
  penalized = penalized_covariates,
  standardize = T,
  fold = 10,
  trace = F
)

stepwise <- penalized(
  Surv(age_below_1$age, age_below_1$dead),
  penalized = penalized_covariates,
  standardize = T,
  steps = 20,
  trace = F
)

plotpath(
  stepwise,
  labelsize = 0.8,
  standardize = T
)

abline(
  v = opt$lambda,
  col = "gray",
  lwd = 2
)

Results with picked predictors

g1_opt_lasso <- penalized(
  Surv(age_below_1$age, age_below_1$dead),
  penalized = penalized_covariates,
  standardize = T,
  lambda1 = opt$lambda,
  trace = F
)

nonzero_varnames <- names(g1_opt_lasso@penalized)[g1_opt_lasso@penalized != 0]
formula <- paste0(" Surv(age, dead) ~ ", paste(nonzero_varnames, collapse = " + "),
                  " + mom_sci + ", "adv_mom")
g1_final_coxph <- coxph(
  as.formula(formula),
  data = age_below_1
)

g1_final_coxph
## Call:
## coxph(formula = as.formula(formula), data = age_below_1)
## 
##                       coef  exp(coef)   se(coef)      z     p
## rain             0.0008125  1.0008128  0.0008521  0.953 0.340
## mom_sci         -0.0095839  0.9904619  0.3438678 -0.028 0.978
## adv_mom         -0.0685234  0.9337716  0.5356914 -0.128 0.898
## rain:mom_sci     0.0007175  1.0007177  0.0008498  0.844 0.398
## rain:adv_mom     0.0014948  1.0014959  0.0013616  1.098 0.272
## mom_sci:adv_mom  0.1602615  1.1738177  0.2661212  0.602 0.547
## 
## Likelihood ratio test=17.09  on 6 df, p=0.008959
## n= 561, number of events= 119

Interaction plot



Age above 1 and below 5

summary

covariate_names <- c("adv_mom", "adv_sib", "rain", "density", "mom_rank", "mom_sci")

covariate_lin_pval <- c(0.0003, 0.998, 0.25, 0.835, 0.41, 0.53)
covariate_non_pval <- c(NA, NA, 0.38, 0.026, 0.65, 0.62)

options(knitr.kable.NA = '')
data.frame(
  list(linear_pval = covariate_lin_pval,
       nonlinear_pval = covariate_non_pval)
) %>%
   `row.names<-`(covariate_names) %>%
  kable(col.names = c("linear p-val", "nonlinear p-val")) %>%
  kable_styling(full_width = FALSE, position = "center", font_size = 18)
linear p-val nonlinear p-val
adv_mom 0.0003
adv_sib 0.9980
rain 0.2500 0.380
density 0.8350 0.026
mom_rank 0.4100 0.650
mom_sci 0.5300 0.620

adv_mom

g2_adv_mom_coxph <- coxph(Surv(age, dead) ~ adv_mom, data = age_above_1_below_5)
g2_adv_mom_coxph
## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom, data = age_above_1_below_5)
## 
##           coef exp(coef) se(coef)     z        p
## adv_mom 0.7654    2.1498   0.2160 3.543 0.000395
## 
## Likelihood ratio test=11.52  on 1 df, p=0.0006876
## n= 440, number of events= 92

adv_sib

g2_adv_sib_coxph <- coxph(Surv(age, dead) ~ adv_sib, data = age_above_1_below_5)
g2_adv_sib_coxph
## Call:
## coxph(formula = Surv(age, dead) ~ adv_sib, data = age_above_1_below_5)
## 
##              coef exp(coef)  se(coef)     z     p
## adv_sib 0.0008429 1.0008433 0.2751020 0.003 0.998
## 
## Likelihood ratio test=0  on 1 df, p=0.9976
## n= 440, number of events= 92

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

Age above 1 below 5: All in one model

Main effects

g2_all <- coxph(
  Surv(age, dead) ~ adv_mom + pspline(density, df = 4),
  data = age_above_1_below_5
)
print(g2_all)
## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom + pspline(density, 
##     df = 4), data = age_above_1_below_5)
## 
##                               coef se(coef)      se2    Chisq   DF       p
## adv_mom                    0.75826  0.21654  0.21625 12.26152 1.00 0.00046
## pspline(density, df = 4), -0.00263  0.01227  0.01226  0.04590 1.00 0.83035
## pspline(density, df = 4),                             9.17130 3.07 0.02875
## 
## Iterations: 3 outer, 11 Newton-Raphson
##      Theta= 0.691 
## Degrees of freedom for terms= 1.0 4.1 
## Likelihood ratio test=24.1  on 5.07 df, p=2e-04
## n= 440, number of events= 92
par(mfrow = c(1, 2))
termplot(g2_all, terms = c(1, 2), se=T, ylabs = "Log-hazard")

Lasso path

penalized_covariates <- model.matrix(
  ~adv_mom + adv_sib + rain + density + mom_rank + mom_sci +
    adv_mom:density,
  data = age_above_1_below_5
) %>% .[, colnames(.) != "(Intercept)"]

opt <- optL1(
  Surv(age_above_1_below_5$age, age_above_1_below_5$dead),
  penalized = penalized_covariates,
  standardize = T,
  fold = 10,
  trace = F
)

stepwise <- penalized(
  Surv(age_above_1_below_5$age, age_above_1_below_5$dead),
  penalized = penalized_covariates,
  standardize = T,
  steps = 20,
  trace = F
)

plotpath(
  stepwise,
  labelsize = 0.8,
  standardize = T
)

abline(
  v = opt$lambda,
  col = "gray",
  lwd = 2
)

Results with picked predictors

g2_opt_lasso <- penalized(
  Surv(age_above_1_below_5$age, age_above_1_below_5$dead),
  penalized = penalized_covariates,
  standardize = T,
  lambda1 = opt$lambda,
  trace = F
)

nonzero_varnames <- names(g2_opt_lasso@penalized)[g2_opt_lasso@penalized != 0]
formula <- paste0(" Surv(age, dead) ~ ", paste(nonzero_varnames, collapse = " + "))
g2_final_coxph <- coxph(
  as.formula(formula),
  data = age_above_1_below_5
)

g2_final_coxph
## Call:
## coxph(formula = as.formula(formula), data = age_above_1_below_5)
## 
##           coef exp(coef) se(coef)     z        p
## adv_mom 0.7654    2.1498   0.2160 3.543 0.000395
## 
## Likelihood ratio test=11.52  on 1 df, p=0.0006876
## n= 440, number of events= 92



Age above 5

summary

covariate_names <- c("adv_mom", "adv_sib", "rain", "density", "mom_rank", "mom_sci")

covariate_lin_pval <- c(0.00173, 0.203, 0.57, 0.865, 0.78, 0.24)
covariate_non_pval <- c(NA, NA, 0.20, 0.032, 0.43, 0.52)

options(knitr.kable.NA = '')
data.frame(
  list(linear_pval = covariate_lin_pval,
       nonlinear_pval = covariate_non_pval)
) %>%
   `row.names<-`(covariate_names) %>%
  kable(col.names = c("linear p-val", "nonlinear p-val")) %>%
  kable_styling(full_width = FALSE, position = "center", font_size = 18)
linear p-val nonlinear p-val
adv_mom 0.00173
adv_sib 0.20300
rain 0.57000 0.200
density 0.86500 0.032
mom_rank 0.78000 0.430
mom_sci 0.24000 0.520

adv_mom

g3_adv_mom_coxph <- coxph(Surv(age, dead) ~ adv_mom, data = age_above_5)
g3_adv_mom_coxph
## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom, data = age_above_5)
## 
##           coef exp(coef) se(coef)     z       p
## adv_mom 0.6881    1.9899   0.2197 3.133 0.00173
## 
## Likelihood ratio test=8.78  on 1 df, p=0.003042
## n= 277, number of events= 110

adv_sib

g3_adv_sib_coxph <- coxph(Surv(age, dead) ~ adv_sib, data = age_above_5)
g3_adv_sib_coxph
## Call:
## coxph(formula = Surv(age, dead) ~ adv_sib, data = age_above_5)
## 
##           coef exp(coef) se(coef)     z     p
## adv_sib 0.3019    1.3525   0.2371 1.273 0.203
## 
## Likelihood ratio test=1.53  on 1 df, p=0.2158
## n= 277, number of events= 110

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

Age above 5: All in one model

Main effects

g3_all <- coxph(
  Surv(age, dead) ~ adv_mom + pspline(density, df = 4),
  data = age_above_5
)
print(g3_all)
## Call:
## coxph(formula = Surv(age, dead) ~ adv_mom + pspline(density, 
##     df = 4), data = age_above_5)
## 
##                               coef se(coef)      se2    Chisq   DF       p
## adv_mom                    0.77232  0.22245  0.22160 12.05363 1.00 0.00052
## pspline(density, df = 4), -0.00205  0.01299  0.01297  0.02502 1.00 0.87431
## pspline(density, df = 4),                            10.88041 3.08 0.01334
## 
## Iterations: 2 outer, 7 Newton-Raphson
##      Theta= 0.634 
## Degrees of freedom for terms= 1.0 4.1 
## Likelihood ratio test=23.9  on 5.07 df, p=2e-04
## n= 277, number of events= 110
par(mfrow = c(1, 2))
termplot(g3_all, terms = c(1, 2), se=T, ylabs = "Log-hazard")

Lasso path

penalized_covariates <- model.matrix(
  ~adv_mom + adv_sib + rain + density + mom_rank + mom_sci +
    adv_mom:density,
  data = age_above_5
) %>% .[, colnames(.) != "(Intercept)"]

opt <- optL1(
  Surv(age_above_5$age, age_above_5$dead),
  penalized = penalized_covariates,
  standardize = T,
  fold = 10,
  trace = F
)

stepwise <- penalized(
  Surv(age_above_5$age, age_above_5$dead),
  penalized = penalized_covariates,
  standardize = T,
  steps = 20,
  trace = F
)

plotpath(
  stepwise,
  labelsize = 0.8,
  standardize = T
)

abline(
  v = opt$lambda,
  col = "gray",
  lwd = 2
)

Results with picked predictors

g3_opt_lasso <- penalized(
  Surv(age_above_5$age, age_above_5$dead),
  penalized = penalized_covariates,
  standardize = T,
  lambda1 = opt$lambda,
  trace = F
)

nonzero_varnames <- names(g3_opt_lasso@penalized)[g3_opt_lasso@penalized != 0]
formula <- paste0(" Surv(age, dead) ~ ", paste(nonzero_varnames, collapse = " + "))
g3_final_coxph <- coxph(
  as.formula(formula),
  data = age_above_5
)

g3_final_coxph
## Call:
## coxph(formula = as.formula(formula), data = age_above_5)
## 
##               coef  exp(coef)   se(coef)      z       p
## adv_mom  0.6906337  1.9949793  0.2213202  3.121 0.00181
## adv_sib  0.4484606  1.5658998  0.2481102  1.808 0.07068
## rain    -0.0006586  0.9993416  0.0007137 -0.923 0.35606
## mom_sci -0.2263512  0.7974380  0.1220062 -1.855 0.06356
## 
## Likelihood ratio test=14.2  on 4 df, p=0.006684
## 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.