離散値/連続値型の割合データのモデリング

　二項分布のGLMでは、「n個のうちy個」というような、離散値/離散値型の割合データに基づくモデリングを行った。けれども、割合データは必ずしも離散値だけで構成されているとは限らない。例えば、時間当たりの鳥の観測数についてモデリングするとしよう。なるべく条件をそろえたいという願望から、類似した気象条件がそろう時間Aの間だけ、鳥の個体数yを観測することに努める。すると、気象条件は人間が操作できないから、日によっては観測できる時間が長かったり、短かったりする。その中で、観察時間中、何個体の鳥を観測できたかをモデリングするのだ。当然だが、観測時間が長ければ長いほど、観測できる鳥の数は増えるだろうから、観測時間を何かしらの形でモデリングしないと、誤った結論を導きかねない。そこで、単位時間当たりの観察数y/Aを解析したくなるものである。この時、分子は離散値だが、分母は連続値である。このようなデータに合う、確率分布を考えられるだろうか。実はこのようなデータは、ちょっとした工夫をすることでポアソン回帰に帰着をして考えることができる。それは、観察時間をオフセットとして扱う、という方法である。

オフセット項：係数を推定しない説明変数

　オフセットoffsetは、モデルの中で係数を推定しない説明変数である。もっと言えば、オフセット項の係数は1で固定している、ということもできるだろう。数式的にオフセットとは何かを下記で示す。線形予測子をβ・xとし、y/A ~ exp(β・x)というモデルにおいてβを推定する。

式変形を行うと、対数リンクのポアソン回帰をした時とほとんど同じ式の形になる。唯一違うのは、Yに対して、線形予測子にlogAという項が追加されていることだ。このlogAをオフセット項と呼ぶ。このモデル式は、yとの関係性を見たい説明変数xと規格化のための調整項としてオフセット項が存在しており、係数を推定するのはxの係数βだけである。オフセット項には推定するべき係数が存在しないため、係数は1で固定されている。ここで、yがポアソン分布に従うことを仮定して、GLMを行うのである。

オフセット項を含むポアソン分布を使ったGLM

　では、オフセットを考慮したGLMについて考察してゆこう。データの生成過程を図示して様子をつかむことにする。 上記のオフセット項の考察から、exp(β・x)に係数Aがついていることに注意しよう。

例えば、一つの説明変数xが存在する時、上記のような状況でデータが生成されていると考える。対数尤度は下記の通りになる。

さて、この対数尤度をパラメータβの関数と考えて、対数尤度を最大にするようなβを求めるのが最尤法だった。対数尤度が説明変数xを含まなければ下記のようにβ0を解析的に解くことができる。

すると、対数尤度を最大にするβ0は(Aの平均)×e^β0 = (yの平均)を満たすことがわかる。リンク関数に立ち返れば、(Aの平均)×e^β0 = λであるので、最尤推定値はポアソン分布の平均を(Aの平均)で割った値を推定していることになる。

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

　ではRを使って、GLMの検出力に関して検討しよう。以下のように、線形予測子の切片の係数は0ではないが、説明変数xの係数は0である場合をまず考える。オフセット項の効果として、log(A)が線形予測子に足されていることに注意する。

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

library(plotn)

library(viridis)

floor2 <- function(x, digits = 0){

  floor(x*10^digits)/(10^digits)

}

ceiling2 <- function(x, digits = 0){

  ceiling(x*10^digits)/(10^digits)

}

b0 <- 0.5

b1 <- 0

n <- 100

x <- runif(n, 0, 30)

A <- runif(n, 40, 60)

y <- rpois(n = n, lambda = exp(b1 * x + b0 + log(A)))

plotn(y ~ x) #図1の図示

plotn(y ~ A) #図2の図示

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

データの生成過程通り、xに対して特に関係のなさそうなデータが生成された。一方、yとAに関しては正の相関がありそうだ。では、このデータをGLMで解析する。ここではガンマ分布と対数リンクを指定するので、poisson(link="log")と指定する。glm関数内に、offset引数があり、そこにオフセットにあたるAを入力する。ただし、オフセットにはリンク関数で変形した値を入力する。つまり、ここではlog(A)である。回帰曲線も描いてみよう。回帰曲線は、ちょっと工夫をして、xを固定した時にAを変化させたときの図とAを固定した時にxを変化させたときの図を作成する。発想は、以前行った重回帰分析の時の図示の仕方と同じだ。

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

fit1 <- glm(y ~ x, family = poisson(link = "log"), offset = log(A)) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"), offset = log(A))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.48823  -0.63318   0.06046   0.62864   2.29207  

## 

## Coefficients:

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

## (Intercept) 4.892e-01  2.184e-02  22.398   <2e-16 ***

## x           5.771e-05  1.304e-03   0.044    0.965    

## ---

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

## 

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

## 

##     Null deviance: 91.846  on 99  degrees of freedom

## Residual deviance: 91.845  on 98  degrees of freedom

## AIC: 717.68

## 

## Number of Fisher Scoring iterations: 4

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

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

plotn(y ~ x) #図3の図示

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

m <- floor2(min(log(A)), digits = 1)

M <- ceiling2(max(log(A)), digits = 1)

for(i in seq(m, M, length = 5)){

  yy <- exp(b0e + b1e * xx + i)

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

}

plotn(y ~ A) #図4の図示

AA <- seq(min(A), max(A), length = 200)

for(i in c(min(x), mean(x), max(x))){

  yy <- exp(b0e + b1e * i + log(AA))

  overdraw(points(AA, yy, type = "l"))

}

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

GLMの結果をsummary関数を通して得られる結果では、線形予測子の切片の推定値は0.49、xの係数の推定値は0.000058であり、初めの設定と大きく変わらないことがわかる。オフセット項の係数は推定されていないことに注目しよう。図3ではyの大きなところの回帰曲線ほど、Aの値が大きいときの回帰曲線であることを示す。図4でも同様にyの大きなところの回帰曲線は、xの値が大きいのであるが、xの係数がほぼ0なので回帰曲線はほぼ重なっている。

　ところで、類似した解析としてオフセット項を考慮しない場合、と、オフセット項を係数を推定する説明変数に加える場合が、考えられる。以下にそれぞれの解析結果を示し、回帰曲線の様子を描く。まずは、オフセット項を考慮しない解析について。

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

fit3 <- glm(y ~ x, family = poisson(link = "log")) #一般化線形モデル

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.2062  -1.0907   0.0029   0.9283   3.3320  

## 

## Coefficients:

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

## (Intercept) 4.361620   0.021887 199.281   <2e-16 ***

## x           0.002074   0.001308   1.586    0.113    

## ---

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

## 

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

## 

##     Null deviance: 206.82  on 99  degrees of freedom

## Residual deviance: 204.31  on 98  degrees of freedom

## AIC: 830.14

## 

## Number of Fisher Scoring iterations: 4

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

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

plotn(y ~ x) #図5の図示

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

yy <- exp(b0e + b1e * xx)

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

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

