区間推定

信頼区間と予測区間

何かの予測が収まる範囲？

　私たちはこれまでにさんざん、検定や回帰などを行ってきた。基本的にそこで取り上げてきたのは、どんな平均値が予測されるか、であった。回帰に対して、そういう印象はないかもしれないが、回帰曲線を引くという行為は、連続値の説明変数に対して被説明変数の平均値がどう推移するかを予測しているという行為に他ならない。このような平均値などの推定を点推定point estimationを呼んでいる。標本平均は母平均の点推定であるわけだが、言わずもがなぴったり一致することはない。であれば、点推定からどれくらいの範囲に母平均がありそうか、を知りたくなってくるだろう。このような範囲に注目した推定を区間推定interval estimationと呼ぶ。本項では、様々な区間推定を行ってみよう。

信頼区間・予測区間

　まず、頻度主義統計の考え方をおさらいしておこう。頻度主義統計では、データを生み出すモデルや構造を決めるパラメータはばらつくことがなく唯一の値が存在すると仮定している。例えば、平均1、分散1の正規分布からデータが生成されている、というような感じで、平均も分散もただ一つの値として固定されており、それが真の値であると考えている。これは当たり前のようにも感じるが、ベイズ統計学では当たり前ではないことを付記しておく。とりあえず、今までの解析では、なにかデータを生成するときは、パラメータを固定して、そこから生成していた。この考え方は頻度主義統計出る。

　このような前提の元、区間推定には2種類ある。それは、信頼区間と予測区間である。ともに点推定を「中心」とした範囲で示される。計算はさておき、まずは用語の説明だけ行ってしまう。

信頼区間confidence interval：

　モデルからデータをn回生成し、n回、それぞれで信頼区間を求めるとする。このとき、m回の信頼区間が真のモデルを含んでいた時、これを m/n × 100 %信頼区間と呼ぶ。

　何やらわかりにくいので具体例を出そう。平均1、分散1の正規分布からデータが生成されているとする。ここからサンプルサイズ3のデータを1000回生成することを考える。

1000回、標本平均と信頼区間を求めたとする。標本平均が母平均1の点推定である。標本平均を中心として信頼区間が求まる。このとき、950個の信頼区間が母平均を含むように範囲を設定したものが、その範囲が95%信頼区間である。もし、800個が母平均を含むような範囲を設定すれば80%信頼区間である。

　繰り返すが、母平均はただ一つの値だから、信頼区間を求めたとき、その信頼区間において、母平均は含まれるか、含まれないかの二択である。だから、95%信頼区間を求めたときに、「この範囲に95%の確率で母平均が含まれる」と言うことはできない。ただ、信頼区間を伸ばしていったときに、母平均が含まれている可能性は高まっているだろう。

予測区間prediction interval：

　モデルからさらにデータを生成したとき、それらのデータのうちa %が含まれる範囲をa %予測区間と呼ぶ。こちらはシンプルに、データが含まれるであろう範囲を予測する。

これらの区間推定はサンプルサイズに強く依存しており、信頼区間の幅は0に近づき、予測区間は仮定した確率分布の分散に近づく。これらの2つの区間推定を駆使することで、今回得られている回帰曲線などの信頼度を測ることができる、と考えて差し支えない。

1標本の区間推定

　まず、1標本の区間推定から始める。以下のようなデータの区間推定をしてみよう。

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

library(plotn)

library(glmmTMB)

y <- c(14.1, 11.1, 11.2, 10.3, 8.6, 9, 10.3, 5.3,

10.9, 11.2, 9.4, 10, 6.3, 11.6, 12, 9.7, 8.3,

11.7, 11, 14.1, 12.5, 9.4, 11.1, 9, 13.2, 10.5,

8.6, 7.8, 10.6, 8.6)

d <- data.frame(y = y)

boxplotn(y, data = d, jitter.method = "swarm")#図1の描画

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

図1 1標本のデータ

線形回帰とt検定を行ってみる。これらの解析結果内に、信頼区間などの情報が入っている。

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

fit1 <- lm(y ~ 1, data = d)

fit2 <- summary(fit1)

fit2

##

## Call:

## lm(formula = y ~ 1, data = d)

##

## Residuals:

## Min 1Q Median 3Q Max

## -4.9467 -1.2467 0.1533 0.9533 3.8533

##

## Coefficients:

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

## (Intercept) 10.247 0.368 27.84 <2e-16 ***

## ---

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

##

## Residual standard error: 2.016 on 29 degrees of freedom

fit3 <- t.test(y, mu = 0)

fit3

##

## One Sample t-test

##

## data: y

## t = 27.842, df = 29, p-value < 2.2e-16

## alternative hypothesis: true mean is not equal to 0

## 95 percent confidence interval:

## 9.493972 10.999362

## sample estimates:

## mean of x

## 10.24667

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

t検定の解析結果内にある、95 percent confidence intervalが95%信頼区間である。これは線形回帰のStd.Errorの値を使うことで下記のように求められる。qt関数はt分布のq値を求める関数である。t分布における95%の範囲で補正したものが、信頼区間なのである。

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

#信頼区間

lc_upr <- coef(fit2)[1] + qt(p = 0.975, df = df.residual(fit1)) * coef(fit2)[2]

