DS2406_7

0. 前回のアンケートから: 16名

• 回帰分析について深い理解。「目から鱗」。
• 身近なデータを用いた実演。
• いよいよついていけるか怪しくなってきた。 ===> 復習の重要性

• 求めたa, bが不適切と判断した場合の対処方法 ===> 今週のテーマ

7. 重回帰分析 : 予測等に使う、連続変量の関係

　ここまでは、構造がシンプルな単回帰分析を例に、この分析方法の考え方や残差の取り方、および、なぜそのように取る必要があるか等を説明した。この節では、より実用場面で登場する重回帰分析まで拡張し、分析の進め方や利用時の注意点を含めて紹介する。いろいろな場面での応用が可能な手法なので、大いに利用していただきたい。

• 2変量以上の説明する変量(説明変量)で 1変量(目的変量)を説明
• 説明変量が複数になる : 単 ===> 重
• 体重を身長と胸囲で説明(回帰)したい : [体重]=a+b[身長]+c[胸囲] : 回帰係数
  • 説明される変量 : 目的変数、従属変数、dependent variable
  • 説明する変量 : 説明変数、独立変数、independent variable

• アイディアは単回帰分析の時と全く同じ。
  • 両者の関係性は直線＝線形を仮定する
  • 回帰直線からのズレ(誤差)のことを回帰分析では「残差」と呼ぶ。
  • 説明される変量(目的変量)と平行に残差を取る。なぜなら、指定された独立変数における測定誤差だから。
  • 残差の分布は正規分布を仮定する
  • 測定の場所によらず、残差は同一の分布であることを仮定する
  • 残差の二乗和を最小にする(最小二乗法)

7.1. 定式化: 式の展開、解法

回帰平面の方程式: \[ \textbf{y} = a_{0} + \sum_{j=1}^{p} a_{j}\textbf{x}^t_{j} \]

測定値と予測値のズレ(残差): \[ \varepsilon_{i} = y_{i} - (a_{0}+\sum_{j=1}^{p}a_{j}x_{ij}) \;\;\;\;\;\; (i=1,2,...,n) \]

残差の2乗和を最小に: \[ S = \sum_{i=1}^{n}{\varepsilon_{i}^2} = \sum_{i=1}^{n} \{ {y_{i} - ( a_{0}+\sum_{j=1}^{p}a_{j}x_{ij}})\}^2 \;\;\; { \to min} \]

行列表記: \[ {\textbf{Y}}={\textbf{X} \textbf{A}} \]

ここで、 \[ \textbf{Y}= \begin{pmatrix} y_{1} \\ y_{2} \\ y_{3} \\ \vdots \\ y_{i} \\ \vdots \\ y_{n} \end{pmatrix}, \;\;\; %%% {\textbf{X}}= \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1j} & \cdots & x_{1p}\\ 1 & x_{21} & x_{23} & \cdots & x_{2j} & \cdots & x_{2p} \\ 1 & x_{31} & x_{32} & \cdots & x_{3j} & \cdots & x_{3p} \\ \vdots & \vdots & \vdots & \ddots & \vdots & & \vdots \\ 1 & x_{i1} & x_{i2} & \cdots & x_{ij} & \cdots & x_{ip} \\ \vdots & \vdots & \vdots & & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nj} & \cdots & x_{np} \\ \end{pmatrix}, \;\;\; %%% \textbf{A}= \begin{pmatrix} a_{0} \\ a_{1} \\ a_{2} \\ \vdots \\ a_{j} \\ \vdots \\ a_{p} \end{pmatrix}. \]

\(\textbf{A}\) を求めるには、左側から\(\textbf{X}^t\)をかけて、

\[ \textbf{X}^t {\textbf{Y}}={\textbf{X}^t {\textbf{X} \textbf{A}}} \]

その後、今度は左側から\(({{\textbf{X}^t {\textbf{X}}}})^{-1}\)をかけて、

\[ \textbf{A} = ({{\textbf{X}^t {\textbf{X}}}})^{-1} \textbf{X}^t {\textbf{Y}} \]

よって、\(\textbf{A}\) を求めるには、\(({{\textbf{X}^t {\textbf{X}}}})\)が逆行列を持つことが条件となる。

7.2. 前準備: 完全データの抽出

　単回帰分析の時も説明したが、 データに欠損値が含まれていると、後々の分析(残差プロット等)時に、 毎回欠損値を除去する手順を組み込む必要が生じて面倒なので、 分析対象の変量だけをまず切り出し、その上で、欠損値を除外(完全データ)しておく。

# データ切り出し
dim(Student24)

## [1] 792   9

StudHWC<-na.omit(Student24[,3:5]) # 欠損値を含むデータを削除
dim(StudHWC)

