source("http://myweb.uiowa.edu/pbreheny/7210/f19/notes/fun.R")
library(survival)

# Slide 22
Data <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/Pike1966.txt")
d <- sum(Data$Death)
v <- sum(Data$Time/365.25)
fit <- survfit(Surv(Time/365.25, Death) ~ 1, Data)
plot(fit, mark.time=FALSE, col=pal(2)[1], lwd=3, bty='n', las=1,
     xlab='Time on study (Years)', ylab='Survival', conf.int=FALSE)
tt <- seq(0, 1, len=199)
lines(tt, pexp(tt, d/v, low=FALSE), lwd=3, col=pal(2)[2])

# Slides 23-27: Setup
lam <- seq(0.5, 3, len=199)
l <- function(ll) {d*log(ll) - ll*v - d*log(d/v) + d}
#l <- function(ll) {d*log(ll) - ll*v}
baseplot <- function() {plot(lam, l(lam), col=pal(2)[1], lwd=3, type='l', bty='n',
                             xlab=expression(lambda), ylab=expression(loglik(lambda)))}
baseplot()

# Slide 23-24: Score
baseplot()
lines(c(0, 2), l(1) + c(-1,1) * (d-v), lwd=3, col='gray')
mtext("Slope: 11")
d - v
sqrt(d)
(d-v)/sqrt(d)

# Slide 25-26: Wald
baseplot()
arrows(1, 0, d/v, 0, lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(hat(theta)-1: " "*0.44))
baseplot()
arrows(1, 0, d/v, 0, lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(hat(theta)-1: " "*0.44))
lines(lam, -(lam-d/v)^2/(2*d/v^2), col=pal(2)[2], lwd=3)
d/v-1
sqrt(d)/v
(d/v-1)/(sqrt(d)/v)

# Slide 27-28: LRT
baseplot()
arrows(d/v, 0, d/v, l(1), lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(loglik(hat(theta))-loglik(1): " "*2.14))
baseplot()
arrows(d/v, 0, d/v, l(1), lwd=3, col='gray', angle=90, code=3, xpd=TRUE)
mtext(expression(loglik(hat(theta))-loglik(1): " "*2.14))
abline(h=-qchisq(.95,1)/2, col=pal(2)[2], lwd=3)
l(1)
pchisq(-2*l(1), 1, low=FALSE)

# Slide 29: Exact
vv <- seq(0, 60, len=199)
plot(vv, dgamma(vv, shape=d, rate=1), type='l', lwd=3, bty='n', col=pal(2)[1],
     xlab='v', ylab='Density', las=1)
abline(v=v, lwd=3, col='gray')
pgamma(v, shape=d, rate=1)

# Comparison of tests
2*pnorm((d-v)/sqrt(d), lower.tail = FALSE)   # Wald/Score
pchisq(-2*l(1), 1, low=FALSE)                # LRT
2*pgamma(v, shape=d, rate=1)                 # Exact