死亡率などの解析

補対数対数回帰

従属変数が 0 か 1 のように 2 つの値をしか取らないモデルは、誤差構造として二項分布が用いられる。たとえば、ある薬剤 A を微生物を培養した培地に追加したとき、 死亡したコロニーの割合は薬剤 A の濃度によってどのように変化するかを調べる場合、誤差構造として二項分布が用いられる。また、年齢が高くなると死亡率はどのように変化するかを調べる際にもよく使われる。

サンプルデータとして、ある生物の年齢を 20 段階に分け、それぞれの段階における死亡率を調べたところ、以下のようになったとする。(※ このデータは確率となっているが、カウントデータがあればカウントデータのままで解析したほうが良い。)

age <- 1:20
mortality <- c(0.000, 0.018, 0.002, 0.005, 0.006,
               0.013, 0.020, 0.027, 0.039, 0.040,
               0.103, 0.145, 0.210, 0.301, 0.454,
               0.689, 0.790, 0.820, 0.984, 1.000)

これをグラフに描くと以下のようになる。

これについてリンク関数をロジスティック関数としてモデル化し、パラメーターを推定する。推定結果を summary 関数で見える。

cloglog.model <- glm(mortality ~ age, family = binomial("cloglog"))
## Warning message:
## In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
summary(cloglog.model)
## Call:
## glm(formula = mortality ~ age, family = binomial("cloglog"))
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.224436  -0.023862  -0.001975   0.025501   0.216775  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -7.0790     2.9962  -2.363   0.0181 *
## x             0.4384     0.1868   2.347   0.0189 *
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13.75289  on 19  degrees of freedom
## Residual deviance:  0.16584  on 18  degrees of freedom
## AIC: 9.2222
## 
## Number of Fisher Scoring iterations: 7

パラメーター β1 = -7.0790、β2 = 0.4384 であるので、このパラメーターでモデルをグラフにすると以下のようになる。

plot(age, mortality, xlim = c(0, 20), ylim = c(0, 1), xlab = "age", ylab = "mortality")
par(new = T)
curve(1 - exp(- exp(0.4384 * x - 7.0790)), type = "l",
      xlim = c(0, 20), ylim = c(0, 1), xlab = "", ylab = "")