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)