統計特論２

東京大学大学院農学生命科学研究科　大森宏

３．分散分析

　分散分析は，ANOVA (Analysis of Variance) と略記されることもある．分散分析は，複数の処理を同時に 行ったときに，処理効果を推定するための最も基本的な手法である．データ全体の持つ情報は，総平方和 にまとめられているが，これを，処理の分散成分（処理平均平方）と誤差の分散成分（誤差平均平方）とに 分離して，その大きさを比較することにより，処理の効果を見積もるものである．

因子と水準

　経済学では価格や成長率，工学では作業時間や故障率，農学では収量や抵抗性など， 調査研究したい特性を形質(character)という．着目した形質に影響を与えると考えられるもの， 例えば，収量では品種，温度，施肥量などを要因または因子(factor)という． 要因の影響を調べるためいくつかの品種を用いたり， 施肥量に段階を設けたりするが，それを水準(level)という．

一元配置（one-way layout）

構造モデル

　t 検定では，2 つの処理平均の比較を行ったが，この節ではこれを拡張して，複数の処理平均の比較を行う 手法を考える．いま，a 水準の処理（treatment）A1，…，Aa，があり， 処理 Ai を行った ni 個の標本， Xi1，…，Xini，が得られた とする．処理 Ai からの標本は，平均 μi = μ + αi，分散 σ2 の正規分布に従うと仮定する．ここで，μ を総平均（grand mean）， αi を処理効果（treatment effect）もしくは，主効果（main effect）と言い， Σiαi = 0，である．ここで，平均 0，分散 σ2 を持つ 誤差項（error term）eij を 導入し，標本の構造モデル，
model
として表現すると，データの持つ構造が理解しやすくなる．

平方和分解

　いま，処理 Ai の標本平均を X-i.， 標本総平均を X-.. とすると，これらは，
model
と計算される．標本全体の持つ情報は，総平方和 ST（Total Sum of Squares）で表現される．これは，
model
のように誤差平方和 Se（Error Sum of Squares）と 処理平方和 SA（Treatment Sum of Squares）とに分解される．これは，積の項が
model
のように 0 となるからである．
　なお，誤差平方和を群内平方和（within groups sum of squares），処理平方和を群間平方和（between groups sum of squares）と呼ぶことも多い．

平方和の期待値

　個々の標本 Xij と処理 Ai の標本平均 X-i.， 標本総平均 X-.. の構造モデルがそれぞれ，
model
のようになるので，誤差平方和 Se と処理平方和 SA の期待値は，それぞれ，
model
model
のように計算できる．

帰無仮説のもとでの平方和の比の分布

　一元配置モデルにおける帰無仮説は，すべての処理効果がない，つまり，

H0：αi = 0，i = 1，…，a，　

である．前節の平方和の期待値から，帰無仮説のもとで，Se/σ2 は自由度 n - a の χ2 分布に従い，SA/σ2 は自由度 a - 1 の χ2 分布に 従うことがわかる．これらの χ2 分布をその自由度で割った比の F 値は，
model
のように自由度 a - 1，n - a の F 分布に従う．
　ここで，MA は，処理平方和 SA を その自由度 a - 1 で割ったもので，処理平均平方（treatment mean square）と呼ばれ，処理平均から求めた 誤差分散 σ2 の推定値である．一方，Me は，誤差平方和 Se を その自由度 n - a で割ったもので，誤差平均平方（error mean square）と呼ばれ，誤差分散の推定値である．
　帰無仮説のもとでは MA と Me はほぼ等しいことが期待されるので，その比 F 値は 1 に近いことが期待される．よって，F 値が大きな値をとるときは帰無仮説が正しくないと考え，帰無仮説を 棄却する．F 値が大きいか小さいかの判断基準が対応する自由度の F 分布で決められる．

分散分析表と F 検定

　一元配置モデルの解析結果は，以下の分散分析表（ANOVA table）にまとめられる．

  変動因   自由度（df） 平方和（S.S.） 平均平方（M.S.）   F 値  
主効果 a - 1 SA MA = SA/(a - 1) MA/Me
誤　差 n - a Se Me = Se/(n - a)  
全　体 n - 1 ST    

　この表から検定統計量 F 値が求められる．そして， 自由度 a - 1，n - a の F 分布の 1 - γ 点（例えば 95 ％点）F（a - 1，n - a）1 - γ より F 値が大きい，すなわち，

F ＞ F（a - 1，n - a）1 - γ

であるとき，帰無仮説を棄却すると，有意水準 γ （例えば 5 ％）の検定が行える．これを，F 検定（F test） という．

対比（contarst）

