# ----------------------------------------- # # TYPE OF FACILITY AND MORTALITY - R SCRIPT # # ----------------------------------------- # # READ TEXT FILES # 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) names(survival)[c(3,4)] <- c("dth","time") survival <- survival[order(survival$time),] # FIGURE 1: STANDARDIZED NONPARAMETRIC AND SMOOTH CUMULATIVE MORTALITY CURVES BY TYPE OF FACILITY # Smooth cumulative mortality curves and their 95% CIs for each type of facility under the two standardizations survival$cm.spl.std1 <- 1 - exp(-exp(survival$lnch.spl.std1)) survival$ll.cm.spl.std1 <- 1 - exp(-exp(survival$lnch.spl.std1 - qnorm(0.975)*survival$se.lnch.spl.std1)) survival$ul.cm.spl.std1 <- 1 - exp(-exp(survival$lnch.spl.std1 + qnorm(0.975)*survival$se.lnch.spl.std1)) survival$cm.spl.std2 <- 1 - exp(-exp(survival$lnch.spl.std2)) survival$ll.cm.spl.std2 <- 1 - exp(-exp(survival$lnch.spl.std2 - qnorm(0.975)*survival$se.lnch.spl.std2)) survival$ul.cm.spl.std2 <- 1 - exp(-exp(survival$lnch.spl.std2 + qnorm(0.975)*survival$se.lnch.spl.std2)) # 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.75,3.5,0,0),tcl=-0.35) plot(survival$time,1-survival$s.km.std2,type="n",xlim=c(0,15),ylim=c(0,1),xlab="",ylab="",axes=F) lines(survival$time[survival$time<=15&survival$factype==0],1-survival$s.km.std2[survival$time<=15&survival$factype==0], type="s",lty=1,lwd=0.75,col="red3") lines(survival$time[survival$time<=15&survival$factype==1],1-survival$s.km.std2[survival$time<=15&survival$factype==1], type="s",lty=1,lwd=0.75,col="darkorange2") lines(survival$time[survival$time<=15&survival$factype==2],1-survival$s.km.std2[survival$time<=15&survival$factype==2], type="s",lty=1,lwd=0.75,col="green2") lines(survival$time[survival$time<=15&survival$factype==3],1-survival$s.km.std2[survival$time<=15&survival$factype==3], type="s",lty=1,lwd=0.75,col="deepskyblue2") lines(survival$time[survival$time<=15&survival$factype==0],survival$cm.spl.std2[survival$time<=15&survival$factype==0], type="l",lty=1,lwd=1.5,col="red3") lines(survival$time[survival$time<=15&survival$factype==1],survival$cm.spl.std2[survival$time<=15&survival$factype==1], type="l",lty=1,lwd=1.5,col="darkorange2") lines(survival$time[survival$time<=15&survival$factype==2],survival$cm.spl.std2[survival$time<=15&survival$factype==2], type="l",lty=1,lwd=1.5,col="green2") lines(survival$time[survival$time<=15&survival$factype==3],survival$cm.spl.std2[survival$time<=15&survival$factype==3], type="l",lty=1,lwd=1.5,col="deepskyblue2") axis(1,at=seq(0,15,5),labels=seq(0,15,5),cex.axis=0.85,font.axis=1,mgp=c(3,0.35,0),lty=1,lwd=0.75) axis(2,at=seq(0,1,0.25),labels=seq(0,1,0.25)*100,cex.axis=0.85,font.axis=1,las=1,adj=1,mgp=c(3,0.65,0),lty=1,lwd=0.75) mtext(side=1,"Follow-up time (years)",cex=1,font=1,line=1.5) mtext(side=2,"Cumulative all-cause mortality (%)",cex=1,font=1,las=3,line=2.25) legend("topleft",inset=0.01,bty="n",y.intersp=1.1, legend=c("Facility","Large-sized public","Medium-sized public/subsidized","Medium-sized private","Small-sized private"), lty=1,lwd=1.5,col=c("white","red3","darkorange2","green2","deepskyblue2"),cex=0.85) dev.off() # TABLE 2: STANDARDIZED DIFFERENCES IN CUMULATIVE MORTALITY AT 2, 5, AND 10 YEARS OF FOLLOW-UP BY TYPE OF FACILITY # Log cumulative hazards and their SEs at 2, 5, and 10 years for each facility type (spline interpolation from values at observed times) lnch.std1 <- matrix(nrow=3,ncol=4,dimnames=list(c("t2","t5","t10"),c("factype0","factype1","factype2","factype3"))) se.lnch.std1 <- lnch.std1 for(i in 0:3){ lnch.std1[,i+1] <- spline(log(survival$time[survival$factype==i]),survival$lnch.spl.std1[survival$factype==i],xout=log(c(2,5,10)),method="fmm")$y se.lnch.std1[,i+1] <- spline(log(survival$time[survival$factype==i]),survival$se.lnch.spl.std1[survival$factype==i],xout=log(c(2,5,10)),method="fmm")$y } lnch.std2 <- matrix(nrow=3,ncol=4,dimnames=list(c("t2","t5","t10"),c("factype0","factype1","factype2","factype3"))) se.lnch.std2 <- lnch.std2 for(i in 0:3){ lnch.std2[,i+1] <- spline(log(survival$time[survival$factype==i]),survival$lnch.spl.std2[survival$factype==i],xout=log(c(2,5,10)),method="fmm")$y se.lnch.std2[,i+1] <- spline(log(survival$time[survival$factype==i]),survival$se.lnch.spl.std2[survival$factype==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=1,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) # TABLE 3: STANDARDIZED DIFFERENCES IN FOLLOW-UP TIMES TO 35%, 50%, AND 65% CUMULATIVE MORTALITY RISKS BY TYPE OF FACILITY # Follow-up times to 35%, 50%, and 65% cumulative mortality risks and their SEs for each type of facility t.std1 <- matrix(c(2.608908,3.159894,3.947797,3.475313, 4.285734,4.912112,7.009877,5.455101, 6.745677,7.412544,11.85726,8.101406), nrow=3,byrow=T,dimnames=list(c("p35","p50","p65"),c("factype0","factype1","factype2","factype3"))) se.t.std1 <- matrix(c(0.332436 ,0.4038666,0.9076167,0.5870994, 0.4842063,0.5889183,1.459509 ,0.7832254, 0.6674207,0.8449553,2.47597 ,1.028859 ), nrow=3,byrow=T,dimnames=list(c("p35","p50","p65"),c("factype0","factype1","factype2","factype3"))) t.std2 <- matrix(c(2.34487 ,3.161648,4.482917,3.823437, 3.978936,4.90547 ,7.601904,5.900443, 6.610855,7.466687,12.16812,8.596556), nrow=3,byrow=T,dimnames=list(c("p35","p50","p65"),c("factype0","factype1","factype2","factype3"))) se.t.std2 <- matrix(c(0.3521731,0.4284832,1.0392 ,0.5852339, 0.5747114,0.6432603,1.491249,0.7689705, 0.9624053,0.946332 ,2.411367,1.033384 ), nrow=3,byrow=T,dimnames=list(c("p35","p50","p65"),c("factype0","factype1","factype2","factype3"))) # Confidence intervals for the differences in follow-up times matrix(paste(format(t.std1[,-1] - t.std1[,1],digits=1,trim=T)," (", format(t.std1[,-1] - t.std1[,1] - qnorm(0.975)*sqrt(se.t.std1[,-1]^2+se.t.std1[,1]^2),digits=0,nsmall=1,trim=T)," to ", format(t.std1[,-1] - t.std1[,1] + qnorm(0.975)*sqrt(se.t.std1[,-1]^2+se.t.std1[,1]^2),digits=0,nsmall=1,trim=T),")",sep=""),nrow=3) matrix(paste(format(t.std2[,-1] - t.std2[,1],digits=1,trim=T)," (", format(t.std2[,-1] - t.std2[,1] - qnorm(0.975)*sqrt(se.t.std2[,-1]^2+se.t.std2[,1]^2),digits=0,nsmall=1,trim=T)," to ", format(t.std2[,-1] - t.std2[,1] + qnorm(0.975)*sqrt(se.t.std2[,-1]^2+se.t.std2[,1]^2),digits=1,nsmall=1,trim=T),")",sep=""),nrow=3) # FIGURE 2: STANDARDIZED DIFFERENCES IN 5-YEAR CUMULATIVE MORTALITY AMONG TYPES OF FACILITY BY SUBGROUP # Log cumulative hazards and their SEs at 5 years of follow-up for each type of facility and resident subgroup lnch.age <- matrix(c(-0.4871136,-0.72503 ,-0.8219972,-0.6253072, 0.1785898, 0.0104031,-0.6689032,-0.2955748), nrow=2,byrow=T,dimnames=list(c("age0","age1"),c("factype0","factype1","factype2","factype3"))) se.lnch.age <- matrix(c(0.1942836,0.1699694,0.3435879,0.2020701, 0.1238823,0.1964779,0.2939996,0.2334796), nrow=2,byrow=T,dimnames=list(c("age0","age1"),c("factype0","factype1","factype2","factype3"))) lnch.sex <- matrix(c(-0.2129048,-0.4777274,-1.082866 ,-0.6158715, -0.1449048,-0.2017469,-0.1765523,-0.5288587), nrow=2,byrow=T,dimnames=list(c("sex0","sex1"),c("factype0","factype1","factype2","factype3"))) se.lnch.sex <- matrix(c(0.1615192,0.1911236,0.2806308,0.188908 , 0.1452927,0.1569496,0.2316074,0.236197), nrow=2,byrow=T,dimnames=list(c("sex0","sex1"),c("factype0","factype1","factype2","factype3"))) lnch.educ <- matrix(c(-0.1256174,-0.2943419,-0.8373717,-1.239264 , -0.1756018,-0.4494317,-0.6506289,-0.3894342), nrow=2,byrow=T,dimnames=list(c("educ0","educ1"),c("factype0","factype1","factype2","factype3"))) se.lnch.educ <- matrix(c(0.1296716,0.1978863,0.4214151,0.7312996, 0.1702239,0.1894238,0.1841399,0.141281 ), nrow=2,byrow=T,dimnames=list(c("educ0","educ1"),c("factype0","factype1","factype2","factype3"))) lnch.stay <- matrix(c( 0.0835501,-0.2024604,-0.4896653,-0.4206716, -0.2825336,-0.4256342,-0.9958541,-0.7373466), nrow=2,byrow=T,dimnames=list(c("stay0","stay1"),c("factype0","factype1","factype2","factype3"))) se.lnch.stay <- matrix(c(0.1925949,0.1536246,0.2402243,0.1774187, 0.1055043,0.1972737,0.3740195,0.3603284), nrow=2,byrow=T,dimnames=list(c("stay0","stay1"),c("factype0","factype1","factype2","factype3"))) lnch.dem <- matrix(c(-0.3546984,-0.4638559,-0.6630431,-0.5471544, 0.4381047, 0.1004327,-0.7840685,-0.6258733), nrow=2,byrow=T,dimnames=list(c("dem0","dem1"),c("factype0","factype1","factype2","factype3"))) se.lnch.dem <- matrix(c(0.1590958,0.1480463,0.2283596,0.1837014, 0.1579198,0.1866224,0.4472981,0.3842935), nrow=2,byrow=T,dimnames=list(c("dem0","dem1"),c("factype0","factype1","factype2","factype3"))) lnch.ndis <- matrix(c(-0.5163615,-0.3418415,-0.8063789,-0.6193898, 0.0982925,-0.2877999,-0.7052239,-0.2931403), nrow=2,byrow=T,dimnames=list(c("ndis0","ndis1"),c("factype0","factype1","factype2","factype3"))) se.lnch.ndis <- matrix(c(0.2581225,0.1594591,0.3270186,0.2204978, 0.1242116,0.193336 ,0.2560003,0.2093569), nrow=2,byrow=T,dimnames=list(c("ndis0","ndis1"),c("factype0","factype1","factype2","factype3"))) lnch.bar <- matrix(c(-0.7006372,-0.6735002,-1.041816 ,-1.453901, -0.158837 ,-0.0283158,-0.4097493,-0.074176), nrow=2,byrow=T,dimnames=list(c("bar0","bar1"),c("factype0","factype1","factype2","factype3"))) se.lnch.bar <- matrix(c(0.1481214,0.2145434,0.3511056,0.3812504, 0.38638 ,0.1442693,0.2629978,0.1522465), nrow=2,byrow=T,dimnames=list(c("bar0","bar1"),c("factype0","factype1","factype2","factype3"))) # Standardized differences in 5-year cumulative mortality risks and their SEs (delta methods) between types of facility for each resident subgroup RD <- se.RD <- p.homo <- NULL mat <- cbind(-1,diag(3)) for(covar in c("age","sex","educ","stay","dem","ndis","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),3)) # Test for homogeneity of risk differences across resident subgroups } names(p.homo) <- c("age","sex","educ","stay","dem","ndis","bar") # Figure 2 var <- c("Age (years)","Sex","Educational level","Length of stay (years)","Dementia","No. of chronic conditions","Functional dependency","Overall") label <- c("65–84",expression("">="85"), "Women","Men", "Less than primary","Primary or more", "<3",expression("">="3"), "No","Yes", "0–2",expression("">="3"), "No/mild","Moderate/severe/total") RD <- rbind(RD,all=c(-0.06107804,-0.1892404,-0.12900663)) # Overall risk differences se.RD <- rbind(se.RD,all=c(0.06336190,0.07380360,0.06525554)) # Overall SEs ci.RD <- matrix(paste(format(100*RD,digits=1,trim=F)," (", format(100*(RD-qnorm(0.975)*se.RD),digits=1,nsmall=1,trim=T)," to ", format(100*(RD+qnorm(0.975)*se.RD),digits=0,nsmall=1,trim=T),")",sep=""),nrow=15) p.homo <- format(p.homo,digits=1) tiff(filename="c:/.../Figure 2.tiff", width=17.75,height=10,units="cm",pointsize=6,res=300) par(mfrow=c(1,3),mar=c(2.25,11,5.25,0),oma=c(0,14.5,0,0.25),tcl=-0.45) plot(RD[,1],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),type="n",xlim=c(-0.50,0.20),ylim=c(-0.50,28),xlab="",ylab="",axes=F) symbols(RD[-15,1],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-15], squares=1/se.RD[-15,1],inches=0.065*max(1/se.RD[-15,1])/max(1/se.RD[-15,]),add=T,lty=1,lwd=1,bg="black") segments((RD-qnorm(0.975)*se.RD)[-c(11,14,15),1],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(11,14,15)], (RD+qnorm(0.975)*se.RD)[-c(11,14,15),1],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(11,14,15)],lty=1,lwd=1.5) arrows((RD-qnorm(0.975)*se.RD)[c(11,14),1],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[c(11,14)], rep(0.20,2),c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[c(11,14)],code=2,length=0.035,lty=1,lwd=1.5) polygon(c((RD-qnorm(0.975)*se.RD)[15,1],RD[15,1],(RD+qnorm(0.975)*se.RD)[15,1],RD[15,1]), c(0,-0.50,0,0.50),border=F,density=-1,col=1) segments(RD[15,1],0.50,RD[15,1],28,lty=2,lwd=0.75) abline(v=0,lty=1,lwd=0.75) axis(1,at=seq(-0.50,0.20,0.10),labels=seq(-50,20,10),cex.axis=0.85/0.66,font.axis=1,mgp=c(3,0.85,0),lty=1,lwd=0.75) mtext(side=2,at=c(29,25,21,17,13,9,5,0),var,cex=0.85,font=1,las=1,adj=0,line=25) mtext(side=2,at=30,"Resident subgroup",cex=0.85,font=1,las=1,adj=0,line=25) mtext(side=2,at=c(28,27,24,23,20,19,16,15,12,11,8,7,4,3)[-c(2,8,12)],label[-c(2,8,12)],cex=0.85,font=1,las=1,adj=0,line=24) mtext(side=2,at=c(28,27,24,23,20,19,16,15,12,11,8,7,4,3)[c(2,8,12)],label[c(2,8,12)],cex=0.85,font=1,las=1,adj=0,line=24.2) mtext(side=2,at=c(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=13.6) mtext(side=2,at=c(26,22,18,14,10,6,2),p.homo,cex=0.85,font=1,las=1,adj=0,line=13.6) mtext(side=2,at=c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),ci.RD[,1],cex=0.85,font=1,las=1,adj=0,line=10) mtext(side=2,at=30.75,"Medium-sized public/subsidized",cex=0.85,font=1,las=1,adj=0.5,line=-4.75) mtext(side=2,at=30,"versus large-sized public facilities",cex=0.85,font=1,las=1,adj=0.5,line=-4.75) mtext(side=2,at=31.75,"Standardized differences in 5-year cumulative all-cause mortality (%)",cex=0.85,font=1,las=1,adj=0.5,line=-35.25) plot(RD[,2],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),type="n",xlim=c(-0.50,0.20),ylim=c(-0.50,28),xlab="",ylab="",axes=F) symbols(RD[-15,2],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-15], squares=1/se.RD[-15,2],inches=0.065*max(1/se.RD[-15,2])/max(1/se.RD[-15,]),add=T,lty=1,lwd=1,bg="black") segments((RD-qnorm(0.975)*se.RD)[-c(10,14,15),2],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(10,14,15)], (RD+qnorm(0.975)*se.RD)[-c(10,14,15),2],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(10,14,15)],lty=1,lwd=1.5) arrows(-0.50,c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[10], (RD+qnorm(0.975)*se.RD)[10,2],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[10],code=1,length=0.035,lty=1,lwd=1.5) arrows((RD-qnorm(0.975)*se.RD)[14,2],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[14], 0.20,c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[14],code=2,length=0.035,lty=1,lwd=1.5) polygon(c((RD-qnorm(0.975)*se.RD)[15,2],RD[15,2],(RD+qnorm(0.975)*se.RD)[15,2],RD[15,2]), c(0,-0.50,0,0.50),border=F,density=-1,col=1) segments(RD[15,2],0.50,RD[15,2],28,lty=2,lwd=0.75) abline(v=0,lty=1,lwd=0.75) axis(1,at=seq(-0.50,0.20,0.10),labels=seq(-50,20,10),cex.axis=0.85/0.66,font.axis=1,mgp=c(3,0.85,0),lty=1,lwd=0.75) mtext(side=2,at=c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),ci.RD[,2],cex=0.85,font=1,las=1,adj=0,line=10) mtext(side=2,at=30.75,"Medium-sized private",cex=0.85,font=1,las=1,adj=0.5,line=-4.75) mtext(side=2,at=30,"versus large-sized public facilities",cex=0.85,font=1,las=1,adj=0.5,line=-4.75) plot(RD[,3],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),type="n",xlim=c(-0.50,0.20),ylim=c(-0.50,28),xlab="",ylab="",axes=F) symbols(RD[-15,3],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-15], squares=1/se.RD[-15,3],inches=0.065*max(1/se.RD[-15,3])/max(1/se.RD[-15,]),add=T,lty=1,lwd=1,bg="black") segments((RD-qnorm(0.975)*se.RD)[-c(5,10,14,15),3],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(5,10,14,15)], (RD+qnorm(0.975)*se.RD)[-c(5,10,14,15),3],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[-c(5,10,14,15)],lty=1,lwd=1.5) arrows(rep(-0.50,2),c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[c(5,10)], (RD+qnorm(0.975)*se.RD)[c(5,10),3],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[c(5,10)],code=1,length=0.035,lty=1,lwd=1.5) arrows((RD-qnorm(0.975)*se.RD)[14,3],c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[14], 0.20,c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0)[14],code=2,length=0.035,lty=1,lwd=1.5) polygon(c((RD-qnorm(0.975)*se.RD)[15,3],RD[15,3],(RD+qnorm(0.975)*se.RD)[15,3],RD[15,3]), c(0,-0.50,0,0.50),border=F,density=-1,col=1) segments(RD[15,3],0.50,RD[15,3],28,lty=2,lwd=0.75) abline(v=0,lty=1,lwd=0.75) axis(1,at=seq(-0.50,0.20,0.10),labels=seq(-50,20,10),cex.axis=0.85/0.66,font.axis=1,mgp=c(3,0.85,0),lty=1,lwd=0.75) mtext(side=2,at=c(28,27,24,23,20,19,16,15,12,11,8,7,4,3,0),ci.RD[,3],cex=0.85,font=1,las=1,adj=0,line=10) mtext(side=2,at=30.75,"Small-sized private",cex=0.85,font=1,las=1,adj=0.5,line=-4.75) mtext(side=2,at=30,"versus large-sized public facilities",cex=0.85,font=1,las=1,adj=0.5,line=-4.75) dev.off() # SECONDARY ANALYSES CLASSIFYING FACILITY ACCORDING TO OWNERSHIP AND PROFIT STATUS # Log cumulative hazards and their SEs at 5 years of follow-up for each type of facility lnch <- c(-0.2178283,-0.2555195,-0.497803,-0.6657686) names(lnch) <- c("factypeb0","factypeb1","factypeb2","factypeb3") se.lnch <- c(0.0840613,0.2563074,0.1714064,0.2463056) names(se.lnch) <- c("factypeb0","factypeb1","factypeb2","factypeb3") # Standardized differences in 5-year cumulative mortality risks by type of facility cm <- 1 - exp(-exp(lnch)) var.cm <- (exp(-exp(lnch))*exp(lnch)*se.lnch)^2 paste(format(100*(cm[-1] - cm[1]),digits=1,nsmall=1,trim=T)," (", format(100*(cm[-1] - cm[1] - qnorm(0.975)*sqrt(var.cm[-1]+var.cm[1])),digits=1,nsmall=1,trim=T)," to ", format(100*(cm[-1] - cm[1] + qnorm(0.975)*sqrt(var.cm[-1]+var.cm[1])),digits=1,nsmall=1,trim=T),")",sep="") # Median survival times and their SEs for each type of facility t <- c(4.269333,4.480515,5.6107,7.296227) names(t) <- c("factypeb0","factypeb1","factypeb2","factypeb3") se.t <- c(0.3957669,1.090405,0.7904961,2.37241) names(se.t) <- c("factypeb0","factypeb1","factypeb2","factypeb3") # Standardized differences in median survival times by type of facility paste(format(t[-1] - t[1],digits=1,trim=T)," (", format(t[-1] - t[1] - qnorm(0.975)*sqrt(se.t[-1]^2+se.t[1]^2),digits=0,nsmall=1,trim=T)," to ", format(t[-1] - t[1] + qnorm(0.975)*sqrt(se.t[-1]^2+se.t[1]^2),digits=1,nsmall=1,trim=T),")",sep="")