『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第8章　単回帰分析(SRA)

Author

SUGINO Isamu, Build with R4.3.1

Published

November 30, 2024

0　全体の章構成

1　単回帰分析の基礎

\[Y=\beta_{0}+\beta_{1}X_{1}+\epsilon,\;\;\;\;\;\;\epsilon\;\sim\;N(0,\sigma_{\epsilon}^2)\] \[\color{blue}y_{i}\color{black}=\color{green}\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}\color{black}+\color{red}e_{i}\color{black}=\color{green}\hat{y}_{i}\color{black}+\color{red}e_{i}\color{black}\] 青が実測値・観測値で，緑の予測値・予測式と赤の残差の和からなると考える。

1-1　最小二乗和法(Least Squares Method)

残差平方和
　ここでは式変形を見易くする為に，$x_{1i}$と書かずに単に$x_{i}$と書く。
\[\sum_{i=1}^{n}\color{red}e_{i}\color{black}^{2}=\sum_{i=1}^{n}(\color{blue}y_{i}\color{black}-\color{green}\hat{y_{i}}\color{black})^2=\sum_{i=1}^{n}(\color{blue}y_{i}\color{black}-\color{green}\hat{\beta_{0}}-\hat{\beta_{1}}x_{i}\color{black})^2\]

式変形による解

単回帰の最小二乗和推定量の計算は単純な式変形で導くことができる。
ポイントは，未知数についての 2 次式であると見て式変形を行うことである。
まずは式の中でより単純な$\hat{\beta_{0}}$の方について整理して平方完成を行い，残った部分を$\hat{\beta_{1}}$について平方完成する。
式変形の途中で，x や y の算術平均，標本分散，標本共分散の定義式を利用する。　 \[\begin{align*} \sum_{i=1}^{n}e_{i}^{2} &=\sum_{i=1}^{n}(-\hat{\beta_{0}}+y_{i}-\hat{\beta_{1}}x_{i})^2\\ &=\sum_{i=1}^{n}\left\{\hat{\beta_{0}}^2-2(y_{i}-\hat{\beta_{1}}x_{i})\hat{\beta_{0}}+(y_{i}-\hat{\beta_{1}}x_{i})^2\right\}\\ &=\sum_{i=1}^{n}\hat{\beta_{0}}^2-2\sum_{i=1}^{n}(y_{i}-\hat{\beta_{1}}x_{i})\hat{\beta_{0}}+\sum_{i=1}^{n}(y_{i}-\hat{\beta_{1}}x_{i})^2\\ &=n\hat{\beta_{0}}^2-2\left\{\sum_{i=1}^{n}y_{i}-\hat{\beta_{1}}\sum_{i=1}^{n}x_{i}\right\}\hat{\beta_{0}}+\sum_{i=1}^{n}(y_{i}-\hat{\beta_{1}}x_{i})^2\\ &=n\hat{\beta_{0}}^2-2n(\bar{y}-\hat{\beta_{1}}\bar{x})\hat{\beta_{0}}+\sum_{i=1}^{n}(y_{i}-\hat{\beta_{1}}x_{i})^2\\ &=n\left\{\hat{\beta_{0}}-(\bar{y}-\hat{\beta_{1}}\bar{x})\right\}^2+\sum_{i=1}^{n}(y_{i}-\hat{\beta_{1}}x_{i})^2-n(\bar{y}-\hat{\beta_{1}}\bar{x})^2 \end{align*}\]

$\hat{\beta_{0}}$の平方完成ができたので，上式の後半の部分を$\hat{\beta_{1}}$について平方完成する。

\[\begin{align*} \sum_{i=1}^{n}(y_{i}-\hat{\beta_{1}}x_{i})^2-n(\bar{y}-\hat{\beta_{1}}\bar{x})^2 &=\sum_{i=1}^{n}(y_{i}^2-2\hat{\beta_{1}}x_{i}y_{i}+\hat{\beta_{1}}^2x_{i}^2)-n\bar{y}^2+2n\hat{\beta_{1}}\bar{x}\bar{y}-n\hat{\beta_{1}}^2\bar{x}^2\\ &=\sum_{i=1}^{n}y_{i}^2-n\bar{y}^2-2\hat{\beta_{1}}\left(\sum_{i=1}^{n}x_{i}y_{i}-n\bar{x}\bar{y}\right)+\hat{\beta_{1}}^2\left(\sum_{i=1}^{n}x_{i}^2-n\bar{x}^2\right)\\ &=n\left(\frac{1}{n}\sum_{i=1}^{n}y_{i}^2-\bar{y}^2\right)-2n\hat{\beta_{1}}\left(\frac{1}{n}\sum_{i=1}^{n}x_{i}y_{i}-\bar{x}\bar{y}\right)+n\hat{\beta_{1}}^2\left(\frac{1}{n}\sum_{i=1}^{n}x_{i}^2-\bar{x}^2\right)\\ &=n(s_{y}^2-2s_{xy}\hat{\beta_{1}}+s_{x}^2\hat{\beta_{1}}^2)=n\left\{s_{x}^2\left(\hat{\beta_{1}}^2-2\frac{s_{xy}}{s_{x}^2}\hat{\beta_{1}}+s_{y}^2\right)\right\}\\ &=n\left\{s_{x}^2\left(\hat{\beta_{1}}-\frac{s_{xy}}{s_{x}^2}\right)^2+s_{y}^2-\left(\frac{s_{xy}}{s_{x}^2}\right)^2\right\} \end{align*}\]

