source('https://myweb.uiowa.edu/pbreheny/7110/f20/notes/fun.R')

# Likelihood
l <- function(p, n1, n0, p1, p0) {
  n1*log(p*p1) + n0*log((1-p)*p0) - (n0+n1)*log(p*p1 + (1-p)*p0)
}
pp <- seq(0, 0.6, len=501)
ll <- l(pp, 20, 40, 1, 0.5)
plotL(pp, exp(ll-max(ll)))

# Pseudo
pseudo_l <- function(p, n1, n0, p1, p0) {
  n1/p1*log(p) + n0/p0*log(1-p)
}
pll <- pseudo_l(pp, 20, 40, 1, 0.5)
plotL(pp, exp(pll-max(pll)), col=pal(2)[1], add=TRUE)
rightlegend(legend=c('True', 'Pseudo'), lwd=3, col=pal(2)[2:1])

# Score
u <- function(p, n1, n0, p1, p0) {
  n1/p - n0/(1-p) - (n0+n1)*(p1-p0)/(p*p1+(1-p)*p0)
}
pp <- seq(0.05, 0.5, len=101)
plot(pp, u(pp, 20, 40, 1, 0.5), type='l', lwd=3, col=pal(2)[2], las=1, bty='n', xlab=expression(pi), ylab='Score', ylim=c(-100, 300))
abline(h=0, lty=2, col='gray')

# Pseudo-score
pseudo_u <- function(p, n1, n0, p1, p0) {
  (n1/p1)/p - (n0/p0)/(1-p)
}
lines(pp, pseudo_u(pp, 20, 40, 1, 0.5), lwd=3, col=pal(2)[1])
toplegend(legend=c('True', 'Pseudo'), lwd=3, col=pal(2)[2:1])

# Sample data
p1 <- 1
p0 <- 0.5
n1 <- 20
n0 <- 40
p <- n1*p0 /(n1*p0 + n0*p1)

# Compare information / variance (true/pseudo/obs/exp)
IT <- n1/p^2 + n0/(1-p)^2 - (n0+n1)*(p1-p0)^2/(p*p1+(1-p)*p0)^2
1/IT
IP <- n1/(p1*p^2) + n0/(p0*(1-p)^2)
1/IP
B <- (1-p1)/p1^2 * n1/p^2 + (1-p0)/p0^2 * n0/(1-p)^2  # Kalbfleish
1/IP + B/IP^2
V <- n1/(p1*p)^2 + n0/(p0*(1-p))^2  # Robust obs
V/IP^2
EIP <- 1/p + 1/(1-p)  # Robust Fisher
EV <- 1/(p1*p) + 1/(p0*(1-p))
EV/EIP^2 * (p1*p + p0*(1-p))/(n0+n1)

# Confidence intervals
SE1 <- sqrt(1/IT)
SE2 <- sqrt(1/IP)
SE3 <- sqrt(V/IP^2)
p + qnorm(c(.025, .975)) * SE1
p + qnorm(c(.025, .975)) * SE2
p + qnorm(c(.025, .975)) * SE3