如何在R中对Beta回归模型的预测值及置信水平进行自助法计算

Great question—Beta regression's confidence intervals do require a bit more care than linear models, since we're dealing with bounded proportional responses and a non-normal distribution. Let's walk through two reliable methods to get those intervals for your amplitude values (40, 50, 60, 70).

假设你已经拟合了Beta回归模型 #

First, I'll assume you're using the betareg package (the standard for Beta regression in R) and have a model like this:

library(betareg)

# Replace df with your actual data frame, prop with your proportional response variable
model <- betareg(prop ~ amplitude, data = df)

方法1：基于渐近正态性的置信区间（内置方法） #

The predict() function in betareg supports confidence intervals directly, but you need to pay attention to the type argument to ensure you get intervals on the proportional response scale (not the link function scale, like logit).

Here's how to generate predictions and confidence intervals for your target amplitude values:

# Create a data frame with the amplitude values you want to predict
new_data <- data.frame(amplitude = c(40, 50, 60, 70))

# Get predictions + 95% confidence intervals on the response scale
preds <- predict(
  model,
  newdata = new_data,
  interval = "confidence",
  type = "response",  # Critical: returns values on the 0-1 proportion scale
  level = 0.95        # Adjust if you want a different confidence level
)

# Convert the result to a data frame for easier plotting/analysis
pred_df <- as.data.frame(preds)
pred_df$amplitude <- new_data$amplitude

The output pred_df will have columns:

• fit: The predicted proportion for each amplitude
• lwr: Lower bound of the 95% confidence interval
• upr: Upper bound of the 95% confidence interval

This method is fast and works well if you have a reasonably large sample size, as it relies on asymptotic normality of the model coefficients.

方法2：Bootstrap置信区间（更稳健） #

If your sample size is small, asymptotic intervals might not be reliable. Bootstrap resampling is a better choice here—it directly simulates the variability of your predictions by repeatedly fitting the model to resampled versions of your data.

We'll use the boot package for this:

library(boot)

# Define a function that fits the model to resampled data and returns predictions
boot_pred_fn <- function(data, indices) {
  # Resample the data using the provided indices
  resampled_data <- data[indices, ]
  # Fit the Beta regression to the resampled data
  resampled_model <- betareg(prop ~ amplitude, data = resampled_data)
  # Return predictions for your target amplitude values
  predict(resampled_model, newdata = new_data, type = "response")
}

# Run 1000 bootstrap iterations (adjust R if you need more precision)
boot_results <- boot(
  data = df,
  statistic = boot_pred_fn,
  R = 1000
)

# Calculate 95% confidence intervals using quantiles of the bootstrap predictions
boot_ci <- apply(boot_results$t, 2, function(x) quantile(x, c(0.025, 0.975)))

# Organize results into a data frame
boot_pred_df <- data.frame(
  amplitude = new_data$amplitude,
  fit = colMeans(boot_results$t),  # Mean of bootstrap predictions
  lwr = boot_ci[1, ],              # Lower CI bound
  upr = boot_ci[2, ]               # Upper CI bound
)

Bootstrap intervals are more robust to violations of asymptotic assumptions, making them ideal for small datasets or cases where your Beta distribution has extreme shape parameters.

把置信区间加到你的现有图里 #

If you're using ggplot2, adding the intervals is straightforward. Here's an example using the asymptotic predictions:

library(ggplot2)

ggplot(df, aes(x = amplitude, y = prop)) +
  geom_point(size = 2, alpha = 0.6) +  # Your original data points
  geom_line(data = pred_df, aes(y = fit), color = "darkblue", linewidth = 1) +  # Predicted line
  geom_ribbon(
    data = pred_df,
    aes(ymin = lwr, ymax = upr),
    fill = "darkblue",
    alpha = 0.2
  ) +  # Confidence interval ribbon
  labs(x = "Amplitude", y = "Proportion") +
  theme_minimal()

Just replace pred_df with boot_pred_df if you want to use the bootstrap intervals instead.

关键注意事项 #

• Link function choice: If you used a non-default link function (e.g., probit instead of logit), the type = "response" argument still ensures you get intervals on the proportional scale—no extra conversion needed.
• Precision parameter: If your model includes a precision component (e.g., betareg(prop ~ amplitude | amplitude, data = df)), both methods above will automatically account for it in the predictions and intervals.
• Confidence level: Adjust the level argument in predict() (for method 1) or the quantiles in quantile() (for method 2) if you need a different confidence level (e.g., 90% instead of 95%).

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