La théorie derrière l'argument des poids dans R lors de l'utilisation de lm ()

Après une année d'études supérieures, ma compréhension des "moindres carrés pondérés" est la suivante: soit , soit matrice de conception , \ boldsymbol \ beta \ in \ mathbb {R} ^ p soit un vecteur de paramètres, \ boldsymbol \ epsilon \ in \ mathbb {R} ^ n soit un vecteur d'erreur tel que \ boldsymbol \ epsilon \ sim \ mathcal {N} (\ mathbf {0}, \ sigma ^ 2 \ mathbf {V}) , où \ mathbf {V} = \ text {diag} (v_1, v_2, \ dots, v_n) et \ sigma ^ 2> 0 . Ensuite, le modèle \ mathbf {y} = \ mathbf {X} \ boldsymbol \ beta + \ boldsymbol \ epsilon$\mathbf{y} \in \mathbb{R}^n$$\mathbf{X}$$n \times p$$\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\epsilon$$ sous les hypothèses est appelé le modèle des "moindres carrés pondérés". Le problème WLS finit par trouver $$\begin{equation} \arg\min_{\boldsymbol \beta}\left(\mathbf{y}-\mathbf{X}\boldsymbol\beta\right)^{T}\mathbf{V}^{-1}\left(\mathbf{y}-\mathbf{X}\boldsymbol\beta\right)\text{.} \end{equation}$$ Supposons $\mathbf{y} = \begin{bmatrix} y_1 & \dots & y_n\end{bmatrix}^{T}$ , $\boldsymbol\beta = \begin{bmatrix} \beta_1 & \dots & \beta_p\end{bmatrix}^{T}$ , et $$\mathbf{X} = \begin{bmatrix} x_{11} & \cdots & x_{1p} \\ x_{21} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots \\ x_{n1} & \cdots & x_{np} \end{bmatrix} = \begin{bmatrix} \mathbf{x}_{1}^{T} \\ \mathbf{x}_{2}^{T} \\ \vdots \\ \mathbf{x}_{n}^{T} \end{bmatrix}\text{.}$$ $\mathbf{x}_i^{T}\boldsymbol\beta\in \mathbb{R}^1$ , donc $$\mathbf{y}-\mathbf{X}\boldsymbol\beta = \begin{bmatrix} y_1-\mathbf{x}_{1}^{T}\boldsymbol\beta \\ y_2-\mathbf{x}_{2}^{T}\boldsymbol\beta \\ \vdots \\ y_n-\mathbf{x}_{n}^{T}\boldsymbol\beta \end{bmatrix}\text{.}$$ Cela donne $$\begin{align} (\mathbf{y}-\mathbf{X}\boldsymbol\beta)^{T}\mathbf{V}^{-1} &= \begin{bmatrix} y_1-\mathbf{x}_{1}^{T}\boldsymbol\beta &y_2-\mathbf{x}_{2}^{T}\boldsymbol\beta & \cdots & y_n-\mathbf{x}_{n}^{T}\boldsymbol\beta \end{bmatrix}\text{diag}(v_1^{-1}, v_2^{-1}, \dots, v_n^{-1}) \\ &= \begin{bmatrix} v_1^{-1}(y_1-\mathbf{x}_{1}^{T}\boldsymbol\beta) &v_2^{-1}(y_2-\mathbf{x}_{2}^{T}\boldsymbol\beta) & \cdots & v_n^{-1}(y_n-\mathbf{x}_{n}^{T}\boldsymbol\beta) \end{bmatrix} \end{align}$$ v_n ^ {- 1} (y_n- \ mathbf {x} _ {n} ^ {T} \ boldsymbol \ beta) \ end {bmatrix} \ end {align} donnant ainsi $$\begin{equation} \arg\min_{\boldsymbol \beta}\left(\mathbf{y}-\mathbf{X}\boldsymbol\beta\right)^{T}\mathbf{V}^{-1}\left(\mathbf{y}-\mathbf{X}\boldsymbol\beta\right) \end{equation} = \arg\min_{\boldsymbol \beta}\sum_{i=1}^{n}v_i^{-1}(y_i-\mathbf{x}^{T}_i\boldsymbol\beta)^2\text{.}$$ $\boldsymbol\beta$ est estimé à l'aide de $$\hat{\boldsymbol\beta} = (\mathbf{X}^{T}\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{V}^{-1}\mathbf{y}\text{.}$$ C'est l'étendue des connaissances que je connais. On ne m'a jamais appris comment $v_1, v_2, \dots, v_n$ , bien qu'il semble que, à en juger par ici , cela habituellement $\text{Var}(\boldsymbol\epsilon) = \text{diag}(\sigma^2_1, \sigma^2_2, \dots, \sigma^2_n)$, ce qui est intuitif. (Donner des poids très variables moins de poids dans le problème WLS, et donner des observations avec moins de variabilité plus de poids.)

