The result of R: lm () is different when using the `weightights` arguments and when using data with manual reversal

Question

To correct heteroskedasticity under error conditions, I run the following weighted least squares regression in R:

#Call: #lm(formula = a ~ q + q2 + b + c, data = mydata, weights = weighting) #Weighted Residuals: # Min 1Q Median 3Q Max #-1.83779 -0.33226 0.02011 0.25135 1.48516 #Coefficients: # Estimate Std. Error t value Pr(>|t|) #(Intercept) -3.939440 0.609991 -6.458 1.62e-09 *** #q 0.175019 0.070101 2.497 0.013696 * #q2 0.048790 0.005613 8.693 8.49e-15 *** #b 0.473891 0.134918 3.512 0.000598 *** #c 0.119551 0.125430 0.953 0.342167 #--- #Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #Residual standard error: 0.5096 on 140 degrees of freedom #Multiple R-squared: 0.9639, Adjusted R-squared: 0.9628 #F-statistic: 933.6 on 4 and 140 DF, p-value: < 2.2e-16 

Where "weighting" is a variable (a function of q ) used to weight observations. q2 just q^2 .

Now, to double-check my results, I manually weigh my variables, creating new weighted variables:

 mydata$a.wls <- mydata$a * mydata$weighting mydata$q.wls <- mydata$q * mydata$weighting mydata$q2.wls <- mydata$q2 * mydata$weighting mydata$b.wls <- mydata$b * mydata$weighting mydata$c.wls <- mydata$c * mydata$weighting 

And run the following regression without a weight parameter and without a constant - since the constant is weighted, column 1 in the original forecast matrix should now be equal to the variable weight:

 Call: lm(formula = a.wls ~ 0 + weighting + q.wls + q2.wls + b.wls + c.wls, data = mydata) #Residuals: # Min 1Q Median 3Q Max #-2.38404 -0.55784 0.01922 0.49838 2.62911 #Coefficients: # Estimate Std. Error t value Pr(>|t|) #weighting -4.125559 0.579093 -7.124 5.05e-11 *** #q.wls 0.217722 0.081851 2.660 0.008726 ** #q2.wls 0.045664 0.006229 7.330 1.67e-11 *** #b.wls 0.466207 0.121429 3.839 0.000186 *** #c.wls 0.133522 0.112641 1.185 0.237876 #--- #Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #Residual standard error: 0.915 on 140 degrees of freedom #Multiple R-squared: 0.9823, Adjusted R-squared: 0.9817 #F-statistic: 1556 on 5 and 140 DF, p-value: < 2.2e-16 

As you can see, the results are similar but not identical. Am I doing something wrong when I manually weigh variables, or does the "weight" parameter do more than just multiply the variables by the weighing vector?

Answer 1

With proper manual weighing you will not see discrepancies.

So the correct way:

 X <- model.matrix(~ q + q2 + b + c, mydata) ## non-weighted model matrix (with intercept) w <- mydata$weighting ## weights rw <- sqrt(w) ## root weights y <- mydata$a ## non-weighted response X_tilde <- rw * X ## weighted model matrix (with intercept) y_tilde <- rw * y ## weighted response ## remember to drop intercept when using formula fit_by_wls <- lm(y ~ X - 1, weights = w) fit_by_ols <- lm(y_tilde ~ X_tilde - 1) 

Although it is generally recommended to use lm.fit and lm.wfit when directly passing to the matrix:

 matfit_by_wls <- lm.wfit(X, y, w) matfit_by_ols <- lm.fit(X_tilde, y_tilde) 

But when using these internal routines lm.fit and lm.wfit , all the input data must be complete cases without NA , otherwise the main procedure C stats:::C_Cdqrls will complain.

If you still want to use the formula interface rather than the matrix, you can do the following:

 ## weight by square root of weights, not weights mydata$root.weighting <- sqrt(mydata$weighting) mydata$a.wls <- mydata$a * mydata$root.weighting mydata$q.wls <- mydata$q * mydata$root.weighting mydata$q2.wls <- mydata$q2 * mydata$root.weighting mydata$b.wls <- mydata$b * mydata$root.weighting mydata$c.wls <- mydata$c * mydata$root.weighting fit_by_wls <- lm(formula = a ~ q + q2 + b + c, data = mydata, weights = weighting) fit_by_ols <- lm(formula = a.wls ~ 0 + root.weighting + q.wls + q2.wls + b.wls + c.wls, data = mydata) 

---

Reproducible example

Let's use the R built-in trees dataset. Use head(trees) to verify this dataset. There is no NA in this dataset. We strive to fit the model:

 Height ~ Girth + Volume 

with arbitrary weights from 1 to 2:

 set.seed(0); w <- runif(nrow(trees), 1, 2) 

We customize this model using weighted regression, either passing weights to lm , or manually converting the data and calling lm without weights:

 X <- model.matrix(~ Girth + Volume, trees) ## non-weighted model matrix (with intercept) rw <- sqrt(w) ## root weights y <- trees$Height ## non-weighted response X_tilde <- rw * X ## weighted model matrix (with intercept) y_tilde <- rw * y ## weighted response fit_by_wls <- lm(y ~ X - 1, weights = w) #Call: #lm(formula = y ~ X - 1, weights = w) #Coefficients: #X(Intercept) XGirth XVolume # 83.2127 -1.8639 0.5843 fit_by_ols <- lm(y_tilde ~ X_tilde - 1) #Call: #lm(formula = y_tilde ~ X_tilde - 1) #Coefficients: #X_tilde(Intercept) X_tildeGirth X_tildeVolume # 83.2127 -1.8639 0.5843 

So, indeed, we see identical results.

Alternatively, we can use lm.fit and lm.wfit :

 matfit_by_wls <- lm.wfit(X, y, w) matfit_by_ols <- lm.fit(X_tilde, y_tilde) 

We can check the coefficients:

 matfit_by_wls$coefficients #(Intercept) Girth Volume # 83.2127455 -1.8639351 0.5843191 matfit_by_ols$coefficients #(Intercept) Girth Volume # 83.2127455 -1.8639351 0.5843191 

Again, the results are the same.