Figure 3

Author

Sarah Urbut

Published

March 7, 2024

Cumulative incidence plot

Age cat specific plots

interaction

Here we break into three age groups and stratify by PRS, asking about resolution by PRS:

Code

dyoung=df[df$phenos.enrollment<55,]
dyoung$pce.age=scale(dyoung$ascvd_10y_accaha_all)
dyoung$pce.r=rank(dyoung$pce.age)/length(dyoung$pce.age)
dyoung$pce.cat=cut(dyoung$pce.r,breaks=c(0,0.20,0.80,1))
dyoung$Interaction=interaction(dyoung$prscat,dyoung$pce.cat)

## now with interaction

fit <- survfit(Surv(phenos.CAD_censor_age,phenos.has_CAD)~Interaction, data=dyoung)
 #fit <- survfit(Surv((phenos.CAD_censor_age), phenos.has_CAD) ~prscat, data=df[df$phenos.enrollment<55,])
g_young=ggsurvplot(fit,
           pval = TRUE, conf.int = TRUE,
        
           fun="cumhaz",
           censor=FALSE,
          
           xlim = c(40, 67),
           ylim=c(0,0.4),
           
            #xlab = "Time since Enrollment",
           xlab = "Time of Event (Age)",
           ylab = "Cumulative Hazard of Incident CAD",
           #ylim=c(0,0.15),
           title = "Enrollment Age <55 yrs",
          legend.title = "Genomic and Clinical Risk Strata",
          ggtheme = theme_classic2(base_size=18),
palette=c(colorRampPalette(c("#DEEBF7", "#9ECAE1", "#3182BD"))(3),
colorRampPalette(c("#E5F5E0", "#A1D99B", "#31A354"))(3),
colorRampPalette(c("#FEE6CE", "#FDAE6B", "#E6550D"))(3)),break.x.by = 5,

          
          legend.labs = c("Low PCE, Low PRS",
                          "Low PCE, Intermediate PRS" ,
                          "Low PCE, High PRS",
                          
                          "Intermediate PCE, Low PRS",
                          "Intermediate PCE, Intermediate PRS",
                          "Intermediate PCE, High PRS",
                          
                          "High PCE, Low PRS",
                          "High PCE, Intermediate PRS",
                          "High PCE, High PRS"))

g_young

Code

###

dmid=df[df$phenos.enrollment>55&df$phenos.enrollment<65,]
dmid$pce.age=scale(dmid$ascvd_10y_accaha_all)
dmid$pce.r=rank(dmid$pce.age)/length(dmid$pce.age)
dmid$pce.cat=cut(dmid$pce.r,breaks=c(0,0.20,0.80,1))
dmid$Interaction=interaction(dmid$prscat,dmid$pce.cat)
### interaction

fit <- survfit(Surv(phenos.CAD_censor_age,phenos.has_CAD)~Interaction, data=dmid)


g_mid=ggsurvplot(fit,
           pval = TRUE, 
           conf.int = TRUE,
           censor=FALSE,
           ##risk.table = F, 
           ##risk.table.col = "strata",
           fun="cumhaz",
           #xlim = c(0, 15), 
           xlim = c(55, 80),
           ylim=c(0,0.4),
          #xlab = "Time since Enrollment",
           xlab = "Time of Event (Age)",
           ylab = "Cumulative Hazard",
           #ylim=c(0,0.15),
           title = "Enrollment Age 55-65 yrs.",
          legend.title = "PRS:PCE Rank",
          ggtheme = theme_classic2(base_size=18),
           #font.family = "Arial",
               #  font.family = "Arial",
          break.x.by = 5,
palette=c(colorRampPalette(c("#DEEBF7", "#9ECAE1", "#3182BD"))(3),
colorRampPalette(c("#E5F5E0", "#A1D99B", "#31A354"))(3),
colorRampPalette(c("#FEE6CE", "#FDAE6B", "#E6550D"))(3)),

          
    
          legend.labs = c("Low PCE, Low PRS",
                          "Low PCE, Intermediate PRS" ,
                          "Low PCE, High PRS",
                          
                          "Intermediate PCE, Low PRS",
                          "PCE Intermediate, Intermediate PRS",
                          "Intermediate PCE, High PRS",
                          
                          "High PCE, Low PRS",
                          "High PCE, Intermediate PRS",
                          "High PCE, High PRS"))
g_mid

Code

###

Code

dold=df[df$phenos.enrollment>65,]
dold$pce.age=scale(dold$ascvd_10y_accaha_all)
dold$pce.r=rank(dold$pce.age)/length(dold$pce.age)
dold$pce.cat=cut(dold$pce.r,breaks=c(0,0.20,0.80,1))
dold$Interaction=interaction(dold$prscat,dold$pce.cat)

### Interaction


fit <- survfit(Surv(phenos.CAD_censor_age, phenos.has_CAD)~Interaction, data=dold) 


