setwd("~/dynamicHRpaper/")
source("code/utils.R")
fig_path="./Figs/Fig2/"
df=readRDS("data/amit_df_with_pc_correction.rds")
dt=df[,c("prs.r","phenos.has_CAD","phenos.enrollment")]
names(dt)=c("x","y","phenos.enrollment")
my <- glm(formula = y~x,data = dt[which(dt$phenos.enrollment<55),],family="binomial")
newdf=data.frame(x=seq(0,1,by=0.01),y=rep(1,101),phenos.enrollment=rep(1,101))
preds <- predict(my, newdata = newdf, type = "response", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- preds$fit + (critval * preds$se.fit)
lwr <- preds$fit - (critval * preds$se.fit)
fit <- preds$fit
predy=data.frame(fit,upr,lwr)
rm(preds)
rm(my)
####################
mm <- glm(formula = y~ x #+ phenos.enrollment,
,data = dt[which(dt$phenos.enrollment<65&dt$phenos.enrollment>55),],family="binomial")
preds <- predict(mm, newdata = newdf, type = "response", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- preds$fit + (critval * preds$se.fit)
lwr <- preds$fit - (critval * preds$se.fit)
fit <- preds$fit
predm=data.frame(fit,upr,lwr)
rm(preds)
rm(mm)
#########
mo<- glm(formula = y~ x #+ phenos.enrollment,
,data = dt[which(dt$phenos.enrollment>65),],family="binomial")
preds <- predict(mo, newdata = newdf, type = "response", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- preds$fit + (critval * preds$se.fit)
lwr <- preds$fit - (critval * preds$se.fit)
fit <- preds$fit
predo=data.frame(fit,upr,lwr)
rm(preds)
rm(mo)
predf=rbind(predy,predm,predo)
## link scale
#predr=exp(predf)/(1+exp(predf))
predr=predf
predr$variable=c(rep("<55",101),rep("55-65",101),rep(">65",101))
predr$x=rep(newdf$x,3)
predr$variable= factor(
predr$variable,
levels = c("<55","55-65",">65"),labels = c(
"<55 years","55-65 years",">65 years"
)
)
predp=predr
gp=ggplot(data=predp, aes(x=x*100, y=fit, ymin=lwr, ymax=upr, color=variable,fill=variable)) +
geom_line(linewidth=2) +
geom_ribbon(alpha=0.3) +
xlab("Polygenic score percentile") +
ylab("Absolute risk of coronary artery disease")+scale_fill_futurama(name="Age")+scale_color_futurama(name="Age")+
# scale_fill_manual(name="Age Group",values=c("<55"="deepskyblue","55-65"="deepskyblue3",">65"="deepskyblue4"))+
# scale_color_manual(name="Age Group",values=c("<55 years"="deepskyblue","55-65 years"="deepskyblue3",">65 years"="deepskyblue4"))+
theme_classic(base_size =20)#+theme(legend.position = "null")
#ggsave(gp,filename = paste0(fig_path,"gprs_abs_withleg.tiff"),dpi=300,height=6,width = 6)
ga=ggplot(data=predp, aes(x=x*100, y=fit, ymin=lwr, ymax=upr, color=variable,fill=variable)) +
geom_line(linewidth=2) +
geom_ribbon(alpha=0.3) +
xlab("Polygenic score percentile") +
ylab("Absolute risk of coronary artery disease")+scale_fill_futurama(name="Age Group")+scale_color_futurama(name="Age")+guides(color="none")+
# scale_fill_manual(name="Age Group",values=c("<55"="deepskyblue","55-65"="deepskyblue3",">65"="deepskyblue4"))+
# scale_color_manual(name="Age Group",values=c("<55 years"="deepskyblue","55-65 years"="deepskyblue3",">65 years"="deepskyblue4"))+
theme_classic(base_size =20)+theme(legend.position = "null")
#ggsave(ga,filename = paste0(fig_path,"gprs_abs.tiff"),dpi=300,height=6,width = 6)
ggplotly(ga)
