library(survival)
library(MASS)

#ADDICTS <- read.table("D:/Courses/survival/Tablman monograph/C4088/Datasets/ADDICTS.txt",sep="\t",header=T)

attach(ADDICTS)

ADDICTS
summary(ADDICTS)

surv <- Surv(Days.survival,Status)

#	Q1:

par(mfcol=c(1,1))
plot(survfit(surv~Clinic),col=c(1,2))
legend(800,1,c("Clinic 1","Clinic 2"),lty=1,lwd=2,col=c(1,2))


#	Q2:

print(survfit(surv~Clinic),show.rmean=T)

#	Q3:

survdiff(surv~Clinic)


#	Q4, alef:

# first input the following R function for estimating the cumulative hazard:

hazard.km <- 
function(data){
## Author: Mara Tableman  Date: 20 November 2002 
## Purpose:  To compute the two types of empirical hazards,
##           and the cumulative hazards along their s.e.'s
## Arguments:  data is survfit object
time <- summary(data)$time
ni <- summary(data)$n.risk
di <- summary(data)$n.event
surv <- summary(data)$surv
stderr <- summary(data)$std.err
hitilde <- di/ni
tau <- diff(time, lag = 1)#length of interval to right of ti
tau[length(tau) + 1] <- NA
hihat <- hitilde/tau
Hhat <-  -log(surv)
Htilde <- cumsum(hitilde)
sqri <- di/(ni^2)
se.Hhat <- stderr/surv
se.Htilde <- sqrt(cumsum(sqri))
hazardtable <- round(data.frame(time, ni, di, hihat, hitilde, Hhat,se.Hhat, Htilde, se.Htilde),4)
} 


# next run the following

surv.1 <- Surv(Days.survival[Clinic==1],Status[Clinic==1])
surv.2 <- Surv(Days.survival[Clinic==2],Status[Clinic==2])

par(mfcol=c(1,2))
h1 <- hazard.km(survfit(surv.1))
h2 <- hazard.km(survfit(surv.2))

xlm <- range(c(h2[,1],h1[,1]))
ylm <- range(c(h2[,8],h1[,8]))

plot(h1[,1],h1[,8],pch=1,lty=1,xlab="time",ylab="log-cumulative hazard",
main="log-Htilde comparison Clinic = 1, 2",log="y",xlim=xlm,ylim=ylm)
lines(h1[,1],h1[,8],col=4)
lines(h2[,1],h2[,8],col=5)
points(h2[,1],h2[,8],pch=2)

plot(h1[,1],h1[,8],pch=1,lty=1,xlab="time",ylab="log-cumulative hazard",
main="log-Htilde comparison Clinic = 1, 2",log="xy",xlim=xlm,ylim=ylm)
lines(h1[,1],h1[,8],col=4)
lines(h2[,1],h2[,8],col=5)
points(h2[,1],h2[,8],pch=2)

#	Q4, bet & gimmel:

(fit.0 <- coxph(surv ~ Prison + Dose))
(fit.1 <- coxph(surv ~ Prison + Dose + Clinic ))



#	Q4, dalet:

par(mfcol=c(1,1))
plot(survfit(fit.1, data.frame(Prison = 0, Dose = 62, Clinic = 1),conf.type="none")) 
lines(survfit(fit.1, data.frame(Prison = 0, Dose = 62, Clinic = 2)),col=2) 

(a1 <- summary(survfit(fit.1, data.frame(Prison = 0, Dose = 62, Clinic = 1))) )
(a2 <- summary(survfit(fit.1, data.frame(Prison = 0, Dose = 62, Clinic = 2))) )

approx(a1$time,a1$surv,xout=500)
approx(a2$time,a2$surv,xout=500)


#	Q4, heh:


summary(coxph(surv ~ Prison + Dose,subset=(Clinic==1)))
summary(coxph(surv ~ Prison + Dose,subset=(Clinic==2)))

dose.1 <- ifelse(Clinic == 1, Dose, 0)
dose.2 <- ifelse(Clinic == 2, Dose, 0)

coxph(surv ~ Prison + Dose * Clinic )
coxph(surv ~ Dose + Prison * Clinic )



#	Q5, alef:

# fit QQ for the 3 parameteric models

surv.fit <- survfit(surv)
surv.vec <-  summary(surv.fit)$surv
time.vec <- summary(surv.fit)$time
p.vec <- 1 - (surv.vec + c(1,surv.vec[-length(surv.vec)]))/2

y.vec <- log(time.vec)

par(mfcol=c(1,3))
x  <- log(qweibull(p.vec,shape=1,scale=1))
plot(x,y.vec,main="Weibull") ; abline(lm(y.vec~x))

