Prediction interval for sum of observations from logistic regression

Question

It is possible make prediction intervals for $Y\mid X=x$ where $X \sim \operatorname{Binom}(n,p)$ and $Y \sim \operatorname{Binom}(m,p)$ as answered in my previous question .

Now I wish to make a prediction interval for $Y$ in a logistic regression settings. Yes, making prediction intervals for single observations in logistic regression is not very useful. Rather, the goal is to create a prediction interval for the aggregate sum.

Here is an example of what I am trying to do. Assume we have a binary response, $z$, and three binary predictors $A,B,C$. The model $\operatorname{logit}(z) \sim A + B + C$ is fit (only main effects) using training data with both the response and predictors available which has $n=100$ observations of various levels of $A,B,C$. Next, someone hands us $m=50$ new rows of data that only contains $A,B,C$ e.g...

$$ \begin{array}{cccl} A & B & C \\[6pt] 1 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 1 & 0 \\ & \vdots \\ & \vdots & & \text{46 rows ommitted} \\ 0 & 0 & 0 \end{array} $$

Our goal is to predict how many $Z_i=1$ for these 50 new observations (i.e. the aggregate sum $Y=\sum_{i=1}^m{Z}_i$). We can get a point estimate by summing the estimated probabilities for each new observation, but how would we generate a prediction interval? Since the sum is not Binomial ($p$ varies), there doesn't seem to be an answer available.

Answer 1

Bootstrapping is indeed attractive. This procedure consists of drawing with replacement from the data, fitting the logistic model to obtain a new parameter estimate $\hat\beta$, and using that to compute upper and lower prediction limits for the sum of new responses. Repeating it will produce the (bivariate) bootstrap distribution of the prediction limits.

Obtaining the prediction limits itself requires some simulation. In this case, with three categorical regressors, that can be done efficiently by noting only eight combinations are possible, and therefore the $50$ new responses can be grouped into eight groups and summed within each group (to produce a Binomial response). Thus, simulating a single sum of all response requires summing just eight Binomial variates. Doing this enough times will pin down the prediction limits with sufficient precision.

Here are some illustrations. The data consist of all eight possible combinations of $A$, $B$, and $C$ (binary coded), each repeated 24 times for a total of 192 observations. The true model is that the logit of the response equals $A-C$. Indeed, the coefficients of $A$ and $C$ are significant in the actual data while the other two coefficients (intercept and $B$) are not. For each of $500$ bootstrap iterations, $200$ sums-of-responses were obtained in order to estimate a $95\%$ two-sided prediction interval $(\text{LPL},\text{UPL}]$. The correct prediction interval, based on the true model, is $(18, 31]$.

The histograms show the marginal distributions of the prediction limits: LPL in blue, UPL in red. The correct limits of $18$ and $31$ are marked by vertical colored lines. The scatterplot indicates the bivariate relationship of the prediction limits. The mean width $\text{UPL}-\text{LPL}$ is $12.2$, plotted as the dark red line.

The agreement of these results with the true values is encouraging. Iteration of this short experiment (using varying random number seeds) provides some intuition about how much can be expected from a dataset like this. Change the datasets to study any other situation in the same manner.

beta <- c(0,1,0,-1) # Model parameters
n.rep <- 24         # Replications of each (A,B,C) combination
alpha <- 0.05       # Predict 1 - alpha of the response
#
# Create data.
#
#set.seed(17)
x <- matrix(rep(as.matrix(expand.grid(A=0:1, B=0:1, C=0:1)), each=n.rep), ncol=3)
colnames(x) <- c("A","B","C")
nvars <- length(beta) - 1

logistic <- function(x) 1 / (1 + exp(-x))
Y <- rbinom(nrow(x), 1, logistic(x %*% beta[-1] + beta[1]))
df <- as.data.frame(cbind(x, Y))
#
# Fit the model.
#
fit <- glm(cbind(Y,1-Y) ~ ., df, family=binomial)
summary(fit)
#
# Generate new data.
#
n.new <- 50
df.new <- as.data.frame(matrix(sample.int(2, 3*n.new, replace=TRUE)-1, ncol=nvars))
names(df.new) <- names(df)[1:nvars]
#
# Obtain the true distribution of the sum of their responses.
#
PL <- function(beta, df.new, alpha, n.iter){
  p <- logistic(as.matrix(df.new) %*% beta[-1] + beta[1])
  probs <- table(p)
  dist <- colSums(matrix(rbinom(length(probs)*n.iter, probs, as.numeric(names(probs))),
                         ncol=n.iter, byrow=FALSE))
  quantile(dist, c(alpha/2, 1-alpha/2))
}
pl <- PL(beta, df.new, alpha=alpha, n.iter=1e4)
#
# Bootstrap a prediction interval.
#
n.iter <- ceiling(10/alpha) # Number of iterations per bootstrap
B <- replicate(5e2, {
  df.boot <- df[sample.int(nrow(df), replace=TRUE), ]
  fit.boot <- glm(cbind(Y,1-Y) ~ ., df.boot, family=binomial)
  PL(coef(fit.boot), df.new, alpha=alpha, n.iter=n.iter)
})
#
# Display the results.
#
par(mfrow=c(1,2))
h <- hist(B, main="Histogram of LPL and UPL", xlab="Sum", border=NA)
hist(B[2,], add=TRUE, col="#80000040", breaks=h$breaks)
hist(B[1,], add=TRUE, col="#00008040", breaks=h$breaks)
abline(v=pl, lwd=2, col=c("Blue", "Red"))

plot(B[1,] + runif(ncol(B), -0.4, 0.4), B[2,] + runif(ncol(B), -0.4, 0.4),
     bty="n",
     xlab="LPL", ylab="UPL", pch=19, col="#00000020", asp=1)
abline(c(mean(B[2,]-B[1,]), 1), col="#a00000", lwd=2)
par(mfrow=c(1,1))