如何在R中移除负二项回归模型（glm.nb()）的不显著因素？

看起来你已经用glm.nb()搭好了负二项回归模型，但里面有些变量和交互项的统计显著性不足。下面我会结合你的模型输出，一步步教你怎么精简模型，同时保证拟合效果。

首先先回顾下你的模型输出：

summary(nb)
Call: glm.nb(formula = LOS ~ patient + age + obesity + race + miscarriage + primigravida + age*procedure, data = db , link = "log", init.theta = 3.193556108)
Deviance Residuals: 
    Min      1Q  Median      3Q     Max 
-2.5905 -0.7482 -0.3547  0.1226  7.7317 

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)               0.992085   0.089697  11.060  < 2e-16 ***
patient                   0.645335   0.034131  -1.435  0.331313    
raceNon-Bumiputera       -0.154894   0.052383  -2.957  0.003107 ** 
ageMiddle                 0.102854   0.090392   1.138  0.255178    
ageElder                  0.055655   0.099996   0.557  0.577815    
obesityYes               -0.250358   0.145008  -1.727  0.084256 .  
miscarriageMiscarried    -0.060868   0.022096  -2.755  0.005875 ** 
primigravidaYoung         0.015143   0.084857   0.178  0.858366    
primigravidamiddle        0.231431   0.105432   2.195  0.028159 *  
primigravidaElder         0.348212   0.125971   2.764  0.005706 ** 
procedureabortion         0.316578   0.432455   2.467  0.341234    
procedurecsection         0.014367   0.006113  -1.673  0.032131 ** 
ageMiddle:procedureabortion  0.035266   0.041567  -1.451  0.034786 ** 
ageElder:procedurecsection   0.658313   0.412243   2.111  0.324342    
ageMiddle:procedureabortion  0.124248   0.233214   1.353  0.143533    
ageElder:procedurecsection   0.236575   0.034353  -1.235  0.013543 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(3.1936) family taken to be 1)

    Null deviance: 5759.6  on 5381  degrees of freedom
Residual deviance: 4852.7  on 5354  degrees of freedom
AIC: 26956

Number of Fisher Scoring iterations: 1

Theta:  3.1936 
Std. Err.:  0.0899 

2 x log-likelihood:  -26897.5170 

方法一：手动逐步剔除（推荐，更可控） #

手动剔除的好处是你可以结合领域知识判断，不会盲目移除有业务意义的变量。步骤如下：

• 先处理最不显著的单变量：看你的输出，patient（p=0.331）、primigravidaYoung（p=0.858）、ageElder（p=0.578）的p值远高于常用的显著性水平（0.05），可以优先移除。
  比如先移除patient，重新拟合模型：

  nb1 <- glm.nb(formula = LOS ~ age + obesity + race + miscarriage + primigravida + age*procedure, 
                data = db, link = "log", init.theta = nb$theta)
  summary(nb1)
  

• 处理不显著的因子水平：对于age变量，ageMiddle（p=0.255）和ageElder都不显著，你可以考虑把这两个水平合并成一个（比如"Non-Young"），或者如果合并后还是不显著，再考虑移除整个age变量。
  合并因子水平的示例：

  db$age <- factor(db$age, levels = c("Young", "Middle", "Elder"), 
                   labels = c("Young", "Non-Young", "Non-Young"))
  db$age <- droplevels(db$age)
  # 重新拟合模型
  nb2 <- glm.nb(formula = LOS ~ age + obesity + race + miscarriage + primigravida + age*procedure, 
                data = db, link = "log", init.theta = nb1$theta)
  summary(nb2)
  

• 处理交互项：注意回归模型的基本原则：如果保留交互项，必须保留对应的主效应。看你的输出，ageElder:procedurecsection（p=0.324）和重复的ageMiddle:procedureabortion（p=0.143）不显著，而对应的主效应ageElder本身也不显著，所以可以先移除这些交互项，再看主效应是否需要保留。
  移除特定交互项的公式示例：

  nb3 <- glm.nb(formula = LOS ~ age + obesity + race + miscarriage + primigravida + 
                  ageMiddle:procedureabortion + ageElder:procedurecsection, # 只保留显著的交互项
                data = db, link = "log", init.theta = nb2$theta)
  summary(nb3)
  

• 每次移除后检查拟合指标：重点看AIC（值越低模型拟合越好）、残差偏差的变化，确保模型没有变差。

方法二：自动逐步回归（快速但需谨慎） #

如果你想快速得到精简模型，可以用R的step()函数，基于AIC准则自动筛选变量：

# 双向逐步回归（既加变量也删变量）
nb_step <- step(nb, direction = "both")
# 查看最终模型
summary(nb_step)

但要注意几个坑：

• 自动方法可能会移除有业务意义但统计不显著的变量，所以最后一定要结合领域知识审核。
• 它会自动遵守"保留交互项则必须保留主效应"的规则，但还是要检查最终模型的公式是否合理。

关键注意事项 #

• 不要只看p值：统计显著性不等于业务重要性，比如obesityYes（p=0.084）接近显著性水平，如果临床肥胖确实和住院时长（LOS）有关，应该保留。
• 检查多重共线性：如果变量之间高度相关，会导致p值偏高，可以用vif()函数（来自car包）检查方差膨胀因子，VIF>5提示有共线性问题，需要先处理。
• 保持theta值的一致性：重新拟合模型时，尽量用上一个模型的theta值作为初始值，加快收敛速度。

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