　処理平均 μi のある群と他の群との違いに特に興味がある場合がある．例えば，処理 1，2，3 の 平均と処理 4，5 の平均の間
model
に差があるかをみたいような場合である．一般に，
model
である比較を対比という．

　対比 C = Σiciμi は C^ = ΣiciX-i. で推定されるが， 帰無仮説（αi = 0）のもとでは， 対比推定量の平均と分散は，

model
となるので，分散 σ2 をその推定量 s2 で置き換えた検定統計量 t は，
model
のように自由度 n - a の t 分布に従うので，対比 C = Σiciμi = 0 の 検定を行うことができる．特に，各処理水準の標本数が ni = m と一定で，対比の係数を Σic2i = 1 と標準化すると，検定統計量は，
model
と簡略化される．

　ところで，別の対比 Σidiμi があったときに， Σici di = 0，となるものを直交対比という． 直交対比の組は同時に検定ができる． 　事前に比較が決められるときは，後述する多重比較による有意確率の補正を行わないことが多いと思われる．

多重比較（multiple comparison）

　分散分析（正確には実験計画）の文脈では，試験設計の段階で帰無仮説の設定が行われる．つまり，検定の内容 が事前に決定されている．このような「先付け」のときは，検定の数がそれほど多くないなら， 複数の検定を行っても有意水準についての補正を行わないのが普通だと思われる．
　しかし，データが得られた後，「後付け」でどの処理間の差が有意であるか調べたい誘惑にかられることが多い． 結果として差が 大きかった処理間で t 検定を繰り返して行うと，たくさんの検定を行うので，たまたま有意になる確率が 名目上の有意水準（たとえば 5 ％）を超えてしまう恐れがある．これが，多重比較である．現在では， コンピュータにより多くの検定を簡単に行うことができるので，以前に比べて多重比較の問題を考慮しなければ ならないと考えられる．

　前節の対比でも，考えられる対比を無原則に行うときは多重比較の問題を考慮しなければならないが，この節 では対比の中でも最も単純な対比較（pairwise comparison），すなわち， 個々の処理水準間の比較のみを取り上げる．
　いま，処理平均 μi と μj の比較を行う場合を考える． 2 つの処理平均の差 μi - μj は， dij = X-i. - X-j. で推定される． 帰無仮説（αi = αj = 0）のもとで，dij の平均と分散は，

model
となるので，分散 σ2 をその推定量 s2 で置き換えた 検定統計量 tij は，
model
のように自由度 n - a の t 分布に従うので dij = 0 の検定を行うことができる．
　有意水準 α'（たとえば 5 ％） の検定は， 自由度 n - a の t 分布の 1 - α'/2（たとえば 97.5 ％) 分位点， t(n - a)1 - α'/2 を用いて，
model
が成り立つとき μi と μj の効果に違いがあると判定される．ここで， LSD（Least squared distance）は最小有意差という量で，以前は，α' = 0.05 として，処理効果のある組み合わせ を見つけるためよく用いられていたが，最近は，多重比較を考慮に入れた有意水準の補正を考えるのが普通 なので，単純な LSD は使用しない方が良いと思われる．
　いま，a 水準の主効果があったとすると，すべての組み合わせは r = a(a - 1)/2 通りあり，「後付け」の 検定を行うときは，全体で r 回の検定を行っていると考えなければならない． R でも 対比較では多重比較による有意確率の補正が簡単に行える．

二元配置（two-way layout）

構造モデル

　2 つの因子 A，B に対し，その水準の数をそれぞれ a，b とする．同じ因子と水準の繰り返し （repetition）を r とする．A 因子の第 i 水準で B 因子の第 j 水準の第 k 番目の標本データ Xijk は，
model
とおける．μ は総平均，αi は因子 A の主効果，βj は因子 B の主効果， (αβ)ij は因子 A と B の交互作用（interaction）で，
model
の制約を満たしている．eijk はモデルで説明できない誤差項である．

各種平方和

　因子 A，B の第 i，j 水準の平均，因子 A の第 i 水準の平均，因子 B の第 j 水準の平均および 標本総平均をそれぞれ，
model
とおく．すると，総平方和 ST，因子 A の平方和 SA， 因子 B の平方和 SB， 交互作用平方和 SA×B，誤差平方和 Se はそれぞれ，
model
と計算される．1 元配置分散分析のときと同様に，
ST = SA + SB + SA×B + Se
という平方和の分解ができる．

分散分析表と F 検定

　二元配置モデルの解析結果は，以下の分散分析表（ANOVA table）にまとめられる．

  変動因   自由度（df） 平方和（S.S.） 平均平方（M.S.）   F 値  
