重回帰分析の係数の求め方

Rのサンプル airquality を使用

R
dat <- na.omit(airquality)
head(dat)
> head(dat)
  Ozone Solar.R Wind Temp Month Day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
7    23     299  8.6   65     5   7
8    19      99 13.8   59     5   8
R
fit <- lm(Solar.R ~ Wind+Temp, data=dat)
summary(fit)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -49.8135    95.5582  -0.521  0.60323   
Wind          0.6478     2.7151   0.239  0.81187   
Temp          2.9331     1.0136   2.894  0.00461 **

下記の係数の推定値を求めてみましょう

R
summary(fit)$coef[1:3] 
> summary(fit)$coef[1:3] 
[1] -49.8135134   0.6478038   2.9331301

\(\bf b=
\left(
\begin{array}{ccc}
-49.8135134 \\
0.6478038 \\
2.9331301
\end{array}
\right)
\)

データフレームを行列に変換して変数行列 \(\bf X\) を作成

(切片用の列を追加するのを忘れずに)

R
dat$b <- c(rep(1, dim(dat)[1]))
dat2 <- dat[ ,c("b",  "Wind", "Temp")]
mat <- as.matrix(dat2)

\(\bf b = (X^TX)^{-1}X^Ty\)より

R
solve(
    t(mat) %*% mat) %*% 
    t(mat) %*% 
    dat$Solar.R
            [,1]
b    -49.8135134
Wind   0.6478038
Temp   2.9331301

コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください

タイトルとURLをコピーしました