Dobson Chap.7 1-4
Yuta NAGANO
11/29/2019
本章では、応答変数が２値で得られる一般化線形モデルを扱う。
7.1 確率分布
- 2値で得られるデータは患者の生死や存在の有無などが挙げられるが、以下では成功(TRUE)と失敗(FALSE)と称する。
- 2値の確率変数として以下を定義する。 \[Z = \left\{\begin{array}{II}1 & TRUE \\ 0 & FALSE \end{array}\right.\] この時、\(Pr(Z = 1) = \pi\)、\(Pr(Z = 0) = 1-\pi\)とする。
- このような\(n\)個の変数\(Z_1, ..., Z_n\)が存在し、互いに独立かつ\(Pr(Z_j = 1) = \pi_j\)であるならば、同時確率は、以下のように表現できる。 \[\prod^n_{j=1}\pi_j^{z_j}(1-\pi_j)^{1-z_j} = exp\biggl[\sum^n_{j=1}z_jlog\Bigl(\frac{\pi_j}{1-\pi_j}\Bigr)+\sum^n_{j=1}log(1-\pi_j)\biggr] \tag{7.1}\] この式からわかるように、指数型分布族に属することがわかる\((\; f(y; \theta) = exp[a(y)b(\theta) + c(\theta) + d(y)]\;)\)。
- もし、\(\pi_j\)が全ての\(j\)について等しい時は、以下を定義する。 \[Y = \sum^n_{j=1}Z_j\] そうすれば、\(Y\)は\(n\)回の施行において成功した回数を表す。このような確率変数\(Y\)は二項分布に従う。 \[Pr(Y = y) = \begin{pmatrix}n\\y\end{pmatrix}\pi^y(1-\pi)^{1-y}\;,\;\;\; y = 0, 1, ..., n \tag{7.2}\]
- より一般的な場合として、\(N\)個の独立な確率変数\(Y_1, ..., Y_N\)が\(N\)個のサブグループにおける成功数を示す場合、すなわち\(Y_i ~ binomial(n_i, \pi)\)ならば、対数尤度関数は以下のように表される。 \[l(\pi_1, ..., \pi_N; y_1, ..., y_N) = \sum^N_{i=1}\biggl[y_ilog\Bigl(\frac{\pi_i}{1-\pi_i}\Bigr)+n_ilog(1-\pi_i)+log\begin{pmatrix} n_i\\y_i \end{pmatrix}\biggr] \tag{7.3}\]
7.2 一般化線形モデル
一般化線形モデルのモチベーション：
各サブグループにおける成功割合\(P_i=Y_i/n_i\)を、そのサブグループを特徴づけるいくつかの説明変数に関して表したい。
- 具体的には\(E(P_i)=E(Y_i)/n_i=\pi_i\)なので、確率\(\pi_i\)を以下のようにモデル化する。 \[g(\pi_i)=\boldsymbol{x}^T_i\boldsymbol{\beta}\] ここで、\(\boldsymbol{x}_i\)は説明変数ベクトルであり、\(\boldsymbol{\beta}\)はパラメータベクトル、\(g\)はリンク関数とである。
- もっとも単純なケースとしては、線形モデルがある。 \[\pi_i=\boldsymbol{x}^T_i\boldsymbol{\beta}\] このようなモデルが用いられる場合もあるが、\(\boldsymbol{x}_i\)の値によっては、\(\boldsymbol{x}_i^T\boldsymbol{\beta}\)が０以下や１以上になってしまうので確率を表すモデルとしては、不適当である。
- \(\pi_i\)を区間[0,1]に制限するために、次のように累積分布関数を用いる方法がよく使われる。 \[\pi_i = \int^{t_i}_{-\infty}f(s)ds\] ここでは、\(f(s)\geq 0\)かつ\(\int^{\infty}_{-\infty}f(s)ds=1\)である。この時の確率密度関数\(f(s)\)は、許容値分布と呼ばれている。よく使われる許容値分布の例は以下で取り上げる。
7.3 用量反応モデル
- ２値データに対して初めて、回帰のようなモデルが使われたのは、バイオアッセイの分野である。
- このようなデータは定性反応データと呼ばれており、分析の目的は失敗率(毒の量に対する死亡率とか)\(\pi\)を用量\(x\)の関数としてモデル化すること。
- もし、\(f(s)\)が区間\([c_1, c_2]\)上の一様分布ならば、 \[f(s) = \left\{\begin{array}{II} \frac{1}{c_2-c_1}&c_1\leq s\leq c_2 \\ 0 & others \end{array}\right.\] となるので、\(\pi\)は\(x\)の関数として \[\pi = \int^x_{c_1} f(s)ds = \frac{x-c_1}{c_1-c_2}, \;\;\; c_1\leq x\leq c_2\] と表される。
- 結局、\(\pi = \beta_1+\beta_2x\)、\(c_1 \leq x \leq c_2\)、ただし、 \[\beta_1 = \frac{-c_1}{c2_c1}, \;\;\; \beta_2=\frac{1}{c_2-c_1}\] というリンク関数\(g\)を恒等関数とした線形モデルになる。
- しかし、\(c_1 \leq x \leq c_2\)という余分な条件があるので、GLMの標準的なパラメータ推定の方法は適用できない。
実用上このモデルは広く使われていない
プロビットモデル
- バイオアッセイのデータに最初に使われたのは、プロビットモデルである。
このモデルでは正規分布を許容値分布として用いる。 \[\begin{align} \pi &= \frac{1}{\sigma\sqrt{2\pi}}\int^x_{-\infty}exp\biggl[ -\frac{1}{2}\Bigl(\frac{s-\mu}{\sigma}\Bigr)^2\biggr] \\ &= \Phi^{-1}\Bigl(\frac{x-\mu}{\sigma}\Bigr)\end{align}\] ただし、\(\Phi\)は標準正規分布\(N(0,1)\)の分布関数である。 - このモデルは\(\beta_1=-\mu/\sigma\)、\(\beta_2=1/\sigma\)とおけば、 \[\Phi^{-1}(\pi) = \beta_1+\beta_2x\] となるので、リンク関数\(g\)を標準正規分布の逆累積分布関数としたGLMである。
ロジスティックモデル
- プロビットモデルと非常に似た計算結果を返してくるが、数値計算上いくぶん扱い安いモデルとしてロジスティックモデルあるいはロジットモデルと呼ばれるものがある。
- ロジスティックモデルでは、許容値分布は、以下で与えられる。 \[f(s)=\frac{\beta_2 exp(\beta_1+\beta_2s)}{[1+exp(beta_1+\beta_2s)]^2}\] これから、以下と計算される。 \[\pi = \int^x_{-\infty}f(s)ds=\frac{exp(\beta_1+\beta_2x)}{1+exp(\beta_1+beta_2x)}\] これより、 \[log\Bigl(\frac{\pi}{1-\pi}\Bigr)=\beta_1+\beta_2x\] となり、リンク関数として\(g(\pi)=log[\pi/(1-\pi)]\)をもったGLMであることがわかる。
- \(log[\pi/(1-\pi)]\)はしばしばロジット関数と呼ばれるが、これは対数オッズという自然な解釈をもっている。
- ロジスティックモデルとプロビットモデルでは、\(f(S)\)の形状はよく似ているが、裾野における確率分布がやや異なる（らしい）。
- 用量反応データのモデル化に用いられるその他の許容値分布としては、極値分布が挙げられる。極値分布は以下の式で与えられる。 \[f(s) = \beta_2 exp[(\beta_1+\beta_2s)-sxp(\beta_1+\beta_2s)]\] これから、以下が計算される。 \[\pi = 1-exp[-exp(\beta_1+\beta_2x)]\] 変形すると、\(log[-log(1-\pi)]=\beta_1+\beta_2x\)となり、この場合のリンク関数\(log[-log(1-\pi)]\)は、complementary log-log関数と呼ばれている。
- このモデルは、\(\pi\)の値が0.5の近傍の時はロジスティックモデルやプロビットモデルによく類似するが、\(\pi\)が0や1に近い時はズレが生じる（らしい）。
7.3.1 例：カブトムシの死亡率
カブトムシに二硫化炭酸ガスを5時間浴びせた時の、ガス濃度と死亡率の関係
(kabuto <- tibble(noudo = c(1.69072, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839),
all = c(59, 60, 62, 56, 63, 59, 62, 60),
death = c(6, 13, 18, 28, 52, 53, 61, 60)))## # A tibble: 8 x 3
## noudo all death
## <dbl> <dbl> <dbl>
## 1 1.69 59 6
## 2 1.72 60 13
## 3 1.76 62 18
## 4 1.78 56 28
## 5 1.81 63 52
## 6 1.84 59 53
## 7 1.86 62 61
## 8 1.88 60 60kabuto %>% mutate(death_ratio = death/all*100) %>%
ggplot(aes(x = noudo, y = death_ratio)) + geom_point() + theme_base(base_family = "HiraKakuPro-W3") +
labs(x = "二硫化炭酸ガスの濃度", y = "カブトムシの死亡率(%)") + xlim(1.65, 1.9) + ylim(10, 100)
- まずは、ロジスティックモデルを当てはめてみる。 \[\pi_i = \frac{exp(\beta_1+\beta_2x_i)}{1+exp(\beta_1+\beta_2x_i)}\] この式から、 \[log\Bigl(\frac{\pi_i}{1-\pi_1}\Bigr) = \beta_1+\beta_2x_i\] および、 \[log(1-\pi_i) = -log[1+exp(\beta_1+\beta_2x_i)]\] が導かれる。これらと式(7.3)より、対数尤度関数は、 \[l = \sum^N_{i=1}\biggl[y_i(\beta_1+\beta_2x_i)-n_ilog[1+exp(\beta_1+\beta_2x_i)]+log\begin{pmatrix} n_i \\ y_i \end{pmatrix} \biggr]\] と表され、かつパラメター\(\beta_1\)と\(\beta_2\)に関するスコア関数は以下のように計算される。 \[\begin{align} U_1 &= \frac{\delta l}{\delta \beta_1} = \sum \Biggl[y_i-n_i\biggl[\frac{exp(\beta_1+\beta_2x_i)}{1+exp(\beta_1+\beta_2x_i)}\biggr]\Biggr] = \sum(y_i-n_i\pi_i) \\ U_2 &= \frac{\delta l}{\delta \beta_2}= \sum \Biggl[y_ix_i-n_ix_i\biggl[\frac{exp(\beta_1+\beta_2x_i)}{1+exp(\beta_1+\beta_2x_i)}\biggr]\Biggr] = \sum x_i(y_i-n_i\pi_i) \end{align}\]
- 情報行列は以下のように計算される。 \[\boldsymbol{\Im} = \begin{bmatrix} \sum n_i\pi_i(1-\pi_i) & \sum n_ix_i\pi_i(1-\pi_i) \\ \sum n_ix_i\pi_i(1-\pi_i) & \sum n_ix_i^2\pi_i(1-\pi_i)\end{bmatrix}\]
- 最尤推定値は、以下の逐次式を解くことで得られる。 \[\boldsymbol{\Im}^{(m-1)}\boldsymbol{b}^{(m)} = \boldsymbol{\Im}^{(m-1)}\boldsymbol{b}^{(m-1)}+\boldsymbol{U}^{(m-1)}\] ここで、上付き添字\((m)\)は\(m\)回目の近似を、\(\boldsymbol{b}\)は推定値ベクトルをそれぞれ意味する。
- 最終的には、6回の反復計算で収束し、以下の結果が得られる。
- \(b_1 = -60.7171, \;\;\;\; s.e.(b_1) = \sqrt{26.840} = 5.18\)
- \(b_2 = 34.27, \;\;\;\; s.e.(b_2) = \sqrt{8.481} = 2.91\)
- 逸脱度は、 \[D = 2\sum^N_{i=1}\biggl[y_ilog\Bigl(\frac{y_i}{\hat{y}_i}\Bigr)+(n-y_i)log\Bigl(\frac{n_i-y_i}{n_i-\hat{y}_i}\Bigr)\biggr]\]
R code
考えてみたけどうまくいかないです。理由はわかっているので、今後アップデートしてきます。
logistic <- function(x, b1, b2){
N <- length(x)
pi <- rep(0, length.out = N)
for(n in 1:N){
pi[n] <- exp(b1+(b2*x[n]))/(1+exp(b1+(b2*x[n])))
}
return(pi)
}
XWX <- function(data, b1, b2){
element_11 <- sum(data$all*logistic(data$noudo, b1, b2)*(1-logistic(data$noudo, b1, b2)))
element_12 <- sum(data$all*logistic(data$noudo, b1, b2)*data$noudo*(1-logistic(data$noudo, b1, b2)))
element_21 <- sum(data$all*logistic(data$noudo, b1, b2)*data$noudo*(1-logistic(data$noudo, b1, b2)))
element_22 <- sum(data$all*logistic(data$noudo, b1, b2)*data$noudo*data$noudo*(1-logistic(data$noudo, b1, b2)))
xwx <- matrix(c(element_11, element_12, element_21, element_22), ncol = 2, byrow = TRUE)
return(xwx)
}
XWz <- function(data, b1, b2){
upper <- sum(data$death*data$all*logistic(data$noudo, b1, b2)*(1-logistic(data$noudo, b1, b2)))
lower <- sum(data$death*data$all*logistic(data$noudo, b1, b2)*data$noudo*(1-logistic(data$noudo, b1, b2)))
xwz <- matrix(c(upper, lower), ncol = 1)
return(xwz)
}
b_new <- function(mat1, mat2){
b_new = solve(mat1) %*% mat2
return(b_new)
}
# b_newを収束するまで反復計算
# 初期値を与える
b1 <- 0
b2 <- 0
N <- 1000 # 適当な数字（反復計算の最大回数）
# 反復計算
for(i in 1:N){
mat1 <- XWX(kabuto, b1 = b1, b2 = b2)
mat2 <- XWz(kabuto, b1 = b1, b2 = b2)
bnew <- b_new(mat1, mat2)
if(abs(b1-bnew[1,]) <= 0.00001 & abs(b2-bnew[2,]) <= 0.00001) break
b1 <- bnew[1, ]
b2 <- bnew[2, ]
}
# 結果を見てみる
bnew## [,1]
## [1,] -347.2693
## [2,] 208.9461# GLMでやってみとく
glm(cbind(death, (all-death)) ~ noudo, data = kabuto, family = binomial) %>% summary()##
## Call:
## glm(formula = cbind(death, (all - death)) ~ noudo, family = binomial,
## data = kabuto)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5940 -0.3945 0.8331 1.2589 1.5938
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -60.721 5.181 -11.72 <2e-16 ***
## noudo 34.272 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.229 on 6 degrees of freedom
## AIC: 41.427
##
## Number of Fisher Scoring iterations: 47.4 一般ロジスティック回帰モデル
- 前節で出てきた、単純な線形ロジスティックモデル\(log[\pi_i/(1-\pi_i)] = \beta_1 + \beta_2 x_i\)は、一般ロジスティック回帰モデル \[logit(\pi_i) = log\biggl(\frac{\pi_i}{1-\pi_i}\biggr) = \boldsymbol{x}^T_i\boldsymbol{\beta}\] の特別な場合に相当する。
- ここに\(\boldsymbol{x}_i\)は連続今日変量の測定値や名義尺度の水準を表すダミー変数から構成されるベクトルである。
- このモデルは、2値反応といくつかの説明変数を含むデータを解析するのに広く使われており、連続反応を解析するために使われる重回帰や分散分析に似た強力な方法である。
- パラメータ\(\boldsymbol{\beta}\)の最尤推定値、したがって確率\(\pi_i = g^{-1}(\boldsymbol{x}_i^T\boldsymbol{\beta})\)の最尤推定値を得るには、対数尤度関数 \[l(\boldsymbol{\pi}; \boldsymbol{y}) \sum^N_{i=1}\biggl[y_ilog\pi_i + (n_i-y_i)log(1-\pi_i) + log\begin{pmatrix}n_i\\y_i\end{pmatrix}\biggr]\] を最大化すれば良い。
- この推定のプロセスは、データが各共変量パターン(説明変数の値が全て同一であるような観測値群)について、グループ化されていてもされていなくても変わらない。
- グループ化されていない場合は、\(Y_i\)は0または1をとる2値変数であり、共変量パターンは別々のものとして扱われる。
- グループ化されていれば、\(Y_i\)は第\(i\)共変量パターンに対する成功数を表し、二項分布に従う。
- 逸脱度は、 \[D = 2\sum^N_{i=1}\biggl[y_ilog\Bigl(\frac{y_i}{\hat{y}_i}\Bigr)+(n_i-y_i)log\Bigl(\frac{n_i-y_i}{n_i-\hat{y}_i}\Bigr)\biggl]\] と表されるが、これは、 \[D = 2\sum o \;log\frac{o}{e}\] という形をとる。 ここで、\(o\)は観測された度数\(y_i\)または\((n_i-y_i)\)を、\(e\)は推定された期待度数\(\hat{y}_i = n_i\hat{\pi}_i\)または、\((n-\hat{y}_i)=(n-n_i\hat{\pi}_i)\)をそれぞれ表している。
- 逸脱度\(D\)は、局外パラメータを１つも含まないので、モデルの適合度の検定は、近似 \[D \sim \chi^2(N-p)\] を利用することによって直接できる。
- N：共変量パターンの総数
- p：推定されるパラメータの個数
- 統計的推測に使われる推定法や標本分布は漸近的な結果に依存している。
7.4.1 例：葯からの不定胚形成
- チョウセンアサガオのデータ
(asagao <- tibble(treatment = c("cont", "cont", "treat", "treat"),
type = c("futei", "yaku", "futei", "yaku"),
"40" = c(55, 102, 55, 76),
"150" = c(52, 99, 50, 81),
"350" = c(57, 108, 50, 90)) %>% gather(key = power, value = n, -treatment, -type) %>%
mutate_at(vars(power), funs(as.numeric(.))))## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.## # A tibble: 12 x 4
## treatment type power n
## <chr> <chr> <dbl> <dbl>
## 1 cont futei 40 55
## 2 cont yaku 40 102
## 3 treat futei 40 55
## 4 treat yaku 40 76
## 5 cont futei 150 52
## 6 cont yaku 150 99
## 7 treat futei 150 50
## 8 treat yaku 150 81
## 9 cont futei 350 57
## 10 cont yaku 350 108
## 11 treat futei 350 50
## 12 treat yaku 350 90asagao %>% spread(key = type, value = n) %>% mutate(ratio = futei/yaku) %>%
ggplot(aes(x = log(power), y = ratio)) + geom_point(aes(color = treatment)) +
theme_base(base_family = "HiraKakuPro-W3") + labs(x = "Log(遠心力)", y = "不定胚の割合") +
scale_color_hue(name = "貯蔵条件", labels = c(cont = "対照群", treat = "処理群"))
- 以下の３つのロジスティックモデルを検討する。
- モデル1：logit \(\pi_{jk} = \alpha_j + \beta_jx_k\)（異なった切片、傾き）
- モデル2：logit \(\pi_{jk} = \alpha_j + \beta x_k\)（異なった切片、共通の傾き）
- モデル3：logit \(\pi_{jk} = \alpha + \beta x_k\)（共通の切片、傾き）
- 最尤推定値は以下のようになる。
- モデル1
- \(\alpha_1 = 0.234(0.628)\)
- \(\alpha_2 - \alpha_1 = 1.977(0.998)\)
- \(\beta_1 = -0.023(0.127)\)
- \(\beta_2 - \beta_1 = -0.319(0.199)\)
- \(D_{model1} = 0.028\)
- モデル2
- \(\alpha_1 = 0.877(0.487)\)
- \(\alpha_2-\alpha_1 = 0.407(0.175)\)
- \(\beta = -0.155(0.097)\)
- \(D_{model2} = 2.619\)
- モデル3
- \(\alpha = 1.021(0.481)\)
- \(\be4ta = -0.148(0.096)\)
- \(D_{model3} = 8.092\)
- モデル1
- 「傾きが処理群と対照群でおなじ」という帰無仮説を検定したい。
- モデル２とモデル１の逸脱度の差\(D_{model2}-D_{model1} = 2.591\)を使う。
- 帰無仮説が正しいなら、自由度1のカイ二乗分布に従うが、p値は0.1~0.2だったので、棄却できない。 →多分、FALSE-Negativeがあるんじゃないのか？