$s_{y}^2=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\bar{y})^2=\frac{1}{n}\sum_{i=1}^{n}y_{i}^2-\bar{y}^2$などを利用するのがポイントである。

　以上より，$\hat{\beta_{0}}=\bar{y}-\hat{\beta_{1}}\bar{x}, \hat{\beta_{1}}=\frac{s_{xy}}{s_{x}^2}=r_{xy}\frac{s_{y}}{s_{x}}$の時に残差平方和は最小値をとる。
これで，標本データから予測式を求める事が出来る。 \[\hat{y_{i}}=\bar{y}-r_{xy}\frac{s_{y}}{s_{x}}\bar{x}+r_{xy}\frac{s_{y}}{s_{x}}x_{i}=r_{xy}\frac{s_{y}}{s_{x}}(x_{i}-\bar{x})+\bar{y}\] 　これは平面における直線の式なので回帰直線とも呼ぶ。
式から分かる様に，この直線は$(\bar{x},\bar{y})$を通る。
また，xとyの標準偏差が等しい時，回帰係数は積率相関係数に一致する。

残差の平均，残差と独立変数の共分散，残差と予測値の共分散

　次に，残差の平均や，残差と独立変数の共分散，残差と予測値の共分散が0になる事について。
　最小二乗和法で残差平方和を最小にした時，残差 e の平均，残差 e と x の共分散，残差 e と予測値の共分散はそれぞれ 0 になる。 $\hat{\beta_{0}}=\bar{y}-\hat{\beta_{1}}\bar{x}, \hat{\beta_{1}}=\frac{s_{xy}}{s_{x}^2}=r_{xy}\frac{s_{y}}{s_{x}}$であるから，

\[\begin{align*} \bar{e} &=\frac{1}{n}\sum_{i=1}^{n}e_{i}=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})=\frac{1}{n}\left(\sum_{i=1}^{n}y_{i}-\sum_{i=1}^{n}\hat{\beta_{0}}-\sum_{i=1}^{n}\hat{\beta_{1}}x_{i}\right)\\ &=\bar{y}-\hat{\beta_{0}}-\hat{\beta_{1}}\bar{x}=\bar{y}-(\bar{y}-\hat{\beta_{1}}\bar{x})-\hat{\beta_{1}}\bar{x}=0\\ s_{xe} &=\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})(e_{i}-\bar{e})=\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})(e_{i}-0)=\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})\\ &=\frac{1}{n}\sum_{i=1}^{n}x_{i}(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})+\bar{x}\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})=\frac{1}{n}\sum_{i=1}^{n}x_{i}(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})+\bar{x}\cdot0\\ &=\frac{1}{n}\sum_{i=1}^{n}(x_{i}y_{i}-\hat{\beta_{0}}x_{i}-\hat{\beta_{1}}x_{i}^2)=\frac{1}{n}\sum_{i=1}^{n}x_{i}y_{i}-\hat{\beta_{0}}\frac{1}{n}\sum_{i=1}^{n}x_{i}-\hat{\beta_{1}}\frac{1}{n}\sum_{i=1}^{n}x_{i}^2\\ &=\frac{1}{n}\sum_{i=1}^{n}x_{i}y_{i}-(\bar{y}-\hat{\beta_{1}}\bar{x})\bar{x}-\hat{\beta_{1}}\frac{1}{n}\sum_{i=1}^{n}x_{i}^2=\frac{1}{n}\sum_{i=1}^{n}x_{i}y_{i}-\bar{x}\bar{y}+\hat{\beta_{1}}\bar{x}^2-\hat{\beta_{1}}\frac{1}{n}\sum_{i=1}^{n}x_{i}^2\\ &=\frac{1}{n}\sum_{i=1}^{n}x_{i}y_{i}-\bar{x}\bar{y}-\hat{\beta_{1}}\bar{x}^2-\hat{\beta_{1}}\frac{1}{n}\sum_{i=1}^{n}x_{i}^2=s_{xy}-\hat{\beta_{1}}\left(\frac{1}{n}\sum_{i=1}^{n}x_{i}^2-\bar{x}^2\right)=s_{xy}-\hat{\beta_{1}}s_{x}^2\\ &=s_{xy}-\frac{s_{xy}}{s_{x}^2}s_{x}^2=0\\ s_{\hat{y}e} &=\frac{1}{n}\sum_{i=1}^{n}(\hat{y_{i}}-\bar{\hat{y}})(e_{i}-\bar{e})=\frac{1}{n}\sum_{i=1}^{n}(\hat{y_{i}}-\bar{\hat{y}})e_{i}\\ &=\frac{1}{n}\sum_{i=1}^{n}(\hat{\beta_{0}}+\hat{\beta_{1}}x_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}\bar{x})(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})=\frac{1}{n}\sum_{i=1}^{n}\hat{\beta_{1}}(x_{i}-\bar{x})(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})\\ &=\hat{\beta_{1}}\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\hat{\beta_{0}}-\hat{\beta_{1}}x_{i})=\hat{\beta_{1}}\cdot0=0 \end{align*}\]

