経営統計入門　期末対策講座　回帰編

はじめに

今回の例題はこちらのデータを使うので、コピーして実際に動かしながら見てみてください

time <- c(1, 3, 4, 4, 7, 2, 5, 2, 6, 9)
score <- c(24, 58, 76, 42, 70, 30, 80, 27, 90, 95)

Code

plot(time, score)

Ch12

回帰

データの傾向を読み取る

データに対してそれらしい直線を引く

→未知の値に対して予測ができる

例

勉強時間が8時間

→94点

赤い回帰直線は、\(\hat{y}=\beta_1 x+\beta_0\)の式で表す

\(\beta_1\): 傾き（xが1増えた時yがどれだけ増えるか）

\(\beta_0\): 切片（xが0の時のyの値）

回帰直線の求め方

最小二乗法

\(\beta_1=\frac{S_{xy}}{S_{xx}}\), \(\beta_0=\bar{y}-b\bar{x}\)

\(S_{xy}=\Sigma(x_i-\bar{x})(y_i-\bar{y})=\Sigma x_iy_i-\frac{(\Sigma x_i)(\Sigma y_i)}{n}\)

\(S_{xx}=\Sigma(x_i-\bar{x})^2=\Sigma x^2_i-\frac{(\Sigma x_i)^2}{n}\)

Code

# データの個数を記録
n <- length(time) # length(score)でも可

# S_xy
S_xy <- sum(time*score)-sum(time)*sum(score)/n
# S_xx
S_xx <- sum(time^2)-sum(time)^2/n

# beta_1
beta_1 <- S_xy/S_xx
# beta_0
beta_0 <- mean(score)-beta_1*mean(time)

beta_0; beta_1

[1] 19.082

[1] 9.329768

Code

# 回帰直線を図示
plot(time, score)
abline(a = beta_1, b = beta_0)

R

# lm(予測したい変数~予測に用いる変数)
lm_fit <- lm(score~time)
lm_fit

Call:
lm(formula = score ~ time)

Coefficients:
(Intercept)         time  
      19.08         9.33  

Coefficients: (Intercept): y切片

time: 傾き（名前はデータの項目名によって変わる）

回帰式: \(y=9.33x+19.08\)

Code

# 回帰直線を図示
plot(time, score)
abline(lm_fit)

練習問題

• 問題
• 解説 (1)
• 解説 (2)

x	y
2	3
2	4
4	6
5	5
3	4

1. 左のようなデータを得た時、データを正しく表現していてかつ正しい回帰直線が引けているグラフはどれか

2. また、回帰直線の傾きと切片は何か

問題のデータを、xとyに分けてRに入れる

x <- c(2, 2, 4, 5, 3)
y <- c(3, 4, 6, 5, 4)

回帰を行う

• 最小二乗法の場合

Code

# データの個数を記録
n <- length(x) # length(score)でも可

# S_xy
S_xy <- sum(x*y)-sum(x)*sum(y)/n
# S_xx
S_xx <- sum(x^2)-sum(x)^2/n

# beta
beta <- S_xy/S_xx
# alpha
alpha <- mean(y)-beta*mean(x)

• Rの場合

Code

# lm(予測したい変数~予測に用いる変数)
lm_fit <- lm(y~x)

ここで得た情報をプロットすると

Code

plot(x, y)

# 最小二乗法の場合
abline(a = alpha, b = beta)
# Rの場合
abline(lm_fit)

よって答えは1

(1)より

• 最小二乗法の場合

Code

alpha; beta

[1] 2.235294

[1] 0.6764706

• Rの場合

Code

lm_fit

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
     2.2353       0.6765  

よって

傾き: 0.6765

切片: 2.2353

回帰の精度

回帰直線によって、データのばらつきをどの程度表現できているかみる

\(S_{yy}=\Sigma y_i^2-\frac{(\Sigma{y_i})^2}{n}\)

# S_yy
S_yy <- sum(score^2)-sum(score)^2/n

\(S_{yy}\) / \(TSS\): yの全体のばらつき（実測値から平均までのばらつき）

\(SSR\): yの平均から回帰直線までのばらつき

\(SSE\): 回帰直線から実測値までのばらつき

\(TSS=SSR+SSE\)

分散分析