## [1] 242   3

colnames(StudHWC)

## [1] "Height" "Weight" "Chest"

StudHWC[1:10,]

##    Height Weight Chest
## 2   145.5   42.0    76
## 3   146.7   41.0    85
## 5   148.0   43.0    80
## 8   150.0   43.0    82
## 9   150.0   46.0    86
## 16  151.7   41.5    80
## 17  152.0   35.0    77
## 23  153.0   46.5    87
## 25  153.0   55.0    78
## 26  153.0   57.0    38

　元々の収集データは772サンプルであったが、 身長と体重と胸囲だけのデータを抽出し、 この中に欠損値を含むサンプルを除外すると、残ったのは239サンプルであった。 よって、欠損値を含んでいたのは533サンプルあったことになる。 変量が多くなると欠損値を含むサンプルが多くなる傾向にある。 以後は、Student24の代わりにStudHWCを使って分析を行う。

7.3. 散布図行列

　既に第4回で紹介したので再掲となるが、 3変量以上の場合に変量間の関係を概観するのに便利な「散布図行列」を再度描画しておく。 3変量の中から2変量ごとを組み合わせて描画した散布図を、 行列形式で配置したものである。 複数変量の関係を視覚的に把握するのに優れた手法であり、 変量間でどの程度関係性がありそうかを視覚的に把握することができて便利である。 例えば、データ分析の初期段階でまずは全体像を掴もうとするような際には有効に機能する。

# 散布図行列
plot(StudHWC)

7.4. プログラム ＆ 計算結果

dim(StudHWC)  # サンプルサイズを改めて確認

## [1] 242   3

# 重回帰分析
RresultHWC<-lm(Weight~Height+Chest, 
               data=StudHWC) # 右辺に説明する変量をプラスマークで追加していく

# 計算結果要約
summary(RresultHWC)

## 
## Call:
## lm(formula = Weight ~ Height + Chest, data = StudHWC)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.019  -4.341  -1.411   3.280  31.561 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -99.70771    9.92295 -10.048  < 2e-16 ***
## Height        0.80174    0.06103  13.136  < 2e-16 ***
## Chest         0.29231    0.04930   5.929 1.06e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.038 on 239 degrees of freedom
## Multiple R-squared:  0.5393, Adjusted R-squared:  0.5354 
## F-statistic: 139.9 on 2 and 239 DF,  p-value: < 2.2e-16

*以下、順に計算結果を見ていくと、

• 回帰係数 : Coefficients の Estimate の数値
  • 定数項、変量に対する係数

    Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
    (Intercept) -99.70771    9.92295 -10.048  < 2e-16 ***
    Height        0.80174    0.06103  13.136  < 2e-16 ***
    Chest         0.29231    0.04930   5.929 1.06e-08 ***

• 決定係数 : Adjusted R-squared
  • 目的変量が説明変量でどの程度説明しているかの割合。
  • 1 に近いほど当てはまりが良いと言える。

    Adjusted R-squared:  0.5354

• 説明変量が予測に役立っているかの検定
  • 回帰に役立っているか : p-value : 小さいと有意(役立っている)

    F-statistic: 139.9 on 2 and 239 DF,  p-value: < 2.2e-16

• 回帰係数がゼロではないかの検定
  • 回帰係数の検定(係数=0 か?) : Pr(>|t|) : 小さいと有意(係数=0ではないと言える)

    Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
    (Intercept) -99.70771    9.92295 -10.048  < 2e-16 ***
    Height        0.80174    0.06103  13.136  < 2e-16 ***
    Chest         0.29231    0.04930   5.929 1.06e-08 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7.5. 残差分析もお忘れなく

• 回帰直線の係数が求まったらそれでおしまい?

• 仮定が3つあったはず
  • 両者の関係性は直線＝線形を仮定する
  • 誤差の分布は正規分布を仮定する
  • 測定の場所によらず、残差は同一の分布であることを仮定する

