Расчет доверительного интервала множественной линейной регрессии (МНК) вручную

Я пытаюсь понять, как вручную рассчитать доверительный интервал множественной линейной регрессии (OLS). Моя проблема в том, что я не знаю, как рассчитать стандартную ошибку для всех отдельных коэффициентов.

Для регрессии только с одной независимой переменной я следовал следующему руководству: http://stattrek.com/regression/slope-confidence-interval.aspx. В этом руководстве представлена ​​следующая формула:

формула

Как оказалось, формула работает. Однако я не до конца понял формулу. Например, почему (-2) вверху формулы. Чтобы проверить правильность, я написал следующий код r, который уже показывает стандартные ошибки:

x<-1:50
y<-c(x[1:48]+rnorm(48,0,5),rnorm(2,150,5))

QR <- rq(y~x, tau=0.5)
summary(QR, se='boot')

LM<-lm(y~x)

alligator = data.frame(
      lnLength = c(3.78, 3.71, 3.73, 3.78),
      lnWeight = c(4.43, 4.38, 4.42, 4.25)
)

alli.mod1 = lm(lnWeight ~ ., data = alligator)

newdata = data.frame(
      lnLength = c(3.78, 3.71, 3.73, 3.78)
)

y_predicted = predict(alli.mod1, newdata, interval="predict")[,1]
length_mean = mean(alligator$lnLength)
> summary(alli.mod1)

 Call:
 lm(formula = lnWeight ~ ., data = alligator)

 Residuals:
        1        2        3        4 
  0.08526 -0.02368  0.03316 -0.09474 

 Coefficients:
             Estimate Std. Error t value Pr(>|t|)
 (Intercept)   7.5279     5.7561   1.308    0.321
 lnLength     -0.8421     1.5349  -0.549    0.638

 Residual standard error: 0.09462 on 2 degrees of freedom
 Multiple R-squared:  0.1308,   Adjusted R-squared:  -0.3038 
 F-statistic: 0.301 on 1 and 2 DF,  p-value: 0.6383

Затем я вручную вычислил SE, используя следующий код r (согласно приведенной выше формуле):

rss = (alligator$lnWeight[1] - y_predicted[1])^2 + 
      (alligator$lnWeight[2] - y_predicted[2])^2 +
      (alligator$lnWeight[3] - y_predicted[3])^2 + 
      (alligator$lnWeight[4] - y_predicted[4])^2

a = sqrt(rss/(length(y_predicted)-2))

b = sqrt((alligator$lnLength[1] - length_mean)^2 + 
         (alligator$lnLength[2] - length_mean)^2 +
         (alligator$lnLength[3] - length_mean)^2 + 
         (alligator$lnLength[4] - length_mean)^2)

a/b
1.534912

В результате получилось то же значение, что и SE сводки (alli.mod1). Итак, подумал я, может быть, это сработает, когда я попробую с двумя переменными. К сожалению, это привело к неправильному ответу. Как показано в приведенном ниже коде:

alligator = data.frame(
    lnLength = c(3.78, 3.71, 3.73, 3.78),
    lnColor = c(2.43, 2.59, 2.46, 2.22),
    lnWeight = c(4.43, 4.38, 4.42, 4.25)
  )

alli.mod1 = lm(lnWeight ~ ., data = alligator)

newdata = data.frame(
  lnLength = c(3.78, 3.71, 3.73, 3.78),
  lnColor = c(2.43, 2.59, 2.46, 2.22)
)


y_predicted = predict(alli.mod1, newdata, interval="predict")[,1]
length_mean = mean(alligator$lnLength)
color_mean = mean(alligator$lnColor)


rss = (alligator$lnWeight[1] - y_predicted[1])^2 + 
      (alligator$lnWeight[2] - y_predicted[2])^2 +
      (alligator$lnWeight[3] - y_predicted[3])^2 + 
      (alligator$lnWeight[4] - y_predicted[4])^2

