predict.glm(, type="terms") 实际上做了啥？

Posted 2023-02-16

【中文标题】predict.glm(, type="terms") 实际上做了啥？

【英文标题】：What does predict.glm(, type="terms") actually do?predict.glm(, type="terms") 实际上做了什么？

【发布时间】：2016-10-24 03:38:38

【问题描述】：

我对 R 中 predict.glm 函数的工作方式感到困惑。 根据帮助，

“术语”选项返回一个矩阵，该矩阵给出模型公式中每个术语在线性预测变量尺度上的拟合值。

因此，如果我的模型的形式为 f(y) = X*beta，那么命令

predict(model, X, type='terms')

预计会产生相同的矩阵 X，乘以 beta 元素。例如，如果我训练以下模型

test.data = data.frame(y = c(0,0,0,1,1,1,1,1,1), x=c(1,2,3,1,2,2,3,3,3))
model = glm(y~(x==1)+(x==2), family = 'binomial', data = test.data)

得到的系数是

beta <- model$coef

设计矩阵是

X <- model.matrix(y~(x==1)+(x==2), data = test.data)

  (Intercept) x == 1TRUE x == 2TRUE
1           1          1          0
2           1          0          1
3           1          0          0
4           1          1          0
5           1          0          1
6           1          0          1
7           1          0          0
8           1          0          0
9           1          0          0

然后乘以它应该看起来像的系数

pred1 <- t(beta * t(X))

  (Intercept) x == 1TRUE x == 2TRUE
1    1.098612  -1.098612  0.0000000
2    1.098612   0.000000 -0.4054651
3    1.098612   0.000000  0.0000000
4    1.098612  -1.098612  0.0000000
5    1.098612   0.000000 -0.4054651
6    1.098612   0.000000 -0.4054651
7    1.098612   0.000000  0.0000000
8    1.098612   0.000000  0.0000000
9    1.098612   0.000000  0.0000000

但是，predict.glm 生成的实际矩阵似乎与此无关。以下代码

pred2 <- predict(model, test.data, type = 'terms')

      x == 1     x == 2
1 -0.8544762  0.1351550
2  0.2441361 -0.2703101
3  0.2441361  0.1351550
4 -0.8544762  0.1351550
5  0.2441361 -0.2703101
6  0.2441361 -0.2703101
7  0.2441361  0.1351550
8  0.2441361  0.1351550
9  0.2441361  0.1351550
attr(,"constant")
[1] 0.7193212

如何解释这样的结果？

【问题讨论】：

似乎在预测terms 时，预测使用了不同的对比，但内置的似乎都不起作用。另外，确认all.equal(rowSums(predict(model, test.data, type = 'terms')) + attributes(predict(model, test.data, type = 'terms'))$constant, predict(model, test.data))

哲元，别这么慌；）

【参考方案1】：

我已经编辑了您的问题，包括获取（原始）模型矩阵、模型系数和您预期的术语预测的“正确”方式。所以你关于如何获得这些的另一个问题已经解决了。下面我将帮助你理解predict.glm()。

---

predict.glm()（实际上是predict.lm()）在进行逐项预测时为每个模型项应用了居中约束。

最初，你有一个模型矩阵

X <- model.matrix(y~(x==1)+(x==2), data = test.data)

但它是居中的，通过删除列意味着：

avx <- colMeans(X)
X1 <- sweep(X, 2L, avx)

> avx
(Intercept)  x == 1TRUE  x == 2TRUE 
  1.0000000   0.2222222   0.3333333 

> X1
  (Intercept) x == 1TRUE x == 2TRUE
1           0  0.7777778 -0.3333333
2           0 -0.2222222  0.6666667
3           0 -0.2222222 -0.3333333
4           0  0.7777778 -0.3333333
5           0 -0.2222222  0.6666667
6           0 -0.2222222  0.6666667
7           0 -0.2222222 -0.3333333
8           0 -0.2222222 -0.3333333
9           0 -0.2222222 -0.3333333

然后使用这个居中的模型矩阵完成逐项计算：

t(beta*t(X1))

  (Intercept) x == 1TRUE x == 2TRUE
1           0 -0.8544762  0.1351550
2           0  0.2441361 -0.2703101
3           0  0.2441361  0.1351550
4           0 -0.8544762  0.1351550
5           0  0.2441361 -0.2703101
6           0  0.2441361 -0.2703101
7           0  0.2441361  0.1351550
8           0  0.2441361  0.1351550
9           0  0.2441361  0.1351550

居中后，不同的项垂直移动到零均值。结果，截距将变为 0。不用担心，通过聚合所有模型项的变化来计算新的截距：

intercept <- as.numeric(crossprod(avx, beta))
# [1] 0.7193212

现在您应该已经看到 predict.glm(, type = "terms") 为您提供了什么。

【讨论】：

有没有办法使用 predict() 获得未居中的术语？