# Quasipoisson illustration
n <- 10000
g <- rexp(n, 1)
x <- runif(n)
y <- rpois(n, g*exp(x))

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

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

# Where do these SEs come from?
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))

# Quasipoisson sim
N <- 1000
covered <- function(fit, b0=1) {
  s <- coef(summary(fit))[2,]
  abs(s[1] - b0)/s[2] < qnorm(0.975)
}
if (!file.exists('cov.rds')) {
  Cov <- SE <- matrix(NA, N, 2, dimnames=list(1:N, c('Poisson', '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]
    saveRDS(Cov, 'cov.rds')
    saveRDS(SE, 'se.rds')
  }
} else {
  Cov <- readRDS('cov.rds')
  SE <- readRDS('se.rds')
}
 
Tab <- cbind(apply(Cov, 2, mean),
             apply(SE, 2, mean))
knitr::kable(Tab, 'latex', booktabs=TRUE, digits=3, col.names=c('Coverage', 'Average SE')) |>
  kableExtra::kable_styling(position = "center")

# 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 function: no cap
xx <- seq(-3, 3, len=99)
col <- c("#FF4E37FF", "#008DFFFF")
plot(xx, xx, type='l', lwd=3, bty='n', las=1, xlab='x', ylab=expression(psi[i](mu*'|'*x)), col=col[1])

# Influence function: Huber
plot(xx, xx, type='l', lwd=2, lty=2, bty='n', las=1, xlab='x', ylab=expression(psi[i](mu*'|'*x)), col=col[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=col[2])