上記のオフセットを考慮しない解析では、切片およびxの係数が大きく推定されている。xの係数は元の50倍程度大きいが、図としてはあまり影響が見られない（まあ、もともと0に近いのでこの変動はかなりあてにならない）。それ以上に、切片項の影響が大きく見える。というのも、オフセット項で考慮されるべき変動は、すべて切片項に押し付けられる形で推定されるからである。Y ~ β・x + log(A)というモデル式からわかるように、Y ~ (β0 + log(A)) + β1 x1 + β2 x2 +……という形に変形できる。つまり、log(A)を解析上無視することは、β0（切片項）を過大推定することに他ならない。

　とはいえ、考慮できなかった説明変数の存在で推定値がずれることはある意味よくあることであるから、注意したほうが良い、としか言えないかもしれない。それ以上に、ある意味で厄介な問題をはらむのが次の、オフセット項を係数を推定する説明変数に含む場合である。状況をそろえるため、log(A)の形で説明変数に組み込む

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

fit5 <- glm(y ~ x + log(A), family = poisson(link = "log")) #一般化線形モデル

fit6 <- summary(fit5)

fit6

## 

## Call:

## glm(formula = y ~ x + log(A), family = poisson(link = "log"))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.45639  -0.63724   0.02484   0.63257   2.36815  

## 

## Coefficients:

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

## (Intercept)  3.243e-01  3.823e-01   0.848    0.396    

## x           -2.094e-05  1.316e-03  -0.016    0.987    

## log(A)       1.042e+00  9.832e-02  10.603   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 206.823  on 99  degrees of freedom

## Residual deviance:  91.658  on 97  degrees of freedom

## AIC: 719.49

## 

## Number of Fisher Scoring iterations: 4

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

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

bAe <- coef(fit6)[3,1]

plotn(y ~ x) #図6の図示

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

m <- floor2(min(log(A)), digits = 1)

M <- ceiling2(max(log(A)), digits = 1)

for(i in seq(m, M, length = 5)){

  yy <- exp(b0e + b1e * xx + bAe*i)

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

}

plotn(y ~ A) #図7の図示

AA <- seq(min(A), max(A), length = 200)

for(i in c(min(x), mean(x), max(x))){

  yy <- exp(b0e + b1e * i + bAe*log(AA))

  overdraw(points(AA, yy, type = "l"))

}

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

オフセット項として組み込んでも、組み込まなくても、推定された係数および図は類似した結果が出ており、どっちでもいい気がしてしまう……が、いくつか点において明確にオフセット項として組み込むべき様相が出ている。

　一つ目は切片項の推定値が不安定になることである。これは標準誤差がオフセット項として組み込んだ時と比較した標準誤差の10倍以上になっていることからわかる。上記の、切片項の過大推定と同様、本来の切片項β0とlog(A)がともに切片項としての役割を担うため、多重共線性と類似の問題が生じ、不安定になる。

　二つ目は（平均的には）AICが無駄に高くなることである。AICの意義は深く紹介しないが、「AICが小さいほど予測性のいいモデルである」という指標であり、発想は罰則項に近い。AICは-2×(最大対数尤度) + 2×(推定したパラメータ数)と計算する。log(A)の係数を推定しなかった時（オフセット項として考慮）はパラメータは2つ、推定した場合はパラメータは3つとなる。さてこの時、log(A)の係数を推定しなくても、推定しても、最大対数尤度がほとんど変わらないとき、推定したパラメータ分だけAICは大きくなる。上記の例でもAICはほぼ2増加しており、これはlog(A)の係数を推定しても最大対数尤度の大きな改善はなかったことを示している。それもそのはずで、そもそもオフセット項の出自を考えれば、係数1から外れた値をとることはほとんどないはずである。それゆえに、1と仮定しないことは無駄に係数を推定していることになる。

　三つ目は二つ目と関連して、自由度を無駄に消費することである。ありがたみを感じにくい要素であるが、この発想は一般化線形混合モデルを考えるときにも関連してくる。例えば、あるモデルの残差の自由度は(サンプルサイズ) - (推定したパラメータ数)である。自由度が大きいほど、パラメータの推定精度が高まるわけだから、無駄なパラメータや興味のないパラメータを推定すればするほど、推定精度は低くなる。

　上記ではオフセット項は1つで済んでいるが、これがもっと多くなれば上記の問題はより顕著になるだろう。例えば、この世のデータは「単位時間あたり」だけでなく、「単位時間・面積あたり」であることもあるだろう。もっと割る数が増える可能性もある。

　では、上記の設定でデータを10000回生成して、検出力を確認してみる。上記のように、（1）オフセット項として考慮、（2）オフセット項を考慮しない、（3）オフセットのするべき項のパラメータ推定をする、という3パターンで解析する。P値だけでなく、係数の推定値も取り出しておこう。

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

b0 <- 0.5

b1 <- 0

n <- 30

p1_1 <- NULL

p1_2 <- NULL

p2_1 <- NULL

p2_2 <- NULL

p3_1 <- NULL

p3_2 <- NULL

p3_3 <- NULL

b1_1 <- NULL

b1_2 <- NULL

b2_1 <- NULL

b2_2 <- NULL

b3_1 <- NULL

b3_2 <- NULL

b3_3 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  A <- runif(n, 40, 60)

  y <- rpois(n = n, lambda = exp(b1 * x + b0 + log(A)))

  fit1 <- glm(y ~ x, family = poisson(link = "log"), offset = log(A)) #一般化線形モデル

  fit2 <- summary(fit1)

  fit3 <- glm(y ~ x, family = poisson(link = "log")) #一般化線形モデル

  fit4 <- summary(fit3)

  fit5 <- glm(y ~ x + log(A), family = poisson(link = "log")) #一般化線形モデル

  fit6 <- summary(fit5)

  p1_1 <- c(p1_1, coef(fit2)[1,4]) #オフセット考慮、切片項のP値

  p1_2 <- c(p1_2, coef(fit2)[2,4]) #オフセット考慮、x係数のP値

  p2_1 <- c(p2_1, coef(fit4)[1,4]) #オフセット考慮なし、切片項のP値

  p2_2 <- c(p2_2, coef(fit4)[2,4]) #オフセット考慮なし、x係数のP値

  p3_1 <- c(p3_1, coef(fit6)[1,4]) #オフセットの係数推定、切片項のP値

  p3_2 <- c(p3_2, coef(fit6)[2,4]) #オフセットの係数推定、x係数のP値

  p3_3 <- c(p3_3, coef(fit6)[3,4]) #オフセットの係数推定、log(A)係数のP値

  b1_1 <- c(b1_1, coef(fit2)[1,1]) #オフセット考慮、切片項

  b1_2 <- c(b1_2, coef(fit2)[2,1]) #オフセット考慮、x係数

  b2_1 <- c(b2_1, coef(fit4)[1,1]) #オフセット考慮なし、切片項

  b2_2 <- c(b2_2, coef(fit4)[2,1]) #オフセット考慮なし、x係数

  b3_1 <- c(b3_1, coef(fit6)[1,1]) #オフセットの係数推定、切片項

  b3_2 <- c(b3_2, coef(fit6)[2,1]) #オフセットの係数推定、x係数

  b3_3 <- c(b3_3, coef(fit6)[3,1]) #オフセットの係数推定、log(A)係数

}