以上から，残差と$x$，残差と予測値$\hat{y}$の共分散ががゼロである事が分かる。共分散がゼロであれば積率相関係数も当然ゼロであり，統計的に独立（無相関）である事が示された。

1-2　回帰係数(regression coefficient)のt検定と区間推定(interval estimation)

\[ \frac{\color{red}\hat{\beta}_{1}\color{black}-\color{blue}\beta_{1}\color{black}}{\frac{\color{red}\hat{\sigma}_{\epsilon}\color{black}}{\sqrt{\color{red}SS_{x_{1}}\color{black}}}}\;\sim\;t_{(n-2)} \] 回帰係数のt検定を行う場合は，（殆どの場合は）$H_{0}: \color{blue}\beta_{1}\color{black}=\color{green}0\color{blue}$として，

\[\frac{\color{red}\hat{\beta}_{1}\color{black}}{\frac{\color{red}\hat{\sigma}_{\epsilon}\color{black}}{\sqrt{\color{red}SS_{x_{1}}\color{black}}}}\;\sim\;t_{(n-2)}\] から，検定統計量tが棄却限界値を超えるかどうかを見る。

母回帰係数の区間推定をする場合には，以下の式を$\color{blue}\beta_{1}\color{black}$を中心として変形する。

\[ t_{0.025(df=n-2)}<\frac{\color{red}\hat{\beta}_{1}\color{black}-\color{blue}\beta_{1}\color{black}}{\frac{\color{red}\hat{\sigma}_{\epsilon}\color{black}}{\sqrt{\color{red}SS_{x_{1}}\color{black}}}}<t_{0.975(df=n-2)} \] \[ \color{red}\hat{\beta}_{1}\color{black}-t_{0.975(df=n-2)}\frac{\color{red}\hat{\sigma}_{\epsilon}\color{black}}{\sqrt{\color{red}SS_{x_{1}}\color{black}}}<\color{blue}\beta_{1}\color{black}<\color{red}\hat{\beta}_{1}\color{black}-t_{0.025(df=n-2)}\frac{\color{red}\hat{\sigma}_{\epsilon}\color{black}}{\sqrt{\color{red}SS_{x_{1}}\color{black}}} \]

一例として$n=296$とすると，t分布の値はそれぞれ，$t_{0.025(df=n-2)}$=-1.9680657，$t_{0.975(df=n-2)}$=1.9680657となる。

　ここはデータフレイム名を簡潔に d01 としよう。世帯年収で幸福度を説明する単回帰モデルの分析結果は以下の通りである。

d01 <- read.csv("practice.csv")

まずは変数の分布を確認

with(d01, table(fincome, useNA = "always"))

fincome
   0   25   75  125  200  300  400  500  600  700  800  925 1125 1375 1750 2500 
   5    1    7    8   17   19   28   36   32   27   23   33   31   13   11    6 
<NA> 
  87

with(d01, table(q1700, useNA = "always"))

q1700
   0    1    2    3    4    5    6    7    8    9   10   11 <NA> 
   2    2    4   17   16   65   42   91   81   41   17    1    5

必要に応じて欠損値処理（方法は色々）

d01$q1700[d01$q1700 == 11] <- NA
with(d01, table(q1700, useNA = "always"))

q1700
   0    1    2    3    4    5    6    7    8    9   10 <NA> 
   2    2    4   17   16   65   42   91   81   41   17    6

lm( ) と summary( )

reg01 <- lm(q1700 ~ fincome, d01)
reg01


Call:
lm(formula = q1700 ~ fincome, data = d01)

Coefficients:
(Intercept)      fincome  
   5.908985     0.001229

sum01 <- summary(reg01)
sum01


Call:
lm(formula = q1700 ~ fincome, data = d01)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1547 -1.1855  0.2309  1.2309  3.9989 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.9089853  0.1917289  30.819  < 2e-16 ***
fincome     0.0012288  0.0002233   5.503 8.13e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.796 on 294 degrees of freedom
  (88 observations deleted due to missingness)
Multiple R-squared:  0.09339,   Adjusted R-squared:  0.09031 
F-statistic: 30.28 on 1 and 294 DF,  p-value: 8.127e-08

オブジェクトに保存されている情報

　それぞれのオブジェクトにどの様な情報が格納されているのか名前を確認しておく。

names(reg01) # lm( )関数の結果

 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "na.action"     "xlevels"       "call"          "terms"        
[13] "model"

names(sum01) # summary( )関数の結果

 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"

reg01$coefficients

(Intercept)     fincome 
5.908985346 0.001228757

sum01$coefficients

               Estimate  Std. Error   t value     Pr(>|t|)
(Intercept) 5.908985346 0.191728906 30.819481 4.386859e-94
fincome     0.001228757 0.000223283  5.503138 8.126713e-08

