被説明変数が離散値に基づく割合データのモデリング

　ポワソン分布のGLMの時のように、他にもカウントデータに基づくモデリングが必要になる場面がある。例えば、合計でn個調べたうち、y個結実していたなどの結実率だ。nもyも共に0以上の整数であり、nが11.3とか、yは6.9とかになったりはしない。けれども、結実率といったら、y/nをモデリングするものなのではないか、と考えるかもしれない。例えば、n = 11, 20, 15で、それぞれy = 2, 7, 6みたいなデータの時に、y/n = 0.181, 0.35, 0.4を解析することを考える。するとこれらの値は連続値になるので正規分布とかを考えるだろう。しかしながら、割り算の解析は統計上微妙な部分があり、あまり推奨されない。それよりも、もっと妥当な解析方法を考えることができる。つまり、n個中y個というような値を記述する確率分布を考えることができるのだ。

二項分布：0以上n以下の離散値をとる確率分布

　結実率のようなカウントデータ/カウントデータの割合を解析するときによく用いられる確率分布が二項分布binomial distributionである。この分布は0以上n以下の離散値（0, 1, 2……, n）を台とする分布である。二項分布の形を決めるパラメータは試行回数nと(成功)確率pである。パラメータの名前は一般的に使われるものを採用しただけなので、いったんはその名前のことは忘れてもらってよい。なお、n = 1のとき、特別にベルヌーイ分布Bernoulli distributionと呼ばれ、ベルヌーイ分布に従う試行をベルヌーイ試行と呼んでいる（追記1）。二項分布はベルヌーイ試行をn回行うときの分布、と言い換えることもできる。二項分布を表す確率質量関数Probability mass functionはnとpをパラメータとして、下記のとおりである。

(n x)はnCx、つまり組み合わせを表す記号である。二項分布の主張は、例えば、表が出る確率pのコインをn回投げて、表がちょうどx回出る確率を記述している。この定義に基づくと、この分布の平均はnp、分散はnp(1 - p)とあらわせる。

続いて分散は下記のとおりである。

つまり、二項分布はnが大きくなれば平均と分散が大きくなり、pが大きくなるほど平均が大きく、0.5に近づくほど分散が大きくなることがわかる。続いて、具体的に乱数を生成してデータの様子を確認しよう。rbinom関数を使うことで、二項分布に従うデータを生成できる。引数nは生成する個数、sizeは上記の試行回数、probは上記の成功確率である。

------------------------------------------------------

library(plotn)

library(viridis)

rbinom(n = 30, size = 30, prob = 0.5)

##  [1] 11 22 15 13 17 18 15 16 13 10 17 17 15 15 13 15 16 12 17 15 13 20 17 16 16

## [26] 15 14 11 15 14

------------------------------------------------------

このように乱数を生成してみるとわかるが、生成されるデータはsize × prob = 15周りが多いことがわかるだろう。続いて、probのパラメータによって確率分布がどうなるかをいくつか確認してみる。試しに、p = 0.1, 0.5, 0.7を指定してみよう。

------------------------------------------------------

dd <- data.frame(X = 0:15)

dd$prob01 <- dbinom(dd$X, size = 15, prob = 0.1)

dd$prob05 <- dbinom(dd$X, size = 15, prob = 0.5)

dd$prob07 <- dbinom(dd$X, size = 15, prob = 0.7)

plotn(dd[,1], dd[,2:4], mode = "m", legend = T, leg.title = "prob",

      leg.lab = c(0.1, 0.5, 0.7), leg.sp = 5, ylab = "確率質量分布")#図1の描画

------------------------------------------------------

p < 0.5では左寄り、p > 0.5では右寄りの分布になることが確認できる。p = 0.5のとき、分布は左右対称である。次はさまざまなnを指定してみよう。

------------------------------------------------------

dd <- data.frame(X = c(0:5,0:10,0:15), 

                 size = c(rep(5, 6), rep(10, 11), rep(15, 16)),

                 prob = c(dbinom(0:5, size = 5, prob = 0.5),

                          dbinom(0:10, size = 10, prob = 0.5),

                          dbinom(0:15, size = 15, prob = 0.5))

)

plotn(prob ~ X, data = dd, group = "size", mode = "m", legend = T, leg.title = "size", 

      leg.lab = c(5, 10, 15), leg.sp = 5, ylab = "確率質量分布")#図2の描画

------------------------------------------------------

nが大きくなるほど、なだらかな山形になることがわかるだろう。実際、nが十分に大きければ二項分布は正規分布への近似が可能になることが知られている。

二項分布を使ったGLM

