Теория весового аргумента в R при использовании lm ()

12

После года, проведенного в аспирантуре, мое понимание «взвешенных наименьших квадратов» таково: пусть , будет некоторой матрицей проектирования, \ boldsymbol \ beta \ in \ mathbb {R} ^ p - вектор параметров, \ boldsymbol \ epsilon \ in \ mathbb {R} ^ n - вектор ошибок, такой что \ boldsymbol \ epsilon \ sim \ mathcal {N} (\ mathbf {0}, \ sigma ^ 2 \ mathbf {V}) , где \ mathbf {V} = \ text {diag} (v_1, v_2, \ dots, v_n) и \ sigma ^ 2> 0 . Тогда модель \ mathbf {y} = \ mathbf {X} \ boldsymbol \ beta + \ boldsymbol \ epsilonyRnXn×pβRpϵRnϵN(0,σ2V)V=diag(v1,v2,,vn)σ2>0

y=Xβ+ϵ
Под предположениями называется модель «взвешенных наименьших квадратов». Проблема WLS заключается в том, чтобы найти
argminβ(yXβ)TV1(yXβ).
Предположим, y=[y1yn]T , β=[β1βp]T и
X=[x11x1px21x2pxn1xnp]=[x1Tx2TxnT].
xiTβR1 , поэтому
yXβ=[y1x1Tβy2x2TβynxnTβ].
Это дает
(yXβ)TV1=[y1x1Tβy2x2TβynxnTβ]diag(v11,v21,,vn1)=[v11(y1x1Tβ)v21(y2x2Tβ)vn1(ynxnTβ)]
v_n ^ {- 1} (y_n- \ mathbf {x} _ {n} ^ {T} \ boldsymbol \ beta) \ end {bmatrix} \ end {align}, что дает
argminβ(yXβ)TV1(yXβ)=argminβi=1nvi1(yixiTβ)2.
β оценивается с использованием
β^=(XTV1X)1XTV1y.
Это та степень знаний, с которыми я знаком. Я никогда не учил , как v1,v2,,vn должен быть выбран, хотя, кажется , что, судя по здесь , что обычно Var(ϵ)=diag(σ12,σ22,,σn2), что имеет интуитивный смысл. (Дайте сильно изменяющимся весам меньший вес в задаче WLS и сделайте наблюдения с меньшей изменчивостью большим весом.)

Что меня особенно интересует, так это то, как Rобрабатываются веса в lm()функции, когда весам назначаются целые числа. От использования ?lm:

Невесовые значения NULLмогут использоваться для указания того, что разные наблюдения имеют разные отклонения (при этом значения в весах обратно пропорциональны отклонениям); или эквивалентно, когда элементы весов являются положительными целыми числами , каждый ответ является средним значением наблюдений за единицу веса (включая случай, когда есть наблюдения равные и данные суммированы).wiyiwiwiyi

Я перечитал этот параграф несколько раз, и для меня это не имеет смысла. Используя структуру, которую я разработал выше, предположим, что у меня есть следующие смоделированные значения:

x <- c(0, 1, 2)
y <- c(0.25, 0.75, 0.85)
weights <- c(50, 85, 75)

lm(y~x, weights = weights)

Call:
lm(formula = y ~ x, weights = weights)

Coefficients:
(Intercept)            x  
     0.3495       0.2834  

Используя структуру, которую я разработал выше, как эти параметры получены? Вот моя попытка сделать это вручную: предполагая, что , у нас есть и выполнение этого в дает (обратите внимание, что обратимость не работает в этом случае, поэтому я использовал обобщенное обратное):V=diag(50,85,75)

[β^0β^1]=([111111]diag(1/50,1/85,1/75)[111111]T)1[111111]Tdiag(1/50,1/85,1/75)[0.250.750.85]
R
X <- matrix(rep(1, times = 6), byrow = T, nrow = 3, ncol = 2)
V_inv <- diag(c(1/50, 1/85, 1/75))
y <- c(0.25, 0.75, 0.85)

library(MASS)
ginv(t(X) %*% V_inv %*% X) %*% t(X) %*% V_inv %*% y

         [,1]
[1,] 0.278913
[2,] 0.278913

Они не соответствуют значениям из lm()вывода. Что я делаю неправильно?

Кларнетист
источник

Ответы:

4

Матрица должна быть не Кроме того, ваш должен быть , а не .X

[101112],
[111111].
V_invdiag(weights)diag(1/weights)
x <- c(0, 1, 2)
y <- c(0.25, 0.75, 0.85)
weights <- c(50, 85, 75)
X <- cbind(1, x)

> solve(t(X) %*% diag(weights) %*% X, t(X) %*% diag(weights) %*% y)
       [,1]
  0.3495122
x 0.2834146
mark999
источник
Спасибо за исправление неправильной матрицы дизайна, особенно! Я довольно ржавый на этом материале. Итак, как последний вопрос, означает ли это, что в предположениях WLS? Var(ϵ)=diag(1/weights)
Кларнетист
Да, хотя веса должны быть пропорциональны только 1 / дисперсии, а не обязательно равны. Например, если вы используете weights <- c(50, 85, 75)/2в своем примере, вы получите тот же результат.
mark999
3

Чтобы ответить на этот вопрос более кратко, регрессия взвешенных наименьших квадратов с использованием weightsin Rделает следующие допущения: предположим, что мы имеем weights = c(w_1, w_2, ..., w_n). Пусть , - матрица дизайна , - вектор параметров, и будет вектором ошибок со средним значением и дисперсионной матрицей , где . Тогда Следуя тем же шагам деривации в исходном сообщении, мы имеем yRnXn×pβRpϵRn0σ2Vσ2>0

V=diag(1/w1,1/w2,,1/wn).
argminβ(yXβ)TV1(yXβ)=argminβi=1n(1/wi)1(yixiTβ)2=argminβi=1nwi(yixiTβ)2
и оценивается с помощью из GLS предположения .β
β^=(XTV1X)1XTV1y
Кларнетист
источник