Ce que je suis particulièrement curieux de savoir, c'est comment Rgère les poids dans la lm()fonction lorsque les poids sont affectés à des entiers. De l'utilisation ?lm:

Les non- NULLpoids peuvent être utilisés pour indiquer que différentes observations ont des variances différentes (les valeurs en poids étant inversement proportionnelles aux variances); ou de manière équivalente, lorsque les éléments de poids sont des entiers positifs , que chaque réponse est la moyenne des observations de poids unitaire (y compris le cas où il y a observations égales à et que les données ont été résumées).$w_i$$y_i$$w_i$$w_i$$y_i$

J'ai relu ce paragraphe plusieurs fois, et cela n'a aucun sens pour moi. En utilisant le cadre que j'ai développé ci-dessus, supposons que j'ai les valeurs simulées suivantes:

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  

En utilisant le cadre que j'ai développé ci-dessus, comment ces paramètres sont-ils dérivés? Voici ma tentative de le faire à la main: en supposant que , nous avons et le faire en donne (notez que l'inversibilité ne fonctionne pas dans ce cas, j'ai donc utilisé un inverse généralisé):$\mathbf{V} = \text{diag}(50, 85, 75)$

$$\begin{align}&\begin{bmatrix} \hat\beta_0 \\ \hat\beta_1 \end{bmatrix} = \\ &\left(\begin{bmatrix} 1 & 1\\ 1 & 1\\ 1 & 1 \end{bmatrix}\text{diag}(1/50, 1/85, 1/75)\begin{bmatrix} 1 & 1\\ 1 & 1\\ 1 & 1 \end{bmatrix}^{T} \right)^{-1}\begin{bmatrix} 1 & 1\\ 1 & 1\\ 1 & 1 \end{bmatrix}^{T}\text{diag}(1/50, 1/85, 1/75)\begin{bmatrix} 0.25 \\ 0.75 \\ 0.85 \end{bmatrix} \end{align}$$ 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

Ceux-ci ne correspondent pas aux valeurs de la lm()sortie. Qu'est-ce que je fais mal?

La matrice doit être pas En outre, vous devriez l'être , non .$X$

$$\begin{bmatrix} 1 & 0\\ 1 & 1\\ 1 & 2 \end{bmatrix},$$ $$\begin{bmatrix} 1 & 1\\ 1 & 1\\ 1 & 1 \end{bmatrix}.$$ 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

Pour répondre à cette question de manière plus concise, la régression des moindres carrés pondérés utilisant weightsin Rfait les hypothèses suivantes: supposons que nous l'avons weights = c(w_1, w_2, ..., w_n). Soit , une matrice de conception, un vecteur de paramètres, et être un vecteur d'erreur avec la moyenne et la matrice de variance , où . Ensuite, suivant les mêmes étapes de la dérivation dans le message d'origine, nous avons $\mathbf{y} \in \mathbb{R}^n$$\mathbf{X}$$n \times p$$\boldsymbol\beta\in\mathbb{R}^p$$\boldsymbol\epsilon \in \mathbb{R}^n$$\mathbf{0}$$\sigma^2\mathbf{V}$$\sigma^2 > 0$

$$\mathbf{V} = \text{diag}(1/w_1, 1/w_2, \dots, 1/w_n)\text{.}$$ $$\begin{align} \arg\min_{\boldsymbol \beta}\left(\mathbf{y}-\mathbf{X}\boldsymbol\beta\right)^{T}\mathbf{V}^{-1}\left(\mathbf{y}-\mathbf{X}\boldsymbol\beta\right)&= \arg\min_{\boldsymbol \beta}\sum_{i=1}^{n}(1/w_i)^{-1}(y_i-\mathbf{x}^{T}_i\boldsymbol\beta)^2 \\ &= \arg\min_{\boldsymbol \beta}\sum_{i=1}^{n}w_i(y_i-\mathbf{x}^{T}_i\boldsymbol\beta)^2 \end{align}$$ et est estimé en utilisant du GLS hypothèses .$\boldsymbol\beta$ $$\hat{\boldsymbol\beta} = (\mathbf{X}^{T}\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{V}^{-1}\mathbf{y}$$