DS2405_1
Hayashi Atsuhiro
2024/10/31
0. 前回のアンケートから: 18名
- 欠損値の存在とその扱い。取捨選択するための工夫が要る。
- ちょうど良いスピード。何とか付いていけそう。
- 講義内容が難しくなった。時間がギリギリ。 ===> 復習の重要性
- 何の目的のために、どの関数を使うのかの理解。
- 課題の期限は、11月20日を考えている。石井先生の講義と被らないように。
1. 回帰分析
今回は、多変量解析の代表的な手法である回帰分析について解説する。
工学系や農学系の実験等を行う領域では頻繁に使用される手法であるが、 日常的な話題の中でも概念は広く利用されているので、取っ付き易い手法ではないだろうか。近年は機械学習等でも取り扱われていると思う。過去のデータからその構造を把握し、新規に測定されたデータに対する予測を 行ないたいと言うときに、回帰分析は有用である。 構造のシンプルな単回帰分析でこの手法の原理を理解し、 複数の説明変量を用いた重回帰分析に拡張する。 残差の取り方や、その二乗和を最少にするという考えは同じである。
1.1. アイディア
散布図にもっともらしい「直線」を当てはめたい。
- 何に使う?
- 何に使える?
- 散布図: 縦軸=体重、横軸=身長
皆さんから収集した体格データを用いて、散布図を描いてみよう。どういう直線が「もっともらしい」と考えるか?
-
[注意1]
散布図を丁寧に観るとデータの稠密な部分と疎な部分があることを理解せよ。
[注意2] ここに掲載している散布図は、視覚的に密度を理解し易くするために、第3回の第3節で紹介した、乱数を用いて描画点をブレさせて表示してある。
- どの様に直線の位置決めをすれば良い?
- 一意に決めるにはどの様なアイディアがある?
- そもそも、どういう考えから導き出された?
-
回帰分析とは?
- なぜ僕らは分布の中央付近を通過する直線が相応しいと感じるのか?
- それはどうやって決めたのか?
1.2. 前準備: 完全データの抽出
データに欠損値が含まれていると、後々の分析(残差プロット等)時に、 毎回欠損値を除去する手順を組み込む必要が生じて面倒なので、 分析対象の変量だけを先に切り出し、その上で、欠損値を除外(完全データ)しておく。
# データの切り出し
colnames(Student24)## [1] "No" "Sex" "Height" "Weight" "Chest"
## [6] "Residence" "Remittance" "Carrier" "Fee"StudTmp1<-Student24[,3:4] # 身長(3番目)と体重(4番目)だけを切り出す
StudTmp1[1:15,]## Height Weight
## 1 145.0 38
## 2 145.5 42
## 3 146.7 41
## 4 148.0 42
## 5 148.0 43
## 6 148.9 NA
## 7 149.0 45
## 8 150.0 43
## 9 150.0 46
## 10 150.0 47
## 11 151.0 42
## 12 151.0 45
## 13 151.0 46
## 14 151.0 50
## 15 151.0 NAdim(StudTmp1)## [1] 792 2元々の身長と体重のデータは 792例だったが、欠損値を含むデータが101例あり、 残ったのは691例である。
StudHW<-na.omit(StudTmp1) # 欠損値を含むデータを削除
dim(StudHW)## [1] 691 2colnames(StudHW)## [1] "Height" "Weight"StudHW[1:15,]## Height Weight
## 1 145.0 38.0
## 2 145.5 42.0
## 3 146.7 41.0
## 4 148.0 42.0
## 5 148.0 43.0
## 7 149.0 45.0
## 8 150.0 43.0
## 9 150.0 46.0
## 10 150.0 47.0
## 11 151.0 42.0
## 12 151.0 45.0
## 13 151.0 46.0
## 14 151.0 50.0
## 16 151.7 41.5
## 17 152.0 35.0以後は、Student24に代えて、欠損値を除去して完全データにしたStudHWを使って分析を行う。
【蛇足: パズルの要領で】上記では、欠損値を含むデータの例数を、101例=792-691として求めたが、 欠損値を含むデータの例数を直接的に数える式を作ってみると、例えば以下のようになる(他にも実現方法はあると思う)。 sum()関数は、本来は引数(カッコ内の式)の合計値を計算する関数だが、論理値 TRUE/FALSEの場合はTRUEを1、FALSEを0と捉えてその総数(=個数)を算出してくれる。 2通りの方法で実現してみたが、それぞれがどのように機能しているか理解できるであろうか? なお、コメントアウト(行頭に＃を付けてある)行は確認のために変数の内容を表示させる行である。表示させたければ、行頭の＃を削除すれば良い。
[実現方法1]
# 欠損値を含むデータの例数カウント(方法1)
StudTmp2<-!is.na(StudTmp1)
# StudTmp2[1:15,]
StudTmp3<-StudTmp2[,1] & StudTmp2[,2]
# StudTmp3[1:15]
# (StudTmp3==F)[1:15]
sum(StudTmp3==F)## [1] 101[実現方法2]
# 欠損値を含むデータの例数カウント(方法2)
StudTmp2<-is.na(StudTmp1)
# StudTmp2[1:15,]
StudTmp3<-StudTmp2[,1] | StudTmp2[,2]
# StudTmp3[1:15]
sum(StudTmp3) ## これと同じ意味 ===> sum(StudTmp3==T)## [1] 101
2. 単回帰分析 : 予測等に使う、連続変量の関係
- 体重を身長で説明(回帰)したい : [体重]=a+b[身長] : 回帰係数
- 関係性は直線＝線形を仮定する
- 説明される変量 : 目的変数、従属変数、dependent variable
-
説明する変量 : 説明変数、独立変数、independent variable
- 誤差の取り方 : 指定された独立変数における測定誤差。だから誤差は垂直方向に分布する。
- 回帰直線からのズレ(誤差)のことを回帰分析では「残差」と呼ぶ。
- 残差の分布は正規分布を仮定する
- 測定の場所によらず、残差は同一の分布であることを仮定する
- どうやって直線を決める? : 予測誤差の2乗和を最小にする
2.1. 定式化
-
[余談] 回帰分析では「2乗和」を最小にすることを考えるが、「絶対値和」とか「符号付き和(1乗和)」を最小にする方法もアイディアとしてはあり得る。 2乗和だと式の展開が楽になる。
-
説明する変数の個数
- 一つ: 単回帰分析: 身長で説明したい : [体重]=a+b[身長]
- 複数: 重回帰分析: 身長と胸囲で説明したい : [体重]=a+b[身長]+c[胸囲]
回帰平面の方程式: \[ \textbf{y} = a + b \textbf{x} \]
測定値と予測値のズレ: \[ \varepsilon_{i} = y_{i} - (a+ b x_{i}) \;\;\;\;\;\; (i=1,2,...,n) \]
ズレの2乗の和を最小に: \[ S = \sum_{i=1}^{n}{\varepsilon_{i}^2} = \sum_{i=1}^{n} \{ {y_{i} - ( a + b x_{i}})\}^2 \;\;\; { \to min} \]
# 回帰分析(指定方法1)
RresultHW <- lm(Weight ~ Height,
data=StudHW) # 線形モデル, Linear Model
RresultHW # 簡略化した計算結果表示##
## Call:
## lm(formula = Weight ~ Height, data = StudHW)
##
## Coefficients:
## (Intercept) Height
## -79.5566 0.8228なお、以下のように記述しても同じ結果が得られる。 こちらの方が表記としては理解し易いかもしれないが、長くなるのが難点か。
# 回帰分析(指定方法2)
RresultTmp1 <- lm(StudHW$Weight ~ StudHW$Height,
data=StudHW) # 線形モデル, Linear Model
RresultTmp1 # 簡略化した計算結果表示- 回帰直線の係数はそれぞれ -79.6と0.823と求められた。
- これを使えば、教室の入り口にメジャーを貼り付けておけば体重を言い当てられるはず。
[演習2.1-1] 試しに自分の体重がうまく予測できているか計算してみよ。
- -79.6+0.823*178=66.9Kg。
- 残差は 86.5-66.9=19.6Kg。
- うーん、メタボ健診に引っかかるわけだ。
- 部分的に係数や残差を取り出して表示することもできる。
coef <-RresultHW$coefficient # 計算結果の取り出し(回帰係数)
coef # 回帰係数## (Intercept) Height
## -79.5565674 0.8228171resid <-RresultHW$residuals # 計算結果の取り出し(残差)
resid[1:20] # 残差(先頭の20例)## 1 2 3 4 5 7
## -1.7519139 1.8366775 -0.1507030 -0.2203652 0.7796348 1.9568176
## 8 9 10 11 12 13
## -0.8659995 2.1340005 3.1340005 -2.6888166 0.3111834 1.3111834
## 14 16 17 18 19 20
## 5.3111834 -3.7647886 -10.5116337 -2.5116337 -1.5116337 -0.5116337
## 21 22
## -5.3344508 -4.3344508- 散布図に回帰直線を付加する方法。
plot(Weight ~ Height, data=StudHW,
xlab = "Height", ylab = "Weight",
main="Scatter Plot of Shintyou and Taijyuu") # 散布図の描画(これまでと別の指定方法)
abline(RresultHW) # 回帰直線を付加
abline(h=seq(40,100,10), lty=3)
abline(v=seq(140,190,5), lty=3)
3. 回帰係数以外の重要な情報
回帰分析について、回帰係数以外にも知っておくべき情報がある。 それらは「summary(RresultHW)」の出力中に表示されている。
summary(RresultHW) # より詳細な計算結果表示##
## Call:
## lm(formula = Weight ~ Height, data = StudHW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.665 -4.671 -1.019 3.437 34.386
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -79.5566 6.1948 -12.84 <2e-16 ***
## Height 0.8228 0.0365 22.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.188 on 689 degrees of freedom
## Multiple R-squared: 0.4244, Adjusted R-squared: 0.4236
## F-statistic: 508.1 on 1 and 689 DF, p-value: < 2.2e-16-
決定係数 : Adjusted R-squared
- 目的変量が説明変量でどの程度説明しているかの割合。
-
1 に近いほど当てはまりが良いと言える。
Adjusted R-squared: 0.4236
-
説明変量が予測に役立っているかの検定
-
回帰に役立っているか : p-value : 小さいと有意(役立っている)
F-statistic: 508.1 on 1 and 689 DF, p-value: < 2.2e-16
-
回帰に役立っているか : p-value : 小さいと有意(役立っている)
-
回帰係数がゼロではないかの検定
-
回帰係数の検定(係数=0 か?) : Pr(>|t|) :
小さいと有意(係数=0ではないと言える)
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -79.5566 6.1948 -12.84 <2e-16 *** Height 0.8228 0.0365 22.54 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-
回帰係数の検定(係数=0 か?) : Pr(>|t|) :
小さいと有意(係数=0ではないと言える)
4. 残差分析もお忘れなく
-
回帰直線の係数が求まったらそれでおしまい?
-
仮定が3つあったはず
- 関係性は直線＝線形を仮定する
- 誤差の分布は正規分布を仮定する
- 測定の場所によらず、残差は同一の分布(正規分布)であることを仮定する
-
仮定が成り立っているかは確認せねば。どうやって?
- もし、残差が何らかの傾向を示しているなら、直線＝線形関係ではないのかもしれないと疑う必要がある。 ==> バナナカーブ
- 測定する場所で残差の傾向が異なるなら、残差の分布が同一とは言えないかもしれない。 ==> 順に広がっているラッパカーブ
- 残差のヒストグラム、箱ひげ図で偏りが判定できる。
-
Q-Q
プロットの斜め直線に乗っていれば正規分布と言える。ズレが大きいようだと検討の余地がある。
- 程度問題ではあるが、少しぐらいの逸脱は許容範囲と言える。
- 残差が正規分布をするということを意識してもらえば良い。
# 残差プロット: 対身長に対して
plot(StudHW$Height, RresultHW$residuals) # 残差プロット(対身長)
abline(v=seq(140,190,5),lty=3) # 点線を追記
abline(h=seq(-10,30,10),lty=3) # 点線を追記
abline(h=0,lty=1) # 原点を追記
# 残差プロット: 対体重に対して
plot(StudHW$Weight, RresultHW$residuals) # 残差プロット(対体重)
abline(v=seq(30,100,5),lty=3) # 点線を追記
abline(h=seq(-10,30,10),lty=3) # 点線を追記
abline(h=0,lty=1) # 原点を追記
# 残差のヒストグラム
hist(RresultHW$residuals, right=F)
abline(h=seq(0,200,25),lty=3) # 点線を追記
# 残差の箱ひげ図
boxplot(RresultHW$residuals, horizontal=T)
abline(v=seq(-10,30,10),lty=3) # 点線を追記
- 残差のヒストグラム、箱ひげ図で偏りが判定できる。
# 残差のQQプロット: 正規確率プロット
qqnorm(RresultHW$residuals) # QQプロット(正規確率プロット)
qqline(RresultHW$residuals, lty=2) # 基準線の描画
abline(h=0, lty=1) # 実線を追記
abline(v=0, lty=1) # 実線を追記
abline(h=seq(-10,30,10),lty=3) # 点線を追記
-
Q-Q
プロットの斜め直線に乗っていれば正規分布と言える。ズレが大きいようだと検討の余地がある。
- [検討事項] 残差が20を超える辺りの(少なくとも)9例ほどは分布を乱しているように見える。===> 外れ値か? 吟味が必要。この後で解説。
5. この分析のまとめ(結果の見方)
-
対象になったのは 691名。
-
説明変量が予測に役立っているか?
-
回帰に役立っているか : p-value : 小さいと有意(役立っている)
-
[この例] 1% 未満(0.01%) なので役に立っていると言える。
-
回帰に役立っているか : p-value : 小さいと有意(役立っている)
-
決定係数 : Adjusted R-squared: ( 相関係数 : R )
- 目的変量が説明変量でどの程度説明しているかの割合。
-
1 に近いほど当てはまりが良いと言える。
-
[この例] 43% 程(約半分)を説明できている。
-
回帰係数 : Estimate
-
[この例] a=-79.6, b=0.823
-
説明変数が予測に役立っているか?
回帰係数の検定(係数=0 か?) : Pr(>|t|) : 小さいと有意(係数=0ではないと言える)-
[この例] 両者とも 1% 未満(0.01%)
なので回帰係数はゼロではない(何らかの意味がある数字と言える)。
-
残差の性質 ===> 正規性 : 残差プロット、残差解析
- 残差(予測誤差)は正規分布をしていると仮定してモデルが構築されている。
- この仮定が覆ると、回帰分析として成立していないことになる。
- 残差が正規分布をしているか確認する必要がある。
- 均等に散らばっているか?
- 傾向はないか?
-
…
[この例] 残差には概ね傾向は見られない。 ただし体重の大きい 9例程度は要確認。場合によっては外れ値として除外も。
[演習5-1] 上述では[体重]を[身長]で説明する回帰式の例で説明した。では、[体重]を[胸囲]で説明するとどうなるか? 各自でやってみよ。
[演習5-2] スマホ月額料金を仕送り額で説明する回帰分析等、興味のある組み合わせで実行してみよ。
6. 有効桁数に注意せよ : どこまでが「意味ある桁」か?
測定精度上回る計算結果は表示させることはできても、意味はない。
プログラムの出力をそのまま書き写さないように。
[重要な注意]
統計ソフトは単なる道具。使いこなすのは各自。
-
12.3 <=== 12.25～12.34
12 <=== 11.5 ～12.4
67.8 <=== 67.75～67.84
68 <=== 67.5 ～68.4
[例2] 日本の観測史上の最高気温は、
2018(平成30)年7月23日に熊谷市で観測された41.1度であり、
最低気温は、1902(明治35)年1月25日に北海道旭川市の-41度であった。===>
-41.0度
[例3]
2001年のイチロー選手の打率は3割5分であった。
2006年は3割3分1厘であった。===> 3割5分0厘
【今週の思い】回帰分析のアイディアを理解し、分析手順を習得せよ。
81. 参考
このページで取り扱ったプログラムだけを抜き出して以下に列挙しておく。
## 1.1. アイディア
tmpRand<-rnorm(2000)
tmpRand[1:40]
StudTmp<-Student24
StudTmp$Height <-Student24$Height+tmpRand[1:772]*0.2
StudTmp$Weight <-Student24$Weight+tmpRand[1001:1772]*0.2
# 乱数で位置をずらした散布図
plot(StudTmp$Height, StudTmp$Weight,
xlab = "Height", ylab = "Weight",
main="Scatter Plot of Shintyou and Taijyuu")
abline(h=seq(40,100,10), lty=3)
abline(v=seq(140,190,5), lty=3)
## # 元データそのものの散布図(点が重なって判読できない)
## plot(Student24$Height, Student24$Weight)
## abline(h=seq(40,100,10), lty=3)
## abline(v=seq(140,190,10), lty=3)
RresultTmp <- lm(Weight ~ Height, data=StudTmp)
## RresultTmp
## summary(RresultTmp)
plot(Weight ~ Height, data=StudTmp,
xlab = "Height", ylab = "Weight",
main="Scatter Plot of Shintyou and Taijyuu")
abline(RresultTmp)
abline(h=seq(40,100,10), lty=3)
abline(v=seq(140,190,5), lty=3)
## 1.2. 前準備(Rでの利用に際して)
# データの切り出し
colnames(Student24)
StudTmp1<-Student24[,3:4] # 身長(3番目)と体重(4番目)だけを切り出す
StudTmp1[1:15,]
dim(StudTmp1)
StudHW<-na.omit(StudTmp1) # 欠損値を含むデータを削除
dim(StudHW)
colnames(StudHW)
StudHW[1:15,]
# 欠損値を含むデータの例数カウント(方法1)
StudTmp2<-!is.na(StudTmp1)
# StudTmp2[1:15,]
StudTmp3<-StudTmp2[,1] & StudTmp2[,2]
# StudTmp3[1:15]
# (StudTmp3==F)[1:15]
sum(StudTmp3==F)
# 欠損値を含むデータの例数カウント(方法2)
StudTmp2<-is.na(StudTmp1)
# StudTmp2[1:15,]
StudTmp3<-StudTmp2[,1] | StudTmp2[,2]
# StudTmp3[1:15]
sum(StudTmp3) ## これと同じ意味 ===> sum(StudTmp3==T)
## 2.1. 定式化
# 回帰分析の表記方法(その1)
RresultHW <- lm(Weight ~ Height,
data=StudHW) # 線形モデル, Linear Model
RresultHW # 簡略化した計算結果表示
# 回帰分析の表記方法(その2)
RresultTmp1 <- lm(StudHW$Weight ~ StudHW$Height,
data=StudHW) # 線形モデル, Linear Model
RresultTmp1 # 簡略化した計算結果表示
coef <-RresultHW$coefficient # 計算結果の取り出し(回帰係数)
coef # 回帰係数
resid <-RresultHW$residuals # 計算結果の取り出し(残差)
resid[1:20] # 残差(先頭の20例)
# 散布図に回帰直線を付加する方法
plot(Weight ~ Height, data=StudHW,
xlab = "Height", ylab = "Weight",
main="Scatter Plot of Shintyou and Taijyuu") # 散布図の描画(これまでと別の指定方法)
abline(RresultHW) # 回帰直線を付加
abline(h=seq(40,100,10), lty=3)
abline(v=seq(140,190,5), lty=3)
## 3. その他の重要な情報
summary(RresultHW) # より詳細な計算結果表示
## 4. 残差分析もお忘れなく
# 残差プロット: 対身長に対して
plot(StudHW$Height, RresultHW$residuals) # 残差プロット(対身長)
abline(v=seq(140,190,5),lty=3) # 点線を追記
abline(h=seq(-10,30,10),lty=3) # 点線を追記
abline(h=0,lty=1) # 原点を追記
# 残差プロット: 対体重に対して
plot(StudHW$Weight, RresultHW$residuals) # 残差プロット(対体重)
abline(v=seq(30,100,5),lty=3) # 点線を追記
abline(h=seq(-10,30,10),lty=3) # 点線を追記
abline(h=0,lty=1) # 原点を追記
# 残差のヒストグラム
hist(RresultHW$residuals, right=F)
abline(h=seq(0,200,25),lty=3) # 点線を追記
# 残差の箱ひげ図
boxplot(RresultHW$residuals, horizontal=T)
abline(v=seq(-10,30,10),lty=3) # 点線を追記
# 残差のQQプロット: 正規確率プロット
qqnorm(RresultHW$residuals) # QQプロット(正規確率プロット)
qqline(RresultHW$residuals, lty=2) # 基準線の描画
abline(h=0, lty=1) # 実線を追記
abline(v=0, lty=1) # 実線を追記
abline(h=seq(-10,30,10),lty=3) # 点線を追記