スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

重回帰分析

前回は単回帰分析の内容とRでの実装方法を勉強したので、今日は重回帰分析を勉強したいと思います。

■ビジネスにおける適応例
重回帰分析では、複数の要因が結果に対して効いているのか、いないのかも判断できるので、モデルを作るだけではなく必要な要因は何なのか、それが結果に対してどのように関係しているのかを明らかにできる。そのため、ビジネスでも複数の要因を説明変数として売り上げに寄与している要因を調べるなどの使い方がされている。

■重回帰分析のおおまかな処理の流れ
①モデルの当てはまりの良さの確認
→自由度調整済み決定係数を確認(クロスセクションデータ:0.5以上、タイムシリーズデータ:0.7以上であれば当てはまりがよい)
※重回帰分析の場合、説明変数を増やせば、決定係数が高くなるという性質があります。そこで、重回帰分析のモデルの当てはまりの良さを図る指標としては自由度調整済み決定係数を使用します。
②各変数の統計的優位性の確認
→各変数のP値を確認(ここで算出されるP値は目的変数と各説明変数でt検定を実施し、「説明変数の影響力は0である」という帰無仮説が棄却される際の危険率)
③各変数の偏回帰係数を解釈
→符号条件が理論的におかしくないかを確認

このながれにそってRのairqualityをサンプルデータとして重回帰分析を実施。
オゾン濃度に対してどの変数が影響しているかを分析してみることにします。

■Rで重回帰
airqualityのデータの中身を確認してみる。
#データの読み込み
data(airquality)
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

早速、オゾンの量を目的変数、残りの太陽光、風、気温を説明変数にして重回帰のモデルを作成。どの説明変数の組み合わせを使用した際に一番制度の高いモデルが作れるのかわからないので、step関数を使用してステップワイズ法でAICが一番小さな値になるモデルを採用することにします。
#オゾン濃度を説明するモデルを作成する
#年月のデータを省いて分析
airq <- airquality[,1:4]
#重回帰分析のモデルを作成
airq.lm <- lm(Ozone~.,data=airq)
#ステップワイズ法でモデルを評価
step(airq.lm)

Start:  AIC=681.7
Ozone ~ Solar.R + Wind + Temp

Df Sum of Sq RSS AIC
<none> 48003 682
- Solar.R 1 2986 50989 686
- Wind 1 11642 59644 704
- Temp 1 19050 67053 717


Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp, data = airq)

Coefficients:
(Intercept) Solar.R Wind Temp
-64.3421 0.0598 -3.3336 1.6521

今回は、すべての説明変数を使用した場合が一番AICが小さいモデルになるみたいですね。
では、ここから先ほどの「重回帰 
分析のおおまかな処理の流れ」に従ってモデルの評価、結果の確認をしたいと思います。
#基本等計量を確認
summary(airq.lm)


Call:
lm(formula = Ozone ~ ., data = airq)

Residuals:
Min 1Q Median 3Q Max
-40.48 -14.22 -3.55 10.10 95.62

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -64.3421 23.0547 -2.79 0.0062 **
Solar.R 0.0598 0.0232 2.58 0.0112 *
Wind -3.3336 0.6544 -5.09 1.5e-06 ***
Temp 1.6521 0.2535 6.52 2.4e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.2 on 107 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.606, Adjusted R-squared: 0.595
F-statistic: 54.8 on 3 and 107 DF, p-value: <2e-16

①モデルの当てはまりの良さの確認
自由度調整済み決定係数の値が0.595>0.5なので、当てはまりの良さ的には問題なさそうです。

②各変数の統計的優位性の確認
Wind、Tempについては、p値が非常に小さな値になっているので有効な変数みたいです。ただ、Solar.Rのp値が0.0112>0.005なので、この変数については、信頼度は薄そうです。とはいえ、ステップワイズ法でAICの確認もしているのでモデルとしてはSolar.Rを説明変数に入れておいた方が当てはまりのよいモデルになるようです。

