# ------------------------------------------ #
# SOCIAL ENGAGEMENT AND MORTALITY - R SCRIPT #
# ------------------------------------------ #

# READ TEXT FILES

# Combined weights and results from standardized Kaplan-Meier methods and spline-based survival models
survival <- read.csv("c:/.../survival.csv",header=T,sep=",",dec=".",stringsAsFactors=F)
names(survival) <- make.names(names(survival),allow_=F)
survival <- survival[order(survival$time),]

# SUPPLEMENTARY FIGURE 1
# Distributions of combined weights based on the three models

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

par(mfrow=c(1,1),mar=c(1.75,2.75,0,0),tcl=-0.35)
boxplot(survival$cweight1,survival$cweight2,survival$cweight3,xlim=c(1-1/3,3+1/3),ylim=c(0,4.1),xlab="",ylab="",axes=F,
        range=1.5,outline=T,boxwex=2/3,boxlty=1,boxlwd=1,medlty=1,medlwd=2,whisklty=1,whisklwd=1,staplelty=1,staplelwd=1,outpch=1,outlwd=1)
   points(1:3,c(mean(survival$cweight1),mean(survival$cweight2),mean(survival$cweight3)),pch=3,lwd=1)
   axis(1,at=1:3,labels=paste("Model",1:3,sep=" "),cex.axis=0.85,font.axis=1,mgp=c(3,0.50,0),lty=1,lwd=1)
   axis(2,at=0:4,labels=0: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,"Combined weight",cex=0.85,font=1,las=3,line=1.5)

dev.off()

# FIGURE 1: FULLY-STANDARDIZED NONPARAMETRIC AND SMOOTH CUMULATIVE MORTALITY CURVES BY BASELINE LEVEL OF SOCIAL ENGAGEMENT

# Fully-standardized smooth cumulative mortality curves and their 95% CIs for each level of social engagement
survival$cm.spl.std3    <- 1 - exp(-exp(survival$lnch.spl.std3))
survival$ll.cm.spl.std3 <- 1 - exp(-exp(survival$lnch.spl.std3 - qnorm(0.975)*survival$se.lnch.spl.std3))
survival$ul.cm.spl.std3 <- 1 - exp(-exp(survival$lnch.spl.std3 + qnorm(0.975)*survival$se.lnch.spl.std3))

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