lc_lwr <- coef(fit2)[1] - qt(p = 0.975, df = df.residual(fit1)) * coef(fit2)[2]

c_upr <- lc_upr

c_lwr <- lc_lwr

c_upr

## [1] 10.99936

c_lwr

## [1] 9.493972

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

そもそも、1標本の線形回帰のEstimateは平均、Std.Errorはデータの標準誤差そのものなので不偏分散をvar関数で計算すれば、下記のように計算できる。信頼区間とは、標準誤差をt分布で補正したものであると言える。なぜt分布で補正するかというと、分散が未知のパラメータであるため、平均の推定値が正規分布ではなくt分布に従うためである。

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

v <- var(y)

sqrt(v/length(d$y))

## [1] 0.3680247

mean(d$y)

## [1] 10.24667

#信頼区間

lc_upr <- mean(d$y) + qt(p = 0.975, df = df.residual(fit1)) * sqrt(v/length(d$y))

lc_lwr <- mean(d$y) - qt(p = 0.975, df = df.residual(fit1)) * sqrt(v/length(d$y))

c_upr <- lc_upr

c_lwr <- lc_lwr

c_upr

## [1] 10.99936

c_lwr

## [1] 9.493972

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

一方、予測区間は下記のように計算できる。

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

#予測区間

lp_upr <- mean(d$y) + qt(p = 0.975, df = df.residual(fit1)) *

sqrt(v*(1 + 1/length(d$y)))

lp_lwr <- mean(d$y) - qt(p = 0.975, df = df.residual(fit1)) *

sqrt(v*(1 + 1/length(d$y)))

p_upr <- lp_upr

p_lwr <- lp_lwr

p_upr

## [1] 14.43749

p_lwr

## [1] 6.055839

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

信頼区間では不偏分散vにサンプルサイズの逆数を掛けていたが、予測区間ではサンプルサイズの逆数に1を足した値を掛けている。これは平方根の中身を分解してみると、「不偏分散と標準誤差の2乗の和」になっている。中心極限定理を思い出すと、標準誤差とは平均値の従う分布だったのだから、もともとの確率分布に依存するばらつきと、平均値のばらつきの両方を考慮した推定になっているということである。

ゆえに、サンプルサイズが大きくなれば、信頼区間は0に近づき、予測区間は分散そのものになる。

　predict関数は信頼区間と予測区間を計算するうえで便利な関数である。線形回帰の場合、下記のようにすれば両方とも計算可能である。

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

pred <- predict(fit1, se.fit = T)

#信頼区間

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

#予測区間

lp_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) *

sqrt(pred$se.fit^2 + pred$residual.scale^2)

lp_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) *

sqrt(pred$se.fit^2 + pred$residual.scale^2)

c_upr <- lc_upr[1]

c_lwr <- lc_lwr[1]

p_upr <- lp_upr[1]

p_lwr <- lp_lwr[1]

res <- pred$fit[1]

c_upr

## 1

## 10.99936

c_lwr

## 1

## 9.493972

p_upr

## 1

## 14.43749

p_lwr

## 1

## 6.055839

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

predict関数の出力内には、$fitに点推定の値、$se.fitに標準誤差、$residual.scaleに標準偏差（不偏分散の平方根）が格納されている。

　あるいは、もっと簡単に下記のようにすれば計算してくれる。

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

pred1 <- predict(fit1, interval = "confidence", level = 0.95)

pred2 <- predict(fit1, interval = "prediction", level = 0.95)

pred1[1,]

## fit lwr upr

## 10.246667 9.493972 10.999362

pred2[1,]

## fit lwr upr

## 10.246667 6.055839 14.437495

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

これらを図示すれば以下のようにできるだろう。

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

lim <- range(c(pred2, d$y))

boxplotn(y, data = d, ylim = lim, jitter.method = "swarm")#図2の描画

overdraw(segments(0.7, pred2[1,2], 0.7, pred2[1,3], lwd = 2),

segments(0.7, pred1[1,2], 0.7, pred1[1,3], lwd = 4),

points(0.7, mean(d$y), cex = 1.5, pch = 16, col = "red"))

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

図2 1標本のデータ。赤い点が平均値、太い線が95%信頼区間、細い線が95%予測区間

線形回帰の区間推定

　続いて、線形回帰に関しての区間推定を行ってみよう。下記のようなデータを考える。

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

x <- c(7, 11, 27, 15, 5, 18, 9, 24, 4, 12,

19, 6, 8, 8, 24, 22, 11, 28, 0, 1,

12, 17, 10, 12, 9, 6, 18, 6, 29, 6)

y <- c(-20, -20, -53, -28, -4, -32, -11, -49, 3, -24,

-31, -11, -20, -12, -45, -45, -18, -48, 6, -3,

-22, -36, -15, -22, -26, -21, -21, -17, -53, -11)

d <- data.frame(x = x, y = y)

plotn(y ~ x, data = d)#図3の描画

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

図3 線形回帰のデータ

線形回帰を行ってみる。

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

fit1 <- lm(y ~ x, data = d) #回帰分析

fit2 <- summary(fit1)

fit2

##

## Call:

## lm(formula = y ~ x, data = d)

##

## Residuals:

## Min 1Q Median 3Q Max

## -10.0640 -3.8539 0.0378 3.2214 12.3430

##

## Coefficients:

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

## (Intercept) 0.2675 1.8326 0.146 0.885

## x -1.8673 0.1216 -15.354 3.64e-15 ***

## ---

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

##

## Residual standard error: 5.297 on 28 degrees of freedom

## Multiple R-squared: 0.8938, Adjusted R-squared: 0.89

## F-statistic: 235.7 on 1 and 28 DF, p-value: 3.64e-15

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

では、この例を使って信頼区間と予測区間を計算してみる。とりあえず、定義に従って下記のように計算できる。もちろん、線形回帰の切片と傾きはxとyの偏差平方和から計算可能なので、それで差し替えてもよい。

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

#回帰曲線を描く範囲

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

#不偏分散（yの実測値 - yの予測値の平方和を自由度で割った）

v <- sum((d$y - (coef(fit2)[1,1] + coef(fit2)[2,1]*d$x))^2)/df.residual(fit1)

#xの偏差平方和

Sxx <- var(d$x) * (length(d$y) - 1)

#予測するxの値と実測のx平均値からのずれ

tmp <- (xx-mean(d$x))^2

#信頼区間

lc_upr <- (coef(fit2)[1,1] + coef(fit2)[2,1]*xx) +

qt(p = 0.975, df = df.residual(fit1)) *

sqrt(v*(1/length(d$y) + tmp/Sxx))

lc_lwr <- (coef(fit2)[1,1] + coef(fit2)[2,1]*xx) -

qt(p = 0.975, df = df.residual(fit1)) *

sqrt(v*(1/length(d$y) + tmp/Sxx))

#予測区間

lp_upr <- (coef(fit2)[1,1] + coef(fit2)[2,1]*xx) +

qt(p = 0.975, df = df.residual(fit1)) *

sqrt(v*(1 + 1/length(d$y) + tmp/Sxx))

lp_lwr <- (coef(fit2)[1,1] + coef(fit2)[2,1]*xx) -

qt(p = 0.975, df = df.residual(fit1)) *

sqrt(v*(1 + 1/length(d$y) + tmp/Sxx))

c_upr <- lc_upr

c_lwr <- lc_lwr

p_upr <- lp_upr

p_lwr <- lp_lwr

head(c_upr)

## [1] 4.021409 3.718508 3.415707 3.113010 2.810417 2.507933

head(c_lwr)

## [1] -3.486473 -3.727796 -3.969219 -4.210745 -4.452376 -4.694115

head(p_upr)

## [1] 11.74833 11.46618 11.18415 10.90222 10.62040 10.33868

head(p_lwr)

## [1] -11.21339 -11.47547 -11.73766 -11.99995 -12.26235 -12.52486

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

図示は後に回して、計算の中身を見てみる。信頼区間の平方根の中身をばらしてみると、以下のような構成になっている。

つまり、データの線形回帰からのずれに関する標準誤差S.E.の2乗と傾きに関する回帰係数の分散の和である。標準誤差に関しては回帰直線の切片に関するばらつき、回帰係数の分散は傾きに関するばらつきととらえられて、この2つが合わさったものが回帰直線全体のばらつきを規定するということである。回帰直線は必ず点（E(x), E(y)）を通るから、回帰直線の中央からずれたxほど大きなずれを生み出すことを(x - E(x))^2を使って表していると言える。

　予測区間は上記の和にさらに不偏分散を足したものが平方根の中身になっている。サンプルサイズを極めて大きくすると、信頼区間は0に、予測区間は正規分布の分散そのものになる。

信頼区間と予測区間の他の計算方法を示しておく。predict関数を使ったとき、1標本の時と同様に$fitに点推定の値、$se.fitに標準誤差（標準誤差の2乗と回帰係数の分散の平方根をとったもの）、$residual.scaleに標準偏差（不偏分散の平方根）が格納されている。

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

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

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lp_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) *

sqrt(pred$se.fit^2 + pred$residual.scale^2)

lp_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) *

sqrt(pred$se.fit^2 + pred$residual.scale^2)

c_upr <- lc_upr

c_lwr <- lc_lwr

p_upr <- lp_upr

p_lwr <- lp_lwr

head(c_upr)

## 1 2 3 4 5 6

## 4.021409 3.718508 3.415707 3.113010 2.810417 2.507933

head(c_lwr)

## 1 2 3 4 5 6

## -3.486473 -3.727796 -3.969219 -4.210745 -4.452376 -4.694115

head(p_upr)

## 1 2 3 4 5 6

## 11.74833 11.46618 11.18415 10.90222 10.62040 10.33868

head(p_lwr)

## 1 2 3 4 5 6

## -11.21339 -11.47547 -11.73766 -11.99995 -12.26235 -12.52486

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

あるいはもっと簡単に下記のようにも計算できる。これを使って図示しよう。

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

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

newdata <- data.frame(x = xx)

pred1 <- predict(fit1, newdata = newdata, interval = "confidence", level = 0.95)

pred2 <- predict(fit1, newdata = newdata, interval = "prediction", level = 0.95)

head(pred1)

## fit lwr upr

## 1 0.267468016 -3.486473 4.021409

