predict.glm(, type="terms") 实际上是做什么的？

Question

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

The "terms" option returns a matrix giving the fitted values of each term in the model formula on the linear predictor scale.

因此，如果我的模型具有 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

如何解读这样的结果？

Answer 1

我已经编辑了你的问题，包括 "correct" 获取（原始）模型矩阵、模型系数和你预期的逐项预测的方法。所以你关于如何获得这些的其他问题已经解决了。下面小编就带大家了解一下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.glm(, type="terms") 实际上是做什么的？

What does predict.glm(, type="terms") actually do?

regression

r

predict

lm

glm