逻辑回归拟合度评估咨询：二项式数据分组方案探讨

Hey there! Great call remembering that raw deviance tests aren't ideal for individual-level binomial logistic regression data—grouping is absolutely the right approach to evaluate fit here. Since you mentioned drafting two grouping methods, let’s walk through the most common (and practical) approaches you’re likely considering, along with their pros, cons, and how to put them into action:

This method slices your continuous MeanTemp variable into intervals (like equal-width or quantile-based bins), groups all observations within each interval, then calculates the actual number of snow days and total observations per group. You can then use this aggregated data to run goodness-of-fit tests (like Pearson’s chi-squared or deviance tests).

• Pros: Super intuitive, directly tied to your key predictor. It lets you easily spot where the model over- or under-predicts snow events across temperature ranges.
• Key considerations:
  • Make sure each bin has enough observations (aim for at least 5 successes/failures per group to avoid unstable expected frequencies). If your sample size is small, don’t split into too many bins.
  • Bin choice matters: Equal-width bins might leave you with sparse groups at extreme temperatures, while quantile-based bins (e.g., deciles, quartiles) balance sample sizes but can artificially split natural temperature patterns.
• Quick R example:

# Split MeanTemp into decile bins
data$temp_bin <- cut(data$MeanTemp, 
                     breaks = quantile(data$MeanTemp, probs = seq(0, 1, 0.1)),
                     include.lowest = TRUE)

# Aggregate actual snow counts and total observations per bin
grouped_data <- aggregate(SnowBinary ~ temp_bin, data = data, 
                          FUN = function(x) c(total = length(x), snow_days = sum(x)))
grouped_data <- do.call(data.frame, grouped_data)

# Fit original model and get predicted probabilities per bin
model <- glm(SnowBinary ~ MeanTemp, data = data, family = binomial)
grouped_data$pred_prob <- predict(model, newdata = grouped_data, type = "response")
grouped_data$expected_snow <- grouped_data$SnowBinary.total * grouped_data$pred_prob

# Run Pearson chi-squared test
chisq.test(x = grouped_data$SnowBinary.snow_days, 
           p = grouped_data$pred_prob, 
           n = grouped_data$SnowBinary.total)

Here, you first fit your initial logistic regression model to get a predicted snow probability for each observation. Then you group these predictions into intervals (again, often quantile-based), and compare actual vs. expected snow counts per group.

• Pros: Focuses directly on the model’s predictive performance, avoiding the subjectivity of choosing how to bin your temperature variable. It’s great for checking if the model’s predicted probabilities align with real outcomes across the entire probability spectrum.
• Key considerations:
  • Results depend on your initial model’s quality—if the initial fit is poor, the grouped test might highlight that, but it’s still useful for diagnosing issues.
  • Watch for sparse groups in extreme probability ranges (e.g., very cold days where snow is almost certain, or warm days where it’s nearly impossible).
• Quick R example:

# Fit initial model
model <- glm(SnowBinary ~ MeanTemp, data = data, family = binomial)
data$pred_prob <- predict(model, type = "response")

# Split predicted probabilities into decile bins
data$prob_bin <- cut(data$pred_prob, 
                     breaks = quantile(data$pred_prob, probs = seq(0, 1, 0.1)),
                     include.lowest = TRUE)

# Aggregate actual counts, total observations, and average predicted probability per bin
grouped_data <- aggregate(cbind(SnowBinary, pred_prob) ~ prob_bin, data = data,
                          FUN = function(x) c(total = length(x), snow_days = sum(x), avg_pred = mean(x)))
grouped_data <- do.call(data.frame, grouped_data)
grouped_data$expected_snow <- grouped_data$SnowBinary.total * grouped_data$pred_prob.avg_pred

# Run deviance test via a grouped binomial model
grouped_model <- glm(cbind(SnowBinary.snow_days, SnowBinary.total - SnowBinary.snow_days) ~ pred_prob.avg_pred,
                     data = grouped_data, family = binomial)
# The deviance statistic here is your goodness-of-fit measure
cat("Deviance statistic:", summary(grouped_model)$deviance, "\n")

• Avoid over-grouping: Too many bins lead to small sample sizes per group, making test results unreliable.
• Pair these tests with visual checks: Plot the actual snow proportion vs. predicted probability per group (a calibration curve) to get a more intuitive sense of model fit.
• If you have a large sample, the Hosmer-Lemeshow test is a pre-packaged version of the predicted-probability grouping method—many stats libraries (like ResourceSelection in R) have built-in functions for this.

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