• 仮定が成り立っているかは確認せねば。どうやって?
  • もし、残差が何らかの傾向を示しているなら、直線＝線形関係ではないのかもしれないと疑う必要がある。 ==> バナナカーブ
  • 測定する場所で残差の傾向が異なるなら、残差の分布が同一とは言えないかもしれない。 ==> 順に広がっているラッパカーブ
  • 残差のヒストグラム、箱ひげ図で偏りが判定できる。
  • Q-Q プロットの斜め直線に乗っていれば正規分布と言える。ズレが大きいようだと検討の余地がある。
    • 程度問題ではあるが、少しぐらいの逸脱は許容範囲と言える。
    • 残差が正規分布をするということを意識してもらえば良い。

  # 残差プロット
  plot(StudHWC$Height, RresultHWC$residuals)  # 対身長
  abline(v=seq(140,190,10), lty=3)     # 点線を追記
  abline(h=seq(-10,30,10), lty=3)     # 点線を追記
  abline(h=0, lty=1)                  # ゼロを実線で追記

  

  plot(StudHWC$Chest, RresultHWC$residuals)  # 対胸囲
  abline(v=seq(30,120,10), lty=3)      # 点線を追記
  abline(h=seq(-10,30,10), lty=3)     # 点線を追記
  abline(h=0, lty=1)                  # ゼロを実線で追記

  

  plot(StudHWC$Weight, RresultHWC$residuals)  # 対体重
  abline(v=seq(30,100,10), lty=3)      # 点線を追記
  abline(h=seq(-10,30,10), lty=3)     # 点線を追記
  abline(h=0, lty=1)                  # ゼロを実線で追記

  

  残差の正規性を確かめてみよう。

  • 残差のヒストグラム、箱ひげ図で偏りが判定できる。
  • 残差のQQプロット

  # 残差のヒストグラム
  hist(RresultHWC$residuals, right=F, ylim=c(0,100))
  abline(h=seq(0,100,20), lty=3)     # 点線を追記

  

  # 残差の箱ひげ図
  boxplot(RresultHWC$residuals, horizontal=T)
  abline(v=seq(-10,30,10), lty=3)      # 点線を追記

  

  # 残差のQQプロット: 正規確率プロット
  qqnorm(RresultHWC$residuals)
  qqline(RresultHWC$residuals, lty=2)
  abline(h=seq(-10,30,10), lty=3)   # 点線を追記
  abline(h=0, lty=1)                # 実線を追記
  abline(v=0, lty=1)               # 実線を追記

  

  • Q-Q プロットの斜め直線に乗っていれば正規分布と言える。ズレが大きいようだと検討の余地がある。
  • 程度問題ではあるが、少しぐらいの逸脱は許容範囲と言える。
  • 残差が正規分布をするということを意識してもらえば良い。

  • [検討事項] 残差が20を超える辺りの5サンプルほどは(少なくとも)分布を乱しているように見える。===> 外れ値か? 吟味が必要。

7.6. この分析のまとめ(結果の見方)

• 対象になったのは 242名。

• 説明変量が予測に役立っているか?
  • 回帰に役立っているか : 「F-statistic」中の「p-value」 : 小さいと有意(役立っている)
    [この例] 1% 未満(0.01%) なので役に立っていると言える。

• 決定係数 : Adjusted R-squared (相関係数 : R)
  • 目的変量が説明変量でどの程度説明しているかの割合。
  • 1 に近いほど当てはまりが良いと言える。
    [この例] 53.5% 程(約半分)を説明できている。

• 回帰係数 : 「Coefficients」 の 「Estimate」
  [この例] a=-99.7(定数項), b=0.802(身長), c=0.292(胸囲)

• 説明変数が予測に役立っているか?
  回帰係数の検定(係数=0 か?) : 「Coefficients」 の 「Pr(>|t|)」 : 小さいと有意(係数=0ではないと言える)
  [この例] 3つの係数とも 1% 未満(0.01%) なので回帰係数はゼロではない(何らかの意味がある数字と言える)。

• 残差の性質 ===> 正規性 : 残差プロット、残差解析
  • 残差(予測誤差)は正規分布をしていると仮定してモデルが構築されている。
  • この仮定が覆ると、回帰分析として成立していないことになる。
  • 残差が正規分布をしているか確認する必要がある。
  • 均等に散らばっているか?
  • 傾向はないか?
  • …
    [この例] 残差には概ね傾向は見られない。 ただし残差が20を超える辺りの5サンプルほどは(少なくとも)分布を乱しているように見える。加えて、今頃になってだが、胸囲が50cm未満の者が4サンプル程度ある。===> 外れ値か? 吟味が必要。

# 計算結果要約(再掲)
summary(RresultHWC)

## 
## Call:
## lm(formula = Weight ~ Height + Chest, data = StudHWC)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.019  -4.341  -1.411   3.280  31.561 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -99.70771    9.92295 -10.048  < 2e-16 ***
## Height        0.80174    0.06103  13.136  < 2e-16 ***
## Chest         0.29231    0.04930   5.929 1.06e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.038 on 239 degrees of freedom
## Multiple R-squared:  0.5393, Adjusted R-squared:  0.5354 
## F-statistic: 139.9 on 2 and 239 DF,  p-value: < 2.2e-16

8. 外れ値への対応: 集団と異なったふるまいをするサンプル

• 吟味する: なぜ外れている? 入力ミス?
• どう対応する: 外す、修正する、数理モデルが不適当?