g_old=ggsurvplot(fit,
           pval = TRUE, conf.int = TRUE,censor=FALSE,
           ##risk.table = F, 
           ##risk.table.col = "strata",
           fun="cumhaz",
           #xlim = c(0, 15), 
           xlim = c(65, 85),
           ylim=c(0,0.4),
            #xlab = "Time since Enrollment",
           xlab = "Time of Event (Age)",
           #ylab = "Cumulative Hazard",
           #ylim=c(0,0.15),
           title = "Enrollment Age >65 yrs",
          legend.title = "PRS:PCE Rank",
          ggtheme = theme_classic2(base_size=18),
           #font.family = "Arial",
           palette=c(colorRampPalette(c("#DEEBF7", "#9ECAE1", "#3182BD"))(3),
colorRampPalette(c("#E5F5E0", "#A1D99B", "#31A354"))(3),
colorRampPalette(c("#FEE6CE", "#FDAE6B", "#E6550D"))(3)),
break.x.by = 5,
       
          legend.labs = c("Low PCE:Low PRS",
                          "Low PCE:Intermediate PRS" ,
                          "Low PCE:High PRS",
                          
                          "Intermediate PCE:Low PRS",
                          "PCE Intermediate:Intermediate PRS",
                          "Intermediate PCE:High PRS",
                          
                          "High PCE:Low PRS",
                          "High PCE:Intermediate PRS",
                          "High PCE:High PRS"))

g_old

alltogether

Code

gg1=ggarrange(g_young$plot,g_mid$plot+labs(y=NULL),g_old$plot+labs(y=NULL),nrow=1,ncol=3,labels=c("A.","B.","C."),common.legend = T,legend="bottom")

gg1

Code

ggsave(gg1,dpi=300,filename = "Figs/Fig3/interaction_time_byage_withPCS.tiff",height=10,width = 20,create.dir = TRUE)

return cumulative hazard for low and high of each age

Code

get_cumI=function(ggsurvplot,stratum){
  stratum_data <- ggsurvplot$data.survplot[ggsurvplot$data.survplot$strata == stratum, ]
  cumulative_hazard <- -log(stratum_data$surv)
  cum_haz_se=stratum_data$std.err/stratum_data$surv
  return(list(ch=cumulative_hazard[length(cumulative_hazard)],se=(stratum_data$std.err/stratum_data$surv)[length(cumulative_hazard)]))
}



get_cumI(g_young,stratum="Interaction=low.(0,0.2]")

$ch
[1] 0.004525493

$se
[1] 0.001113627

Code

get_cumI(g_young,stratum="Interaction=high.(0.8,1]")

$ch
[1] 0.1416114

$se
[1] 0.008012507

Code

get_cumI(g_young,stratum="Interaction=high.(0.8,1]")$ch-get_cumI(g_young,stratum="Interaction=high.(0.8,1]")$se*1.96

[1] 0.1259069

Code

get_cumI(g_young,stratum="Interaction=high.(0.8,1]")$ch+get_cumI(g_young,stratum="Interaction=high.(0.8,1]")$se*1.96

[1] 0.1573159

Code

get_cumI(g_mid,stratum="Interaction=low.(0,0.2]")

$ch
[1] 0.008308347

$se
[1] 0.002048986

Code

get_cumI(g_mid,stratum="Interaction=high.(0.8,1]")

$ch
[1] 0.2710532

$se
[1] 0.07265383

Code

get_cumI(g_old,stratum="Interaction=low.(0,0.2]")

$ch
[1] 0.04624608

$se
[1] 0.02588439

Code

get_cumI(g_old,stratum="Interaction=low.(0,0.2]")$ch-get_cumI(g_old,stratum="Interaction=low.(0,0.2]")$se*1

[1] 0.02036169

Code

get_cumI(g_old,stratum="Interaction=low.(0,0.2]")$ch+get_cumI(g_old,stratum="Interaction=low.(0,0.2]")$se*1

[1] 0.07213047

Code

get_cumI(g_old,stratum="Interaction=low.(0,0.2]")

$ch
[1] 0.04624608

$se
[1] 0.02588439

Code

get_cumI(g_old,stratum="Interaction=high.(0.8,1]")$ch-get_cumI(g_old,stratum="Interaction=high.(0.8,1]")$se*1

[1] 0.1797599

Code

get_cumI(g_old,stratum="Interaction=high.(0.8,1]")$ch+get_cumI(g_old,stratum="Interaction=high.(0.8,1]")$se*1

[1] 0.5728429

Code

get_cumI=function(data,stratum){
  fit <- survfit(Surv(phenos.CAD_censor_age,phenos.has_CAD)~1, data=data[data$Interaction==stratum])
  cumulative_hazard <- fit$cumhaz[length(fit$cumhaz)]
  cum_haz_se=fit$std.err[length(fit$cumhaz)]
  return(list(ch=cumulative_hazard,se=cum_haz_se))
}


get_cumI(dyoung,stratum="low.(0,0.2]")

$ch
[1] 0.004524879

$se
[1] 0.001108599