par(mfrow=c(1,1),mar=c(2.60,3.25,0,0),tcl=-0.35)
plot(survival$time,1-survival$s.km.std3,type="n",xlim=c(0,10),ylim=c(0,0.75),xlab="",ylab="",axes=F)
   lines(survival$time[survival$socpart==0],1-survival$s.km.std3[survival$socpart==0],
         type="s",lty=1,lwd=0.75,col="red3")
   lines(survival$time[survival$socpart==1],1-survival$s.km.std3[survival$socpart==1],
         type="s",lty=1,lwd=0.75,col="green2")
   lines(survival$time[survival$socpart==2],1-survival$s.km.std3[survival$socpart==2],
         type="s",lty=1,lwd=0.75,col="deepskyblue2")
   lines(survival$time[survival$socpart==0],survival$cm.spl.std3[survival$socpart==0],
         type="l",lty=1,lwd=1.5,col="red3")
   lines(survival$time[survival$socpart==1],survival$cm.spl.std3[survival$socpart==1],
         type="l",lty=1,lwd=1.5,col="green2")
   lines(survival$time[survival$socpart==2],survival$cm.spl.std3[survival$socpart==2],
         type="l",lty=1,lwd=1.5,col="deepskyblue2")
   axis(1,at=seq(0,10,2),labels=seq(0,10,2),cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   axis(2,at=seq(0,0.75,0.25),labels=seq(0,0.75,0.25)*100,cex.axis=0.85,font.axis=1,las=1,adj=1,mgp=c(3,0.65,0),lty=1,lwd=1)
   mtext(side=1,"Follow-up (years)",cex=0.85,font=1,line=1.35)
   mtext(side=2,"Cumulative all-cause mortality (%)",cex=0.85,font=1,las=3,line=2)
   legend("topleft",inset=0.02,bty="n",y.intersp=1.2,
          legend=c("Level of social engagement at baseline","Low/null","Moderate","High"),
          lty=1,lwd=1.5,col=c("white","red3","green2","deepskyblue2"),cex=0.85)

dev.off()

# TABLE 2: STANDARDIZED DIFFERENCES IN CUMULATIVE MORTALITY AT 2, 5, AND 10 YEARS OF FOLLOW-UP BY BASELINE LEVEL OF SOCIAL ENGAGEMENT

# Log cumulative hazards and their SEs at 2, 5, and 10 years for each level of social engagement (spline interpolation from values at observed times)
lnch.std1 <- matrix(nrow=3,ncol=3,dimnames=list(c("t2","t5","t10"),c("socpart0","socpart1","socpart2")))
se.lnch.std1 <- lnch.std1
for(i in 0:2){
   lnch.std1[,i+1] <- spline(log(survival$time[survival$socpart==i]),survival$lnch.spl.std1[survival$socpart==i],xout=log(c(2,5,10)),method="fmm")$y
   se.lnch.std1[,i+1] <- spline(log(survival$time[survival$socpart==i]),survival$se.lnch.spl.std1[survival$socpart==i],xout=log(c(2,5,10)),method="fmm")$y
   }

lnch.std2 <- matrix(nrow=3,ncol=3,dimnames=list(c("t2","t5","t10"),c("socpart0","socpart1","socpart2")))
se.lnch.std2 <- lnch.std2
for(i in 0:2){
   lnch.std2[,i+1] <- spline(log(survival$time[survival$socpart==i]),survival$lnch.spl.std2[survival$socpart==i],xout=log(c(2,5,10)),method="fmm")$y
   se.lnch.std2[,i+1] <- spline(log(survival$time[survival$socpart==i]),survival$se.lnch.spl.std2[survival$socpart==i],xout=log(c(2,5,10)),method="fmm")$y
   }

lnch.std3 <- matrix(nrow=3,ncol=3,dimnames=list(c("t2","t5","t10"),c("socpart0","socpart1","socpart2")))
se.lnch.std3 <- lnch.std3
for(i in 0:2){
   lnch.std3[,i+1] <- spline(log(survival$time[survival$socpart==i]),survival$lnch.spl.std3[survival$socpart==i],xout=log(c(2,5,10)),method="fmm")$y
   se.lnch.std3[,i+1] <- spline(log(survival$time[survival$socpart==i]),survival$se.lnch.spl.std3[survival$socpart==i],xout=log(c(2,5,10)),method="fmm")$y
   }

# Confidence intervals for the differences in cumulative mortality by delta methods
cm.std1 <- 1 - exp(-exp(lnch.std1))
var.cm.std1 <- (exp(-exp(lnch.std1))*exp(lnch.std1)*se.lnch.std1)^2
matrix(paste(format(100*(cm.std1[,-1] - cm.std1[,1]),digits=0,nsmall=1,trim=T)," (",
             format(100*(cm.std1[,-1] - cm.std1[,1] - qnorm(0.975)*sqrt(var.cm.std1[,-1]+var.cm.std1[,1])),digits=1,nsmall=1,trim=T)," to ",
             format(100*(cm.std1[,-1] - cm.std1[,1] + qnorm(0.975)*sqrt(var.cm.std1[,-1]+var.cm.std1[,1])),digits=1,nsmall=1,trim=T),")",sep=""),nrow=3)

cm.std2 <- 1 - exp(-exp(lnch.std2))
var.cm.std2 <- (exp(-exp(lnch.std2))*exp(lnch.std2)*se.lnch.std2)^2
matrix(paste(format(100*(cm.std2[,-1] - cm.std2[,1]),digits=0,nsmall=1,trim=T)," (",
             format(100*(cm.std2[,-1] - cm.std2[,1] - qnorm(0.975)*sqrt(var.cm.std2[,-1]+var.cm.std2[,1])),digits=1,nsmall=1,trim=T)," to ",
             format(100*(cm.std2[,-1] - cm.std2[,1] + qnorm(0.975)*sqrt(var.cm.std2[,-1]+var.cm.std2[,1])),digits=1,nsmall=1,trim=T),")",sep=""),nrow=3)

cm.std3 <- 1 - exp(-exp(lnch.std3))
var.cm.std3 <- (exp(-exp(lnch.std3))*exp(lnch.std3)*se.lnch.std3)^2
matrix(paste(format(100*(cm.std3[,-1] - cm.std3[,1]),digits=0,nsmall=1,trim=T)," (",
             format(100*(cm.std3[,-1] - cm.std3[,1] - qnorm(0.975)*sqrt(var.cm.std3[,-1]+var.cm.std3[,1])),digits=1,nsmall=1,trim=T)," to ",
             format(100*(cm.std3[,-1] - cm.std3[,1] + qnorm(0.975)*sqrt(var.cm.std3[,-1]+var.cm.std3[,1])),digits=1,nsmall=1,trim=T),")",sep=""),nrow=3)

# STANDARDIZED DIFFERENCES IN MEDIAN SURVIVAL TIMES BY BASELINE LEVEL OF SOCIAL ENGAGEMENT

# Median survival times and their SEs for each level of social engagement
medt.std1 <- c(5.213409,5.769626,8.719002)
names(medt.std1) <- c("socpart0","socpart1","socpart2")
se.medt.std1 <- c(0.5811262,0.6167421,0.8345822)
names(se.medt.std1) <- c("socpart0","socpart1","socpart2")

medt.std2 <- c(5.118691,5.779498,8.544968)
names(medt.std2) <- c("socpart0","socpart1","socpart2")
se.medt.std2 <- c(0.5495559,0.6218628,0.8744135)
names(se.medt.std2) <- c("socpart0","socpart1","socpart2")

medt.std3 <- c(5.422458,5.838873,8.418187)
names(medt.std3) <- c("socpart0","socpart1","socpart2")
se.medt.std3 <- c(0.6329756,0.6448675,0.9399465)
names(se.medt.std3) <- c("socpart0","socpart1","socpart2")

# Confidence intervals for the differences in median survival times
paste(format(medt.std1[-1] - medt.std1[1],digits=1,trim=T)," (",
      format(medt.std1[-1] - medt.std1[1] - qnorm(0.975)*sqrt(se.medt.std1[-1]^2+se.medt.std1[1]^2),digits=0,nsmall=1,trim=T)," to ",
      format(medt.std1[-1] - medt.std1[1] + qnorm(0.975)*sqrt(se.medt.std1[-1]^2+se.medt.std1[1]^2),digits=0,nsmall=1,trim=T),")",sep="")

paste(format(medt.std2[-1] - medt.std2[1],digits=1,trim=T)," (",
      format(medt.std2[-1] - medt.std2[1] - qnorm(0.975)*sqrt(se.medt.std2[-1]^2+se.medt.std2[1]^2),digits=0,nsmall=1,trim=T)," to ",
      format(medt.std2[-1] - medt.std2[1] + qnorm(0.975)*sqrt(se.medt.std2[-1]^2+se.medt.std2[1]^2),digits=0,nsmall=1,trim=T),")",sep="")

paste(format(medt.std3[-1] - medt.std3[1],digits=1,trim=T)," (",
      format(medt.std3[-1] - medt.std3[1] - qnorm(0.975)*sqrt(se.medt.std3[-1]^2+se.medt.std3[1]^2),digits=0,nsmall=1,trim=T)," to ",
      format(medt.std3[-1] - medt.std3[1] + qnorm(0.975)*sqrt(se.medt.std3[-1]^2+se.medt.std3[1]^2),digits=0,nsmall=1,trim=T),")",sep="")

# FIGURE 2: STANDARDIZED DIFFERENCES IN 5-YEAR CUMULATIVE MORTALITY AMONG BASELINE LEVELS OF SOCIAL ENGAGEMENT BY SUBGROUP

# Log cumulative hazards and their SEs at 5 years of follow-up for each baseline level of social engagement and subgroup
lnch.all <- c(-0.4659166,-0.5358698,-1.103881)
names(lnch.all) <- c("socpart0","socpart1","socpart2")
se.lnch.all <- c(0.1359842,0.1318469,0.2704908)
names(se.lnch.all) <- c("socpart0","socpart1","socpart2")

lnch.age <- matrix(c(-0.8302066,-0.7410141,-1.918386 ,
                      0.1789821,-0.2590429,-0.9619553),
                   nrow=2,byrow=T,dimnames=list(c("age0","age1"),c("socpart0","socpart1","socpart2")))
se.lnch.age <- matrix(c(0.2103195,0.2016724,0.3328082,
                        0.1722267,0.1849641,0.3208147),
                      nrow=2,byrow=T,dimnames=list(c("age0","age1"),c("socpart0","socpart1","socpart2")))

lnch.sex <- matrix(c(-0.6542245,-0.5362902,-1.237453 ,
                     -0.0542673,-0.563398 ,-0.5816566),
                   nrow=2,byrow=T,dimnames=list(c("sex0","sex1"),c("socpart0","socpart1","socpart2")))
se.lnch.sex <- matrix(c(0.1796125,0.1829337,0.4535605,
                        0.2249564,0.1811729,0.3301233),
                      nrow=2,byrow=T,dimnames=list(c("sex0","sex1"),c("socpart0","socpart1","socpart2")))

lnch.owner <- matrix(c(-0.3880301,-0.5644763,-0.9528314,
                       -0.4840522,-0.3455065,-2.16895  ),
                     nrow=2,byrow=T,dimnames=list(c("owner0","owner1"),c("socpart0","socpart1","socpart2")))
se.lnch.owner <- matrix(c(0.1635307,0.1421271,0.2981138,
                          0.2859218,0.2823356,0.4323213),
                        nrow=2,byrow=T,dimnames=list(c("owner0","owner1"),c("socpart0","socpart1","socpart2")))

lnch.size <- matrix(c(-0.4609295,-0.5821285,-1.12553  ,
                      -0.5258389,-0.5076232,-0.5301442),
                    nrow=2,byrow=T,dimnames=list(c("size0","size1"),c("socpart0","socpart1","socpart2")))
se.lnch.size <- matrix(c(0.2184854,0.1982668,0.3763036,
                         0.2324401,0.1654763,0.507702 ),
                       nrow=2,byrow=T,dimnames=list(c("size0","size1"),c("socpart0","socpart1","socpart2")))

lnch.bar <- matrix(c(-0.8784455,-0.8410527,-1.823073,
                     -0.2838791,-0.3580084,-0.844191),
                   nrow=2,byrow=T,dimnames=list(c("bar0","bar1"),c("socpart0","socpart1","socpart2")))
se.lnch.bar <- matrix(c(0.2973433,0.2251123,0.3420723,
                        0.1584453,0.1523986,0.3066228),
                      nrow=2,byrow=T,dimnames=list(c("bar0","bar1"),c("socpart0","socpart1","socpart2")))

# Standardized differences in 5-year cumulative mortality risks and their SEs (delta methods) between levels of social engagement for each subgroup
cm.all <- 1 - exp(-exp(lnch.all))
var.cm.all <- (exp(-exp(lnch.all))*exp(lnch.all)*se.lnch.all)^2
rd.all <- cm.all[-1] - cm.all[1]
se.rd.all <- sqrt(var.cm.all[-1]+var.cm.all[1])

RD <- se.RD <- p.homo <- NULL
mat <- cbind(-1,diag(2))
for(covar in c("age","sex","owner","size","bar")){
   lnch <- get(paste("lnch.",covar,sep=""))
   se.lnch <- get(paste("se.lnch.",covar,sep=""))
   cm <- 1 - exp(-exp(lnch))
   var.cm <- (exp(-exp(lnch))*exp(lnch)*se.lnch)^2
   RD <- rbind(RD,cm[,-1]-cm[,1])
   se.RD <- rbind(se.RD,sqrt(var.cm[,-1]+var.cm[,1]))
   b <- cm[2,]-cm[1,]
   V <- diag(var.cm[2,]+var.cm[1,])
   p.homo <- c(p.homo,1-pchisq(t(mat%*%b)%*%solve(mat%*%V%*%t(mat))%*%(mat%*%b),2))  # Test for homogeneity of risk differences across subgroups
   }
names(p.homo) <- c("age","sex","owner","size","bar")

RD <- rbind(RD,overall=rd.all)
se.RD <- rbind(se.RD,overall=se.rd.all)

# Figure 2
var <- c("Age (years)","Sex","Facility ownership","Facility size (beds)","Functional dependency","Overall")
label <- c("65–84",expression("">="85"),
           "Women","Men",
           "Public/subsidized","Private",
           "< 300",expression("">="300"),
           "No","Mild/moderate")

est.RD <- matrix(format(100*RD,digits=0,nsmall=1,trim=T),nrow=11)
ci.RD <- matrix(paste(" (",format(100*(RD-qnorm(0.975)*se.RD),digits=1,nsmall=1,trim=T)," to ",
                           format(100*(RD+qnorm(0.975)*se.RD),digits=1,nsmall=1,trim=T),")",sep=""),nrow=11)
p.homo <- format(p.homo,digits=2)

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

par(mfrow=c(1,2),mar=c(1.6,7.5,3.1,0),oma=c(0,10,0,0.25),tcl=-0.35)
plot(RD[,1],c(20,19,16,15,12,11,8,7,4,3,0),type="n",xlim=c(-0.50,0.20),ylim=c(-0.50,20),xlab="",ylab="",axes=F)
   symbols(RD[-11,1],c(20,19,16,15,12,11,8,7,4,3,0)[-11],
           squares=1/se.RD[-11,1],inches=0.065*max(1/se.RD[-11,1])/max(1/se.RD[-11,]),add=T,lty=1,lwd=1,bg="black")
   segments((RD-qnorm(0.975)*se.RD)[-c(6,11),1],c(20,19,16,15,12,11,8,7,4,3,0)[-c(6,11)],
            (RD+qnorm(0.975)*se.RD)[-c(6,11),1],c(20,19,16,15,12,11,8,7,4,3,0)[-c(6,11)],lty=1,lwd=1.5)
   arrows((RD-qnorm(0.975)*se.RD)[6,1],c(20,19,16,15,12,11,8,7,4,3,0)[6],
          0.20,c(20,19,16,15,12,11,8,7,4,3,0)[6],code=2,length=0.035,lty=1,lwd=1.5)
   polygon(c((RD-qnorm(0.975)*se.RD)[11,1],RD[11,1],(RD+qnorm(0.975)*se.RD)[11,1],RD[11,1]),
           c(0,-0.50,0,0.50),border=F,density=-1,lty=1,lwd=1,col=1)
   segments(RD[11,1],0.50,RD[11,1],20,lty=2,lwd=1)
   abline(v=0,lty=1,lwd=1)
   axis(1,at=seq(-0.50,0.20,0.10),labels=seq(-50,20,10),cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   mtext(side=2,at=c(21,17,13,9,5,0),var,cex=0.85,font=1,las=1,adj=0,line=17)
   mtext(side=2,at=22,"Baseline subgroup",cex=0.85,font=1,las=1,adj=0,line=17)
   mtext(side=2,at=c(20,19,16,15,12,11,8,7,4,3)[-c(2,8)],label[-c(2,8)],cex=0.85,font=1,las=1,adj=0,line=16)
   mtext(side=2,at=c(20,19,16,15,12,11,8,7,4,3)[c(2,8)],label[c(2,8)],cex=0.85,font=1,las=1,adj=0,line=16.1)
   mtext(side=2,at=c(18,14,10,6,2)-0.05,expression(paste(italic(P)," for heterogeneity = ")),cex=0.85,font=1,las=1,adj=1,line=9.1)
   mtext(side=2,at=c(18,14,10,6,2),p.homo,cex=0.85,font=1,las=1,adj=0,line=9.1)
   mtext(side=2,at=c(20,19,16,15,12,11,8,7,4,3,0),est.RD[,1],cex=0.85,font=1,las=1,adj=1,line=5)
   mtext(side=2,at=c(20,19,16,15,12,11,8,7,4,3,0),ci.RD[,1],cex=0.85,font=1,las=1,adj=0,line=5)
   mtext(side=2,at=22,"Moderate vs. low/null social engagement at baseline",cex=0.85,font=1,las=1,adj=0.5,line=-2.6)
   mtext(side=2,at=23,"Standardized difference in 5-year cumulative all-cause mortality (%)",cex=0.85,font=1,las=1,adj=0.5,line=-12)
plot(RD[,2],c(20,19,16,15,12,11,8,7,4,3,0),type="n",xlim=c(-0.50,0.20),ylim=c(-0.50,20),xlab="",ylab="",axes=F)
   symbols(RD[-11,2],c(20,19,16,15,12,11,8,7,4,3,0)[-11],
           squares=1/se.RD[-11,2],inches=0.065*max(1/se.RD[-11,2])/max(1/se.RD[-11,]),add=T,lty=1,lwd=1,bg="black")
   segments((RD-qnorm(0.975)*se.RD)[-c(2,6,8,11),2],c(20,19,16,15,12,11,8,7,4,3,0)[-c(2,6,8,11)],
            (RD+qnorm(0.975)*se.RD)[-c(2,6,8,11),2],c(20,19,16,15,12,11,8,7,4,3,0)[-c(2,6,8,11)],lty=1,lwd=1.5)
   arrows(rep(-0.50,2),c(20,19,16,15,12,11,8,7,4,3,0)[c(2,6)],
          (RD+qnorm(0.975)*se.RD)[c(2,6),2],c(20,19,16,15,12,11,8,7,4,3,0)[c(2,6)],code=1,length=0.035,lty=1,lwd=1.5)
   arrows((RD-qnorm(0.975)*se.RD)[8,2],c(20,19,16,15,12,11,8,7,4,3,0)[8],
          0.20,c(20,19,16,15,12,11,8,7,4,3,0)[8],code=2,length=0.035,lty=1,lwd=1.5)
   polygon(c((RD-qnorm(0.975)*se.RD)[11,2],RD[11,2],(RD+qnorm(0.975)*se.RD)[11,2],RD[11,2]),
           c(0,-0.50,0,0.50),border=F,density=-1,lty=1,lwd=1,col=1)
   segments(RD[11,2],0.50,RD[11,2],20,lty=2,lwd=1)
   abline(v=0,lty=1,lwd=1)
   axis(1,at=seq(-0.50,0.20,0.10),labels=seq(-50,20,10),cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=1)
   mtext(side=2,at=c(20,19,16,15,12,11,8,7,4,3,0),est.RD[,2],cex=0.85,font=1,las=1,adj=1,line=5)
   mtext(side=2,at=c(20,19,16,15,12,11,8,7,4,3,0),ci.RD[,2],cex=0.85,font=1,las=1,adj=0,line=5)
   mtext(side=2,at=22,"High vs. low/null social engagement at baseline",cex=0.85,font=1,las=1,adj=0.5,line=-2.7)

dev.off()