二値および計数データの回帰分析
Hiroyoshi Iwata
2025/01/10
必要なパッケージ
この授業では、農学関係のデータを集めたagridatを使う。また、data.frameを編集するのに便利なツールがそろうtidyverseも使う。
require(agridat)
require(tidyverse)ロジスティック回帰
共変量\(x_i\)に依存する応答\(y_i\)がカテゴリ変数で、2つの異なる結果をもつ場合を考える。そのような応答の例として、罹病と非罹病、成功と失敗、勝ちと負けなどが考えられる。こうした2値の応答では、簡単のために、応答の種類によらず、応答の値として0、1という仮の数値が割りふられる。ここで、興味があることは、観察された共変量のもとで、\(y=1\)が得られる確率をモデル化することである。
このようなモデル化に、単回帰モデル\(y_i = \beta_0 + \beta_1 x_i + \epsilon_i\)を用いることは、適切ではない。応答\(y_i\)は0, 1の2つの値しかとることができないが、共変量\(x_i\)の値には制約がないために、\(\beta_0 + \beta_1 x_i\)の値にも制約がなく、単回帰モデルは0より小さい、あるいは、1より大きな値を与える可能性がある。そこで、変換された尺度のもとで確率モデルを組み立てることを考える。
具体的には、\(y\)が\(Ey = P(y=1)=p\)のベルヌーイ分布に従うとする。すると、\(E(y_i | X = x_i) = P(y_i = 1 | X = x_i) = p_i\)は0から1の値をとるため、\(p_i\)を累積分布関数\(F\)により、\(F(\beta_0 + \beta_1 x_i)\)とモデル化すれば、\(\beta_0 + \beta_1 x_i\)がどのような値をとっても、\(p_i\)すなわち\(P(y_i = 1 | X = x_i)\)は0から1の値をとる。このとき、\(p_i\)と\(\beta_0 + \beta_1 x_i\)の関係は、累積分布関数の逆関数\(F^{-1}\)を用いて、 \[ F^{-1}(p_i) = \beta_0 + \beta_1 x_i \] と表せる。
累積分布関数\(F\)は単調増加するので、原理的には、どのような\(F\)でも、確率\(p\)と共変量\(y\)を結びつける関数（リンク（link）関数とよばれる）として利用できる。しかし、多くの場合は、ロジスティック分布、正規分布、相補的log-log（complementary log-log）分布の累積分布関数がリンク関数として用いられ、それぞれ、ロジスティック回帰、プロビット回帰、clog-log回帰とよばれるモデルとなる。このうち最もよく用いられているのがロジスティック回帰である。ロジスティック回帰では、\(F^{-1}(p) = \log \frac{p}{1-p}\)をロジットとよび、\(logit(p)\)と記す。
ロジスティック回帰がよく用いられるのは、回帰係数の解釈が容易なためである。回帰係数\(\beta_1\)は、2値の応答である\(y\)に対する共変量の効果を表すが、ロジスティック回帰では、その値を\(\{y=1\}\)という事象に関する対数オッズとして解釈できる。具体例をもとにこれを見ていく。
ロジスティック回帰の当てはめ
ロジスティック回帰における基本的な統計モデルは、応答\(y_i\)（値は、0または1）がパラメータ\(p_i\)のベルヌーイ分布\(Ber(p_i)\)として、 \[ y_i \sim Ber(p_i) \\ \text {logit}(p_i) = \log \frac{p_i}{1-p_i}=\beta_0+\beta_1x_i, \:\:i=1,...,n \] と表される。
同じ共変量の値\(x_i\)に対して、複数の計測データがある場合は、応答をベルヌーイ分布ではなく、二項分布\(Bin(n_i, p_i)\)からの標本として表現するほうが便利である。すなわち、\(x_i\)に対する応答を\(n_i\)標本分まとめて、 \[ y_i \sim Bin(n_i, p_i) \\ \text {logit}(p_i) = \log \frac{p_i}{1-p_i}=\beta_0+\beta_1x_i, \:\:i=1,...,n \\ \sum_{i=1}^k n_i = n \] と表せる。原理的には、二項分布応答のモデルを、\(y_i\)個分、応答1の標本をならべ、\(n_i - y_i\)個分、応答0の標本を並べることで、ベルヌーイ分布のモデルとして表せる。
上の2つのロジスティックモデルのパラメータ\(\beta_0, \beta_1\)は、線形回帰に用いられる最小二乗法では求められない。モデル係数の推定には、非線形方程式を解く必要があり、数値計算を繰り返して解が得られる。例えば、Newton-Raphson法などが、非線形尤度方程式の解法に用いられ、母数\(\beta_0, \beta_1\)の推定量\(b_0, b_1\)を計算できる。モデルパラメータの推定方法の詳細は、ここでは説明しない。
一度、パラメータ\(\beta_0, \beta_1\)が推定されれば、確率\(p_i = P(y = 1| x = x_i)\)は、logit関数の逆関数（ロジスティック関数）を用いて以下のように計算できる。 \[ \hat p_i = \frac {\exp(b_0 + b_1 x_i)}{1 + \exp(b_0 + b_1 x_i)}=\frac 1 {1+\exp(-b_0 - b_1x_i)} \] これは、\(z_i = \log \frac{p_i}{1-p_i}\)として、その両辺の指数をとり、整理すると \[ \exp(z_i) = \frac{p_i}{1-p_i} \\ \exp(z_i) - \exp(z_i) p = p \\ p = \frac{\exp(z_i)}{\exp(z_i) + 1} = \frac{1}{1 + exp(-z_i)} \] として導出できる。
ここまでは、共変量が1つの場合について説明してきたが、共変量が\(p-1\)個の場合は、 \[ \boldsymbol{x_i'b}=l_i=b_0+b_1x_{i1}+b_2x_{i2}+\cdots+b_{p-1}x_{i,p-1} \] で\(b_0+b_1x_i\)を置き換える。ここで、\(\boldsymbol{x_i'}\)は、サイズ\(n\times p\)の行列\(\boldsymbol{X}\)の\(i\)番目の行であり、\(\boldsymbol{b}\)は回帰係数の推定値\(\boldsymbol{b}=(b_0,b_1,...,b_{p-1})'\)のベクトルである。
このとき、ロジスティック回帰モデルは、 \[ \hat p_i = \frac {\exp(\boldsymbol{x_i'b})}{1 + \exp(\boldsymbol{x_i'b})}=\frac 1 {1+\exp(-\boldsymbol{x_i'b})} \] と表される。
例 帝王切開における感染症
帝王切開を受ける母親は、感染症、過度の出血、血栓、分娩後の痛みの増加、入院期間の延長、および、著しく長い回復期間などのリスクがある。ここでの例は、帝王切開による感染症に関するものである（Fahrmeir and Tutz、1996、Vidakovic 2011より引用）。関心のある応答変数は、感染の発生または非発生であり、それに対して、それぞれ2つの水準をもつ3つの共変量による影響をモデル化する。
- noplan - 帝王切開分娩が計画的なものか（0）、計画されていなかったものか（1）
- riskfac - 糖尿病、肥満、帝王切開の経験など、母親に危険因子が存在するか（1）、存在しないか（0）
- antibio - 予防として抗生物質が与えられたか（1）、与えられなかったか（0）
これらの共変量のもとで整理された結果は以下のとおり。
| Planned | No Plan | ||||
|---|---|---|---|---|---|
| Infection | Infection | ||||
| yes | no | yes | no | ||
| Antibiotics | |||||
| Risk factor yes | 1 | 17 | 11 | 87 | |
| Risk factor no | 0 | 2 | 0 | 0 | |
| No Antibiotics | |||||
| Risk factor yes | 28 | 30 | 23 | 3 | |
| Risk factor no | 8 | 32 | 0 | 9 |
infection <- c(1, 11, 0, 0, 28, 23, 8, 0)
noinfection <- c(17, 87, 2, 0, 30, 3, 32, 9)
total <- c(18, 98, 2, 0, 58, 26, 40, 9)
noplan <- c(0, 1, 0, 1, 0, 1, 0, 1)
riskfac <- c(1, 1, 0, 0, 1, 1, 0, 0)
antibio <- c(1, 1, 1, 1, 0, 0, 0, 0)
model <- glm(cbind(infection, noinfection) ~ noplan + riskfac + antibio, family = binomial)
summary(model)##
## Call:
## glm(formula = cbind(infection, noinfection) ~ noplan + riskfac +
## antibio, family = binomial)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8926 0.4124 -4.589 4.45e-06 ***
## noplan 1.0720 0.4254 2.520 0.0117 *
## riskfac 2.0299 0.4553 4.459 8.25e-06 ***
## antibio -3.2544 0.4813 -6.761 1.37e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 83.491 on 6 degrees of freedom
## Residual deviance: 10.997 on 3 degrees of freedom
## AIC: 36.178
##
## Number of Fisher Scoring iterations: 5推定に用いられたモデルは、 \[ \log \frac{P(\mbox{infection})}{P(\mbox{no infection})} = \beta_0 + \beta_1 \cdot \mbox{noplan} + \beta_2 \cdot \mbox{riskfac} + \beta_3 \cdot \mbox{antibio} \] となる。glm関数で計算した結果、\(\beta\)の推定値は、\(b_0 = -1.8926, b_1 = 1.0720, b_2 = 2.0299, \beta_3 = -3.2544\)と求められる。
\(\beta\)の推定値は、先述したようにオッズ比\(\frac{P(\mbox{infection})}{P(\mbox{no infection})}\) に注目すると解釈しやすくなる。上のモデル式の両辺の対数をとると、オッズ比は、 \[ \frac{P(\mbox{infection})}{P(\mbox{no infection})} = \exp(\beta_0) \cdot \exp(\beta_1\cdot \mbox{noplan} )\cdot \exp(\beta_2 \cdot \mbox{riskfac} )\cdot \exp(\beta_3 \cdot \mbox{antibio}) \] と表される。これを求められた\(\beta\)の推定値を当てはめて表すと、 \[ \frac{P(\mbox{infection})}{P(\mbox{no infection})} = \exp(-1.8926)\cdot \exp(1.0720\cdot \mbox{noplan} )\cdot \exp(2.0299 \cdot \mbox{riskfac} )\cdot \exp(-3.2544 \cdot \mbox{antibio}) \] となる。
例えば、抗生物質が与えられたとき（antibio=1）は、与えられなかったとき（antibio=0)に比べ、感染のオッズの推定値\(\frac{P(\mbox{infection})}{P(\mbox{no infection})}\)が、\(\exp(-3.25) = 0.0388\)倍になる。ただし、こうした推論は、モデルがよく当てはまっている場合にのみ有効であることに注意する。また、他のモデル（例えば、プロビットやclog-log）を用いた場合には、推定されるリスク因子の影響は異なってくる可能性がある。
なお、先述した結果では、\(\beta\)の標準誤差と、\(z\)値が示されている。これは、Waldの\(Z\)統計量とよばれるもので、正規分布で近似され、\(\beta\)の仮説検定に用いられます。この点については、後ほど説明する。
逸脱度（deviance）は、通常の回帰分析における\(R^2\)などと同様に、適合度の尺度で、この例では自由度3の\(\chi^2\)分布に従う。自由度は、全グループ数の7（1組合せは欠測のため）から推定されたパラメータ数の4（\(p\)）を減じて、3と計算される。
それでは、逸脱度を検定してみよう。
dev <- 10.9967
pval <- 1 - pchisq(dev, 7-4)
pval## [1] 0.01174373検定の結果得られた\(p値\)の値から、逸脱度が有意であることが分かる。つまり、このモデルの当てはまりは良くない。
そもそもなぜ回帰モデルが必要なのか、ということを疑問に感じるかもしれない。興味ある確率は、相対的な頻度としても推定できるためである。例えば、(noplan = 0, riskfac = 1, antibio = 1)のとき、感染の相対的頻度は\(1/18 = 0.0556\)となり、モデルから推定される\(\hat p = 0.0424\)と大きく異ならない。しかし、モデル化することには2つの利点がある。1つ目の利点は、モデルは観察されなかった場合、例えば（noplan = 1, riskfac = 0, antibio = 1）、についても感染の確率を予測できる点である。2つ目の利点は、感染（\(y = 1\)）が見られなかった場合についても、モデルを使うことで他のデータの「力を借りて」推定できる点である。例えば、(noplan = 1, riskfac = 0, antibio = 0)のとき、感染の確率を0とするよりも、モデルベースの推定値\(\hat p = 0.3056\)を提示するほうがより現実的である。
なお、次に示す図は、観察値とモデルから推定される感染割合を示す。共変量の三つ組み（noplan, riskfac, antibio）に対して、x軸のコードは、1 = (0,1,1), 2=(1,1,1), 3=(0,0,1), 4=(1,0,1), 5=(0,1,0), 6=(1,1,0), 7=(0,0,0), 8=(1,0,0)とする。なお、4=(1,0,1)は観察値が得られていない組合せだが、モデルベースであれば感染確率を求めることができる。
obs.p <- infection / total
model.p <- fitted(model)
matplot(cbind(obs.p, model.p), xlab = "Code for a triple", ylab = "Probability of Infection", pch = 1, ylim = c(0, 1))
legend("topleft", legend = c("Obs", "Model"), pch = 1, col = 1:2)
当てはめられたロジスティック回帰モデルの評価
線形回帰の当てはまりの良さを評価する尺度として、通常の回帰分析では、\(R^2, F, MSE\)などが用いられまたが、これらの尺度はロジスティック回帰においては適当ではない。ロジステック回帰で用いられるうちの、いくつかの尺度について説明する。
モデルパラメータ\(\beta_0, \beta_1, ...\)の有意性は、Wald検定とよばれる方法によって検定される。統計量\(Z_i = \frac {b_i}{s(b_i)}\)は、係数\(\beta_i\)が0のとき、正規分布で近似される。また、統計量\(W_i = \frac {b_i^2}{s^2(b_i)}\)は、自由度1の\(\chi^2\)分布に従う。\(|Z_i|\)や\(W_i\)が大きな値を示すと\(H_0: \beta_i = 0\)が棄却される。ここで、\(s(b_i), s^2(b_i)\)は、それぞれ、\(b_i\)の標本標準偏差と標本分散である。
適合度の尺度として慣習的に用いられているのは、逸脱度（deviance）で、 \[ D = -2 \log \frac {\mbox{likelihood of the fitted model}}{\mbox{likelihood of the saturated model}} \] と定義される。ここで、上述したロジスティック回帰では、\(y_i\)はクラス\(i\)の1の数で、\(n_i - y_i\)は0の数であり、尤度（likelihood）は、二項分布モデルのもとで\(y_i\)および\(n_i - y_i\)を観察する確率（に比例する値）であり、 \[ L = \prod_{i=1}^k \left(\begin{array}{cc} n_i\\y_i\end{array}\right) p_i^{y_i} (1 - p_i)^{n_i - y_i} \] と計算される。ここで、尤度の対数（対数尤度、log-likelihood）をとると、 \[ LL = \sum_{i=1}^k \left[y_i \log p_i + (n_i - y_i) \log(1-p_i) + \log \left(\begin{array}{cc} n_i\\y_i\end{array}\right) \right] \] となる。 このとき、逸脱度は、 \[ D = -2 \sum_{i = 1}^k \biggl[ y_i \log \biggl(\frac{\hat y_i}{y_i} \biggr) + (n_i - y_i) \log \biggl(\frac{n_i - \hat y_i}{n_i - y_i} \biggr) \biggr] \\ = 2 \sum_{i = 1}^k \biggl[ y_i \log \biggl(\frac{y_i} {\hat y_i}\biggr) + (n_i - y_i) \log \biggl(\frac{n_i - y_i}{n_i - \hat y_i} \biggr) \biggr] \] として計算される。なお、\(\hat y_i = n_i \hat p_i\)で\(y_i\)に対してモデルで推定された値（当てはめた値）である。飽和モデル（saturated model）とは、\(p_i\)を\(\hat p_i = y_i / n_i\)として推定し、\(\hat y_i = y_i\)として、完全に観察値に合わせたときのモデルである。逸脱度統計量\(D\)は、先にも述べたのように自由度\(k - p\)の\(\chi^2\)分布にしたがう。ここで、\(k\)はクラス/グループの数、\(p\)はモデルパラメータの数である。 当てはまりの良さ（適合度）の尺度\(G\)は、帰無モデル（切片のみのモデル）と現在検討中のモデルの逸脱度の違いとして定義される。\(G\)は自由度\(p-1\)の\(\chi^2\)分布に従い、\(G\)が小さな値をとったときに棄却される。\(G\)が小さな値をとることは、共変量を加えても逸脱度が改良されなかったことを示唆する。
個々のサンプルについての評価を行うためには、個別に残差を計算する必要がある。例えば、残差逸脱度という尺度があり、二項分布モデルでは \[ r_i^D = \mbox{sign}(y_i - \hat y_i) \sqrt {2 \left[y_i \log \left(\frac{y_i}{\hat y_i}\right) + (n_i - y_i)\log \left(\frac{n_i - y_i}{n_i - \hat y_i}\right)\right]}, \: i=1,...,k \] と表される。逸脱度は、分散分析のように、残差逸脱度の平方和として分割できる。すなわち、\(D = \sum(r_i^D)^2\)、ここで、\((r_i^D)^2\)は\(i\)番目のクラスの逸脱度への寄与を表す。残差逸脱度を例えば、\(y\)に対してにプロットすることで、\(y\)の大きさとの関連や、外れ値がないかを確認することができる。
例 ジャガイモの疫病
Johnsonら（1996）は、ジャガイモの疫病（potato blight）を前年の疫病発生の有無とその年の気象条件から、病気の発生を予測するモデルを構築した。ここでは、同じデータを用いて同様のモデルを構築してみる。
まずは、データを読み込む。
data("johnson.blight") # agridatの中に含まれているデータjohnson.blightをロードする
head(jb <- johnson.blight) # jbという短い名前のオブジェクトに代入しておく## year area blight rain.am rain.ja precip.m
## 1 1970 0 0 8 1 5.84
## 2 1971 0 0 9 4 6.86
## 3 1972 0 0 9 6 47.29
## 4 1973 0 0 6 1 8.89
## 5 1974 50 1 16 6 7.37
## 6 1975 810 1 10 7 5.08病気の発生（blight）は、発生面積（area, ha）が0より大きいときに1になっています。また、その年の4月5月の降雨日数（rain.am, days）、7月8月の降雨日数（rain.ja, days）、および、5月の温度が5度以上のときの降雨量（precip.m, mm）のデータもあります。このデータをもとに、ロジスティック回帰を行う。
この例では、ベルヌーイ分布型のロジスティック回帰となる。そのため、独立変数を2列にする必要はない。
model <- glm(blight ~ rain.am + rain.ja + precip.m, data = jb, family = binomial)
summary(model)##
## Call:
## glm(formula = blight ~ rain.am + rain.ja + precip.m, family = binomial,
## data = jb)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.37994 2.85978 -2.231 0.0257 *
## rain.am 0.55574 0.23073 2.409 0.0160 *
## rain.ja 0.21433 0.17084 1.255 0.2096
## precip.m -0.08568 0.07841 -1.093 0.2745
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 34.617 on 24 degrees of freedom
## Residual deviance: 20.029 on 21 degrees of freedom
## AIC: 28.029
##
## Number of Fisher Scoring iterations: 5\(\beta\)の推定値と標準誤差、Wald検定の\(p\)値が与えられている。切片は有意にゼロから隔たっている。また、変量rain.am（4月5月の降雨日数）も有意である。
これは信頼区間にも反映されており、\(\beta_0,\beta_1\)だけが0を含んでいない。同様にオッズ比である\(\exp(\beta_0), \exp(\beta_1)\)も1を含まない。
confint(model)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -13.5098823 -1.83263904
## rain.am 0.1934111 1.14976021
## rain.ja -0.1019369 0.60070028
## precip.m -0.2906831 0.03141216exp(confint(model))## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 1.357478e-06 0.1599908
## rain.am 1.213382e+00 3.1574357
## rain.ja 9.030866e-01 1.8233952
## precip.m 7.477526e-01 1.0319107次の図は、観察された疫病の有無（0または1）とそのロジスティック回帰への当てはめを表す。
lin <- predict(model, type = "link")
x <- seq(-6, 6, 0.1)
plot(x, exp(x) / (1 + exp(x)), col = "red", type = "l", xlab = "Linear Predictor", ylab = "Probability of potato blight")
points(lin, jb$blight)
次に、\(\hat p_i\)を得て、逸脱度を計算する。
phat <- predict(model, type = "response")
devi <- -2 * sum(jb$blight * log(phat + .Machine$double.eps) + (1 - jb$blight) * log(1 - phat + .Machine$double.eps))
devi # directly## [1] 20.02893deviance関数を用いて逸脱度を表示できる。また、modelからも簡単に引き出せる。
deviance(model) # glm output## [1] 20.02893model$deviance## [1] 20.02893\(G\)の値を計算してみる。
model0 <- glm(blight ~ 1, data = jb, family = binomial)
model0$deviance - model$deviance## [1] 14.58842なお、model内に帰無モデルの逸脱度の情報が既にふくまれている。
model$null.deviance - model$deviance## [1] 14.58842残差（ordinary, deviance）を計算して、\(\hat p\)に対してプロットしてみます。
op <- par(mfrow = c(1,2))
# ordinary
ro <- jb$blight - phat
plot(phat, ro)
abline(h = 0)
# deviance
rdev <- sign(jb$blight - phat) * sqrt(-2 * jb$blight * log(phat + .Machine$double.eps) - 2 * (1 - jb$blight) * log(1 - phat + .Machine$double.eps))
plot(phat, rdev)
abline(h = 0)
par(op)全サンプルの残差逸脱度の和が、逸脱度になることを確認してみます。
# model deviance
sum(rdev^2)## [1] 20.02893deviance(model)## [1] 20.02893プロビット回帰と相補log-log（Complementary log-log）回帰
ロジスティック回帰では、累積分布関数\(F\)として \[ \hat p_i = F(l_i) = \frac {\exp(l_i)}{1 + \exp(l_i)} \] を用いて、モデルの線形部分\(l_i = b_0 + b_1 x_{i1} + \cdots + b_{p-1}x_{i,p-1}\)と確率\(p_i\)を結びつけた。
プロビット回帰では、累積分布関数\(F\)として標準正規分布の累積分布関数 \[ \hat p_i = \Phi (l_i) \] が用いられる。プロビット回帰の欠点は、近似法や数値アルゴリズムはあるものの、リンク関数となる累積分布の逆関数\(\Phi^{-1}\)が陽的に表現できないことである。
相補log-log回帰では、累積分布関数 \[ F(l_i) = 1 - \exp(-\exp(l_i)) \] が用いられる。
線形部分\(l_i\)がプロビット回帰に当てはめらられと、その確率はそれぞれ以下のように推定できる。 \[ \hat p_i = \Phi(l_i) \] また、相補log-log回帰に当てはめられると、その確率はそれぞれ以下のように推定できる。 \[ \hat p_i = 1 - \exp(-\exp(l_i)) \]
Rでは、プロビットや相補log-logリンクを以下の例に示すように用いることができます。
例 Blissのデータ
様々な濃度の二硫化炭素ガスに5時間さらされた後に死亡した甲虫の数を計測した実験結果がある（Bliss 1935、Vidakovic 2011より引用）。このデータセットを用いて、3つの手法を比較してみる。
dose <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
beetle <- c(59, 60, 62, 56, 63, 59, 62, 60)
killed <- c(6, 13, 18, 28, 52, 53, 61, 60)
# Logit
model.logit <- glm(cbind(killed, beetle - killed) ~ dose, family = binomial)
summary(model.logit)##
## Call:
## glm(formula = cbind(killed, beetle - killed) ~ dose, family = binomial)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -60.717 5.181 -11.72 <2e-16 ***
## dose 34.270 2.912 11.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.202 on 7 degrees of freedom
## Residual deviance: 11.232 on 6 degrees of freedom
## AIC: 41.43
##
## Number of Fisher Scoring iterations: 4# Probit
model.probit <- glm(cbind(killed, beetle - killed) ~ dose, family = binomial(link = "probit"))
summary(model.probit)##
## Call:
## glm(formula = cbind(killed, beetle - killed) ~ dose, family = binomial(link = "probit"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -34.935 2.648 -13.19 <2e-16 ***
## dose 19.728 1.487 13.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.20 on 7 degrees of freedom
## Residual deviance: 10.12 on 6 degrees of freedom
## AIC: 40.318
##
## Number of Fisher Scoring iterations: 4# Clog-log
model.cloglog <- glm(cbind(killed, beetle - killed) ~ dose, family = binomial(link = "cloglog"))
summary(model.cloglog)##
## Call:
## glm(formula = cbind(killed, beetle - killed) ~ dose, family = binomial(link = "cloglog"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -39.572 3.240 -12.21 <2e-16 ***
## dose 22.041 1.799 12.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.2024 on 7 degrees of freedom
## Residual deviance: 3.4464 on 6 degrees of freedom
## AIC: 33.644
##
## Number of Fisher Scoring iterations: 4相補log-log回帰の残差逸脱度が最も低い（モデルがよく当てはまっている）。
これらを、グラフを描いて比較をしてみる。
plot(dose, killed / beetle, xlim = c(1.5, 2), ylim = c(0,1))
df <- data.frame(dose = seq(1.5, 2, 0.01))
p.logit <- predict(model.logit, df, type = "response")
lines(df$dose, p.logit, col = "red")
p.probit <- predict(model.probit, df, type = "response")
lines(df$dose, p.probit, col = "blue")
p.cloglog <- predict(model.cloglog, df, type = "response")
lines(df$dose, p.cloglog, col = "green")
曲線からの逸脱が、相補log-logモデル（緑色の曲線）で最も小さい。
ポアソン回帰
ポアソン回帰は、多数の試行の中で稀に起こる事象の回数\(y = {0,1,2,3,...}\)をモデル化するための手法である。典型的な例としては、異常な有害事象、事故、まれな病気の発生、特定の期間中の機器の故障などの回数を予測するモデルが考えられる。ポアソン分布に従う確率変数\(Y \sim Poi(\lambda)\)の確率質量関数（probability mass function）は \[ f(y) = P(Y = y) = \frac{\lambda^y}{y!} \exp(-\lambda), \:\: y = 0,1,2,3,... \] であり、平均と分散は\(\lambda\) (\(\lambda>0\))となる。
ここで、\(n\)個の計数\(y_i, \:\:i=1,...,n\)が観察されており、それぞれの計数が共変量\(x_i, \:i=1,...,n\)の値に対応して変化しているとする。典型的なポアソン回帰は、以下のように定式化できる。すなわち、 \[ y_i \sim Poi(\lambda_i) \] \[ \log(\lambda_i) =\beta_0+\beta_1 x_i,\:i=1,...,n \] なお、\(\lambda_i\)と線形部分\(\beta_0+\beta_1 x_i\)の関係については、\(\lambda_i\)が正になるのであれば、様々な形に拡張できる。例えば、\(\lambda_i\)は\(p-1\)の共変量と\(p\)個のパラメータをもつ線形表現にリンクすることもできる。すなわち、 \[ \log(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_{p-1} x_{i,p-1}, \:\:i = 1,...,n \]
上記のモデルでは、\(\lambda_i = E(y_i|X = x_i)\)が期待値であり、その対数が\(\log E(y_i|X = x_i) = \beta_0 + \beta_1 x_i\)として表される。このとき、共変量 \(x_i\)が1単位分大きくなり\(x_i + 1\)になると、 \[ \log E(y_i|X = x_i + 1) = \beta_0 + \beta_1 x_i + \beta_1 = \log E(y_i|X = x_i) + \beta_1 \] となる。したがって、\(\beta_1\)は、共変量が1増加したときの、対数比率（log rate）の増加量として解釈することができる。逆に、\(\exp(\beta_1)\)は、比率の割合 \[ \exp(\beta_1) = \frac {E(y_i|x_i + 1)}{E(y_i|x_i)} \] となる。
モデルによって推定される平均応答は、\(\hat y_i=\exp(b_0 + b_1 x_i)\)と表される（ここで、\(b_0, b_1\)は、\(\beta_0, \beta_1\)の推定値）。
モデルの逸脱度（deviance）\(D\)は以下のように定義される。 \[ D = -2 \sum_{i=1}^n \biggl[y_i \log \frac {\hat y_i}{y_i} - (\hat y_i - y_i)\biggr] \\ = 2 \sum_{i=1}^n \biggl[y_i \log \frac {y_i}{\hat y_i} - (y_i - \hat y_i)\biggr] \] ここで、\(y_i = 0\)のとき\(y_i \log y_i = 0\)となる。ロジスティック回帰のように、逸脱度はモデルの当てはまりの良さの指標であり、ポアソン回帰では、自由度\(n-p\)の\(\chi^2\)分布に従う。
また、残差逸脱度（deviance residuals）は、 \[ r_i^{dev} = \mbox{sign}(y_i - \hat y_i) \sqrt{2 \left[y_i \log \frac {y_i}{\hat y_i} - (y_i - \hat y_i)\right]} \] と定義され、やはり\(D = \sum_{i = 1}^n (r_i^{dev})^2\)を満たす。
例 コットンボールの数
agridatパッケージに含まれるSilviaら（2012）の綿の実のデータ（silvia.cotton）をもとに、綿の実（cotton boll）の数と、人工的な適葉（artificial defoliation）の割合、および、生育ステージ（stage）の関係をポアソン回帰でモデル化してみる。
まずは、データを読み込む。
data("silva.cotton")
sc <- silva.cotton図を描いて生育ステージごとに変数間の関係を確認してみましょう。
tmp <- sc %>% filter(stage == "vegetative")
boxplot(tmp$bolls ~ tmp$defoliation)
tmp <- sc %>% filter(stage == "flowerbud")
boxplot(tmp$bolls ~ tmp$defoliation)
tmp <- sc %>% filter(stage == "blossom")
boxplot(tmp$bolls ~ tmp$defoliation)
tmp <- sc %>% filter(stage == "boll")
boxplot(tmp$bolls ~ tmp$defoliation)
tmp <- sc %>% filter(stage == "bollopen")
boxplot(tmp$bolls ~ tmp$defoliation)
適葉率が上がると、綿の実の数も減っていくことが多いことがわかる。ただし、生育ステージによってその影響は変化する。
今回は綿の実が実ったときのデータをもとに、ポアソン回帰を行ってみる。ロジスティック回帰と同様にglm関数を用いて回帰を行うことができる。
sc.boll <- sc %>% filter(stage == "boll")
model <- glm(bolls ~ defoliation, data=sc.boll, family = poisson(link = "log"))
summary(model)##
## Call:
## glm(formula = bolls ~ defoliation, family = poisson(link = "log"),
## data = sc.boll)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.541567 0.119712 12.88 < 2e-16 ***
## defoliation -0.007244 0.002229 -3.25 0.00115 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 29.824 on 49 degrees of freedom
## Residual deviance: 19.035 on 48 degrees of freedom
## AIC: 173.38
##
## Number of Fisher Scoring iterations: 4\(\beta_0\)も\(\beta_1\)もいずれも有意である。\(b_1 = -0.007\)と推定されており、共変量が1増加すると対数比率が0.007減少することがわかる。
次に、\(\beta_0\)も\(\beta_1\)の信頼区間を計算してみる。
confint(model)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 1.29933095 1.769004235
## defoliation -0.01165423 -0.002906803\(beta_0, beta_1\)のいずれの信頼区間も0をまたいでいない。
次に、\(y\)の推定値を求めて、残差逸脱度を計算してみます。
yhat <- fitted(model) # get estimated values
y <- sc.boll$bolls
rdev <- sign(y - yhat) * sqrt(2 * y * log(y/(yhat + .Machine$double.eps))-2*(y - yhat))
# .Machine$double.epsはyがゼロのときのための対策です。
plot(y, rdev, main = "deviance")
abline(h = 0)
残差逸脱度の二乗和が、モデルの逸脱度になっていることを確認します。
sum(rdev^2)## [1] 19.03451(D <- model$deviance)## [1] 19.03451\(\chi^2\)検定を用いて逸脱度の検定をしてみる。
pval <- 1 - pchisq(D, length(y) - 2)
pval## [1] 0.9999424逸脱度は有意ではなく、モデルによく当てはまっていることがわかる。
最後にGの値を計算し、Gの検定をしてみる。
(G <- model$null.deviance - model$deviance) # the same as G, difference of deviance## [1] 10.78944pval <- 1 - pchisq(G, 1)
pval## [1] 0.001020808モデルが帰無モデル（切片だけで共変量のないモデル）である\(H_0\)のもとでは、統計量\(G\)は\(df = p - 1\)の自由度をもつ。この場合、\(df = 1\)である。この検定が有意であることから、共変量defoliationは有意にモデルに寄与していることが分かる。
参考文献
以下は、本テキストを作成するにあたり参考にした。より詳しく学ぶためにもお薦めの本です。
- Vidakovic B (2011) Statistics for Bioengineering Sciences With MATLAB and WinBUGS Support. Springer.