在R与Stan中提取Logistic回归模型generated quantities块预测值并与真实观测对比的方法验证及简便实现咨询

I want to extract the predictions from the generated quantities block in my Stan fit and compare them to the true observed values, but I haven't found a straightforward solution. Here's my implementation for a simple logistic regression model:

library(rstan)
library(tidyverse)
library(boot)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
T <- 40
set.seed(123)
x <- sort(runif(T, 0, 10))
alpha <- 1
beta <- 0.2
logit_p <- alpha + beta * x
p <- inv.logit(logit_p)
y <- rbinom(T, 1, p)
model_code <- "
data {
  int<lower=0> N;
  vector[N] x;
  int<lower=0,upper=1> y[N];
}
parameters {
  real alpha;
  real beta;
}
model {
  y ~ bernoulli_logit(alpha + beta * x);
}
generated quantities {
  vector[N] z;
  for (n in 1:N) z[n] = bernoulli_logit_rng(alpha + beta * x[n]);
}"
model_data <- list(
  N = T,
  x = x,
  y = y
)
stan_run <- stan(
  data = model_data,
  model_code = model_code
)
posterior <- rstan::extract(stan_run)
df <- as.data.frame(posterior$z)
df <-df %>% summarise(across(everything(.), ~ ifelse(length(.[which(. == 1)]) > length(.[which(. == 0)]), 1, 0)))

I'm not sure if the above implementation is correct, and I'm wondering if there's a more direct and simpler way to do this.

Great question! Your current approach is technically correct, but there are cleaner, more idiomatic ways to extract and summarize these posterior predictions for comparison with your observed values. Let's break down the improvements:

1. Simplify Your Majority-Vote Calculation #

Your existing code for getting majority-vote predictions works, but it can be streamlined using vectorized operations (faster and more readable for this task):

# Extract posterior draws of your generated z values
posterior_z <- posterior$z

# Calculate the proportion of draws where z == 1 for each observation
pred_probs <- colMeans(posterior_z)

# Get majority-vote predictions (matches your original logic)
pred_majority <- ifelse(pred_probs > 0.5, 1, 0)

This achieves the exact same result as your dplyr pipeline but cuts out unnecessary steps, which is especially helpful for larger datasets.

2. Use Stan's Built-in posterior_predict() (No Custom Generated Quantities Needed) #

If you don't need the custom z variable in your Stan code, you can skip defining it in the generated quantities block entirely. Stan's posterior_predict() function automatically generates posterior predictive draws for your outcome variable based on the model's likelihood:

# Generate posterior predictive draws directly from the fit
posterior_pred <- posterior_predict(stan_run)

# Calculate prediction probabilities and majority votes
pred_probs <- colMeans(posterior_pred)
pred_majority <- ifelse(pred_probs > 0.5, 1, 0)

This is a more straightforward workflow because you don't have to manually code the predictive draw logic in your Stan model—posterior_predict() handles that for you.

3. Compare Predictions to True Observations #

Once you have your predictions, combining them with the true values for analysis is easy. Here's how to create a comparison dataframe and calculate basic metrics:

comparison_df <- tibble(
  x = x,
  true_y = y,
  pred_success_prob = pred_probs,
  pred_majority = pred_majority
)

# Check the first few rows of the comparison
head(comparison_df)

# Calculate overall prediction accuracy
accuracy <- mean(comparison_df$pred_majority == comparison_df$true_y)
cat(paste0("Prediction Accuracy: ", round(accuracy * 100, 2), "%\n"))

4. Visualize the Comparison #

A plot helps you intuit how well the model's predictions align with reality. This plot shows true values, majority-vote predictions, and the model's predicted probability of success:

ggplot(comparison_df, aes(x = x)) +
  geom_point(aes(y = true_y), color = "navy", alpha = 0.7, size = 2) +
  geom_point(aes(y = pred_majority), color = "firebrick", alpha = 0.7, shape = 1, size = 2) +
  geom_line(aes(y = pred_success_prob), color = "forestgreen", linewidth = 1) +
  labs(
    x = "Predictor Variable x",
    y = "Outcome (1 = Success)",
    title = "True Values vs. Posterior Predictions",
    subtitle = "Navy = True Observation, Red = Majority Vote, Green = Predicted Success Probability"
  ) +
  theme_minimal()

Quick Validation of Your Original Code #

For what it's worth, your original method does produce correct majority-vote predictions—the simplified version above just makes the code more efficient and readable. Using posterior probabilities (instead of just majority votes) also gives you more insight into the model's uncertainty for each prediction, which is often more useful than a binary yes/no.

内容的提问来源于stack exchange，提问作者Amin Shn