histn(p1_1, xlab = "P value", ylab = "頻度") #図8の図示

try(histn(p2_1, xlab = "P value", ylab = "頻度"))

## Error in hist.default(x = xx, breaks = breaks, plot = F) : 

##   invalid number of 'breaks'

histn(p3_1, xlab = "P value", ylab = "頻度") #図9の図示

histn(p1_2, xlab = "P value", ylab = "頻度") #図10の図示

histn(p2_2, xlab = "P value", ylab = "頻度") #図11の図示

histn(p3_2, xlab = "P value", ylab = "頻度") #図12の図示

histn(p3_3, xlab = "P value", ylab = "頻度") #図13の図示

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

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

sum(p1_1 < 0.05) #オフセット考慮、切片項のP値

## [1] 10000

sum(p1_2 < 0.05) #オフセット考慮、x係数のP値

## [1] 499

sum(p2_1 < 0.05) #オフセット考慮なし、切片項のP値

## [1] 10000

sum(p2_2 < 0.05) #オフセット考慮なし、x係数のP値

## [1] 1742

sum(p3_1 < 0.05) #オフセットの係数推定、切片項のP値

## [1] 1074

sum(p3_2 < 0.05) #オフセットの係数推定、x係数のP値

## [1] 503

sum(p3_3 < 0.05) #オフセットの係数推定、log(A)係数のP値

## [1] 9995

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

シミュレーションをしてみると、データ生成過程から当然であるがオフセット項として考慮するモデルが最も正しい推定になっていることがわかる。オフセット項として考慮する場合、切片項は検出力100%、x係数の危険率5%程度である。オフセット項を考慮しない場合、切片項は検出力100%、x係数の危険率17%程度であり、xの係数の推定はあてにならないことがわかる。オフセット項の係数を推定する場合、切片項は検出力10%、x係数の危険率5%程度であり、こちらは切片項の推定があてにならない。続いて、具体的に推定された係数について考察しよう。

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

histn(b1_1, xlab = "estimate", ylab = "頻度") #図14の図示

histn(b2_1, xlab = "estimate", ylab = "頻度") #図15の図示

histn(b3_1, xlab = "estimate", ylab = "頻度") #図16の図示

histn(b1_2, xlab = "estimate", ylab = "頻度") #図17の図示

histn(b2_2, xlab = "estimate", ylab = "頻度") #図18の図示

histn(b3_2, xlab = "estimate", ylab = "頻度") #図19の図示

histn(b3_3, xlab = "estimate", ylab = "頻度") #図20の図示

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

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

mean(b1_1) #オフセット考慮、切片項の平均

## [1] 0.4995244

var(b1_1) #オフセット考慮、切片項の分散

## [1] 0.001739038

mean(b2_1) #オフセット考慮なし、切片項の平均

## [1] 4.411001

var(b2_1) #オフセット考慮なし、切片項の分散

## [1] 0.003567823

mean(b3_1) #オフセットの係数推定、切片項の平均

## [1] 0.4895007

var(b3_1) #オフセットの係数推定、切片項の分散

## [1] 0.5207213

mean(b1_2) #オフセット考慮、x係数の平均

## [1] 1.476596e-05

var(b1_2) #オフセット考慮、x係数の分散

## [1] 5.807758e-06

mean(b2_2) #オフセット考慮なし、x係数の平均

## [1] 2.075538e-05

var(b2_2) #オフセット考慮なし、x係数の分散

## [1] 1.206363e-05

mean(b3_2) #オフセットの係数推定、x係数の平均

## [1] 1.756517e-05

var(b3_2) #オフセットの係数推定、x係数の分散

## [1] 6.033597e-06

mean(b3_3) #オフセットの係数推定、log(A)係数の平均

## [1] 1.002483

var(b3_3) #オフセットの係数推定、log(A)係数の分散

## [1] 0.03382327

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

まず特に見比べるべきは、オフセット考慮とオフセット係数推定の場合の、切片項の平均値と分散である。これらの平均値は約0.5で一致するが、分散が大きく変わっており、オフセット係数推定すると切片項の推定が不安定になることがわかる。オフセット係数に関しては、やはり平均すれば1となり、その分散も小さいことから、わざわざ未知のパラメータとして推定する意義は薄い。

　では次は説明変数xの係数が0でない場合を考えよう。以下のようにデータを生成し、図示する。

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

b0 <- 0.5

b1 <- 0.03

n <- 100

x <- runif(n, 0, 30)

A <- runif(n, 40, 60)

y <- rpois(n = n, lambda = exp(b1 * x + b0 + log(A)))

plotn(y ~ x) #図21の図示

plotn(y ~ A) #図22の図示

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

データの生成過程通り、xとyは強く相関しそうだ。一方、yとAに関しては一見すればあまり関係がないようにも見える。このように、xとyの関係が強くなるとyとAの関係は弱く見えるようになる。では、このデータをGLMで解析する。回帰曲線も描いてみよう。

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

fit1 <- glm(y ~ x, family = poisson(link = "log"), offset = log(A)) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"), offset = log(A))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -3.03001  -0.62745  -0.06347   0.72691   3.08377  

## 

## Coefficients:

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

## (Intercept) 0.456701   0.019537   23.38   <2e-16 ***

## x           0.032898   0.001008   32.63   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 1238.2  on 99  degrees of freedom

## Residual deviance:  130.1  on 98  degrees of freedom

## AIC: 803.5

## 

## Number of Fisher Scoring iterations: 4

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

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

plotn(y ~ x) #図23の図示

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

m <- floor2(min(log(A)), digits = 1)

M <- ceiling2(max(log(A)), digits = 1)

for(i in seq(m, M, length = 5)){

  yy <- exp(b0e + b1e * xx + i)

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

}

plotn(y ~ A) #図24の図示

AA <- seq(min(A), max(A), length = 200)

for(i in c(min(x), mean(x), max(x))){

  yy <- exp(b0e + b1e * i + log(AA))

  overdraw(points(AA, yy, type = "l"))

}

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

GLMの結果をsummary関数を通して得られる結果では、線形予測子の切片の推定値は0.46、xの係数の推定値は0.033であり、初めの設定と大きく変わらないことがわかる。図23および24では、yの大きなところの回帰曲線ほど、それぞれAの値が大きいとき、xの値が大きいときの回帰曲線を示す。

　類似した解析としてオフセット項を考慮しない場合、と、オフセット項を係数を推定する説明変数に加える場合についても解析し、回帰曲線の様子を描く。まずは、オフセット項を考慮しない解析について。

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

fit3 <- glm(y ~ x, family = poisson(link = "log")) #一般化線形モデル

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.3507  -1.2139   0.3207   1.0358   3.5132  