　では、二項分布を使ったGLMについて考察してゆこう（こういう解析の場合、二項回帰とはあまり言わない気がする……）。二項分布を使う場合、Rのプログラム上では、リンク関数はデフォルトでロジット関数logit functionを使うことになっている。ロジット関数の逆関数はロジスティック関数logistic functionである。ロジスティック関数とロジット関数は聞きなれないと思うので、まずはこれらの関数について紹介する。パラメータβと説明変数xを下記のように定義した時、ロジスティック関数は下記のような関数である。

ロジスティック関数はxが正負の無限大で定義されるが、yは0から1までの値しかとらない。xの係数が正ならばxが小さくなればyは0に近づき、xが大きくなればyは1に近づく。xの係数が負ならばxが小さくなればyは1に近づき、xが大きくなればyは0に近づく。この性質から、ロジスティック関数は確率を記述するのに向いている。このロジスティック関数を線形予測子について解くと、それはロジット関数になっている。

リンク関数は元の関数の式を線形予測子に変換する式、リンク関数の逆関数が元の関数の形であるから、二項分布とロジットリンクをGLMで使うということは、元の関数はロジスティック関数を仮定していることになる。二項分布のGLMでは成功確率pを推定するのでロジスティック関数との相性が良い。これまでと異なり、分布の平均値を推定しているわけではないことに注意が必要だ（後述するように二項分布の平均的な確率を推定している）。ではこのような状況において、尤度を最大化するようなパラメータβを求める。データの生成過程を図示して様子をつかむことにする。 

例えば、一つの説明変数xが存在する時、上記のような状況でデータが生成されていると考える。二項分布の成功確率pはロジットリンクであることから、説明変数に対してロジスティック関数の挙動を示す。予想される確率に対して、被説明変数はその確率と同時に得られている試行回数niをもとに二項分布に従ってデータが生成されている。そう考えると、xiとniが得られた時のyiが生じる条件付き確率は上記のように表現できる。今回はデータセットがm組（nは試行回数の文字で使ってしまったので）あると考えて、それらすべての観測値xi、yiについて掛け合わせたものが尤度なのだから、尤度は下記の通りになる。

そして対数尤度は下記のようになり、それの式を整理する。

さて、この対数尤度をパラメータβの関数と考えて、対数尤度を最大にするようなβを求めるのが最尤法だった。しかし、この関数の微分係数が0となるようなβは一般に解析的には解けない。なので以前紹介したような近似などを用いて数値的に解くことになる。対数尤度が説明変数xを含まなければ下記のようにβ0を解析的に解くことができる。

すると、対数尤度を最大にするp = 1/(1+exp(b0))は(yの総和)/(nの総和)を満たすことがわかる。つまり、平均的な確率を計算している。この推定値に基づけば、(nの平均)×pは丁度、yの平均を推定していることになる。

一般化線形モデルの検出力

　ではRを使って、GLMの検出力に関して検討しよう。以下のように、線形予測子の切片の係数は0ではないが、説明変数xの係数は0である場合をまず考える。リンク関数はロジットリンクを想定するので、yを生成するときの引数probはロジスティック関数で指定する。試行回数sizeも20～30でランダムに生成されるとしよう。データを生成した後で、描画してみよう。

------------------------------------------------------

b0 <- 2

b1 <- 0

m <- 30

x <- runif(m, 0, 30)

n <- sample(20:30, size = m, replace = T)

y <- rbinom(n = m, size = n, prob = 1/(1+exp(-(b1 * x + b0))))

plotn(y/n ~ x)#図4の描画

------------------------------------------------------

データの生成過程通り、xに対して特に関係のなさそうなデータが生成された。では、このデータをGLMで解析する。ここでは二項分布とロジットリンクを指定するので、binomial(link="logit")と指定する。二項分布を指定する場合、被説明変数の設定が他と比べて違うことに注意する。被説明変数はyの他に、試行回数nから成功回数yを引いた値n-yをともに組み込む。そのために、cbind(y, n-y)としている。こうすることで試行回数nも解析に組み込めるのだ。併せて、回帰曲線も描いてみよう。

------------------------------------------------------

