###################################################
### chunk number 1: 
###################################################
library(lattice)
source("formats.R")
source("coxph.disparity.R")
options(width=65,show.signif.stars=FALSE)
options(SweaveHooks=list(fig=function() {
ltheme <- canonical.theme(color = FALSE)     ## in-built B&W theme
ltheme$strip.background$col <- "transparent" ## change strip bg
ltheme$plot.symbol$col <- "black"
ltheme$axis.text$cex <- 1.2
ltheme$par.xlab.text$cex <- 1.2
ltheme$par.ylab.text$cex <- 1.2
lattice.options(default.theme = ltheme)      ## set as default
}))
#library(xtable)
#library(lattice)
#ltheme <- canonical.theme(color = FALSE)     ## in-built B&W theme
#ltheme$strip.background$col <- "transparent" ## change strip bg
#ltheme$plot.symbol$col <- "black"
#ltheme$axis.text$cex <- 1.2
#ltheme$par.xlab.text$cex <- 1.2
#ltheme$par.ylab.text$cex <- 1.2
#lattice.options(default.theme = ltheme)      ## set as default
#ps.options(pointsize=14)


###################################################
### chunk number 2: 
###################################################
data(leuk,package="MASS")
feigl <- leuk
feigl <- transform(feigl,wbc=wbc/1000,ag=factor(ag,labels=c("+","-")))
feigl$ag <- relevel(feigl$ag,ref="-")


###################################################
### chunk number 3:  eval=FALSE
###################################################
## data(feigl,package="SMIR")


###################################################
### chunk number 4: feigl.plot1 eval=FALSE
###################################################
## print(xyplot(time~wbc,data=feigl))


###################################################
### chunk number 5: 
###################################################
print(xyplot(time~wbc,data=feigl))


###################################################
### chunk number 6: feigl.plot2 eval=FALSE
###################################################
## print(xyplot(time~wbc,data=feigl,xlab="wbc (log scale)",groups=ag,
##              ylab="time (log scale)", pch=c(3,1),col="black",
##              scales=list(x=list(log=TRUE,
##                            at=(xv <- c(2,5,10,20,40,60,80)),
##                            labels=xv),
##                y=list(log=TRUE,
##                  at=(yv <- c(2,5,10,25,50,100,150)),
##                  labels=yv))))


###################################################
### chunk number 7: 
###################################################
print(xyplot(time~wbc,data=feigl,xlab="wbc (log scale)",groups=ag,
             ylab="time (log scale)", pch=c(3,1),col="black",
             scales=list(x=list(log=TRUE,
                           at=(xv <- c(2,5,10,20,40,60,80)),
                           labels=xv),
               y=list(log=TRUE,
                 at=(yv <- c(2,5,10,25,50,100,150)),
                 labels=yv))))


###################################################
### chunk number 8: feigl.glm1
###################################################
feigl <- transform(feigl,lwbc=log(wbc))
feigl.glm1 <- glm(time~ag*lwbc,data=feigl,family=Gamma(link="log"))
summary(feigl.glm1,dispersion=1)


###################################################
### chunk number 9: feigl.glm3
###################################################
feigl.glm2 <- update(feigl.glm1,.~.-ag:lwbc)
print(coef(summary(feigl.glm2,dispersion=1)),digits=5)
feigl <- transform(feigl,z=ifelse(ag=="+",0,lwbc))
feigl.glm3 <- update(feigl.glm2,.~ag+z)
anova(feigl.glm2,feigl.glm3)
print(coef(summary(feigl.glm3,dispersion=1)),3)


###################################################
### chunk number 10: feigl.glm4
###################################################
feigl.glm4 <- update(feigl.glm1,family=Gamma)
summary(feigl.glm4,dispersion=1)


###################################################
### chunk number 11: feigl.glm5
###################################################
feigl.glm5 <- update(feigl.glm4,.~.-ag:lwbc)
summary(feigl.glm5,dispersion=1)


###################################################
### chunk number 12: feigl.glm6
###################################################
feigl.glm6 <- update(feigl.glm1,
family=Gamma(link="identity"))
#family=quasi(link="identity",variance="mu^2"))
summary(feigl.glm6,dispersion=1)


###################################################
### chunk number 13: feigl.glm7
###################################################
feigl.glm7 <- update(feigl.glm6,.~.-ag:lwbc)


###################################################
### chunk number 14: 
###################################################
feigl <- transform(feigl, u = time/fitted(feigl.glm1))


###################################################
### chunk number 15: 
###################################################
feigl <- transform(feigl, 
s=1-rank(u)/(length(u)+1))
#seq(1,length(u))/(length(u)+1))


###################################################
### chunk number 16: feigl.survival.plot eval=FALSE
###################################################
## print(xyplot(s~u,xlab="u",ylab="S(u)",data=feigl,
## panel=function(x,y){
## panel.xyplot(x,y)
## panel.curve(exp(-x),lwd=1.5)}))


###################################################
### chunk number 17: 
###################################################
print(xyplot(s~u,xlab="u",ylab="S(u)",data=feigl,
panel=function(x,y){
panel.xyplot(x,y)
panel.curve(exp(-x),lwd=1.5)}))


###################################################
### chunk number 18: 
###################################################
source('NPL.bands.R')


###################################################
### chunk number 19:  eval=FALSE
###################################################
## library("SMIR")


###################################################
### chunk number 20: feigl.bounds.plot eval=FALSE
###################################################
## print(xyplot((1-lower)+(1-upper)~x,
## data=NPL.bands(feigl$u),xlab="u",ylab="S(u)",
## panel=function(x,y){
## panel.xyplot(x,y,pch=19)
## panel.curve(exp(-x))
## }))


###################################################
### chunk number 21: 
###################################################
print(xyplot((1-lower)+(1-upper)~x,
data=NPL.bands(feigl$u),xlab="u",ylab="S(u)",
panel=function(x,y){
panel.xyplot(x,y,pch=19)
panel.curve(exp(-x))
}))