## 

## Coefficients:

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

## (Intercept) 4.376625   0.019533  224.06   <2e-16 ***

## x           0.032851   0.001008   32.59   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 1351.96  on 99  degrees of freedom

## Residual deviance:  247.66  on 98  degrees of freedom

## AIC: 921.05

## 

## Number of Fisher Scoring iterations: 4

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

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

plotn(y ~ x) #図25の図示

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

yy <- exp(b0e + b1e * xx)

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

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

上記のオフセットを考慮しない解析では、切片の係数が大きく推定され、xの係数はそれほど変わらない。今回もAを考慮しない分、すべて切片項に押し付けられる形で推定されているため、切片項が大きくなっている。

　オフセット項の係数を推定すると同だろうか。

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

fit5 <- glm(y ~ x + log(A), family = poisson(link = "log")) #一般化線形モデル

fit6 <- summary(fit5)

fit6

## 

## Call:

## glm(formula = y ~ x + log(A), family = poisson(link = "log"))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.69204  -0.56573  -0.02565   0.76650   2.92985  

## 

## Coefficients:

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

## (Intercept) 1.127194   0.296095   3.807 0.000141 ***

## x           0.032880   0.001008  32.611  < 2e-16 ***

## log(A)      0.829233   0.075273  11.016  < 2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 1351.96  on 99  degrees of freedom

## Residual deviance:  124.97  on 97  degrees of freedom

## AIC: 800.36

## 

## Number of Fisher Scoring iterations: 4

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

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

bAe <- coef(fit6)[3,1]

plotn(y ~ x) #図26の図示

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

m <- floor2(min(log(A)), digits = 1)

M <- ceiling2(max(log(A)), digits = 1)

for(i in seq(m, M, length = 5)){

  yy <- exp(b0e + b1e * xx + bAe*i)

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

}

plotn(y ~ A) #図27の図示

AA <- seq(min(A), max(A), length = 200)

for(i in c(min(x), mean(x), max(x))){

  yy <- exp(b0e + b1e * i + bAe*log(AA))

  overdraw(points(AA, yy, type = "l"))

}

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

オフセット項として組み込んでも、組み込まなくても、推定された係数および図は類似した結果が出ているが、後述のようにシミュレートしてみると、オフセット項の整数推定が必ずしも必要ないことがわかる。

　では、上記の設定でデータを10000回生成して、検出力を確認してみる。上記のように、（1）オフセット項として考慮、（2）オフセット項を考慮しない、（3）オフセットのするべき項のパラメータ推定をする、という3パターンで解析する。P値だけでなく、係数の推定値も取り出しておこう。

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

b0 <- 0.5

b1 <- 0.03

n <- 30

p1_1 <- NULL

p1_2 <- NULL

p2_1 <- NULL

p2_2 <- NULL

p3_1 <- NULL

p3_2 <- NULL

p3_3 <- NULL

b1_1 <- NULL

b1_2 <- NULL

b2_1 <- NULL

b2_2 <- NULL

b3_1 <- NULL

b3_2 <- NULL

b3_3 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  A <- runif(n, 40, 60)

  y <- rpois(n = n, lambda = exp(b1 * x + b0 + log(A)))

  fit1 <- glm(y ~ x, family = poisson(link = "log"), offset = log(A)) #一般化線形モデル

  fit2 <- summary(fit1)

  fit3 <- glm(y ~ x, family = poisson(link = "log")) #一般化線形モデル

  fit4 <- summary(fit3)

  fit5 <- glm(y ~ x + log(A), family = poisson(link = "log")) #一般化線形モデル

  fit6 <- summary(fit5)

  p1_1 <- c(p1_1, coef(fit2)[1,4]) #オフセット考慮、切片項のP値

  p1_2 <- c(p1_2, coef(fit2)[2,4]) #オフセット考慮、x係数のP値

  p2_1 <- c(p2_1, coef(fit4)[1,4]) #オフセット考慮なし、切片項のP値

  p2_2 <- c(p2_2, coef(fit4)[2,4]) #オフセット考慮なし、x係数のP値

  p3_1 <- c(p3_1, coef(fit6)[1,4]) #オフセットの係数推定、切片項のP値

  p3_2 <- c(p3_2, coef(fit6)[2,4]) #オフセットの係数推定、x係数のP値

  p3_3 <- c(p3_3, coef(fit6)[3,4]) #オフセットの係数推定、log(A)係数のP値

  b1_1 <- c(b1_1, coef(fit2)[1,1]) #オフセット考慮、切片項

  b1_2 <- c(b1_2, coef(fit2)[2,1]) #オフセット考慮、x係数

  b2_1 <- c(b2_1, coef(fit4)[1,1]) #オフセット考慮なし、切片項

  b2_2 <- c(b2_2, coef(fit4)[2,1]) #オフセット考慮なし、x係数

  b3_1 <- c(b3_1, coef(fit6)[1,1]) #オフセットの係数推定、切片項

  b3_2 <- c(b3_2, coef(fit6)[2,1]) #オフセットの係数推定、x係数

  b3_3 <- c(b3_3, coef(fit6)[3,1]) #オフセットの係数推定、log(A)係数

}

histn(p1_1, xlab = "P value", ylab = "頻度") #図28の図示

try(histn(p2_1, xlab = "P value", ylab = "頻度"))

## Error in hist.default(x = xx, breaks = breaks, plot = F) : 

##   invalid number of 'breaks'

histn(p3_1, xlab = "P value", ylab = "頻度") #図29の図示

histn(p1_2, xlab = "P value", ylab = "頻度") #図30の図示

histn(p2_2, xlab = "P value", ylab = "頻度") #図31の図示

histn(p3_2, xlab = "P value", ylab = "頻度") #図32の図示

histn(p3_3, xlab = "P value", ylab = "頻度") #図33の図示

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

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

sum(p1_1 < 0.05) #オフセット考慮、切片項のP値

## [1] 10000

sum(p1_2 < 0.05) #オフセット考慮、x係数のP値

## [1] 10000

sum(p2_1 < 0.05) #オフセット考慮なし、切片項のP値

## [1] 10000

sum(p2_2 < 0.05) #オフセット考慮なし、x係数のP値

## [1] 10000

sum(p3_1 < 0.05) #オフセットの係数推定、切片項のP値

## [1] 1456

sum(p3_2 < 0.05) #オフセットの係数推定、x係数のP値

## [1] 10000

sum(p3_3 < 0.05) #オフセットの係数推定、log(A)係数のP値

## [1] 10000

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

シミュレーションをしてみると、データ生成過程から当然であるがオフセット項として考慮するモデルが最も正しい推定になっていることがわかる。オフセット項の係数を推定する場合、切片項は検出力14%程度であり、切片項の推定があてにならない。続いて、具体的に推定された係数について考察しよう。

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

histn(b1_1, xlab = "estimate", ylab = "頻度") #図34の図示

histn(b2_1, xlab = "estimate", ylab = "頻度") #図35の図示

