Получение неправильных бета-версий при выполнении регрессии OLS в R
Мой первый вопрос здесь. Эта проблема украла дни из моей жизни. Я знаю, это не так важно, но в то же время: мне нужно знать! Я знаю, что есть много хороших формул для регрессии. Но когда я пытаюсь сделать это, используя старую арифметику, просто для того, чтобы разобраться в этом, я получаю нелепые ответы на бета-версии.
Предполагается, что бета-вектор (X'X)^(-1)X'y (где X - матрица регрессоров, а y - вектор ответов). Я приведу один пример (и то, что он не подходит для OLS, не имеет значения - я просто хочу b:s здесь):
X <- matrix(1:10)
y <- matrix(2:11)
b <- (t(X) %*% X)^(-1) %*% t(X) %*% y
Что дает b = 1,142857, а итоговое значение (lm(y~X)) дает бета = 1 и точку пересечения 1. Я добавляю константу к X, чтобы получить точку пересечения: X <-cbind (X, 1) и результаты I get is b = (2.324675,14.5), который вообще не имеет смысла. Что я здесь не так делаю?
3 ответа
Здесь есть две проблемы. Первое - это проблема обозначений. Степень -1 в формуле фактически указывает на обратную матрицу. Это рассчитывается с solve в R а не с ^-1, что указывает на поэлементные взаимные.
Затем вам нужно создать матрицу дизайна, которая на самом деле содержит перехват.
X <- matrix(1:10)
y <- matrix(2:11)^2
coef(lm(y~X))
#(Intercept) X
# -21 13
X <- cbind(1, X)
solve(t(X) %*% X) %*% t(X) %*% y
# [,1]
#[1,] -21
#[2,] 13
Очевидно, что вы не должны делать эту инверсию матрицы в реальных приложениях (и R lm не делает этого).
Проблема с использованием ^(-1) для обратного. Это не работает так для матриц. solve используется для получения обратной матрицы: https://www.statmethods.net/advstats/matrix.html
# use solve
b <- solve(t(X) %*% X) %*% t(X) %*% y
# fit model without intercept
m <- lm(y~-1+X)
summary(m)
# same coefficients
b
m$coefficients
# with intercept
X2 <- cbind(rep(1, 10), X)
b2 <- solve(t(X2) %*% X2) %*% t(X2) %*% y
m2 <- lm(y~+X)
summary(m2)
b2
m2$coefficients
X <- cbind(1, matrix(1:10))
b<-solve(t(X)%*%X)%*%t(X)%*%y
https://www.rdocumentation.org/packages/Matrix/versions/0.3-26/topics/solve.Matrix