MSR: yの平均から回帰直線までの自由度1に対する平均のばらつき

MSE: 回帰直線までから実測値までの自由度1に対する平均のばらつき

Source	df	SS	MS	F
Regression	1	\(SSR=\frac{(S_{xy})^2}{S_{xx}}\)	\(MSR=\frac{SSR}{1}\)	\(F=\frac{MSR}{MSE}\)
Error	n-2	\(SSE=S_{yy}-\frac{(S_{xy})^2}{S_{xx}}\)	\(MSE=\frac{SSE}{n-2}\)
Total	n-1	\(SST=S_{yy}\)

この表に従って最終的にF値を求め、この回帰直線が未知の値に対しても予測できるかを調べる→

\(H_0: 有効性=0\)（この回帰直線は役に立たない）

\(H_a: 有効性\neq0\)（この回帰直線は役に立たないことはない）

練習問題

• 問題
• ヒント
• 解答

勉強時間とテストの得点を表すデータを使って、分散分析表を完成させなさい

Source	df	SS	MS	F
Regression
Error
Total

\(S_{xx}=\) 6.8

\(S_{xy}=\) 4.6

Code

# データの個数を記録
n <- length(score) # length(score)でも可

# S_xy
S_xy <- sum(time*score)-sum(time)*sum(score)/n
# S_xx
S_xx <- sum(time^2)-sum(time)^2/n
# S_yy
S_yy <- sum(score^2)-sum(score)^2/n

# SSR
SSR <- S_xy^2/S_xx
# SSE
SSE <- S_yy-S_xy^2/S_xx
# SST
SST <- S_yy

# MSR
MSR <- SSR/1
# MSE
MSE <- SSE/(n-2)

# F
F_value <- MSR/MSE

Source	df	SS	MS	F
Regression	\(1\)	\(SSR=\frac{(S_{xy})^2}{S_{xx}}=\) 4883.2007	\(MSR=\frac{SSR}{1}=\) 4883.2007	\(F=\frac{MSR}{MSE}=\) 24.3491
Error	\(n-2=\) 8	\(SSE=S_{yy}-\frac{(S_{xy})^2}{S_{xx}}=\) 1604.3993	\(MSE=\frac{SSE}{n-2}=\) 200.5499
Total	\(n-1=\) 9	\(SST=S_{yy}=\) 6487.6

決定係数

回帰直線が、データに対してどれほど当てはまりがいいか

\(決定係数r^2=\frac{SSR}{SST}\)

この値は相関係数の二乗と同じ値になる

# 相関係数の二乗として計算
cor(time, score)^2

[1] 0.7526976

# 回帰から得られる情報として算出
lm_fit <- lm(score~time)
summary(lm_fit)

Call:
lm(formula = score ~ time)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.401 -10.069  -6.077  13.434  19.599 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   19.082      9.282   2.056  0.07384 . 
time           9.330      1.891   4.934  0.00114 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.16 on 8 degrees of freedom
Multiple R-squared:  0.7527,    Adjusted R-squared:  0.7218 
F-statistic: 24.35 on 1 and 8 DF,  p-value: 0.001143

Multiple R-squaredが決定係数を示す

回帰中のt検定

２つのデータの間に関係性があるかを検定する

もし傾き\(\beta_1\)が0だったらどうなるのか

傾きが0だとxがどんな値になろうともyの値が一定

→xとyの間に関係がない

\(H_0:\beta_1=0\)を棄却できれば、\(\beta_1\neq0\)ではない

→xとyの間に関係があるといえる

回帰中で使うt検定のdfは\(n-2\)

勉強時間とテストの点数の回帰

Call:
lm(formula = score ~ time)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.401 -10.069  -6.077  13.434  19.599 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   19.082      9.282   2.056  0.07384 . 
time           9.330      1.891   4.934  0.00114 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.16 on 8 degrees of freedom
Multiple R-squared:  0.7527,    Adjusted R-squared:  0.7218 
F-statistic: 24.35 on 1 and 8 DF,  p-value: 0.001143

t valueが\(\beta_1=0\)とした時のt値を表す

棄却域に入らない

→傾きは0ではない/２つのデータの間には関係性がある

練習問題

• 問題
• 解説 (1)
• 解説(2)