## 2 -0.004643808 -3.727796 3.718508

## 3 -0.276755633 -3.969219 3.415707

## 4 -0.548867457 -4.210745 3.113010

## 5 -0.820979281 -4.452376 2.810417

## 6 -1.093091106 -4.694115 2.507933

head(pred2)

## fit lwr upr

## 1 0.267468016 -11.21339 11.74833

## 2 -0.004643808 -11.47547 11.46618

## 3 -0.276755633 -11.73766 11.18415

## 4 -0.548867457 -11.99995 10.90222

## 5 -0.820979281 -12.26235 10.62040

## 6 -1.093091106 -12.52486 10.33868

lim <- range(c(pred2, d$y))

plotn(y ~ x, data = d, ylim = lim)#図4の描画

overdraw(polygon(c(xx, rev(xx)), c(pred1[,2], rev(pred1[,3])),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(pred2[,2], rev(pred2[,3])),

border = F, col = "#00000030"),

points(xx, pred1[,1], type = "l", col = "red"))

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

図4 線形回帰のデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

一般化線形モデルの区間推定：ポアソン回帰

　次は一般化線形モデルに関して信頼区間と予測区間を求めていこう。ただし、予測区間に関しては推定値の誤差を考慮しないものとする。というのも、そのあたりを考慮し始めるとかなり難しいようで、シンプルな計算ができるわけではないようだ。ciToolsパッケージを使う方法やyahoo知恵袋にブートストラップで計算する方法が紹介されていたりする。

　上記の線形回帰で示してきたように、サンプルサイズが十分大きければ、従う確率分布そのものになることがわかる。実際、自分でブートストラップで計算してみても、元の確率分布とあまり変わらないので、もとの確率分布の分散の影響がかなり大きいので、予測区間としてはそれで十分だろうと思われる。

　では以下のようなデータを解析する。

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

x <- c(13, 10.7, 8.1, 5.2, 22.9, 5.7, 20.1, 22.2, 4.7,

8.1, 30, 3.3, 19.6, 17.1, 3.9, 21.6, 0.9, 11.7,

4.9, 5.8, 29, 1.4, 26.1, 8.3, 24.7, 8.1, 15.8,

14.6, 29.2, 21.1, 20.2, 25.4, 22.5, 28, 17, 24.4,

10.1, 28.1, 11.3, 21.8, 9, 18.3, 13.7, 2.1, 11.3,

16.9, 0.7, 23.9, 18.5, 11.9)

y <- c(1, 0, 0, 0, 5, 1, 1, 2, 3, 0, 7, 0,

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

3, 1, 2, 1, 11, 2, 6, 9, 2, 5, 1, 4,

2, 3, 1, 5, 0, 4, 2, 0, 0, 4, 0, 3, 1, 0)

d <- data.frame(x = x, y = y)

plotn(y ~ x, data = d)#図5の描画

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

図5 ポアソン回帰のデータ

ではGLMによる解析を行う。

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

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

fit2 <- summary(fit1)

fit2

##

## Call:

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

##

## Deviance Residuals:

## Min 1Q Median 3Q Max

## -1.4696 -0.9610 -0.4634 0.4594 2.3778

##

## Coefficients:

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

## (Intercept) -1.16373 0.31234 -3.726 0.000195 ***

## x 0.10425 0.01366 7.634 2.27e-14 ***

## ---

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

##

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

##

## Null deviance: 124.484 on 49 degrees of freedom

## Residual deviance: 53.302 on 48 degrees of freedom

## AIC: 157.67

##

## Number of Fisher Scoring iterations: 5

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

そして、信頼区間と予測区間は下記のように計算、図示できる。同様にpredict関数が使えるが、glmやglmmTMBオブジェクトに対して、予測区間を計算してくれることはない。また、信頼区間に関しても計算に工夫が必要となる。リンク関数の逆関数による変換の前、つまり線形予測子の時点で推定値の誤差を考慮し、その後に変換を行うことで信頼区間を求めている。

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

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

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

#信頼区間

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_upr <- fit1$family$linkinv(lc_upr)#fit1オブジェクトのこの場所にlink関数の逆関数がある

c_lwr <- fit1$family$linkinv(lc_lwr)#fit1オブジェクトのこの場所にlink関数の逆関数がある

#回帰曲線

res <- fit1$family$linkinv(pred$fit)

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

予測区間はデータが従う確率分布、ポアソン分布の95%の範囲を表示することで代用している。離散分布なので滑らかな曲線にならないことに注意。無理やり指数関数で回帰して滑らかにつないでもよいだろうが、それだと予測区間 = データが実現する範囲であるにもかかわらず、実際にはその値が実現することがないため、意味のない値となるだろう。信頼区間のほうは回帰曲線が通る範囲なので、滑らかな曲線で間違いではない。

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

#予測区間

p_upr <- qpois(p = 0.975, lambda = res)

p_lwr <- qpois(p = 0.025, lambda = res)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図6の描画

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)),

border = F, col = "#00000030"),

points(xx, res, type = "l", col = "red"))

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

図6 ポアソン回帰のデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

一般化線形モデルの区間推定：ガンマ回帰

　では以下のようなデータを解析する。

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