###################################################
### chunk number 22: feigl.bounds.log.plot eval=FALSE
###################################################
## print(xyplot(log(1-lower)+log(1-upper)~x,
## data=NPL.bands(feigl$u),xlab="u",ylab="log S(u)",
## panel=function(x,y){
## panel.xyplot(x,y,pch=19)
## panel.curve(-x)
## }))


###################################################
### chunk number 23: 
###################################################
print(xyplot(log(1-lower)+log(1-upper)~x,
data=NPL.bands(feigl$u),xlab="u",ylab="log S(u)",
panel=function(x,y){
panel.xyplot(x,y,pch=19)
panel.curve(-x)
}))


###################################################
### chunk number 24: feigl.bounds.log.log.plot eval=FALSE
###################################################
## print(xyplot((-log(-log(1-lower)))+(-log(-log(1-upper)))~log(x),
## data=NPL.bands(feigl$u),xlab="log u",ylab="-log[-log S(u)]",
## panel=function(x,y){
## panel.xyplot(x,y,pch=19)
## panel.curve(-x)
## }))


###################################################
### chunk number 25: 
###################################################
print(xyplot((-log(-log(1-lower)))+(-log(-log(1-upper)))~log(x),
data=NPL.bands(feigl$u),xlab="log u",ylab="-log[-log S(u)]",
panel=function(x,y){
panel.xyplot(x,y,pch=19)
panel.curve(-x)
}))


###################################################
### chunk number 26: feigl.leverage.plot eval=FALSE
###################################################
## print(dotplot(influence(feigl.glm1)$hat,type='h',
## scales=list(y=list(cex=0.8)),
## ylab="observation number", 
## xlab="leverage",
## panel=function(...){
## panel.dotplot(...)
## panel.abline(v=2*4/33,lty=2)}))


###################################################
### chunk number 27: 
###################################################
print(dotplot(influence(feigl.glm1)$hat,type='h',
scales=list(y=list(cex=0.8)),
ylab="observation number", 
xlab="leverage",
panel=function(...){
panel.dotplot(...)
panel.abline(v=2*4/33,lty=2)}))


###################################################
### chunk number 28: feigl.fitted.log.reciprocal.plots
###################################################
feigl.fitted.base <- xyplot(time~wbc,data=feigl,groups=ag,pch=c(3,1),
              scales=list(
                x=list(log="e",at=(xv <- c(1,2,5,10,20,50,100)),
                            labels=xv),
                y=list(log="e",at=(yv <- c(1,2,5,10,20,50,100,200,400)),
                  labels=yv)))
print(update(feigl.fitted.base,
              panel=function(...){
                panel.xyplot(...)
                panel.curve(4.7303-1.0176-0.3044*x,lwd=1.5,)
                panel.curve(4.7303-0.3044*x,lwd=1.5)
                panel.curve(log(1/(0.005795+0.006105*x)),lty=3,lwd=1.5)
                panel.curve(log(1/(0.005795+0.034415+0.006105*x)),lty=3,lwd=1.5)}))


###################################################
### chunk number 29: 
###################################################
feigl.fitted.base <- xyplot(time~wbc,data=feigl,groups=ag,pch=c(3,1),
              scales=list(
                x=list(log="e",at=(xv <- c(1,2,5,10,20,50,100)),
                            labels=xv),
                y=list(log="e",at=(yv <- c(1,2,5,10,20,50,100,200,400)),
                  labels=yv)))
print(update(feigl.fitted.base,
              panel=function(...){
                panel.xyplot(...)
                panel.curve(4.7303-1.0176-0.3044*x,lwd=1.5,)
                panel.curve(4.7303-0.3044*x,lwd=1.5)
                panel.curve(log(1/(0.005795+0.006105*x)),lty=3,lwd=1.5)
                panel.curve(log(1/(0.005795+0.034415+0.006105*x)),lty=3,lwd=1.5)}))


###################################################
### chunk number 30: feigl.boxcox
###################################################
library(MASS)
boxcox(time~ag*wbc,data=feigl,lambda=seq(-2,2,0.5),plotit=FALSE)


###################################################
### chunk number 31: feigl.boxcox2
###################################################
boxcox(time~ag*lwbc,data=feigl,lambda=seq(-2,2,0.5),plotit=FALSE)


###################################################
### chunk number 32: 
###################################################
-2*logLik(feigl.glm1)


###################################################
### chunk number 33: 
###################################################
feigl.lognormal <- update(feigl.glm1,log(time)~.,family=gaussian)


###################################################
### chunk number 34: 
###################################################
feigl.lognormal1 <- update(feigl.lognormal,.~.-ag:lwbc)


###################################################
### chunk number 35: feigl.lognormal.exp.fitted.plot eval=FALSE
###################################################
## print(update(feigl.fitted.base,
##               panel=function(...){
##                 panel.xyplot(...)
##                 panel.curve(4.8022-0.9883-0.5709*x,lty=3,lwd=1.5)
##                 panel.curve(4.8022-0.5709*x,lty=3,lwd=1.5)
##                 panel.curve(4.7303326-0.3044004*x,lty=1,lwd=1.5)
##                 panel.curve(4.7303326-1.0176452-0.3044004*x,lty=1,lwd=1.5)}))


###################################################
### chunk number 36: 
###################################################
print(update(feigl.fitted.base,
              panel=function(...){
                panel.xyplot(...)
                panel.curve(4.8022-0.9883-0.5709*x,lty=3,lwd=1.5)
                panel.curve(4.8022-0.5709*x,lty=3,lwd=1.5)
                panel.curve(4.7303326-0.3044004*x,lty=1,lwd=1.5)
                panel.curve(4.7303326-1.0176452-0.3044004*x,lty=1,lwd=1.5)}))