fit1 <- glm(cbind(y, n-y) ~ x, family = binomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.4342  -0.9218   0.1589   0.9229   1.6425  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  2.27702    0.20342  11.193   <2e-16 ***

## x           -0.01770    0.01186  -1.492    0.136    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for binomial family taken to be 1)

## 

##     Null deviance: 28.099  on 29  degrees of freedom

## Residual deviance: 25.879  on 28  degrees of freedom

## AIC: 110.12

## 

## Number of Fisher Scoring iterations: 4

b0e <- coef(fit2)[1,1]

b1e <- coef(fit2)[2,1]

plotn(y/n ~ x)

xx <- seq(min(x), max(x), length = 200)

yy <- 1/(1+exp(-(b0e + b1e * xx)))

overdraw(points(xx, yy, type = "l")) #図5の描画

------------------------------------------------------

GLMの結果をsummary関数を通して得られる結果の解釈も線形回帰の時とそれほど変わらない。今回の場合、線形予測子の切片の推定値は2.28、xの係数の推定値は-0.018であり、初めの設定と大きく変わらないことがわかる。Devianceなどの簡単な説明は前回参照。AICは当該ページ参照。

　では、上記の設定でデータを10000回生成して、検出力を確認してみる。

------------------------------------------------------

b0 <- 2

b1 <- 0

m <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(m, 0, 30)

  n <- sample(20:30, size = m, replace = T)

  y <- rbinom(n = m, size = n, prob = 1/(1+exp(-(b1 * x + b0))))

  fit1 <- glm(cbind(y,n-y) ~ x, family = binomial(link = "logit")) #一般化線形モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

}

histn(p1, xlab = "P value", ylab = "頻度") #図6の描画

histn(p2, xlab = "P value", ylab = "頻度") #図7の描画

sum(p1 < 0.05)

## [1] 10000

sum(p2 < 0.05)

## [1] 461

------------------------------------------------------

すると切片の係数の検出力は100%であることがわかる。一方で、xの係数の危険率は5%程度に収まっている。このように、適切な解析が行われていることがわかるだろう。

　次はxの係数が存在している場合を考えてみる。

------------------------------------------------------

b0 <- -3

b1 <- 0.25

m <- 30

x <- runif(m, 0, 30)

n <- sample(20:30, size = m, replace = T)

y <- rbinom(n = m, size = n, prob = 1/(1+exp(-(b1 * x + b0))))

plotn(y/n ~ x)#図8の描画

------------------------------------------------------

今回はxが増加するとともに、yも増加している傾向が見て取れる。ではGLMで解析しよう。

------------------------------------------------------

fit1 <- glm(cbind(y, n-y) ~ x, family = binomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

## 

## Call:

## glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.7099  -0.6758  -0.2318   0.6175   1.6820  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept) -2.94880    0.23573  -12.51   <2e-16 ***

## x            0.24309    0.01758   13.83   <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: 342.955  on 29  degrees of freedom

## Residual deviance:  20.263  on 28  degrees of freedom

## AIC: 113.81

## 

## Number of Fisher Scoring iterations: 4

b0e <- coef(fit2)[1,1]

b1e <- coef(fit2)[2,1]

plotn(y/n ~ x)

xx <- seq(min(x), max(x), length = 200)

yy <- 1/(1+exp(-(b0e + b1e * xx)))

overdraw(points(xx, yy, type = "l")) #図9の描画

------------------------------------------------------

線形予測子の切片の推定値は-2.95、xの係数の推定値は0.24であり、これも初めの設定と大きく変わらないことがわかる。回帰曲線は線形回帰の時と異なり、ロジスティック関数となっている。

　では、上記の設定でデータを10000回生成して、検出力を確認してみる。

------------------------------------------------------

b0 <- -3

b1 <- 0.25

m <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(m, 0, 30)

  n <- sample(20:30, size = m, replace = T)

  y <- rbinom(n = m, size = n, prob = 1/(1+exp(-(b1 * x + b0))))

  fit1 <- glm(cbind(y,n-y) ~ x, family = binomial(link = "logit")) #一般化線形モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)[1,4])

  p2 <- c(p2, coef(fit2)[2,4])

}

histn(p1, xlab = "P value", ylab = "頻度") #図10の描画

histn(p2, xlab = "P value", ylab = "頻度") #図11の描画

sum(p1 < 0.05)

## [1] 10000

sum(p2 < 0.05)

## [1] 10000

------------------------------------------------------

すると切片の係数の検出力は100%、xの係数の検出力は100%である。xの係数の設定を変えたことで、xの係数が有意に検出されるようになった。

一般化線形モデルを使った解析

　では、今度はデータを与えられた時の解析の方針などを考えてゆこう。例えば、次のようなデータが与えられたとする。nが調べた総数で、yはそのうちいくつが結実していたか、というようなデータだとしよう。データの図示はまず最初にやろう。

------------------------------------------------------

