source("http://myweb.uiowa.edu/pbreheny/7210/f18/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)

# Slide 31: Bayes, reference prior
lam <- sort(c(seq(0, 3, len=199), 1))
par(mar=c(4.5,5,0.1,0.1))
plot(lam, dgamma(lam, 1, 0), col='gray', lwd=3, type='l', bty='n', ylim=c(0, 1.75),
     xlab=expression(lambda), ylab=expression(p(lambda)))
par(mar=c(4.5,5,0.1,0.1))
plot(lam, dgamma(lam, 1, 0), col='gray', lwd=3, type='l', bty='n', ylim=c(0, 1.75),
     xlab=expression(lambda), ylab=expression(p(lambda)))
lines(lam, dgamma(lam, d+1, v), col=pal(2)[2], lwd=3)
pgamma(1, d+1, v)

# Slide 31: Bayes, informative prior
par(mar=c(4.5,5,0.1,0.1))
plot(lam, dgamma(lam, 3, 3), col='gray', lwd=3, type='l', bty='n', ylim=c(0, 1.75),
     xlab=expression(lambda), ylab=expression(p(lambda)))
lines(lam, dgamma(lam, 2+d, 2+v), col=pal(2)[2], lwd=3)
pgamma(1, d+3, v+3)