histn(b3_1, xlab = "estimate", ylab = "頻度") #図36の図示

histn(b1_2, xlab = "estimate", ylab = "頻度") #図37の図示

histn(b2_2, xlab = "estimate", ylab = "頻度") #図38の図示

histn(b3_2, xlab = "estimate", ylab = "頻度") #図39の図示

histn(b3_3, xlab = "estimate", ylab = "頻度") #図40の図示

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

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

mean(b1_1) #オフセット考慮、切片項の平均

## [1] 0.5001033

var(b1_1) #オフセット考慮、切片項の分散

## [1] 0.001358293

mean(b2_1) #オフセット考慮なし、切片項の平均

## [1] 4.412519

var(b2_1) #オフセット考慮なし、切片項の分散

## [1] 0.003392623

mean(b3_1) #オフセットの係数推定、切片項の平均

## [1] 0.4977224

var(b3_1) #オフセットの係数推定、切片項の分散

## [1] 0.3252525

mean(b1_2) #オフセット考慮、x係数の平均

## [1] 0.0299807

var(b1_2) #オフセット考慮、x係数の分散

## [1] 3.757601e-06

mean(b2_2) #オフセット考慮なし、x係数の平均

## [1] 0.0299447

var(b2_2) #オフセット考慮なし、x係数の分散

## [1] 1.059196e-05

mean(b3_2) #オフセットの係数推定、x係数の平均

## [1] 0.02997579

var(b3_2) #オフセットの係数推定、x係数の分散

## [1] 3.927864e-06

mean(b3_3) #オフセットの係数推定、log(A)係数の平均

## [1] 1.000607

var(b3_3) #オフセットの係数推定、log(A)係数の分散

## [1] 0.02111381

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

まず特に見比べるべきは、オフセット考慮とオフセット係数推定の場合の、切片項の平均値と分散である。これらの平均値は約0.5で一致するが、分散が大きく変わっており、オフセット係数推定すると切片項の推定が不安定になることがわかる。オフセット係数に関しては、やはり平均すれば1となり、その分散も小さいことから、わざわざ未知のパラメータとして推定する意義は薄い。

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

　では、データを与えられた時の解析の方針などを考えてゆこう。例えば、次のようなデータが与えられたとする。今回はy/Aを解析したいと考えている。データの図示はまず最初にやろう。

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

A <- c(50.7, 51.8, 43.1, 51.1, 49.7, 55.2, 50.8, 50.4, 53, 

       47.2, 46.9, 41.7, 46.6, 59.2, 55.1, 48.7, 48.4, 53.1,

       46.5, 51.2, 51.1, 43.2, 58.2, 58.8, 54, 56, 51.5,

       47.2, 54, 54.1, 55.4, 45.7, 58.9, 50.2, 58.2, 52.3, 

       43.6, 40.1, 53.8, 47.2, 48.8, 51.2, 48.7, 40.1, 43.5,

       50.6, 54.1, 55.3, 45.7, 43.9, 46.6, 51.3, 42.2, 40.4,

       45.3, 47.5, 42.8, 53.5, 47.3, 41, 47.4, 51.5, 41.5, 

       42.9, 41.7, 49.1, 45.1, 41.8, 50.7, 50.4, 45.6, 59.1, 

       54.7, 55.8, 42.5, 44.6, 50.6, 57.3, 42, 50.2, 47.5, 

       49.2, 59, 45.6, 40.1, 48.9, 48, 59.8, 46.5, 54.8, 42.2, 

       52.5, 47, 56.2, 58.8, 53.5, 56.7, 46.8, 58.3, 41)

y <- c(16, 27, 18, 25, 20, 23, 9, 24, 21, 19, 21, 15, 18, 18,

       23, 24, 15, 20, 21, 20, 24, 13, 20, 23, 17, 22, 17, 13,

       19, 10, 24, 12, 19, 15, 24, 17, 12, 10, 23, 17, 24, 9,

       22, 16, 10, 9, 27, 19, 11, 25, 16, 18, 15, 12, 19, 28,

       12, 19, 17, 22, 17, 17, 13, 11, 10, 20, 12, 19, 28, 20,

       26, 19, 21, 14, 7, 15, 26, 25, 10, 12, 17, 21, 30, 14,

       14, 14, 19, 24, 12, 22, 14, 25, 27, 13, 25, 18, 14, 14, 18, 22)

plotn(y ~ A) #図41の図示

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

注意点として、あくまでyとAとの関係にはあまり興味がないということである。Aの効果を考慮（あるいは除去？）した時に、yを説明する切片項は有意に0から異なっているか、に興味があるのである。図41を見ればわかるようにyはAが大きくなるにつれて大きな値を示す傾向がある。GLMで解析してみると下記のようになる。

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

fit1 <- glm(y ~ 1, family = poisson(link = "log"), offset = log(A)) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ 1, family = poisson(link = "log"), offset = log(A))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.51635  -0.76789  -0.09994   0.68444   2.32381  

## 

## Coefficients:

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

## (Intercept) -1.00232    0.02346  -42.73   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 129.01  on 99  degrees of freedom

## Residual deviance: 129.01  on 99  degrees of freedom

## AIC: 601.15

## 

## Number of Fisher Scoring iterations: 4

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

すると切片の係数は有意に0と異なっていることがわかる。今回の場合、exp(切片の係数)は(yの合計)/(Aの合計)と一致する。見方を変えれば、(yの平均) = exp(切片の係数)×(Aの平均)でもある。

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

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

exp(b0e)

## [1] 0.3670262

sum(y)/sum(A)

## [1] 0.3670262

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

続いて、対数尤度を計算するように、下記のように関数を定義する。

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

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

  sapply(b0t, function(tmp) sum(log(dpois(x = y, lambda = A*exp(tmp)))))

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

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

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

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

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, A, 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")) #図42の図示

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

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

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

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, A, 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")) #図43の図示

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

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

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

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

1/sqrt(exp(b0e)*sum(A))

## [1] 0.0234597

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

逸脱度残差Deviance residualsは下記の通り。

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

d <- 2*(log(dpois(x = y, lambda = y)) - log(dpois(x = y, lambda = A*exp(b0e))))

quantile(sign(y-A*exp(b0e))*sqrt(d))

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

## -2.51634884 -0.76789195 -0.09993724  0.68443544  2.32381391

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

ここで求めたdの総和は残差逸脱度Residual devianceで、残差逸脱度はヌル逸脱度Null devianceと一致する。

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

sum(d)

## [1] 129.0091

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

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

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

-2*(f(y, A, b0e) -  sum(log(dpois(y, lambda = y))))

## [1] 129.0091

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

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

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

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

## [1] 601.1525

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

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

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

x <- c(13, -8.7, 1.1, 15.7, -8.6, -2.4, 7.5, -8.5, 15.1, -5.9, 

       8.9, 11, 2.6, 13.5, -5.5, -7, 8.2, 5.4, 12.3, -3.8, 3.8, 

       -6.1, -7.7, 1.9, 16.4, -0.4, -8.7, -5.3, 18.4, 9.3, 12.9, 

       17.5, 17, 7.5, 17.2, 15.3, -5.2, -2.2, 5.1, -0.9, -6.2, 

       7, -7.4, 9.8, -4, 3.7, -2.4, 13.2, -2.3, 18.2)

