分类变量水平的 p 值在泊松回归中代表什么?

机器算法验证 r 假设检验 分类数据 广义线性模型 泊松回归
2022-03-26 17:11:28

我有一个密度不同的泊松模型:

set.seed(1)
df = data.frame(density = 1:5, events = rpois(2000, 1:5))

如果我对此进行回归,我会得到截距约为log(3),这是有道理的,因为 3 是 1:5 的平均值。

glm(events ~ 1, df, family = poisson)  # returns 1.089

但现在假设我想读回密度系数:

glm(events ~ as.factor(density), df, family = poisson)

(为简单起见,我已将其用作density字段的 ID 及其密度。)我希望系数为density[i]log(3-i)因为截距仍为 3。但是,截距似乎并不保持 3 - 在此情况下,截距设置为log(1)在玩这个时,似乎 glm 将截距设置为第一个因素的系数。

现在我开始想知道 glm 回归中的 p 值表示什么。零假设density[i]是否与截距(aka density[1])相同?还是那样density[i] = mean(density)

2个回答

模型系数是基于数据框如何在密度因子水平中生成对比的估计对比。看看这个:

fit <- glm(events ~ as.factor(density), df, family = poisson)

model.matrix(fit)

要查看如何估计这些对比,请将 GLM 作为对象存储在工作区中。这种情况下的截距现在density是等于 1 时的平均对数率(即 1 的对数,即接近 0)。每个参数,例如第一个,标记as.factor(density)2为当密度等于 2 与密度等于 1 时比较事件的对数相对率。

由于中心极限定理,这些模型参数中的每一个都具有已知的极限渐近分布。这方面的理论很好理解,但有点先进。有关结果的说明,请参阅 McCullagh & Nelder,“广义线性模型”。基本上,与线性回归一样,广义线性模型中的自然参数在研究重复下收敛到正态分布。因此,我们可以计算零假设下的极限分布,并计算观察模型系数与实验获得的不一致或更不一致的概率。这与通常的解释非常相似p-从 OLS 模型参数或列联表的简单 Pearson 检验或 t 检验获得的值。

请注意,如果您删除 的as.factor编码density,您将估计一个平均对数相对率,比较density相差 1 个单位的值,并且当 时,截距将被内插为对数事件率density=0,这可能是也可能不是无用的数量。您生成的数据中的对数相对比率不是恒定的,因此模型效应将代表“平均效应”。

例如:

   ## the actual relative rates comparing subsequent density values
relRates <- exp(diff(log(1:5)) 

modelFit <- glm(events ~ density, data=df, family=poisson)
   ## model based relative rate, weighted by random data
exp(coef(modelFit))[2] 
   ## the approximate average log relative rate, converted to relative rate
exp(mean(log(relRates))  

R,默认情况下,使用参考单元编码(我在这里解释:regression-based-for-example-on-days-of-week)。请注意,这称为在 中使用“处理对比R(加州大学洛杉矶分校的统计帮助网站描述了许多类型的编码方案。)

正如@COOLSerdash 所说,您指定的因子水平的 p 值正在测试它们的截距等于参考类别的截距的零假设。您的参考水平的截距只是(Intercept)在输出中调用。在那里测试的零假设是数据来自真实值为0. (我在这里指的是“截距”,而您指的是“均值”,但请记住,如果您没有连续变量,即只有y~1,则截距 = 均值。)

让我们快速看一下泊松模型:

ln(E(Y))=ln(λ)=β0+β1X
λ=exp(β0+β1X)
既然你没有X, 只是ln(λ)=β0. 如果β0=0, 然后λ=exp(0)=1. 因此,这就是针对参考类别进行测试的内容。

我们也可以看这个R

> set.seed(1234)
> y = rpois(10000, 1)
> summary(glm(y~1, family="poisson"))

...
Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.0009995  0.0099950     0.1     0.92
...
> exp(0.0009995)
[1] 1.001
> mean(y)
[1] 1.001

> y2 = rpois(1000, 2)
> summary(glm(y2~1, family="poisson"))

...
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.72900    0.02196   33.19   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
...
> exp(0.72900)
[1] 2.073007
> mean(y2)
[1] 2.073