.

.

データ解析
第三回「回帰分析」
鈴木 大慈
理学部情報科学科
西八号館 W707 号室
s-taiji@is.titech.ac.jp

1 / 23

今日の講義内容
回帰分析 (lm)lm 関数返り値の解釈
回帰係数の有意性検定
AIC によるモデル選択

2 / 23線形回帰モデル
ガウス・マルコフモデル:

y = β1 x1 + β2 x2 + · · · + βd xd + βd+1 + ϵ(ϵ ∼ N(0, σ 2 ))y ∈ R: 従属変数x ∈ Rd : 説明変数
βi (i = 1, . . . , d):（偏）回帰係数, βd+1 : 切片 (切片項は省略することもある)
単回帰 (d = 1)y = β1 x1 + β2 + ϵ
重回帰 (d ≥ 1)

y = β1 x1 + β2 x2 + · · · + βd xd + βd+1 + ϵ

3 / 23

204060area
801004 / 23

2e+07

4e+07

price
6e+07

8e+07

最小二乗法n サンプル (観測値)：(yi , xi ) ∈ R × Rd (i = 1, . . . , n)
  y1— x⊤
ϵ1
1 — 1
 

..
..  ∈ Rn×(d+1) , ϵ =  ..  ∈ Rn ,Y =  ...  ∈ Rn , X = ..
.yn
ϵn— x⊤—
1nとおき，β ∗ を真の回帰係数 (これを推定したい) とすると,

Y = X β ∗ + ϵ.
最小二乗推定量 (最尤推定量)：β̂ = (X ⊤ X )−1 X ⊤ Y .

不偏推定量Cramer-Rao の下限を達成 (一様最小分散不偏推定量)
5 / 23最尤推定量の従う分布.
X 固定のもと，β̂ の分布は次式で与えられる：(()−1 )√1 ⊤∗2n(β̂ − β ) ∼ N 0, σX X.
n.
(導出)β̂ = (X ⊤ X )−1 X ⊤ Y = (X ⊤ X )−1 X ⊤ (X β ∗ + ϵ)√√⇒ n(β̂ − β ∗ ) = n(X ⊤ X )−1 X ⊤ ϵ.

(1)ϵi ∼ N(0, σ 2 ) より，右辺は平均 0 の多変量正規分布. その分散共分散行列は
E[n(X ⊤ X )−1 X ⊤ ϵϵ⊤ X (X ⊤ X )−1 ] = n(X ⊤ X )−1 X ⊤ E[ϵϵ⊤ ]X (X ⊤ X )−1 =()−1σ 2 n1 X ⊤ X.
□
これより回帰係数の信頼区間が求まるが，σ 2 を知らない場合はそれも推定する必要がある (次ページ)．
6 / 23最尤推定量の信頼区間
ϵ̂ を ϵ の推定量として
ϵ̂ = Y − X β̂ = X (β ∗ − β̂) + ϵ = (I − X (X ⊤ X )−1 X ⊤ )ϵとおく．すると，S = ( n1 X ⊤ X )−1 としたとき，.

.

√√n(β̂j − βj∗ )/ Sjj√
∼ t(n − (d + 1))∥ϵ̂∥2 /(n − (d + 1))(自由度 n − (d + 1) の t 分布).

これより，βj の信頼区間や検定が可能になる．
例: βj∗ = 0 の検定．もし√√n|β̂j |/ Sjj√
≥ tα (n − (d + 1))∥ϵ̂∥2 /(n − (d + 1))ならば，βj∗ = 0 なる帰無仮説を棄却する．ただし，tα は t 分布の上側 α 分位点．

(βj∗ = 0 にしては β̂j が大きすぎ)
7 / 23

導出. まず，1

ϵ̂ = Y − X β̂ = X (β ∗ − β̂) + ϵ = (I − X (X ⊤ X )−1 X ⊤ )ϵに気をつける．これは平均 0 の多変量正規分布で，かつβ̂ − β ∗ = (X ⊤ X )−1 X ⊤ ϵ (Eq. (1)) と独立である．独立性は，
E[ϵ̂(β̂ − β ∗ )⊤ ] = O であることと，ϵ̂ と β̂ − β ∗ の同時分布が正規分布であることより，それらは独立であることがわかる (チェックせよ)．

. 次に，(I − X (X ⊤ X )−1 X ⊤ ) は n − (d + 1) 次元部分空間への射影なので，∥ϵ̂∥2 /σ 2 は自由度 n − (d + 1) の χ2 分布に従う．
√.3 また，最尤推定量の分布より √n(β̂j − β ∗ )/ Sjj ∼ N(0, σ 2 ) である．
j√√n(β̂j −βj∗ )/ Sjj
4. よって， √は互いに独立な正規分布と χ2 分布の比なので t 分22∥ϵ̂∥ /(n−(d+1))

布となる．σ 2 は分母と分子で打ち消す．

□8 / 23

