code for understanding firth's bias correction for logit models

soren.R

library(ggplot2)
library(MASS)
library(logistf)

set.seed(12345678)
############ FAKE DATA ############
rho <- 0.0
k <- 3
sigma <- matrix(rho, k, k); diag(sigma) <- 5
n <- 500
b0 <- 0
b1 <- 0.5
n.sims <- 1
X <- cbind(1, mvrnorm(n, rep(0, k), Sigma = sigma))
b <- c(b0, b1, rep(0.2, k - 1))
p <- plogis(X%*%b)
# Gen y, fit the model
y <- rbinom(n, 1, p)
model <- glm(y ~ X - 1, family = "binomial")
summary(model)
model2 <- logistf(y ~ X - 1)
summary(model2)

# eta, pi, and W
eta <- X%*%model$coef
pi <- exp(X%*%model$coef)/(1 + exp(X%*%model$coef))
W <- as.vector((pi*(1-pi)))*diag(length(X[,1]))
# also W
# W <- exp(X%*%model$coef)/(1 + exp(X%*%model$coef))*
#   (1-exp(X%*%model$coef)/(1 + exp(X%*%model$coef)))
# Information.
I.b <- t(X)%*%W%*%X
I.b.inv <- solve(I.b)
sqrt(diag(I.b.inv)) # standard errors!
cbind(sqrt(diag(I.b.inv)), summary(model)$coef[,2])
# ‘‘Asymptotic variance of X/hat matrix’’
H <- sqrt(W)%*%X%*%solve(t(X)%*%W%*%X)%*%t(X)%*%sqrt(W)
# Lowercase h = diagonal of H
h <- d.H <- diag(H)
# Quantity of interest: bias correcting vector W.xi
W.xi <- d.H*(pi - 0.5)
# Can actually compute bias
bias <- solve(t(X)%*%W%*%X)%*%t(X)%*%W.xi
bias
# Estimates the same (bias correcting versus penalized)
cbind(summary(model)$coef[,1] - bias, summary(model2)$coef)
# Let’s call that rounding/approximation
data.plot <- data.frame(X[,2], X[,3], y, eta, pi, diag(W), d.H, W.xi)
names(data.plot) <- c("X2", "X3", "y", "eta", "pi", "diagW", "dH", "Wxi")
aa <- ggplot(data.plot, aes(eta, y))
aa + geom_point(aes(color = abs(Wxi)), size = 3)
bb <- ggplot(data.plot, aes(eta, pi))
bb + geom_point(aes(color = abs(Wxi)), size = 3)
cc <- ggplot(data.plot, aes(X2, eta))
cc + geom_point(aes(color = abs(Wxi)), size = 3)
dd <- ggplot(data.plot, aes(X3, eta))
dd + geom_point(aes(color = abs(Wxi)), size = 3)

# plot the correction again pi
ee <- ggplot(data.plot, aes(pi, Wxi)) + geom_point(size = 3); ee

# plot the correction again pi
ff <- ggplot(data.plot, aes(eta, Wxi)) + geom_point(size = 3); ff

# plot the second derivative of the link function
x <- seq(-3, 3, length.out = 100)
d <- numeric(length(x))
for (i in 1:length(x)) {
  d[i] <- grad(dlogis, x[i])
}
qplot(plogis(x), d)