###################################################
### chunk number 37: linear.hazard.model eval=FALSE
###################################################
## summary(glm(feigl$w~model.matrix(feigl.glm1)*feigl$time-1,
## family=poisson(link="identity")))


###################################################
### chunk number 38: feigl.glm8
###################################################
feigl <- transform(feigl,w=rep(1,length(time)), lt=log(time))
feigl.glm8 <- glm(w~offset(lt)+ag*lwbc,data=feigl,family=poisson)
summary(feigl.glm8)


###################################################
### chunk number 39: reciprocal.hazard.model eval=FALSE
###################################################
## model.matrix(feigl.glm1)/feigl$time
## summary(glm(feigl$w~model.matrix(feigl.glm1)/feigl$time-1,
## family=poisson(link="reciprocal")))


###################################################
### chunk number 40: 
###################################################
plot(0:1,0:1,type="n",axes=FALSE,xlab="",ylab="Individual")
axis(2,at=seq(0,8)/8,label=c("","n",":","i",":",3,2,1,""),tck=FALSE,las=2,lty=0,font=3)
segments(0,0,0,1.1/8)
segments(0,1/8,0,5/8,lty=3)
segments(0,4.9/8,0,1)
axis(3,at=c(0:5)/6,labels=c(0,expression(a[1]),expression(a[2]),
expression(a[3]),expression(a[j]),expression(a[N])),tck=0.01,pos=1,lty=0,font=4)
segments(c(0:5)/6,6/6,c(1:6)/6,6/6,lty=c(1,1,1,3,3,1))
segments((1:5)/6,1,(1:5)/6,1-1/30)
mtext("Time point",side=3,line=3)
segments((0:2+0.1)/6,7/8,(1:3-0.1)/6,7/8)
text(1:3/6,7/8,c("0","0","X"))
segments((0:4+0.1)/6,6/8,(1:5-0.1)/6,6/8,lty=c(1,1,1,3,3))
text(c(1:5)/6,6/8,0)
segments((0:1+0.1)/6,5/8,(1:2-0.1)/6,5/8)
text(c(1:2)/6,5/8,0)
segments((0:4+0.1)/6,3/8,(1:5-0.1)/6,3/8,lty=c(1,1,1,3,3))
text(c(1:5)/6,3/8,c("0","0","0","0","X"))
segments((0:1+0.1)/6,1/8,(1:2-0.1)/6,1/8)
text(c(1:2)/6,1/8,c("0","X"))


###################################################
### chunk number 41:  eval=FALSE
###################################################
## data(gehan,package="SMIR")


###################################################
### chunk number 42: 
###################################################
data(gehan,package="MASS") 
gehan$treat <-relevel(gehan$treat,ref="control")


###################################################
### chunk number 43: gehan.table
###################################################
xtabs(~cens+time,data=gehan,subset=treat=="6-MP")


###################################################
### chunk number 44: gehan
###################################################
library(survival)
summary(survfit(Surv(time,cens),data=gehan,subset=treat=="6-MP",type="kap"))


###################################################
### chunk number 45: gehan.surv.plot1 eval=FALSE
###################################################
## library(survival)
## gehan.surv <- with(gehan,Surv(time,cens))
## gehan.survfit <- survfit(gehan.surv~treat,data=gehan,type="kap",conf.type="log-log")
## #plot(gehan.survfit,xlab="t",ylab="S(t)",lty=c(3,1),legend.text=levels(gehan$treat))
## print(xyplot(surv~time,data=summary(gehan.survfit),groups=gehan$treat,
## scales=list(y=list(at=seq(0,1,by=0.1))),
## ylab="S(u)",pch=19))


###################################################
### chunk number 46: 
###################################################
library(survival)
gehan.surv <- with(gehan,Surv(time,cens))
gehan.survfit <- survfit(gehan.surv~treat,data=gehan,type="kap",conf.type="log-log")
#plot(gehan.survfit,xlab="t",ylab="S(t)",lty=c(3,1),legend.text=levels(gehan$treat))
print(xyplot(surv~time,data=summary(gehan.survfit),groups=gehan$treat,
scales=list(y=list(at=seq(0,1,by=0.1))),
ylab="S(u)",pch=19))


###################################################
### chunk number 47: gehan.surv.plot2 eval=FALSE
###################################################
## print(xyplot(-log(-log(surv))~log(time),
## data=summary(gehan.survfit),groups=gehan$treat,
## scales=list(y=list(at=pretty(-log(-log(seq(0,1,by=0.1)))))),
## ylab="-log[-log S(u)]",xlab="log time",pch=19))


###################################################
### chunk number 48: 
###################################################
print(xyplot(-log(-log(surv))~log(time),
data=summary(gehan.survfit),groups=gehan$treat,
scales=list(y=list(at=pretty(-log(-log(seq(0,1,by=0.1)))))),
ylab="-log[-log S(u)]",xlab="log time",pch=19))


###################################################
### chunk number 49: gehan.glm
###################################################
gehan.glm <- glm(cens~offset(log(time))+treat,data=gehan,family=poisson)
round(coef(summary(gehan.glm)),3)


###################################################
### chunk number 50: gehan.bounds eval=FALSE
###################################################
## gehan$u <- fitted(gehan.glm)
## gehan.survfit <- survfit(Surv(u,cens)~1,data=gehan)
## print(xyplot((1-lower)+(1-upper)~x,
## data=NPL.bands(summary(gehan.survfit)$time[-19]),
## xlab="u",ylab="S(u)",
## panel=function(x,y){
## panel.xyplot(x,y,pch=20,col=1)
## panel.curve(exp(-x))
## panel.xyplot(summary(gehan.survfit)$time[-19],summary(gehan.survfit)$surv[-19])
## }))