主効果 A a - 1 SA MA = SA/(a - 1) MA/Me
主効果 B b - 1 SB MB = SB/(b - 1) MB/Me
交互作用 (a - 1)(b - 1) SA×B MA×B = SA×B/(a - 1)(b - 1) MA×B/Me
誤　差 ab(r - 1) Se Me = Se/(n - a)  
全　体 n - 1 ST    

実験計画法

に続く．

実験計画法（追加）

４．回帰分析

相関

　標本（サンプル）に対し，２つの変数 x，y が測定されているとする． たとえば，x が身長(m)であり，y が体重(kg)である． 大きさ n の標本（サンプル）に対し，２つの変数の組のデータが，

（x1 ，y1 ）， （x2 ，y2 ）， …， （xn ，yn ）
であったとする． 変数間の関連性の強さを測る量として共分散（Covariance），Cov[x ，y ] がある． これは，変数に対する平均を，
として，
と定義される．

　共分散は測定単位により大きさが変わるので，これをおのおのの変数の標本分散， Var[x]，Var[y]，

で標準化したものがピアソン（Peason）の積率 相関係数 r，-1≦r≦1，であり，

sokan

と定義される．これは，変数間の線形的関係の強さ， （x が大きいと y も大きく，x が小さいと y も小さい，） を測る指標で，｜r｜＝1 のときは， 変数 x，y は完全な直線関係にあり，r ＝0 のときは，線形的な関係 がない．r が 1 に近いときは，正の相関関係があるといい， r が -1 に近いときは，負の相関関係があるという．

データ散布図と相関係数
sokan sokan
sokan sokan

相関データダウンロード

# 相関係数の R スクリプト
r <- read.csv("r1.csv") # csv データ読み込み
x <- r[,1]; y <- r[,2] #
plot(x,y, xlim=c(2,8), ylim=c(0,6)) # x，y の散布図
title(main="相関係数 r = 0.23") #
cov(x, y) # x と y の共分散
cor(x, y) # x と y　の相関係数
x <- 10*x #
plot(x,y, xlim=c(2,8), ylim=c(0,6)) # x，y の散布図
cov(x, y) # x と y の共分散
cor(x, y) # x と y　の相関係数

直線回帰

モデル

　２つの変数 x ，y に対し，y の値が x の値の動きにつれて 線形的に変化すると仮定される，つまり，

y ＝ a ＋ b x

という関係が成り立っていると考えられる場合である．これを y の x に 対する直線回帰といい，a ，b を回帰係数という． また，変数 y を従属変数，目的変数といい，変数 x を独立変数，説明変数と いう．

最小２乗法

　データに最もよくあてはまる直線回帰式を得るには，データ点 （xi ，yi ）， と回帰による推定点， （xi ，y^i ）， y^i ＝ a ＋ b xi ， の間の距離の２乗和 S が最小になるような回帰係数 a ，b を求める．つまり，

minsqure

を最小化する a ，b を求める問題に帰着する．これを最小２乗法という．

　これは，S を a ，b で偏微分して 0 とおくことによって得られる．つまり，

partial

の連立方程式を a ，b で解けばよい．これより，回帰係数の推定値 a^， b^ が，

coef

と得られる．

回帰式の統計モデル

　推定された直線回帰式がどの程度現実のデータに適合しているかを調べるために， 回帰式が従う統計モデルを考える．標本の格データ点， （xi ，yi ）， が，

yi ＝ a ＋ b xi ＋ ei ， ei 〜 N( 0，σ2 )

であると仮定する．ei は誤差（error），あるいは， 残差（residual）で，直線回帰 式では説明がつかない部分を表し，これが互いに独立に平均 0，分散 σ2 の正規分布に従うと仮定する．誤差の大きさが大きいときは，直線回帰式ではデータが説明できない と考える．

残差分散と回帰係数の標準誤差

　回帰で説明がつかない残差平方和 Se は，

minsqure
で求められる．これの自由度は n−2 であるので（２つの回帰係数分の自由度を除く），回帰の 残差（誤差）分散の推定値は，
minsqure
で求められる．

　一般に，Var(yi ) ＝ σ2 であるとき，その定数（c） 倍の分散は，

Var(cyi ) ＝ c2σ2， Var(Σici yi ) ＝ Σici 2 σ2

であり，従属変数 y のデータ yi は，

yi 〜 N( a ＋ b xi ，σ2 )

と分布するので，回帰係数の推定値 b^ の分散は，
minsqure
となる．この分散の平方根を回帰係数推定値 b^ の標準誤差という．

回帰係数の標準誤差による t 検定

目的変数 y が説明変数 x との回帰関係にないという 帰無仮説，

H0：b ＝ 0，