③各変数の偏回帰係数を解釈
オゾン濃度には気温がプラスの影響を、風速がマイナスの影響を与えていることが分かりますね。


■参考
ビジネスデータ分析研究室
重回帰分析
Rと重回帰分析
http://www.aoni.waseda.jp/abek/document/regression-2.html
スポンサーサイト

統計学の復習(擬似相関と層別解析)

■擬似相関と偏相関係数
擬似相関とは、変数Aと変数B、変数Aと変数Cの間に相関関係があるとき、実際には相関関係がないにも関わらず変数Bと変数Cの間にも相関関係が生じること。擬似相関が生じていると判断を誤ってしまうことがあるので注意が必要です。
擬似相関については、このサイトにいい例とRのコードが記載されてました。擬似相関、統計的消去について非常にわかりやすくかかれているので参考にしてください。
上記サイトの例では、変数Aと変数Bの残差、変数Aと変数Cの残差をそれぞれ算出して残差同士の相関を確認していますが、偏相関係数を使用すると統計的消去を行い、直接変数間の相関関係を算出することが出来ます。
上記サイトで使用しているデータをお借りして偏相関係数を算出。

#上記サイトで使用しているデータをお借りします。
data <- read.table("http://dl.dropboxusercontent.com/u/432512/20130708/SchoolChildrensMaths.txt", header = TRUE)
#偏相関係数を使用するため、ソースの読み込みを行う
source("http://aoki2.si.gunma-u.ac.jp/R/src/partial_cor.R", encoding="euc-jp")
#相関係数を算出
cor(data)
#偏相関係数を算出
partial.cor(data)

▼結果
> cor(data)

AMA YEARS HGT
AMA 1.0000000 0.9915810 0.9799811
YEARS 0.9915810 1.0000000 0.9889419
HGT 0.9799811 0.9889419 1.0000000

> partial.cor(data)

Var 1 Var 2 Var 3
Var 1 NA 0.7598976 -0.03306099
Var 2 0.75989764 NA 0.66762606
Var 3 -0.03306099 0.6676261 NA

普通に相関係数を算出した場合だとすべての変数間で非常に高い相関値が算出されているのに対して、偏相関係数を算出した場合は擬似相関を除いた相関値が出力されていることがわかります。

■層別解析
お次は層別解析。層別解析とは、データをグループ別に分けて解析すること。解析対象のデータに正確が異なるいくつかの部分集団を含んでいる場合全体では相関関係は現れないが、グループ別に分析すると相関関係が洗われることがあるため、データの相関関係をみる際にはあらかじめ注意して行う必要がある。
Rではby()関数を使用することで、層別解析を行うことが出来る。

#データの読み込み
data<-read.csv("http://fileman.rakurakuhp.net/UserFiles/40164/File/1199162168.csv")
#変数のアタッチ
attach(data)
#全データの要約統計量の算出
summary(data)
#性別ごとの要約統計量の算出
by(data, SEX, summary)

▼結果
> summary(data)
     SUBJID    SEX         AGE       ARMCD     WEIGHT         HEIGHT   
Min. :101 F:15 Min. :20.0 A:8 Min. : 8.0 Min. :146
1st Qu.:108 M:16 1st Qu.:27.5 B:8 1st Qu.:65.0 1st Qu.:158
Median :116 Median :35.0 C:8 Median :75.0 Median :167
Mean :116 Mean :35.0 D:7 Mean :69.9 Mean :169
3rd Qu.:124 3rd Qu.:42.5 3rd Qu.:78.0 3rd Qu.:178
Max. :131 Max. :50.0 Max. :98.0 Max. :198

