对列中的每个值执行泊松回归

Question

我有一个长格式数据框，我正在对其执行泊松回归。

 'data.frame':  20000 obs. of  6 variables:
 $ cal_y  : int  2008 2008 2008 2008 2008 ...
 $ age_y  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ gender : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ cal_m  : int  9 7 8 1 6 11 2 10 3 4 ...
 $ n_outcome: int  276 187 164 352 229 250 332 267 348 291 ...
 $ n_atrisk : int  4645 4645 4645 4645 4645 4645 4645 4645 4645 4645 ...

glm(n_outcome ~ factor(cal_y) + factor(cal_m) + gender + offset(log(n_atrisk)),
data = df, family =poisson)

我想知道每个 age_y 的结果 n_outcome 的曝光系数 cal_y 并且最好能够将这些信息汇总到一个 df 中。

我试过几个误入歧途的 lapply() 和 tapply() 版本。目前，我最好的解决方案是手动执行此操作：

glm(n_outcome ~ factor(cal_y) + factor(cal_m) + gender + offset(log(n_atrisk)),
data = filter(df, age_y >= 0, age_y <1), family =poisson)

但这都很乏味 (range(age_y) = 0 105)，结果不容易合并到一个新的 df 中，我不确定在统计之前对数据进行子集化是否正确执行回归。

欢迎任何指点、评论或帮助。

Answer 1

你可以用 dplyr 和我的 broom 包来做到这一点：

library(dplyr)
library(broom)

results <- df %>%
  group_by(age_y) %>%
  do(tidy(glm(n_outcome ~ factor(cal_y) + factor(cal_m) + gender + offset(log(n_atrisk)),
              data = ., family =poisson)))

这是有效的，因为 group_by 和 do 对每个 age_y 值执行回归，然后 tidy 将每个回归转换为可以重新组合的数据框。

有关更多信息，请参阅 broom and dplyr 插图。

Answer 2

你在上面的评论中描述的是如果你使用效果编码与均值模型会发生什么。两者是等效的，只是表示对参数的不同约束（留一与总和为 0）。但这似乎在歪曲你的想法。如果您使用 age_y 作为分类编码，您将获得适当的输出

通过使用子集回归，您不会在每个模型中包含所有可用信息，这是统计上无效的方法。这也增加了 I 类错误率。您应该在模型中使用所有可用信息。因此，这是正确的规范：

# this is the default way that R handles things, leaving one level out.
glm(n_outcome ~ factor(age_y) + factor(cal_y) + factor(cal_m) + gender + offset(log(n_atrisk)),
data = df, family= poisson(link = "log"))

# In contrast, this will provide an estimate for each level of
# factor(age_y) where the test-statistic is if the coefficient
# is statistically different from zero.
glm(n_outcome ~ 0 + factor(age_y) + factor(cal_y) + factor(cal_m) + gender + offset(log(n_atrisk)),
data = df, family= poisson(link = "log"))

以上是正确的，除非您期望每个 age_y 对 cal_y、gender、and/or n_atrisk 产生不同的效果，您可能需要以估计这些变量与 age_y 的相互作用（仍然可以在单个模型中指定）。

如果你想测试age_y的水平是否彼此不同，你可以测试它们的对比。

也就是说，使用 (... 0 + factor(age_y) ...)——这听起来像你的偏好——为你提供了 age_y 相对于零假设的每个水平的信息，但它不提供之间的统计检验说 age_y:level1 vs age_y:level2。要做这个测试，你需要测试他们的对比。

对列中的每个值执行泊松回归

Perform poisson regression for each value in column

r

poisson