二元Logistic回归模型选择基础及数据分析技术问询

Great question! Let’s walk through the entire workflow tailored to your dataset (1158 observations, binary outcome infection, 2 binary predictors, 2 integer-truncated continuous predictors) step by step:

1. Model Selection Basics #

Univariate Screening #

Start by checking individual predictor-outcome relationships to identify variables worth including in your multivariate model:

• For binary predictors (like gender): Use chi-squared tests (or Fisher’s exact test if small cell counts exist) to compare infection rates across groups.
• For continuous predictors (like age): Use t-tests (if normally distributed) or Mann-Whitney U tests to compare distributions between infected/uninfected groups, or fit univariate logistic regressions to get crude odds ratios.

Multivariate Model Building #

Skip over-reliance on automated stepwise selection (e.g., step() in R)—it often leads to overfitting and ignores clinical/theoretical context. Instead:

• Start with a theory-driven model including predictors you have reason to believe link to infection.
• If you want to explore variable selection without overfitting, use penalized regression like LASSO logistic regression.

2. Fitting the Logistic Regression Model #

In R, use the glm() function with the binomial family to fit your model:

# Replace var2 with your second continuous/binary predictor name
model <- glm(infection ~ age + gender + var1 + var2, 
             data = data, 
             family = binomial(link = "logit"))

summary(model)

• Exponentiate coefficients with exp(coef(model)) to get adjusted odds ratios, and check p-values for statistical significance.

3. Assessing Model Goodness of Fit #

a. Deviance Goodness of Fit Test #

This compares your model to a "perfect fit" saturated model:

# Deviance test
pchisq(model$deviance, model$df.residual, lower.tail = FALSE)

A non-significant p-value (p > 0.05) suggests good fit. Note: This test can be overly sensitive with large sample sizes like yours, so pair it with other checks.

b. Hosmer-Lemeshow Test #

Splits observations into 10 groups by predicted probability and compares observed vs. expected infection counts:

library(ResourceSelection)
hoslem.test(data$infection, fitted(model), g = 10)

Again, a non-significant p-value indicates the model fits well. Watch out for uneven group sizes, which can skew results.

c. Residual Analysis #

Check for outliers or poorly fitted observations with deviance residuals:

plot(model, which = 1) # Residuals vs fitted values
plot(model, which = 2) # Q-Q plot of residuals

Residuals with absolute values >3 are potential influential points to investigate.

4. Evaluating Predictive Quality #

a. AUC-ROC Curve #

Measures how well the model distinguishes between infected and uninfected cases. An AUC of 0.5 is no better than chance; >0.7 is acceptable, >0.8 is strong:

library(pROC)
roc_obj <- roc(data$infection, fitted(model))
plot(roc_obj, main = "ROC Curve", print.auc = TRUE)

b. Calibration Curve #

Verifies if predicted probabilities match real-world observed frequencies. Use the rms package:

library(rms)
cal <- calibration(infection ~ fitted(model), data = data)
plot(cal, main = "Calibration Curve")

A curve close to the 45-degree line means your predictions are well-calibrated.

c. Confusion Matrix & Classification Metrics #

Pick a probability cutoff (often 0.5, adjust based on whether you prioritize sensitivity or specificity) to generate metrics:

predicted <- ifelse(fitted(model) > 0.5, 1, 0)
conf_matrix <- table(Actual = data$infection, Predicted = predicted)

# Calculate key metrics
sensitivity <- conf_matrix[2,2]/sum(conf_matrix[2,])
specificity <- conf_matrix[1,1]/sum(conf_matrix[1,])
accuracy <- sum(diag(conf_matrix))/sum(conf_matrix)

cat("Sensitivity:", round(sensitivity, 2), "\n")
cat("Specificity:", round(specificity, 2), "\n")
cat("Accuracy:", round(accuracy, 2), "\n")

Quick Bonus Checks #

• Multicollinearity: Use car::vif(model) to check variance inflation factors—values >5 indicate potential issues.
• Overfitting: Validate with 10-fold cross-validation using the caret package:

library(caret)
train_control <- trainControl(method = "cv", number = 10)
cv_model <- train(infection ~ age + gender + var1 + var2, 
                  data = data, 
                  method = "glm", 
                  family = binomial, 
                  trControl = train_control)
print(cv_model)

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