目次

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

・（第1種）ベータ分布：0～1の連続値をとる確率分布

・ベータ分布を使ったGLM（ベータ回帰）

・検出力

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

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

　これまで二項分布やオフセット項を考慮したGLMで割合データの統計モデリングを行ってきた。これらでも取り扱うことができないものが、連続値/連続値に基づく0～1の範囲の割合データである。具体的には、連続的な確率を取り扱う場合や、面積/面積のパーセンテージとか、が例として上がる。こういう場合はベータ分布を使って解析するのが便利だろう。

（第1種）ベータ分布：0～1の連続値をとる確率分布

　何かの連続的確率や連続値ベースのパーセンテージを表現するのに使える確率分布がベータ分布beta distributionである。この分布は0より大きく1より小さいの連続値を台とする。ベータ分布の形を決めるパラメータは形状母数α、βである。ベータ分布を表す確率密度関数Probability density functionはαとβをパラメータとして、下記のとおりである（形状の挙動は極めて複雑なので後述）。

Bはベータ関数であり、Γはガンマ関数である。ガンマ関数が階乗の一般化であったことを思い出すと、ベータ関数は組み合わせcombination計算の一般化であることがわかる。この定義に基づくと、この分布の平均はα/(α+β)、分散はαβ/((α+β)^2(α+β+1))とあらわせる。

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

ところで、実は後述するベータ分布を使ったGLMではα、βを推定するのではない。E(x) = α/(α+β) = μ、α+β=φとして推定しているため、以下のようにパラメータを置き換える。このうち、φはdispersion parameterと呼ばれ、二項分布におけるsize parameterに相当し、推定精度の指標ともなっている。

つまり、ベータ分布はφが大きくなれば分散が大きくなり、μが0.5に近づくほど分散が大きくなることがわかる。続いて、具体的に乱数を生成してデータの様子を確認しよう。rbeta関数を使うことで、ベータ分布に従うデータを生成できる。引数nは生成する個数、shape1と2は形状母数である。なお、通常のglm関数ではベータ分布を使った解析をできないので、ここでglmmTMBパッケージをロードする。

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

library(plotn)

library(viridis)

library(glmmTMB)

rbeta(n = 30, shape1 = 2, shape = 5)

##  [1] 0.29525745 0.44163978 0.31709819 0.18317000 0.46227996 0.32093776

##  [7] 0.18720403 0.38981954 0.43107867 0.17634799 0.46469506 0.34541792

## [13] 0.19842048 0.31316652 0.48775302 0.03240684 0.14264789 0.45756899

## [19] 0.11146020 0.06105098 0.28977653 0.52851361 0.04978864 0.32121917

## [25] 0.29895673 0.14286188 0.22934342 0.17423185 0.21498687 0.32249203

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

このように乱数を生成してみるとわかるが、生成されるデータはshape1/(shape1+shape2) = 2/7 = 0.28周りが多いことがわかるだろう。続いて、shaoeのパラメータによって確率分布がどうなるかをいくつか確認してみる。試しに、shape1 = 2で固定し、shape2には1, 1.5, 4を指定してみよう。

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

dd <- data.frame(X = seq(0.001, 0.999, length = 1000))

dd$shape2.1 <- dbeta(dd$X, shape1 = 2, shape2 = 1)

dd$shape2.15 <- dbeta(dd$X, shape1 = 2, shape2 = 1.5)

dd$shape2.4 <- dbeta(dd$X, shape1 = 2, shape2 = 4)

plotn(dd[,1], dd[,2:4], mode = "m", type = "l", 

      legend = T, leg.title = "shape", pch.leg = NA, lty.leg = 1,

      leg.lab = c("α = 2, β = 1", "α = 2, β = 1.5", "α = 2, β = 4"), 

      leg.sp = 8, ylab = "確率密度分布")#図1の描画

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

極めて特徴的なのは、形状母数が1の時で、この時は直線になる。次はshape1 = 1で固定し、shape2には1, 1.5, 4を指定してみよう。

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

dd <- data.frame(X = seq(0.001, 0.999, length = 1000))

dd$shape1.1 <- dbeta(dd$X, shape1 = 1, shape2 = 1)

dd$shape1.15 <- dbeta(dd$X, shape1 = 1, shape2 = 1.5)

dd$shape1.4 <- dbeta(dd$X, shape1 = 1, shape2 = 4)

plotn(dd[,1], dd[,2:4], mode = "m", type = "l", 

      legend = T, leg.title = "shape", pch.leg = NA, lty.leg = 1,

      leg.lab = c("α = 1, β = 1", "α = 1, β = 1.5", "α = 1, β = 4"), 

      leg.sp = 8, ylab = "確率密度分布")#図2の描画

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