• 一部のサンプルを除外して分析してみよう
  残差分析の結果、残差が20を超えている5サンプル程について吟味したところ、 一般的な学生の体型と異なっていると判断し、 これらを除外して分析をしてみよう。

• 検討の結果、以下の処置を取ることとした。
  • 体重が85Kgを超える者を除外
  • 胸囲が50cmに満たない者を除外

8.1. 該当サンプルの除外

元のExcelファイルに戻って除外する手もあるが、Rの中で実現するには、 不等号や論理演算子、否定(!, Not)等をうまく組み合わせて、対象とするサンプルを除外する。 それぞれの式の意味している事項が理解できるであろうか? もし理解できないようなら、式の中を中から順に分割して出力してみれば理解の促進になるであろう。

# 外れ値と考えたサンプルを抽出
StudTaiOver85<-StudHWC$Weight>85
StudHWC[StudTaiOver85,]                # 体重が85超えのサンプル一覧

##     Height Weight Chest
## 366  169.3   88.5    94
## 539  173.0   90.0    90
## 628  176.9   90.1    98
## 670  178.0  100.0   112
## 684  179.0   87.0   102
## 735  182.0   90.0   100
## 741  183.0   87.0   100

StudKyoLess50<-StudHWC$Chest<50
StudHWC[StudKyoLess50,]                # 胸囲が50未満のサンプル一覧

##     Height Weight Chest
## 26     153     57  38.0
## 355    169     60  30.5
## 538    173     84  46.0
## 599    175     72  45.0

StudHWC[StudTaiOver85|StudKyoLess50,]  # 体重が85超え または 胸囲が50未満のサンプル一覧

##     Height Weight Chest
## 26   153.0   57.0  38.0
## 355  169.0   60.0  30.5
## 366  169.3   88.5  94.0
## 538  173.0   84.0  46.0
## 539  173.0   90.0  90.0
## 599  175.0   72.0  45.0
## 628  176.9   90.1  98.0
## 670  178.0  100.0 112.0
## 684  179.0   87.0 102.0
## 735  182.0   90.0 100.0
## 741  183.0   87.0 100.0

# 一旦変数に置き直すと利用が簡単になる(理解がしやすくなる、実現していることは上記と同じ)
StudOutlier<-StudTaiOver85|StudKyoLess50
StudHWC[StudOutlier,]                 # 体重が85超え または 胸囲が50未満のサンプル

##     Height Weight Chest
## 26   153.0   57.0  38.0
## 355  169.0   60.0  30.5
## 366  169.3   88.5  94.0
## 538  173.0   84.0  46.0
## 539  173.0   90.0  90.0
## 599  175.0   72.0  45.0
## 628  176.9   90.1  98.0
## 670  178.0  100.0 112.0
## 684  179.0   87.0 102.0
## 735  182.0   90.0 100.0
## 741  183.0   87.0 100.0

「体重が85超え または 胸囲が60未満」の11サンプルを特定。 元々は242サンプルあったが、11サンプルを除外したので231サンプルになっている。 以下では変数名をStudHWC2とした。

StudHWC2<-StudHWC[!StudOutlier,]
StudHWC2[1:10,]

##    Height Weight Chest
## 2   145.5   42.0    76
## 3   146.7   41.0    85
## 5   148.0   43.0    80
## 8   150.0   43.0    82
## 9   150.0   46.0    86
## 16  151.7   41.5    80
## 17  152.0   35.0    77
## 23  153.0   46.5    87
## 25  153.0   55.0    78
## 32  154.4   44.0    75

dim(StudHWC2)

## [1] 231   3

dim(StudHWC)

## [1] 242   3

8.2. 改めて重回帰分析

では、抽出したデータ(228サンプル)「StudHWC2」を対象に、改めて重回帰分析を行ってみよう。

dim(StudHWC2)  # サンプルサイズを改めて確認

## [1] 231   3

# 重回帰分析
RresultHWC2<-lm(Weight~Height+Chest, data=StudHWC2)

# 計算結果要約
summary(RresultHWC2)

## 
## Call:
## lm(formula = Weight ~ Height + Chest, data = StudHWC2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.911  -3.404  -0.978   3.232  17.344 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -98.63620    7.71208  -12.79   <2e-16 ***
## Height        0.66740    0.04701   14.20   <2e-16 ***
## Chest         0.53132    0.04993   10.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.2 on 228 degrees of freedom
## Multiple R-squared:  0.6698, Adjusted R-squared:  0.6669 
## F-statistic: 231.2 on 2 and 228 DF,  p-value: < 2.2e-16

*以下、順に計算結果を見ていくと、