A <- c(59.9, 30, 33.1, 49.3, 20.6, 69.1, 26.3, 77.2, 64.7, 39.3,

       73, 28.5, 66.3, 28, 35.9, 37.5, 28.4, 72.3, 79.9, 21.8,

       53.2, 40.7, 20.5, 43.8, 68.6, 29.2, 39.9, 75.7, 32.2, 

       47.7, 55.1, 20.6, 47.5, 78.2, 55.5, 61.2, 79.6, 47.6, 

       74.1, 60, 23, 20.5, 54.1, 78.2, 61.4, 58, 53.5, 76.5, 30.1, 67.5)

y <- c(21, 76, 27, 8, 47, 96, 15, 174, 17, 71, 36, 11, 51, 8, 64,

       83, 8, 43, 23, 34, 40, 72, 40, 40, 7, 28, 90, 127, 3, 26,

       12, 2, 10, 58, 6, 13, 129, 61, 56, 47, 44, 7, 109, 29, 

       83, 29, 62, 24, 46, 12)

plotn(y ~ x)#図44の図示

plotn(y ~ A)#図45の図示

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

一見、図44を見る限り、通常のポアソン回帰でも問題はなさそうだが、図45を見ると弱いながらもyとAにも正の相関が見て取れる。前述の問題同様、yとAの関係にはあまり興味がない。Aによるyの変動を考慮しつつ、yとxの関係をオフセット項を利用した回帰で明らかにしよう。併せて、回帰曲線を描く。

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

fit1 <- glm(y ~ x, family = poisson(link = "log"), offset = log(A)) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"), offset = log(A))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.44027  -0.53818  -0.06146   0.50649   3.02609  

## 

## Coefficients:

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

## (Intercept)  0.004386   0.021859   0.201    0.841    

## x           -0.097161   0.003019 -32.182   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 1339.206  on 49  degrees of freedom

## Residual deviance:   47.568  on 48  degrees of freedom

## AIC: 313.42

## 

## Number of Fisher Scoring iterations: 4

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

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

plotn(y ~ x)#図46の図示

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

m <- floor2(min(log(A)), digits = 1)

M <- ceiling2(max(log(A)), digits = 1)

j <- 1

for(i in seq(m, M, length = 5)){

  yy <- exp(b0e + b1e * xx + i)

  overdraw(points(xx, yy, type = "l", col = .default_col[j]))

  j <- j+1

}

plotn(y ~ A)#図47の図示

AA <- seq(min(A), max(A), length = 200)

j <- 1

for(i in c(min(x), mean(x), max(x))){

  yy <- exp(b0e + b1e * i + log(AA))

  overdraw(points(AA, yy, type = "l", col = .default_col[j]))

  j <- j+1

}

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

線形予測子の切片は-1.16、xの係数は0.10と推定され、回帰曲線は指数関数の形状となった。

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

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

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

  mapply(function(tmp1, tmp2) {

    sum(log(dpois(x = y, lambda = A*exp(tmp1 + tmp2*x))))

    }, b0t, b1t)

}

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, A = A)

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)#図48の図示

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")#図49の図示

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

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

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

d <- 2*(log(dpois(x = y, lambda = y)) - log(dpois(x = y, lambda = A*exp(b0e + b1e*x))))

quantile(sign(y-A*exp(b0e+b1e*x))*sqrt(d))

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

## -2.44026787 -0.53817568 -0.06145837  0.50648997  3.02608530

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

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

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

sum(d)

## [1] 47.56835

-2*(f(x, y, A, b0e, b1e) - sum(log(dpois(y, lambda = y))))

## [1] 47.56835

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

ヌル逸脱度を計算するためには改めて、説明変数がない場合のGLMを行って、係数を求めて計算をする方法もある。しかし、説明変数がないときE(y) = exp(b0)×E(A)であり、A exp(b0) = λだったことを思い出せば、すぐにヌルモデルのb0を計算できて、以下のようにヌル逸脱度を計算できる。

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

dn <- 2*(log(dpois(x = y, lambda = y)) - log(dpois(x = y, lambda = A*mean(y)/mean(A))))

sum(dn)

## [1] 1339.206

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

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

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

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

## [1] 313.4213

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

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

追記：二項分布との比較

　二項分布のGLMでも割り算値を取り扱うことができたが、オフセット項との推定の違いはあるだろうか。二項分布に従うように生成した、以下のデータを解析することで考察しよう。

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

x <- c(15.2, 11.9, 21.4, 3.8, 21.5, 9.7, 8.2, 12.5, 4.2, 13.3, 12.1,

       25.3, 9.2, 4.4, 25, 22.4, 10.8, 26.6, 3.2, 5.6, 5, 6.7, 29.4,

       3.9, 7.3, 23.9, 26.5, 1.2, 0.9, 13.6, 17.8, 23.8, 12.1, 19.1,

       26, 14.2, 29.4, 15.1, 9.4, 12.7, 20.1, 5.6, 10, 10.7, 13.9, 

       17.7, 25.9, 26.4, 7.9, 17.6, 27, 29.6, 7.6, 5.5, 21.8, 20.8,

       2.8, 21.9, 20, 13.5, 13.2, 1.8, 2.5, 18.1, 19.6, 5.3, 11.3, 

       3.6, 20.6, 16.7, 29.2, 17.5, 24.9, 3.8, 24.6, 3.9, 26.2, 9.8,

       9.1, 27.5, 18.2, 18, 6.7, 19.5, 19.8, 26.3, 4.7, 12.2, 1.7, 

       9.4, 1.2, 3.8, 17.2, 2.2, 1.4, 6, 12.9, 27.3, 14.8, 19.9)

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

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

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

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

       12, 19, 17, 14, 20, 19, 14, 11, 13, 19, 19, 20, 11, 18, 14, 

       17, 16, 19, 13, 12, 10, 11, 11, 17, 20, 14, 14, 14, 10, 18, 

       20, 19, 11, 17, 12, 19, 13, 20, 11, 18)

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

       4, 4, 0, 2, 0, 0, 0, 4, 9, 3, 1, 0, 4, 0, 0, 3, 0, 1, 4, 0,

       0, 1, 0, 4, 2, 0, 0, 0, 1, 1, 0, 0, 1, 3, 0, 0, 6, 0, 0, 0, 

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

       0, 0, 2, 0, 1, 0, 5, 0, 6, 2, 11, 7, 0, 7, 7, 5, 0, 0, 1, 0)

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

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

まず二項分布を仮定した解析は下記の通り。

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

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  

## -2.1554  -0.6768  -0.2790   0.4429   2.5762  

## 

## Coefficients:

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

## (Intercept)  0.28868    0.14489   1.992   0.0463 *  