回帰の有意性: F 検定
そもそも回帰に意味がないかもしれない．回帰式に説明力があることを検定．

H0 : β ∗1 = β ∗2 = · · · = β ∗d = 0vsH1 : それ以外(H0 でも切片の非ゼロは許す).
帰無仮説のもとで，F =.ただし，Ȳ =∑ni=1 yin∥X β̂ − 1Ȳ ∥2 /d∼ F (d, n − d − 1).
∥ϵ̂∥2 /(n − (d + 1))．

導出は確率統計第二の講義を参照．また標準的な数理統計の教科書には大体書いてある．
2)：
決定係数 (R 2 ) と自由度調整済み決定係数 (RAR2 = 1 −∥Y − X β̂∥2,∥Y − 1Ȳ ∥22=1−RA∥Y − X β̂∥2 /(n − d − 1)
∥Y − 1Ȳ ∥2 /(n − 1)

決定係数が大きければ帰無仮説は棄却され，回帰式には説明力があることになる．
9 / 23

AIC による変数選択
重回帰分析では説明変数を追加するごとに残差は小さくなってゆく．
それでいいのか？
答え：必ずしも良くはない．
観測済みデータによく当てはまっても，未観測のデータへの当てはまりが悪くなることがある．
過適合 (overﬁt) という．
過適合を避けるにはサンプル数に比して複雑なモデルは使うべきではない．
.

赤池情報量規準, AIC (Akaike Information Criterion).
予測精度が一番良いモデルを選択するための規準
d 次元パラメトリックモデル Md の AIC は次式で定義される:
AIC = −2 log(p({Xi }ni=1 |θ̂Md )) + 2m.
.ただし，θ̂Md は Md 上の最尤推定量．
言い換えればAIC = −2 × 最尤推定量の対数尤度 + 2 × モデルの次元である．AIC が最小になるように説明変数の数を選ぶ．

10 / 23

AIC の解釈次元が高い (複雑な) モデルであるほど負の対数尤度が小さくなるが，それにモデルの次元 (複雑さ) を罰則として加えることで過適合を防ぐ．
モデルの次元とデータへの当てはまりのトレードオフ．

AIC は真の分布との KL-divergence を漸近展開した１次項である．
予測誤差 (真との KL-divergence) を最小化．変数を選択するための客観的規準 (予測誤差) を一つ持ってきたところがミソ．
(導出は省略)11 / 23

中古マンション価格データによる実演

国土交通省が公開している不動産取引価格情報から世田谷区の中古マンション取
引価格データ (平成２５年度第３四半期分) を取得．ここから一部を抜粋した
データで回帰分析をやってみる．
http://www.land.mlit.go.jp/webland/download.html

従属変数：価格
説明変数：最寄り駅からの距離（徒歩），延床面積，建物の構造，建ぺい率，容積率，建築年，最寄り駅に急行が止まるか

12 / 23

ダミー変数
説明変数「建物の構造」は「RC」と「SRC」という文字列 (カテゴリカル変数)．
これを {0, 1} に置き換える．ダミー変数という．
RC → 1, SRC → 0 とした．

13 / 23線形回帰分析のための関数
書式:
lm(formula, data, subset, weights, na.action,method = ”qr”, model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
singular.ok = TRUE, contrasts = NULL, oﬀset, ...)
引数，大事な部分だけ抜粋
formula
当てはめられるモデルのシンボリックな記述式
data
モデル中の変数を含むオプションのデータフレーム
subsetオプションのベクトルで，当てはめ過程で使われる観測値の部分集合を指定
weights
当てはめ過程で使われるオプションの重みベクトルna.action データが NA を含むときそれを処理する関数．

14 / 23

formula の書式
y˜xy ˜ x1 + x2
y ˜ x1 * x2
y˜x-1
y ˜ 1 + x + I(xˆ2)
y˜x—zy ˜ . , data = データ名y = β1 x + b + ϵ
y = β1 x1 + β2 x2 + b + ϵ
交互作用項を含んだモデル
y = β1 x1 + β2 x2 + β3 x1x2 + b + ϵ
切片 (定数) 項を除外する．y = β1 x + ϵ多項式回帰: y = b + β1 x + β2 x 2 + ϵz で条件付けしたときの y の x への単回帰．
あるデータに目的変数 y と説明変数 x1 , . . . , xd
∑dが入っている場合，y = b + j=1 βj xj + ϵ とする．

例：

x <- seq(-10,10); y<- 3*x + rnorm(21);lm(y ~ x) # y を x に回帰
x <- seq(-10,10); z <- seq(-10,10)^2; y<- 3*x + 4*z + rnorm(21);datain <- data.frame(x=x,z=z,y=y)lm(y ~ x, data=datain) # y を x に回帰
lm(y ~ ., data=datain) # y を y 以外の変数に回帰
15 / 23データの取得とプロット

x <- read.csv("setagaya_manshion.csv")
sman = data.frame(price=x[[1]],walk=x[[3]],area=x[[5]],str=ifelse(x[[6]]=="ＲＣ",1,0),kenpei=x[[7]],youseki=x[[8]],tikunen=x[[9]],kyuko=x[[10]])#構造はダミー変数 0-1 に置き換える．
plot(sman)#散布図のプロット010

200.00.6

100

4000.00.6

walk80010202e+07

price

str80

0.00.620area

youseki

2000

100

400

4060

kenpei0.00.61970

tikunenkyuko

16 / 23回帰分析関数 (lm)

2e+07

4e+07

price

6e+07

8e+07

sman.lm <- lm(price ~ area,data=sman) #回帰分析はこの一行で OK！
plot(sman$area,sman$price, xlab="area",ylab="price") #結果をプロッ
トabline(sman.lm , lwd=1 , col="blue")

20

4060

80

100

area

17 / 23

分析の要約 (summary)> summary(sman.lm)
Call:lm(formula = price ~ area, data = sman)
Residuals:
Min
1Q-16082568 -5634345

Median-789511

3Q4806399

Max16822060

Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept) -1029167
1519818 -0.677
0.5area673025
26408 25.485<2e-16 ***--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05

