# -------------------------------------------------------- #
# ASSOCIATION OF FALLS WITH LONG-TERM MORTALITY - R SCRIPT #
# -------------------------------------------------------- #

# -------------------------------------------------
# FIGURE 1: ADJUSTED NON-PARAMETRIC SURVIVAL CURVES

# READ TEXT FILE

figure1 <- read.csv("c:/.../Figure1.csv",
                    header=T,sep=",",dec=".",stringsAsFactors=F)
names(figure1) <- c("surv.cox","time",names(figure1)[3])
figure1 <- figure1[order(figure1$time),]

# FIGURE 1
# Adjusted non-parametric survival curves were obtained from a stratified Cox model at weighted mean values of covariates

cdf <- matrix(nrow=3,ncol=5,dimnames=list(c("no","non-severe","severe"),1:5))
for(i in 1:3){
   for(j in 1:5){
      cdf[i,j] <- 1 - min(figure1$surv.cox[figure1$fallcat==(i-1)&figure1$time<=j])
      }
   }

# tiff(filename="c:/.../Figure 1.tiff",
#      width=8.5,height=7.5,units="cm",pointsize=6,res=300)

pdf(file="c:/.../Figure 1.pdf",
    width=8.5/2.54,height=7.5/2.54,pointsize=6)

par(mfrow=c(1,1),mar=c(6.5,3.5,0,0),oma=c(0,0,0,0),tcl=-0.35)
plot(figure1$time,1-figure1$surv.cox,type="n",xlim=c(0,5),ylim=c(0,0.60),xlab="",ylab="",axes=F)
   lines(figure1$time[figure1$fallcat==0],1-figure1$surv.cox[figure1$fallcat==0],type="s",lty=1,lwd=1,col=gray(0.80))
   lines(figure1$time[figure1$fallcat==1],1-figure1$surv.cox[figure1$fallcat==1],type="s",lty=1,lwd=1,col=gray(0.50))
   lines(figure1$time[figure1$fallcat==2],1-figure1$surv.cox[figure1$fallcat==2],type="s",lty=1,lwd=1,col="black")
   axis(1,at=0:5,labels=0:5,cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   axis(2,at=seq(0,0.6,0.1),labels=seq(0,0.6,0.1)*100,cex.axis=0.85,font.axis=1,mgp=c(3,0.65,0),las=1,adj=1,lty=1,lwd=1)
   mtext(side=1,"Follow-up time (years)",cex=0.85,font=1,line=1.5)
   mtext(side=2,"Cumulative all-cause mortality (%)",cex=0.85,font=1,las=3,line=2.25)
   legend("bottomright",inset=0.025,bty="n",legend=c("No fall","Non-severe fall","Severe fall"),
          cex=0.85,lty=1,lwd=1,col=c(gray(0.80),gray(0.50),"black"),y.intersp=1.25)
   mtext(side=1,at=1:5,format(100*cdf[1,],digits=3),cex=0.85,font=1,line=2.75)
   mtext(side=1,at=-0.7,"No fall",cex=0.85,font=1,adj=0,line=2.75)
   mtext(side=1,at=1:5,format(100*cdf[2,],digits=3),cex=0.85,font=1,line=4)
   mtext(side=1,at=-0.7,"Non-severe fall",cex=0.85,font=1,adj=0,line=4)
   mtext(side=1,at=1:5,format(100*cdf[3,],digits=3),cex=0.85,font=1,line=5.25)
   mtext(side=1,at=-0.7,"Severe fall",cex=0.85,font=1,adj=0,line=5.25)

dev.off()

# -------------------------------------------------------------
# FIGURE 2: TIME-DEPENDENT EFFECT OF FALL SEVERITY ON MORTALITY

# READ TEXT FILE

figure2 <- read.csv("c:/.../Figure2.csv",
                    header=T,sep=",",dec=".",stringsAsFactors=F)
names(figure2) <- c("time","lnt","lnt.s1")
figure2 <- figure2[order(figure2$time),]

# SMOOTH EFFECT OF FALL SEVERITY ON MORTALITY AS A FUNCTION OF TIME

# Parameter estimates from a spline-based parametric survival model stratified by severity of fall
b1 <- c(-0.14121697,0.03668751,-0.04190252)
names(b1) <- c("cons","lnt","lnt.s1")
V1 <- matrix(c( 0.17739271,-0.0769419 , 0.02277862,
               -0.0769419 , 0.04834324,-0.01394953,
                0.02277862,-0.01394953, 0.01482787),
             nrow=3,dimnames=list(c("cons","lnt","lnt.s1"),c("cons","lnt","lnt.s1")))

b2 <- c(-0.06816066,0.28408735,0.28447667)
names(b2) <- c("cons","lnt","lnt.s1")
V2 <- matrix(c( 0.3114961 ,-0.24219956,-0.16416681,
               -0.24219956, 0.24694128, 0.14181629,
               -0.16416681, 0.14181629, 0.11852101),
             nrow=3,dimnames=list(c("cons","lnt","lnt.s1"),c("cons","lnt","lnt.s1")))

# FIGURE 2

# Natural cubic splines for log time with a single internal knot at the 50th percentile
mat <- cbind(1,figure2$lnt,figure2$lnt.s1)

# Time-dependent log cumulative hazard ratios for non-severe and severe fall compared with no fall
lnCHR1 <- as.numeric(mat%*%b1)
se.lnCHR1 <- sqrt(diag(mat%*%V1%*%t(mat)))
ll.lnCHR1 <- lnCHR1 - qnorm(0.975)*se.lnCHR1
ul.lnCHR1 <- lnCHR1 + qnorm(0.975)*se.lnCHR1

lnCHR2 <- as.numeric(mat%*%b2)
se.lnCHR2 <- sqrt(diag(mat%*%V2%*%t(mat)))
ll.lnCHR2 <- lnCHR2 - qnorm(0.975)*se.lnCHR2
ul.lnCHR2 <- lnCHR2 + qnorm(0.975)*se.lnCHR2

# Cumulative hazard ratios at 1 to 5 years of follow-up (spline interpolation)
CHR1 <- format(exp(spline(figure2$lnt,lnCHR1,xout=log(1:5),method="fmm")$y),digits=2)
ci.CHR1 <- paste("(",format(exp(spline(figure2$lnt,ll.lnCHR1,xout=log(1:5),method="fmm")$y),digits=2),"–",
                     format(exp(spline(figure2$lnt,ul.lnCHR1,xout=log(1:5),method="fmm")$y),digits=3),")",sep="")

CHR2 <- format(exp(spline(figure2$lnt,lnCHR2,xout=log(1:5),method="fmm")$y),digits=3)
ci.CHR2 <- paste("(",format(exp(spline(figure2$lnt,ll.lnCHR2,xout=log(1:5),method="fmm")$y),digits=2),"–",
                     format(exp(spline(figure2$lnt,ul.lnCHR2,xout=log(1:5),method="fmm")$y),digits=3),")",sep="")

# Figure 2
# tiff(filename="c:/.../Figure 2.tiff",
#      width=8.5,height=13.75,units="cm",pointsize=6,res=300)

pdf(file="c:/.../Figure 2.pdf",
    width=8.5/2.54,height=13.75/2.54,pointsize=6)

par(mfrow=c(2,1),mar=c(2,3.5,0.5,1.25),oma=c(5.5,0,0.25,0),tcl=-0.35)
plot(figure2$time,lnCHR1,type="n",xlim=c(0,5),ylim=log(c(0.25,4)),xlab="",ylab="",axes=F)
   lines(figure2$time[figure2$time>=0.5],lnCHR1[figure2$time>=0.5],lty=1,lwd=2)
   lines(rev(figure2$time[figure2$time>=0.5]),rev(ll.lnCHR1[figure2$time>=0.5]),lty=2,lwd=2)
   lines(rev(figure2$time[figure2$time>=0.5]),rev(ul.lnCHR1[figure2$time>=0.5]),lty=2,lwd=2)
   abline(h=0,lty=2,lwd=1)
   axis(1,at=0:5,labels=0:5,cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   axis(2,at=log(c(0.25,0.5,1,2,4)),labels=c(0.25,0.5,1,2,4),cex.axis=0.85,font.axis=1,mgp=c(3,0.65,0),las=1,adj=1,lty=1,lwd=1)
   mtext(side=2,"Cumulative hazard ratio for all-cause mortality",cex=0.85,font=1,las=3,line=2.25)
   mtext(side=3,"Non-severe fall vs. no fall",cex=0.85,font=1,line=-0.5)
   text(5,log(0.35),labels=expression(paste(italic(P)," for proportional hazards = 0.94")),cex=0.85,font=1,adj=1)
plot(figure2$time,lnCHR2,type="n",xlim=c(0,5),ylim=log(c(0.25,4)),xlab="",ylab="",axes=F)
   lines(figure2$time[figure2$time>=0.5],lnCHR2[figure2$time>=0.5],lty=1,lwd=2)
   lines(rev(figure2$time[figure2$time>=0.5]),rev(ll.lnCHR2[figure2$time>=0.5]),lty=2,lwd=2)
   lines(rev(figure2$time[figure2$time>=0.5]),rev(ul.lnCHR2[figure2$time>=0.5]),lty=2,lwd=2)
   abline(h=0,lty=2,lwd=1)
   axis(1,at=0:5,labels=0:5,cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   axis(2,at=log(c(0.25,0.5,1,2,4)),labels=c(0.25,0.5,1,2,4),cex.axis=0.85,font.axis=1,mgp=c(3,0.65,0),las=1,adj=1,lty=1,lwd=1)
   mtext(side=1,"Follow-up time (years)",cex=0.85,font=1,line=1.5)
   mtext(side=2,"Cumulative hazard ratio for all-cause mortality",cex=0.85,font=1,las=3,line=2.25)
   mtext(side=3,"Severe fall vs. no fall",cex=0.85,font=1,line=-0.5)
   text(5,log(0.35),labels=expression(paste(italic(P)," for proportional hazards = 0.70")),cex=0.85,font=1,adj=1)
   mtext(side=1,at=-0.75,"Non-severe vs. no",cex=0.85,font=1,adj=0,line=3)
   mtext(side=1,at=1:5,CHR1,cex=0.85,font=1,line=3)
   mtext(side=1,at=1:5,ci.CHR1,cex=0.85,font=1,line=4)
   mtext(side=1,at=-0.75,"Severe vs. no",cex=0.85,font=1,adj=0,line=5)
   mtext(side=1,at=1:5,CHR2,cex=0.85,font=1,line=5)
   mtext(side=1,at=1:5,ci.CHR2,cex=0.85,font=1,line=6)

dev.off()

# ----------------------------------------------------------
# FIGURE 3: STRATIFIED EFFECTS OF FALL SEVERITY ON MORTALITY

# EFFECTS OF FALL SEVERITY ON MORTALITY BY SUBGROUPS

var <- c("Age (years)",
         "Sex",
         "Type of facility",
         "Dementia",
         "No. of chronic conditions",
         "Use of antidepressants",
         "No. of prescribed medications",
         "Urinary incontinence",
         "Functional dependency",
         "Overall")

label <- c("65–84","", # See below
           "Women","Men",
           "Public/subsidized","Private",
           "No","Yes",
           "0–3","", # See below
           "No","Yes",
           "0–4","", # See below
           "No","Mild/severe",
           "Independent","Mild/moderate/severe")

# Parameter estimates and SEs for severity of fall obtained from Cox models with interactions
b1 <- c( 0.1070561,-0.2218006,
        -0.2103577, 0.3407182,
        -0.4670683, 0.5254737,
        -0.0519604,-0.1292472,
        -0.1376056,-0.0741748,
        -0.3034225, 0.5778953,
        -0.1949541,-0.017002 ,
        -0.127971 ,-0.1318905,
        -1.135898 ,-0.0714127,
        -0.0847297) # Overall
se1 <- c(0.3055091,0.3336101,
         0.2672442,0.3583815,
         0.2739167,0.3454897,
         0.2722872,0.3490546,
         0.3165214,0.2872005,
         0.2703734,0.4824233,
         0.3864045,0.245663 ,
         0.475844 ,0.2399392,
         0.8491262,0.232007 ,
         0.2269232) # Overall

b2 <- c( 0.0945624, 0.4210839,
         0.1690932, 0.6913788,
         0.2555455, 0.3103749,
         0.2237501, 0.4164196,
         0.6269168, 0.043405 ,
         0.3927113,-0.2806308,
         0.6986755,-0.0958862,
         0.0902771, 0.3088016,
         0.7415685, 0.2471979,
         0.3079675) # Overall
se2 <- c(0.5099837,0.457834 ,
         0.4253105,0.2458494,
         0.401406 ,0.4555522,
         0.4823074,0.3389422,
         0.2599456,0.5170424,
         0.33437  ,1.15862  ,
         0.2810515,0.5942471,
         0.6067363,0.3913703,
         0.1680568,0.3232823,
         0.3070507) # Overall

p.int <- c(0.6996,
           0.2735,
           0.0691,
           0.9436,
           0.6340,
           0.2541,
           0.3417,
           0.9530,
           0.0757)

HR1 <- paste(format(exp(b1),digits=2)," (",
             format(exp(b1-qt(0.975,41)*se1),digits=1),"–",
             format(exp(b1+qt(0.975,41)*se1),digits=3),")",sep="")

HR2 <- paste(format(exp(b2),digits=2)," (",
             format(exp(b2-qt(0.975,41)*se2),digits=1),"–",
             format(exp(b2+qt(0.975,41)*se2),digits=3),")",sep="")

# FIGURE 3

# tiff(filename="c:/.../Figure 3.tiff",
#      width=12.5,height=12.5,units="cm",pointsize=6,res=300)

pdf(file="c:/.../Figure 3.pdf",
    width=12.5/2.54,height=12.5/2.54,pointsize=6)

par(mfrow=c(1,2),mar=c(1.75,6.25,2.75,0),oma=c(0,10,0,0.25),tcl=-0.35)
plot(b1,c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),type="n",xlim=log(c(0.25,4)),ylim=c(0,36),xlab="",ylab="",axes=F)
   symbols(b1[-19],c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-19],
           squares=1/se1[-19],inches=0.09*max(1/se1[-19])/max(1/se1[-19],1/se2[-19]),add=T,lty=1,lwd=1,bg="black")
   segments((b1-qt(0.975,41)*se1)[-c(12,17,19)],c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(12,17,19)],
            (b1+qt(0.975,41)*se1)[-c(12,17,19)],c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(12,17,19)],lty=1,lwd=1.5)
   arrows(log(0.25),4,(b1+qt(0.975,41)*se1)[17],4,code=1,length=0.035,angle=30,lty=1,lwd=1.5)
   arrows((b1-qt(0.975,41)*se1)[12],15,log(4),15,code=2,length=0.035,angle=30,lty=1,lwd=1.5)
   polygon(c((b1-qt(0.975,41)*se1)[19],b1[19],(b1+qt(0.975,41)*se1)[19],b1[19]),
           c(0,-0.50,0,0.50),border=T,density=-1,lty=1,lwd=1,col=1)
   segments(b1[19],0.50,b1[19],36,lty=2,lwd=1)
   abline(v=0,lty=1,lwd=1)
   axis(1,at=log(c(0.25,0.5,1,2,4)),labels=c(0.25,0.5,1,2,4),cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   mtext(side=2,at=c(37,33,29,25,21,17,13,9,5,0),var,cex=0.85,font=1,las=1,adj=0,line=15.75)
   mtext(side=2,at=38,"Subgroup",cex=0.85,font=1,las=1,adj=0,line=15.75)
   mtext(side=2,at=c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-19],label,cex=0.85,font=1,las=1,adj=0,line=14.75)
   mtext(side=2,at=c(35,19,11),expression("">=85,"">=4,"">=5),cex=0.85,font=1,las=1,adj=0,line=14.85)
   mtext(side=2,at=c(34,30,26,22,18,14,10,6,2)-0.05,expression(paste(italic(P)," for heterogeneity = ")),cex=0.85,font=1,las=1,adj=1,line=7.90)
   mtext(side=2,at=c(34,30,26,22,18,14,10,6,2),format(p.int,digits=1),cex=0.85,font=1,las=1,adj=0,line=7.90)
   mtext(side=2,at=c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),HR1,cex=0.85,font=1,las=1,adj=0,line=5.25)
   mtext(side=2,at=38,"Non-severe fall vs. no fall",cex=0.85,font=1,las=1,adj=0.5,line=-3.75)
   mtext(side=2,at=39,"Hazard ratio for all-cause mortality (95% CI)",cex=0.85,font=1,las=1,adj=0.5,line=-13.50)