## x           -0.23806    0.01895 -12.560   <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: 404.343  on 99  degrees of freedom

## Residual deviance:  92.404  on 98  degrees of freedom

## AIC: 231.06

## 

## Number of Fisher Scoring iterations: 5

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

nをオフセットとし、ポアソン分布を仮定した解析は下記の通り。

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

fit3 <- glm(y ~ x, family = poisson(link = "log"), offset = log(n)) #一般化線形モデル

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"), offset = log(n))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.8537  -0.7428  -0.3483   0.3825   2.4110  

## 

## Coefficients:

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

## (Intercept) -0.34153    0.11124   -3.07  0.00214 ** 

## x           -0.19450    0.01599  -12.16  < 2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 343.997  on 99  degrees of freedom

## Residual deviance:  78.086  on 98  degrees of freedom

## AIC: 233.09

## 

## Number of Fisher Scoring iterations: 5

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

では、これらのデータをもとに回帰曲線を作図しよう。

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

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

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

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

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

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

overdraw(points(xx, yy1, type = "l", col = "red"))#二項分布仮定の回帰曲線

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

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

yy2 <- exp(b0e + b1e * xx)

overdraw(points(xx, yy2, type = "l", col = "blue"))#オフセット仮定の回帰曲線

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

当然、正しいモデルは二項分布仮定のほうだが、これを見る限りどちらでも大きな差はない（ように見える）。このように、一見うまくいっているように見えるのは、y/nが下限付近しか存在していないためである。二項分布で推定されるロジスティック関数を微分するとわかるが、ロジスティック関数の下限周りの微分係数は、ほとんど指数関数の微分係数と同じになるので、類似した見た目の推定値が算出されるのである。

　次のようなデータの時に問題は顕著に表れる。

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

x <- c(9.8, 14.7, 20.4, 8.3, 11.4, 19.9, 17.1, 15.7, 19, 4.9, 28.8,

       23.6, 20.8, 26.2, 16, 20, 9, 12.3, 0.6, 26.9, 19.7, 26.2, 7.8,

       24, 9.5, 21.7, 18.4, 22.7, 10.2, 11.9, 2.9, 10, 10.2, 22.1, 

       8.1, 7.4, 1.6, 18.1, 12.6, 12.3, 23.4, 28.5, 15.2, 13.2, 22.3, 

       28.3, 7.6, 1.7, 22.9, 13.5, 25.7, 27.6, 7, 1.7, 18.7, 25.8, 3.6,

       11.6, 23.4, 1.1, 1.4, 22.5, 24.2, 2.8, 7.9, 5.1, 29.6, 12.3, 

       29, 8.4, 27.3, 23.4, 4.1, 10.8, 7.9, 19.4, 14, 3.6, 10.7, 7.4, 

       20, 21, 6, 20.4, 22, 25.9, 4.1, 29.1, 27.8, 3.7, 21, 2.4, 16,

       5.7, 22.4, 15.1, 11.6, 25.2, 24, 20)

n <- c(11, 17, 12, 11, 17, 10, 13, 15, 14, 17, 12, 14, 14, 11, 11, 

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

       19, 12, 11, 18, 18, 12, 17, 13, 11, 13, 11, 16, 20, 18, 10, 

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

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

       13, 10, 18, 16, 14, 15, 15, 12, 20, 20, 19, 18, 11, 20, 12, 

       15, 15, 18, 16, 12, 14, 12, 14, 20, 18)

y <- c(8, 8, 1, 9, 11, 4, 6, 9, 2, 16, 1, 1, 1, 0, 3, 2, 12, 5, 20, 

       0, 4, 0, 10, 0, 15, 1, 0, 1, 9, 9, 19, 8, 8, 2, 16, 12, 16, 

       3, 7, 7, 0, 0, 7, 11, 3, 0, 10, 15, 0, 11, 0, 0, 16, 12, 1, 

       0, 19, 10, 0, 13, 20, 0, 0, 11, 16, 16, 0, 7, 0, 15, 0, 1, 17, 

       18, 13, 1, 3, 17, 8, 12, 1, 1, 12, 1, 1, 0, 18, 0, 0, 11, 1,

       15, 3, 15, 0, 4, 10, 1, 2, 3)

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

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

このデータの特徴は上限があることだ。オフセット項を含みポアソン分布を仮定するGLMは、あくまでポアソン回帰であるので、割合でありながらも割合は無限大に大きくなることを仮定している。とりあえず、解析しよう。二項分布仮定では下記の通り。

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

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  

## -2.8255  -0.7580  -0.3068   0.6440   2.3517  

## 

## Coefficients:

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

## (Intercept)  4.65341    0.23897   19.47   <2e-16 ***

## x           -0.33203    0.01585  -20.95   <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: 1234.9  on 99  degrees of freedom

## Residual deviance:  106.8  on 98  degrees of freedom

## AIC: 270.91

## 

## Number of Fisher Scoring iterations: 5

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

nをオフセットとし、ポアソン分布を仮定した解析は下記の通り。

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

fit3 <- glm(y ~ x, family = poisson(link = "log"), offset = log(n)) #一般化線形モデル

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"), offset = log(n))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.1488  -1.2778  -0.3515   0.8813   2.5659  

## 

## Coefficients:

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

## (Intercept)  0.501490   0.059767   8.391   <2e-16 ***

## x           -0.113700   0.005709 -19.914   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 660.48  on 99  degrees of freedom

## Residual deviance: 150.52  on 98  degrees of freedom

## AIC: 438.01

## 

## Number of Fisher Scoring iterations: 5

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

では、これらのデータをもとに回帰曲線を作図しよう。

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

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

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

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

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

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

overdraw(points(xx, yy1, type = "l", col = "red"))#二項分布仮定の回帰曲線

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

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

yy2 <- exp(b0e + b1e * xx)

overdraw(points(xx, yy2, type = "l", col = "blue"))#オフセット仮定の回帰曲線

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

オフセット項を仮定する解析は、全体として大外れだが、特にy/nが上限に近づいたあたりの推定値は大きく外れてしまう。このように、オフセット項＆ポアソン分布の解析は、必ずしも二項分布の解析の代替とはならないことを覚えておこう。

　二項分布っぽいデータに関して、ポアソン分布を使った解析法として、上記とは異なる方法を久保先生が紹介されていた。こちらも取り扱って考察してみよう。下記の式に基づけば、二項回帰で推定される係数βbはポアソン回帰で各カテゴリA、Bについて個別に求めたβAおよびβBの差に相当すると考えられる（この手法を"Poisson Trick"と呼んでいる文献を見かけた）。

そこで下記のようにポアソン回帰して、上記と一致するか考える。

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

fit5 <- glm(y ~ x, family = poisson(link = "log")) #一般化線形モデル

fit6 <- summary(fit5)

fit6

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.3886  -1.2657  -0.2348   0.4402   3.3875  

## 

## Coefficients:

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

## (Intercept)  3.279790   0.060788   53.95   <2e-16 ***

## x           -0.121279   0.005876  -20.64   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 700.10  on 99  degrees of freedom

