データ解析第三回「回帰分析」鈴木　大慈理学部情報科学科西八号館W707号室s-taiji@is.titech.ac.jp
。回帰分析(Im) 。 Im関数返り値の解釈●回帰係数の有意性検定● AICによるモデル選択
ガウス・マルコフモデル: y = β1x1 + β2x2 + · · + βdxd + βd+1 + c    (ε ~ N(0,σ²)) y ∈ R:従属変数x ∈ R :説明変数β i (i = 1, . . ,d ): (偏)回帰係数, β +1:切片(切片項は省略することもある)単回帰(d = 1) y = β1×1 + β2 + c重回帰(d > 1) y = β1x1 + β2x2 + · · + βaxa + βd+1 + c
nサンプル(観測値) : (y ; , × ;) ∈ R × R d (i = 1, . . , n)
最小二乗推定量(最尤推定量) :
●不偏推定量● Cramer-Raoの下限を達成(一様最小分散不偏推定量)
X固定のもと，βの分布は次式で与えられる:
(導出)
⇒ √ n(β − β*) = √ n(X x)−¹X ε.
ε; ～ N(0.σ²)より，右辺は平均0の多変量正規分布.その分散共分散行列はE[n(XTX)-1XTecTX(XTX)-1] = n(XTX)-1XTE[ecT]X(XTX)-1 = σ² (≒XTX)−¹.これより回帰係数の信頼区間が求まるが，σ²を知らない場合はそれも推定する必要がある(次ベージ).
6 / 1