1. 勉強時間とテストの得点を示すデータを使って、この２つの要素間に関係があるか99%信頼区間を用いてt検定しなさい。

2. 作成した回帰式とデータの当てはまりを示す決定係数を算出しなさい。

帰無仮説\(H_0:\beta_1=0\)とした時のt値を算出する

lm_fit <- lm(score~time)
summary(lm_fit)

Call:
lm(formula = score ~ time)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.401 -10.069  -6.077  13.434  19.599 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   19.082      9.282   2.056  0.07384 . 
time           9.330      1.891   4.934  0.00114 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.16 on 8 degrees of freedom
Multiple R-squared:  0.7527,    Adjusted R-squared:  0.7218 
F-statistic: 24.35 on 1 and 8 DF,  p-value: 0.001143

このデータのサンプル数は\(n-2=10-2=8\)なので、自由度\(df=9\)のt分布で99%信頼区間となる棄却域を求める

# 信頼係数
a <- 0.99

# 棄却域のはじまるt値を求める
qt((1-a)/2, 8)

[1] -3.355387

timeのt valueは4.934で、棄却域は3.36からなので、このt値は棄却域に入る

よってtimeは99%信頼区間においてテストの得点に関係があると示された

• \(\frac{SSR}{SST}\)を使う場合

# データの個数を記録
n <- length(time) # length(score)でも可

# S_xy
S_xy <- sum(time*score)-sum(time)*sum(score)/n
# S_xx
S_xx <- sum(time^2)-sum(time)^2/n

# SSR
SSR <- S_xy^2/S_xx
# SST
SST <- sum(score^2)-sum(score)^2/n

# r_2
r_2 <- SSR/SST

r_2

[1] 0.7526976

• 相関係数を使う場合

# 相関係数を計算
cor(time, score)^2

[1] 0.7526976

# 回帰を行う
lm_fit <- lm(score~time)
# 詳細な情報を取得
summary(lm_fit)

Call:
lm(formula = score ~ time)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.401 -10.069  -6.077  13.434  19.599 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   19.082      9.282   2.056  0.07384 . 
time           9.330      1.891   4.934  0.00114 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.16 on 8 degrees of freedom
Multiple R-squared:  0.7527,    Adjusted R-squared:  0.7218 
F-statistic: 24.35 on 1 and 8 DF,  p-value: 0.001143

Multiple R-squared: 0.7525

Ch13

重回帰分析

複数の項目を使って、値を予測する

今までの勉強時間とテストの得点のデータに加え、授業の出席割合のデータを加えます

attend <- c(0.7, 0.78, 0.9, 0.72, 0.83, 0.64, 0.95, 0.3, 1, 1)

3軸のグラフで回帰を当てはめると面になる

回帰を当てはめ、その回帰がデータに対してどの程度当てはまっているのか、その予測がどれほど正確なのかを検証していく

多重回帰の求め方

• R

lm_fit <- lm(score~time+attend)
lm_fit

Call:
lm(formula = score ~ time + attend)

Coefficients:
(Intercept)         time       attend  
    -14.476        5.542       63.743  

Coefficients: (Intercept)(切片): -14.476

time(傾き 1): 5.542

attendance(傾き 2): 63.743

\(y=\beta_1 x_1 + \beta_2 x_2 + \beta_0\)

多重回帰における仮説設定

回帰ではまず予測に使用している項目が予測したい項目に対して影響があるかを調べる

使用しているすべての項目を一括りにして、予測する項目に影響があるかをみる

→\(H_0:\beta_1=\beta_2=0\)

この仮説が棄却される（この項目の中の一つ、またはそれ以上の項目は予測する項目に影響を与えている）

→個々の変数に焦点を絞り検定をする

分散分析表

Source	df	SS	MS	F
Regression	予測に使う変数の数: \(x\)	\(SSR=SST-SSE\)	\(MSR=\frac{SSR}{df}\)	\(F=\frac{MSR}{MSE}\)
Error	Total df-x	\(SSE=SST-SSR\)	\(MSE=\frac{SSE}{Total df-x}\)
Total	n-1	\(SST=SSR+SSE\)

Regressionのdfは予測につかっている変数の数が入る

→Ch12では、テストの点数の予測に勉強時間のみをつかっていたため1が入る

→今回は予測に勉強時間と出席率の２つを使っているため、2が入る