sum01$df

[1]   2 294   2

定数項と回帰係数のt検定

定数項と回帰係数のt検定の両側有意確率を求めておく。

# オブジェクトsum01に含まれる情報からt統計量を計算し，t検定を行う。
t0 <- sum01$coefficients[,1]/sum01$coefficients[,2]; t0

(Intercept)     fincome 
  30.819481    5.503138

2 * pt(t0, df = sum01$df[2], lower.tail = F) # 両側t検定の有意確率

 (Intercept)      fincome 
4.386859e-94 8.126713e-08

先のsummary(lm( ))と結果を見比べて，正しく計算出来ている事が確認出来る。

信頼区間の求め方

　定数項や回帰係数の信頼区間の求め方を示す。
最初は定義に従って自分で式を組み立てて求める方法，次にconfint( )関数によって自動的に求める方法である。

# 自由度と，t分布の限界値を確認
sum01$df[2]; qt(.975, df=sum01$df[2])

[1] 294

[1] 1.968066

# 95%信頼区間下限
sum01$coefficients[,1] - qt(.975, df = sum01$df[2]) * sum01$coefficients[,2]

 (Intercept)      fincome 
5.5316502652 0.0007893216

# 95%信頼区間上限
sum01$coefficients[,1] + qt(.975, df = sum01$df[2]) * sum01$coefficients[,2]

(Intercept)     fincome 
6.286320427 0.001668193

# 信頼区間を求める関数
confint(reg01, level = .95)

                   2.5 %      97.5 %
(Intercept) 5.5316502652 6.286320427
fincome     0.0007893216 0.001668193

対応する箇所をよく見比べて欲しい。
　母回帰係数の95%信頼区間が.00079から.00167であると云う事は，有意水準5%で母回帰係数のt検定を行った場合，ゼロ仮説「母回帰係数は.00079である」から「母回帰係数は.00167である」までのゼロ仮説はどれも棄却されないと云う事でもある。
区間に0を含んでいないので，「母回帰係数は0である」は棄却される。
これは先の回帰係数のt検定の結果の通りである。

1-3　回帰モデルのF検定

分散分析表の作成

　分散分析表を作成する為に，anova( )関数を使用する。

# lm( )の結果であるreg01をanova( )に与える。
ano01 <- anova(reg01); ano01

names(ano01) # 格納されている情報を確認

[1] "Df"      "Sum Sq"  "Mean Sq" "F value" "Pr(>F)"

ano01$"Sum Sq" # 偏差平方和だけを取り出す

[1]  97.72471 948.70435

sum(ano01$"Sum Sq") # 総平方和を計算する

[1] 1046.429

総平方和が，従属変数の偏差平方和である事を確認する。
独立変数か従属変数にNAを含むケースを全て除外して完備ケースだけにしないと結果が一致しなくなるので注意。
完備ケースだけのデータフレイム d02 を作成する。

d02 <- d01[complete.cases(d01$fincome, d01$q1700), c("fincome", "q1700")]
sum((d02$q1700 - mean(d02$q1700))^2)

[1] 1046.429

総平方和が従属変数の偏差平方和である事が確認出来た。

単回帰分析の分散分析表
	平方和SS	自由度df	平均平方MS	F値	p値
回帰モデル	97.72	1	97.725	30.285	8.126713e-08
残差	948.7	294	3.227
全体	1046.43	295

回帰係数のt検定とモデルのF検定の同値性

単回帰分析においては，回帰係数のt検定とモデルのF検定は同値である事を示す。

# 回帰係数のt検定の結果の一部
names(sum01)

 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"

sum01$coefficients

               Estimate  Std. Error   t value     Pr(>|t|)
(Intercept) 5.908985346 0.191728906 30.819481 4.386859e-94
fincome     0.001228757 0.000223283  5.503138 8.126713e-08

(t0 <- sum01$coefficients[2, 3]) # t統計量

[1] 5.503138

t0^2 # t統計量の二乗

[1] 30.28453

# モデルのF検定の結果の一部
names(ano01)

[1] "Df"      "Sum Sq"  "Mean Sq" "F value" "Pr(>F)"

ano01$"F value"[1] # F統計量

[1] 30.28453

sum01$coefficients[2, 4] # t検定のp値

[1] 8.126713e-08

ano01$"Pr(>F)"[1] # F検定のp値

[1] 8.126713e-08

1-4　決定係数(coefficient of determinance)，分散説明率(proportion of variance explained)，誤差減少率(PRE)

　決定係数（分散説明率）を幾通りかの方法で確認してみよう。

平方和の比

ano01$"Sum Sq"[1] # モデル平方和

[1] 97.72471

sum(ano01$"Sum Sq") # 総平方和

[1] 1046.429

ano01$"Sum Sq"[1]/sum(ano01$"Sum Sq") # 分散説明率

[1] 0.09338876

回帰分析の結果のR^2

sum01 # Multiple R-squared と Adjusted R-squared に注目