> by(data, SEX, summary)
SEX: F
SUBJID SEX AGE ARMCD WEIGHT HEIGHT
Min. :102 F:15 Min. :21 A:0 Min. :58.0 Min. :146
1st Qu.:109 M: 0 1st Qu.:28 B:8 1st Qu.:70.5 1st Qu.:162
Median :116 Median :35 C:0 Median :76.0 Median :167
Mean :116 Mean :35 D:7 Mean :76.6 Mean :168
3rd Qu.:123 3rd Qu.:42 3rd Qu.:82.5 3rd Qu.:176
Max. :130 Max. :49 Max. :98.0 Max. :186
--------------------------------------------------------
SEX: M
SUBJID SEX AGE ARMCD WEIGHT HEIGHT
Min. :101 F: 0 Min. :20.0 A:8 Min. : 8.0 Min. :148
1st Qu.:108 M:16 1st Qu.:27.5 B:0 1st Qu.:55.5 1st Qu.:158
Median :116 Median :35.0 C:8 Median :67.0 Median :167
Mean :116 Mean :35.0 D:0 Mean :63.7 Mean :169
3rd Qu.:124 3rd Qu.:42.5 3rd Qu.:76.0 3rd Qu.:180
Max. :131 Max. :50.0 Max. :87.0 Max. :198


by()関数一発で、性別ごとの要約統計量を得ることが出来ました。これは便利!

■参考
偏相関係数
ほくそ笑む
学びing

統計学の復習(度数分布とヒストグラム)

久しぶりに統計学の教科書を開いたらいろいろと忘れていたので、一から復習をすることにします。
統計学の本をもとにRを使って実践しつつ基本をおさらいしていきます。

教科書はこの本。


■基本用語
ヒストグラム:度数分布表をもとに作成した柱状グラフ。調査や実験に使用するデータセットが手に入ったらまずは基礎統計量ヒストグラムを作成して全体のデータの分布状況を確認する。
階級値:階級を代表する値のこと。階級の上限値と下限値の中間値を階級値とするのが一般的
相対度数:全体の大きさを1とした時の各階級に属する観測値の個数の全体中での割合。正規かされた値なので、データの大きさが異なる複数のデータの分布の比較を行う時に有効。

■スタージェスの公式
度数分布表やヒストグラムを作成する時に注意するべき点は階級数階級幅の問題。特に階級数は多すぎても少なすぎても得られるヒストグラムから読み取れることを変わってきてしまうので、身長に階級数を決定する必要がある。階級数をどうやって決めるかはルール化されていないが、スタージェスの公式というものを使えばデータセットから最適な階級数を算出してくれる。

Rの標準のデータセット「airmiles」を使って最適な階級数を確認。
#データセットの設定
data(airmiles)
airmiles
#スタージェスの公式による階級数の確認
nclass.Sturges(airmiles)

▼結果
[1] 6

上記の結果から「airmiles」のデータセットだと、6刻みで階級数を決定するのがいいらしい。


■Rでのヒストグラムの作成方法
Rではヒストグラムの横幅はデフォルトでは スタージェスアルゴリズムによって決められるらしい。横幅を変更する場合は、breaks オプションに分割数、横幅のサイズまたはアルゴリズムを指定する。アルゴリズムを指定するとき、Sturges、Scott、FD、Freedman-Diacoins のを設定する。


#ヒストグラムの作成
hist(airmiles)

▼結果
histgram of airmiles

なるほど、ちゃんとスタージェスの公式で算出した通り6分割でヒストグラムが作成されているのがわかります。

■参考
バイオスタティスティクス
R tips
ほくそ笑む

線形単回帰分析

今日は線形単回帰分析。

線形回帰分析:説明変数(x)を用いて目的変数(y)を説明する統計モデル(y=f(x))をデータから求める分析手法。
例えば、「車のスピード(x)」と「ブレーキをかけた後に止まるまでの距離(y)」の関係を数式(y=f(x))で表すことが出来る。この関係を数式で表すことにより、時速50kmで車が走ってきたときに約何mくらいで止まるのかという予測が出来るようになる。

■処理ステップ


