「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

Rによるデータサイエンス(10/XX)

 Rについて

をベースに学んでいく。
 今回は線形回帰分析(PP.133-146)を扱う。

11. 線形回帰分析

 データ間の関係の規則を関数式で表し、入力に相当する変数X*1を説明変数、出力に相当する変数Yを目的変数として



\begin{aligned}
Y=f(X)
\end{aligned}


という形式の統計モデルおよびそれをデータから求める分析方法を回帰分析という。説明変数が1つならば単回帰分析、複数ならば重回帰分析という。

11.1 単回帰分析

 単回帰モデル



\begin{aligned}
Y=\beta_0+\beta_1 X+\varepsilon, 
\end{aligned}


では最小二乗法により



\begin{aligned}
\hat{\beta}_1&=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})(Y_i-\bar{Y})}}{\displaystyle{\sum_{i=1}^{n}(X_i-\bar{X})^2}}},\\
\hat{\beta}_0&=\bar{Y}-\hat{\beta}_1 \bar{X}
\end{aligned}


と推定するのが代表的である。

11.2 Rパッケージ

 \mathrm{R}では線形回帰分析のために\mathrm{lm}がある。

# lm
lm(formula, data, weights, subset, na.action,...)
   \mathrm{formula} モデル形式を指定。Ex. y~x(y=ax+b), y~x-1(y=ax)
   \mathrm{data} 回帰分析に用いるデータセット(データフレーム)
   \mathrm{weights} 説明変数にウェイトを付ける場合のウェイトベクトル
   \mathrm{subsets} データセットの中の一部分を用いる場合に指定
   \mathrm{na.action} 欠損値の扱い

 出力は変数に格納するとリスト形式で得られ、以下のような構成になっている*2

   \mathrm{coef} 回帰係数 \mathrm{coef(\cdot)}
   \mathrm{fitted} 推定に用いたデータの予測値 \mathrm{fitted(\cdot)}
   \mathrm{deviance} 残差平方和 \mathrm{deviance(\cdot)}
   \mathrm{anova} 回帰係数の分散分析 \mathrm{anova(\cdot)}
   \mathrm{predict} 新たなデータに対する予測値 \mathrm{predict(\cdot)}
   \mathrm{print} \mathrm{summary}よりも簡単な結果要約 \mathrm{print(\cdot)}
   \mathrm{summary} 結果の要約 \mathrm{summary(\cdot)}
   \mathrm{plot} 回帰診断プロット \mathrm{plot(\cdot)}

11.3 Rでの実行例

11.3.1 単回帰の実行例
library("ggplot2")

##########################################
### dataとしてcarsを用いて単回帰を実施 ###
##########################################
data(cars)

# 分布
g <- ggplot(data = cars,mapping = aes(x = speed, y = dist)) + geom_point()+ theme_minimal() + 
  ggtitle("carsの散布図") + 
  theme(plot.title = element_text(hjust = 0.5),legend.position = "bottom",
        legend.title=element_text(size = 7),
        legend.text=element_text(size = 7))
g <-  g + geom_vline(xintercept = 0,linetype = 1,size = 0.5, alpha = 0.7)
g <-  g + geom_hline(yintercept = 0,linetype = 1,size = 0.5, alpha = 0.7)
plot(g)

# 相関の有無を相関係数でも確認
print(cor(x = cars$speed,y = cars$dist,method = "pearson")) # Pearsonの積率相関
print(cor(x = cars$speed,y = cars$dist,method = "kendall")) # Kendallの順位相関


####################
### 単回帰を実施 ###
####################
lm_cars <- lm(formula = formula(dist ~ speed),data = cars)

# 結果の要約
summary(lm_cars)