## Residual deviance: 148.25  on 98  degrees of freedom

## AIC: 435.74

## 

## Number of Fisher Scoring iterations: 5

fit7 <- glm((n - y) ~ x, family = poisson(link = "log")) #一般化線形モデル

fit8 <- summary(fit7)

fit8

## 

## Call:

## glm(formula = (n - y) ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.6094  -1.1400  -0.3152   0.7383   2.9719  

## 

## Coefficients:

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

## (Intercept) 0.283651   0.115490   2.456    0.014 *  

## x           0.096565   0.005144  18.772   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 598.82  on 99  degrees of freedom

## Residual deviance: 169.78  on 98  degrees of freedom

## AIC: 501.82

## 

## Number of Fisher Scoring iterations: 5

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

以下のように作図する。

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

plotn(y/n ~ x)#図54の作図

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

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

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

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

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

b0A <- coef(fit6)[1,1]

b1A <- coef(fit6)[2,1]

b0B <- coef(fit8)[1,1]

b1B <- coef(fit8)[2,1]

b0e2 <- b0A - b0B

b1e2 <- b1A - b1B

yy2 <- 1/(1+exp(-(b0e2 + b1e2 * xx)))

overdraw(points(xx, yy2, type = "l", col = "blue"))

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

傾向は類似するが、一致しない。おそらくこれは、データの生成を二項分布から行ったことが原因で、2つのカテゴリが独立でないためであると考えられる。2つのカテゴリを独立にポアソン分布から生成すると、よく一致するようになる。

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

b0A <- 3

b1A <- -0.15

b0B <- -4

b1B <- 0.2

n <- 100

x <- runif(n, 0, 30)

x <- round(x, digits = 1)

yA <- rpois(n = n, lambda = exp(b0A+b1A*x))

yB <- rpois(n = n, lambda = exp(b0B+b1B*x))

n <- yA+yB

plotn(yA ~ x)#図55の作図

plotn(yB ~ x)#図56の作図

plotn(yA/n ~ x)#図57の作図

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

以下のように解析し、回帰曲線を描く。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.9343  -0.5740  -0.2289   0.7565   2.4104  

## 

## Coefficients:

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

## (Intercept)  2.90717    0.33558   8.663   <2e-16 ***

## x           -0.29824    0.02977 -10.019   <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: 367.002  on 94  degrees of freedom

## Residual deviance:  85.642  on 93  degrees of freedom

## AIC: 125.93

## 

## Number of Fisher Scoring iterations: 6

fit5 <- glm(yA ~ x, family = poisson(link = "log")) #一般化線形モデル

fit6 <- summary(fit5)

fit6

## 

## Call:

## glm(formula = yA ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.1635  -0.5224  -0.2845   0.2078   2.2716  

## 

## Coefficients:

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

## (Intercept)  2.10719    0.13630   15.46   <2e-16 ***

## x           -0.22190    0.02166  -10.25   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 281.078  on 99  degrees of freedom

## Residual deviance:  71.265  on 98  degrees of freedom

## AIC: 194.25

## 

## Number of Fisher Scoring iterations: 5

fit7 <- glm(yB ~ x, family = poisson(link = "log")) #一般化線形モデル

fit8 <- summary(fit7)

fit8

## 

## Call:

## glm(formula = yB ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.3396  -1.1090  -0.2628   0.5132   3.5204  

## 

## Coefficients:

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

## (Intercept) -1.043215   0.210459  -4.957 7.16e-07 ***

## x            0.104595   0.009329  11.211  < 2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 288.45  on 99  degrees of freedom

## Residual deviance: 132.54  on 98  degrees of freedom

## AIC: 339.49

## 

## Number of Fisher Scoring iterations: 5

plotn(yA/n ~ x)#図58の作図

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

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

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

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

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

b0A <- coef(fit6)[1,1]

b1A <- coef(fit6)[2,1]

b0B <- coef(fit8)[1,1]

b1B <- coef(fit8)[2,1]

b0e2 <- b0A - b0B

b1e2 <- b1A - b1B

yy2 <- 1/(1+exp(-(b0e2 + b1e2 * xx)))

overdraw(points(xx, yy2, type = "l", col = "blue"))

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

それでも完全に一致するわけではない。また、二項分布でも、下記のような、カテゴリカル解析のような場合は一致する。

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

b0 <- 1

b1 <- 2

m <- 50

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

y1 <- rbinom(n = m, size = n1, prob = 1/(1+exp(-b0)))

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

y2 <- rbinom(n = m, size = n2, prob = 1/(1+exp(-b0-b1)))

d <- data.frame(group = c(rep("A",m), rep("B",m)),

                n = c(n1, n2), y = c(y1, y2))

fit1 <- glm(cbind(y, n-y) ~ group, data = d, family = binomial(link = "logit"))

fit2 <- summary(fit1)

fit2

## 

## Call:

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

##     data = d)

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.50327  -0.65424   0.04388   1.25831   2.82233  

## 

## Coefficients:

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

## (Intercept)  1.02118    0.06325   16.14   <2e-16 ***

## groupB       2.13290    0.15635   13.64   <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: 386.42  on 99  degrees of freedom

## Residual deviance: 118.23  on 98  degrees of freedom

## AIC: 363.81

## 

## Number of Fisher Scoring iterations: 5

fit5 <- glm(y ~ group, data = d, family = poisson(link = "log")) #一般化線形モデル

fit6 <- summary(fit5)

fit6

## 

## Call:

## glm(formula = y ~ group, family = poisson(link = "log"), data = d)

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -1.43457  -0.62459   0.02401   0.47918   1.34113  

## 

## Coefficients:

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

## (Intercept)  2.93810    0.03255  90.272  < 2e-16 ***

## groupB       0.23578    0.04354   5.415 6.14e-08 ***

## ---

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

## 

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

## 

##     Null deviance: 80.618  on 99  degrees of freedom

## Residual deviance: 51.096  on 98  degrees of freedom

## AIC: 544.07

## 

## Number of Fisher Scoring iterations: 4

fit7 <- glm((n - y) ~ group, data = d, family = poisson(link = "log")) #一般化線形モデル

fit8 <- summary(fit7)

fit8

## 

## Call:

## glm(formula = (n - y) ~ group, family = poisson(link = "log"), 

##     data = d)

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.78678  -1.16403  -0.01987   0.53651   2.41242  

## 

## Coefficients:

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

## (Intercept)  1.91692    0.05423   35.35   <2e-16 ***

## groupB      -1.89712    0.15016  -12.63   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 346.99  on 99  degrees of freedom

## Residual deviance: 107.75  on 98  degrees of freedom

## AIC: 371.07

## 

## Number of Fisher Scoring iterations: 5

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

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

b0e

## [1] 1.021181

b1e

## [1] 2.132895

b0A <- coef(fit6)[1,1]

b1A <- coef(fit6)[2,1]

b0B <- coef(fit8)[1,1]

b1B <- coef(fit8)[2,1]