x <- c(10.8, 14.2, 15.8, 7.1, 8.9, 22.4, 20.9, 18.1, 4.4, 4, 4.2,

11.7, 27.6, 19.2, 14.6, 29.4, 17.8, 3.2, 0.9, 0.8, 14.5,

16.6, 19.8, 23.6, 4.5, 7.7, 27.5, 18.5, 9, 17.9, 4, 28.3,

5, 1.1, 23.6, 22.8, 8.3, 27.5, 11.1, 5.7, 8.5, 5.2, 3.8,

12, 0, 16.8, 22.5, 6, 1.2, 24)

y <- c(0.66, 2.14, 1.16, 0.54, 0.69, 2.85, 2.86, 2.13, 0.37, 0.36,

0.42, 0.59, 6.52, 2.94, 1.54, 4.94, 1.25, 0.7, 0.41, 0.27,

1.67, 2.41, 3.76, 3.11, 0.37, 0.39, 5.98, 1.83, 0.77, 2.25,

0.92, 4.61, 0.39, 0.48, 3.34, 6.32, 1.06, 7.43, 1.85, 0.49,

1.11, 0.49, 0.38, 1.63, 0.26, 1.3, 3.39, 0.48, 0.34, 1.85)

d <- data.frame(x = x, y = y)

plotn(y ~ x, data = d)#図7の描画

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

図7 ガンマ回帰のデータ

ではGLMによる解析と信頼区間、予測区間の図示を合わせて行ってしまおう。

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

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

fit2 <- summary(fit1)

fit2

##

## Call:

## glm(formula = y ~ x, family = Gamma(link = "log"), data = d)

##

## Deviance Residuals:

## Min 1Q Median 3Q Max

## -0.6813 -0.2807 -0.1275 0.1917 0.7045

##

## Coefficients:

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

## (Intercept) -1.133524 0.088991 -12.74 <2e-16 ***

## x 0.104872 0.005682 18.46 <2e-16 ***

## ---

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

##

## (Dispersion parameter for Gamma family taken to be 0.1206351)

##

## Null deviance: 44.9458 on 49 degrees of freedom

## Residual deviance: 5.4295 on 48 degrees of freedom

## AIC: 55.971

##

## Number of Fisher Scoring iterations: 4

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

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_upr <- fit1$family$linkinv(lc_upr)

c_lwr <- fit1$family$linkinv(lc_lwr)

res <- fit1$family$linkinv(pred$fit)

p_upr <- qgamma(p = 0.975, shape = 1/fit2$dispersion, rate = 1/(fit2$dispersion*res))

p_lwr <- qgamma(p = 0.025, shape = 1/fit2$dispersion, rate = 1/(fit2$dispersion*res))

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図8の描画

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)),

border = F, col = "#00000030"),

points(xx, res, type = "l", col = "red"))

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

図8 ガンマ回帰のデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

一般化線形モデルの区間推定：二項回帰

　では以下のようなデータを解析する。

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

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)

d <- data.frame(x = x, n = n, y = y, rate = y/n)

plotn(rate ~ x, data = d)#図9の描画

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

図9 二項回帰のデータ

ではGLMによる解析と信頼区間、予測区間の図示を合わせて行ってしまおう。二項分布には、試行回数パラメータをあらかじめ決めておく必要があるが、今回は得られたデータの算術平均を丸めたもので代用している。

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

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

fit2 <- summary(fit1)

fit2

##

## Call:

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

## data = d)

##

## 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

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

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_upr <- fit1$family$linkinv(lc_upr)

c_lwr <- fit1$family$linkinv(lc_lwr)

res <- fit1$family$linkinv(pred$fit)

p_upr <- qbinom(p = 0.975, size = round(mean(n)), prob = res)/round(mean(n))

p_lwr <- qbinom(p = 0.025, size = round(mean(n)), prob = res)/round(mean(n))

lim <- range(c(c_upr, c_lwr, d$rate))

plotn(rate ~ x, data = d, ylim = lim)#図10の描画

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)),

border = F, col = "#00000030"),

points(xx, res, type = "l", col = "red"))

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

図10 二項回帰のデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

一般化線形モデルの区間推定：負の二項回帰

　では以下のようなデータを解析する。

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

x <- c(22.6, 23.1, 26.4, 5.8, 20.8, 10.8, 14.1, 24.9, 26.6,

11.6, 15.3, 1.9, 2.6, 13.6, 8.1, 27.7, 21.9, 10.7,

28.1, 4.8, 27.2, 1.2, 12.1, 19, 24.2, 22.3, 28.6,

21.4, 21, 20.4, 7.2, 27.6, 19.4, 17.5, 25.9, 6.9,

21.2, 27.6, 28.9, 22.7, 0.2, 15, 23.9, 6.3, 3.6,

1.8, 4.9, 10.5, 20, 21.9)

y <- c(12, 14, 5, 92, 20, 60, 28, 10, 14, 52, 26, 127, 93,

13, 91, 8, 14, 47, 13, 118, 12, 112, 53, 35, 5, 9,

8, 17, 31, 13, 80, 7, 13, 23, 8, 49, 12, 17, 2, 14,

134, 57, 7, 120, 52, 141, 96, 62, 39, 12)