# Call: 適用したlm()
# Residuals: 残差の特性値
# Coefficients: 回帰結果
#   Intercept:回帰切片, speed: 回帰係数
#   Estimate: 推定値, Std.Error: 標準誤差, 
#   t value:t値, Pr(>|t|): p値
# Residual standard error: 残差の標準誤差, xx degrees of freedom:自由度
# Multiple R-squared: 決定係数
# Adjusted R-squared: 修正済決定係数
# F-statistic: F値(すべての回帰係数が0であるという帰無仮説の検定統計量)

# より簡素な要約
print(lm_cars)

# 残差
residuals(lm_cars)

# 回帰係数
coefficients(lm_cars) # 推定したモデル:(dist) = -17.579 + 3.932 * (speed)

predict(object = lm_cars,newdata = data.frame("speed" = 0:25))

# 回帰直線を追加
g2 <- g + geom_abline(slope = coefficients(lm_cars)[2],intercept = coefficients(lm_cars)[1],colour = "red")
plot(g2)

# 回帰診断図
op <- par(mfrow = c(2,2))
plot(lm_cars)
par(op)


回帰診断図(左上:残差と予測値,右上:残差の正規Q-Qプロット

左下:残差の平方根プロット,右下:残差と影響力プロット)

11.3.2 重回帰の実行例
library("ggplot2")
library("GGally")

#####################
### 線形重回帰分析 ###
######################

# データセットairqualityを利用
data("airquality")

# 散布図行列
ggpairs(airquality[,c(-5,-6)],aes_string(alpha = 0.5))

####################
### 重回帰を実施 ###
####################
lm_airquality <- lm(formula = formula("Ozone ~ Solar.R + Wind + Temp"),data = airquality)

# 要約
summary(lm_airquality)

# より簡素な要約
print(lm_airquality)

# 残差
residuals(lm_airquality)

# 回帰係数
coefficients(lm_airquality) # 推定したモデル:(dist) = -17.579 + 3.932 * (speed)

# 回帰診断図
op <- par(mfrow = c(2,2))
plot(lm_airquality)
par(op)





11.3.3 交互作用モデル

 交互作用を考慮した


\begin{aligned}
(\mathrm{Ozone})=&\beta_0+\beta_1(\mathrm{Solar.R}) +\beta_2(\mathrm{Wind}) +\beta_3(\mathrm{Temp} )+\beta_4(\mathrm{Solar.R})\cdot(\mathrm{Wind})\\&+\beta_5(\mathrm{Solar.R})\cdot(\mathrm{Wind})+\beta_6(\mathrm{Wind})\cdot(\mathrm{Temp})+\varepsilon
\end{aligned}

モデルを考えてみる。
 とはいえ、フルモデルが良いものになるとは限らない(実際の推定結果を見てもp値が少し高い)。そこで\mathrm{AIC}でモデル選択をしてみる。

################################
### 交互作用を追加したモデル ###
################################

lm_airquality2 <- lm(formula = formula("Ozone ~ (Solar.R + Wind + Temp)^2"),data = airquality)

summary(lm_airquality2)

# AICを用いて変数選択を行うことにする
lm_airquality2_AIC <- step(lm_airquality2) # 最もAICの小さい、Solar.R:Windを除いたモデルを選択する

# 最終的なモデル
lm_airquality_f <- lm(formula = formula("Ozone ~ Solar.R + Wind + Temp + Solar.R:Temp + Wind:Temp"),
                      data = airquality)
summary(lm_airquality_f)


モデル選択結果(左:フルモデル,右:選択モデル)

補足 スペック情報

エディション Windows 10 Home
バージョン 20H2
プロセッサ Intel(R) Core(TM) i5-1035G4 CPU @ 1.10GHz 1.50 GHz
実装 RAM 8.00 GB
システムの種類 64 ビット オペレーティング システム、x64 ベース プロセッサ
R バージョン 3.6.3 (2020-02-29)
RStudio バージョン 1.2.5033

*1:ベクトルでも可。

*2:そのリストを以下の関数の引数にすることでも取得できる。

プライバシーポリシー お問い合わせ