Call:
lm(formula = q1700 ~ fincome, data = d01)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1547 -1.1855  0.2309  1.2309  3.9989 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.9089853  0.1917289  30.819  < 2e-16 ***
fincome     0.0012288  0.0002233   5.503 8.13e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.796 on 294 degrees of freedom
  (88 observations deleted due to missingness)
Multiple R-squared:  0.09339,   Adjusted R-squared:  0.09031 
F-statistic: 30.28 on 1 and 294 DF,  p-value: 8.127e-08

従属変数の観測値と予測値の相関係数の二乗

names(reg01)

 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "na.action"     "xlevels"       "call"          "terms"        
[13] "model"

R0 <- cor(d02$q1700, reg01$"fitted.values"); R0

[1] 0.3055957

R0^2

[1] 0.09338876

従属変数の観測値と独立変数の相関係数の二乗

R1 <- cor(d02$q1700, d02$fincome); R1

[1] 0.3055957

R1^2

[1] 0.09338876

lm( )関数をanova( )関数に引き渡した結果であるano01の情報から分散説明率を計算した結果は，lm( )関数をsummary( )関数に引き渡した結果であるsum01が表示する分散説明率(Multiple R-squared)に一致している。
自由度調整済み分散説明率(Adjusted R-squared)は無調整の分散説明率より小さい。これは母集団における分散説明率の推定値である。
（完備ケースの）従属変数の観測値とその予測値(fitted.values)の積率相関係数は約.3056，その二乗は分散説明率に一致する。
従属変数（観測値）と独立変数の積率相関係数R1はR0に一致する。予測値は独立変数の一次変換であり，積率相関係数は一次変換に関して（絶対値）不変である事から理解出来る。

lm( )$model からスマートに計算

　実はlm( )関数の結果であるreg01は，“model”という名前で，分析に用いられた独立変数と従属変数のデータを格納しているので，これを用いる方が相関係数の計算はスマートに出来る。

reg01$model

cor(reg01$model$fincome, reg01$model$q1700) # 独立変数と従属変数

[1] 0.3055957

cor(reg01$model$fincome, reg01$model$q1700)^2

[1] 0.09338876

cor(reg01$"fitted.values", reg01$model$q1700) # 従属変数の予測値と観測値

[1] 0.3055957

cor(reg01$"fitted.values", reg01$model$q1700)^2

[1] 0.09338876

2　回帰係数についての検討

2-1　相関係数と回帰係数

　標本相関係数と標本回帰係数の関係を確認する。

reg01$coefficients[2] # 標本回帰係数

    fincome 
0.001228757

# 標本相関係数R1とx,yの標準偏差から計算
R1*sd(d02$q1700)/sd(d02$fincome)

[1] 0.001228757

# 分散は偏差平方の平均である事を利用
R1*sqrt(mean((d02$q1700 - mean(d02$q1700))^2))/sqrt(mean((d02$fincome - mean(d02$fincome))^2))

[1] 0.001228757

選抜効果を描くグラフのスクリプト

　ここで，本文の様に選抜効果を確認するグラフを描くスクリプトを紹介する。
設定した相関係数にぴたりと一致するヴェクトル対の作成法などは単回帰分析や合成変数の分散・共分散について深い理解を必要とするので難しく，理解出来なくても支障無いが，グラフの描き方などは知っておくと色々と便利な工夫が含まれている。

r  <- .60 # 選抜前の相関係数を任意に設定
n  <- 400 # ケース数を任意に設定
y  <- rnorm(n, mean = 0, sd = 1)
x0 <- rnorm(n, mean = 0, sd = 1)
x  <- scale(x0)
e0 <- lm(y ~ x)$residuals
e <- scale(e0)

u <- sqrt(1-r^2)*e + r*x
cor(x, u) # 設定した相関係数にぴたりと一致するヴェクトル対が出来た

     [,1]
[1,]  0.6

X <- x*10 + 50 # 偏差値化
U <- u*10 + 50
cor(X, U) # 相関係数は当然変わっていない

     [,1]
[1,]  0.6

b <- lm(U ~ X)$coefficients[2] # 回帰係数をbに格納

# 選抜のない散布図
plot(X, U, 
     bty  = "n", 
     family = "serif", # Macの場合は要変更
     main = sprintf("相関係数%.2fの散布図(n=%d)", r, n),
     xlab = "", 
     ylab = "", 
     xlim = c(10, 90), 
     ylim = c(10, 90))
abline(lm(U ~ X), lty = 2) # 回帰直線を点線で書き込む
text(25, 70, sprintf("回帰係数%.2f", b))

# 選抜のある散布図
X1 <- X[X >= 50]; U1 <- U[X >= 50]
r1 <- cor(X1, U1); n1 <- length(X1)
b1 <- lm(U1 ~ X1)$coefficients[2]

plot(X1, U1, 
     bty  = "n", 
     family = "serif",
     main = sprintf("相関係数%.2fの散布図(n=%d)", r1, n1),
     xlab = "", 
     ylab = "", 
     xlim = c(10, 90), 
     ylim = c(10, 90))
abline(lm(U1 ~ X1), lty=2) # 回帰直線を点線で書き込む
text(25, 70, sprintf("回帰係数%.2f", b1))

