This document contains the meta-analysis of the surivival models, using time.age since baseline as the metric of interest and neuroticism by conscientiousness as the predictor of interest.

Code

The following packages were used to generate this table:

library(papaja)
library(tidyverse)
library(metafor)
library(knitr)
library(kableExtra)
library(here)

The files needed for this table are available at osf.io/mzfu9 in the Individual Study Output folder.

First we load the individual study analysis objects.

study.names = c("EAS", "HRS", "LBC", "LBLS", "MAP", "MAS", "MIDUS", "NAS", "OATS", "ROS", "SLS", "WLS")

lapply(here(paste0("mortality/study output/", study.names, "_survival_output.Rdata")), load, .GlobalEnv)
meta.data.time.age <- data.frame()
n <- 0

for(i in study.names){
  n <- n+1
  x <- get(paste0(i, "_survival_output"))
  meta.data.time.age[n, "study"] <- i
  meta.data.time.age[n, "coef"] <- x$age.interaction$coef["z.neur:z.con", "coef"]
  meta.data.time.age[n, "se"] <- x$age.interaction$coef["z.neur:z.con", "se(coef)"]
  meta.data.time.age[n, "n"] <- x$age.interaction$ntotal
  
  meta.data.time.age[n, "low.c_30.est"] <- x$age.interaction$simple.HR$low.c_30["z.neur", "coef"]
  meta.data.time.age[n, "low.c_30.se"] <- x$age.interaction$simple.HR$low.c_30["z.neur", "se(coef)"]
  meta.data.time.age[n, "hgh.c_30.est"] <- x$age.interaction$simple.HR$hgh.c_30["z.neur", "coef"]
  meta.data.time.age[n, "hgh.c_30.se"] <- x$age.interaction$simple.HR$hgh.c_30["z.neur", "se(coef)"]
  
  meta.data.time.age[n, "low.c_70.est"] <- x$age.interaction$simple.HR$low.c_70["z.neur", "coef"]
  meta.data.time.age[n, "low.c_70.se"] <- x$age.interaction$simple.HR$low.c_70["z.neur", "se(coef)"]
  meta.data.time.age[n, "hgh.c_70.est"] <- x$age.interaction$simple.HR$hgh.c_70["z.neur", "coef"]
  meta.data.time.age[n, "hgh.c_70.se"] <- x$age.interaction$simple.HR$hgh.c_70["z.neur", "se(coef)"]
  
  #calculate hazard ratios at low and high C
  low.c_30.HR <- exp(x$age.interaction$simple.HR$low.c_30["z.neur", "coef"])
  low.c_30.HR.lwr <- exp(x$age.interaction$simple.HR$low.c_30["z.neur", "coef"] - 1.96*x$age.interaction$simple.HR$low.c_30["z.neur", "se(coef)"])
  low.c_30.HR.upr <- exp(x$age.interaction$simple.HR$low.c_30["z.neur", "coef"] + 1.96*x$age.interaction$simple.HR$low.c_30["z.neur", "se(coef)"])
  meta.data.time.age[n, "low.c_30.CI"] <- paste0(printnum(low.c_30.HR)," [",printnum(low.c_30.HR.lwr),", ",printnum(low.c_30.HR.upr),"]")

  hgh.c_30.HR <- exp(x$age.interaction$simple.HR$hgh.c_30["z.neur", "coef"])
  hgh.c_30.HR.lwr <- exp(x$age.interaction$simple.HR$hgh.c_30["z.neur", "coef"] - 1.96*x$age.interaction$simple.HR$hgh.c_30["z.neur", "se(coef)"])
  hgh.c_30.HR.upr <- exp(x$age.interaction$simple.HR$hgh.c_30["z.neur", "coef"] + 1.96*x$age.interaction$simple.HR$hgh.c_30["z.neur", "se(coef)"])
  meta.data.time.age[n, "hgh.c_30.CI"] <- paste0(printnum(hgh.c_30.HR),
                                                  " [",printnum(hgh.c_30.HR.lwr),", ",
                                                  printnum(hgh.c_30.HR.upr),"]")
  
  low.c_70.HR <- exp(x$age.interaction$simple.HR$low.c_70["z.neur", "coef"])
  low.c_70.HR.lwr <- exp(x$age.interaction$simple.HR$low.c_70["z.neur", "coef"] - 1.96*x$age.interaction$simple.HR$low.c_70["z.neur", "se(coef)"])
  low.c_70.HR.upr <- exp(x$age.interaction$simple.HR$low.c_70["z.neur", "coef"] + 1.96*x$age.interaction$simple.HR$low.c_70["z.neur", "se(coef)"])
  meta.data.time.age[n, "low.c_70.CI"] <- paste0(printnum(low.c_70.HR)," [",printnum(low.c_70.HR.lwr),", ",printnum(low.c_70.HR.upr),"]")

  hgh.c_70.HR <- exp(x$age.interaction$simple.HR$hgh.c_70["z.neur", "coef"])
  hgh.c_70.HR.lwr <- exp(x$age.interaction$simple.HR$hgh.c_70["z.neur", "coef"] - 1.96*x$age.interaction$simple.HR$hgh.c_70["z.neur", "se(coef)"])
  hgh.c_70.HR.upr <- exp(x$age.interaction$simple.HR$hgh.c_70["z.neur", "coef"] + 1.96*x$age.interaction$simple.HR$hgh.c_70["z.neur", "se(coef)"])
  meta.data.time.age[n, "hgh.c_70.CI"] <- paste0(printnum(hgh.c_70.HR)," [",printnum(hgh.c_70.HR.lwr),", ",printnum(hgh.c_70.HR.upr),"]")

}

