library(survival)
attach(aml)


surv.fit <- survfit(Surv(time[x=="Maintained"],status[x=="Maintained"]),type="kaplan-meier")


aml.surv  <-  Surv(aml$time,aml$status)
aml1.surv <-  Surv(aml$time[aml$x=="Maintained"],   aml$status[aml$x=="Maintained"])
aml0.surv <-  Surv(aml$time[aml$x=="Nonmaintained"],aml$status[aml$x=="Nonmaintained"])

group <- ifelse(aml$x=="Nonmaintained",rep(0,23),rep(1,23))


#------- 2 sample exponential fit --------

expo.fit <- survreg(aml.surv~group,dist="expo")
summary(expo.fit)
expo.fit$coeff
expo.fit$var

sum(aml$status == 1 & group==1)
sum(aml$status == 1 & group==0)
sum(aml$time[group==1])
sum(aml$time[group==0])

summary(survreg(aml.surv~1,dist="expo"))


#------- 2 sample K-M  loglogsitic fit --------


loglog.fit1 <- survreg(aml1.surv~1,dist="loglogistic")
loglog.fit0 <- survreg(aml0.surv~1,dist="loglogistic")
loglog.fit  <- survreg(aml.surv~group,dist="loglogistic")

summary(loglog.fit)
summary(loglog.fit1)
summary(loglog.fit0)


plot(survfit(aml1.surv~1),conf.int=F,main="Loglogistic fit")
lines(survfit(aml0.surv~1),conf.int=F)

x.ser <- seq(0,161,length=30)
lines(x.ser,1-plogis(log(x.ser),location=3.515,scale=0.542),lty=3)
lines(x.ser,1-plogis(log(x.ser),location=2.9,scale=0.497),lty=3)

lines(x.ser,1-plogis(log(x.ser),location=2.899+0.604,scale=0.513),lty=3,col=2)
lines(x.ser,1-plogis(log(x.ser),location=2.899,scale=0.513),lty=3,col=2)

legend(80,.9,c("K-M (AML 0/1)","Loglog (AML 0/1)","Loglog (AML 0 & 1)"),lty=c(1,3,3),col=c(1,1,2))



# ---- 2 sample loglogistic medians ----
 
med1hat <- predict(loglog.fit1,type="uquantile",p=0.5,se.fit=T)
med0hat <- predict(loglog.fit0,type="uquantile",p=0.5,se.fit=T)
medhat  <- predict(loglog.fit,type="uquantile",p=0.5,se.fit=T)

CI.med1 <- exp(med1hat$fit[1]+c(0,-1.96,1.96)*med1hat$se.fit[1])
CI.med0 <- exp(med0hat$fit[1]+c(0,-1.96,1.96)*med0hat$se.fit[1])

names(CI.med1)<-c("median","LCL","UCL") ; names(CI.med0)<-c("median","LCL","UCL")

CI.med1 ; CI.med0

exp(medhat$fit[1]+c(0,-1.96,1.96)*medhat$se.fit[1])
exp(medhat$fit[13]+c(0,-1.96,1.96)*medhat$se.fit[13])



#----  2 Sample  loglogistic and Weibull regression model -----


summary(survReg(Surv(aml$weeks,aml$status)~aml$group,dist="loglogistic"))
summary(survReg(Surv(aml$weeks,aml$status)~aml$group,dist="weibull"))


# -----  4 nested models ----
 
weib.fit0 <- survReg(aml.surv~1,dist="weibull")
weib.fit1 <- survReg(aml.surv~aml$group,dist="weibull")
weib.fit20 <- survReg(aml1.surv~1,dist="weibull")
weib.fit21 <- survReg(aml2.surv~1,dist="weibull")
 summary(weib.fit0)
summary(weib.fit1)
summary(weib.fit20)
summary(weib.fit21)

anova(weib.fit0,weib.fit1,test="Chisq")