n <- c(18, 10, 16, 15, 11, 17, 20, 14, 18, 16, 15, 11, 

       19, 17, 18, 19, 17, 19, 16, 15, 16, 14, 13, 12, 

       10, 14, 20, 14, 10, 20, 19, 15, 11, 17, 12, 11, 

       17, 11, 16, 13, 14, 13, 11, 17, 17, 14, 13, 13, 

       19, 12, 11, 17, 18, 13, 12, 18, 19, 15, 15, 11, 

       14, 20, 12, 19, 13, 12, 17, 14, 17, 20, 10, 20, 

       18, 15, 13, 18, 14, 16, 14, 17, 11, 11, 13, 11, 

       20, 20, 17, 14, 14, 20, 14, 17, 20, 14, 13, 17, 

       15, 15, 20, 13)

y <- c(2, 2, 2, 2, 4, 8, 2, 4, 5, 1, 5, 1, 6, 4, 8, 8, 

       5, 4, 3, 6, 5, 5, 3, 4, 2, 6, 4, 2, 5, 5, 6, 2, 

       1, 3, 2, 3, 5, 2, 3, 6, 7, 3, 3, 2, 4, 3, 2, 4, 

       7, 4, 3, 4, 7, 4, 5, 5, 4, 2, 5, 6, 3, 7, 4, 6, 

       5, 5, 3, 5, 5, 6, 2, 4, 2, 1, 2, 5, 1, 7, 3, 5, 

       2, 5, 4, 3, 0, 5, 6, 4, 6, 8, 3, 6, 5, 3, 4, 3, 

       5, 3, 4, 3)

histn(y/n)#図12の描画

------------------------------------------------------

データは0以上の離散値データであり、カウントデータの割合を解析しようとしている。平均的な割合は0.3程度であろうか。こういう時は二項分布を解析に使うとよい。一旦は練習なので、説明変数がないデータの解析をしてみることにする。これはデータが正規分布の時であれば、平均が有意に0から異なるかを調べる1標本のt検定に相当する。今回の場合は、切片の係数が0から有意に異なるかを考えているわけだから、切片推定値がlog((yの総和)/(n-yの総和))であることを考慮すれば、(yの総和)と(n-yの総和)が有意に異なるかを検定していると言える。GLMで解析してみると下記のようになる。

------------------------------------------------------

fit1 <- glm(cbind(y, n-y) ~ 1, family = binomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = cbind(y, n - y) ~ 1, family = binomial(link = "logit"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.5103  -0.6798   0.0569   0.5841   1.9629  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept) -1.01947    0.05811  -17.54   <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: 101.97  on 99  degrees of freedom

## Residual deviance: 101.97  on 99  degrees of freedom

## AIC: 387.3

## 

## Number of Fisher Scoring iterations: 4

------------------------------------------------------

すると切片の係数は有意に0と異なっていることがわかる。今回の場合は二項分布の平均は下記のように4.03であり、確かに(nの平均)/(1+exp(切片の係数))と一致している。

------------------------------------------------------

mean(y)

## [1] 4.03

b0e <- coef(fit2)[1,1]

mean(n)*1/(1+exp(-b0e))

## [1] 4.03

------------------------------------------------------

さてここで尤度のことを思い出そう。今回推定された係数は尤度を最大にするように選ばれているはずである。対数尤度を計算するために下記のように関数を定義する。

------------------------------------------------------

f <- function(y, n, b0t){

  sapply(b0t, function(tmp) sum(log(dbinom(x = y, size = n, prob = 1/(1 + exp(-tmp))))))

} #あるb0に対応する対数尤度を計算

------------------------------------------------------

今回得られた最尤推定値b0e = -1.01947から±1の範囲の連続値を生成し、それらに対応する対数尤度を計算して、図示してみる。

------------------------------------------------------

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, n, seq(b0min, b0max, length = 2000)) #対数尤度の計算

plotn(seq(b0min, b0max, length = 2000), logL, xlab = "b0", ylab = "logL", type = "l")

overdraw(abline(v = b0e),

         axis(side = 1, at = b0e, labels = "^b0")) #図13の図示

------------------------------------------------------

この図を見るとわかるように、確かに推定値b0eで最大の値を示していることがわかるだろう。さらに狭い範囲で図示すると次のようになる。

------------------------------------------------------

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, n, seq(b0min, b0max, length = 2000)) #対数尤度の計算

plotn(seq(b0min, b0max, length = 2000), logL, xlab = "b0", ylab = "logL", type = "l")

overdraw(abline(v = b0e),

         axis(side = 1, at = b0e, labels = "^b0")) #図14の図示

------------------------------------------------------