.

0.1

1

Residual standard error: 7594000 on 128 degrees of freedom
Multiple R-squared: 0.8354,
Adjusted R-squared: 0.8341
F-statistic: 649.5 on 1 and 128 DF, p-value: < 2.2e-16
18 / 23

要約の意味 (1)

Call:lm(formula = price ~ area, data = sman)呼び出し式Residuals:
Min1QMedian-16082568 -5634345-789511⊤
残差 yi − yˆi = yi − β̂ x の要約統計量

3Q4806399

Max16822060

19 / 23

要約の意味 (2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept) -1029167
1519818 -0.677
0.5
area673025
26408 25.485<2e-16 ***--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 .
1
回帰係数の推定量，標準偏差，t-統計量，p-値(Intercept) は切片. t-統計量の意味は前のスライドを参照．

0.1



ここの p-値がある有意水準より小さければ，その回帰係数が 0 であるという帰無仮説は棄却される．つまり，その変数は y を説明するための情報を持っているということ．
この例だと，切片項は有意でなく，area は有意である (たとえば有意水準
α = 0.05 のとき)．
20 / 23

要約の意味 (3)Residual standard error: 7594000 on 128 degrees of freedom
残差の標準偏差 (σ の推定量)．この場合，サンプル数 130 で変数は 1 変数+切片なので自由度は 130-2=128.
F-statistic: 649.5 on 1 and 128 DF, p-value: < 2.2e-16
β = 0 の検定における F 統計量．自由度は (d, n − d − 1) (前のスライドを参照)．
この場合は回帰は有意であると言える (α = 0.05)．
Multiple R-squared: 0.8354,Adjusted R-squared: 0.8341

決定係数と自由度調整済み決定係数 (Adjusted R-squared)．Adjusted のほうが重要．大きいほど相関が強い．

21 / 23回帰分析の診断par(mfrow=c(2,2)) #一つの Figure に 2x2 の図を出力
plot(sman.lm)

3e+075e+072
139−2−101

Theoretical Quantiles

Scale−Location

Residuals vs Leverage

5e+07

Fitted values

7e+0701243
45−1−20.5

1.0Standardized residuals

58(下左) Scale-Location プロット:規準化した残差の絶対値の平方根．この値は予測値の値に対して大体一定であるべき．そうでない場合，場所によってノイズの分散が異なることになる．

2Fitted values

3e+07(上左) 残差プロット:
予測値 ŷi vs 残差 yi − ŷi(上右) 残差の Q-Q プロット:
正規分布との Q-Q プロット．対角線に近ければより正規分布に近い．残差の正規性はガウスマルコフモデルの仮定．

0−27e+07583996

1e+0796 58−1Standardized residuals

0e+00

391e+070.0

Standardized residuals

Normal Q−Q5896

1.5−2e+07

Residuals

2e+07

Residuals vs Fitted

Cook’s distance(下右) Cook 距離のプロット:
個々のデータが推定に与える影響を表した距離 (そのデータがない場合とある場合の予測値の変化量). 大きいと外れ値の可能性がある．0.5 を超えると「大きい」とされる．

0.00 0.01 0.02 0.03 0.04 0.05
Leverage22 / 23重回帰分析
築年数も含めて分析してみる．
> sman.lm3 <- lm(price ~ area + tikunen, data=sman)> (AIC(sman.lm3))[1] 4443.881
AIC(·) で AIC が計算できる．
すべての説明変数を用いて回帰分析






sman.lmall <- lm(price ~., data=sman)
sman.lmAIC <- step(sman.lmall)summary(sman.lmAIC)

step(·) で AIC 最小のモデル (説明変数の組) を探索．walk + area + tikunen の三変数モデルが採用された．

23 / 23