を考えてみよう． 回帰係数 b の推定値 b^ の分散は，
minsqure
と推定できるので，b^ の標準偏差（標準誤差）は， s b と推定 される．これより，回帰係数をその標準誤差で割った t 値が，帰無仮説のもとで，
minsqure
のように，自由度 n−2 の t 分布に従うことを利用して回帰係数の検定が行える．すなわち， 自由度 n−2 の t 分布の 97.5％点を t0 とすると，

|t | ＞ t0 → 帰無仮 説を有意水準 5 ％で棄却（回帰関係が有意に認められる）
|t | ≦ t0 → 帰無仮説を棄却しない（回帰関係が認められない）

と定式化できる．

分散分析

平方和分解

　回帰式により， 従属変数 y のデータ yi は，

yi ＝ y^i ＋ (yi − y^i ) ＝ 回帰値 ＋ 残差

のように分解される．この分解に対応して従属変数データの総平方和 ST は，

ST ＝ Σi （y i − y- ) 2 ＝ Σi （y^i − y- ) 2 ＋ Σi （y i − y^i ) 2 ＝ SR ＋ Se
総平方和 ＝ 回帰平方和 ＋ 残差平方和

のように分解される．これを平方和の分解という．この分解に対応して自由度は，
n−1 ＝ 1 ＋ n−2
と分解される．

決定係数（重相関係数の２乗）

　データが直線回帰式でよく説明できるのは，回帰平方和が大きく，残差平方和 が小さい場合である．総平方和のうち回帰平方和で説明される割合を決定係数，もしくは 重相関係数の２乗といい，

minsqure
で定義される．なお，重相関係数 R とは，データ y i と回帰値 y^i との間の相関係数である．これより， 以下の分散分析表ができる．

回帰分析の分散分析表
変動因平方和自由度平均平方F 値
回帰 SR 1 SR F ＝ SR/se2
残差 Se n−2 se2 ＝ Se/n−2 　
全体 ST n−1 　 　

F 検定

従属変数 y が説明変数 x の回帰関係にないという 帰無仮説，

H0：b ＝ 0，

を考える．帰無仮説のもとでは，回帰平均平方 SR と残差分散 se2 がともに誤差 σ2 の不偏推定量になるので， その比 F 値が，

F ＝ SR/se2 〜 F（1，n−2），

という F 分布に従うことを利用して検定ができる．すなわち，分子，分母 自由度が 1，n−2 である F 分布 F（1，n−2）の95％点を F0 とすると，

F ＞ F0 → 帰無仮説を有意水準 5 ％で棄却（回帰関係が有意に認められる）
F ≦ F0 → 帰無仮説を棄却しない（回帰関係が認められない）

と定式化できる．

# 回帰分析の詳細の R スクリプト 
anova(mbaf.lm)  			# 分散分析表  
n <- nrow(mbaf); n  		# データ数  
st <- var(mbaf[,3])*(n-1); st  #　総平方和 
sr <- var(fitted(mbaf.lm))*(n-1); sr  #　回帰平方和 
se <- var(resid(mbaf.lm))*(n-1); se  #　残差平方和 
s <- sqrt(se/(n-2)); s  		#　残差標準偏差 
sx <- var(mbaf[,2])*(n-1); sx  	#　x の総平方和 
b <- cov(mbaf[,2],mbaf[,3])*(n-1)/sx; b  #　回帰係数 
sb <- s/sqrt(sx); sb  		#　回帰係数bの標準誤差 
r2 <- sr/st; r2  		# 重相関係数の２乗 
r <- cor(mbaf[,3], fitted(mbaf.lm)); r  # 目的変数と回帰推定値の相関 
fv <- sr/(se/(n-2)); fv  		# F値 
tv <- b/sb; tv  			# t値

回帰式の信頼区間

回帰係数の信頼区間

　回帰係数の標準誤差 s b を用いて， 回帰係数 b の信頼区間がつくれる．すなわち， 自由度 n−2 の t 分布の 97.5％点を t0 とすると， 回帰係数 b の 95％信頼区間の幅 d は，d ＝ t0 s b となるので， 95％信頼区間は，

b^ − t0 s b ＜ b ＜ b^ ＋ t0 s b

となる．
回帰直線の信頼区間
　データから推定された回帰直線は，データの平均（x-，y-） を通るので，

Y = y- + b^ （x - x-）

とおける．これより，Y の分散は，
minsqure
となる．誤差分散 σ2 は， データの残差分散 se2 で推定されるので，Y の標準誤差は，
minsqure
となり，これが自由度 n - 2 を持つ．よって，推定回帰式にこの標準誤差の t0 倍を加えたものが，回帰式の 95％信頼幅となる．
回帰予測値の信頼幅
　回帰式から得られる予測値 y~ は，回帰式に誤差項が加わって，