• 回帰係数 : Coefficients の Estimate の数値
  • 定数項、変量に対する係数

    Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
    (Intercept) -98.63620    7.71208  -12.79   <2e-16 ***
    Height        0.66740    0.04701   14.20   <2e-16 ***
    Chest         0.53132    0.04993   10.64   <2e-16 ***

• 決定係数 : Adjusted R-squared
  • 目的変量が説明変量でどの程度説明しているかの割合。
  • 1 に近いほど当てはまりが良いと言える。

    Adjusted R-squared:  0.6669

• 説明変量が予測に役立っているかの検定
  • 回帰に役立っているか : p-value : 小さいと有意(役立っている)

    F-statistic: 231.2 on 2 and 228 DF,  p-value: < 2.2e-16

• 回帰係数がゼロではないかの検定
  • 回帰係数の検定(係数=0 か?) : Pr(>|t|) : 小さいと有意(係数=0ではないと言える)

    Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
    (Intercept) -98.63620    7.71208  -12.79   <2e-16 ***
    Height        0.66740    0.04701   14.20   <2e-16 ***
    Chest         0.53132    0.04993   10.64   <2e-16 ***
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.3. 改めて残差分析もお忘れなく

# 残差プロット
plot(StudHWC2$Height, RresultHWC2$residuals)  # 対身長
abline(v=seq(140,190,10), lty=3)    # 点線を追記
abline(h=seq(-10,30,5), lty=3)      # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

plot(StudHWC2$Chest, RresultHWC2$residuals)  # 対胸囲
abline(v=seq(30,120,10), lty=3)     # 点線を追記
abline(h=seq(-10,30,5), lty=3)      # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

plot(StudHWC2$Weight, RresultHWC2$residuals)  # 対体重
abline(v=seq(30,100,10), lty=3)     # 点線を追記
abline(h=seq(-10,30,5), lty=3)      # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

残差の正規性を確かめてみよう。

• 残差のヒストグラム、箱ひげ図で偏りが判定できる。
• 残差のQQプロット

# 残差のヒストグラム
hist(RresultHWC2$residuals, right=F, ylim=c(0,100))
abline(h=seq(0,100,20), lty=3)     # 点線を追記

# 残差の箱ひげ図
boxplot(RresultHWC2$residuals, horizontal=T)
abline(v=seq(-10,15,5), lty=3)      # 点線を追記

• 残差のヒストグラム、箱ひげ図で偏りが判定できる。

# 残差のQQプロット: 正規確率プロット
qqnorm(RresultHWC2$residuals)
qqline(RresultHWC2$residuals, lty=2)
abline(h=seq(-10,30,5), lty=3)   # 点線を追記
abline(h=0, lty=1)               # 実線を追記
abline(v=0, lty=1)               # 実線を追記

外れ値を除外する前後で残差のQQプロットが改善されたことを確認せよ。

8.4. 改めての分析のまとめ(結果の見方)

• 対象になったのは 231名。
  [この例] 一般的な学生の体格と異なる11サンプルを除外した。

• 説明変量が予測に役立っているか?
  • 回帰に役立っているか : 「F-statistic」中の「p-value」 : 小さいと有意(役立っている)
    [この例] 1% 未満(0.01%) なので役に立っていると言える。

• 決定係数 : Adjusted R-squared (相関係数 : R)
  • 目的変量が説明変量でどの程度説明しているかの割合。
  • 1 に近いほど当てはまりが良いと言える。
    [この例] 66.7% (3分の2)を説明できている。単回帰分析のときよりもより多く説明できている。

• 回帰係数 : 「Coefficients」 の 「Estimate」
  [この例] a=-98.6(定数項), b=0.667(身長), c=0.531(胸囲)

• 説明変数が予測に役立っているか?
  回帰係数の検定(係数=0 か?) : 「Coefficients」 の 「Pr(>|t|)」 : 小さいと有意(係数=0ではないと言える)
  [この例] 3つとも 1% 未満(0.01%) なので回帰係数はゼロではない(何らかの意味がある数字と言える)。

• 残差の性質 ===> 正規性 : 残差プロット、残差解析
  • 残差(予測誤差)は正規分布をしていると仮定してモデルが構築されている。
  • この仮定が覆ると、回帰分析として成立していないことになる。
  • 残差が正規分布をしているか確認する必要がある。
  • 均等に散らばっているか?
  • 傾向はないか?
  • …
    [この例] 残差には概ね傾向は見られない。QQプロットもほぼ斜めライン上に絡まるように推移している。

# 計算結果要約(再掲)
summary(RresultHWC2)

## 
## Call:
## lm(formula = Weight ~ Height + Chest, data = StudHWC2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.911  -3.404  -0.978   3.232  17.344 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -98.63620    7.71208  -12.79   <2e-16 ***
## Height        0.66740    0.04701   14.20   <2e-16 ***
## Chest         0.53132    0.04993   10.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.2 on 228 degrees of freedom
## Multiple R-squared:  0.6698, Adjusted R-squared:  0.6669 
## F-statistic: 231.2 on 2 and 228 DF,  p-value: < 2.2e-16