###################################################
### chunk number 51: 
###################################################
gehan$u <- fitted(gehan.glm)
gehan.survfit <- survfit(Surv(u,cens)~1,data=gehan)
print(xyplot((1-lower)+(1-upper)~x,
data=NPL.bands(summary(gehan.survfit)$time[-19]),
xlab="u",ylab="S(u)",
panel=function(x,y){
panel.xyplot(x,y,pch=20,col=1)
panel.curve(exp(-x))
panel.xyplot(summary(gehan.survfit)$time[-19],summary(gehan.survfit)$surv[-19])
}))


###################################################
### chunk number 52: gehan.log.survivor.bounds eval=FALSE
###################################################
## print(xyplot(log(1-lower)+log(1-upper)~x,
## data=NPL.bands(summary(gehan.survfit)$time[-19]),
## xlab="u",ylab="log S(u)",
## panel=function(x,y){
## panel.xyplot(x,y,pch=20,col=1)
## panel.curve(-x)
## stmp <- summary(gehan.survfit)
## panel.xyplot(stmp$time[-19],log(stmp$surv[-19]),col="black")}))


###################################################
### chunk number 53: 
###################################################
print(xyplot(log(1-lower)+log(1-upper)~x,
data=NPL.bands(summary(gehan.survfit)$time[-19]),
xlab="u",ylab="log S(u)",
panel=function(x,y){
panel.xyplot(x,y,pch=20,col=1)
panel.curve(-x)
stmp <- summary(gehan.survfit)
panel.xyplot(stmp$time[-19],log(stmp$surv[-19]),col="black")}))


###################################################
### chunk number 54: gehan.log.log.survivor.bounds eval=FALSE
###################################################
## print(xyplot((-log(-log(1-lower)))+(-log(-log(1-upper)))~log(x),
## data=NPL.bands(summary(gehan.survfit)$time[-19]),
## xlab="log u",ylab="-log[-log S(u)]",
## panel=function(x,y){
## panel.xyplot(x,y,col="black",pch=20)
## panel.curve(-x)
## stmp <- summary(gehan.survfit)
## panel.xyplot(log(stmp$time[-19]),-log(-log(stmp$surv[-19])),col="black")}))


###################################################
### chunk number 55: 
###################################################
print(xyplot((-log(-log(1-lower)))+(-log(-log(1-upper)))~log(x),
data=NPL.bands(summary(gehan.survfit)$time[-19]),
xlab="log u",ylab="-log[-log S(u)]",
panel=function(x,y){
panel.xyplot(x,y,col="black",pch=20)
panel.curve(-x)
stmp <- summary(gehan.survfit)
panel.xyplot(log(stmp$time[-19]),-log(-log(stmp$surv[-19])),col="black")}))


###################################################
### chunk number 56: gamma.functions eval=FALSE
###################################################
## print(xyplot(1~1,xlim=c(-0.1,20),ylim=c(-0.01,0.4),ylab="density",xlab="t",
## panel=function(...){
## for (s in c(0.5,1,2,5)) panel.curve(dgamma(x,shape=s,rate=s/5),from=0.01,to=20,lty=s*2)}))


###################################################
### chunk number 57: 
###################################################
print(xyplot(1~1,xlim=c(-0.1,20),ylim=c(-0.01,0.4),ylab="density",xlab="t",
panel=function(...){
for (s in c(0.5,1,2,5)) panel.curve(dgamma(x,shape=s,rate=s/5),from=0.01,to=20,lty=s*2)}))


###################################################
### chunk number 58: gamma.hazard.functions eval=FALSE
###################################################
## print(xyplot(1~1,xlim=c(-0.2,20),ylim=c(-0.02,1),ylab="hazard",xlab="t",
## panel=function(...){
## for (s in c(0.5,1,2,5)) panel.curve(dgamma(x,shape=s,rate=s/5)/(1-pgamma(x,shape=s,rate=s/5))
## ,from=0.01,to=20,lty=s*2)}))


###################################################
### chunk number 59: 
###################################################
print(xyplot(1~1,xlim=c(-0.2,20),ylim=c(-0.02,1),ylab="hazard",xlab="t",
panel=function(...){
for (s in c(0.5,1,2,5)) panel.curve(dgamma(x,shape=s,rate=s/5)/(1-pgamma(x,shape=s,rate=s/5))
,from=0.01,to=20,lty=s*2)}))


###################################################
### chunk number 60: feigl.exp
###################################################
feigl.gamma <- glm(time~ag*lwbc,data=feigl,family=Gamma(link="log"))
summary(feigl.gamma)


###################################################
### chunk number 61: feigl.weibull
###################################################
library(MASS)
round(coef(summary(feigl.gamma,dispersion=gamma.dispersion(feigl.gamma))),4)
round(-2*logLik(feigl.gamma)[1],2)


###################################################
### chunk number 62: gehan.gamma
###################################################
gehan.gamma <- survreg(Surv(time,cens)~treat,data=gehan)
summary(gehan.gamma)
(gehan.gamma1 <- update(gehan.gamma,.~1))
round(exp(confint(gehan.gamma)[2,]))


###################################################
### chunk number 63: double.modelling
###################################################
library(dglm)
feigl.double.glm <- dglm(time~ag*lwbc,dformula=~ag*lwbc,family=Gamma(link="log"),
data=feigl,method="ml",dlink="log")
print(summary(feigl.double.glm),digits=3)


###################################################
### chunk number 64: feigl.dglm1
###################################################
feigl.double.glm1 <- update(feigl.double.glm,.~ag+lwbc,dformula=~lwbc)
print(summary(feigl.double.glm1),digits=3)


###################################################
### chunk number 65: 
###################################################
feigl.double.glm2 <- update(feigl.double.glm,.~ag+lwbc,dformula=~1)
print(summary(feigl.double.glm2),digits=3)


