求GLM及Probit模型中基于多向聚类稳健VCov绘制置信区间的方法

Got it, let's break this down—since predict.glm() doesn't natively let you pass a custom variance-covariance matrix (like the one from multiwaycov::cluster.vcov()), we need to calculate prediction standard errors manually, then build our confidence intervals from there. This approach works for probit models (and all GLMs, really) with multi-way clustered SEs. Here's a step-by-step guide with code:

Step 1: Fit Your Probit Model & Get Clustered VCOV #

First, fit your probit model as usual, then generate the two-way clustered variance-covariance matrix using multiwaycov::cluster.vcov().

# Load required packages
library(multiwaycov)
library(ggplot2)
library(dplyr)

# Example data (replace with your actual dataset)
set.seed(123)
df <- tibble(
  y = rbinom(n = 500, size = 1, prob = 0.3),
  x1 = rnorm(500),
  x2 = rnorm(500),
  firm = sample(1:50, 500, replace = TRUE),
  year = sample(2010:2020, 500, replace = TRUE)
)

# Fit probit model
probit_model <- glm(y ~ x1 + x2, data = df, family = binomial(link = "probit"))

# Generate two-way clustered VCOV (cluster on firm and year)
vcov_clustered <- cluster.vcov(probit_model, cluster = df[, c("firm", "year")])

Step 2: Prepare New Data & Calculate Predictions #

Next, create a dataset for which you want predictions (e.g., varying x1 while holding x2 at its mean). Then compute the linear predictions (on the probit link scale) and their standard errors using the clustered VCOV.

# Create new data for predictions (adjust based on your variables)
new_data <- tibble(
  x1 = seq(min(df$x1), max(df$x1), length.out = 100),
  x2 = mean(df$x2, na.rm = TRUE) # Hold x2 at its mean
)

# Generate design matrix for new data (matches the model's formula)
X_new <- model.matrix(formula(probit_model), data = new_data)

# Calculate linear predictions (link scale: probit)
linear_pred <- X_new %*% coef(probit_model)

# Calculate standard errors for linear predictions
se_link <- sqrt(diag(X_new %*% vcov_clustered %*% t(X_new)))

Step 3: Convert to Response Scale & Build Confidence Intervals #

Probit predictions are on the link scale, so we use pnorm() to convert them to probabilities. We also convert the link-scale confidence intervals to the probability scale.

# Build prediction dataframe with 95% CIs
pred_df <- new_data %>%
  mutate(
    pred_prob = pnorm(linear_pred), # Convert link-scale to probability
    ci_low = pnorm(linear_pred - 1.96 * se_link),
    ci_high = pnorm(linear_pred + 1.96 * se_link)
  )

Step 4: Plot Predictions & Confidence Intervals #

Finally, use ggplot2 to visualize the predicted probabilities and clustered confidence intervals.

ggplot(pred_df, aes(x = x1, y = pred_prob)) +
  geom_line(color = "#2c3e50", linewidth = 1) +
  geom_ribbon(aes(ymin = ci_low, ymax = ci_high), alpha = 0.2, fill = "#3498db") +
  labs(
    x = "Predictor X1",
    y = "Predicted Probability of Y=1",
    title = "Probit Predictions with Two-Way Clustered 95% Confidence Intervals"
  ) +
  theme_minimal()

Generalizing to Other GLMs #

This approach works for any GLM—you just need to swap the inverse link function:

• For logit models: Use plogis() instead of pnorm()
• For Poisson models: Use exp() instead of pnorm()

The core logic stays the same: calculate linear predictions on the link scale, compute SEs using your clustered VCOV, then convert everything to the response scale for plotting.

Key Notes #

• Ensure your new data's design matrix matches the original model (same variables, same coding for factors, etc.)
• Clustered SEs account for dependence in your data, so double-check that you're passing the correct cluster variables to cluster.vcov()
• We calculate CIs on the link scale first because variance is more stable there—converting to the response scale after avoids biased intervals.

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