Die Theorie hinter dem Gewichtungsargument in R bei Verwendung von lm ()

11

Nach einem Jahr in der Graduiertenschule verstehe ich die "gewichteten kleinsten Quadrate" wie folgt: Sei , eine Entwurfsmatrix, \ boldsymbol \ beta \ in \ mathbb {R} ^ p ist ein Parametervektor, \ boldsymbol \ epsilon \ in \ mathbb {R} ^ n ist ein Fehlervektor, so dass \ boldsymbol \ epsilon \ sim \ mathcal {N} (\ mathbf {0}, \ sigma ^ 2 \ mathbf {V}) , wobei \ mathbf {V} = \ text {diag} (v_1, v_2, \ dots, v_n) und \ sigma ^ 2> 0 . Dann das Modell \ mathbf {y} = \ mathbf {X} \ boldsymbol \ beta + \ boldsymbol \ epsilonyRnXn×pβRpϵRnϵN(0,σ2V)V=diag(v1,v2,,vn)σ2>0

y=Xβ+ϵ
unter den Annahmen wird das Modell "gewichtete kleinste Quadrate" genannt. Das WLS-Problem besteht darin,
argminβ(yXβ)TV1(yXβ).
Angenommen, y=[y1yn]T , β=[β1βp]T und
X=[x11x1px21x2pxn1xnp]=[x1Tx2TxnT].
xiTβR1 , also
yXβ=[y1x1Tβy2x2TβynxnTβ].
Dies ergibt
(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} ergibt also
argminβ(yXβ)TV1(yXβ)=argminβi=1nvi1(yixiTβ)2.
β wird unter Verwendung von
β^=(XTV1X)1XTV1y.
Dies ist der Umfang des Wissens, mit dem ich vertraut bin. Mir wurde nie beigebracht, wie v1,v2,,vn sollten, obwohl es nach hier zu urteilen scheint , dass normalerweise Var(ϵ)=diag(σ12,σ22,,σn2), was intuitiv Sinn macht. (Geben Sie stark variablen Gewichten weniger Gewicht im WLS-Problem und geben Sie Beobachtungen mit weniger Variabilität mehr Gewicht.)

Was mich besonders interessiert, ist, wie RGewichte in der lm()Funktion behandelt werden, wenn Gewichte Ganzzahlen zugewiesen werden. Von der Verwendung ?lm:

Nichtgewichte NULLkönnen verwendet werden, um anzuzeigen, dass unterschiedliche Beobachtungen unterschiedliche Abweichungen aufweisen (wobei die Werte in Gewichten umgekehrt proportional zu den Abweichungen sind); oder äquivalent, wenn die Elemente von Gewichten positive ganze Zahlen sind, ist jede Antwort der Mittelwert von Einheitsgewichtsbeobachtungen (einschließlich des Falls, dass es Beobachtungen gibt, die gleich und die Daten zusammengefasst wurden).wiyiwiwiyi

Ich habe diesen Absatz mehrmals gelesen und es macht für mich keinen Sinn. Angenommen, ich habe unter Verwendung des oben entwickelten Frameworks die folgenden simulierten Werte:

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  

Wie werden diese Parameter unter Verwendung des oben entwickelten Frameworks abgeleitet? Hier ist mein Versuch, dies von Hand zu tun: Unter der Annahme von haben wir und dies in ergibt (beachten Sie, dass die Invertierbarkeit in diesem Fall nicht funktioniert, daher habe ich eine verallgemeinerte Umkehrung verwendet):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

Diese stimmen nicht mit den Werten aus der lm()Ausgabe überein . Was mache ich falsch?

Klarinettist
quelle

Antworten:

3

Die Matrix sollte nicht Auch sollten Sie nicht sein .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
quelle
Vielen Dank, dass Sie insbesondere die falsche Designmatrix geklärt haben! Ich bin ziemlich rostig auf diesem Material. dies als letzte Frage, dass in den WLS-Annahmen? Var(ϵ)=diag(1/weights)
Klarinettist
Ja, obwohl die Gewichte nur proportional zu 1 / Varianz sein müssen, nicht unbedingt gleich. Wenn Sie weights <- c(50, 85, 75)/2beispielsweise in Ihrem Beispiel verwenden, erhalten Sie das gleiche Ergebnis.
Mark999
2

Um dies genauer zu beantworten, geht die gewichtete Regression der kleinsten Quadrate unter Verwendung von weightsin von Rfolgenden Annahmen aus: Nehmen wir an, wir haben weights = c(w_1, w_2, ..., w_n). Sei , eine Entwurfsmatrix, ein Parametervektor und ist ein Fehlervektor mit dem Mittelwert und der Varianzmatrix , wobei . Dann ist Nach den gleichen Schritten der Ableitung im ursprünglichen Beitrag haben wir 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
und werden unter Verwendung von geschätzt vom GLS Annahmen .β
β^=(XTV1X)1XTV1y
Klarinettist
quelle