weib.fit20$loglik
 loglik3 <- weib.fit20$loglik[2] + weib.fit21$loglik[2]

loglik3
lrt23 <- -2 * (weib.fit1$loglik[2] -loglik3)

lrt23


1-pchisq(lrt23,1)


#----   Motorette example ---

library(MASS)

motorette <- motors[11:40,]
x      <- 1000 / (273.2+motorette$temp)

motorette.surv <- Surv(motorette$time,motorette$cens)
weib.fit <- survreg(motorette.surv~x,dist="weibull")

summary(weib.fit)
weib.fit$var
predict(weib.fit,newdata=data.frame(x=2.3),se.fit=T,type="quantile",p=c(.05,.50,.95))

mu <- sum(weib.fit$coeff[1] + weib.fit$coeff[2]*2.3)
sig  <- weib.fit$scale
exp(mu) * ( -log(1-.5))^sig




weib.fit <- survreg(motorette.surv~x,dist="weibull")
lglg.fit <- survreg(motorette.surv~x,dist="loglogistic")
lgnm.fit <- survreg(motorette.surv~x,dist="lognorm")

#-------------  parametric regression in log time  ----------------------------------------

par(mfcol=c(1,2))
ind <- motorette[,3] == 1
plot(x[ind],log(motorette[ind,2]),xlab="x",ylab="log survival time (hours)",ylim=range(log(motorette[,2])),
main="Log-time model")
points(x[!ind],log(motorette[!ind,2]),pch=3)

x.ser <- seq(from=min(x),to=max(x),length=10)
weib.pred <- predict(weib.fit,newdata=data.frame(x=x.ser),se.fit=F,type="uquantile",p=c(.25,.50,.75))
lglg.pred <- predict(lglg.fit,newdata=data.frame(x=x.ser),se.fit=F,type="uquantile",p=c(.25,.50,.75))
lgnm.pred <- predict(lgnm.fit,newdata=data.frame(x=x.ser),se.fit=F,type="uquantile",p=c(.25,.50,.75))
lines(x.ser,weib.pred[,1],lty=2) 
lines(x.ser,weib.pred[,2],lty=1)  
lines(x.ser,weib.pred[,3],lty=2)

lines(x.ser,lglg.pred[,1],lty=2) 
lines(x.ser,lglg.pred[,2],lty=1) 
lines(x.ser,lglg.pred[,3],lty=2)

lines(x.ser,lgnm.pred[,1],lty=2) 
lines(x.ser,lgnm.pred[,2],lty=1)  
lines(x.ser,lgnm.pred[,3],lty=2)



#-------------  parametric regression in time  ----------------------------------------

plot(x[ind],motorette[ind,2],xlab="x",ylab="survival time (hours)",ylim=range(motorette[2]),
main="Time model")
points(x[!ind],motorette[!ind,2],pch=3)

weib.pred <- predict(weib.fit,newdata=data.frame(x=x.ser),se.fit=F,type="quantile",p=c(.25,.50,.75))
lglg.pred <- predict(lglg.fit,newdata=data.frame(x=x.ser),se.fit=F,type="quantile",p=c(.25,.50,.75))
lgnm.pred <- predict(lgnm.fit,newdata=data.frame(x=x.ser),se.fit=F,type="quantile",p=c(.25,.50,.75))
lines(x.ser,weib.pred[,1],lty=2) ;
lines(x.ser,weib.pred[,2],lty=1) ; 
lines(x.ser,weib.pred[,3],lty=2)

lines(x.ser,lglg.pred[,1],lty=2) ;
lines(x.ser,lglg.pred[,2],lty=1) ; 
lines(x.ser,lglg.pred[,3],lty=2)

lines(x.ser,lgnm.pred[,1],lty=2) ;
 lines(x.ser,lgnm.pred[,2],lty=1) ; 
lines(x.ser,lgnm.pred[,3],lty=2)




#------------------   proportional odds property of logistic regression 