線形回帰分析は大まかに下記の5ステップで実施。
1.説明変数と目的変数の相関関係を確認
2.回帰式の算出
3.回帰式の精度を確認
4.回帰係数の検定
5.信頼区間と予測区間の算出

で、線形回帰分析の概念自体は特に難しい話ではないので、実際にRに標準で用意されているデータセットcarsを使って線形回帰分析に挑戦。

1.説明変数と目的変数の相関関係を確認


無相関のデータに対して線形回帰をやっても意味がないので、「車のスピード(x)」と「ブレーキをかけた後に止まるまでの距離(y)」相関係数を算出。
cor(cars$speed, cars$dist)


▼cor(cars$speed, cars$dist)の結果
0.8068949

0.8なので、「車のスピード(x)」と「ブレーキをかけた後に止まるまでの距離(y)」の間には強い正の相関があることを確認。

2.回帰式の算出


cars.lm <- lm(dist~speed, data=cars)
summary(cars.lm)

lm()は線形回帰分析の関数で「dist~speed」は目的変数distを説明変数speedで説明することを表しています。「data = cars」は分析対象のデータセットはcarsであることを表しています。

▼summary(cars.lm)の結果
Call:
lm(formula = dist ~ speed, data = cars)

Residuals:
Min 1Q Median 3Q Max
-29.069 -9.525 -2.272 9.215 43.201

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791 6.7584 -2.601 0.0123 *
speed 3.9324 0.4155 9.464 1.49e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438
F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12

「Residuals」は残差の最小値、第1四分位数、中央値、第3四分位数、最大値がそれぞれ表示されています。ちなみに、残差は回帰直線とデータとの誤差のこと。
t値とp値は係数が回帰モデルにやくだっているかどうかに関する統計量で、p値が大きいほどその係数が役にたっていないことを表している。いろいろと出てきているが一番知りたいのは作成した予測モデルがどれだけフィットしているか。で、それを表すのが決定係数、調整済み決定係数。これが1に近づくほどフィットしているとこを表します。
この回帰直線はy=3.9324x-17.5791のモデルで表現できるようです。

値だけでみててもよくわからんので、とりあえず散布図を作成。

plot(cars) # 散布図の作成
abline(cars.lm, lwd=2) # 回帰直線の作成

▼散布図と回帰直線
Rplot01.png
それらしい回帰直線が引かれている。

3.回帰式の精度を確認


残差を視覚的に分析するために、回帰診断図を表示する。
par(mfrow=c(2,2)) # 2x2のマトリックスで回帰診断図を表示
plot(cars.lm) # 回帰診断図を描画

▼回帰診断図
回帰診断図

よくわからない図が出現。。。図の見方を調べてみる。

■残差とフィット値のプロット(左上)
横軸が予測値、縦軸が残差をあらわしており、残差の全体像を外観するときに使用。

■残差の正規Q-Qプロット(右上)
データの正規性を確認するための図。Q-Qプロットはデータが正規分布に従うと点が直線上に並ぶ。回帰分析では、残差が標準正規分布に従うことを仮定しているため、Q-Qプロットを使用することで

■残差の平方根プロット(左下)
残差の変動状況を考察するための図。

■残差と影響力プロット(右下)
モデルの当てはまりへの影響力を測るための図。横軸は梃子値で、縦軸は標準化した残差。点線でクックの距離0.5を示している。クックの距離が0.5を超えると影響力あり、1を超えると特異に大きい。


■参考
分かって使う統計学 -相関と回帰分析-


多次元尺度構成法(MDS)

今日は多次元尺度構成法(Multi-Dimensional Scaling)のお勉強。

MDSは多次元データを2次元もしくは3次元上に可視化するための手法。


【処理ステップ】
●距離を求める。
●座標値を求める
●2~3次元上で個体を配置する(散布図を作成する)。
●信頼性について考察する
プロフィール

HitTips

Author:HitTips
FC2ブログへようこそ!

最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QR
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。