y~ = y- + b^ （x - x-）+ e

となるので，その分散は，
minsqure
となる．先ほどと同様に，誤差標準偏差を残差標準偏差で置き換えると，回帰予測値 y~ の 標準誤差は，
minsqure
となり，回帰式の 95％信頼幅の外側に回帰予測値の 95％信頼幅が描ける．
minsqure

# 回帰式の信頼幅の R スクリプト 
new <- data.frame(Entexam=seq(400,700,by=2))  		# 予測したい範囲の定義  
Grac <- predict(mbaf.lm, new, interval="confidence", level=0.95)  # 回帰推定値（回帰直線）の 95 ％信頼幅  
Grap <- predict(mbaf.lm, new, interval="prediction", level=0.95)  # 回帰予測値の 95 ％信頼幅  
matplot(new$Entexam, cbind(Grac,Grap[,-1]), lty=c(1,2,2,3,3), type="l", 
col=c("blue","blue","blue","red","red"), xlab="入試得点", ylab="初年度成績")  
points(mbaf[,2:3])  					# データの表示 
title(main="回帰式と予測値の 95 ％信頼幅")  #
# 信頼幅を詳しく計算 
n <- nrow(mbaf); n  				# データ数  
sx <- var(mbaf$Entexam)*(n-1); sx  		# x の偏差平方和  
se <- summary(mbaf.lm)$sigma; se  		# 残差標準偏差  
b <- coef(mbaf.lm)[2]; b  			# 回帰係数  
mx <- mean(mbaf$Entexam)  			# x の平均  
my <- mean(mbaf$Grade1st)  			# y の平均  
t0 <- qt(0.975, df=(n-2)); t0  		# 自由度 n - 2 の t 分布 97.5％点  
x11()				# 新しいグラフウィンドウで表示
plot(mbaf[,2:3], xlim=c(400,700), ylim=c(2,4))  # データ散布図  
abline(mbaf.lm, col="blue")  				# 回帰直線  
sr <- 1/n + (new$Entexam-mx)^2/sx  #  
y1 <- my + b*(new$Entexam-mx) + t0*se*sqrt(sr)  # 回帰直線 95％信頼幅上限  
y2 <- my + b*(new$Entexam-mx) - t0*se*sqrt(sr)  # 回帰直線 95％信頼幅下限  
points(new$Entexam, y1, type="l", lty=2, col="blue")  #  
points(new$Entexam, y2, type="l", lty=2, col="blue")  #  
yp1 <- my + b*(new$Entexam-mx) + t0*se*sqrt(sr+1)  # 回帰予測値 95％信頼幅上限  
yp2 <- my + b*(new$Entexam-mx) - t0*se*sqrt(sr+1)  # 回帰予測値 95％信頼幅下限  
points(new$Entexam, yp1, type="l", lty=3, col="red")  #  
points(new$Entexam, yp2, type="l", lty=3, col="red")  #  
title(main="回帰式と予測値の 95 ％信頼幅")

課題：回帰式と予測値の 99 ％信頼幅をつくれ．

重回帰分析

統計モデル
　説明変数が 2 つ以上になった場合である．いま，p 個の説明変数 x1， x2，…，xp，により目的変数 y が，

y = b0 + b1x1 + b2x2 + … + bpxm + e

と表現できるとする．ここで，b0b1，…，bp は 回帰係数，e は，誤差である．
　ここで，n（n ＞ p）個のデータがあり，目的変数が y = (y1，…，yn)'，j 番目の説明変数 xj の値が， x1 = (x1j，…，xnj)' であったとする． 回帰係数ベクトルを b = (b0，…，bm)'，誤差ベクトル を e = (e1，…，en)' とすると，データが満たす構造 モデルは，
minsqure
minsqure
となる．ここで，誤差ベクトルは，平均 0，分散共分散行列 σ2I の n 次元 正規分布に従うとする．これは，n 個の誤差 ei が互いに独立に平均 0 分散 σ2 の正規分布に従うと考えたとき（普通行う仮定）の多変量表現である．
　このモデルのもとでの残差平方和は，
mvreg1
となる．最小２乗法は，S を最小にする回帰係数ベクトル b を求めることである． これは，S を b で偏微分して 0 とおくと，
mvreg2
となるので，まとめると，正規方程式
mvreg3
が得られる．n > p であるときは，この方程式の解として，回帰係数の推定値は，
minsqure
と求められる．
多変量の線形関数の平均と分散
　p 変量確率変数 x の平均を μ，分散共分散行列を Σ とすると，任意の p×p スカラー 行列 A に対し，