###################################################
### chunk number 66: 
###################################################
library(lattice)
print(xyplot(1:20~1:20,ylim=c(0,0.4),xlab="t",ylab="density",
panel=function(...){
panel.curve(dweibull(x,shape=0.5,scale=5),from=0.01,to=20)
panel.curve(dweibull(x,shape=1,scale=5),from=0.01,to=20,add=TRUE,lty='dotted')
panel.curve(dweibull(x,shape=2,scale=5),from=0.01,to=20,add=TRUE,lty='dotdash')
panel.curve(dweibull(x,shape=5,scale=5),from=0.01,to=20,add=TRUE,lty='dashed')}))


###################################################
### chunk number 67: 
###################################################
print(xyplot(0:20~0:20,ylim=c(-0.01,1.2),xlab="t",ylab="hazard",
panel=function(...){
panel.curve(dweibull(x,shape=0.5,scale=5)/(1-pweibull(x,shape=0.5,scale=5)),
from=0.01,to=20)
panel.curve(dweibull(x,shape=1,  scale=5)/(1-pweibull(x,shape=1,  scale=5)),
from=0.01,to=20,lty='dotted')
panel.curve(dweibull(x,shape=2,  scale=5)/(1-pweibull(x,shape=2,  scale=5)),
from=0.01,to=20,lty='dotdash')
panel.curve(dweibull(x,shape=5,  scale=5)/(1-pweibull(x,shape=5,  scale=5)),
from=0.01,to=20,lty='dashed')
}))


###################################################
### chunk number 68: gehan.weibull
###################################################
library(survival)
gehan.weibull <- survreg(Surv(time,cens)~treat,data=gehan,dist="weibull")
summary(gehan.weibull)
(gehan.exp <- survreg(Surv(time,cens)~treat,data=gehan,dist="exp"))
#print(gehan.exp)
(gehan.weibull1 <- update(gehan.weibull,.~1))
#summary(gehan.weibull1)


###################################################
### chunk number 69: feigl.weibull
###################################################
library(survival)
(feigl.exp <- survreg(Surv(time)~ag*lwbc,data=feigl,dist="exp"))
#print(feigl.exp)
feigl.weibull <- survreg(Surv(time)~ag*lwbc,data=feigl)
summary(feigl.weibull)
feigl.lognormal <- survreg(Surv(time)~ag*lwbc,data=feigl,dist="lognormal")
#feigl.lognormal$loglik[2]*(-2)


###################################################
### chunk number 70: 
###################################################
curve(exp(x-exp(x)),from=-5,to=3,xlab="y",ylab="density")


###################################################
### chunk number 71: gehan.extreme.glm
###################################################
library(survival)
summary(gehan.extreme <- survreg(Surv(time,cens)~treat,data=gehan,dist="extreme"))
summary(gehan.extreme1 <- update(gehan.extreme,.~1))


###################################################
### chunk number 72: 
###################################################
summary(feigl.extreme <- survreg(Surv(time)~ag*lwbc,data=feigl,dist="extreme"))


###################################################
### chunk number 73: gehan_logsurv eval=FALSE
###################################################
## ntimes <- (gehan$time*exp(-gehan.weibull$linear.predictors))^
## (1/gehan.weibull$scale)
## gehankm <- survfit(Surv(ntimes,gehan$cens))
## changes <- with(gehankm,c(-diff(n.risk),rev(n.event)[1]))
## #
## surv <- rev(gehankm$surv[gehankm$n.event!=0])[-1]
## riset <- rev(gehankm$n.risk[gehankm$n.event!=0])[-1]
## r <- riset*surv
## times <- rev(gehankm$time[gehankm$n.event!=0])[-1]
## tmp <- data.frame(NPL.bands(r,riset))
## print(xyplot(log(lower) + log(upper)~times,
##       data=tmp,
##              ylab="log S(u)", xlab="u", 
##              panel=function(x,y){
##                panel.xyplot(x,y,pch=20)
##                panel.points(x,log(surv))
##                panel.curve(-x)}))


###################################################
### chunk number 74: 
###################################################
ntimes <- (gehan$time*exp(-gehan.weibull$linear.predictors))^
(1/gehan.weibull$scale)
gehankm <- survfit(Surv(ntimes,gehan$cens))
changes <- with(gehankm,c(-diff(n.risk),rev(n.event)[1]))
#
surv <- rev(gehankm$surv[gehankm$n.event!=0])[-1]
riset <- rev(gehankm$n.risk[gehankm$n.event!=0])[-1]
r <- riset*surv
times <- rev(gehankm$time[gehankm$n.event!=0])[-1]
tmp <- data.frame(NPL.bands(r,riset))
print(xyplot(log(lower) + log(upper)~times,
      data=tmp,
             ylab="log S(u)", xlab="u", 
             panel=function(x,y){
               panel.xyplot(x,y,pch=20)
               panel.points(x,log(surv))
               panel.curve(-x)}))


###################################################
### chunk number 75: gehan_loglogsurv eval=FALSE
###################################################
## print(xyplot(-log(-log(lower)) + -log(-log(upper))~log(times),
##       data=tmp,
##              ylab="-log [- log S(u)]", xlab="u", 
##              panel=function(x,y){
##                panel.xyplot(x,y,pch=20)
##                panel.points(log(times),-log(-log(surv)))
##                panel.curve(-x)}))


###################################################
### chunk number 76: 
###################################################
print(xyplot(-log(-log(lower)) + -log(-log(upper))~log(times),
      data=tmp,
             ylab="-log [- log S(u)]", xlab="u", 
             panel=function(x,y){
               panel.xyplot(x,y,pch=20)
               panel.points(log(times),-log(-log(surv)))
               panel.curve(-x)}))