a = sqrt(rss/(length(y_predicted)-2))

b = sqrt((alligator$lnColor[1] - color_mean)^2 + 
         (alligator$lnColor[2] - color_mean)^2 +
         (alligator$lnColor[3] - color_mean)^2 + 
         (alligator$lnColor[4] - color_mean)^2)

b1 = sqrt((alligator$lnLength[1] - length_mean)^2 + 
          (alligator$lnLength[2] - length_mean)^2 +
          (alligator$lnLength[3] - length_mean)^2 + 
          (alligator$lnLength[4] - length_mean)^2)

> summary(alli.mod1)
Call:
lm(formula = lnWeight ~ ., data = alligator)

Residuals:
        1         2         3         4 
 0.006725 -0.041534  0.058147 -0.023338 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -3.5746     8.8650  -0.403    0.756
lnLength      1.6569     2.1006   0.789    0.575
lnColor       0.7140     0.4877   1.464    0.381

Residual standard error: 0.07547 on 1 degrees of freedom
Multiple R-squared:  0.7235,    Adjusted R-squared:  0.1705 
F-statistic: 1.308 on 2 and 1 DF,  p-value: 0.5258

> a/b
             1 
     0.2009918 
> a/b1
             1 
     0.8657274

Есть ли общий подход, которому я мог бы следовать, чтобы вычислить стандартную ошибку?


r linear-regression confidence-interval least-squares standard-error
person Danique    schedule 11.11.2017    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Я бы посоветовал немного почитать о OLS, включая множественную регрессию. Есть несколько свободно доступных источников информации; одной из отправных точек может быть веб-сайт STAT 501 штата Пенсильвания. Вы можете найти вывод формулы для стандартных ошибок OLS на слайдах 8 и 9 документа эти слайды из MIT Open Courseware.

По сути, стандартная ошибка - это квадратный корень из дисперсии. Как вы можете видеть на слайдах, которые я связал, формула для ковариационной матрицы коэффициентов составляет 2 (X 'X) -1, где 2 - это дисперсия условия ошибки. Тогда дисперсия каждого j равна j -й диагонали этой матрицы. Поскольку нам неизвестно истинное значение 2, мы оцениваем его, как вы это делали выше - мы извлекаем квадратный корень из суммы квадратов ошибок, деленной на n - p. , где p - количество независимых переменных (включая / + точку пересечения) - в простой регрессии p = 2. Однако, в то время как в случае простой регрессии диагонали (X 'X) -1 можно найти по знаменателю вашей формулы наверху, этого не будет в множественная регрессия; вам нужно будет заняться матричной алгеброй. К счастью, в R это очень просто:

# First we make the example data
alligator = data.frame(
    lnLength = c(3.78, 3.71, 3.73, 3.78),
    lnColor = c(2.43, 2.59, 2.46, 2.22),
    lnWeight = c(4.43, 4.38, 4.42, 4.25)
)
# Then we use lm() for a check on our answers later
alli.mod1 = lm(lnWeight ~ ., data = alligator)
# Find the sum of the squared residuals
rss <- sum(alli.mod1$residuals^2)
# And use that to find the estimate of sigma^2, commonly called S
S <- sqrt(rss / (length(alli.mod1$residuals) - length(alli.mod1$coefficients)))
# Make the X matrix; a column of 1s for the intercept and one for each variable
X <- cbind(rep(1, nrow(alligator)), alligator$lnLength, alligator$lnColor)
# We can multiply matrices using %*%, transpose them with t(),
# and invert them with solve(); so we directly apply the formula above with:
std.errors <- S * sqrt(diag(solve(t(X) %*% X)))
# Now we check our answers:
summary(alli.mod1)$coefficients[ , 2] # the second column is the std. errors
# (Intercept)    lnLength     lnColor 
#   8.8650459   2.1005738   0.4876803 
std.errors
# [1] 8.8650459 2.1005738 0.4876803
person duckmayr    schedule 12.11.2017