E[Ax] = AE[x] = Aμ， Var[Ax] = AVar[x]A' = AΣA'，

が成り立つ．
回帰係数の分散
　重回帰モデルの回帰係数ベクトルの推定量 b^ の分散は，
mreg
となる．残差分散 σ2 の推定値を se2 とすると， 回帰係数の分散は，se2(X'X)-1 で推定される．
ハット行列とてこ比
　さて，回帰推定値ベクトル y^ は，
と書ける．ここで，H は，説明変数データベクトル y から 回帰推定値ベクトル y^ を生成する行列なので，ハット行列（hat matrix）と 呼ばれている．ハット行列の対角成分 hii が大きいと i 番目のデータは推定値に大きな 影響を及ぼすので，これを「てこ比（leverage）」と呼んでいる．てこ比の平均的な値は，(p + 1)/n なので， たとえば，2(p + 1)/n より大きなてこ比をもつようなデータには注意を払う必要がある．すなわち， このデータの存在が回帰係数の推定に大きな影響を与えているからである．
標準化残差
　回帰残差推定量ベクトル e^ は，
とハット行列を用いて書ける．これより，その分散は，
となる．なお，I - H がべき等行列であることに注意せよ．従って，個々の 回帰残差推定量 e^i の分布は，
となる．この ri を標準化残差（standarized residual）と言う．
Cook 距離
　個々のデータが回帰係数の推定値に与える影響を調べるため，i 番目のデータを除いたときの回帰係数 推定値を b^-i として，これが全データを用いたときの回帰係数 推定値 b^ との違いの大きさで測る．これが Cook 距離で，以下のように 定義される．
mreg
ここで，y^-i = Xb^-i で ある．また，多少の計算により Cook 距離は標準化残差とてこ比で計算され，最後の式のように表せる．
Cook 距離が 0.5 を超えるとそのデータは影響が「大きい」とされ，1 を超えると影響が「特に大きい」と される．

回帰診断

　回帰分析を行ったときに，データが回帰モデルによくフィットしているかや，データの中に回帰モデル からはずれたも（異常値）がないか，などを調べた方がより安全である．これを 回帰診断（regression diagnostics）という． R などの統計ソフトが普及する以前は，手間がかかるので回帰診断まで行うことはあまりやらなかったが， 現在では手軽にできるので，行うのが普通になってきていると思われる．このため，回帰診断の考え方や 着目点などを理解する必要がでてきたと言える． 　回帰診断は，以下の２点からなる．
diag diag
diag diag
# 回帰診断の R スクリプト 
mba <- read.csv("mbagrade.csv")  		# データ読み込み　 
mbaf <- mba[mba[,1]=="F",]  		# 女性データのみ抽出 　 
mbaf.lm <- lm(Grade1st ~ Entexam, data=mbaf) 			# 回帰
summary(mbaf.lm)  				# 結果表示 　 
plot(mbaf.lm)				# 回帰診断はたったこれだけ．（４枚のグラフが出てくる）
#
# 回帰診断を詳しくみる
# てこ比
lev <- hatvalues(mbaf.lm); lev		# てこ比
X <- model.matrix(mbaf.lm)			# 説明変数行列
H <- X %*% solve(t(X) %*% X) %*% t(X)	# ハット行列
diag(H)					# てこ比の計算による導出
# 標準化残差
rstd <- rstandard(mbaf.lm); rstd		# 標準化残差
se <- summary(mbaf.lm)$sigma		# 残差標準偏差
resid(mbaf.lm)/(se*sqrt(1-lev))		# 標準化残差の計算による導出
# Cook 距離
cookd <- cooks.distance(mbaf.lm); cookd	# Cook 距離
b <- coef(mbaf.lm); b 		# 回帰係数
lm1 <- lm(Grade1st[-1] ~ Entexam[-1], data=mbaf)
b1 <- coef(lm1); b1 			# １番目のデータを除いた回帰係数
t(b-b1) %*% t(X) %*% X %*% (b-b1)/se^2/2	# １番目のデータの Cook 距離
rstd^2*lev/(1-lev)/2			# Cook 距離の計算による導出
# 回帰推定値 - 残差プロット
x0 <- fitted(mbaf.lm)			# 回帰推定値
y0 <- resid(mbaf.lm)			# 回帰残差
plot(x0, y0, type="n", xlab="回帰推定値", ylab="残差")	# 
text(x0, y0, rownames(mbaf), cex=0.8)
abline(h=0, lty=3)
title(main="回帰推定値 - 残差プロット")
# 標準化残差 Q - Q プロット
a <- qqnorm(rstd, type="n", main="")		# 標準化残差の Q - Q プロット（正規分布からのずれ）
text(a, rownames(mbaf), cex=0.8)
qqline(rstd, col="red")
title("標準化残差正規 Q-Q プロット")
# S - L プロット
yr <- sqrt(abs(rstd))			# 標準化残差絶対値の平方根
plot(x0, yr, type="n", xlab="回帰推定値", ylab="残差絶対値の平方根")	# 
text(x0, yr, rownames(mbaf), cex=0.8)
title(main="S - L プロット")
# てこ比 - 標準化残差 プロット
plot(lev, rstd, type="n", xlim=c(0,0.5), xlab="てこ比", ylab="標準化残差")
text(lev, rstd, rownames(mbaf), cex=0.8)
abline(h=0, lty=3)
abline(v=0, lty=3)
title(main="てこ比 - 標準化残差 プロット")
xx <- seq(0, 0.5, by=0.01)
yc <- .5*2*(1-xx)/xx			# Cook 距離 0.5 となる標準化残差
yc1 <- sqrt(yc)
yc2 <- -yc1
points(xx, yc1, type="l", lty=2, col="red")		# Cook 距離 0.5 線
points(xx, -yc1, type="l", lty=2, col="red")
課題：47 番のデータは，回帰モデルから大きく離れていることがわかった．このデータを 異常値として取り除いた回帰分析を行なえ．

重回帰分析の例

モデル間の尤度比検定

モデル選択と AIC

　確率モデルのパラメータ推定には，通常，最尤法が用いられる．しかしながら，重回帰分析などで説明変数の 個数（パラメータ数）を決めようとすると，一般に，パラメータ数が多いほどデータへのモデルの当てはまり （fitting）が良くなるので，最尤法でパラメータ数を決めるとパラメータ数の多いモデルが「良い」とされて しまう．パラメータ数の多いモデルは，パラメータの値を推定したデータにはよく当てはまるが， 同様の状況から得られた別のデータへの当てはまりが悪くなることが知られている．このような現象を解釈 しすぎ（over fitting）という．
　これを避けるには，できるだけ単純なモデルを考えるのがよいとされている．これを実現するモデル選択 の基準として，

X = (モデルのデータへの当てはまり) + (モデルの複雑さへのペナルティ)

の形式のものがいくつか提案されている．この中で有名な基準の一つが AIC (Akaike Information Criterion) である．AIC は，

AIC = - 2×(モデルの最大対数尤度) + 2×(モデルの自由パラメータ数)

と定義される．モデルの最大対数尤度は，確率モデルの最尤推定値を確率モデルに代入したときの尤度の対数 を取ったものであり，モデルのデータへの当てはまりのよさを評価している．モデルの自由パラメータ数 は，モデルの複雑さの尺度の一つで，パラメータ数の少ないモデルほど単純でよいものと考えられる． 結局，AIC の小さなモデルがよいとされる．

　k 個のパラメータ θ を持つ回帰モデル

y = f(x ; θ) + e
において，残差 e が正規分布に従うモデルでは，n 個のデータから得られた 残差分散の最尤推定値を v2 とすると， 回帰モデルの最大対数尤度は，
l = -(1/2) [n log 2π + n log v2 + n ]
となる．これより，回帰モデルでの AIC は，
AIC = (n log 2π + n log v2 + n ) + 2k
となり，AIC の小さな回帰モデルがよいとされる．

説明変数の選択

すべての組み合わせの AIC
　全変数を使って重回帰分析を行ったが，経度（x2）の回帰係数の有意確率が小さくないので， 経度の情報は気温を説明するのに必要ないかも知れない．これは日本では，緯度が高く標高が高いほど気温 が低いと考えられることとも一致している．モデル選択の方法として AIC を利用してみる．
　3 個の説明変数があるので，説明変数の組み合わせは 23 = 8 通りある．このすべての組合わせ に対して AIC の値を計算し，最も小さな値をもつモデルを採用することにする．

# 変数選択の R スクリプト
kion0.lm <- lm(気温 ~ 1, data=kion) # 説明変数無し回帰
anova(kion0.lm) #
kion11.lm <- lm(気温 ~ 緯度, data=kion) # 説明変数：緯度，回帰
anova(kion11.lm) #
kion12.lm <- lm(気温 ~ 経度, data=kion) # 説明変数：経度，回帰
anova(kion12.lm) #
kion13.lm <- lm(気温 ~ 標高, data=kion) # 説明変数：標高，回帰
anova(kion13.lm) #
kion21.lm <- lm(気温 ~ 緯度 + 経度, data=kion) # 説明変数：緯度，経度，回帰
anova(kion21.lm) #
kion22.lm <- lm(気温 ~ 経度 + 標高, data=kion) # 説明変数：経度，標高，回帰
anova(kion22.lm) #
kion23.lm <- lm(気温 ~ 緯度 + 標高, data=kion) # 説明変数：緯度，標高，回帰
anova(kion23.lm) #
kion3.lm <- lm(気温 ~ 緯度 + 経度 + 標高, data=kion) # 説明変数：緯度，経度，標高，回帰
anova(kion3.lm) #
AIC(kion0.lm) # 説明変数無し回帰 AIC
AIC(kion11.lm) # 説明変数：緯度，回帰 AIC
AIC(kion12.lm) # 説明変数：経度，回帰 AIC
AIC(kion13.lm) # 説明変数：標高，回帰 AIC
AIC(kion21.lm) # 説明変数：緯度，経度，回帰 AIC
AIC(kion22.lm) # 説明変数：経度，標高，回帰 AIC
AIC(kion23.lm) # 説明変数：緯度，標高，回帰 AIC
AIC(kion3.lm) # 説明変数：緯度，経度，標高，回帰 AIC

変数増減法の AIC
　説明変数の数が多くなり，すべての組み合わせを調べることが大変な場合には，変数増加法，変数減少法， その組み合わせである変数増減法がある．ここでは，変数増減法もやってみる．

# 変数選択の R スクリプト
library(MASS) # MASS ライブラリィー読み込み
null <- lm(気温 ~ 1, kion) # 説明変数無し
full <- lm(気温 ~ 緯度 + 経度 + 標高, kion) # 説明変数３つ
result <- stepAIC(null, scope=list(lower=null, upper=full), data=kion) # 変数増減法
summary(result) # 結果表示

多項式回帰

　目的変数 y と説明変数 x との関係が直線関係では説明できないようなときには， x の多項式で回帰することを考える．p 次までの多項式を考えると，

y = b0 + b1x + b2x2 + … + bpxp + e

というモデルとなる．ここで，x2 = x2，…， xp = xp と変数変換を行えば通常の重回帰分析と同じである．

ブートストラップ法による信頼区間の構成

　回帰残差が正規分布するというモデルのもとでは，回帰係数推定量が t 分布に従う ことを利用して，回帰係数の信頼区間を構成することができた．
　回帰残差の正規性が疑われる場合や，正規性が成立しない場合 の信頼区間を構成する方法の一つとしてブートストラップがある．ブートストラップ法は， 大きさ n の標本（サンプル）から復元抽出 を許して新しく大きさ n の標本を生成する（リサンプリング）手法で，生成された標本ごとに回帰係数の 推定値が得られる．このブートストラップサンプルによる回帰式はサンプルごとに微妙に異なる．
　1000回程度のリサンプリングを行えば，回帰係数が1000個得られるので，この分布から 回帰係数の信頼区間を構成することができる．

多変量解析

に続く．

冬休みレポート

ミヤマクワガタの相対成長解析

データダウンロード

成長段階の異なる 47 頭のミヤマクワガタのパーツ別の重量データ(g)．パーツは， 頭部（WHEAD），前胸部（WTHORAX），中胸〜腹部（WABDOM），交尾器（WGENI）からなる．

　形態形質 x，y に対して，初期成長指数を a，成長比を b と したときの生物形態の相対成長式，
logy = loga + b logx
において，中胸〜腹部（WABDOM）の重量(g)を体サイズの指標 x としたとき，y として 頭部（WHEAD）の重量(g)，前胸部（WTHORAX）の重量(g)，交尾器（WGENI）の重量(g) の３通りを考えたときの相対成長式の回帰係数をそれぞれ求める．
　体の成長に合わせてパーツも比例して大きくなれば，成長比 b = 1 となるはずである． 以下の問に答えよ．

参考文献（古い順）

  1. Introduction to the Theory of Statistics, Mood, A. M., Graubill, F. A. & Boes, D. C., 1974, McGRAW-HILL
  2. 工学のためのデータサイエンス入門−フリーな統計環境Rを用いたデータ解析−，間瀬茂ら，2004， 数理工学社
  3. 実践生物統計学−分子から生態まで−（第 1 章，第 2 章）， 東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
  4. The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
  5. R で学ぶデータマインニング I −データ解析の視点から−，熊谷悦生・船尾暢男，2007，九天社
  6. R で学ぶデータマインニング II −シミュレーションの視点から−，熊谷悦生・船尾暢男，2007，九天社
  7. 生物統計学入門，上村賢治・高野泰・大森宏，2008，オーム社
Copyright (C) 2008, Hiroshi Omori. 最終更新日：2013年4月16日