êをεの推定量として
とおく.すると, Ｓ＝(≒XTX)−¹としたとき, √ n ( β − β ∗ )/ √ S jj = ~ t(n − (d + 1))   (自由度n − (d + 1)のt分布). √||ê||²/( n − (d + 1))これより，βの信頼区間や検定が可能になる.例: β* = 0の検定.もし
ならば，βi = 0なる帰無仮説を棄却する．ただし，t。はt分布の上側α分位点. (β* = 0にしてはβiが大きすぎ) 7 / 1 ●まず, ê = Y – X̃B = X( β * – β ) + c = (I – X( X̃X̃X̃X̃X̃X̃)εに気をつける. これは平均０の多変量正規分布で，かつβ − β* = (X T X)− X T ε (Eq. (1))と独立である．独立性は, E[ê(β – β*) ] = 0であることと, ê è – β*の同時分布が正規分布であることより，それらは独立であることがわかる(チェックせよ). ●次に， (I − X(X T X)− X T )はn − (d + 1)次元部分空間への射影なので, ||ê||²/σ²は自由度n − (d + 1)のχ²分布に従う. ☻また，最尤推定量の分布より√π(βj − βj )/ √5jj ~ N(0,σ²)である. ●よって， √π(β;−β;)/ √5,布となる．σ²は分母と分子で打ち消す． 口
そもそも回帰に意味がないかもしれない．回帰式に説明力があることを検定.
帰無仮説のもとで,
導出は確率統計第二の講義を参照。また標準的な数理統計の教科書には大体書いてある.決定條数(R²)と自由度調整済み決定條数(R²) :
決定係数が大きければ帰無仮説は棄却され，回帰式には説明力があることになる. 9 / 1
重回帰分析では説明変数を追加するごとに残差は小さくなってゆく.それでいいのか？答え：必ずしも良くはない. 。観測済みデータによく当てはまっても，未観測のデータへの当てはまりが悪くなることがある. 。過適合(overfit)という. ●適適合を避けるにはサンプル数に比して複雑なモデルは使うべきではない.赤池情報量規準, AIC (Akaike Information Criterion)予測精度が一番良いモデルを選択するための規準d次元パラメトリックモデルMaのAICは次式で定義される: AIC = −2 log(ρ({X;};=1| δ Mg )) + 2m.ただし,  θ μ」は M』上の最光推定量.言い換えればAIC = −2 ×最尤推定量の対数尤度+2 ×モデルの次元である． AICが最小になるように説明変数の数を選ぶ. ●次元が高い(複雑な)モデルであるほど負の対数尤度が小さくなるが，それにモデルの次元(複雑さ)を罰則として加えることで過適合を防ぐ. ●モデルの次元とデータへの当てはまりのトレードオフ, ● AICは真の分布とのKL-divergenceを漸近展開した1次項である. ●予測誤差(真とのKL-divergence)を最小化、変数を選択するための客観的規準(予測誤差)を一つ持ってきたところがミソ. (導出は省略)
国土交通省が公開している不動産取引価格情報から世田谷区の中古マンション取引価格データ(平成２５年度第３四半期分)を取得．ここから一部を抜粋したデータで回帰分析をやってみる. http://www.land.mlit.go.jp/webland/download.html従属変数：価格説明変数：最寄り駅からの距離（徒歩）,延床面積，建物の構造，建べい率，容積率，建築年，最寄り駅に急行が止まるか説明変数「建物の構造」は「RC」と「SRC」という文字列(カテゴリカル変数).これを{0,1}に置き換える. ダミー変数という. RC → 1, SRC → 0とした.
書式: Im(formula, data, subset, weights, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, ...)
別: 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以外の変数に回帰
x <- read.csv("setagaya_manshion.csv") sman = data.frame(price=x[[1]],walk=x[[3]],area=x[[5]], str=ifelse(x[[6]]=="RC",1,0),kenpei=x[[7]],youseki=x[[8]], tikunen=x[[9]],kyuko=x[[10]]) #機造はダミー変数0-1に置き換える. plot(sman)   #散布図のプロット
sman.lm <- lm(price ~ area,data=sman)  #回帰分析はこの一行でOK！ plot(sman$area,sman$price, xlab="area",ylab="price")  #結果をプロッL abline(sman.lm , lwd=1 , col="blue")
> summary(sman.lm) Call: lm(formula = price ~ area, data = sman) Residuals: Min       10 Median 30 Max -789511   4806399  16822060 −16082568   −5634345 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1029167    1519818  -0.677      0.5 area            673025      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
Call: lm(formula = price ~ area, data = sman)呼び出し式Residuals: Min 10 Median 30 Max −16082568   −5634345 4806399  16822060 -789511残差y; − y; = y; − β ' xの要約統計量
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1029167    1519818  -0.677      0.5 673025 26408  25.485   <2e-16 *** area --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ 1回帰係数の推定量，標準偏差，t-統計量，p-値(Intercept)は切片. t-統計量の意味は前のスライドを参照.ここのp-植がある有意水準より小さければ，その回帰係数が0であるという帰無仮説は棄却される．つまり，その変数はyを説明するための情報を持っているということ.この例だと，切片項は有意でなく，areaは有意である(たとえば有意水準α = 0.05のとき).
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のほうが重要.大きいほど相関が強い.
par(mfrow=c(2,2)) #一つのFigureに2x2の図を出力plot(sman.lm) (上左) 残差プロット:子潤値ý, vs残差y, - ŷ,残差のQ-Qブロット:正規分布とのQ-Qブロット,対角線に近ければより正規分布に近い。残差の正規性はガウスマルコフモデルの仮定.下左） Scale-Locationプロット:規準化した残差の絶対値の平方根.この値は予測値の値に対して大体一定であるべき.そうでない場合,場所によってノイズの分散が異なることになる. Cook距離のプロット:個々のデータが推定に与える影響を表した距離(そのデータがない場合とある場合の予測値の変化量).大きいと外れ値の可能性がある．0.5を超えると「大きい」とされる.
predict関数にinterval="confidence"なるオプションを入れることで予測値( β T X )の信頼区間を得ることができる. area.plot <- seq(min(sman$area)*0.9,max(sman$area)*1.1,by=1) #信頼区間を計算する範囲の設定sman.con <- predict(sman.lm, data.frame(area=area.plot),interval="confidence") #位钢区間sirca (m 2)
築年数も含めて分析してみる. > sman.lm3 <- lm(price ~ area + tikunen, data=sman) > (AIC(sman.lm3)) [1] 444.881 AIC(۰)でAICが計算できる.すべての説明変数を用いて回帰分析sman.lmall <- lm(price ~., data=sman) sman.lmAIC <- step(sman.lmall) summary(sman.lmAIC) step(-)でAIC最小のモデル(説明変数の組)を探索、walk + area + tikunenの三変数モデルが採用された.
最後に残差の正規性を検定、正規性は棄却されない. > shapiro.test(sman.lmAIC$residuals) Shapiro-Wilk normality test data:  sman.lmAIC$residuals W = 0.9874, p-value = 0.2806
qqnorm(sman.lmAIC$residual) QQプロットからも残差の分布が正規分布に近いことが確認できる.