モモノキ＆ナノネと一緒にRで重回帰分析をやってみよう
定番サンプルデータセットAirqualityを使った重回帰分析の練習
ナノネ、統計解析用フリーソフト『R』の使い方を練習するよ。今回は定番サンプルデータの『Airquality』を使って重回帰分析をやってみよう。
『Airquality』を使った分析は、演習教材でもよくみかけるけど、どんなデータだったっけ？
Rに標準で入っているデータフレームで、ニューヨーク市の大気環境を観測したデータだよ。1973年5月から9月まで154日間の環境観測値で、サンプルデータに含まれている変数は6つ。Ozone（ppb単位のオゾン濃度）、Solar.R（セントラルパークの太陽放射強度、Langley単位）、Wind（LaGuarddia空港の平均風速、マイル/時単位）、Temp（日最高気温、華氏単位）、それとMonth（月）とDay（日）の合計6カラムで構成されたデータセットだよ。
New York Air Quality Measurements
Daily air quality measurements in New York, May to September 1973.
Daily readings of the following air quality values for May 1, 1973 (a Tuesday) to September 30, 1973.
- Ozone: Mean ozone in parts per billion from 1300 to 1500 hours at Roosevelt Island.
- Solar.R: Solar radiation in Langleys in the frequency band 4000–7700 Angstroms from 0800 to 1200 hours at Central Park.
- Wind: Average wind speed in miles per hour at 0700 and 1000 hours at LaGuardia Airport.
- Temp: Maximum daily temperature in degrees Fahrenheit at LaGuardia Airport.
Ozone（オゾン量）
at Roosevelt Island
Solar.R（日射量）
at Central Park
Wind（風速）
at LaGuardia Airport
Temp（気温）
at LaGuardia Airport
今回はOzone、Solar.R、Wind、Temp、4つの量的変数を分析に使うよ。ニューヨーク市のオゾン濃度を、太陽放射強度、平均風速、最高気温によって説明する重回帰モデルを実践してみよう。
ラジャー。
まずは分析用データの準備。サンプルデータの読み込みとデータチェックから始めるよ。
df_0 <- airquality
str(df_0)
変数MonthとDayは使わないから、除去してOKだよね。
df_1 <- df_0[,1:4] # Month、Dayを除く
print(head(df_1))
summary(df_1)
OzoneとSolar.Rは欠損データがあるみたいだけど、どうする？
今回はNAを含む行データは全部除去しちゃおう。na.omit()関数を使えば簡単に処理できるよ。
df_2 <- na.omit(df_1) # NA含む行データを全て除去
summary(df_2)
print(head(df_2))
NAデータの除去完了。
ナノネ、データの数値をよくみてみると。。。
オゾンと太陽光強度はいいとしても．．．風量と温度は普段馴染みの単位と違うから、観測値が少しイメージしにくいかも。
馴染みの単位に換算してもいいかなぁ？
ナノネがイメージしやすいように値を変換してもいいよ。今回は練習だし、換算しても分析に大きな影響なさそうだし。
df <- df_2
# Wind(風速)の単位変換 マイル/時 -> メートル/秒
df[, "Wind"] <- round(df_2[, "Wind"] * 1609.344 / 3600, 1) # 1マイルは1609.344メートル
# Temp(温度)の単位変換 華氏F -> 摂氏C
df[, "Temp"] <- round((df_2[, "Temp"] - 32) * (5/9), 1)
print(head(df)) # 今回分析に使う単位に変化後のデータ Wind(メートル/秒)、Temp（摂氏℃)
無事に単位変換できた。風速はメートル毎秒(m/s)、温度は摂氏（℃）に変えてみた。
グラフにプロットしてデータのおおまかな傾向を確認してみよう。
plot(df)
Windは負の相関、Tempは正の相関。Solarはわかりにくいけど、弱い正の相関かな？
相関係数の値も確認してみて。cor()関数で表示できるよ。
print(cor(df)) # 相関係数
Ozoneとの相関係数は、Solar.R（0.35）、Wind（-0.61）、Temp（0.70）だ。グラフのイメージとだいたい一緒。
もっとデータの特性を調べたいところだけど、細かなチェックは省略。さっそく重回帰分析を実践してみよう。
重回帰モデルはlm()関数で実行できるよ。単回帰のときと同じだね。引数にモデル式とデータを指定すれば取りあえずOKだよ。
モデル式は、目的変数 ~ 説明変数1 + 説明変数2 + 説明変数3、みたいな感じで表現するよ。
res <- lm(Ozone ~ Solar.R + Wind + Temp, data=df) #線形回帰モデル
# res_1 <- lm(Ozone ~ ., data=df) # 目的変数を除いた残り全変数を説明変数に使う場合は.(ドット)でもOK
今回のモデル式は、Ozone ~ Solar.R + Wind + TempでOKだね。
lm()関数の出力結果を確認してみよう。
summary(res) # モデルの要約
数値がいろいろ出てきたけど．．．。結果はどう見るの？
一番上Callに適用したモデル式が表示されているよ。
今回はOzone ~ Solar.R + Wind + Temp。
2番目Residualsは残差の統計量。目的変数Onzeについて、観測値とモデル予測値の差（ズレ）だよ。値が小さいほど、モデルの精度が良い指標になるよ。
3番目Coefficientsは（偏）回帰係数の結果だよ。Estimate列に各説明変数の偏回帰係数の値が出力されているよ。
- (Intercept) -11.58258
- Solar.R 0.05953
- Wind -7.44281
- Temp 2.97684
だとすると、回帰式は．．．．こんな感じ？？
Ozone = -11.58258 + 0.05953＊Soloar.R -7.44281＊Wind + Temp＊2.97684
うん、そうだね。線形回帰モデルは解析結果を簡単に数式で表現できるから、なにかと応用ができて便利だね。
それとEstimateの右側にある3列は係数の統計量で、標準誤差、T値、P値と並んでいるよ。P値の横を見ると今回の3つ説明変数は、どれも＊付きで検定結果は有意になっているよ。検定については、帰無仮説がどうこう話し出すと長くなるので今回省略、あとで詳しい統計教材を調べてみてね。
ラジャー、統計検定はあとで自習する。
最後の3行2段目にあるR-squaredはモデルの決定係数（0～1.0範囲）で、簡単に言うと説明変数を使ってどれくらい目的変数を説明できているかの指標だよ。今回の場合、6割位（あまり良くないけど）を説明できるモデルと言えるね。
また、その横にあるAdjusted R-squaredは調整済みの決定係数で、モデルに使っている説明変数の個数で補正した決定係数の値だよ。決定係数は説明変数の数を増やせば増やすほど必然的に向上しちゃうから、説明変数の個数にペナルティを与えて計算した値だよ。重回帰分析ではAdjusted R-squaredの方を重視するから覚えておいてね。
重回帰分析ではAdjusted R-squaredを要チェックね。
重回帰モデルが一応作れたから、続いて『多重共線性』のチェックと『変数選択』もやってみよう。
多重共線性？？
多重共線性（(multicolinearity）はマルチコとも呼ばれて、説明変数同土の相関に関係することだよ。説明変数同士に強い相関があると問題があって、モデルの偏回帰係数を正しく計算きないんだ。もし多重共線性が見つかった場合は、強い相関関係にあるどちらか一方の変数をモデルから取り除く必要があるんだ。
同じ意味合いの説明変数を複数いれてはダメってことだね。マルチコさんのチェック方法は？
今回はcarパッケージのvif()関数を使って確認してみよう。
vif()関数の値が10以上だと、多重共線性に問題ありと判断できるよ。
# install.packages("car")
library(car) # carパッケージを使用
print(vif(res)) # VIFの結果が10以上だと多重共線性に問題あり、2以下ベター？？
vif()関数の結果だと、マルチコは今回問題は無さそう。
次は変数選択だよ。
たしか、モデルに入れる説明変数の種類を変えて最適な組み合わせを選ぶヤツだよね。
今回はstep()関数を使って変数選択をやってみよう。引数にモデルを指定するだけでOK。デフォルトは変数を減らす方向の探索動作で、AICの値を基準に最適なモデルを探してくれるよ。
AICは赤池情報量規準（Akaike's Information Criterion）で、モデルの良し悪しを評価する指標のひとつ。AIC値は相対的に値が小さいほど良いモデルといえて、精度と複雑さ（パラメータ数にペナルティを与える）を元に算出されているよ。
# 変数選択 step関数で最もAICが低くなるモデルを選択
res_best <- step(res) # デフォルトは変数を減らす方向でサーチ
res_best
1ステップで終わっちゃったみたいだけど？
そうだね。出力されたAICの値をみてみると、最初のステップでどの変数を抜いてもAICが改善されていないから、説明変数を3つ含んだ元のモデルがベストってことだね。
ちなみに見かたは、＜none＞はそのステップで変数を何も抜かない比較モデル、－付き変数名はその変数を抜いたモデルの結果を示しているよ。
ナノネ、AIC()関数を使って、選択されたモデルのAIC値を確認してみて。
round(AIC(res_best),2) # AIC() -> 998.72 # AiC確認
# step関数で出力されるAIC値は、内部的にextractAIC() 関数を利用
# モデル選択（相対比較でよい）が目的のため、計算量を減らして高速化。
# round(extractAIC(res_best)[2],2) # extractAIC() -> 681.94 step()関数のAIC値と同じ
AICは998.94。step()関数で出力されているAIC値（681.94）と全然違うね。
step()関数は変数選択が目的なので、AICは相対比較さえできれば必要十分。step()関数のAICは計算量を減らして処理速度を優先しているよ。正確なAIC値は、AIC()関数などで確認するようにしてね。
なるほど、覚えておく。
結果をグラフにプロットしてみよう。目的変数Ozoneについて、実際の観測値とモデル予測値を単純にプロットするとこんな感じ。
plot(x=res$fitted.value, y=res$model$Ozone,
main="Airquality Ozone (obs vs lm)",
xlab="lm(formula = Ozone ~ Solar.R + Wind + Temp)",
ylab="observation data",
xlim=c(-50, 150), ylim=c(-50,150))
abline(0, 1, col="red")
abline(h=0, v=0, col="gray")
う～ん．．．フィティング感はチョット微妙。
決定係数の結果からは6割を説明できるモデルらしいから、グラフ的にはこんなもんじゃない？手持ちデータに過剰適合させても予測には役立たないモデルになっちゃうし、チューニングする際もバランスが重要だね。
なにごともバンランス重視だね。
ちなみにこのモデルの欠点として、オゾン量の値がマイナスになっちゃうケースがあるね。たとえば負の相関にある風(Wind）が強く吹いたときとか。lm()関数とは違う非線形モデルを使えば対処できそうだけど、今回は非線形は扱わないのでこのままで。
せっかくだから、説明変数（Solar.R、Wind、Temp）の影響とかもグラフに表現できる？
いっぺんに３つ変数の影響をグラフで表現するのは難しいかも。変数が増えると二次元のプロットではうまくデータを表現できないからね。
値を変化させる変数を絞って（3つの説明変数の内2つは平均値で固定）、小分けのグラフをプロットしてみよう。
# 説明変数の平均値計算
print(Solar.R_mean <- mean(df$Solar.R)) # -> 184.8018 Solar.R平均
print(Wind_mean <- mean(df$Wind)) # -> 4.43964 Wind平均
print(Temp_mean <- mean(df$Temp)) # -> 25.43784 Wind平均
# Solar.R、Wiind、Tempがそれぞれの平均値のとき、モデルによるOzone予測値
pred_mean <- predict(res_best, list(Solar.R=Solar.R_mean, Wind=Wind_mean, Temp=Temp_mean))
round(pred_mean, 1) # Ozone -> 42.1
# summary(df$Ozone)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.0 18.0 31.0 42.1 62.0 168.0
参考までに説明変数が3つとも平均値のとき、Ozone予測値は42.1の結果。
# 各説明変数の影響をグラフにプロット（説明変数３つの内、２つは平均値で固定）
N <- 100
S <- seq(min(df$Solar.R), max(df$Solar.R), length=N)
W <- seq(min(df$Wind), max(df$Wind), length=N)
T <- seq(min(df$Temp), max(df$Temp), length=N)
pred_mean_WT <- predict(res_best, list(Solar.R=S, Wind=rep(Wind_mean, N), Temp=rep(Temp_mean, N)))
pred_mean_RT <- predict(res_best, list(Solar.R=rep(Solar.R_mean, N), Wind=W, Temp=rep(Temp_mean, N)))
pred_mean_RW <- predict(res_best, list(Solar.R=rep(Solar.R_mean, N), Wind=rep(Wind_mean, N), Temp=T))
par(mfrow=c(2,2))
plot(x=S, y=pred_mean_WT, # Ozone vs Solar.R
main="Ozone vs Solar.R\n(Wind,Temp@mean)",
xlab="Solar.R",
ylab="Ozone(pred.)",
ylim=c(0, 100), type="l", lwd=3, col="orange")
plot(x=W, y=pred_mean_RT, # Ozone vs Wind
main="Ozone vs Wind\n(Solar.R,Temp@mean)",
xlab="Wind",
ylab="Ozone(pred.)",
ylim=c(0, 100), type="l", lwd=3, col="skyblue")
plot(x=T, y=pred_mean_RW, # Ozone vs Temp
main="Ozone vs Temp\n(Solar.R,Wind@mean)",
xlab="Temp",
ylab="Ozone(pred.)",
ylim=c(0, 100), type="l", lwd=3, col="orange")
オゾン量は、日射が強くなると少し増加、風が強くふくと減少、気温が上昇すると増加だね。オゾン量の増減に一番影響しているのはどれだろう？（風か温度のどっちかみたいだけど）
各説明変数の偏回帰係数でわかるかも？絶対値的には風（Wind）が一番影響強そう？？
# 偏回帰係数の確認(標準化前)
print(round(res_best$coefficients, 4))
# print(round(coef(res_best), 4)) # coef()を利用する場合
残念ながら単位がバラバラだから、そのままの偏回帰係数では比較できないんだ。でも、標準化した偏回帰係数を求めることで説明変数の影響度を比較することができるようになるよ。
標準偏回帰係数はどうやって求めるの？
標準化（平均0、標準偏差1［分散1］）したデータを使ってモデルを作成すれば、標準偏回帰係数を求めることができるよ。データの標準化はscale()関数で簡単に実行できるよ。
# 標準偏回帰係数を求める
df_std <- data.frame(scale(df)) # 分析用データを全て標準化
# ※1 説明変数のみ標準化する場合（目的変数そのまま）
# df_std2 <- data.frame(Ozone=df[,1], scale(df[,2:4]))
print(head(df_std)) # 平均0、標準偏差1(分散1)に標準化
# print(head(df_std2)) # ※1 説明変数のみ標準化した場合
標準化したデータを準備して、モデルを再作成。
res_std <- lm(Ozone ~ Solar.R + Wind + Temp, data=df_std) # 標準化データを適用
# res_std2 <- lm(Ozone ~ Solar.R + Wind + Temp, data=df_std2) # ※1 説明変数のみ標準化した場合
print(round(res_std$coefficients, 4)) # 標準化偏回帰係数
# print(round(res_std2$coefficients, 4)) # ※1 偏回帰係数の絶対値は違うが、比率は同じ
標準偏回帰係数が求まった。オゾン量の増減は、温度（Temp）の影響が一番強いみたい。
上ではモデル作成用のデータを標準化することで標準偏回帰係数を求めたけど、下のように手動計算で求めることもできるから参考にしてみてね。
# 標準化偏回帰係数を手動で計算する場合（各変数の標準偏差を利用する）
coef_std_manual <-
res_best$coefficients *
c(0, sd(df[,"Solar.R"]), sd(df[,"Wind"]), sd(df[,"Temp"])) / sd(df[, "Ozone"]) # 0は切片用
print(round(coef_std_manual, 4)) # lm()結果と同じ標準化回帰回帰係数の値が得られる
# coef_std_manual2 <-
# res_best$coefficients *
# c(1, sd(df[,"Solar.R"]), sd(df[,"Wind"]), sd(df[,"Temp"]))
# print(round(coef_std_manual2, 4)) # lm()結果と同じ標準化回帰回帰係数の値が得られる
最後におまけでもう一つ練習。交互作用をモデルに組み込んでみよう。交互作用とは要因と要因の掛け合わせ効果で、ある要因の影響で他の要因の効果が異なってくるような作用だよ。
説明変数同士の相乗効果で影響が強まったり弱まったりするヤツだね。
交互作用を利用するのに何か準備が必要？
交互作用を利用するときは、説明変数の値を中心化する必要があるみたいだよ（データから平均値を引いて平均0にするだけ）。
# 説明変数を中心化（平均0に）
df_mean0 <- data.frame(Ozone = df$Ozone,
Solar.R = df$Solar.R - mean(df$Solar.R),
Wind = df$Wind - mean(df$Wind),
Temp = df$Temp - mean(df$Temp))
print(colMeans(df_mean0))
説明変数は全て平均0に中心化した。準備OK。
今回は詳細なチェック作業などは省略で。
交互作用を組み込むには、モデル式で＋を使っていたところを＊に変えるだけでOK。
例）目的変数 ~ 説明変数1 ＊ 説明変数2 ＊ 説明変数3
あとはさっきと同じように処理してみて。
# 試しに交互作用をモデル式に入れてみる
res2 <- lm(Ozone ~ Solar.R * Wind * Temp, data=df_mean0) # *で交互作用の組み込み
res2_best <- step(res2) # 変数選択
交互作用を入れたモデルを作成して変数選択も完了。結果では説明変数同士が：（コロン）繋がれて表示されている箇所が交互作用みたい。step()関数で最終的に選ばれたのは、Ozone ~ Solar.R + Wind + Temp + Solar.R:Temp + Wind:Tempで。交互作用項としてSolar.R:Temp、Wind:Tempだね。
summary(res2_best) # 交互作用ありのモデル要約
AIC(res2_best) # -> 978.7 # 交互作用ありのAIC値
# 比較（交互作用なしのモデル）
# res_ref <- lm(Ozone ~ Solar.R + Wind + Temp, data=df_mean0)
# res_ref_best <- step(res_ref)
# AIC(res_ref) # -> 998.9 # 交互作用なしのAIC値
決定係数とAICの値が少しだけ向上したみたい。もっと詳しく調べる必要がありそうだけど、今回は割愛。
# プロット（相互作用ありのモデルで観測値vs予測値）
plot(x=res2$fitted.value, y=res2$model$Ozone,
main="Airquality Ozone: obs vs lm*",
xlab="lm(formula = Ozone ~ Solar.R + Wind + Temp + Solar.R:Temp + Wind:Temp)",
ylab="observation data",
xlim=c(0, 150), ylim=c(0,150))
abline(0, 1, col="red")
相互作用を入れたモデルでも最初と同じようにプロットしてみた。チョットだけ変わったかな。
キホン的な重回帰分析の練習は以上で。
またね！
0 件のコメント :
コメントを投稿