summary(lglg.fit)
sig <- lglg.fit$scale ; alpha <- 1/sig

mu <- predict(lglg.fit,newdata=data.frame(x=2.0),se.fit=F,type="lp") ; lambda <- exp(-mu)
surv.t.x2.0 <- 1 / (1+(lambda*300)^alpha)
odds.x2.0 <- surv.t.x2.0  / (1 - surv.t.x2.0)

 
mu <- predict(lglg.fit,newdata=data.frame(x=2.1),se.fit=F,type="lp") ; lambda <- exp(-mu)
surv.t.x2.1 <- 1 / (1+(lambda*300)^alpha)
odds.x2.1 <- surv.t.x2.1  / (1 - surv.t.x2.1)
 
mu <- predict(lglg.fit,newdata=data.frame(x=2.2),se.fit=F,type="lp") ; lambda <- exp(-mu)
surv.t.x2.2 <- 1 / (1+(lambda*300)^alpha)
odds.x2.2 <- surv.t.x2.2  / (1 - surv.t.x2.2)

odds.x2.0 / odds.x2.1


odds.x2.1 / odds.x2.2

exp(-0.1*lglg.fit$coef[2]*alpha)

 
#------------------   proportional times property of logistic regression 

summary(lglg.fit)

T.0.3.x2.0 <- predict(lglg.fit,newdata=data.frame(x=2.0),se.fit=F,type="quantile",p=0.3)
T.0.3.x2.1 <- predict(lglg.fit,newdata=data.frame(x=2.1),se.fit=F,type="quantile",p=0.3)

T.0.8.x2.1 <- predict(lglg.fit,newdata=data.frame(x=2.0),se.fit=F,type="quantile",p=0.8)
T.0.8.x2.2 <- predict(lglg.fit,newdata=data.frame(x=2.1),se.fit=F,type="quantile",p=0.8)

T.0.3.x2.0 / T.0.3.x2.1

T.0.8.x2.1 / T.0.8.x2.2


exp(-0.1*lglg.fit$coef[2])

#------------------   PH property of weibull regression 

summary(weib.fit)
sig <- weib.fit$scale ; alpha <- 1/sig


mu <- predict(weib.fit,newdata=data.frame(x=2.0),se.fit=F,type="lp") ; lambda <- exp(-mu)
haz.x2.0 <- lambda * alpha * (lambda*400)^(alpha-1)

mu <- predict(weib.fit,newdata=data.frame(x=2.1),se.fit=F,type="lp") ; lambda <- exp(-mu)
haz.x2.1 <- lambda * alpha * (lambda*400)^(alpha-1)

mu <- predict(weib.fit,newdata=data.frame(x=2.2),se.fit=F,type="lp") ; lambda <- exp(-mu)
haz.x2.2 <- lambda * alpha * (lambda*400)^(alpha-1)

haz.x2.1 / haz.x2.0


haz.x2.2 / haz.x2.1
 exp(-0.1*weib.fit$coef[2]*alpha)



 


#--------------------   Residuals


par(mfcol=c(1,3))
plot(x[ind],motorette[ind,2],xlab="x",ylab="survival time (hours)",ylim=range(motorette[,2]),
main="Time model")
points(x[!ind],motorette[!ind,2],pch=3)

lines(x.ser,weib.pred[,1],lty=2) 
lines(x.ser,weib.pred[,2],lty=1)  
lines(x.ser,weib.pred[,3],lty=2)

weib.res <- resid(weib.fit,type="deviance")
plot(x[ind],weib.res[ind],ylab="deviance residuals",xlab="x",ylim=range(weib.res))
points(x[!ind],weib.res[!ind],pch=3)
abline(0,0)

plot((1:30)[ind],weib.res[ind],ylab="deviance residuals",xlab="index",ylim=range(weib.res),xlim=c(1,30))
points((1:30)[!ind],weib.res[!ind],pch=3)
abline(0,0)