###################################################
### chunk number 77: 
###################################################
plot(0:1,0:1,type="n",axes=FALSE,xlab="",ylab="Individual")
axis(2,at=seq(0,8)/8,label=c("","n",":","i",":",3,2,1,""),tck=FALSE,las=2,lty=0,font=3)
segments(0,0,0,1.1/8)
segments(0,1/8,0,5/8,lty=3)
segments(0,4.9/8,0,1)
axis(3,at=c(0:5)/6,labels=c(0,expression(a[1]),expression(a[2]),
expression(a[3]),expression(a[j]),expression(a[N])),tck=0.01,pos=1,lty=0,font=4)
segments(c(0:5)/6,6/6,c(1:6)/6,6/6,lty=c(1,1,1,3,3,1))
segments((1:5)/6,1,(1:5)/6,1-1/30)
mtext("Time point",side=3,line=3)
segments((0:2+0.1)/6,7/8,(1:3-0.1)/6,7/8)
text(1:3/6,7/8,c("0","0","X"))
segments((0:4+0.1)/6,6/8,(1:5-0.1)/6,6/8,lty=c(1,1,1,3,3))
text(c(1:5)/6,6/8,0)
segments((0:1+0.1)/6,5/8,(1:2-0.1)/6,5/8)
text(c(1:2)/6,5/8,0)
segments((0:4+0.1)/6,3/8,(1:5-0.1)/6,3/8,lty=c(1,1,1,3,3))
text(c(1:5)/6,3/8,c("0","0","0","0","X"))
segments((0:1+0.1)/6,1/8,(1:2-0.1)/6,1/8)
text(c(1:2)/6,1/8,c("0","X"))
## plot(0:1,0:1,type="n",axes=FALSE,xlab="",ylab="Individual",cex.axis=1.5)
## axis(2,at=c(0,1,2,3,4,5,6)/6,label=c("","n","i",3,2,1,""),tck=FALSE,las=2,lty=0,font=3,cex.axis=1.2)
## segments(0/6,c(0:5)/6,0/6,c(1:6)/6,lty=c(1,3,3,1,1))
## axis(3,at=c(0:5)/6,labels=c(0,expression(a[1]),expression(a[2]),
## expression(a[3]),expression(a[j]),expression(a[N])),tck=0.01,pos=1,lty=0,font=3,cex.axis=1.2)
## mtext("Time point",side=3,line=3,cex=1.5)
## mtext("...",side=2,at=c(1.5,2.5)/6,line=1,cex=1.5)
## mtext("...",side=3,at=c(3.5,4.6)/6,line=0,cex=1.5)
## segments(c(3.48,3.52)/6,0.99,c(3.5,3.54)/6,1.01)
## segments(c(4.48,4.52)/6,0.99,c(4.5,4.54)/6,1.01)
## segments(c(0:5)/6,6/6,c(1:6)/6,6/6,lty=c(1,1,1,3,3,1))
## segments((1:5)/6,1,(1:5)/6,1-1/30)
## #
## segments((c(0,1,2)+0.1)/6,5/6,(c(1,2,2.5)-0.1)/6,5/6)
## text(c(1,2,2.5)/6,5/6,c("0","0","X"))
## segments((0:4+0.1)/6,4/6,(1:5-0.1)/6,4/6,lty=c(1,1,1,3,3))
## text(c(1:5)/6,4/6,"0")
## segments((0:2+0.1)/6,3/6,(c(1,2,2.3)-0.1)/6,3/6)
## text(c(1,2,2.3)/6,3/6,"0")
## segments((c(0,1,2,3,3.5)+0.1)/6,2/6,(c(1,2,3,3.5,4)-0.1)/6,2/6,lty=c(1,1,1,3,1))
## text(c(1,2,3,3.5,4)/6,2/6,c("0","0","0","0","X"))
## segments((0:1+0.1)/6,1/6,(c(1,1.5)-0.1)/6,1/6)
## text(c(1,1.5)/6,1/6,c("0","X"))


###################################################
### chunk number 78: gehan.cph
###################################################
library(survival)
gehan.cph <- coxph(Surv(time,cens)~treat,data=gehan,method='breslow')
print(gehan.cph)


###################################################
### chunk number 79: coxph.disparity
###################################################
#tind <- function(t,ddt){
#    ind <- NULL
#    for ( i in 1:length(t)) ind <- c(ind,sum(ddt<=t[i]))
#    ind
#}
gz <- basehaz(gehan.cph,centered=TRUE)
gz <- transform(gz,
                haz = diff(c(0,hazard)),
                tint = diff(c(0,time)))
ti <- apply(outer(gehan$time,gz$time,">="),1,sum)
#ti <- tind(gehan$time,gz$time)
gehan.cph.disparity <- -2*sum(gehan$cens*(log(gz$haz[ti]/gz$tint[ti]) +
                   gehan.cph$linear)-
       gz$hazard[ti]*exp(gehan.cph$linear))


###################################################
### chunk number 80: gehan.baseline.plot eval=FALSE
###################################################
## plot(-log(hazard)~time,data=basehaz(gehan.cph,centered=FALSE),log="x",
## xlab="time (log scale)", ylab= "-log[-log(S(t))]",las=1)


###################################################
### chunk number 81: 
###################################################
plot(-log(hazard)~time,data=basehaz(gehan.cph,centered=FALSE),log="x",
xlab="time (log scale)", ylab= "-log[-log(S(t))]",las=1)


###################################################
### chunk number 82:  eval=FALSE
###################################################
## data(feigl,package="SMIR")


###################################################
### chunk number 83:  eval=FALSE
###################################################
## feigl <- transform(feigl,lwbc=log(wbc))
## (feigl.cph <- coxph(Surv(time)~ag*lwbc,data=feigl,method='breslow'))
## library(SMIR)
## round(coxph.disparity(feigl.cph),2)


###################################################
### chunk number 84: feigl.cph
###################################################
feigl <- transform(feigl,lwbc=log(wbc))
(feigl.cph <- coxph(Surv(time)~ag*lwbc,data=feigl,method='breslow'))
round(coxph.disparity(feigl.cph),2)