形状母数が両方1になると、一様分布になる。次は、shape1 = 0.5で固定し、shape2には1, 1.5, 4を指定してみよう。

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

dd <- data.frame(X = seq(0.001, 0.999, length = 1000))

dd$shape05.1 <- dbeta(dd$X, shape1 = 0.5, shape2 = 1)

dd$shape05.15 <- dbeta(dd$X, shape1 = 0.5, shape2 = 1.5)

dd$shape05.4 <- dbeta(dd$X, shape1 = 0.5, shape2 = 4)

plotn(dd[,1], dd[,2:4], mode = "m", type = "l", 

      ylim = c(0, 5),

      legend = T, leg.title = "shape", pch.leg = NA, lty.leg = 1,

      leg.lab = c("α = 0.5, β = 1", "α = 0.5, β = 1.5", "α = 0.5, β = 4"), 

      leg.sp = 8, ylab = "確率密度分布")#図3の描画

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

最後は、shape1 = 0.5で固定し、shape2にも1以下の値を指定してみよう。

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

dd <- data.frame(X = seq(0.001, 0.999, length = 1000))

dd$shape05.01 <- dbeta(dd$X, shape1 = 0.5, shape2 = 0.1)

dd$shape05.05 <- dbeta(dd$X, shape1 = 0.5, shape2 = 0.5)

dd$shape05.09 <- dbeta(dd$X, shape1 = 0.5, shape2 = 0.9)

plotn(dd[,1], dd[,2:4], mode = "m", type = "l", 

      ylim = c(0, 5),

      legend = T, leg.title = "shape", pch.leg = NA, lty.leg = 1,

      leg.lab = c("α = 0.5, β = 0.1", "α = 0.5, β = 0.5", "α = 0.5, β = 0.9"), 

      leg.sp = 8, ylab = "確率密度分布")#図4の描画

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

両方のshapeパラメータが1より小さくなると、今までの確率分布とは異なる、中央がくぼんだような分布となる。このようにベータ分布は柔軟な確率の表現をすることが可能である。ここで、α+β=φを固定して、αとβの割合を変化させたときの挙動を見てみる。

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

phi <- 5

dd <- data.frame(X = seq(0.001, 0.999, length = 1000), 

                 shape02 = NA, 

                 shape08 = NA, 

                 shape10 = NA, 

                 shape20 = NA, 

                 shape30 = NA, 

                 shape40 = NA, 

                 shape45 = NA)

rep <- 1

for(i in c(0.2, 0.8, 1, 2, 3, 4, 4.5)){

  rep <- rep + 1

  dd[,rep] <- dbeta(dd$X, shape1 = phi - i, shape2 = i)

}

plotn(dd[,1], dd[,2:8], mode = "m", type = "l", 

      ylim = c(0, 5),

      legend = T, leg.title = "shape", pch.leg = NA, lty.leg = 1,

      leg.lab = c("α = 4.8, β = 0.2", 

                  "α = 4.2, β = 0.8", 

                  "α = 4.0, β = 1.0", 

                  "α = 3.0, β = 2.0", 

                  "α = 2.0, β = 3.0", 

                  "α = 1.0, β = 4.0", 

                  "α = 0.5, β = 4.5"), 

      leg.sp = 8, ylab = "確率密度分布")#図5の描画

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

このようにφを固定した時は、ちょうど二項分布の連続値化したような挙動となっている。

ベータ分布を使ったGLM（ベータ回帰）

　では、ベータ分布を使ったGLMについて考察してゆこう。ベータ分布を使う場合、Rのプログラム上では、リンク関数はデフォルトでロジット関数logit functionを使うことになっている。ロジット関数の逆関数はロジスティック関数logistic functionである。詳しくは二項分布のGLMの項で確認してほしい。。

リンク関数は元の関数の式を線形予測子に変換する式、リンク関数の逆関数が元の関数の形であるから、ベータ分布とロジットリンクをGLMで使うということは、元の関数はロジスティック関数を仮定していることになる。ベータ分布のGLMでは0<y<1を満たすyについて推定するのでロジスティック関数との相性が良い。ではこのような状況において、尤度を最大化するようなパラメータβを求める。データの生成過程を図示して様子をつかむことにする。 

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

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

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

