# Quantifying the uncertainty of aggregated model predictions

Say I have a binary response variable, Y, that I model using a logistic model with four predictors, A, B, C and D. To make matters concrete, imagine that Y = 1 designates a respondent registering support for something, and 0 an absence of support.

Having estimated the relevant parameters on some sample, S, I then want to see what proportion of 1s (i.e., support) I likely would have seen, had all observations in S taken on a particular value on A. Assume conditions for causal inference are satisfied for A, so that changing its value can be thought of as a (hypothetical) intervention.

So I create a "new" sample, S*, identical, to S, save for each observation taking on the desired value on A. I then use the fitted model to “predict” the probability of Y = 1 for each observation in that sample. Taking the mean of those predictions I get an estimated proportion of support under the relevant intervention.

My question is: how should I quantify the uncertainty of that estimate? I can think of three ways, but am not sure which one (if any) makes sense:

1. Resample from the predicted probabilities of the model and bootstrap a confidence interval for the relevant mean that way.
2. Calculate a confidence interval for the prediction made on each observation (the probability that respondent 1 registers support, etc.) like here, and then create a confidence interval for mean support by taking the means of the upr and lwr values of the individual predictions.
3. Resample from S* to fit a large number of models, generate predictions on each model, and then bootstrap a confidence interval for the relevant mean from these predictions.

Any advice here would be greatly appreciated.

If I am understanding you correctly, you want to create a confidence interval around the proportion of observations that would have resulted in $$Y = 1$$ given $$A = a$$. Your logistic regression already provides us with $$P(Y = 1 | A = a, B = b, C = c, D = d)$$. You are justified in taking the mean as a point estimate of $$P(Y = 1|A = a)$$ by the law of total probability. Because the proportion of observations where $$Y = 1$$ is logically equivalent to the probability that an observation will yield $$Y = 1$$ you could use the normal approximation computed on the predicted probabilities themselves to create a confidence interval.

Here is some R code that would do the trick

# fake data for demonstration
# 20 samples drawn from a uniform distribution between 0 and 1
predicted_probs <- runif(n = 50, min = 0, max = 1)

# estimate of the average probability
global_prob_est <- mean(predicted_probs)

# standard error of the estimate
global_prob_se <- sd(predicted_probs)/length(predicted_probs)

# 90% CI using the 5th and 95th percentiles of a normal distribution.
qnorm(p = (0.05, 0.95), mean = global_prob_est, sd = global_prob_se)

• Thanks so much, David! More straightforward than I thought! Just one question: on your proposed way of proceeding, it seems we capture none of the uncertainty in the model estimates (i.e., the individual predictions), as given by the SEs on my proposed (2). For example, say we have two models generating identical vectors of estimates (your predicted_probs), but with massively different SEs (on one model, they're all low; on the other, they're all high). On your proposal, the CIs for both vectors would be the same. On proposal (2), they would not be. Or am I way off? – kh_one Aug 2 '20 at 11:17
• (1/2) Ah, I see. Yes, my proposal takes the model fit as a plug-in estimator (in the same sense as using a fitted model in a bootstrap). I recognize the challenge as being with respect to the propagation of uncertainty from the standard errors of the logistic regression. This would be straight forward in a Bayesian context because the uncertainty is treated as a random variable versus a non-probabilistic interval. – David Nelson Aug 2 '20 at 13:33
• (2/2) You could use the predict function to get the se on each interval and then compute a confidence interval for each point as you have proposed in (2); however the means of the lower and upper boundaries would not necessarily have the same coverage as the intervals at each point e.g., the means of the lower and upper boundaries of 90% intervals do not necessarily capture the "true parameter" (in this case, a vector of probabilities) in 90% of repeated uses. – David Nelson Aug 2 '20 at 13:40
• Thanks again, David. Yes, that was my thinking with (2). It’s the closest I can think of to what I need, but for the reasons you outline re coverage I’m not sure it’s close enough. That said, it does of course capture something about the uncertainty involved. If you end up thinking of a better way, I’d be happy to hear! – kh_one Aug 2 '20 at 13:58