二項分布とロジットリンクを仮定した時の対数尤度関数の二次導関数は以下のように求まる。それが正規分布の確率密度分布の二次導関数と一致することを使って、b0の分散を求めることができる。

以上から、対数尤度関数の頂点付近は分散1/(nの総和)(1/exp (b0) + exp(b0) + 2)の正規分布で近似できることがわかる。それゆえ、b0の標準偏差はこの値の平方根である。実際、下記のように計算すれば上記のGLMの結果である0.05811と一致する。

------------------------------------------------------

sqrt((1/sum(n))*(1/exp(-b0e)+exp(-b0e)+2))

## [1] 0.0581089

------------------------------------------------------

また、前回同様にWald統計量は推定値を推定値の標準誤差で割ったものであり、Wald検定はWald統計量が標準正規分布に従うことを利用して行える。

------------------------------------------------------

b0e/sqrt((1/sum(n))*(1/exp(-b0e)+exp(-b0e)+2))

## [1] -17.54405

sign2 <- function(x){

  if(x < 0){

    return(0)

  } else {

    return(1)

  }

}

pnorm(b0e/sqrt((1/sum(n))*(1/exp(-b0e)+exp(-b0e)+2)), mean = 0, sd = 1, lower.tail = !sign2(b0e)) * 2

## [1] 6.605363e-69

coef(fit2)[1,4]

## [1] 6.605255e-69

------------------------------------------------------

　逸脱度残差の二項分布およびロジットリンクの場合の計算は下記のように行う。

------------------------------------------------------

d <- 2*(log(dbinom(x = y, size = n, prob = y/n)) - log(dbinom(x = y, size = n, prob = 1/(1+exp(-b0e)))))

quantile(sign(y/n - 1/(1+exp(-b0e)))*sqrt(d))

##          0%         25%         50%         75%        100% 

## -3.51034935 -0.67983641  0.05690119  0.58407342  1.96285087

------------------------------------------------------

残差逸脱度は逸脱度残差の総和なので

------------------------------------------------------

sum(d)

## [1] 101.9716

------------------------------------------------------

あるいは残差逸脱度は上で求めた対数尤度の計算関数を使って、下記のように計算もできる。

------------------------------------------------------

-2*(f(y, n, b0e) -  sum(log(dbinom(x = y, size = n, prob = y/n))))

## [1] 101.9716

------------------------------------------------------

さらにf(y, n, b0e)は最大の対数尤度であるが、この値を-2倍し、推定したパラメータの数×2だけ、値を足したものがAICである。

------------------------------------------------------

-2*f(y, n, b0e)+2*1

## [1] 387.2971

------------------------------------------------------

　続いて、説明変数xを含んだ時の解析を紹介する。例えば、次のようなデータである。図示も併せて行おう。

------------------------------------------------------

x <- c(23.8, 20.7, 29.1, 0.9, 1.1, 8.1, 12.1, 16.9, 20, 12.5, 25.1, 

       10, 0.8, 2.8, 23.9, 18.2, 27.9, 2.3, 3.3, 14.9, 26.5, 8.8, 

       14.8, 24.7, 3.8, 0.2, 12.9, 10.2, 11.6, 26.2, 1.6, 25.2, 3.9,

       25.6, 1.5, 21.3, 16.7, 6.7, 11.2, 9.2, 7, 8.6, 21.2, 23.2, 

       24.1, 24, 18.5, 23.6, 17.3, 4.2)

n <- c(14, 15, 19, 10, 14, 11, 16, 20, 14, 14, 13, 15, 15, 14, 12,

       20, 16, 12, 14, 18, 18, 14, 20, 17, 15, 17, 17, 10, 11, 16,

       14, 15, 13, 12, 16, 11, 12, 19, 10, 14, 10, 11, 17, 12, 12,

       16, 19, 17, 14, 19)

y <- c(3, 4, 0, 6, 13, 4, 3, 3, 3, 9, 0, 5, 8, 6, 2, 3, 0, 8, 10, 

       4, 1, 5, 3, 0, 12, 14, 4, 4, 4, 1, 8, 1, 10, 0, 11, 1, 3, 

       12, 5, 5, 7, 5, 1, 0, 1, 0, 2, 1, 3, 10)

plotn(y/n ~ x)#図15の図示

------------------------------------------------------

nに対するyの割合を解析したいと考えており、ともに離散値をとるので二項分布を真っ先に思い浮かべるべきだろう。では解析をしてみよう。回帰を行い、回帰曲線を描く。

------------------------------------------------------