x  <- qlogis(p.vec,location=0,scale=1)
plot(x,y.vec,main="Loglog") ; abline(lm(y.vec~x))

x  <- qnorm(p.vec,mean=0,sd=1)
plot(x,y.vec,main="Lognorm") ; abline(lm(y.vec~x))


#	Q5, bet:

lglg.fit <- survreg(surv ~ Dose + Prison + Clinic,dist="loglogistic")
summary(lglg.fit)
lglg.res <- resid(lglg.fit,type="deviance")

par(mfcol=c(1,2))
plot(1:238,lglg.res,ylab="Deviance residuals",xlab="Index")
lines(lowess(1:238,lglg.res),col=5) 
abline(0,0,col=2) 

plot(log(predict(lglg.fit)),lglg.res,ylab="Deviance residuals",xlab="Pred.")
lines(lowess(log(predict(lglg.fit)),lglg.res),col=5) 
abline(0,0,col=2) 


#  What could have been discovered ??

lglg.fit <- survreg(surv ~ Prison + Clinic,dist="loglogistic")
#lglg.fit <- survreg(surv ~ Prison + Clinic + Dose,dist="loglogistic")
lglg.res <- resid(lglg.fit,type="deviance")

par(mfcol=c(1,1))
plot(Dose,lglg.res,ylab="Deviance residuals",xlab="Dose")
lines(lowess(Dose,lglg.res),col=5) 
abline(0,0,col=2) 



#	Q5, gimmel:


lglg.fit <- survreg(surv ~ Prison + Dose + Clinic,dist="loglogistic")

predict(lglg.fit,newdata=data.frame(Prison = 0, Dose = 62, Clinic = 1),se.fit=T,type="quantile",p=c(0.25,0.75,1-0.49))
predict(lglg.fit,newdata=data.frame(Prison = 0, Dose = 62, Clinic = 2),se.fit=T,type="quantile",p=c(0.25,0.75,1-0.77))

par(mfcol=c(1,1))
plot(survfit(fit.1, data.frame(Prison = 0, Dose = 62, Clinic = 1),conf.type="none"),lty=2) 
lines(survfit(fit.1, data.frame(Prison = 0, Dose = 62, Clinic = 2)),col=2,lty=2) 

logis.mu1 <- predict(lglg.fit,newdata=data.frame(Prison = 0, Dose = 62, Clinic = 1),type="lp")
logis.mu2 <- predict(lglg.fit,newdata=data.frame(Prison = 0, Dose = 62, Clinic = 2),type="lp")
logis.sig <- lglg.fit$scale

x.ser <- seq(0,900,length=400)
lines(x.ser,1-plogis(log(x.ser),location=logis.mu1,scale=logis.sig),lty=1)
lines(x.ser,1-plogis(log(x.ser),location=logis.mu2,scale=logis.sig),lty=1,col=2)

abline(h=0.5,lty=1,col=3)
abline(h=0.25,lty=1,col=3)

lglg.surv1 <- 1-plogis(log(x.ser),location=logis.mu1,scale=logis.sig)
lglg.surv2  <- 1-plogis(log(x.ser),location=logis.mu2,scale=logis.sig)

approx(lglg.surv1,x.ser,xout=0.5)$y
approx(lglg.surv1,x.ser,xout=0.25)$y
approx(lglg.surv2,x.ser,xout=0.5)$y
approx(lglg.surv2,x.ser,xout=0.25)$y

exp(qlogis(0.5,location=logis.mu1,scale=logis.sig))
exp(qlogis(0.75,location=logis.mu1,scale=logis.sig))
exp(qlogis(0.5,location=logis.mu2,scale=logis.sig))
exp(qlogis(0.75,location=logis.mu2,scale=logis.sig))


#	Q5, dalet:

#  Odds ratio at 500 and 300

abline(v=500,lty=1,col=4)
abline(v=300,lty=1,col=4)

(s500.1 <- approx(x.ser,lglg.surv1,xout=500)$y)
(s500.2 <- approx(x.ser,lglg.surv2,xout=500)$y)

(s300.1 <- approx(x.ser,lglg.surv1,xout=300)$y)
(s300.2 <- approx(x.ser,lglg.surv2,xout=300)$y)

(odds.1 <- s500.1 / (1-s500.1))
(odds.2 <- s500.2 / (1-s500.2))

odds.2 / odds.1

(s300.2 / (1-s300.2)) / (s300.1 / (1-s300.1))

alpha <- 1/lglg.fit$scale
exp(lglg.fit$coef[4]*alpha)  # Odds Ratio  =  exp( - alpha * (x.14 - x.24) * b.4)