[参考: 出力の読み方の詳細] 統計検定のための、Rの出力結果からわかること（回帰分析編）

[演習8.4-1] 上述の重回帰分析について男性のみ、女性のみに分離して分析してみよ。

[演習8.4-2] 先週の単回帰分析について残差の大きいサンプルを除外して分析してみよ。

【今週の思い】回帰分析を成立させる背景やアイディアを理解し、その仮定を逸脱しないように分析を進める手順を習得せよ。

9. 4つの尺度と回帰分析

• 比率尺度(ratio scale) : 長さ、重さ。間隔尺度に加えて絶対原点が存在。

• 間隔尺度(interval scale) : セ氏やカ氏で測定された温度。原点は任意。

• 名義尺度(nominal scale) : 性別、職業、支持政党。分類の項目。
• 順序尺度(ordinal scale) : 成績の優・良・可・不可。順序関係あり。「優と良の差」が「可と不可の差」に等しいわけではない。

• 量的データ : 間隔尺度、比率尺度。連続変量とも。
• 質的データ : 名義尺度、順序尺度。離散変量とも。

• 離散変量を回帰分析で用いるには注意が必要である。 線形関係と言えるのか? 割り当てた値を変更したら影響力が変化する。

10. [要約: 回帰分析] 解析する上での注意点

• 2変量以上の説明する変量(説明変量)で 1変量(目的変量)を説明したい
  • 説明される変量 : 目的変数、従属変数、dependent variable
  • 説明する変量 : 説明変数、独立変数、independent variable

• 関係性は直線＝線形を仮定する
• 回帰直線からのズレ(誤差)のことを回帰分析では「残差」と呼ぶ。
• 説明される変量(目的変量)と平行に残差を取る。なぜなら、指定された独立変数における測定誤差だから。
• 残差の分布は正規分布を仮定する
• 測定の場所によらず、残差は同一の分布であることを仮定する

• 残差の二乗和を最小にする(最小二乗法)

• 残差分析を行なおう
  • 残差に傾向はないか? :
    　　　一様性、バナナカーブ、ラッパ状、… ===> 残差プロット
  • 残差の正規性の確認 : Q-Q plot

• 外れ値に対して
  • 外れ値と思われるサンプルを検討する必要あり。場合によっては除外も。
    • 外れ値を持つ個体を除外する : 標準体型の者だけに、…
    • 対象集団を単質のものに絞る : 男性だけ、18歳だけ、…
    • ただし、恣意的な除外は謹むべきである。

  • 「外れ値」と「異常値」は区別して取り扱う・識別する : 「外れて」いても「異常」かどうかは別
    • 外れ値 : 大多数のデータから大きく離れた値
    • 異常値 : 通常ではあり得ない入力ミス等の異常な値

• また、 変量数が多いデータを対象とした場合は、 より影響力の大きい変量だけでモデルを構築する必要に迫られることも多い。 そのような場合は「変数選択」の手法を用いると良い。 一方、変量の組み合わせによって相関係数がどのように変化するかを観るために 全組み合わせを順に発生させて振る舞いを眺める「総当たり法」を用いることもある。 これらは実データを分析する際に試して欲しい。

• 一般的には説明変量を多くするほど説明力 \(\textbf{R}^2\) は上がるが、 「数値計算上の最適モデル」と「その分野の知識からの最適モデル」には違いがあることを知っておくこと。
  • 計算結果としてベターであっても、回帰式・モデルを理解する上で好ましいとは言えない場合も有り得る。

• 上記以外に「多重共線性(multi-colinearity)」の問題もある。 逆行列が取れなくなる(ランク落ち)。 数値計算上、解が求まったかのような振る舞いをするものの不安定な解。 そのような場合は、組み込んでいる説明変数間の関係を見直す必要がある。

  
  \[ {\textbf{X}}= \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1k} & x_{1l} & x_{1m} & \cdots & x_{1p}\\ 1 & x_{21} & x_{23} & \cdots & x_{2k} & x_{2l} & x_{2m} & \cdots & x_{2p} \\ 1 & x_{31} & x_{32} & \cdots & x_{3k} & x_{3l} & x_{3m} & \cdots & x_{3p} \\ \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{i1} & x_{i2} & \cdots & x_{ik} & x_{il} & x_{im} & \cdots & x_{ip} \\ \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nk} & x_{nl} & x_{nm} & \cdots & x_{np} \\ \end{pmatrix} \]
  

  ここで、例えば

  \[ a) \; x_{im} = x_{ik} \;\;\;\;\;\; (i=1,2,...,n) \\ \]

  \[ b) \; x_{im} = x_{ik} + x_{il} \;\;\;\;\;\; (i=1,2,...,n) \\ \]

  \[ c) \; x_{im} = \alpha \; x_{ik} \;\;\;\;\;\; \alpha は実数, \; (i=1,2,...,n) \\ \]

  等だったら、\(\textbf{X}\)はランク落ちする。 すなわち、逆行列が取れなくなり、多重共線性(multi-colinearity)が起こる。 また、既に投入した変量と相関が極めて高い変量を投入した場合も同様の事態が起きる。 このような場合は、投入した変量を見直す必要がある。

