Teoria stojąca za argumentem wag w R przy użyciu lm ()

12

Po roku nauki w szkole, moje rozumienie „ważonych najmniejszych kwadratów” jest następujące: niech , będzie jakaś macierzą projektową, \ boldsymbol \ beta \ in \ mathbb {R} ^ p być wektorem parametrów, \ boldsymbol \ epsilon \ in \ mathbb {R} ^ n być wektorem błędu takim, że \ boldsymbol \ epsilon \ sim \ mathcal {N} (\ mathbf {0}, \ sigma ^ 2 \ mathbf {V}) , gdzie \ mathbf {V} = \ text {diag} (v_1, v_2, \ dots, v_n) i \ sigma ^ 2> 0 . Następnie model \ mathbf {y} = \ mathbf {X} \ boldsymbol \ beta + \ boldsymbol \ epsilonyRnXn×pβRpϵRnϵN(0,σ2V)V=diag(v1,v2,,vn)σ2>0

y=Xβ+ϵ
zgodnie z założeniami nazywany jest modelem „ważonych najmniejszych kwadratów”. Problemem WLS jest znalezienie
argminβ(yXβ)TV1(yXβ).
Załóżmy, że y=[y1yn]T , β=[β1βp]T i
X=[x11x1px21x2pxn1xnp]=[x1Tx2TxnT].
xiTβR1 , więc
yXβ=[y1x1Tβy2x2TβynxnTβ].
To daje
(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} dając w ten sposób
argminβ(yXβ)TV1(yXβ)=argminβi=1nvi1(yixiTβ)2.
β szacuje się za pomocą
β^=(XTV1X)1XTV1y.
Jest to zakres wiedzy, którą znam. Nigdy nie uczono mnie, jak wybrać v1,v2,,vn , chociaż wydaje się, że sądząc po tym , że zwykle Var(ϵ)=diag(σ12,σ22,,σn2), co ma intuicyjny sens. (Daj bardzo zmienne wagi mniej wagi w problemie WLS i daj obserwacjom o mniejszej zmienności większą wagę.)

Szczególnie mnie ciekawi, jak Robsługuje wagi w lm()funkcji, gdy wagi są przypisane jako liczby całkowite. Z użycia ?lm:

Długoterminowe NULLwagi mogą być stosowane w celu wskazania, że różne spostrzeżenia mają różne odchylenia (z wartościami wagi są odwrotnie proporcjonalne do odchylenia); lub równoważnie, gdy elementy wag są dodatnimi liczbami całkowitymi , że każda odpowiedź jest średnią z obserwacji masy jednostkowej (w tym przypadku, gdy obserwacje są równe a dane zostały podsumowane).wiyiwiwiyi

Przeczytałem ten akapit kilka razy i nie ma to dla mnie sensu. Korzystając ze środowiska, które opracowałem powyżej, załóżmy, że mam następujące symulowane wartości:

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  

W jaki sposób opracowano te parametry przy użyciu wyżej opracowanego środowiska? Oto moja próba zrobienia tego ręcznie: zakładając, że , mamy i zrobienie tego w daje (zauważ, że odwracalność nie działa w tym przypadku, więc użyłem uogólnionej odwrotności):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

Nie pasują one do wartości z danych lm()wyjściowych. Co ja robię źle?

Klarnecista
źródło

Odpowiedzi:

4

Matryca powinna wyglądać następująco: „ nie Należy również pamiętać, powinno być , nie .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
źródło
Dziękujemy za usunięcie nieprawidłowej matrycy projektowej, szczególnie! Jestem dość zardzewiały na tym materiale. Czy jako ostatnie pytanie oznacza to, że w założeniach WLS? Var(ϵ)=diag(1/weights)
Klarnecista
Tak, chociaż wagi muszą być tylko proporcjonalne do 1 / wariancji, niekoniecznie równe. Na przykład, jeśli użyjesz weights <- c(50, 85, 75)/2w swoim przykładzie, otrzymasz ten sam wynik.
mark999
3

Aby odpowiedzieć na to bardziej zwięźle, regresja metodą najmniejszych kwadratów ważona za pomocą weightsin Rprzyjmuje następujące założenia: załóżmy, że mamy weights = c(w_1, w_2, ..., w_n). Niech , będzie macierzą projektową , być wektorem parametrów, a być wektorem błędu o średniej i macierzy wariancji , gdzie . Następnie Postępując zgodnie z tymi samymi krokami wyprowadzenia w oryginalnym poście, mamy 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
i szacuje się za pomocą z GLS założenia .β
β^=(XTV1X)1XTV1y
Klarnecista
źródło