回归分析：从逻辑回归到负二项回归的探索

### 回归分析：从逻辑回归到负二项回归的探索 #### 1. 逻辑回归的数据考量 逻辑回归模型基于观测值相互独立的假设运行。若观测值为连续年份，该假设便会被违背。偏差残差和其他诊断统计量有助于验证模型，检测诸如连接函数指定错误等问题，可参考 `LogisticDx` 包获取更多信息。 一般而言，逻辑回归模型要求每个预测变量至少对应 10 个事件，这里的事件指响应中出现频率较低类别的观测值。以死刑案例为例，死亡在响应中属较少出现的类别，数据库中有 68 个死刑判决，据此规则，最多允许 6 - 7 个预测变量。 回归系数通过最大似然法估计。由于不存在获取这些最大似然估计的封闭数学形式，R 采用优化算法。有时可能会收到算法未收敛的错误信息，这意味着无法找到合适的解决方案，可能由预测变量过多、事件过少等多种原因导致。 #### 2. 模型拟合优度 评估模型性能的一个指标是整体模型的显著性。相应的似然比检验用于判断给定模型是否比仅含截距的零模型拟合效果显著更好。 要获取检验结果，需查看输出中的残差偏差，它衡量观测和拟合对数似然函数最大值之间的差异。逻辑回归遵循最大似然原则，目标是最小化偏差残差之和，此残差与线性回归中的原始残差类似，线性回归的目标是最小化残差平方和。 零偏差表示仅含截距的模型对响应的预测效果。判断模型时，需将残差偏差与零偏差进行比较，二者差值服从卡方分布，可使用 `lmtest` 包中的检验函数： ```R library(lmtest) lrtest(binom.model.1) ``` 似然比类似于线性回归模型中的 F 检验，它能揭示模型是否显著，但无法说明拟合优度，而在线性回归中，调整后的 R 平方可用于描述拟合优度。逻辑回归模型虽无等效统计量，但已开发出几种伪 R 平方，其值通常在 0 到 1 之间，值越高表示拟合越好。可使用 `BaylorEdPsych` 包中的 `PseudoR2` 函数计算该值： ```R library(BaylorEdPsych) PseudoR2(binom.model.1) ``` 不过要注意，伪 R 平方不能像普通最小二乘（OLS）的 R 平方那样解释，且存在一些已知问题，但能提供大致情况。在此例中，模型的解释能力较低，考虑到在如犯罪判决这样复杂过程的建模中仅使用了两个预测变量，这并不意外。 #### 3. 模型比较 对于嵌套逻辑回归模型，可使用似然比检验（如 `lmtest` 库中的 `lrtest` 函数）比较残差偏差的差异。 ```R lrtest(binom.model.0, binom.model.1) ``` 对于非嵌套模型，可使用 AIC，在逻辑回归模型中，AIC 是标准输出的一部分，无需单独调用 AIC 函数。例如，`binom.model.1` 的 AIC 低于 `binom.model.0`，且差值大于 2，不可忽略。 #### 4. 计数数据模型 逻辑回归仅能处理二元响应。若有计数数据，如特定时间段或地理区域内的死亡或失败次数，可使用泊松或负二项回归。这类数据在处理聚合数据时尤为常见，聚合数据以不同类别中的事件数量形式呈现。 #### 5. 泊松回归 泊松回归模型是广义线性模型，以对数为连接函数，假设响应服从泊松分布。泊松分布只取整数值，适用于计数数据，如固定时间段内发生的事件，前提是事件较为罕见，如每天硬盘故障的数量。 以下以 2013 年硬盘数据集为例，该数据集可从 [https://docs.backblaze.com/public/hard-drive-data/2013_data.zip](https://docs.backblaze.com/public/hard-drive-data/2013_data.zip) 下载，经过处理和简化。原始数据库中的每条记录对应一个硬盘的每日快照，感兴趣的故障变量可为 0（硬盘正常）或 1（硬盘故障前的最后一天）。 可能影响故障出现的潜在预测因素包括： - `model`：硬盘制造商指定的型号编号 - `capacity_bytes`：硬盘容量（字节） - `age_month`：硬盘平均使用月数 - `temperature`：硬盘温度 - `PendingSector`：表示不稳定扇区是否出现的逻辑值（给定硬盘在给定日期等待重新映射） 将原始数据集按这些变量进行聚合，`freq` 变量表示给定类别中的记录数，加载最终清理和聚合后的数据集： ```R dfa <- readRDS('SMART_2013.RData') ``` 查看不同型号硬盘的故障数量： ```R (ct <- xtabs(~model+failure, data=dfa)) ``` 去除没有故障的硬盘型号： ```R dfa <- dfa[dfa$model %in% names(which(rowSums(ct) - ct[, 1] > 0)),] ``` 使用 `ggplot2` 包绘制按型号编号的对数刻度直方图，快速了解故障数量： ```R library(ggplot2) ggplot(rbind(dfa, data.frame(model='All', dfa[, -1] )), aes(failure)) + ylab("log(count)") + geom_histogram(binwidth = 1, drop=TRUE, origin = -0.5) + scale_y_log10() + scale_x_continuous(breaks=c(0:10)) + facet_wrap( ~ model, ncol = 3) + ggtitle("Histograms by manufacturer") + theme_bw() ``` 使用型号编号作为预测变量，拟合泊松回归模型，使用 `glm` 函数，设置 `family = 'poisson'`，默认情况下对预期对数计数建模，使用对数连接。考虑到数据库中每个观测对应不同数量的硬盘，使用 `offset` 函数处理不同组大小： ```R poiss.base <- glm(failure ~ model, offset(log(freq)), family = 'poisson', data = dfa) summary(poiss.base) ``` 解释系数，型号编号是离散预测变量，使用多个虚拟变量表示。默认输出中不显示参考类别，但可随时查询： ```R contrasts(dfa$model, sparse = TRUE) ``` 结果显示参考类别为 `HGST`，虚拟变量将每个型号与 `HGST` 硬盘进行比较。例如，`Hitachi` 的系数为 1.77，意味着 `Hitachi` 硬盘的预期对数计数比 `HGST` 硬盘大约大 1.77。也可计算指数来表示比率： ```R exp(1.7666) ``` 一般来说，X 增加一个单位，Y 乘以 `exp(b)`。 确定模型的显著性，将当前模型与无预测变量的零模型比较，识别残差偏差和零偏差的差异，期望差异足够大，相应的卡方检验显著： ```R lrtest(poiss.base) ``` 模型看似显著，但需确定模型假设是否可能不成立。泊松回归有独立性假设，即事件相互独立，在硬盘故障案例中该假设成立。另一个重要假设是响应服从泊松