ここで、ψ（プサイ）はディガンマ関数であり、対数ガンマ関数logΓ(x)を微分したもの、つまりψ(x) = Γ'(x)/Γ(x)である。非常に複雑な形（工夫すればもう少しだけ簡単になりそうだけど）であるが、後述するように、上記の式を満たすようなβ0となる。

　見た目だけなら下記のようにもう少しだけ簡単にできる（ディガンマ関数の性質を利用）。

二項分布の回帰の時の結果のアナロジーで考えると、β0の推定値は、log((yの総和)/(1-yの総和))から大きく外れることはないと期待できる。実際は、手計算すると異なっているがいったんはこれくらいの荒い推定値で考えてもいいかもしれない。

検出力

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

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

b0 <- -1

b1 <- 0

phi <- 5

n <- 30

x <- runif(n, 0, 30)

y <- rbeta(n = n, shape1 = phi/(1+exp(-b0-b1*x)), shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

plotn(y ~ x)#図6の描画

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

データの生成過程通り、xに対して特に関係のなさそうなデータが生成された。では、このデータをGLMで解析する。今回は、デフォルトのglm関数ではベータ分布を扱えないので、glmmTMBパッケージ内の、glmmTMB関数とbeta_family関数を使う。ここではベータ分布とロジットリンクを指定するので、beta_family(link="logit")と指定する。併せて、回帰曲線も描いてみよう。

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

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

fit2 <- summary(fit1)

fit2

##  Family: beta  ( logit )

## Formula:          y ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    -32.7    -28.5     19.3    -38.7       27 

## 

## 

## Dispersion parameter for beta family (): 5.87 

## 

## Conditional model:

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

## (Intercept) -0.65337    0.31076  -2.102   0.0355 *

## x           -0.04479    0.01976  -2.267   0.0234 *

## ---

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

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

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

phi_e <- fit2$sigma

plotn(y ~ x) #図7の描画

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

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

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

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

GLMの結果をsummary関数を通して得られる結果の解釈も線形回帰の時とそれほど変わらない。今回の場合、線形予測子の切片の推定値は-0.653、xの係数の推定値は-0.045であった。

　では、上記の設定でデータを1000回生成して、検出力を確認してみる。glmmTMBは計算がやや遅いので、1000回程度でとどめる。

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

b0 <- -1

b1 <- 0

phi <- 5

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- rbeta(n = n, shape1 = phi/(1+exp(-b0-b1*x)), 

             shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

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

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

  fit2 <- summary(fit1)

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

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

}

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

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

sum(p1 < 0.05)

## [1] 890

sum(p2 < 0.05)

## [1] 71

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

すると切片の係数の検出力は89%であることがわかる。一方で、xの係数の危険率は7%程度である。さらにシミュレートの必要があるかもしれないが、これくらいなら許容範囲かもしれない。

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

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

b0 <- -3

b1 <- 0.2

phi <- 5

n <- 30

x <- runif(n, 0, 30)

y <- rbeta(n = n, shape1 = phi/(1+exp(-b0-b1*x)), shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

plotn(y ~ x)#図10の描画

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

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

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

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

fit2 <- summary(fit1)

fit2

##  Family: beta  ( logit )

## Formula:          y ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    -55.8    -51.6     30.9    -61.8       27 

## 

## 

## Dispersion parameter for beta family (): 5.86 

## 

## Conditional model:

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

## (Intercept) -2.68506    0.36491  -7.358 1.86e-13 ***

## x            0.17948    0.02061   8.710  < 2e-16 ***

## ---

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

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

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

phi_e <- fit2$sigma

plotn(y ~ x) #図11の描画

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

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

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

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

線形予測子の切片の推定値は-2.69、xの係数の推定値は0.18であり、これも初めの設定と大きく変わらないことがわかる。

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

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

b0 <- -3

b1 <- 0.2

phi <- 5

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- rbeta(n = n, shape1 = phi/(1+exp(-b0-b1*x)), 

             shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

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

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

  fit2 <- summary(fit1)

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

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

}

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

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

sum(p1 < 0.05)

## [1] 1000

sum(p2 < 0.05)

## [1] 1000

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

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

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

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

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

y <- c(0.7, 0.55, 0.87, 0.83, 0.86, 0.83, 0.49, 0.65, 0.98, 

       0.52, 0.57, 0.69, 0.93, 0.61, 0.66, 0.93, 0.8, 0.77, 

       0.75, 0.81, 0.81, 0.49, 0.73, 0.83, 0.67, 0.86, 0.69,

       0.62, 0.7, 0.89, 0.67, 0.78, 0.71, 0.86, 0.36, 0.67, 

       0.78, 0.87, 0.63, 0.83, 0.66, 0.51, 0.83, 0.66, 0.72, 

       0.72, 0.57, 0.83, 0.88, 0.81, 0.83, 0.78, 0.89, 0.87, 

       0.84, 0.79, 0.5, 0.72, 0.76, 0.85, 0.58, 0.85, 0.81, 

       0.61, 0.82, 0.54, 0.57, 0.75, 0.63, 0.79, 0.49, 0.76,

       0.82, 0.87, 0.81, 0.7, 0.72, 0.92, 0.65, 0.8, 0.65, 

       0.65, 0.59, 0.73, 0.71, 0.81, 0.49, 0.96, 0.71, 0.51, 

       0.7, 0.88, 0.59, 0.26, 0.84, 0.77, 0.78, 0.8, 0.43, 0.82)

histn(y)#図14の描画

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

データは0より大きく、1より小さいの連続値データである。平均的な割合は0.8程度であろうか。こういう時はベータ分布を解析に使うとよい。一旦は練習なので、説明変数がないデータの解析をしてみることにする。これはデータが正規分布の時であれば、平均が有意に0から異なるかを調べる1標本のt検定に相当する。今回の場合は、切片の係数が0から有意に異なるかを考えているわけだから、前述の荒い推定に基づいて切片推定値がlog((yの総和)/(1-yの総和))であることを考慮すれば、(yの総和)と(1-yの総和)が有意に異なるかを検定していると言える。GLMで解析してみると下記のようになる。

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

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

fit2 <- summary(fit1)

fit2

##  Family: beta  ( logit )

## Formula:          y ~ 1

## 

##      AIC      BIC   logLik deviance df.resid 

##   -120.7   -115.5     62.3   -124.7       98 

## 

## 

## Dispersion parameter for beta family (): 9.79 

## 

## Conditional model:

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

## (Intercept)  0.96280    0.06786   14.19   <2e-16 ***

## ---

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

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

すると切片の係数は有意に0と異なっていることがわかる。今回の場合はベータ分布の平均をシンプルに示すことはできない。だが、上述の議論で、logit(y)の平均は、推定値を代入したディガンマ関数の差に等しい。この時、計算のためにdispersion parameterが必要で、それはfit2$sigmaに格納されている。この値は、glm関数でガンマ回帰をした時と異なり、最尤推定されている。

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

mean(logit(y))

## [1] 1.086624

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

phi_e <- fit2$sigma

digamma(phi_e/(1+exp(-b0e)))-digamma(phi_e*exp(-b0e)/(1+exp(-b0e)))

## [1] 1.086624

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

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

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

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

  sapply(b0t, function(tmp) sum(log(dbeta(x = y, shape1 = p/(1+exp(-tmp)), 

                                          shape2 = p*(1-1/(1+exp(-tmp)))

                                          )

                                    )

                                )

         )

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

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

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

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

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

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

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

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

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

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

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

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

glmmTMBの出力は、glmと異なり、対数尤度と逸脱度が表示される（まあ、こっちのほうが素直な気はする）。下記のように上記の関数を使って計算できる。併せてAICも計算しよう。ガンマ回帰の時と同様に、dispersion parameterが線形予測子以外に推定されているので、逸脱度への加算量は2×2である。

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

f(y, phi_e, b0e)

## [1] 62.34569

-2*(f(y, phi_e, b0e))

## [1] -124.6914

-2*f(y, phi_e, b0e) + 2*2

## [1] -120.6914

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

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

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

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)

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

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

では解析をしてみよう。解析ののち、回帰曲線を描く。

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

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

fit2 <- summary(fit1)

fit2

##  Family: beta  ( logit )

## Formula:          y ~ x

## 

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

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

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

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

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

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

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

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

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

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

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

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

  mapply(function(tmp1, tmp2) {

    sum(log(dbeta(x = y, shape1 = p/(1+exp(-tmp1-tmp2*x)), 

                  shape2 = p*(1-1/(1+exp(-tmp1-tmp2*x))))))

    }, b0t, b1t)

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

phi_e <- fit2$sigma

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, p = phi_e)

#対数尤度の計算

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

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

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

上記の関数を使って対数尤度、逸脱度を計算できる。併せてAICも計算しよう。ガンマ回帰の時と同様に、dispersion parameterが線形予測子以外に推定されているので、逸脱度への加算量は2×3である。

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

f(x, y, phi_e, b0e, b1e)

## [1] 75.13232

-2*f(x, y, phi_e, b0e, b1e)

## [1] -150.2646

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

## [1] -144.2646

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