11. 回帰分析にまつわる話題、誤用?!: 【提示資料】

• 構造がシンプルで理解し易い分、“変な”利用も散見される
• 内挿(観測点の内側)はまだしも、外挿(観測点の外側)を予測するのは難しい or 無謀。
• 予測範囲でデータの構造が一定という”条件”が必要。
• 線形で表現できる関係なのか? 線形の関係式が有効なのか? 非線型?
• 自分で判断できる能力を持つこと。疑うこと。

• [例1] 人間の成長曲線 (構造変化)
• [例2] 将来のプログラマ必要数予測 (外挿) : 21世紀(?)には国民全員がプログラマ (’80s)
• [例3] 将来予測：直線回帰 (外挿)
  • オリンピック 100m 走の男女記録 : 2156年には女性の方が速い (2004.09.30) :
    Japan Journal LTD の記事 , 朝日新聞 の記事
    　　　[究極の命題!] 100m に 0.00秒 要する(!?)ようになるのは何時?

12. 【余談】回帰分析と最小二乗法

• 測定誤差が正規分布していることを仮定して、回帰分析はモデルが構築されている。
• 最小二乗法(Least Squares Method)は予測変量の誤差を最小とする数値計算手法であり、求める線や平面に構造を持ち込んでいるわけでなない。
• 最近話題のAIでも回帰分析に基づいた推定等が行われているが、 その際に、「誤差が正規分布をする」との仮定は設けられていないことがある。 その場合は「回帰分析」と言わずに「最小二乗法」と言った方が正確ではないのか? 【私見】

13. 最終レポート

【最終レポート(後期前半)】
　以下の事項について、WordファイルもしくはPDFファイルでレポートを作成し、Moodleから提出下さい。なお、提出フォームは11月14日(木) 昼から利用可能となります。

1.【必須項目】
　かねてからお伝えしていたように、各自が興味を持って収集したデータに対して統計手法を適用し、何らかの知見を引き出して報告してください。統計手法として、本講義で紹介した統計手法には限定しません。また、対象とするデータは1つ以上いくつでもかまいません。

◎記載内容: 以下に挙げるような項目を含めて作成すること。

1. 所属学部名、学籍番号、氏名
2. データ内容の説明
3. どのような点に興味を持ったか
4. 自分の解析目的
5. 何を知りたいためにどのような手法を使ったのか
6. 得られた知見と考察
7. その他、気付いたこと

2.【必須項目】
　本講義を受講することによって「統計」に抱くイメージが変化したかを述べよ。変化した場合 or しない場合の各々で、講義を聞き終えた現状でどのように感じているか、また今後自分として統計に対してどのように取り組みたい/取り組みたくないかを説明せよ。

3.【任意項目(コメントがあれば嬉しいな)】 講義方法、講義の進め方
　講義内容の感想だけでなく、講義方法等に付いて気になった点や感想、改善希望点をお聞かせください。

【期日】11月21日(木) 朝まで (それ以降は受け付けなくなるので要注意)
【注意】単位取得を考えている者は期日までに提出すること。

82. 参考

このページで取り扱ったプログラムだけを抜き出して以下に列挙しておく。

# ディレクトリの移動。必須ではない。個々人の設定に応じて。
setwd("D:/home_sub3/R_Dir")  # ホームディレクトリに移動(Set Working Directory)
getwd()           # 現在のディレクトリ位置を表示
list.files()      # ファイル・ディレクトリ一覧を表示
setwd("KougiDS24")  # ディレクトリを移動
list.files()      # ファイル・ディレクトリ一覧を表示

# データの読み込み
Student24<-read.csv("StudAll24b.csv", skip=5,
                     header=TRUE, na.strings="NA")
dim(Student24)
colnames(Student24)
Student24[1:5,]

## 7.2. 前準備: 完全データの抽出
# データ切り出し
dim(Student24)
StudHWC<-na.omit(Student24[,3:5]) # 欠損値を含むデータを削除
dim(StudHWC)
colnames(StudHWC)
StudHWC[1:10,]