d <- data.frame(x = x, y = y)

plotn(y ~ x, data = d)#図11の描画

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

図11 負の二項回帰のデータ

ではGLMによる解析と信頼区間、予測区間の図示を合わせて行ってしまおう。glmmTMBオブジェクトをpredict関数に通すときは、typeに関していくつか候補がある。linkにすればリンク関数による変換後（つまり線形予測子）の推定値の形で計算される。responseにすればリンク関数による返還前の推定値である。また、予測区間が一見すると滑らかな曲線に見えるが、全体的に値が大きく、そう見えてるだけである。拡大するとガタガタしている。

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

fit1 <- glmmTMB(y ~ x, data = d, family = nbinom2(link = "log")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## Family: nbinom2 ( log )

## Formula: y ~ x

## Data: d

##

## AIC BIC logLik deviance df.resid

## 379.1 384.8 -186.5 373.1 47

##

## Dispersion parameter for nbinom2 family (): 10.8

##

## Conditional model:

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

## (Intercept) 5.051185 0.101263 49.88 <2e-16 ***

## x -0.104092 0.005979 -17.41 <2e-16 ***

## ---

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

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

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T, type = "link")

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_upr <- fit1$modelInfo$family$linkinv(lc_upr)

c_lwr <- fit1$modelInfo$family$linkinv(lc_lwr)

res <- fit1$modelInfo$family$linkinv(pred$fit)

p_upr <- qnbinom(p = 0.975, size = fit2$sigma, prob = fit2$sigma/(fit2$sigma + res))

p_lwr <- qnbinom(p = 0.025, size = fit2$sigma, prob = fit2$sigma/(fit2$sigma + res))

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図12の描画

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)),

border = F, col = "#00000030"),

points(xx, res, type = "l", col = "red"))

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

図12 負の二項回帰のデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

一般化線形モデルの区間推定：ベータ回帰

　では以下のようなデータを解析する。

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

x <- c(2.1, 20.8, 14.6, 1.9, 18.1, 16.9, 10.8, 9.1, 27.8,

21.4, 15.7, 2.2, 29.7, 9.5, 17.9, 9.6, 23.1, 8.7,

8.2, 12.3, 25.6, 22, 1.4, 22.4, 16.2, 6.4, 26.7,

25.8, 9.8, 25.9, 22.7, 15.5, 2.2, 4.6, 8.1, 10.6,

1.5, 21.3, 2.4, 13, 8.8, 8.4, 24.6, 0, 21.8, 27.4,

9.9, 19.6, 25.8, 22.9)

y <- c(0.663, 0.115, 0.413, 0.609, 0.132, 0.199, 0.305, 0.405,

0.021, 0.019, 0.249, 0.525, 0.009, 0.289, 0.159, 0.476,

0.067, 0.454, 0.407, 0.291, 0.082, 0.1, 0.689, 0.082,

0.07, 0.543, 0.012, 0.03, 0.443, 0.057, 0.084, 0.226,

0.532, 0.476, 0.367, 0.331, 0.591, 0.072, 0.593, 0.305,

0.473, 0.547, 0.136, 0.84, 0.076, 0.028, 0.307, 0.204, 0.008, 0.071)

d <- data.frame(x = x, y = y)

plotn(y ~ x, data = d)#図13の描画

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

図13 ベータ回帰のデータ

ではGLMによる解析と信頼区間、予測区間の図示を合わせて行ってしまおう。

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

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

fit2 <- summary(fit1)

fit2

## Family: beta ( logit )

## Formula: y ~ x

## Data: d

##

## AIC BIC logLik deviance df.resid

## -144.3 -138.5 75.1 -150.3 47

##

## Dispersion parameter for beta family (): 37

##

## Conditional model:

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

## (Intercept) 0.898835 0.101678 8.84 <2e-16 ***

## x -0.145935 0.007592 -19.22 <2e-16 ***

## ---

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

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

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T, type = "link")

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_upr <- fit1$modelInfo$family$linkinv(lc_upr)

c_lwr <- fit1$modelInfo$family$linkinv(lc_lwr)

res <- fit1$modelInfo$family$linkinv(pred$fit)

p_upr <- qbeta(p = 0.975, shape1 = fit2$sigma*res, shape2 = fit2$sigma*(1-res))

p_lwr <- qbeta(p = 0.025, shape1 = fit2$sigma*res, shape2 = fit2$sigma*(1-res))

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図14の描画

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)),

border = F, col = "#00000030"),

points(xx, res, type = "l", col = "red"))

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

図14 ベータ回帰のデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

一般化線形モデルの区間推定：ベータ-二項回帰

　では以下のようなデータを解析する。

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

x <- c(5.6, 19.9, 3.1, 2.6, 0.3, 5, 13.7, 27.3, 19.8, 9.5, 8,

9.4, 22.8, 5.9, 29, 20.5, 22.3, 14.1, 3.8, 1.4, 2.9,

7.7, 23.5, 17.1, 7, 0.3, 28.6, 26.3, 8.9, 3.4, 3.2, 9.2,

5, 27.7, 5.8, 25.4, 1, 17.4, 28.4, 10.2, 8.9, 18.7, 1.8,

7, 3.5, 24.5, 0.1, 18.1, 16.1, 4.5)

