Rを使っての重回帰分析をggplotで図示しながらやる
手を動かしながら学ぶ ビジネスに活かすデータマイニングの勉強第2回
- 作者: 尾崎隆
- 出版社/メーカー: 技術評論社
- 発売日: 2014/08/22
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
ggplotと格闘しつつ，きれいな図を作りながら勉強したのでその記録を残す．
今回は4章の重回帰分析を行った．
回帰分析
目的変数を説明変数で説明できるモデルを推定すること
例： アイスクリームの売上(目的変数) と 気温(説明変数)
y_i = \beta_0 + \beta{1i} x{1i} + \beta{2i} x{2i} + \cdots + \beta{ni} x{ni} + \epsilon_i
y : 目的変数
x : 説明変数
\beta : 偏回帰係数(説明変数の重み)
n: n番目の説明変数
i: i番目に測定されたデータ
\betaを求めることが重回帰分析の目的
\betaの正負，大きさで説明変数と目的変数の関係がわかる
オゾン濃度とその他要因の関係を調べる
元データ，月日は使わないので除外する．
> head(airquality)
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
5 NA NA 14.3 56 5 5
6 28 NA 14.9 66 5 6
> airq <- airquality[,1:4]
> head(airq)
Ozone Solar.R Wind Temp
1 41 190 7.4 67
2 36 118 8.0 72
3 12 149 12.6 74
4 18 313 11.5 62
5 NA NA 14.3 56
6 28 NA 14.9 66
ggplotで図示するためにmeltする
> mairq <- melt(airq)
> head(mairq)
variable value
1 Ozone 41
2 Ozone 36
3 Ozone 12
4 Ozone 18
5 Ozone NA
6 Ozone 28
んでこれだけではx軸の情報が表にないので追加する
> mairq$num <- c(1:153)
> head(mairq)
variable value num
1 Ozone 41 1
2 Ozone 36 2
3 Ozone 12 3
4 Ozone 18 4
5 Ozone NA 5
6 Ozone 28 6
.
.
154 Solar.R NA 1
156 Solar.R 28 2
.
.
> p <- ggplot(mairq, aes(x=num, y=value, fill=variable, color=variable)) + geom_line()
> p
できた図がこれ
他のバージョンで図示する．
> p + facet_grid(variable ~ .)
> p + facet_grid(variable ~ ., scales = "free_y")
各項目ごとに値がバラけている場合には分けてみると見やすい．
次に回帰分析lm()をやる．
> airq.lm <- lm(Ozone~., airq)
Ozoneを目的変数に説明変数は.(dot)
この意味を調べるとr - Meaning of ~. (tilde dot) argument? - Stack Overflowより，
The dot means "any columns from data that are otherwise not used". とデータフレームの中で他に使っていないカラムとのことである．
結果をそのまま，summary関数(詳細やまとめを表示する)を通してみると
> airq.lm
Call:
lm(formula = Ozone ~ ., data = airq)
Coefficients:
(Intercept) Solar.R Wind Temp
-64.34208 0.05982 -3.33359 1.65209
> summary(airq.lm)
Call:
lm(formula = Ozone ~ ., data = airq)
Residuals:
Min 1Q Median 3Q Max
-40.485 -14.219 -3.551 10.097 95.619
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -64.34208 23.05472 -2.791 0.00623 **
Solar.R 0.05982 0.02319 2.580 0.01124 *
Wind -3.33359 0.65441 -5.094 1.52e-06 ***
Temp 1.65209 0.25353 6.516 2.42e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 21.18 on 107 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.6059, Adjusted R-squared: 0.5948
F-statistic: 54.83 on 3 and 107 DF, p-value: < 2.2e-16
(Intercept): 切片 今回は負の切片を持つ
Solar.R : オゾンに対し正の相関を持つ 今回は0.06なのだがこれは大きいのだろうか小さいのだろうか?
Wind: オゾンに対し負の相関を持つ
Temp: オゾンに対し正の相関を持つ
まとめると，日照強度(Solar.R)と気温(Temp)が高いほどオゾン濃度(Ozone)が増し，
風速(Wind)が大きいほどオゾン濃度が下がることが判明した．
ビールの売上と花火，気温，CMの関係を調べる
元データとggplotでの図示
> head(beer)
Revenue CM Temp Firework
1 47.14347 141 31 2
2 36.92363 144 23 1
3 38.92102 155 32 0
4 40.46434 130 28 0
5 51.60783 161 37 0
6 32.87875 154 27 0
> mbeer <- melt(beer)
> mbeer$x <- c(1:30)
> head(mbeer)
variable value x
1 Revenue 47.14347 1
2 Revenue 36.92363 2
3 Revenue 38.92102 3
4 Revenue 40.46434 4
5 Revenue 51.60783 5
6 Revenue 32.87875 6
> p <- ggplot(mbeer, aes(x=x, y=value, fill=variable, color=variable)) + geom_line()
> p
> p + facet_grid(variable ~ .)
> p + facet_grid(variable ~ ., scales = "free_y")
y軸をフリーにし，分けてみると見やすい．
花火は直接的に影響ありそうに見える．テレビCMも少しタイミングが遅れて関係しているように見える．
次に重回帰分析を行う．
> beer.lm <- lm(Revenue ~ ., beer)
> summary(beer.lm)
Call:
lm(formula = Revenue ~ ., data = beer)
Residuals:
Min 1Q Median 3Q Max
-6.028 -3.038 -0.009 2.097 8.141
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.23377 12.40527 1.389 0.17655
CM -0.04284 0.07768 -0.551 0.58602
Temp 0.98716 0.17945 5.501 9e-06 ***
Firework 3.18159 0.95993 3.314 0.00271 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.981 on 26 degrees of freedom
Multiple R-squared: 0.6264, Adjusted R-squared: 0.5833
F-statistic: 14.53 on 3 and 26 DF, p-value: 9.342e-06
CMは小さな負の値，つまり薄い負の相関になっているが，図を見る限りタイミングが遅れて効果があるように見える． 重回帰分析は，直線的な予測をするモデルなので，タイミングの遅れには対応できないのかもしれない． このような場合にはCMの効果があるまでの間を把握し，データをもっと大きく区切ることが有効だと思う． 例えばCMの効果が出でるまでに3日かかるのであれば，1日単位のデータを結合して，3日単位のデータにすればいいと思う． この辺りはRでやるのは大変なのでスクリプト言語でやったほうが簡単だろう．
相関係数と重回帰分析
相関係数を求めて，図示する
> beercor <- cor(beer)
> beercor
Revenue CM Temp Firework
Revenue 1.00000000 -0.07355843 0.67461456 0.43715932
CM -0.07355843 1.00000000 0.06244943 -0.12010466
Temp 0.67461456 0.06244943 1.00000000 0.04295698
Firework 0.43715932 -0.12010466 0.04295698 1.00000000
> mbeercor <- melt(beercor)
> mbeercor
X1 X2 value
1 Revenue Revenue 1.00000000
2 CM Revenue -0.07355843
3 Temp Revenue 0.67461456
4 Firework Revenue 0.43715932
5 Revenue CM -0.07355843
6 CM CM 1.00000000
7 Temp CM 0.06244943
8 Firework CM -0.12010466
9 Revenue Temp 0.67461456
10 CM Temp 0.06244943
11 Temp Temp 1.00000000
12 Firework Temp 0.04295698
13 Revenue Firework 0.43715932
14 CM Firework -0.12010466
15 Temp Firework 0.04295698
16 Firework Firework 1.00000000
> p <- ggplot(mbeercor, aes(X1, X2)) +
> geom_tile(aes(fill=value)) +
> geom_text(aes(fill=value, label = round(value,3)),colour="white")
> p
Revenueと他関数の関係は
Revenue CM Temp Firework
Revenue 1.00000000 -0.07355843 0.67461456 0.43715932
となる．
一方で正規化した重回帰分析の値は
> beer.lm$coefficients
(Intercept) CM Temp Firework
17.23377359 -0.04283877 0.98715569 3.18159007
となる．
あれ? 結構値が違う気がするのだけど…教科書と違う…
なんにせよ，重回帰分析のほうが最小二乗法を使っているため正確らしい．