## 7.3. 散布図行列
# 散布図行列
plot(StudHWC)

## 7.4. プログラム ＆ 計算結果
dim(StudHWC)  # サンプルサイズを改めて確認
# 重回帰分析
RresultHWC<-lm(Weight~Height+Chest, 
               data=StudHWC) # 右辺に説明する変量をプラスマークで追加していく

# 計算結果要約
summary(RresultHWC)

## 7.5. 残差分析もお忘れなく
# 残差プロット
plot(StudHWC$Height, RresultHWC$residuals)  # 対身長
abline(v=seq(140,190,10), lty=3)    # 点線を追記
abline(h=seq(-10,30,10), lty=3)     # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

plot(StudHWC$Chest, RresultHWC$residuals)  # 対胸囲
abline(v=seq(30,120,10), lty=3)     # 点線を追記
abline(h=seq(-10,30,10), lty=3)     # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

plot(StudHWC$Weight, RresultHWC$residuals)  # 対体重
abline(v=seq(30,100,10), lty=3)     # 点線を追記
abline(h=seq(-10,30,10), lty=3)     # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

# 残差のヒストグラム
hist(RresultHWC$residuals, right=F, ylim=c(0,100))
abline(h=seq(0,100,20), lty=3)     # 点線を追記

# 残差の箱ひげ図
boxplot(RresultHWC$residuals, horizontal=T)
abline(v=seq(-10,30,10), lty=3)   # 点線を追記

# 残差のQQプロット: 正規確率プロット
qqnorm(RresultHWC$residuals)
qqline(RresultHWC$residuals, lty=2)
abline(h=seq(-10,30,10), lty=3)   # 点線を追記
abline(h=0, lty=1)                # 実線を追記
abline(v=0, lty=1)               # 実線を追記

## 7.6. この分析のまとめ(結果の見方)
# 計算結果要約(再掲)
summary(RresultHWC)

## 8.1. 該当サンプルの除外
# 外れ値と考えたサンプルを抽出
StudTaiOver85<-StudHWC$Weight>85
StudHWC[StudTaiOver85,]                # 体重が85超えのサンプル一覧

StudKyoLess50<-StudHWC$Chest<50
StudHWC[StudKyoLess50,]                # 胸囲が50未満のサンプル一覧

StudHWC[StudTaiOver85|StudKyoLess50,]  # 体重が85超え または 胸囲が50未満のサンプル一覧

# 一旦変数に置き直すと利用が簡単になる(理解がしやすくなる、実現していることは上記と同じ)
StudOutlier<-StudTaiOver85|StudKyoLess50
StudHWC[StudOutlier,]                 # 体重が85超え または 胸囲が50未満のサンプル
StudHWC2<-StudHWC[!StudOutlier,]
StudHWC2[1:10,]
dim(StudHWC2)
dim(StudHWC)

## 8.2. 改めて重回帰分析
dim(StudHWC2)  # サンプルサイズを改めて確認

# 重回帰分析
RresultHWC2<-lm(Weight~Height+Chest, data=StudHWC2)

# 計算結果要約
summary(RresultHWC2)

## 8.3. 改めて残差分析もお忘れなく
# 残差プロット
plot(StudHWC2$Height, RresultHWC2$residuals)  # 対身長
abline(v=seq(140,190,10), lty=3)    # 点線を追記
abline(h=seq(-10,30,5), lty=3)      # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

plot(StudHWC2$Chest, RresultHWC2$residuals)  # 対胸囲
abline(v=seq(30,120,10), lty=3)     # 点線を追記
abline(h=seq(-10,30,5), lty=3)      # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

plot(StudHWC2$Weight, RresultHWC2$residuals)  # 対体重
abline(v=seq(30,100,10), lty=3)     # 点線を追記
abline(h=seq(-10,30,5), lty=3)      # 点線を追記
abline(h=0, lty=1)                  # ゼロを実線で追記

# 残差のヒストグラム
hist(RresultHWC2$residuals, right=F, ylim=c(0,100))
abline(h=seq(0,100,20), lty=3)     # 点線を追記

# 残差の箱ひげ図
boxplot(RresultHWC2$residuals, horizontal=T)
abline(v=seq(-10,15,5), lty=3)      # 点線を追記

# 残差のQQプロット: 正規確率プロット
qqnorm(RresultHWC2$residuals)
qqline(RresultHWC2$residuals, lty=2)
abline(h=seq(-10,30,10), lty=3)   # 点線を追記
abline(h=0, lty=1)                # 実線を追記
abline(v=0, lty=1)               # 実線を追記

## 8.4. 改めての分析のまとめ(結果の見方)
# 計算結果要約(再掲)
summary(RresultHWC2)