シミュレイションなので，実行するたびに多少異なった結果になる（乱数の種seedを指定すれば再現可能）。
回帰係数が殆ど変化しない場合もあれば，大きく変化する場合もある。

2-2　標準化回帰係数(standardized regression coefficient)

　このサポートウェブでは既に活用して来た手法だが，変数をそれぞれ単独で欠損値除外するのではなく，使用する変数全てが有効であるケースだけを残さないと，多変量解析の結果が正しく再現出来ない恐れがある。
こうした分析は完備ケース分析(complete-case analysis)といい，欠損分析の観点からは決して問題が無いとは言えないが，初歩的分析では広く使用されている。
　標準化（偏）回帰係数を求める場合にも，完備ケースだけを選び出した後で独立変数と従属変数を標準化しないと正しい結果が得られない。
　まずは，完備ケースだけになっていないデータフレイム d01 について，独立変数と従属変数を単純にそれぞれ標準化して単回帰分析を行った結果を示す。

d01$z_fincome <- scale(d01$fincome)
d01$z_q1700   <- scale(d01$q1700)
summary(lm(z_q1700 ~ z_fincome, d01))


Call:
lm(formula = z_q1700 ~ z_fincome, data = d01)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2550 -0.6269  0.1221  0.6510  2.1148 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.05886    0.05522   1.066    0.287    
z_fincome    0.30388    0.05522   5.503 8.13e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.95 on 294 degrees of freedom
  (88 observations deleted due to missingness)
Multiple R-squared:  0.09339,   Adjusted R-squared:  0.09031 
F-statistic: 30.28 on 1 and 294 DF,  p-value: 8.127e-08

　標準化されているので独立変数も従属変数も平均は0の筈であり，回帰直線は点（xの平均，yの平均）を通る筈なので，この場合原点を通る筈である。
しかし定数項（切片）が誤差の範囲で0に一致しない。
また，xの標準偏差とyの標準偏差は共に1である筈なので，回帰係数は積率相関係数に一致し，その値は.3056であると既に求められているが，ここでの回帰係数は.3039でやはりずれている。
これらから，適切な変換がなされていない事が分かる。
　正しくは，完備ケースだけにしたデータフレイム d02 を用いる。

d02$z_fincome <- scale(d02$fincome)
d02$z_q1700   <- scale(d02$q1700)
summary(lm(z_q1700 ~ z_fincome, d02))


Call:
lm(formula = z_q1700 ~ z_fincome, data = d02)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2679 -0.6294  0.1226  0.6535  2.1232 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.640e-16  5.544e-02   0.000        1    
z_fincome    3.056e-01  5.553e-02   5.503 8.13e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9538 on 294 degrees of freedom
Multiple R-squared:  0.09339,   Adjusted R-squared:  0.09031 
F-statistic: 30.28 on 1 and 294 DF,  p-value: 8.127e-08

切片は（誤差の範囲で）0となり，回帰係数も上で求めた積率相関係数に一致している。

発展1　2群の母平均の差のt検定と単回帰分析

　性別と幸福度のスチューデントのt検定の結果を確認しておこう。

var.test(q1700 ~ sex, d01) # 一応，等分散性の検定


    F test to compare two variances

data:  q1700 by sex
F = 1.3234, num df = 160, denom df = 216, p-value = 0.05555
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.9933775 1.7743088
sample estimates:
ratio of variances 
          1.323409

t.test(q1700 ~ sex, d01, var.equal=T)


    Two Sample t-test

data:  q1700 by sex
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 -0.9337912 -0.1673909
sample estimates:
mean in group 1 mean in group 2 
       6.366460        6.917051

　次に，男性ダミー変数を作成して単回帰分析を行う。

\[\begin{align*} y_{i}=\beta_{0}+\beta_{1}d_{i}+e_{i} \end{align*}\]

d01$male <- c(1, 0)[d01$sex]
summary(lm(q1700 ~ male, d01))


Call:
lm(formula = q1700 ~ male, data = d01)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3665 -1.3665  0.0829  1.0829  3.6335 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.9171     0.1272  54.385  < 2e-16 ***
male         -0.5506     0.1949  -2.825  0.00498 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.874 on 376 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.02079,   Adjusted R-squared:  0.01818 
F-statistic: 7.982 on 1 and 376 DF,  p-value: 0.004977

　F統計量は上のt統計量の二乗に一致し，両者の有意確率も.004977で一致している。t統計量の自由度とF統計量の自由度2も勿論一致している。
　単回帰分析のCoefficientsのEstimateの(Intercept)はt検定の出力における女性の平均値6.917051と（四捨五入の範囲で）一致している。男性ダミー変数の係数である-.5506は，t検定の出力における男性平均マイナス女性平均の-.550591に（四捨五入の範囲で）一致する。よく見るとその回帰係数のt検定のt統計量と有意確率も，上の二群の母平均の差のt検定におけるt統計量と有意確率に一致している。単回帰であるから，回帰係数のt検定とモデルのF検定は実質的に同じであり，そしてそれらは二群の母平均差のt検定と実質的に同じなのである。数式表現については本文を参照して欲しい。