###################################################
### chunk number 85: feigl.baseline.plot eval=FALSE
###################################################
## plot(-log(hazard)~time,data=basehaz(feigl.cph,centered=FALSE),log="x",
## xlab="time (log scale)", ylab= "-log[-log(S(t))]",las=1)


###################################################
### chunk number 86: 
###################################################
plot(-log(hazard)~time,data=basehaz(feigl.cph,centered=FALSE),log="x",
xlab="time (log scale)", ylab= "-log[-log(S(t))]",las=1)


###################################################
### chunk number 87: feigl.cph1
###################################################
print(feigl.cph1 <- update(feigl.cph,.~.-ag:lwbc),digits=4)
#summary(feigl.cph1)
anova(feigl.cph1,feigl.cph)


###################################################
### chunk number 88: gehan.logistic
###################################################
library(survival)
(gehan.logistic <- survreg(Surv(time,cens)~treat,data=gehan,dist="logistic"))
round(gehan.logistic$loglik[2],2)
#summary(gehan.logistic)


###################################################
### chunk number 89: gehan.loglogistic
###################################################
gehan.loglogistic <- update(gehan.logistic,dist="loglogistic")
summary(gehan.loglogistic)
round(gehan.loglogistic$loglik[2],2)
(gehan.loglogistic.null <- update(gehan.loglogistic,.~1))
#summary(gehan.loglogistic.null)


###################################################
### chunk number 90: 
###################################################
library(survival)
(gehan.normal <- survreg(Surv(time,cens)~treat,data=gehan,dist="gaussian"))
round(gehan.normal$loglik[2],2)
#summary(gehan.normal)
gehan.lognormal <- update(gehan.normal,dist="lognormal")
summary(gehan.lognormal)
round(gehan.lognormal$loglik[2],2)


###################################################
### chunk number 91: gehan.lognormal1
###################################################
gehan.lognormal1 <- update(gehan.lognormal,.~1)
print(anova(gehan.lognormal1,gehan.lognormal),digits=4)
#summary(gehan.lognormal1)


###################################################
### chunk number 92:  eval=FALSE
###################################################
## data(prentice,package="SMIR")
## library(survival)


###################################################
### chunk number 93: 
###################################################
library(survival)
prentice <- read.csv('Data/prentice.csv',header=TRUE)
prentice <- transform(prentice,
                      Type = factor(prentice$type),
                      status = status/10 )


###################################################
### chunk number 94: prentice.survivor.fn eval=FALSE
###################################################
## prentice.survfit <- survfit(Surv(time,censor)~Type,data=prentice)
## print(xyplot(-log(-log(surv))~log(time),groups=strata,
##              summary(prentice.survfit),
##        ylab=expression("-log[-log S(t)]"),xlab="log t"))


###################################################
### chunk number 95: 
###################################################
prentice.survfit <- survfit(Surv(time,censor)~Type,data=prentice)
print(xyplot(-log(-log(surv))~log(time),groups=strata,
             summary(prentice.survfit),
       ylab=expression("-log[-log S(t)]"),xlab="log t"))


###################################################
### chunk number 96: prentice.cph1
###################################################
prentice <- transform(prentice,lmfd=log(mfd),
Prior=factor(prior),Treat=factor(treat),Type=factor(type))
library(survival)
prentice.type1 <- subset(prentice,Type=="1")
prentice.cph1 <- coxph(Surv(time,censor)~status+lmfd+age+Prior+Treat,
data=prentice.type1,method='breslow')
print(summary(prentice.cph1)$coef,digits=3)


###################################################
### chunk number 97: type1.survivor.fn eval=FALSE
###################################################
## #with(survfit(prentice.cph1),plot(-log(-log(surv))~log(time)))
## plot(-log(hazard)~time,data=basehaz(prentice.cph1,centered=FALSE),log="x",
## xlab="time (log scale)", ylab= "-log(-log S)")


###################################################
### chunk number 98: 
###################################################
#with(survfit(prentice.cph1),plot(-log(-log(surv))~log(time)))
plot(-log(hazard)~time,data=basehaz(prentice.cph1,centered=FALSE),log="x",
xlab="time (log scale)", ylab= "-log(-log S)")


###################################################
### chunk number 99: prentice.cph2
###################################################
prentice.cph2 <- update(prentice.cph1,.~.^2)
#
#coxph(Surv(time,censor)~(status+lmfd+age+Prior+Treat)^2,
#data=prentice.type1)
#summary(prentice.cph2)
#drop1(prentice.cph2,test="Chisq")
prentice.cph3 <- update(prentice.cph2,.~status+lmfd+Prior+Prior:status)
#summary(prentice.cph3)
#summary(prentice.weibull2 <- survreg(Surv(time,censor)~(status+lmfd+age+Prior+Treat)^2,data=prentice.type1))


###################################################
### chunk number 100: prentice.weibull
###################################################
prentice.weibull <- survreg(Surv(time,censor)~lmfd+status*Prior,data=prentice.type1)
print(summary(prentice.weibull),digits=3)
round(prentice.weibull$loglik[2],digits=2)
#anova(prentice.weibull,prentice.cph2)


###################################################
### chunk number 101: prentice.t234.cph
###################################################
(prentice.t234.cph <- coxph(Surv(time,censor)~status+lmfd+age+Prior+Type+Treat,
data=prentice,subset=Type!="1",method='breslow'))


###################################################
### chunk number 102: 
###################################################
with(survfit(prentice.t234.cph),plot(-log(-log(surv))~time,log="x",
xlab="time (log scale)", ylab= "-log(-log S)"))


###################################################
### chunk number 103: 
###################################################
(prentice.t234.cph20 <- update(prentice.t234.cph,.~.^2))