n <- c(26, 27, 30, 23, 20, 27, 21, 30, 24, 21, 30, 28, 23, 28,

27, 30, 30, 22, 26, 30, 21, 29, 29, 22, 26, 24, 28, 20,

20, 26, 28, 28, 24, 27, 28, 26, 26, 22, 26, 27, 29, 20,

21, 25, 22, 30, 25, 20, 29, 24)

y <- c(15, 5, 7, 16, 17, 10, 2, 0, 6, 7, 15, 6, 3, 8, 0, 5, 6,

5, 12, 22, 20, 13, 1, 2, 9, 17, 0, 2, 15, 14, 18, 6, 18,

0, 14, 3, 22, 5, 6, 10, 14, 3, 18, 16, 18, 0, 14, 2, 15, 19)

d <- data.frame(x = x, n = n, y = y, rate = y/n)

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

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

図15 ベータ-二項回帰のデータ

ではGLMによる解析と信頼区間、予測区間の図示を合わせて行ってしまおう。ベータ-二項分布に関する関数はデフォルトで用意されていないので、今回は自前で関数を定義して使ってみる。

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

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

fit2 <- summary(fit1)

fit2

## Family: betabinomial ( logit )

## Formula: cbind(y, n - y) ~ x

## Data: d

##

## AIC BIC logLik deviance df.resid

## 258.6 264.4 -126.3 252.6 47

##

## Dispersion parameter for betabinomial family (): 11.6

##

## Conditional model:

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

## (Intercept) 0.99295 0.17604 5.641 1.7e-08 ***

## x -0.13911 0.01528 -9.104 < 2e-16 ***

## ---

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

rbeta_binom <- function(n, size, mu, phi){

rbinom(n = n, size = size, prob = rbeta(n = n, shape1 = mu*phi, shape2 = (1-mu)*phi))

}

dbeta_binom <- function(x, size, mu, phi, log = F){

f <- function(x, size, mu, phi){

choose(size, x)*(gamma(mu*phi + x)/gamma(mu*phi))*(gamma((1-mu)*phi + size - x)/gamma(((1-mu)*phi)))*(gamma(phi)/gamma(phi+size))

}

d <- mapply(f, x = x, size = size, mu = mu, phi = phi)

if(log == T) d <- log(d)

return(d)

}

pbeta_binom <- function(q, size, mu, phi, lower.tail = T, log.p = F){

f <- function(q, size, mu, phi){

tmp <- sum(dbeta_binom(x = 0:floor(q), size, mu, phi))

tmp[q == size] <- 1

return(tmp)

}

p <- mapply(f, q = q, size = size, mu = mu, phi = phi)

if(lower.tail == F) p <- 1 - p

if(log.p == T) p <- log(p)

return(p)

}

qbeta_binom <- function(p, size, mu, phi, lower.tail = T, log.p = F){

if(log.p == T) p <- exp(p)

f <- function(p, size, mu, phi, lower.tail){

tmp <- pbeta_binom(0:size, size, mu, phi, lower.tail)

if(lower.tail == T){

sum(tmp <= p)

} else {

size + 1 - sum(tmp <= p)

}

q <- mapply(f, p = p, size = size, mu = mu, phi = phi, lower.tail = lower.tail)

q[q > size] <- size

return(q)

}

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

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T, type = "link")

lc_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

lc_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_upr <- fit1$modelInfo$family$linkinv(lc_upr)

c_lwr <- fit1$modelInfo$family$linkinv(lc_lwr)

res <- fit1$modelInfo$family$linkinv(pred$fit)

p_upr <- qbeta_binom(p = 0.975, size = round(1/mean(1/n)), mu = res, phi = fit2$sigma)/round(1/mean(1/n))

p_lwr <- qbeta_binom(p = 0.025, size = round(1/mean(1/n)), mu = res, phi = fit2$sigma)/round(1/mean(1/n))

lim <- range(c(p_upr, p_lwr, d$rate))

plotn(rate ~ x, data = d, ylim = lim)#図16の描画

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)),

border = F, col = "#00000030"),

points(xx, res, type = "l", col = "red"))

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

図16 ベータ-二項回帰のデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

ハードルモデルの区間推定

　次にハードルモデルに関して、区間推定を行ってみる。2つのモデルの組み合わせのようなモデルであり、こちらもぱっと調べた限りで信頼区間や予測区間を計算する方法が出てくるわけではない。私なりに考えて出した結果を紹介することにする。試しにシミュレーションしてみると大外れではないようなので、一つの案としてとらえてほしい。

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

x <- c(20.3, 26.1, 17.6, 6.6, 24.7, 7.9, 3.5, 11.6, 7.9, 21.3,

9.7, 12.9, 28.1, 16.3, 15.2, 25.3, 13.6, 2.1, 0.3, 26.5,

0.7, 10.4, 16, 3.6, 29.3, 12, 20.9, 28.1, 23.1, 25.8, 9.1,

13.2, 18.9, 16.6, 0.3, 18, 9.8, 1.1, 11.3, 4.1, 3.6, 15,

8.2, 23.8, 7, 21.8, 7.4, 4.2, 12.3, 20.7, 16.5, 2.1, 5.9,

2.9, 5.6, 22.7, 26.5, 22.3, 9.2, 7.4, 1.7, 18.6, 9.1, 2.5,

19.6, 11.9, 7.5, 25.9, 1.1, 20.5, 5, 14.1, 22.4, 10, 15.9,

15.2, 21.8, 28.2, 26.4, 4.6, 5.9, 21.6, 11.6, 21.9, 29.8,

6.4, 22, 1.1, 24.2, 12.6, 26.5, 22.1, 26.6, 23.8, 5.5, 2.6,

4.2, 19.3, 8.5, 6.7)