等分散正規性など，残差に要求される前提について

　以上で，Studentのt検定と，ダミー変数一個を用いた単回帰分析が同値である事を示した。Studentのt検定は，Welchの検定とは異なり，等分散正規性の条件を必要とするものであった（正規性の前提については，中心極限定理があるので，大標本の場合には統計的検定は頑健(robust)であるとの議論も多いが）。従って単回帰分析でも，残差について等分散正規性が要請されているのである。その他，各測定値が独立である（系列相関やクラスタとしての類似性が無い）と云う条件も満たしている必要がある。
　ここでは，回帰分析による検定・推定が妥当である為のこれらの条件のチェックの方法として，残差プロット(Residuals plot)による回帰診断(regression diagnostics)を紹介する。
　本章の1-2の，世帯年収で11段階幸福度を説明する単回帰の結果について，残差による診断を例示する。この単回帰の lm( ) の結果は既にreg01というオブジェクトとして保存してあるのでそれを使用する。

par(mfrow=c(2,2)) # 4つの診断図が出力されるので，作図領域を2行2列に分割しておく
plot(reg01) # 4つの診断図を出力する

このグラフの見方について，詳しくは以下のサイトなどを参照すると良い。
University of Virginia Library, Research Data Services + Sciences, ‘Understanding Diagnostic Plots for Linear Regression Analysis’(2019年9月4日閲覧)

　簡単に言えば，まず①左上の’Residuals vs Fitted’プロットで，各予測値(fitted value)ごとの残差が，適切に線型パタンを示しているかどうかをチェックする。大体水平に近い赤い線が引かれれば問題無い。
　次に右上の②’Normal Q-Q’プロット（正規確率プロット）で，残差の正規性が満たされているかどうかをチェックする。ほぼ45°線上に並んでいれば，理論的な正規分布からのずれが大きくないのでOKである。
　左下の③’Scale-Location’では，等分散性のチェックが行われている。どの予測値に対しても，同じ様に「標準化残差の平方根」が分布していれば問題ない。
　最後に右下の④’Residuals vs Leverage’プロットでは，全体の結果に対する影響力の大きいケース(high-leverage)で大きな標準化残差を持つものがあると，それだけで全体の分析結果を大きく左右する（つまり問題を孕んだ外れ値の可能性が高い）事を示す。Leverageの高い右側の方で，Standardized residualが0から遠い点に，番号が付される。これらは，全体に対して影響力の高い点であり，それを含めた分析と除外した分析で結果が少なからず変化する事を意味している。
　以上の様にして，残差と予測値に関連があるか否か①，残差が正規性を充たすか②，残差が等分散性を充たすか③，問題を孕んだ外れ値が無いか④をチェックする事が出来る。

発展2　偏相関係数(partial correlation coefficient)

　ゼロ次の積率相関係数行列と，偏相関係数行列を出力する。R Commander(Rcmdr)をinstallしていない場合はインストールから行う。library(Rcmdr)でR Commanderを有効化した時に，更に追加のパッケイジを求められた場合には，指示の通りにインストールすればよい。

# install.packages("Rcmdr", repos="https://cran.ism.ac.jp")
# library(Rcmdr)
cor(cbind("edu"=d01$edu, "income"=d01$income, "book"=d01$q0200), use="complete")

             edu     income       book
edu    1.0000000 0.27521938 0.38115890
income 0.2752194 1.00000000 0.09382867
book   0.3811589 0.09382867 1.00000000

partial.cor(cbind("edu"=d01$edu, "income"=d01$income, "book"=d01$q0200), use="complete")


 Partial correlations:
           edu   income     book
edu    0.00000  0.26016  0.37125
income 0.26016  0.00000 -0.01246
book   0.37125 -0.01246  0.00000

 Number of observations: 338

本文脚注でも述べたが，bookは順序尺度なので本来は間隔尺度に変換する方が良いだろうが，ここでは単純化の為にそのまま使用している。
　教育年数eduを統制した場合の年収incomeと蔵書bookの偏相関係数とは，incomeをeduで説明する回帰分析の残差と，bookをeduで説明する回帰分析の残差との相関係数である。その事を以下で示そう。

# 三つの変数が全て有効であるケースに限定する準備（完備ケース分析）
boolean <- complete.cases(d01$edu, d01$income, d01$q0200)
x <- d01$income[boolean]; y <- d01$q0200[boolean]; z <- d01$edu[boolean]

# 教育年数で年収を，教育年数で蔵書数を説明する回帰分析
summary(lmx <- lm(x ~ z)); summary(lmy <- lm(y ~ z))


Call:
lm(formula = x ~ z)

Residuals:
    Min      1Q  Median      3Q     Max 
-449.02 -237.13  -62.13  150.98 2050.98 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -298.531    124.401  -2.400   0.0169 *  
z             46.722      8.904   5.248 2.74e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 326.9 on 336 degrees of freedom
Multiple R-squared:  0.07575,   Adjusted R-squared:  0.07299 
F-statistic: 27.54 on 1 and 336 DF,  p-value: 2.737e-07