###################################################
### chunk number 104: prentice.t234.cph2
###################################################
(prentice.t234.cph2 <- update(prentice.t234.cph,
.~status*Prior+Treat+I(Type=="4")))
round(coxph.disparity(prentice.t234.cph2),2)


###################################################
### chunk number 105: prentice.t234.weibull
###################################################
prentice.t234.weibull <- survreg(Surv(time,censor)~status*Prior+Treat+I(Type=="4"),
data=prentice,dist="weibull",subset=Type!="1")
summary(prentice.t234.weibull)
round(prentice.t234.weibull$loglik[2],2)


###################################################
### chunk number 106: 
###################################################
stan <- read.table("Data/stan.txt",skip=9,nrows=65,sep="")
names(stan) <- c("id","za","zb","age","surg","acc","died","surv","nmm","hla","mm","rej")
stan$surv[stan$surv==0] <- 0.5


###################################################
### chunk number 107:  eval=FALSE
###################################################
## data(stan,package="SMIR")


###################################################
### chunk number 108: stan.cph
###################################################
library(survival)
(stan.cph <- coxph(Surv(surv,died)~age+surg+mm,data=stan,method="breslow"))


###################################################
### chunk number 109: stan.log.hazard
###################################################
#with(survfit(stan.cph),plot(-log(-log(surv))~log(time)))
plot(-log(hazard)~time,basehaz(stan.cph,center=F),
     log="x",las=1,xlab="time (log scale)",ylab="-log(-log S)")


###################################################
### chunk number 110: 
###################################################
#with(survfit(stan.cph),plot(-log(-log(surv))~log(time)))
plot(-log(hazard)~time,basehaz(stan.cph,center=F),
     log="x",las=1,xlab="time (log scale)",ylab="-log(-log S)")


###################################################
### chunk number 111: 
###################################################
stan <- transform(stan, drej = rej*died)
(stan.reject.cph <- coxph(Surv(surv,drej)~age+surg+mm,data=stan,method="breslow"))


###################################################
### chunk number 112: 
###################################################
stan.reject.cph2 <- update(stan.reject.cph,.~.-surg)
#1.0342/sqrt(3.548)


###################################################
### chunk number 113: stan.transplant.rejection eval=FALSE
###################################################
## plot(-log(hazard)~time,basehaz(stan.reject.cph,centered=FALSE),
##      log="x",las=1,xlab="time (log scale)",ylab="-log(-log S)")


###################################################
### chunk number 114: 
###################################################
plot(-log(hazard)~time,basehaz(stan.reject.cph,centered=FALSE),
     log="x",las=1,xlab="time (log scale)",ylab="-log(-log S)")


###################################################
### chunk number 115: stan.nrej.cph
###################################################
stan <- transform(stan, nrej = (1-rej)*died)
(stan.nrej.cph <- coxph(Surv(surv,nrej)~age+surg+mm,data=stan,method="breslow"))


###################################################
### chunk number 116: stan.nrej.hazard eval=FALSE
###################################################
## plot(-log(hazard)~time,data=basehaz(stan.nrej.cph,centered=FALSE),
##      log="x",las=1,xlab="time (log scale)" )


###################################################
### chunk number 117: 
###################################################
plot(-log(hazard)~time,data=basehaz(stan.nrej.cph,centered=FALSE),
     log="x",las=1,xlab="time (log scale)" )


###################################################
### chunk number 118: stan.nrej.null
###################################################
(stan.nrej.null <- coxph(Surv(surv,nrej)~1,data=stan,method="breslow"))
round(coxph.disparity(stan.nrej.null),2)


###################################################
### chunk number 119: stan.weibull
###################################################
stan.weibull <- survreg(Surv(surv,nrej)~1,data=stan,subset=surv>0,dist="weibull")
summary(stan.weibull)
round(-2*stan.weibull$loglik[2],2)


###################################################
### chunk number 120: 
###################################################
(stan.weibull3 <- update(stan.weibull,.~age+surg+mm))


###################################################
### chunk number 121: gehan.discrete
###################################################
gehand <- data.frame(d = c(7,6,4,1,1,2,0,4,1,2,0,2),
                    r = c(21,14,8,4,3,2,21,20,13,12,7,7),
                    month = gl(6,1),
                    treat = gl(2,6,labels=c("control","6-MP")))
library(gnm)
gehand.gnm <- gnm(cbind(d,(r-d))~1,eliminate=month,data=gehand,
family=binomial(link="cloglog"),verbose=F)
(gehand.gnm1 <- update(gehand.gnm,.~.+treat))
anova(gehand.gnm,gehand.gnm1)


###################################################
### chunk number 122: feigl.discrete
###################################################
#feigl <- read.csv("Data/feigl.csv")
#library(survival)
#feigl.survfit <- survfit(Surv(time,w),data=feigl)
#j <- 1
Nn <- sum(feigl$time)
I <- nrow(feigl)
Feigl <- feigl[rep(1:I,feigl$time),]
#Feigl$y <- rep(0,Nn)
vP <- NULL
for ( i in 1:I)
vP <- c(vP,seq(1,feigl$time[i],by=1))
Feigl <- transform(Feigl,
                   y = ifelse(time==vP,1,0),
                   lwbc=log(wbc),
                   Period=factor(vP))
feigl.discrete.glm <- glm(y~Period+ag*lwbc-1,data=Feigl,family=binomial(link="cloglog"))
print(summary(feigl.discrete.glm)$coef[-c(1:156),],digits=3)


###################################################
### chunk number 123: feigld.glm1
###################################################
feigl.discrete.glm1 <- update(feigl.discrete.glm,.~.-ag:lwbc)
print(summary(feigl.discrete.glm1)$coef[-c(1:156),],digits=4)
#round(feigld.glm1$deviance,2)
anova(feigl.discrete.glm1,feigl.discrete.glm,test="Chisq")