meta.results.time.age <- rma(yi = coef, 
                    sei = se, 
                    ni = n,  
                    measure = "RR",
                    slab = study, 
                    data = meta.data.time.age)
#meta-analyze HR at low and high c

meta.low30.time.age <- rma(yi = low.c_30.est, 
                    sei = low.c_30.se, 
                    ni = n,  
                    measure = "RR",
                    slab = study, 
                    data = meta.data.time.age)

meta.hgh30.time.age <- rma(yi = hgh.c_30.est, 
                    sei = hgh.c_30.se, 
                    ni = n,  
                    measure = "RR",
                    slab = study, 
                    data = meta.data.time.age)

meta.low70.time.age <- rma(yi = low.c_70.est, 
                    sei = low.c_70.se, 
                    ni = n,  
                    measure = "RR",
                    slab = study, 
                    data = meta.data.time.age)

meta.hgh70.time.age <- rma(yi = hgh.c_70.est, 
                    sei = hgh.c_70.se, 
                    ni = n,  
                    measure = "RR",
                    slab = study, 
                    data = meta.data.time.age)


#find plot limits
max.ci.mod1 = max(exp(meta.data.time.age$coef+1.96*meta.data.time.age$se))
min.ci.mod1 = min(exp(meta.data.time.age$coef-1.96*meta.data.time.age$se))
range.mod1 = max.ci.mod1-min.ci.mod1

lower.mod1 = min.ci.mod1-(range.mod1*7)
upper.mod1 = max.ci.mod1+(range.mod1*1.5)

cex.set = .7

#estimate position of extra information
pos.mod1 = min.ci.mod1-lower.mod1
pos.mod1 = pos.mod1/15
pos.mod1 = c(lower.mod1+3*pos.mod1,
             lower.mod1+5*pos.mod1,
             lower.mod1+7.5*pos.mod1,
             lower.mod1+10.5*pos.mod1,
             lower.mod1+13*pos.mod1)

#extract weighted values at -1 and +1 sd
meta_est_low.c_30_mod1 = printnum(exp(meta.low30.time.age$b))
meta_lwr_low.c_30_mod1 = printnum(exp(meta.low30.time.age$ci.lb))
meta_upr_low.c_30_mod1 = printnum(exp(meta.low30.time.age$ci.ub))
meta_low.c_30_ci_mod1 = paste0(meta_est_low.c_30_mod1, " [",
                           meta_lwr_low.c_30_mod1, ", ",
                           meta_upr_low.c_30_mod1, "]")