Code

get_cumI(dyoung,stratum="high.(0.8,1]")

$ch
[1] 0.1415855

$se
[1] 0.006954524

Code

get_cumI(dmid,stratum="low.(0,0.2]")

$ch
[1] 0.008306283

$se
[1] 0.002032033

Code

get_cumI(dmid,stratum="high.(0.8,1]")

$ch
[1] 0.2695373

$se
[1] 0.05540406

Code

get_cumI(dold,stratum="low.(0,0.2]")

$ch
[1] 0.04594303

$se
[1] 0.0247146

Code

get_cumI(dold,stratum="high.(0.8,1]")

$ch
[1] 0.3675943

$se
[1] 0.1349052

Code

## to get standard erro

# https://dominicmagirr.github.io/post/2022-01-18-be-careful-with-standard-errors-in-survival-survfit/

### http://www.sthda.com/english/wiki/survival-analysis-basics


# Get unique levels of the strata variable
strata_levels <- unique(dyoung$Interaction)

# Loop through each level of the strata
for (strata_level in c("low.(0,0.2]","high.(0.8,1]")) {
  # Subset the survival object for the current strata level
  
strata_surv <- survfit(Surv(phenos.CAD_censor_age,phenos.has_CAD)~1, data=dyoung[dyoung$Interaction==strata_level,])
  
  # Calculate the cumulative hazard
  cumhaz <- -log(strata_surv$surv)[length(strata_surv$surv)]
  lower=-log(strata_surv$upper)[length(strata_surv$surv)]
  upper=-log(strata_surv$lower)[length(strata_surv$surv)]
  

  # Print the results for the current strata level
  cat("Strata Level:", strata_level, "\n")
  cat("Cumulative Hazard:", cumhaz, "\n")
  cat("Cumulative Hazard Confidence Interval:", c(lower,upper), "\n\n")
}

Strata Level: low.(0,0.2] 
Cumulative Hazard: 0.004525493 
Cumulative Hazard Confidence Interval: 0.002352679 0.006698307 

Strata Level: high.(0.8,1] 
Cumulative Hazard: 0.1416114 
Cumulative Hazard Confidence Interval: 0.1279808 0.155242

Code

# Loop through each level of the strata
for (strata_level in c("low.(0,0.2]","high.(0.8,1]")) {
  # Subset the survival object for the current strata level
  
strata_surv <- survfit(Surv(phenos.CAD_censor_age,phenos.has_CAD)~1, data=dold[dold$Interaction==strata_level,])
  
  # Calculate the cumulative hazard
  cumhaz <- -log(strata_surv$surv)[length(strata_surv$surv)]
  fit=strata_surv
  lower=fit$cumhaz[length(fit$cumhaz)]-1.96*fit$std.chaz[length(fit$cumhaz)]
  upper=fit$cumhaz[length(fit$cumhaz)]+1.96*fit$std.chaz[length(fit$cumhaz)]
  lower2=-log(strata_surv$upper)[length(strata_surv$surv)]
  upper2=-log(strata_surv$lower)[length(strata_surv$surv)]
  

  # Print the results for the current strata level
  cat("Strata Level:", strata_level, "\n")
  cat("Cumulative Hazard:", cumhaz, "\n")
  cat("Cumulative Hazard Confidence Interval symmetric:", c(lower,upper), "\n\n")
  cat("Cumulative Hazard Confidence Interval log:", c(lower2,upper2), "\n\n")
}

Strata Level: low.(0,0.2] 
Cumulative Hazard: 0.04624608 
Cumulative Hazard Confidence Interval symmetric: -0.001941962 0.09382802 

Cumulative Hazard Confidence Interval log: 0 0.0946858 

Strata Level: high.(0.8,1] 
Cumulative Hazard: 0.3763014 
Cumulative Hazard Confidence Interval symmetric: 0.1199381 0.6152506 

Cumulative Hazard Confidence Interval log: 0.1118921 0.6407107

Code

#Among individuals with CAD <55 years, XXX (XXX%) had high (>20th # percentile) PRS and XXX (XXX%) had high (>20% estimated risk) PCE. 

dyoungevent=df[df$phenos.CAD_censor_age<55&df$phenos.has_CAD==1,]
sum(dyoungevent$prs.r>0.80)

[1] 429

Code

t.test(dyoungevent$prs.r>0.80)


    One Sample t-test

data:  dyoungevent$prs.r > 0.8
t = 26.625, df = 1084, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.3662531 0.4245304
sample estimates:
mean of x 
0.3953917

Code

#0.3935484 
sum(dyoungevent$ascvd_10y_accaha_all>20)

[1] 32

Code

t.test(dyoungevent$ascvd_10y_accaha_all>20)


    One Sample t-test

data:  dyoungevent$ascvd_10y_accaha_all > 20
t = 5.7395, df = 1084, p-value = 1.232e-08
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.01941036 0.03957582
sample estimates:
 mean of x 
0.02949309

Code

## 0.02949309