Call:
lm(formula = y ~ z)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5265 -0.9817  0.0183  0.5632  3.5632 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.83212    0.50364  -1.652   0.0994 .  
z            0.27241    0.03605   7.557 3.95e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.323 on 336 degrees of freedom
Multiple R-squared:  0.1453,    Adjusted R-squared:  0.1427 
F-statistic: 57.11 on 1 and 336 DF,  p-value: 3.947e-13

# それぞれの単回帰分析の残差同士の積率相関係数
cor(lmx$residuals, lmy$residuals)

[1] -0.01245901

第8章の【練習問題】の解答

1) の答え

　ここではデータフレイム名は d01 とするが，最初に data01 として読みこんでいる場合はその様に指定する。

full <- complete.cases(d01$edu, d01$income)
edu <- d01$edu[full]; income <- d01$income[full]
cor.test(edu, income)


    Pearson's product-moment correlation

data:  edu and income
t = 5.2475, df = 336, p-value = 2.737e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1736406 0.3710037
sample estimates:
      cor 
0.2752194

　標本相関係数は.275程度（95%信頼区間は.174～.371）であり，ゼロ仮説「母相関係数は0である」は有意水準0.1%でも棄却される。

2) の答え

summary(lm(income~edu))


Call:
lm(formula = income ~ edu)

Residuals:
    Min      1Q  Median      3Q     Max 
-449.02 -237.13  -62.13  150.98 2050.98 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -298.531    124.401  -2.400   0.0169 *  
edu           46.722      8.904   5.248 2.74e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 326.9 on 336 degrees of freedom
Multiple R-squared:  0.07575,   Adjusted R-squared:  0.07299 
F-statistic: 27.54 on 1 and 336 DF,  p-value: 2.737e-07

　1)の相関係数のt検定と，この2)の回帰係数のt検定は一致している。t統計量5.2475の二乗は27.53626となり，F統計量27.54に一致する。自由度や有意確率も一致している。また，相関係数.2752194の二乗は.07574572で，分散説明率.07575に一致する。

3) の答え

names(cor.test(edu, income))

[1] "statistic"   "parameter"   "p.value"     "estimate"    "null.value" 
[6] "alternative" "method"      "data.name"   "conf.int"

cor.test(edu, income)$estimate

      cor 
0.2752194

sd(edu); sd(income)

[1] 2.000035

[1] 339.5313

cor.test(edu, income)$estimate*sd(income)/sd(edu)

     cor 
46.72197

lm(income ~ edu)$coefficients[2] # 回帰係数

     edu 
46.72197

cor.test(edu, income)$estimate^2 # 相関係数の二乗

       cor 
0.07574571

summary(lm(income ~ edu))$"r.squared" #分散説明

[1] 0.07574571

4) の答え

必要最小限なら下記のスクリプトで実行できる。

t.test(income ~ sex, d01, var.equal=T)


    Two Sample t-test

data:  income by sex
t = 11.307, df = 340, p-value < 2.2e-16
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 296.3890 421.2202
sample estimates:
mean in group 1 mean in group 2 
       554.9296        196.1250

d01$female <- c(0, 1)[d01$sex]
summary(lm(income ~ female, d01))


Call:
lm(formula = income ~ female, data = d01)

Residuals:
    Min      1Q  Median      3Q     Max 
-554.93 -171.13  -54.93  103.87 1945.07 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   554.93      24.27   22.87   <2e-16 ***
female       -358.80      31.73  -11.31   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 289.2 on 340 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared:  0.2733,    Adjusted R-squared:  0.2711 
F-statistic: 127.9 on 1 and 340 DF,  p-value: < 2.2e-16

だが，それぞれの対応関係を確認する為にはやはり一旦オブジェクトに格納した方が便利である。

t01 <- t.test(income ~ sex, d01, var.equal=T)
r01 <- summary(lm(income ~ female, d01))
names(t01); names(r01)

 [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
 [6] "null.value"  "stderr"      "alternative" "method"      "data.name"

 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"

t01$statistic^2 # t値の二乗

       t 
127.8565

r01$fstatistic # F値

   value    numdf    dendf 
127.8565   1.0000 340.0000

diff(t01$estimate) # 標本平均の差

mean in group 2 
      -358.8046

r01$coefficients[2, 1] # 回帰係数

[1] -358.8046

t01$"conf.int" # 平均値差の95%信頼区間

[1] 296.3890 421.2202
attr(,"conf.level")
[1] 0.95

confint(lm(income ~ female, d01))[2, ] # 回帰係数の95%信頼区間

    2.5 %    97.5 % 
-421.2202 -296.3890

0 全体の章構成

1 単回帰分析の基礎

1-1 最小二乗和法(Least Squares Method)

1-2 回帰係数(regression coefficient)のt検定と区間推定(interval estimation)

1-3 回帰モデルのF検定

1-4 決定係数(coefficient of determinance)，分散説明率(proportion of variance explained)，誤差減少率(PRE)

2 回帰係数についての検討

2-1 相関係数と回帰係数

2-2 標準化回帰係数(standardized regression coefficient)

発展1 2群の母平均の差のt検定と単回帰分析