meta_est_hgh.c_30_mod1 = printnum(exp(meta.hgh30.time.age$b))
meta_lwr_hgh.c_30_mod1 = printnum(exp(meta.hgh30.time.age$ci.lb))
meta_upr_hgh.c_30_mod1 = printnum(exp(meta.hgh30.time.age$ci.ub))
meta_hgh.c_30_ci_mod1 = paste0(meta_est_hgh.c_30_mod1, " [",
                           meta_lwr_hgh.c_30_mod1, ", ",
                           meta_upr_hgh.c_30_mod1, "]")

meta_est_low.c_70_mod1 = printnum(exp(meta.low70.time.age$b))
meta_lwr_low.c_70_mod1 = printnum(exp(meta.low70.time.age$ci.lb))
meta_upr_low.c_70_mod1 = printnum(exp(meta.low70.time.age$ci.ub))
meta_low.c_70_ci_mod1 = paste0(meta_est_low.c_70_mod1, " [",
                           meta_lwr_low.c_70_mod1, ", ",
                           meta_upr_low.c_70_mod1, "]")

meta_est_hgh.c_70_mod1 = printnum(exp(meta.hgh70.time.age$b))
meta_lwr_hgh.c_70_mod1 = printnum(exp(meta.hgh70.time.age$ci.lb))
meta_upr_hgh.c_70_mod1 = printnum(exp(meta.hgh70.time.age$ci.ub))
meta_hgh.c_70_ci_mod1 = paste0(meta_est_hgh.c_70_mod1, " [",
                           meta_lwr_hgh.c_70_mod1, ", ",
                           meta_upr_hgh.c_70_mod1, "]")
#forest plot

meta.data.time.age$study = gsub("LBC", "LBC1936", meta.data.time.age$study)


forest(meta.results.time.age,
       xlim = c(lower.mod1, upper.mod1),
       cex = cex.set,
       slab = meta.data.time.age$study,
       transf=exp,
       refline = 1,
       ilab = meta.data.time.age[,c("n",
                                "low.c_30.CI",
                                "hgh.c_30.CI",
                                "low.c_70.CI",
                                "hgh.c_70.CI")],
       ilab.xpos = pos.mod1)

# #additional text
text(pos.mod1[1], length(study.names)+1.5,
     c("N"), cex = cex.set)
text(pos.mod1[2], length(study.names)+1.5,
     c("Low Con"), cex = cex.set)
text(pos.mod1[3], length(study.names)+1.5,
     c("High Con"), cex = cex.set)
text(pos.mod1[4], length(study.names)+1.5,
     c("Low Con"), cex = cex.set)
text(pos.mod1[5], length(study.names)+1.5,
     c("High Con"), cex = cex.set)
text((pos.mod1[2]+pos.mod1[3])/2, length(study.names)+2.5,
     c("Age 30"), cex = cex.set)
text((pos.mod1[4]+pos.mod1[5])/2, length(study.names)+2.5,
     c("Age 70"), cex = cex.set)

text((pos.mod1[3]+pos.mod1[4])/2, length(study.names)+3,
     "Neuroticism Hazard Ratio", cex = cex.set)

text(lower.mod1, length(study.names)+1.5,
     "Study", cex = cex.set, pos = 4)
text(upper.mod1, length(study.names)+1.5,
     "Interaction Hazard Ratio [95% CI]", cex = cex.set, pos=2)
text(pos.mod1[2:5], -1,
     c(meta_low.c_30_ci_mod1, meta_hgh.c_30_ci_mod1, 
       meta_low.c_70_ci_mod1, meta_hgh.c_70_ci_mod1), cex = cex.set)
text(lower.mod1, -2, pos=4, cex=cex.set,
     bquote(paste("(Q = ",
     .(formatC(meta.results.time.age$QE, digits=2, format="f")),
     ", df = ", .(meta.results.time.age$k - meta.results.time.age$p),
     ", p = ", .(formatC(meta.results.time.age$QEp, digits=2, format="f")),
     "; ", I^2, " = ",
     .(formatC(meta.results.time.age$I2, digits=1, format="f")), "%)")))