从 GLM 中删除多个因子预测变量的截距仅适用于模型中的第一个因子

机器算法验证 广义线性模型 分类数据 二项分布 分类编码 截距
2022-04-06 09:11:21

我在 R 中使用 logit 链接函数运行二项式逻辑回归。我的响应是阶乘 [0/1],我有两个多级阶乘预测变量 - 我们称它们为ab在哪里a有 4 个因子水平 (a1,a2,a3,a4)b有 9 个因子水平(b1,b2,,b9). 所以:

mod <- glm(y ~ a + b, family = binomial(logit), data=pretend)
summary(mod)

然后模型输出将显示有关模型的所有信息以及系数。

汇总输出中缺少 a 和 b(a1 和 b1)的因子水平。我知道它受限于模型的“截距”。我已经读过,如果我想删除截距项并查看这些因子水平的估计值,我可以将 -1 添加到模型公式中,即:

mod2 <- glm(y ~ a + b - 1, family=binomial(logit), data=pretend)
summary(mod2)

在新模型 (mod2) 中,截距项消失了,变量 a 的因子水平 a1 在系数列表中给出。但是,变量 b 的因子水平 b1 仍然缺失,并且鉴于不再有截距项,那么我该如何解释该因子水平的优势比呢?

有人可以向我解释如何获得 b1 的系数以及为什么会这样吗?

2个回答

@kjetil b halvorsen 很好地概述了这里的主要思想。让我补充几点。

对于分类变量,抑制截距导致水平意味着编码,而不是默认的参考水平编码。我在这里更详细地解释了这一点: 逻辑回归如何具有因子预测因子而没有截距?

可以使用具有多个分类变量的级别表示编码,但本质上您必须适应完整的交互。在您的情况下,您只想拟合加法模型 ( y~a+b);如前所述,这是您无法做到的。

如果您致力于使用水平手段编码,则该过程相当简单。您首先创建一个新的单个变量作为各种分类变量的所有可能级别的笛卡尔积(组合)。例如,代替原来的两个分类变量(a,具有 4 个级别,和b,具有 9 个级别),您将有一个具有 36 个级别的变量 ( a1b1, a1b2, a1b3, a1b4, a1b5, a1b6, a1b7, a1b8, a1b9, a2b1, ..., a4b9)。然后,您使用新变量使用级别均值编码(即抑制截距)来拟合您的模型:

mod <- glm(y~0+ab, family=binomial(logit), data=pretend)
summary(mod)

再次注意,这等价于glm(y~a*b, ...); 只是输出会以不同的方式呈现。

如您所见,通过删除截距为因子的每个级别获取参数的技巧仅在只有一个因子时才有效。你可以通过计算自由度来理解为什么:让因子aa水平,因素bb水平。然后因素aa1自由度,这意味着指标矩阵a用,用 a 表示的列1在该行所在级别的每一行中,都有排名a1. 同样的因素bb1自由程度。截距有一个自由度。所以模型公式 a+b(这真的是 a+b+1) 有1+a1+b1=a+b1自由程度。去除截距(模型公式 a+b1) 表示相同的模型,只是参数化发生了变化。所以它也必须有a+b1自由程度。1表明不可能有a+b参数,因此其中一个因素仍然必须比水平数少一个参数。

这就解释了你所看到的。但是您仍然可以获得缺失水平的系数b,它应该是零,简单地说。(取决于您使用的对比度)。

为了更清楚地说明这一点,让我们看一个例子。我将使用 R 作为矩阵代数。为了根据因子制作设计矩阵(用 R 的说法“模型矩阵”),我们需要定义对比函数。我使用默认值:

> options("contrasts")
$contrasts
        unordered           ordered 
"contr.treatment"      "contr.poly" 

首先,我们为一个简单的、完全交叉的设计制作了两个因素:

a  <- factor(rep(letters[1:3], 3))
b  <- factor(rep(letters[1:3], each=3))

然后为它们中的每一个设计矩阵:

> X1 <- model.matrix( ~ a-1)
> X2 <- model.matrix( ~b-1)
> X1
  aa ab ac
1  1  0  0
2  0  1  0
3  0  0  1
4  1  0  0
5  0  1  0
6  0  0  1
7  1  0  0
8  0  1  0
9  0  0  1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"

> X2
  ba bb bc
1  1  0  0
2  1  0  0
3  1  0  0
4  0  1  0
5  0  1  0
6  0  1  0
7  0  0  1
8  0  0  1
9  0  0  1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$b
[1] "contr.treatment"

它们中的每一个分别是满级的:

library(MASS)
library(Matrix)  
> Matrix::rankMatrix(X1)
[1] 3
attr(,"method")
[1] "tolNorm2"
attr(,"useGrad")
[1] FALSE
attr(,"tol")
[1] 1.998401e-15
> Matrix::rankMatrix(X2)
[1] 3
attr(,"method")
[1] "tolNorm2"
attr(,"useGrad")
[1] FALSE
attr(,"tol")
[1] 1.998401e-15

但是当组合在一起时,就会出现排名不足,因此它们必须具有“共同”的一个维度:

rankMatrix(cbind(X1, X2))
[1] 5
attr(,"method")
[1] "tolNorm2"
attr(,"useGrad")
[1] FALSE
attr(,"tol")
[1] 1.998401e-15

为了识别公共维度,我们使用Null()package 中的函数MASS,计算空空间:

 Null(t(cbind(X1, X2)))
           [,1]
[1,] -0.4082483
[2,] -0.4082483
[3,] -0.4082483
[4,]  0.4082483
[5,]  0.4082483
[6,]  0.4082483

是的,共同维度是常数向量。