y <- c(12.08, 24.79, 6.79, 0, 21.08, 4.6, 7.08, 3.08, 6.09,

2.83, 0, 0, 8.07, 2.14, 1.61, 6.74, 12.34, 0, 2.15,

8.37, 2.12, 0, 9.55, 0, 6.14, 0, 3.98, 38.71, 0, 9.86,

0, 1.69, 3.07, 2.96, 0, 0, 0, 0, 4.08, 5.39, 4.66, 14.13,

3.68, 4.72, 0, 8.96, 12.16, 0, 0, 7.01, 11.48, 0, 0, 0,

0, 3.06, 15.42, 8.49, 0, 3.67, 0.91, 0, 0, 0, 13.79, 9.93,

1.74, 13.25, 4.38, 0, 0, 12.58, 8.58, 0, 0, 7.27, 2.08,

13.02, 4.46, 0.87, 0, 8.13, 1.61, 2.26, 14.38, 1.33, 15.57,

0, 5.52, 0, 12.63, 10.44, 12.28, 4.85, 0, 0, 2.48, 9.22, 8.68, 0)

d <- data.frame(x = x, y = y)

plotn(y ~ x, data = d)#図17の描画

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

図17 ハードルモデルのデータ

ではGLMによる解析と信頼区間、予測区間の図示を合わせて行ってしまおう。

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

fit1 <- glmmTMB(y ~ x, data = d, family = ziGamma(link = "log"), ziformula = ~ x) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## Family: Gamma ( log )

## Formula: y ~ x

## Zero inflation: ~x

## Data: d

##

## AIC BIC logLik deviance df.resid

## 482.8 495.8 -236.4 472.8 95

##

## Dispersion estimate for Gamma family (sigma^2): 0.408

##

## Conditional model:

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

## (Intercept) 1.129165 0.175778 6.424 1.33e-10 ***

## x 0.049990 0.009299 5.376 7.61e-08 ***

## ---

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

##

## Zero-inflation model:

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

## (Intercept) 1.07043 0.43109 2.483 0.013 *

## x -0.14031 0.03299 -4.253 2.11e-05 ***

## ---

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

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

newdata <- data.frame(x = xx)

pred1 <- predict(fit1, newdata = newdata, se.fit = T, type = "link")#ガンマ回帰の信頼区間

pred2 <- predict(fit1, newdata = newdata, se.fit = T, type = "zlink")#ベルヌーイ回帰の信頼区間

crit <- 0.025

#ガンマ回帰の信頼区間

lc_upr1 <- pred1$fit + qt(p = 1 - crit, df = df.residual(fit1)) * pred1$se.fit

lc_lwr1 <- pred1$fit - qt(p = 1 - crit, df = df.residual(fit1)) * pred1$se.fit

#ベルヌーイ回帰の信頼区間

lc_upr2 <- pred2$fit - qt(p = 1 - crit, df = df.residual(fit1)) * pred2$se.fit

lc_lwr2 <- pred2$fit + qt(p = 1 - crit, df = df.residual(fit1)) * pred2$se.fit

c_upr1 <- fit1$modelInfo$family$linkinv(lc_upr1)

c_lwr1 <- fit1$modelInfo$family$linkinv(lc_lwr1)

c_upr2 <- 1/(1+exp(lc_upr2))

c_lwr2 <- 1/(1+exp(lc_lwr2))

#両方を考慮した信頼区間

c_upr <- c_upr1 * c_upr2

c_lwr <- c_lwr1 * c_lwr2

res1 <- fit1$modelInfo$family$linkinv(pred1$fit)

res2 <- 1/(1+exp(pred2$fit))

res <- res1 * res2

#ガンマ回帰の予測区間

p_upr1 <- qgamma(p = 1-crit, shape = fit2$sigma^(-2), rate = fit2$sigma^(-2)/res1)

p_lwr1 <- qgamma(p = crit, shape = fit2$sigma^(-2), rate = fit2$sigma^(-2)/res1)

#ベルヌーイ回帰の予測区間

p_upr2 <- qbinom(p = 1-crit, size = 1, prob = res2)

p_lwr2 <- qbinom(p = crit, size = 1, prob = res2)

#両方を考慮した予測区間

p_upr <- p_upr1 * p_upr2

p_lwr <- p_lwr1 * p_lwr2

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図18の描画

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)),

border = F, col = "#00000030"),

polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)),

border = F, col = "#00000030"),

points(xx, res, type = "l", col = "red"))

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

図18 ハードルモデルのデータ。赤い線が回帰直線、濃い領域が95%信頼区間、薄い領域が95%予測区間

一般化線形混合モデルの区間推定

　（予定）

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