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

n <- 100
g <- rexp(n, 1)
x <- runif(n)
y <- rpois(n, g*exp(x))

# Pois
fit1 <- glm(y ~ x, family=poisson())
summary(fit1)

# Quasi
fit2 <- glm(y ~ x, family=quasipoisson())
summary(fit2)

# Derive our own
X <- model.matrix(fit1)
eta <- drop(X %*% coef(fit1))
mu <- exp(eta)
phi <- sum((y-mu)^2/mu) / (n-2)
w <- mu/phi
vcov(fit2)
solve(crossprod(X, w*X))

# Small sim
N <- 1000
covered <- function(fit, b0=1) {
  s <- coef(summary(fit))[2,]
  abs(s[1] - b0)/s[2] < qnorm(0.975)
}
Cov <- SE <- matrix(NA, N, 2, dimnames=list(1:N, c('Poission', 'Quasi')))
for (i in 1:N) {
  n <- 100
  g <- rexp(n, 1)
  x <- runif(n)
  y <- rpois(n, g*exp(x))
  fit1 <- glm(y ~ x, family=poisson())
  fit2 <- glm(y ~ x, family=quasipoisson())
  Cov[i,1] <- covered(fit1)
  Cov[i,2] <- covered(fit2)
  SE[i, 1] <- coef(summary(fit1))[2,2]
  SE[i, 2] <- coef(summary(fit2))[2,2]
}
Tab <- cbind(apply(Cov, 2, mean),
             apply(SE, 2, mean))

# Sandwich estimation
library(sandwich)
fit1 <- glm(y ~ x, family=poisson())
fit2 <- glm(y ~ x, family=quasipoisson())
vcov(fit1)
vcov(fit2)
vcovCL(fit1)
vcovCL(fit2)

# Influence functions
xx <- seq(-3, 3, len=99)
plot(xx, xx, type='l', lwd=3, bty='n', las=1, xlab='x', ylab=expression(psi[i](mu*'|'*x)), col=pal(2)[1])
plot(xx, xx, type='l', lwd=2, lty=2, bty='n', las=1, xlab='x', ylab=expression(psi[i](mu*'|'*x)), col=pal(2)[1])
lines(xx, sign(xx)*pmin(abs(xx), 1), type='l', lwd=3, bty='n', las=1, xlab='x', ylab=expression(psi[i](x)), col=pal(2)[2])