図1 modelについての診断プロット
<< 目次へ
最終更新日 ： 2016.9.17
計数データとは言いかえれば個々の要素をカウントしたデータであるといえる。例えば、人数や件数といったものがこれである。今まで扱ってきた連続型のデータと大きく異なる点を以下に列挙してみよう（連続量に対して計数データは離散量と呼ばれている）。
以上の2点を考慮すると、計数データを連続型として扱った場合にうまくない結果が得られてしまう。それは小数点を含むデータが結果として得られるが、そのような結果は事実上在りえない結果である。もう1つはマイナスの値を結果として得られてしまうことが在りえるが、これも既に述べたように実際上在りえないことである。
また統計学の理論として、計数データは誤差が正規分布に従わないということが挙げられる。したがって、もし計数データを解析する際に「誤差が正規分布に従う」と仮定して解析してしまうと、その解析結果の妥当性は低いものとなる（というか、そもそも不適切である）。
ところが、実際にいくつかの論文を見てみると、計数データであるにもかかわらず重回帰モデルを当てはめて解析していることがある。これはどういうことであろうか。結論をいってしまうと、重回帰分析は目的変数が間隔尺度のデータ（人数とか件数は間隔尺度である）であれば問題ないのである。いってみれば便宜的に重回帰モデルを当てはめて解析しているわけである。これは、5件法で回答を求めた質問紙を分析する際に、本来的には順序尺度であるのだが、間隔尺度として分析することに似ているかもしれない。
さて、ここではまず真のカテゴリカル型データとはどのようなものであるかを説明しておこう。以前に「カテゴリカル型データは、分類のために適当な数字を割り当てたものにすぎない」と述べた。例えば、小学校を1、中学校を2、高校を3といったように任意の数字を割り当て、これを水準と呼んだのであった（水準はカテゴリ、属性、レベルなどともいわれるのであった）。これも間違いなくカテゴリカル型のデータであるが、同じカテゴリカル型データといっても計数データとなると少し定義が異なるのである。
少し混乱してしまうかもしれないが、一言でいえばカテゴリカル型である計数データとは独立な個々の要素の計数である、ということである（だから、重複してカウントされたデータはもはや計数データとはいえない）。そして、連続型であるデータは正規分布に従うと仮定された。それに対してカテゴリカル型データが従うと仮定される分布はポアソン分布や二項分布である。前者のポアソン分布は今回紹介する計数データであり、後者の二項分布は次章で紹介する2値データの解析で用いられるものである。
先ほどから計数データと舌打っているが、要するにこれは分割表の解析である。古典的には分割表でまとめられたデータはカイ自乗検定（独立性の検定）やフィッシャーの正確確率検定（Fisher's Exact Test）によって分析される。これについては一般線形モデル（General Liner Model）を拡張した一般化線形モデル（Genelalized Liner Model）を用いて解析することが適切であるが、一般線形モデルを用いて解析することもできる。
分割表でまとめられたデータを一般線形モデルとして扱うためには、
ということを理解しておく必要がある。ここでいういくつかの問題とは、いうまでもなく誤差の正規性や分散の均一性のことである。
試しに表1のデータを一般線形モデルを用いて解析してみよう。表1のデータを一般線形モデルとして扱うためには、表2のようにデータを置き換えてやる必要がある。
| b1 | b2 | b3 | 合計 | |
| a1 | 16 | 12 | 36 | 64 |
| a2 | 12 | 5 | 20 | 37 |
| a3 | 15 | 11 | 24 | 50 |
| a4 | 9 | 2 | 1 | 12 |
| 合計 | 52 | 30 | 81 |
| A | B | DAT |
| 1 | 1 | 16 |
| 1 | 2 | 12 |
| 1 | 3 | 36 |
| 2 | 1 | 12 |
| 2 | 2 | 5 |
| 2 | 3 | 20 |
| 3 | 1 | 15 |
| 3 | 2 | 11 |
| 3 | 3 | 24 |
| 4 | 1 | 9 |
| 4 | 2 | 2 |
| 4 | 3 | 1 |
表2のデータについて、従属変数をDAT、独立変数をAとBとして解析してみる（式1）。
DAT = A + B ・・・ （式1）
> model <- lm(DAT ~ A + B)
> summary(aov(model))
Df Sum Sq Mean Sq F value Pr(>F)
A 3 488.92 162.97 4.0268 0.06922 .
B 2 327.17 163.58 4.0419 0.07732 .
Residuals 6 242.83 40.47
すると、A要因はDATに影響を及ぼしているとはいえない（p = 0.06922）、B要因もDATに影響を及ぼしているとはいえない（p = 0.07732）。しかし、この結果は妥当なものといえるだろうか。それを確かめるためには標準化残差プロットやQ-Qプロットを見てみればよい（図1）。
> par(mfrow=c(2,2)) # プロットウィンドウを4分割 > plot(model, which=1:4)
図1 modelについての診断プロット
左下の図（Scale-Location）が標準化された残差のプロットである。簡単におさらいすると、この標準化残差プロットは、描かれているデータポイントが均一にばらついていることが望ましい。もし横軸の左側にデータポイントが集中しており、横軸の右側（横軸の値が大きくなるにつれて）データポイントが広く散布していくようだと、もはや分散の均一性が保たれているとはいえない。今回の場合は何とも言えないような状態である。
一方、右上の図（Normal Q-Q）は誤差の正規性を表すものである。点線で描かれている直線上にピッタリとデータポイントが乗っかっている状態が好ましい。今回のように明らかに直線状から外れているデータポイントが多い場合は誤差の正規性が保障されていないと判断できる。
続いて、今度は変数DATを平方根変換して解析した場合をみてみよう（式2）。
sqrt(DAT) = A + B ・・・ （式2）
> model2 <- lm(sqrt(DAT) ~ A + B)
> summary(aov(model2))
Df Sum Sq Mean Sq F value Pr(>F)
A 3 12.3656 4.1219 4.9055 0.04699 *
B 2 4.5566 2.2783 2.7115 0.14492
Residuals 6 5.0415 0.8403
先ほどの解析ではA要因は有意でなかったが、今度の解析では有意である（p = 0.04699）。B要因については今度の解析においても有意ではなかった（p = 0.14492）。この結果、すなわち平方根変換を行ったことによって、分散の均一性や誤差の正規性がどれほど改善されたかを見てみよう（図2）。
> par(mfrow=c(2,2)) > plot(model2, which=1:4)
図2 model2についての診断プロット
分散の均一性については図1の左下の図と比べても明らかに改善されたというようには見えない（もともと分散の均一性については問題がなかったといえる）。右上のQ-Qプロットは図1と図2とを見比べると、図2（すなわち平方根変換をした場合）では改善されているのが見てとれる。図2のほうでは明らかに直線上付近にデータポイントが集まっているのが確認できる。
以上のことから、平方根変換を行うことによって分散の均一性と誤差の正規性の問題を解決できるということが分かっただろう。ただし、データがポアソン分布に従っているならば、式2についての分散分析表よりも検出力の高い検定ができるので、その方法を紹介することにしよう。
その方法とは、変数AとBによって説明される平方和に基づくポアソン分散検定と呼ばれるものである。ポアソン分散検定に用いられる検定等計量であるχ2値は式3によって計算することができる。
χ2 = 処理の自由度 * 処理の平均平方（MS） / 0.25 ・・・ （式3）
実際に計算してみると、
Aに対するχ2 ：
χ2 = 3 * 4.1219 / 0.25 = 49.4628
Bに対するχ2 ：
χ2 = 2 * 2.2783 / 0.25 = 18.2264
となる。このχ2による片側検定のp値はRの関数pchisq()を用いて求めることができる。
> 1 - pchisq(49.4628, 3) [1] 1.039674e-10 > 1 - pchisq(18.2264, 2) [1] 0.0001102015
すると、A要因についてはp = 0.0000、B要因については0.0001となるのでいずれの要因も有意であると判断することができる。つまり、AとBは共にDATに影響を及ぼす要因であるといえるわけである。いうまでもなく、このような解析では独立変数（AとB）が従属変数（DAT）の変動にどれほど関与しているかという主張ができるのであって、独立性の検定のように「AとBには関連がある」というような主張をするのとは少し異なるものとなる。
独立性の検定とは、見方を変えれば「2変数の交互作用を検定する」ことである。もし仮にAとBとの間に有意な交互作用が認められるのであれば、Aという条件のもとでBはDATに影響を及ぼす（逆にBという条件のもとでAはDATに影響を及ぼす）といえるので、結果としてこれはAとBには関連があると主張していることになるというわけである。
かといって、交互作用項を含めた
DAT = A + B + A*B
といったモデル解析をすることとは別問題である。確かに交互作用項を検定していることにもなるが、そもそも分割表でまとめられたデータでは交互作用項を含めたモデルを解析することはできない。なぜなら、一言でいってしまえばデータ数（より正確にいえば反復数）が少ないため計算ができないからである。
ところで、分割表でまとめられたデータを一般線形モデルとして扱うということは各要因の主効果を検定するということであった。これを古典的な独立性の検定によって再現することもできる。そのためには、表1にある合計のセルにある度数を利用すればよい。
> a <- c(64, 37, 50, 12) # Aの合計のセルにあるデータ
> b <- c(52, 30, 81) # Bの合計のセルにあるデータ
> chisq.test(a) # aについて検定
Chi-squared test for given probabilities
data: a
X-squared = 35.9939, df = 3, p-value = 7.511e-08
> chisq.test(b) # bについて検定
Chi-squared test for given probabilities
data: b
X-squared = 24.0859, df = 2, p-value = 5.886e-06
これは式2のモデル解析の結果と同じ結果を導くことになる。ただし、計算方法が違うので計算されたp値は異なる。
| 式2のモデル解析 | A | p = 1.039674e-10 = 0.0000000001039674 |
| B | p = 0.0001102015 | |
| 合計セルの独立性の検定 | A | p = 7.511e-08 = 0.00000007511 |
| B | p = 5.886e-06 = 0.000005886 |
また、このような分割表にまとめられたデータ（計数データ）を解析するための方法として、対数線形モデルによる解析がある。このページでは詳しく解説しないが、日本語で読める文献としては次のようなものが挙げられる。【1】はこちらのwebサイトから無料で観覧することができる。
【1】 松田紀之『質的情報の多変量解析』朝倉書店（1988）
【2】 渡邉ほか訳『カテゴリカルデータ解析入門』サイエンティスト社（2005）
【3】 田中ほか訳『一般化線形モデル入門 原著第2版』共立出版（2008）
さて分割表にまとめられたデータは従属変数を平方根変換することによって、一般線形モデルとして扱うことができるということを紹介した。また、それ以外の方法として対数線形モデルというものがあるということにも触れた。今度はその対数線形モデルを用いて解析してみよう。ただし、数学的な詳しい解説は控えることにするので、それについては上で挙げた参考文献を読むとよいだろう。
ここでは最も基本的な二元分割表の解析を紹介し、続いて三元分割表と多元分割表の解析例を紹介していくことにしよう。ただそのまえにいずれの分割表の解析においても、共通して知っておくべき3つのモデルがあるということについて解説しておく。
対数線形モデルを用いて分割表を解析する際には、常に1)飽和モデル、2)加法モデル、3)最小モデルの3つのモデルが考えられる。飽和モデルとは、モデルの中に考えられる全ての交互作用項を含めたモデルである。i行j列のセルの観測度数をμijとし、i行の行効果をAi、j列の列効果をBjとすると飽和モデルは次のように表される。
logμij = λ + Ai + Bj + AiBj
次に加法モデルについてだが、これは独立モデルともよばれる。なぜなら、このモデルが成り立つとき、古典的な独立性の検定における帰無仮説と同じく、2つの変数（要因）は互いに独立であるということができるからである。加法モデルは次式によって表現されるが、要するにこれは交互作用項を含めない場合のモデルである。
logμij = λ + Ai + Bj
最後に最小モデルについて説明する。これは文字通り最も小さなモデルという意味で、モデルには定数項λしかない状態である。最小モデルにおけるλは前セルの平均値となる。したがって、全てのセルにおける観測度数が全く同じであることを意味している。
logμij = λ
それでは二元分割表のデータを実際に解析してみよう。なお、二元分割表というのは表の形式のことで、AとBの2つの要因があるということである。
上に挙げた参考文献を読みながら学習を進められるように、文献【2】のp203-表6.1のデータをRで解析してみよう。この場合、性別（Sex）という要因と生死観（LorD）という要因があるわけである。
| belive | not belive | |
| male | 435 | 147 |
| female | 375 | 134 |
# 観測度数Countと性別Sexと生死観LorDを用意する
> Count <- c(435, 147, 375, 134)
> Sex <- c("male", "male", "female", "female")
# factor()を使ってカテゴリカル型であることを指定する
# 引数levelsは水準の順位を指定するものである
> Sex <- factor(Sex, levels=c("male", "female"))
> LorD <- c("belive", "not belive", "belive", "not belive")
# Sexの場合と同様に、カテゴリカル型であることを指定している
> LorD <- factor(LorD, levels=c("belive", "not belive"))
# 各変数の中身を確認してみる
> Count
[1] 435 147 375 134
> Sex
[1] male male female female
Levels: male female
> LorD
[1] belive not belive belive not belive
Levels: belive not belive
# glm()を用いてモデル解析
# SexとLorDを*でつないでいるのは、飽和モデルであることを表す
# 引数familyではポアソン分布に従うことを仮定している
> model <- glm(Count ~ Sex * LorD, family=poisson)
> summary(model) # 係数表の出力
Call:
glm(formula = Count ~ Sex * LorD, family = poisson)
Deviance Residuals:
[1] 0 0 0 0
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.07535 0.04795 126.711 <2e-16 ***
Sexfemale -0.14842 0.07047 -2.106 0.0352 *
LorDnot belive -1.08491 0.09540 -11.372 <2e-16 ***
Sexfemale:LorDnot belive 0.05583 0.13868 0.403 0.6873
---
係数表を分かりやすいように数式でまとめてみる。
よりこの形に近い結果を得るためには、dummy.coef()という関数を利用すればよい。
> dummy.coef(model)
Full coefficients are
(Intercept): 6.075346
Sex: male female
0.00000 -0.14842
LorD: belive not belive
0.000000 -1.084913
Sex:LorD: male:belive female:belive male:not belive
0.00000000 0.00000000 0.00000000
(Intercept):
Sex:
LorD:
Sex:LorD: female:not belive
0.05582722
さらに、得られた予測式から当てはめ値の対数を求めるためにはpredict()という関数を用いるとよい。また、関数exp()によって、当てはめ値を計算することができる。
> predict(model)
1 2 3 4
6.075346 4.990433 5.926926 4.897840
> exp(predict(model))
1 2 3 4
435 147 375 134
ご覧のとおり、飽和モデルで予測された当てはめ値は観測度数と完全に一致していることが分かる。ただし、このモデルが統計的に有意であるかどうかは別問題であることに注意してほしい。実際、上の係数表において、交互作用項に対するp値から明らかに有意ではないことがわかる（p = 0.6873）。
結論をいってしまうと、加法モデルが適切である。したがって、関数update()を使ってより適切なモデルを当てはめることを考える必要がある。
> model2 <- update(model, ~. -Sex:LorD)
> summary(model2)
Call:
glm(formula = Count ~ Sex + LorD, family = poisson)
Deviance Residuals:
1 2 3 4
0.1394 -0.2377 -0.1494 0.2524
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.06865 0.04512 134.488 <2e-16 ***
Sexfemale -0.13402 0.06069 -2.208 0.0272 *
LorDnot belive -1.05868 0.06923 -15.291 <2e-16 ***
---
すると今度はすべての項が有意であると認められるようになった。これでモデル解析は終了となる。先ほどと同じように、当てはめ値と当てはめ値の対数を求めてみると次のようになる。
> predict(model2) # 当てはめ値の対数
1 2 3 4
6.068655 5.009975 5.934632 4.875953
> exp(predict(model2)) # 当てはめ値
1 2 3 4
432.099 149.901 377.901 131.099
ここでは練習のために、さらにモデルからSexとLorDを除いて最小モデルを作成してみよう。
> model3 <- update(model2, ~. -Sex -LorD)
> summary(model3)
Call:
glm(formula = Count ~ 1, family = poisson)
Deviance Residuals:
1 2 3 4
9.034 -8.353 5.855 -9.329
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.60856 0.03027 185.3 <2e-16 ***
---
最小モデルから当てはめ値を求めてみる。
> predict(model3)
1 2 3 4
5.608556 5.608556 5.608556 5.608556
> exp(predict(model3))
1 2 3 4
272.75 272.75 272.75 272.75
先にも述べたように、最小モデルによって予測された当てはめ値（272.75）は前セルの平均値である。
> mean(Count) [1] 272.75
ところでこれらの対数線形モデルが受容されるという検定の統計量はどのようにして求められるのだろうか。それは次式によって求められる。
検定統計量X2は自由度(i-1)(j-1)のカイ自乗分布に従うという意味である。そして、この検定の帰無仮説は独立性の検定と同じで、「A要因とB要因は独立である（関連がない）」というものである。これに対する対立仮説は「A要因とB要因は独立でない（関連がある）」ということになる。
以下の結果をみると、飽和モデルがもっとも独立であるということについて有意であるといえる。ただし、前述したように、もともと誤ったモデルを当てはめて、それによって得られた当てはめ値から検定統計量を計算しても、それは無意味なことである。
今回は練習のために全てのモデルについて検定を行ったが、加法モデルが最も適したモデルであるわけであるから、加法モデルについて検定を行えばよいわけである。そして、その結果はp = 0.69であり、これは加法モデルが観測度数のデータによく当てはまっているといえるのである。
# 飽和モデル：検定統計量を求めてX1に代入 > X1 <- sum((Count - exp(predict(model))) ^ 2 / exp(predict(model))) > X1 [1] 2.893842e-25 > 1 - pchisq(X1, 1) # X1からp値を求める [1] 1 # 加法モデル > X2 <- sum((Count - exp(predict(model2))) ^ 2 / exp(predict(model2))) > X2 [1] 0.162084 > 1 - pchisq(X2, 1) [1] 0.6872451 # 最小モデル > X3 <- sum((Count - exp(predict(model3))) ^ 2 / exp(predict(model3))) > X3 [1] 263.4088 > 1 - pchisq(X3, 1) [1] 0
オッズ比
文献【2】にはオッズ比についての解説もされているが（p29〜）、これについてはメタアナリシスのメタアナリシスの基本用語のページで解説しているので、そちらを参照してもらいたい。ここではオッズ比が1であるとき、その分割表の2要因は互いに独立であることを示している、ということを覚えてもらいたい。
三元分割表とは要因が3つある場合の分割表のことである。文献【2】のp211-表6.3をRのglm()で解析するためには、次のようなデータフレームを用意する必要がある。
| 飲酒歴 | 喫煙歴 | マリファナ使用歴 | 観測度数 |
| 1 | 1 | 1 | 911 |
| 1 | 1 | 2 | 538 |
| 1 | 2 | 1 | 44 |
| 1 | 2 | 2 | 456 |
| 2 | 1 | 1 | 3 |
| 2 | 1 | 2 | 43 |
| 2 | 2 | 1 | 2 |
| 2 | 2 | 2 | 279 |
Rにデータを取り込む方法はいくつかあるが、ここでは1つ1つコマンドしていく方法を1回だけ紹介するので、あとの方法はユーザー自身で工夫してやってもらいたい。
# 飲酒
> A <- rep(c(1, 2), c(4, 4))
> A <- factor(A, levels=c(1, 2), labels=c("Yes", "No"))
# 喫煙
> C <- rep(c(1, 2, 1, 2), c(2, 2, 2, 2))
> C <- factor(C, levels=c(1, 2), labels=c("Yes", "No"))
# マリファナ
> M <- rep(1:2, 4)
> M <- factor(M, levels=c(1, 2), labels=c("Yes", "No"))
# 観測度数
> Count <- c(911, 538, 44, 456, 3, 43, 2, 279)
# データフレームの作成
> Table6.3 <- data.frame(
+ alcohol = A, cigarette = C, marijuana = M, COUNT = Count)
# glm()を使って、飽和モデルから解析を開始する
# データフレームから分析する場合は、引数dataに作成したデータフレームをしていすること
> model <- glm(COUNT ~ alcohol * cigarette * marijuana, family=poisson, data=Table6.3)
> summary(model) # 係数表の出力
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.81454 0.03313 205.682 < 2e-16
alcoholNo -5.71593 0.57830 -9.884 < 2e-16
cigaretteNo -3.03035 0.15435 -19.633 < 2e-16
marijuanaNo -0.52668 0.05437 -9.686 < 2e-16
alcoholNo:cigaretteNo 2.62489 0.92583 2.835 0.00458
alcoholNo:marijuanaNo 3.18927 0.59962 5.319 1.04e-07
cigaretteNo:marijuanaNo 2.86499 0.16696 17.159 < 2e-16
alcoholNo:cigaretteNo:marijuanaNo -0.58951 0.94236 -0.626 0.53160
さて、上の係数表からalcohol:cigarette:marijuanaの3次の交互作用項は明らかに有意でないことが確認できる。そのため、update()を使ってこの3次の交互作用項を除いたモデルを解析してみる。
> model2 <- update(model, ~. -alcohol:cigarette:marijuana)
> summary(model2)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.81387 0.03313 205.699 < 2e-16 ***
alcoholNo -5.52827 0.45221 -12.225 < 2e-16 ***
cigaretteNo -3.01575 0.15162 -19.891 < 2e-16 ***
marijuanaNo -0.52486 0.05428 -9.669 < 2e-16 ***
alcoholNo:cigaretteNo 2.05453 0.17406 11.803 < 2e-16 ***
alcoholNo:marijuanaNo 2.98601 0.46468 6.426 1.31e-10 ***
cigaretteNo:marijuanaNo 2.84789 0.16384 17.382 < 2e-16 ***
すると、今度は全ての項が有意であるので、これ以上モデルを節約する必要はない。結論としてalcohol:cigarette（AC）とalocohol:marijuana（AM）とcigarette:marijuana（CM）といった3つの2次交互作用項を含めたモデルが最適であるといえる。これは教科書でいうところのモデル（AC, AM, CM）である。
ここで当てはめ値を計算してみよう。
> atehame <- exp(predict(model2))
> atehame
1 2 3 4 5 6
910.383170 538.616830 44.616830 455.383170 3.616830 42.383170
7 8
1.383170 279.616830
さらに二元分割表の解析でも用いた、X2統計量を次式によって計算してみよう。復習になるが、これはモデルの当てはめが良いかどうかを検定するためのものである。
# X^2統計量の計算 > sum((Count - atehame)^2 / atehame) [1] 0.4011005 # X^2に対するp値の計算 > 1 - pchisq(0.4011005, 1) [1] 0.5265215
この場合の帰無仮説は
H0 ： 三元分割表の期待度数はある対数モデルを満たす
というものである。ある対数モデルというのは、今回の例でいえばモデル（AC, AM, CM）のことである。また、この場合、帰無仮説を採択できることが望ましいので、計算されたp値は1に近いほど（より大きいほど）よい。
飽和モデルについては、そもそも自由度が0なので検定を行うことさえできない。ちなみに、自由度はsummary()によって得られる結果の一番下のほうに出力されている
# これはmodelのヤツ
Residual deviance: 2.6690e-13 on 0 degrees of freedom
# これはmodel2のヤツ
Residual deviance: 0.37399 on 1 degrees of freedom
という部分を見ればよい。正確な用語でいうと、残差逸脱度と残差自由度という。model2の残差逸脱度は0.37399であり、残差自由度は1である。また、残差自由度よりも残差逸脱度がかなり大きな場合は過分散を疑ったほうがよい（過分散については文献【2】p125〜127を参照のこと）。
少し寄り道的な話をしよう。もし自由度がいくつか分からないときは、解析したモデルにいくつのパラメータを推定したのかを考えればよい。今回の例だと、A、C、Mいずれも2つの水準（Yes / No）をもっていた。したがって、全自由度は2 * 2 * 2 = 8である。そして推定したパラメータは（A, C, M, AC, AM, CM）という6つの項と定数項（Intercept）の計7つである。自由度は全自由度から推定したパラメータ数を引いたものになる。すなわち、8 - 7 = 1となるわけである。
この自由度のルール（決まりごと）は全ての解析において共通していえる基本なので、これを覚えておけば、場合場合によって異なる自由度の計算公式を暗記しておかなくても問題ないだろう。
ところで当てはまりの良さを判断するためには標準化残差を見てみるという方法もある。標準化残差について説明する前に、ピアソン残差とピアソンの適合度統計量X2の関係について解説しておこう。ピアソン残差eiは、
によって計算される。そしてこのピアソン残差とピアソン適合度検定量には
という関係が成り立つ。
実際に計算してみると同じであることがわかる。
# X^2統計量の計算 > sum((Count - atehame)^2 / atehame) [1] 0.4011005 # ピアソン残差の2乗和 > a <- (Count - atehame) / sqrt(atehame) > sum(a^2) [1] 0.4011005
さて、標準化残差の定義式をみてみよう。標準化残差のことを、別の言い方でスチューデント化残差などともいう。次式のhiはてこ比とよばれるものであるが、てこ比の計算公式はここには載せることはしない。
自分の手で計算しなくても、Rには標準化残差を計算してくれる関数が用意されている。それはrstudent()という関数である。
> rstudent(model2)
1 2 3 4 5 6
0.6333248 -0.6333251 -0.6333561 0.6333246 -0.6384745 0.6332908
7 8
0.6061680 -0.6333257
この標準化された残差の絶対値が2を越える場合は当てはまりが悪いといえる。ちなみに、この標準化残差のプロットは
> plot(model2, which=3)
とコマンドすることによって出力される。
4因子以上の分散分析のことを多因子分散分析（あるいは多要因分散分析）といったが、分割表の解析の場合も4元以上の分割表のことを多元分割表という。一般にある1つの実験において4つ以上の要因を組み込むのは、解析が非常に複雑になるためできるだけ避けるようにするものである。せいぜい3要因までが限度であると認識しておいた方がよいかもしれない。
特に分割表の解析においては、多元分割表の場合、モデルの構築が非常に複雑で困難であることを肝に銘じるべきである。もっとも困難であるかどうかは、取り扱う問題や得られたデータセットにもよるし、同じデータを解析する場合であっても分析者の基本的な技術や経験によるものであることに変わりはない。
多元分割表の解析方法とその結果の解釈の仕方の基本は今まで取り扱ってきた内容と同じなので、ここでは分析例のみを提示することにしよう。用いるデータセットは文献【2】のp220-表6.8である。下表はそのデータフレームである。
| G | L | S | I | Count |
| 1 | 1 | 1 | 1 | 7287 |
| 1 | 1 | 1 | 2 | 996 |
| 1 | 1 | 2 | 1 | 11587 |
| 1 | 1 | 2 | 2 | 759 |
| 1 | 2 | 1 | 1 | 3246 |
| 1 | 2 | 1 | 2 | 973 |
| 1 | 2 | 2 | 1 | 6134 |
| 1 | 2 | 2 | 2 | 757 |
| 2 | 1 | 1 | 1 | 10381 |
| 2 | 1 | 1 | 2 | 812 |
| 2 | 1 | 2 | 1 | 10969 |
| 2 | 1 | 2 | 2 | 380 |
| 2 | 2 | 1 | 1 | 6123 |
| 2 | 2 | 1 | 2 | 1084 |
| 2 | 2 | 2 | 1 | 6693 |
| 2 | 2 | 2 | 2 | 513 |
Rによる解析例 ：
# 飽和モデルから解析を開始する
> model <- glm(Count ~ G * L * S * I, family=poisson)
> summary(model)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.89385 0.01171 759.214 < 2e-16 ***
Gmale 0.35389 0.01528 23.156 < 2e-16 ***
Lsuburb -0.80867 0.02110 -38.322 < 2e-16 ***
Syes 0.46379 0.01495 31.021 < 2e-16 ***
Iyes -1.99010 0.03378 -58.909 < 2e-16 ***
Gmale:Lsuburb 0.28074 0.02655 10.574 < 2e-16 ***
Gmale:Syes -0.40870 0.02027 -20.159 < 2e-16 ***
Lsuburb:Syes 0.17263 0.02636 6.550 5.75e-11 ***
Gmale:Iyes -0.55813 0.04969 -11.232 < 2e-16 ***
Lsuburb:Iyes 0.78531 0.04977 15.779 < 2e-16 ***
Syes:Iyes -0.73554 0.05045 -14.580 < 2e-16 ***
Gmale:Lsuburb:Syes -0.13872 0.03457 -4.013 5.99e-05 ***
Gmale:Lsuburb:Iyes 0.03153 0.06993 0.451 0.6521
Gmale:Syes:Iyes -0.07889 0.08121 -0.971 0.3314
Lsuburb:Syes:Iyes -0.15191 0.07325 -2.074 0.0381 *
Gmale:Lsuburb:Syes:Iyes 0.12918 0.11225 1.151 0.2498
---
# G:L:S:Iの項が有意ではないので取り除く
> model2 <- update(model, ~. -G:L:S:I)
> summary(model2)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.89244 0.01166 762.737 < 2e-16 ***
Gmale 0.35628 0.01514 23.527 < 2e-16 ***
Lsuburb -0.80410 0.02071 -38.831 < 2e-16 ***
Syes 0.46609 0.01482 31.447 < 2e-16 ***
Iyes -1.97844 0.03209 -61.662 < 2e-16 ***
Gmale:Lsuburb 0.27351 0.02579 10.607 < 2e-16 ***
Gmale:Syes -0.41291 0.01994 -20.705 < 2e-16 ***
Lsuburb:Syes 0.16551 0.02561 6.463 1.03e-10 ***
Gmale:Iyes -0.58350 0.04459 -13.087 < 2e-16 ***
Lsuburb:Iyes 0.75989 0.04460 17.038 < 2e-16 ***
Syes:Iyes -0.76170 0.04511 -16.886 < 2e-16 ***
Gmale:Lsuburb:Syes -0.12646 0.03288 -3.846 0.000120 ***
Gmale:Lsuburb:Iyes 0.08176 0.05469 1.495 0.134914
Gmale:Syes:Iyes -0.01144 0.05603 -0.204 0.838208
Lsuburb:Syes:Iyes -0.09685 0.05547 -1.746 0.080803 .
---
# G:S:Iも有意ではない(p = 0.838)ので取り除く
> model3 <- update(model2, ~. -G:S:I)
> summary(model3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.89220 0.01160 766.542 < 2e-16 ***
Gmale 0.35669 0.01501 23.762 < 2e-16 ***
Lsuburb -0.80445 0.02064 -38.968 < 2e-16 ***
Syes 0.46648 0.01470 31.738 < 2e-16 ***
Iyes -1.97646 0.03057 -64.654 < 2e-16 ***
Gmale:Lsuburb 0.27400 0.02568 10.671 < 2e-16 ***
Gmale:Syes -0.41363 0.01963 -21.072 < 2e-16 ***
Lsuburb:Syes 0.16601 0.02549 6.512 7.4e-11 ***
Gmale:Iyes -0.58782 0.03923 -14.984 < 2e-16 ***
Lsuburb:Iyes 0.76046 0.04448 17.096 < 2e-16 ***
Syes:Iyes -0.76616 0.03946 -19.416 < 2e-16 ***
Gmale:Lsuburb:Syes -0.12723 0.03266 -3.895 9.8e-05 ***
Gmale:Lsuburb:Iyes 0.08148 0.05466 1.491 0.1360
Lsuburb:Syes:Iyes -0.09773 0.05529 -1.768 0.0771 .
---
# G:L:Iも有意ではない(p = 0.136)ので取り除く
> model4 <- update(model3, ~. -G:L:I)
> summary(model4)
Call:
glm(formula = Count ~ G + L + S + I + G:L + G:S + L:S + G:I +
L:I + S:I + G:L:S + L:S:I, family = poisson)
Deviance Residuals:
1 2 3 4 5 6 7
-0.05787 0.15683 -0.20767 0.82005 0.14031 -0.25537 0.24697
8 9 10 11 12 13 14
-0.69640 0.04850 -0.17305 0.21372 -1.12571 -0.10203 0.24320
15 16
-0.23596 0.86406
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.89453 0.01148 774.713 < 2e-16 ***
Gmale 0.35273 0.01477 23.886 < 2e-16 ***
Lsuburb -0.81181 0.02007 -40.449 < 2e-16 ***
Syes 0.46504 0.01465 31.740 < 2e-16 ***
Iyes -1.99575 0.02783 -71.720 < 2e-16 ***
Gmale:Lsuburb 0.28566 0.02447 11.674 < 2e-16 ***
Gmale:Syes -0.41151 0.01957 -21.028 < 2e-16 ***
Lsuburb:Syes 0.17069 0.02532 6.741 1.57e-11 ***
Gmale:Iyes -0.54594 0.02729 -20.007 < 2e-16 ***
Lsuburb:Iyes 0.80160 0.03490 22.966 < 2e-16 ***
Syes:Iyes -0.76173 0.03933 -19.366 < 2e-16 ***
Gmale:Lsuburb:Syes -0.13363 0.03239 -4.126 3.70e-05 ***
Lsuburb:Syes:Iyes -0.10800 0.05486 -1.968 0.049 *
---
# 全ての項が有意になったが、念のためL:S:Iを除いて解析してみる
> model5 <- update(model4, ~. -L:S:I)
> summary(model5)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.89195 0.01142 778.756 < 2e-16 ***
Gmale 0.35365 0.01477 23.944 < 2e-16 ***
Lsuburb -0.80411 0.01966 -40.891 < 2e-16 ***
Syes 0.46962 0.01447 32.450 < 2e-16 ***
Iyes -1.97446 0.02547 -77.507 < 2e-16 ***
Gmale:Lsuburb 0.28274 0.02441 11.584 < 2e-16 ***
Gmale:Syes -0.41328 0.01956 -21.134 < 2e-16 ***
Lsuburb:Syes 0.15752 0.02441 6.453 1.09e-10 ***
Gmale:Iyes -0.54483 0.02727 -19.982 < 2e-16 ***
Lsuburb:Iyes 0.75806 0.02697 28.105 < 2e-16 ***
Syes:Iyes -0.81710 0.02765 -29.551 < 2e-16 ***
Gmale:Lsuburb:Syes -0.12858 0.03228 -3.984 6.78e-05 ***
---
さて、model4は次のようなモデルである。
Count ~ G + L + S + I + G:L + G:S + L:S + G:I + L:I + S:I + G:L:S + L:S:I
一方でmodel5は次のようなモデルである。
Count ~ G + L + S + I + G:L + G:S + L:S + G:I + L:I + S:I + G:L:S
両者の違いは3次交互作用項であるL:S:Iを含むか含まないかという点である。どちらがより良いモデルであるかを比較してみよう。まずAICによるモデル比較をしてみよう。AICの値は小さいほど当てはまりの良いモデルであることを示している。RではAIC()という関数が用意されている。
# model4のAIC > AIC(model4) [1] 183.0509 # model5のAIC > AIC(model5) [1] 184.9239
どうやらAICを比較した上ではmodel4のほうが当てはまりの良いモデルだといえそうである。では標準化残差を見比べてみた場合ではどうだろうか。この標準化残差の絶対値に2以上のものがあると当てはまりが良くないことを示すのであった。この際、ベクトルの最小値と最大値を返してくれるrange()を使うと都合が良い。
# model4の標準化残差のレンジ > range(rstudent(model4)) [1] -1.578777 1.571241 # model5の標準化残差のレンジ > range(rstudent(model5)) [1] -2.483799 2.471588
これをみるとmodel5には2を越える標準化残差の値が存在することが分かる。したがって、model5よりもmodel4のほうが当てはまりの良いモデルであるようである。だが、ここでさらにX2検定統計量を計算してみよう。
# 当てはめ値の計算 > m4.ate <- exp(predict(model4)) > m5.ate <- exp(predict(model5)) # X^2統計量の計算 > m4.X <- sum((Count - m4.ate)^2 / m4.ate) > m5.X <- sum((Count - m5.ate)^2 / m5.ate) # p値の計算 > 1 - pchisq(m4.X, 3) # model4の場合 [1] 0.3105178 > 1 - pchisq(m5.X, 4) # model5の場合 [1] 0.1122669
結局、model4のp値（p = 0.311）よりもmodel5のp値（p = 0.112）の方が小さく、model4の方が当てはまりが良さそうであることが分かる。
モデル選択における注意点
ここで紹介したモデル選択の方法はp値をもとに有意でない項を取り除いていくというものであった。たしかにこれは基本的な方法ではあるが、必ずしも適切な方法であるとはいえない。なぜならば、統計的に有意な結果が必ずしも実用面で有用であるとは限らないからである。というのも、たとえ統計学的に最適であると判断されるようなモデルであっても、その結果が複雑すぎてどう解釈してよいか分からないということが起こり得るからである。仮に分析者自信が納得できたとしても、それが他の多くの研究者に納得してもらえなければ意味がない。
したがって、モデル選択を行う際には常に統計的有意性と実用的有意性について考える必要がある。要するに、適合度検定の結果（p値）のみにすべての判断を委ねるというようなことは不適切である、ということである。
試しにmodel4の結果を自分でどう解釈できるかを考えてみたらよいだろう（教科書の説明を一緒に読みながら）。そうすると、おそらくここで述べたことの意味が分かるはずである。
対数線形モデルの解説の最後にちょっとした補足を加えておこう。文献【2】p203の表6.1には、パラメータを推定する際にどの変数（正確にはダミー変数）をベースラインとして解析するかによって、推定される数値が異なるという結果がまとめられている（そうはいっても、これらの推定値は実質的に全く同じ結果である）。これについてRのコマンド方法を解説しておく。
Rのglm()では何も指定しなければSet2のような結果が返される。つまり、最初の変数をベースラインとして計算しているわけである。
> Set2 <- glm(Count ~ Sex + LorD, family=poisson)
> dummy.coef(Set2)
Full coefficients are
(Intercept): 6.068655
Sex: male female
0.0000000 -0.1340224
LorD: belive not belive
0.000000 -1.058680
続いてSet1についてだが、これは最後の変数をベースラインとして計算している。これはSASで採用されている方法である。RではSASの真似もできる。そして、それにはglm()の引数contrastsで次のように指定する必要がある。
> Set1 <- glm(Count ~ Sex + LorD, family=poisson, contrasts=list(Sex="contr.SAS", LorD="contr.SAS"))
> dummy.coef(Set1)
Full coefficients are
(Intercept): 4.875953
Sex: male female
0.1340224 0.0000000
LorD: belive not belive
1.058680 0.000000
そしてSet3は各変数の合計をベースラインとして解析する方法である。これはMinitabで採用されている計算方法で、この場合、定数項は全セルの平均値と一致する。
> Set3 <- glm(Count ~ Sex + LorD, family=poisson, contrasts=list(Sex="contr.sum", LorD="contr.sum"))
> dummy.coef(Set3)
Full coefficients are
(Intercept): 5.472304
Sex: male female
0.06701122 -0.06701122
LorD: belive not belive
0.5293398 -0.5293398
# 当てはめ値の対数の平均
> (6.069 + 5.010 + 5.935 + 4.876) / 4
[1] 5.4725
【1】 松田紀之『質的情報の多変量解析』朝倉書店（1988）
【2】 渡邉ほか訳『カテゴリカルデータ解析入門』サイエンティスト社（2005）
【3】 田中ほか訳『一般化線形モデル入門 原著第2版』共立出版（2008）
【4】 野間口・菊池訳『統計学：Rを用いた入門書』共立出版（2008）
【5】 野間口・野間口訳『一般線形モデルによる生物科学のための現代統計学-あなたの実験をどのように解析するか-』共立出版（2007）
【6】 間瀬ほか著『工学のためのデータサイエンス入門 フリーな統計環境Rを用いたデータ解析』数理工学社（2004）
