Hosmer-Lemeshow模型构建：自动化识别系数变化超20%的变量

Hey there! Let's tackle this coefficient comparison step cleanly—those nested parentheses can get messy fast, especially when you're working with larger models. Here's a straightforward, efficient approach to identify variables where the coefficient changed by more than 20% between your full and reduced models:

Step 1: Extract Coefficients & Align Common Variables #

First, we'll pull the coefficient estimates from both models, then keep only the variables that exist in both models (since your reduced model drops some predictors). This avoids mismatches in row/column positions:

# Extract coefficient estimates (exclude intercept for predictor comparison)
coef_full <- coef(Iris.mod1)[-1]  # Remove intercept term
coef_reduced <- coef(Iris.mod2)[-1]

# Get variables present in both models
common_vars <- intersect(names(coef_full), names(coef_reduced))

# Subset coefficients to only shared predictors
coef_full_common <- coef_full[common_vars]
coef_reduced_common <- coef_reduced[common_vars]

Step 2: Calculate Percentage Change #

Next, compute the absolute percentage change between the reduced and full model coefficients. We use absolute value to capture both increases and decreases beyond 20%:

# Calculate absolute percentage change: |(reduced_coef - full_coef)/full_coef| * 100
pct_change <- abs((coef_reduced_common - coef_full_common) / coef_full_common) * 100

Step 3: Identify Variables Exceeding 20% Change #

Finally, filter for variables where the percentage change crosses your 20% threshold:

# Get variables with >20% coefficient change
vars_over_20 <- names(pct_change[pct_change > 20])

# Print results
cat("Variables with coefficient change >20%:\n")
print(vars_over_20)

Why This Works Better #

• No nested bracket chaos: Breaking the process into separate steps makes the code readable and easy to debug.
• Handles variable mismatches: Using intersect() ensures we only compare predictors that exist in both models, which is critical when your reduced model drops variables.
• Efficient for large datasets: All operations are vectorized (no loops), so it’ll run quickly even with your 93 variables and 1.7M rows—extracting coefficients is a lightweight task compared to fitting the GLMs themselves.

Full Example with Your Iris Data #

Putting it all together with your sample code:

# Prepare binomial iris dataset
iris_sub <- subset(iris, Species != 'virginica')
iris_sub$Species <- as.numeric(ifelse(iris_sub$Species == 'setosa', 1, 0))

# Build full and reduced models
Iris.mod1 <- glm(Species~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width, 
                 data = iris_sub, family = binomial())
Iris.mod2 <- glm(Species~Sepal.Length+Petal.Length, 
                 data = iris_sub, family = binomial())

# Extract and compare coefficients
coef_full <- coef(Iris.mod1)[-1]
coef_reduced <- coef(Iris.mod2)[-1]
common_vars <- intersect(names(coef_full), names(coef_reduced))
coef_full_common <- coef_full[common_vars]
coef_reduced_common <- coef_reduced[common_vars]

pct_change <- abs((coef_reduced_common - coef_full_common) / coef_full_common) * 100
vars_over_20 <- names(pct_change[pct_change > 20])

cat("Variables with >20% coefficient change:\n")
print(vars_over_20)

Notes for Your Large Dataset #

• This code works directly with your actual models—just replace Iris.mod1 and Iris.mod2 with your full dataset model objects.
• To see the exact percentage change values (not just variable names), print the pct_change vector directly: print(pct_change).

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