Как вычисляется перехват в GLM fit?

Я читал код , используемый R для подгонки обобщенной линейной модели (GLM), поскольку исходный код R находится в свободном доступе. Используемый алгоритм называетсяитеративно взвешенными наименьшими квадратами (IRLS), что является достаточно документированным алгоритмом. Для каждой итерации существует вызов функции Фортрана для решения взвешенной задачи наименьших квадратов.

С точки зрения конечного пользователя, например, для логистической регрессии вызов в R выглядит просто вот так: 

y <- rbinom(100, 1, 0.5)
x <- rnorm(100)
glm(y~x, family=binomial)$coefficients

И если вы не хотите использовать перехват , любой из этих вызовов подходит: 

glm(y~x-1, family=binomial)$coefficients
glm(y~x+0, family=binomial)$coefficients
Однако я не могу понять, каким образом формула , т. е. y~x или y~x-1, имеет смысл в коде и понимается как использование перехвата или нет. Я искал часть кода, где столбец единиц был бы привязан к x, но, похоже, его нет.

Спасибо.

PS: насколько я прочитал, булев перехват, который появляется в функции под названием glm.fit, не совпадает с перехватом, о котором я говорю. И то же самое для смещения .

Документация о glm и glm.fit находитсяздесь .

1 2
formula

1 ответ:

Вы, вероятно, искали не в том месте. Обычно model.matrix() вызывается первым в функциях подгонки: 

> D <- data.frame(x1=1:4, x2=4:1)
> model.matrix(~ x1 + x2, D)
  (Intercept) x1 x2
1           1  1  4
2           1  2  3
3           1  3  2
4           1  4  1
attr(,"assign")
[1] 0 1 2
> model.matrix(~ x1 + x2 -1 , D)
  x1 x2
1  1  4
2  2  3
3  3  2
4  4  1
attr(,"assign")
[1] 1 2
>

И это выход model.matrix(), который передается в Fortran. Это относится к lm() и другим модельным слесарям.

Для glm() он отличается и называется только model.frame(), который не добавляет столбец перехвата. Почему это так, имеет отношение к разнице между обобщенными линейными моделями и стандартными линейными моделями и выходит за рамки об этом посте.

6