fit1 <- glm(cbind(y,n-y) ~ x, family = binomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.8133  -0.8522  -0.1149   0.5389   2.5933  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  1.12326    0.16157   6.952 3.59e-12 ***

## x           -0.15543    0.01246 -12.475  < 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: 271.179  on 49  degrees of freedom

## Residual deviance:  51.929  on 48  degrees of freedom

## AIC: 172.65

## 

## Number of Fisher Scoring iterations: 4

b0e <- coef(fit2)[1,1]

b1e <- coef(fit2)[2,1]

plotn(y/n ~ x) #図16の図示

xx <- seq(min(x), max(x), length = 200)

yy <- 1/(1+exp(-(b0e + b1e * xx)))

overdraw(points(xx, yy, type = "l", col = "red"))

------------------------------------------------------

線形予測子の切片は1.12、xの係数は-0.15と推定され、回帰曲線はロジスティック関数の形状となった。

　では、今回も対数尤度の定義から、2変数の場合の対数尤度を求めて図示してみよう。推定値b0e = 1.12、b1e = -0.15周りの対数尤度は下記のように計算できる。

------------------------------------------------------

f <- function(x, y, n, b0t, b1t){

  mapply(function(tmp1, tmp2) {

    sum(log(dbinom(x = y, size = n, prob = 1/(1+exp(-tmp1-tmp2*x)))))

    }, b0t, b1t)

}#あるb0, b1に対応する対数尤度を計算

db0 <- 0.01

db1 <- 0.001

b0min <- b0e - db0

b0max <- b0e + db0

b1min <- b1e - db1

b1max <- b1e + db1

logL <- outer(seq(b0min, b0max, length = 200), 

              seq(b1min, b1max, length = 200), 

              f, x = x, y = y, n = n)

#対数尤度の計算

lev <- unique(c(seq(min(logL), max(logL)- 0.01, length = 10), 

                seq(max(logL)-0.01, max(logL)- 0.0001, length = 10), 

                seq(max(logL) - 0.0005, max(logL)-0.00001, length = 5)))

image(seq(b0min, b0max, length = 200), seq(b1min, b1max, length = 200), logL, 

      col = magma(400), xlab = "b0", ylab = "b1", las = 1)#図17の図示

contour(seq(b0min, b0max, length = 200), seq(b1min, b1max, length = 200), logL,

        levels = lev, add = T)#等高線

abline(v = b0e)

abline(h = b1e)

axis(side = 1, at = b0e, labels = "^b0")

axis(side = 2, at = b1e, labels = "^b1")

------------------------------------------------------

すると、確かに最尤推定値のb0eとb1eがその周りで最も高い尤度となっていることがわかる。以下のような図示の仕方もある。

------------------------------------------------------

persp(seq(b0min, b0max, length = 200), seq(b1min, b1max, length = 200), logL, 

      xlab = "b0", ylab = "b1", zlab = "logL",

      theta=120, phi=20, expand=0.5, ticktype="detailed")#図18の図示

------------------------------------------------------

　逸脱度残差は次のように計算する。

------------------------------------------------------

d <- 2*(log(dbinom(x = y, size = n, prob = y/n)) - log(dbinom(x = y, size = n, prob = 1/(1+exp(-b0e-b1e*x)))))

quantile(sign(y/n - 1/(1+exp(-b0e-b1e*x)))*sqrt(d))

##         0%        25%        50%        75%       100% 

## -1.8133328 -0.8521702 -0.1148924  0.5388882  2.5932706

------------------------------------------------------

また上記と同様にこのdの総和が残差逸脱度であり、対数尤度の計算関数を使っても同様の値を得られる。

------------------------------------------------------

sum(d)

## [1] 51.92879

-2*(f(x, y, n, b0e, b1e) -  sum(log(dbinom(x = y, size = n, prob = y/n))))

## [1] 51.92879

------------------------------------------------------

説明変数がないときp = (yの総和)/(nの総和)であったことを思い出せば、すぐにヌルモデルのb0を計算できて、以下のようにヌル逸脱度を計算できる。

------------------------------------------------------

dn <- 2*(log(dbinom(x = y, size = n, prob = y/n)) - log(dbinom(x = y, size = n, prob = sum(y)/sum(n))))

sum(dn)

## [1] 271.1789

------------------------------------------------------

f(x, y, n, b0e, b1e)は最大の対数尤度であり、AICの定義から-2×最大対数尤度+2×推定したパラメータの数を求めれば、

------------------------------------------------------

-2*f(x, y, n, b0e, b1e) + 2*2

## [1] 172.6504

------------------------------------------------------

以上のようにGLMを使うことでカウントデータに基づく割合データの解析を行うことができた。

追記1：ベルヌーイ分布的なGLM

　前述の通り、ベルヌーイ分布Bernoulli distributionとは二項分布においてn = 1のときの分布である。確率質量関数は下記の通りで、xは0か1をとり、x = 0のときf(x) = p、x = 1のときf(x) = 1 - pとなる。

今回の二項分布によるGLMは、線形回帰やそれ以外のGLMと異なり、被説明変数の指定方法がちょっと特殊であった。二項分布はサンプルサイズで分布の形が変わる。それゆえ、モデリングにはサンプルサイズnも被説明変数yとともに必要になるのである。では、二項分布をサンプルサイズに依存しない分布に置き換えれば、GLMの指定方法は今までと同じになるはずだ。この時に、ベルヌーイ分布が登場する。ベルヌーイ分布はn = 1の二項分布で、サンプルサイズに依存しない分布だ。

　二項分布の代わりにベルヌーイ分布を使うのであれば、当然、被説明変数も変換が必要である。二項分布ではyは0からnまで取れたが、ベルヌーイ分布では0か1しかとれない。二項分布におけるn個中y個のような値は、ベルヌーイ分布では1がy個、0がn - y個に相当する。例えば、二項分布では(y, n) = (3, 10)であるとき、その値に相当するベルヌーイ分布での値はy = (1, 1, 1, 0, 0, 0, 0, 0, 0, 0)である。また、二項分布が説明変数xに紐づいていることもあるだろう。例えば、二項分布において、(x, y, n) = (3.8, 6, 10)であれば、ベルヌーイ分布ではx = (3.8, 3.8, 3.8, ……, 3.8)と3.8が10個並ぶベクトルとy = (1,1,1,1,1,1,0,0,0,0)のベクトルとなる。

こう見ると、そもそも二項分布とは、0か1しかとらない被説明変数yをサンプルサイズでまとめている分布であると考えられる。

　さて、この考え方が正しいことを示すために、ベルヌーイ分布的な発想でGLMを行っても結果が一致することを示そう。上記と同じデータを使う。

------------------------------------------------------

x <- c(23.8, 20.7, 29.1, 0.9, 1.1, 8.1, 12.1, 16.9, 20, 12.5, 25.1, 

       10, 0.8, 2.8, 23.9, 18.2, 27.9, 2.3, 3.3, 14.9, 26.5, 8.8, 

       14.8, 24.7, 3.8, 0.2, 12.9, 10.2, 11.6, 26.2, 1.6, 25.2, 3.9,

       25.6, 1.5, 21.3, 16.7, 6.7, 11.2, 9.2, 7, 8.6, 21.2, 23.2, 

       24.1, 24, 18.5, 23.6, 17.3, 4.2)

n <- c(14, 15, 19, 10, 14, 11, 16, 20, 14, 14, 13, 15, 15, 14, 12,

       20, 16, 12, 14, 18, 18, 14, 20, 17, 15, 17, 17, 10, 11, 16,

       14, 15, 13, 12, 16, 11, 12, 19, 10, 14, 10, 11, 17, 12, 12,

       16, 19, 17, 14, 19)

y <- c(3, 4, 0, 6, 13, 4, 3, 3, 3, 9, 0, 5, 8, 6, 2, 3, 0, 8, 10, 

       4, 1, 5, 3, 0, 12, 14, 4, 4, 4, 1, 8, 1, 10, 0, 11, 1, 3, 

       12, 5, 5, 7, 5, 1, 0, 1, 0, 2, 1, 3, 10)

------------------------------------------------------

今このデータは二項分布的な形になっているので、ベルヌーイ分布的な形に展開する。

------------------------------------------------------

d2 <- data.frame(x = unlist(mapply(function(x, n){rep(x, n)}, x, n)),

                 y = unlist(mapply(function(y, n){c(rep(1, y), rep(0,n-y))}, y, n)))

head(d2, 30)

##       x y

## 1  23.8 1

## 2  23.8 1

## 3  23.8 1

## 4  23.8 0

## 5  23.8 0

## 6  23.8 0

## 7  23.8 0

## 8  23.8 0

## 9  23.8 0

## 10 23.8 0

## 11 23.8 0

## 12 23.8 0

## 13 23.8 0

## 14 23.8 0

## 15 20.7 1

## 16 20.7 1

## 17 20.7 1

## 18 20.7 1

## 19 20.7 0

## 20 20.7 0

## 21 20.7 0

## 22 20.7 0

## 23 20.7 0

## 24 20.7 0

## 25 20.7 0

## 26 20.7 0

## 27 20.7 0

## 28 20.7 0

## 29 20.7 0

## 30 29.1 0

------------------------------------------------------

このように、各xの値ごとに0と1が並んだデータフレームに変換した。x = 23.8を確認すると、yには0が11個、1が3個並んでいる。これは元データでx = 23.8のとき、y = 3、n = 14に対応している。さて、このデータを図示しよう。とはいえ、yは0か1しかとらないので、点がかぶりまくっている。なので、点の色の濃度を調整して、たくさんの値がかぶっているほど濃く見えるようにする。

------------------------------------------------------

plotn(y ~ x, data = d2, cex = 3, col.dot = "#00000010")#図19の図示

------------------------------------------------------

この図は本質的には図15と同じものである。xが大きくなるにつれて、y = 1の点が薄くなり、y = 0の点が濃くなるので、y = 0を示すものが増えていくことがわかるだろう。では、この形のデータをGLMで解析する。既にデータの形はベルヌーイ分布的に変換してあるので、glm関数へのインプットは工夫する必要はない。familyは二項分布でよい。ベルヌーイ分布も二項分布の一種だからである。

------------------------------------------------------

fit1 <- glm(y ~ x, data = d2,  family = binomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ x, family = binomial(link = "logit"), data = d2)

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.6622  -0.7280  -0.3773   0.8295   2.4675  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  1.12326    0.16157   6.952 3.59e-12 ***

## x           -0.15543    0.01246 -12.474  < 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: 909.57  on 733  degrees of freedom

## Residual deviance: 690.32  on 732  degrees of freedom

## AIC: 694.32

## 

## Number of Fisher Scoring iterations: 5

------------------------------------------------------

確かに推定値は二項分布で解析した場合と一致している。回帰曲線も描いておこう。

------------------------------------------------------

b0e <- coef(fit2)[1,1]

b1e <- coef(fit2)[2,1]

plotn(y ~ x, data = d2, cex = 3, col.dot = "#00000010")#図20の図示

xx <- seq(min(x), max(x), length = 200)

yy <- 1/(1+exp(-(b0e + b1e * xx)))

overdraw(points(xx, yy, type = "l", col = "red"))

------------------------------------------------------

ゆえに、上記のようなデータの変換は推定値には影響を及ぼさず、二項分布による回帰の実態は被説明変数が0/1をとるデータの解析に他ならない。しかし、データの図示などを見てもわかると思うが、二項分布的に解釈したほうがわかりやすいだろう。

　それと、自由度、逸脱度（対数尤度）、AICは、二項分布で解析した時と異なっている。今一度、考えると、二項分布において、(x, y, n) = (3.8, 6, 10)であれば、ベルヌーイ分布ではx = (3.8, 3.8, 3.8, ……, 3.8)と3.8が10個並ぶベクトルとy = (1,1,1,1,1,1,0,0,0,0)のベクトルとなるが、本来はx = 3.8における1つの観察値であったものが、10倍になっている。言い方を変えれば、ベルヌーイ分布の0/1値はx = 3.8における擬反復となっている。それゆえ、推定値に影響を与えないとはいえ、今回のデータは二項分布で考えるのが適切だ。もちろん、n個中y個のようなデータではなくて、1回の実験につき食べた/食べない、逃げた/逃げない、のような二値化できるデータなら、ベルヌーイ分布的に解釈をして問題はない。例えば、温度が上昇すると動きやすくなって捕食成功しやすくなる生物を考えよう。各温度につき、獲物と捕食者を1匹ずつ容器に入れて一定時間後に捕食したか否かを計測する。当然、捕食が必ずしも起こるとは限らない。こんな時、温度に対する捕食率 = 捕食成功回数/実験数をモデリングしたくなるが、本項最初の議論と同様、こういうときは1回の施行につき、成功(1)/失敗(0)をベルヌーイ分布を使ってモデリングするのがスマートだ。

追記2：ポアソン分布と二項分布の関係

　ポアソン分布と二項分布は、ともにカウントデータを記述できる確率分布ということでよく似ている。実際、二項分布にある仮定を課すことでポアソン分布を導出することが可能である。その仮定はnpは一定であるとし、nが極めて大きく、pが極めて小さいときに、ポアソン分布となる。そのことを確かめよう。

ここでnについて極限をとる。

この計算から何が言えるかというと、ポアソン分布は試行回数が∞のとき、ごく低頻度で起こる現象を記述しているということである。それゆえに、特にポアソン分布は0付近での解析に威力を発揮するといってよいだろう。