plot(b2,c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),type="n",xlim=log(c(0.25,4)),ylim=c(0,36),xlab="",ylab="",axes=F)
   symbols(b2[-19],c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-19],
           squares=1/se2[-19],inches=0.09*max(1/se2[-19])/max(1/se1[-19],1/se2[-19]),add=T,lty=1,lwd=1,bg="black")
   segments((b2-qt(0.975,41)*se2)[-c(12,19)],c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(12,19)],
            (b2+qt(0.975,41)*se2)[-c(12,19)],c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(12,19)],lty=1,lwd=1.5)
   arrows(log(0.25),15,log(4),15,code=3,length=0.035,angle=30,lty=1,lwd=1.5)
   polygon(c((b2-qt(0.975,41)*se2)[19],b2[19],(b2+qt(0.975,41)*se2)[19],b2[19]),
           c(0,-0.50,0,0.50),border=T,density=-1,lty=1,lwd=1,col=1)
   segments(b2[19],0.50,b2[19],36,lty=2,lwd=1)
   abline(v=0,lty=1,lwd=1)
   axis(1,at=log(c(0.25,0.5,1,2,4)),labels=c(0.25,0.5,1,2,4),cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   mtext(side=2,at=c(36,35,32,31,28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),HR2,cex=0.85,font=1,las=1,adj=0,line=5.25)
   mtext(side=2,at=38,"Severe fall vs. no fall",cex=0.85,font=1,las=1,adj=0.5,line=-3.75)

dev.off()