はじめに

Rで線形回帰を実施すると、下のように出力されます。

> mtcars # Rにもともと入っているデータセット
                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2
> # 線形回帰を実施（目的変数=mpg、説明変数=それ以外）
> result = lm(mpg~cyl+disp+hp+drat+wt+qsec+vs+am+gear+carb, mtcars)
> summary(result) # 線形回帰の結果を表示

Call:
lm(formula = mpg ~ cyl + disp + hp + drat + wt + qsec + vs + 
    am + gear + carb, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4506 -1.6044 -0.1196  1.2193  4.6271 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 12.30337   18.71788   0.657   0.5181  
cyl         -0.11144    1.04502  -0.107   0.9161  
disp         0.01334    0.01786   0.747   0.4635  
hp          -0.02148    0.02177  -0.987   0.3350  
drat         0.78711    1.63537   0.481   0.6353  
wt          -3.71530    1.89441  -1.961   0.0633 .
qsec         0.82104    0.73084   1.123   0.2739  
vs           0.31776    2.10451   0.151   0.8814  
am           2.52023    2.05665   1.225   0.2340  
gear         0.65541    1.49326   0.439   0.6652  
carb        -0.19942    0.82875  -0.241   0.8122  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.65 on 21 degrees of freedom
Multiple R-squared:  0.869,	Adjusted R-squared:  0.8066 
F-statistic: 13.93 on 10 and 21 DF,  p-value: 3.793e-07

このsummary(result)の結果の導出方法を、関数lmを使わないRプログラムと数式で示します。

Coefficients:Estimate

summary(result)のCoefficients:Estimateの項には偏回帰係数（回帰係数）が記載されています。下の方法で求めることができます。

> X = mtcars[, colnames(mtcars)!='mpg']
> Intercept = 1
> X = as.matrix(cbind(Intercept, X))
> y = mtcars[, 'mpg']
> beta = solve(t(X) %*% X) %*% t(X) %*% y 
> beta
                 [,1]
Intercept 12.30337416
cyl       -0.11144048
disp       0.01333524
hp        -0.02148212
drat       0.78711097
wt        -3.71530393
qsec       0.82104075
vs         0.31776281
am         2.52022689
gear       0.65541302
carb      -0.19941925

このbetaは、summary(result)のCoefficients:Estimateと一致します。数式にすると、次式を計算したことになります。

\hat{\boldsymbol{\beta}} = (X^TX)^{-1}X^T\boldsymbol{y}

Residuals

summary(result)のResidualsの項には残差の最小値、1/4分位数、1/2分位数(=中央値)、3/4分位数、最大値が記載されています。下の方法で求めることができます。

> ### Residuals
> epsilon = y - X %*% beta
> quantile(epsilon)
        0%        25%        50%        75%       100% 
-3.4506441 -1.6044018 -0.1196051  1.2192678  4.6270942

このquantile(epsilon)は、summary(result)のResidualsに一致します。epsilonに関して数式にすると、次式を計算したことになります。

\boldsymbol{\epsilon} = \boldsymbol{y} - X\hat{\boldsymbol{\beta}}

Residual standard error

summary(result)のResidual standard errorの項には残差の標準偏差が記載されています。下の方法で求めることができます。

> n = nrow(X)
> p = ncol(X)
> sigma2 = sum(epsilon ** 2) / (n - p)
> sqrt(sigma2)
[1] 2.650197
> n - p
[1] 21

このsqrt(sigma2)は、summary(result)のResidual standard errorに一致します。n-pは、 degrees of freedomに一致します。sigma2に関して数式にすると、次式を計算したことになります。

\hat{\sigma}^2 = \frac{||\boldsymbol{\epsilon}||^2}{n-p}

Coefficients:Std. Error

summary(result)のCoefficients:Std. Errorの項には偏回帰係数の標準偏差が記載されています。下の方法で求めることができます。

> var_beta = diag(sigma2 * solve(t(X) %*% X))
> sqrt(var_beta)
  Intercept         cyl        disp          hp        drat          wt        qsec 
18.71788443  1.04502336  0.01785750  0.02176858  1.63537307  1.89441430  0.73084480 
         vs          am        gear        carb 
 2.10450861  2.05665055  1.49325996  0.82875250

このsqrt(var_beta)は、summary(result)のCoefficients:Std. Errorに一致します。var_betaに関して数式にすると、次式を計算したことになります。

V(\boldsymbol{\hat{\beta}}) = V[(X^TX)^{-1}X^T\boldsymbol{y}] = \hat{\sigma}^2(X^TX)^{-1}

Coefficients:t value

summary(result)のCoefficients:t valueの項には「帰無仮説：偏回帰係数=0 v.s. 対立仮説：偏回帰係数≠0」の検定で使用するt検定量が記載されています。下の方法で求めることができます。

> t_value = beta / sqrt(var_beta)
> t_value
                [,1]
Intercept  0.6573058
cyl       -0.1066392
disp       0.7467585
hp        -0.9868407
drat       0.4813036
wt        -1.9611887
qsec       1.1234133
vs         0.1509915
am         1.2254035
gear       0.4389142
carb      -0.2406258

このt_value は、summary(result)のCoefficients:t valueに一致します。数式にすると、次式を計算したことになります。

t = \frac{\hat{\boldsymbol{\beta}}}{\sqrt{V(\hat{\boldsymbol{\beta}})}}

ただし、ここでの割り算はベクトルの成分同士の割り算を表します。

Coefficients:Pr(>|t|)

summary(result)のCoefficients:Pr(>|t|)の項には上記仮説検定におけるp値が記載されています。この値が小さいほど、偏回帰係数が0であるという仮説が強く否定されます。下の方法で求めることができます。

> pt(abs(t_value), n - p, lower.tail = F) * 2
                [,1]
Intercept 0.51812440
cyl       0.91608738
disp      0.46348865
hp        0.33495531
drat      0.63527790
wt        0.06325215
qsec      0.27394127
vs        0.88142347
am        0.23398971
gear      0.66520643
carb      0.81217871

このpt(abs(t_value), n - p, lower.tail = F) * 2は、summary(result)のCoefficients:Pr(>|t|)valueに一致します。数式にすると、次式を計算したことになります。

P(|\beta| > |t|) = P(\beta < -|t|) +  P(\beta > |t|) = 2P(\beta > |t|)

Multiple R-squared

summary(result)のMultiple R-squaredの項には決定係数が記載されています。1に近いほどモデルがデータによく適合することを示します。下の方法で求めることができます。

> R2 = 1 - sum(epsilon ** 2) / sum((y - mean(y)) ** 2)
> R2
[1] 0.8690158

このR2 は、summary(result)のMultiple R-squaredに一致します。数式にすると、次式を計算したことになります。

R^2 = 1 - \frac{||\boldsymbol{y} - \hat{\boldsymbol{y}}||^2}{||\boldsymbol{y} - \bar{\boldsymbol{y}}||^2}
= 1 - \frac{||\boldsymbol{\epsilon}||^2}{||\boldsymbol{y} - \bar{\boldsymbol{y}}||^2}

Adjusted R-squared

summary(result)のAdjusted R-squaredの項には自由度調整済み決定係数が記載されています。通常の決定係数は、説明変数を増やせば増やすほど1に近づくという性質を持っておりモデルの良し悪しの比較ができません。自由度調整済み決定係数は、自由度による調整を行いモデルの良し悪しの比較をできるようにしたものです。下の方法で求めることができます。

> aR2 = 1 - (sum(epsilon ** 2) / (n - p)) / (sum((y - mean(y)) ** 2) / (n - 1))
> aR2
[1] 0.8066423

このaR2は、summary(result)のAdjusted R-squaredに一致します。数式にすると、次式を計算したことになります。

\tilde{R^2} = 1 - \frac{||\boldsymbol{y} - \hat{\boldsymbol{y}}||^2/(n-p)}{||\boldsymbol{y} - \bar{\boldsymbol{y}}||^2/(n - 1)} = 1 - \frac{||\boldsymbol{\epsilon}||^2/(n-p)}{||\boldsymbol{y} - \bar{\boldsymbol{y}}||^2/(n - 1)}

F-statistic

summary(result)のF-statisticの項には、「帰無仮説：この回帰モデルは有意である v.s. 対立仮説：この回帰モデルは有意ではない」の検定で使用するF検定量が記載されています。下の方法で求めることができます。

F_value = (sum((X %*% beta - mean(y)) ** 2) / (p - 1)) / (sum((y - X %*% beta) ** 2) / (n - p))
F_value
[1] 13.93246

このF_valueは、summary(result)のF-statisticに一致します。数式にすると、次式を計算したことになります。

F = \frac{||\hat{\boldsymbol{y}} - \bar{\boldsymbol{y}}||^2/(p-1)}{||{\boldsymbol{y}} - \hat{\boldsymbol{y}}||^2/(n - p)}

p-value

summary(result)のp-valueの項には上記仮説検定におけるp値が記載されています。この値が小さいほど、回帰モデルが役に立たないものであるという仮説が強く否定されます。下の方法で求めることができます。

> pf(F_value, p - 1, n - p, lower.tail = FALSE)
[1] 3.793152e-07

このpf(F_value, p - 1, n - p, lower.tail = FALSE)は、summary(result)のp-valueに一致します。