setwd("C:/Users/kyra/Desktop/DBSB Survival/Part 1") library(prodlim) library(sas7bdat) library(survival) drop=read.sas7bdat("dropout.sas7bdat") head(drop) ?prodlim ?Hist drop$ind=ifelse(drop$state > 0,1,0) km=prodlim(Hist(week,ind)~drug,data=drop) summary(km,surv=FALSE) plot(km,type="cuminc") plot(km,type="cuminc",ylim=c(0,0.4),background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend="FALSE") axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) legend("topleft",legend=c("Drug","Placebo"),col=c("red","black"),lwd=2,bty="n") aj=prodlim(Hist(week,state)~drug,data=drop) summary(aj) par(mfrow=c(2,2)) plot(aj,cause=1,ylim=c(0,0.4),background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE) axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) plot(aj,cause=2,ylim=c(0,0.4),background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE) axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) plot(aj,cause=3,ylim=c(0,0.4),background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE) axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) drop$ind1=ifelse(drop$state==1,1,0) drop$ind2=ifelse(drop$state==2,1,0) drop$ind3=ifelse(drop$state==3,1,0) c1=prodlim(Hist(week,ind1)~drug,data=drop) summary(c1,surv=FALSE) c2=prodlim(Hist(week,ind2)~drug,data=drop) summary(c2,surv=FALSE) c3=prodlim(Hist(week,ind3)~drug,data=drop) summary(c3,surv=FALSE) ##plots by cause of drop-out par(mfrow=c(1,3)) plot(c1,type="cuminc",background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE,ylim=c(0,0.2), plot.main="CIF cause=1") axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) plot(aj,cause=1,add=T,lty=2,confint=FALSE,col=c("green","purple")) legend("topleft",legend=c("1-KM drug","1-KM Placebo","AJ drug","AJ placebo"), col=c("red","black","purple","green"),lty=c(1,1,2,2),bty="n",lwd=2) plot(c2,type="cuminc",background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE,ylim=c(0,0.2), plot.main="CIF cause=2") axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) plot(aj,cause=2,add=T,lty=2,confint=FALSE,col=c("green","purple")) legend("topleft",legend=c("1-KM drug","1-KM Placebo","AJ drug","AJ placebo"), col=c("red","black","purple","green"),lty=c(1,1,2,2),bty="n",lwd=2) plot(c3,type="cuminc",background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE,ylim=c(0,0.2), plot.main="CIF cause=3") axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) plot(aj,cause=3,add=T,lty=2,confint=FALSE,col=c("green","purple")) legend("topleft",legend=c("1-KM drug","1-KM Placebo","AJ drug","AJ placebo"), col=c("red","black","purple","green"),lty=c(1,1,2,2),bty="n",lwd=2) ##stacked plots drug=prodlim(Hist(week,state)~1,data=drop[drop$drug==1,]) placebo=prodlim(Hist(week,state)~1,data=drop[drop$drug==0,]) par(mfrow=c(1,2)) plot(drug,cause="stacked",background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE,ylim=c(0,0.4),xlim=c(0,52), plot.main="CIF Drug") axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0)) plot(placebo,cause="stacked",background=FALSE,axes=FALSE,atrisk=FALSE,confint=FALSE,legend=FALSE,ylim=c(0,0.4),xlim=c(0,52), plot.main="CIF placebo") axis(side=1,pos=c(0,0)) axis(side=2,pos=c(0,0))