［R］重み付き線形最小二乗法を行う

lm関数を使う。weightsオプションに重みを指定すると、その重みを考慮した計算を行う。

以下の例では、与えられた8組の測定値を使用して、線形最小二乗法により回帰モデルy = b1 + b2 * x + b3 * x ^ 2（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。

> # 計算の準備
> x <- c(1, 2, 3, 5, 6, 7, 8, 9)
> y <- c(1, 4, 10, 25, 32, 49, 64, 75)  # 測定値（図の黒丸）
> w <- rep(1, length(x))
> w[c(3, 5)] <- 100
> print(data.frame(x, y, w))
  x  y   w
1 1  1   1
2 2  4   1
3 3 10 100
4 5 25   1
5 6 32 100
6 7 49   1
7 8 64   1
8 9 75   1
> 
> # 重みを考慮しない計算（図の青実線）
> r1 <- lm(y ~ 1 + x + I(x ^ 2))
> print(r1)
Call:
lm(formula = y ~ 1 + x + I(x^2))
Coefficients:
(Intercept)            x       I(x^2)  
    -0.6731       0.7042       0.8792  
> sm <- summary(r1)
> print(sm)
Call:
lm(formula = y ~ 1 + x + I(x^2))
Residuals:
       1        2        3        4        5        6        7        8 
 0.08963 -0.25225  0.64741  0.17136 -3.20434  1.66150  2.76888 -1.88219 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.6731     2.8545  -0.236  0.82295   
x             0.7042     1.3677   0.515  0.62859   
I(x^2)        0.8792     0.1351   6.507  0.00128 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.225 on 5 degrees of freedom
Multiple R-squared:  0.9954,    Adjusted R-squared:  0.9936 
F-statistic: 546.6 on 2 and 5 DF,  p-value: 1.399e-06
> print(r1)                              # lm関数による計算結果の概要
Call:
lm(formula = y ~ 1 + x + I(x^2))
Coefficients:
(Intercept)            x       I(x^2)  
    -0.6731       0.7042       0.8792  
> print(sm$coef)                         # 推定した回帰係数,その標準誤差,t値,p値
              Estimate Std. Error    t value    Pr(>|t|)
(Intercept) -0.6730552  2.8545448 -0.2357837 0.822953999
x            0.7041995  1.3677286  0.5148679 0.628589629
I(x^2)       0.8792277  0.1351261  6.5067205 0.001280624
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 24.75789
> print(sum((y - fitted(r1)) ^ 2))
[1] 24.75789
> print(sm$df)                           # モデルのランク,自由度,回帰係数の数
[1] 3 5 3
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 4.951578
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 2.225214
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            (Intercept)          x      I(x^2)
(Intercept)   8.1484260 -3.5141134  0.31168334
x            -3.5141134  1.8706816 -0.18061352
I(x^2)        0.3116833 -0.1806135  0.01825906
> 
> # 重みを考慮した計算（図の赤実線）
> r2 <- lm(y ~ 1 + x + I(x ^ 2), weights = w)
> print(r2)
Call:
lm(formula = y ~ 1 + x + I(x^2), weights = w)
Coefficients:
(Intercept)            x       I(x^2)  
      8.024       -2.828        1.143  
> sm <- summary(r2)
> print(sm)
Call:
lm(formula = y ~ 1 + x + I(x^2), weights = w)
Weighted Residuals:
     1      2      3      4      5      6      7      8 
-5.340 -2.942  1.701  2.535 -2.128  4.754  5.433 -0.173 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   8.0244     3.7465   2.142  0.08512 . 
x            -2.8277     1.8037  -1.568  0.17773   
I(x^2)        1.1432     0.1958   5.839  0.00208 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.542 on 5 degrees of freedom
Multiple R-squared:  0.9966,    Adjusted R-squared:  0.9953 
F-statistic: 734.5 on 2 and 5 DF,  p-value: 6.701e-07
> print(r2)                              # lm関数による計算結果の概要
Call:
lm(formula = y ~ 1 + x + I(x^2), weights = w)
Coefficients:
(Intercept)            x       I(x^2)  
      8.024       -2.828        1.143  
> print(sm$coef)                         # 推定した回帰係数,その標準誤差,t値,p値
             Estimate Std. Error   t value    Pr(>|t|)
(Intercept)  8.024356  3.7464592  2.141851 0.085117824
x           -2.827713  1.8037101 -1.567721 0.177732033
I(x^2)       1.143186  0.1957718  5.839381 0.002083779
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 103.1604
> print(sum(w * (y - fitted(r2)) ^ 2))
[1] 103.1604
> print(sm$df)                           # モデルのランク,自由度,回帰係数の数
[1] 3 5 3
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 20.63208
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 4.542255
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            (Intercept)          x      I(x^2)
(Intercept)  14.0359563 -6.6722356  0.70866822
x            -6.6722356  3.2533702 -0.35090188
I(x^2)        0.7086682 -0.3509019  0.03832661
> 
> # 図の描画
> plot(x, fitted(r1), type = "l", col = "blue", xlab = "x", ylab = "y")
> lines(x, fitted(r2), col = "red")
> points(x, y, pch = 20)