目次

・被説明変数が正の連続値データのモデリング

・ガンマ分布：正の連続値をとる確率分布

・ガンマ分布を使ったGLM（ガンマ回帰）

・検出力

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

被説明変数が正の連続値データのモデリング

　これまではポアソン分布と二項分布のような離散値のモデリングをおこなってきた。本項では、連続値データを扱ってみる。例えば、何かの長さや重さをモデリングすることを考えてみよう。これらのデータは1、2、3のような離散値をとらないことは明白であろう。一方で、これらは連続値であっても、-1.5とか0とかのような0以下の値も取らない。以上から、こういうデータは正の連続値をとるような確率分布で記述できることが望ましい。正規分布は負の∞まで値をとるので現実的ではない。

ガンマ分布：正の連続値をとる確率分布

　そこで正の連続値データを解析するときによく用いられる分布がガンマ分布Gamma distributionである。この確率分布は正の連続値を台とする。0の値はとらないことに注意しよう。ガンマ分布の形を決めるパラメータは形状母数shape parameter (k)と比母数 rate parameter (λ)である。あるいは比母数の逆数1/λ＝尺度母数 scale parameter (θ)を用いることもある。ガンマ分布を表す確率密度関数Probability density functionはk、λ、θをパラメータとして、下記のとおりである。

まず特徴的なのはk ≦ 1の時、ガンマ分布は0付近の確率密度が高く、軸にベタっと張り付いたような形状になる（ただし、0の値はとらない）。k > 1の時は軸から離れた一山形分布になる。そして、ともに正の∞方向に裾が重い分布になっている。ここでΓ（ガンマ）とはガンマ関数のことを指し、以下のように定義される。引数が自然数nの場合は階乗と同じ値になるが、ガンマ関数は実数でも定義されている（もっと言えば複素数でも定義されている）。つまりは、階乗の一般化と考えてよいだろう。

式の形はいかついが、慣れてしまえば以下のような階乗計算と類似した特徴が出てくることがわかる。

ガンマ関数の特徴はこれくらいにして、ガンマ分布の特徴に戻るとガンマ分布は平均k/λ、分散k/(λ^2)となる。以下のように計算できる。

つまり、kが大きくなるほど平均や分散は大きくなり、λが大きくなるほど平均や分散は小さくなることがわかる。また、分散 = 平均/λであり、平均が大きくなれば分散も大きくなる。続いて、具体的に乱数を生成してデータの様子を確認しよう。rgamma関数を使うことで、ガンマ分布に従うデータを生成できる。引数nは生成する個数、shapeは形状母数、rateは比母数である。

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

library(plotn)

library(viridis)

rgamma(n = 30, shape = 1.1, rate = 2)

##  [1] 0.008165513 0.113495935 0.361143033 0.118370205 0.743159890 0.286758591

##  [7] 0.508853059 0.595750096 0.107773895 0.437215513 0.319024899 0.576628757

## [13] 1.060932064 0.809948275 1.549132817 0.777316862 0.044942524 0.158238673

## [19] 0.232579043 0.367149078 0.688513754 0.161376490 1.125810372 0.343815438

## [25] 0.026769745 0.299617486 0.230764487 0.719233764 0.521872731 0.005557594

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

このように乱数を生成してみるとわかるが、生成されるデータはshape/rate = 0.55周りが多いことがわかるだろう。続いて、shapeのパラメータによって確率分布がどうなるかをいくつか確認してみる。試しに、k = 0.9, 1.0, 1.1を指定してみよう。

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

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

dd$shape09 <- dgamma(dd$X, shape = 0.9, scale = 1)

dd$shape10 <- dgamma(dd$X, shape = 1, scale = 1)

dd$shape11 <- dgamma(dd$X, shape = 1.1, scale = 1)

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

      leg.lab = c(0.9, 1.0, 1.1), pch.leg = NA, lty.leg = 1, leg.sp = 5, ylab = "確率密度分布")#図1の描画

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

さらにshapeのパラメータを大きくするとどうなるか確認してみる。k = 1, 10, 30, 60を指定してみよう。

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

dd <- data.frame(X = seq(0.01, 100, length = 1000))

dd$shape1 <- dgamma(dd$X, shape = 1, scale = 1)

dd$shape10 <- dgamma(dd$X, shape = 10, scale = 1)

dd$shape30 <- dgamma(dd$X, shape = 30, scale = 1)

dd$shape60 <- dgamma(dd$X, shape = 60, scale = 1)

plotn(dd[,1], dd[,2:5], mode = "m", type = "l", legend = T, leg.title = "shape",

      leg.lab = c(1, 10, 30, 60), pch.leg = NA, lty.leg = 1, leg.sp = 5, ylab = "確率密度分布")#図2の描画

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

次にscaleパラメータを変えるとどうなるか確認してみる。θ = 1/10, 1/30, 1/60を指定してみよう。

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

dd <- data.frame(X = seq(0.01, 2, length = 1000))

dd$scale10 <- dgamma(dd$X, shape = 10, scale = 1/10)

dd$scale30 <- dgamma(dd$X, shape = 10, scale = 1/30)

dd$scale60 <- dgamma(dd$X, shape = 10, scale = 1/60)

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

      leg.lab = c("1/10", "1/30", "1/60"), pch.leg = NA, 

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

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

以上のように形状母数や尺度母数が大きくなるほど、ガンマ分布は左右対称な分布に近づく。実際、ある程度大きくなれば正規分布への近似は可能であり、平均値が大きいデータであれば必ずしもガンマ分布を使ってモデリングをする必要はない。ガンマ分布が威力を発揮するのは、ポアソン分布と同様に0付近のデータである。

ガンマ分布を使ったGLM（ガンマ回帰）

　では、ガンマ分布を使ったGLMについて考察してゆこう。このような回帰分析はガンマ回帰Gamma regressionとも呼ばれる。ガンマ分布を使う場合、Rのプログラム上では、リンク関数はデフォルトで逆数関数inverse functionを使うことになっている。これはつまり、1/xというように逆数にする変換で、逆数関数の逆関数は同様に逆数関数である。ただし、逆数関数はガンマ分布と相性が良い部分もある一方、逆数関数に不連続点があり、不連続点付近では無限大に発散するなど、実数全体の予測に向かない不便な点もある。ガンマ分布の性質上、実数全体で定義される指数関数とも相性が良いので、リンク関数に対数関数を指定するのも手である。

逆数リンクの場合、下記のようなモデル式になる。

対数リンクの場合、下記のようなモデル式になる。

では、尤度を最大化するようなパラメータβを求めるために、尤度関数を求める。データの生成過程を図示して様子をつかむことにする。 

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

そして対数尤度は下記のようになり、それの式を整理する。まずは逆数リンクを考える。

対数尤度が説明変数xを含まなければ下記のようにβ0を解析的に解くことができる。

対数リンクの場合は下記のとおりである。

対数尤度が説明変数xを含まなければ下記のようにβ0を解析的に解くことができる。

以上のように、逆数リンクおよび対数リンクの時、それぞれ対数尤度を最大にするb0は1/b0 = (yの平均)あるいはexp(b0) = (yの平均)を満たすことがわかる。

検出力

　ではRを使って、GLMの検出力に関して検討しよう。以下のように、線形予測子の切片の係数は0ではないが、説明変数xの係数は0である場合をまず考える。shapeは一定とを考える。今回はshape = 1とする。リンク関数は対数リンクを想定するので、yを生成するときの引数rateは、shape/rate = exp(線形予測子)であるので、rate = shape/exp(線形予測子)と指定する。データを生成した後で、描画してみよう。

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

b0 <- 1

b1 <- 0

shape <- 1

n <- 30

x <- runif(n, 0, 30)

y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

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

データの生成過程通り、xに対して特に関係のなさそうなデータが生成された。ガンマ分布の右側の裾が重いので、外れ値が生成しやすいことも注目に値する。では、このデータをGLMで解析する。ここではガンマ分布と対数リンクを指定するので、Gamma(link="log")と指定する。併せて、回帰曲線も描いてみよう。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.9331  -0.8113  -0.2619   0.3709   1.9392  

## 

## Coefficients:

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

## (Intercept) 0.935967   0.340218   2.751   0.0103 *

## x           0.003714   0.019239   0.193   0.8483  

## ---

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

## 

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

## 

##     Null deviance: 36.325  on 29  degrees of freedom

## Residual deviance: 36.286  on 28  degrees of freedom

## AIC: 125.89

## 

## Number of Fisher Scoring iterations: 6

plotn(y ~ x)

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

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

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

b0 <- 1

b1 <- 0

shape <- 1

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

  fit2 <- summary(fit1)

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

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

}

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

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

sum(p1 < 0.05)

## [1] 7243

sum(p2 < 0.05)

## [1] 681

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

予測通り、切片の検出力は72.4%、xの係数の危険率は6.8%程度になった。やや、危険率が高いだろうか。ガンマ分布の特性として知っておいてもよいかもしれない。

　次にshape = 10としてデータを生成して解析してみよう。

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

b0 <- 1

b1 <- 0

shape <- 10

n <- 30

x <- runif(n, 0, 30)

y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

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

データの生成過程通り、xに対して特に関係のなさそうなデータが生成された。では、このデータをGLMで解析する。併せて、回帰曲線も描いてみよう。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -0.52128  -0.21719  -0.05666   0.25231   0.40938  

## 

## Coefficients:

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

## (Intercept) 0.862429   0.109988   7.841 1.53e-08 ***

## x           0.007906   0.005874   1.346    0.189    

## ---

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

## 

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

## 

##     Null deviance: 2.1548  on 29  degrees of freedom

## Residual deviance: 2.0320  on 28  degrees of freedom

## AIC: 68.372

## 

## Number of Fisher Scoring iterations: 4

plotn(y ~ x)

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

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

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

b0 <- 1

b1 <- 0

shape <- 1

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

  fit2 <- summary(fit1)

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

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

}

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

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

sum(p1 < 0.05)

## [1] 10000

sum(p2 < 0.05)

## [1] 499

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

予測通り、切片の検出力は100%、xの係数の危険率は5%程度になった。

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

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

b0 <- 1

b1 <- 0.05

shape <- 1

n <- 30

x <- runif(n, 0, 30)

y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

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

今回はxが増加するとともに、yも増加している……ような気がする。では、このデータをGLMで解析する。併せて、回帰曲線も描いてみよう。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.1724  -0.8158  -0.3550   0.2193   1.6546  

## 

## Coefficients:

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

## (Intercept)  1.29272    0.34482   3.749  0.00082 ***

## x            0.03827    0.01912   2.001  0.05516 .  

## ---

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

## 

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

## 

##     Null deviance: 27.630  on 29  degrees of freedom

## Residual deviance: 25.203  on 28  degrees of freedom

## AIC: 178.4

## 

## Number of Fisher Scoring iterations: 7

plotn(y ~ x)

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

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

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

b0 <- 1

b1 <- 0.05

shape <- 1

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

  fit2 <- summary(fit1)

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

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

}

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

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

sum(p1 < 0.05)

## [1] 7310

sum(p2 < 0.05)

## [1] 6339

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

今回は、切片の検出力は73.1%、xの係数の危険率は63.3%程度になった。

　次にshape = 10としてデータを生成して解析してみよう。

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

b0 <- 1

b1 <- 0.05

shape <- 10

n <- 30

x <- runif(n, 0, 30)

y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

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

今回はxが増加するとともに、yも増加している傾向が見て取れる。では、このデータをGLMで解析する。併せて、回帰曲線も描いてみよう。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -0.51451  -0.29884  -0.05075   0.16503   0.58826  

## 

## Coefficients:

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

## (Intercept) 1.044342   0.109786   9.513 2.87e-10 ***

## x           0.051101   0.006147   8.313 4.81e-09 ***

## ---

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

## 

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

## 

##     Null deviance: 9.1311  on 29  degrees of freedom

## Residual deviance: 2.7558  on 28  degrees of freedom

## AIC: 126.4

## 

## Number of Fisher Scoring iterations: 4

plotn(y ~ x)

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

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

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

b0 <- 1

b1 <- 0.05

shape <- 10

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

  y <- rgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0))

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

  fit2 <- summary(fit1)

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

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

}

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

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

sum(p1 < 0.05)

## [1] 10000

sum(p2 < 0.05)

## [1] 10000

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

すると切片の係数の検出力は100%、xの係数の検出力は100%である。全体を通してみてみると、k ≦ 1のときは、検出感度が低く、危険率が高くなりがちなので、データのばらつき方の様子は注意が必要である。

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

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

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

y <- c(4.78, 1.41, 0.24, 2.56, 2.74, 0.05, 2.42, 0.74, 2.1, 

       0.08, 1.39, 1.53, 0.39, 7.54, 0.31, 3.01, 2.85, 1.22, 

       3.19, 2.96, 4.03, 1.29, 2.59, 3.72, 3.04, 0.001, 1.15, 0.4,

       1.13, 1.11, 1.32, 0.83, 1.47, 0.13, 2.5, 1.52, 0.54,

       4.84, 4.43, 3.45, 2.74, 3.02, 0.9, 0.69, 0.36, 1, 6.39, 

       1.08, 0.38, 1.13, 0.32, 4.23, 0.86, 1.1, 2.18, 6.48, 0.22, 

       4.27, 7.49, 7.33, 0.86, 2.32, 3.28, 2.2, 3.44, 8.45, 1.68, 

       0.61, 2.07, 0.59, 2.76, 2.23, 0.8, 1.65, 0.01, 1.47, 0.68, 

       3.65, 0.58, 0.48, 6.05, 0.8, 0.51, 2.37, 0.51, 0.81, 1.42, 

       1.29, 0.42, 1.57, 0.64, 2.02, 4.1, 0.1, 1.44, 1.35, 5.14, 3.29, 0.26, 1.78)

histn(y)#図20の描画

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

データは0より大きいの連続値データである。データは0周りに分布し、かつ特に0付近の頻度が大きくなっている。このような分布はk ≦ 1のときのガンマ分布の特徴である。一旦は練習なので、説明変数がないデータの解析をしてみることにする。これはデータが正規分布の時であれば、平均が有意に0から異なるかを調べる1標本のt検定に相当する。今回の場合は、切片の係数が0から有意に異なるか、を考えている。ガンマ分布の対数リンクで考えるならば、ガンマ分布の平均がexp(0) = 1から有意に異なるかを調べていると考えてもよいだろう。ここでもし、逆数リンクを考えるのだとしたら、ガンマ分布の平均が1/0ではない……、というようなかなり怪しい解析をすることになるので注意が必要である。この辺りはやはり対数リンクのほうが解析しやすいと思う。GLMで解析してみると下記のようになる。GLMで解析してみると下記のようになる。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ 1, family = Gamma(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.6455  -0.9384  -0.3410   0.3750   1.8156  

## 

## Coefficients:

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

## (Intercept)  0.73645    0.09019   8.166 1.06e-12 ***

## ---

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

## 

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

## 

##     Null deviance: 107.33  on 99  degrees of freedom

## Residual deviance: 107.33  on 99  degrees of freedom

## AIC: 352.16

## 

## Number of Fisher Scoring iterations: 5

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

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

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

mean(y)

## [1] 2.08851

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

exp(b0e)

## [1] 2.08851

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

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

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

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

  sapply(b0t, function(tmp) sum(log(dgamma(x = y, shape = s, rate = s/exp(tmp)))))

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

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

ここでこの自作関数の引数sは形状母数k (shape parameter)の推定値を入力する。glm関数では形状母数の推定値は最尤推定されていない（！）（Faraway 2006, Extending the Linear Model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models）。fit2内に格納されているdispersion parameterの逆数が形状母数の推定値であり、以下のような値である。

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

shape_e <- 1/fit2$dispersion

shape_e

## [1] 1.229482

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

dispersion parameterは、yの実測値とyの推定値の差をyの推定値で標準化したものの二乗和（残差二乗和にあたる）を自由度（サンプルサイズ - 推定されたパラメータ数）で除したものである。以下のように計算できて、一致することがわかるだろう。

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

fit2$dispersion

## [1] 0.8133505

sum(((y - exp(b0e))/(exp(b0e)))^2)/((length(y) - 1))

## [1] 0.8133505

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

ただし、dispersion parameterの逆数で推定される形状母数は"粗い"と考えられており、MASSパッケージのgamma.shape関数を使えば、最尤推定で形状母数を推定できる。

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

library(MASS)

gamma.shape(fit1)

##                 

## Alpha: 1.0675107

## SE:    0.1336891

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

ここでAlphaが形状母数の推定値である。ネタばらしをすれば、今回のデータは形状母数k = 1としてデータを生成しており、gamma.shape関数のほうがよりよい推定値を与えていることがわかる。しかし、glm関数内ではdispersion parameterの逆数を使った推定値で解析をしているので、ここではそれに従う。

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

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

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

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

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

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

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

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

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

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

実は、ここで私はずっと躓いている。b0が最尤推定値であることはグラフからも正しいのであるが、尤度の計算が合わないのである。b0の最尤推定自体はうまくいっているし、尤度の計算自体も定義にのっとっているだけなので、内部的に形状母数の値の使い方が何やら異なるような予感がしている。というのも、尤度関数を計算していくとわかるが、形状母数はパラメータの推定値に影響を与えない定数項にあたるからである。下記は対数リンクの場合である。尤度関数をパラメータβ0およびβ1の関数とみなせば、尤度関数をβ0およびβ1で偏微分した時に、偏導関数 = 0となるようなβ0およびβ1を求めることに他ならないから、最終的に偏導関数から形状母数は消えてしまい、β0およびβ1の推定に影響を与えないことがわかる。

しかし詳細がわからない。また何かわかれば追記する。

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

以上から、対数尤度関数の頂点付近は分散exp(b0)/(k×(yの総和))の正規分布で近似できることがわかる。それゆえ、b0の標準偏差はこの値の平方根である。実際、下記のように計算すれば上記のGLMの結果である0.09019と一致する。この結果を考えても、形状母数の推定値はやはりdispersion parameterの逆数だと思うのだが。

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

1/sqrt(shape_e*sum(y)*exp(-b0e))

## [1] 0.09018595

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

また、前回同様にWald統計量は推定値を推定値の標準誤差で割ったものである。以前にちょっと紹介したがglm関数ではガンマ分布の時、Wald統計量が残差自由度のt分布に従うとしてp値を計算しているのでポアソン分布や二項分布の時と異なり、以下のように計算できる。

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

b0e*sqrt(shape_e*sum(y)*exp(-b0e))

## [1] 8.165916

sign2 <- function(x){

  if(x < 0){

    return(0)

  } else {

    return(1)

  }

}

pt(b0e*sqrt(shape_e*sum(y)*exp(-b0e)), df = fit1$df.residual, lower.tail = !sign2(b0e)) * 2

## [1] 1.057494e-12

coef(fit2)[1,4]

## [1] 1.057494e-12

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

　逸脱度残差のガンマ分布および対数リンクの場合の計算は下記のように行う。ポアソン分布や二項分布の時と違い、形状母数による補正が入っている。より適切に表現するなら、dispersion parameterを掛けることによる補正である（ポアソン分布や二項分布では1と仮定されており、正規分布のglmでもdispersion parameterを掛ける）。

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

d <- 2*(log(dgamma(x = y, shape = shape_e, rate = shape_e/y)) - log(dgamma(x = y, shape = shape_e, rate = shape_e/exp(b0e))))/shape_e

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

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

## -3.6454588 -0.9384414 -0.3409545  0.3749768  1.8156162

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

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

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

sum(d)

## [1] 107.3321

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

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

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

-2*(f(y, shape_e, b0e) -  sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/y)))) /shape_e

## [1] 107.3321

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

さらにf(y, shape_e, b0e)は（本当なら）最大の対数尤度であるが、求まる値からずれている。一応、この値を-2倍し、推定したパラメータの数×2だけ、値を足したものがAICである。なお、ガンマ分布の場合、b0以外に、形状母数も推定されるパラメータとして計上されているので、2×2 = 4だけ足す。

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

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

## [1] 352.3467

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

　続いて、次のようなデータを考えよう。

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

y <- c(2.53, 3.47, 4.47, 2.65, 2.93, 2.86, 2.72, 1.33, 1.83, 

       2.78, 2.61, 3.13, 3.39, 5.43, 5.33, 2.02, 3.12, 2.36, 

       3.11, 1.1, 3.87, 4.29, 2.62, 1.93, 3.22, 2.8, 3.56, 1.77, 

       3, 1.6, 2.86, 1.92, 4.18, 2.69, 2.27, 2.67, 2.85, 2.25, 

       3.47, 2.25, 4.25, 3.49, 2.95, 3.51, 3.06, 2.81, 5.58, 

       1.5, 2.66, 3.76, 2.81, 2.41, 3.54, 2.24, 2.36, 3.09, 

       2.1, 2.46, 2.17, 3.17, 3.69, 1.89, 4.55, 2.03, 1.33, 

       2.78, 2.38, 2.53, 2.19, 1.89, 2.5, 3.66, 2.77, 4.83,

       1.88, 1.59, 2.48, 1.56, 1.86, 3.01, 4.43, 3.44, 3.41,

       3.63, 2.17, 2.8, 2.91, 3.65, 2.38, 1.33, 1.51, 2.43, 

       2.95, 4.2, 3.32, 3.43, 3.07, 2.43, 2.15, 3.1)

histn(y)#図23の描画

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

直上と同様、データは0より大きいの連続値データである。データは左右非対称な分布で、かつ0から少し離れたところの頻度が大きくなっている。このような分布はk > 1のときのガンマ分布の特徴である。GLMで解析してみると下記のようになる。GLMで解析してみると下記のようになる。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ 1, family = Gamma(link = "log"))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -0.82295  -0.22943  -0.02225   0.18536   0.75501  

## 

## Coefficients:

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

## (Intercept)  1.04837    0.03206    32.7   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 10.251  on 99  degrees of freedom

## Residual deviance: 10.251  on 99  degrees of freedom

## AIC: 261.14

## 

## Number of Fisher Scoring iterations: 4

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

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

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

mean(y)

## [1] 2.853

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

exp(b0e)

## [1] 2.853

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

fit2内に格納されているdispersion parameterの逆数が形状母数の推定値であり、以下のような値である。

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

shape_e <- 1/fit2$dispersion

shape_e

## [1] 9.728722

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

dispersion parameterは以下のように計算できて、一致することがわかるだろう。

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

fit2$dispersion

## [1] 0.1027884

sum(((y - exp(b0e))/(exp(b0e)))^2)/((length(y) - 1))

## [1] 0.1027884

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

MASSパッケージのgamma.shape関数を使えば、最尤推定で形状母数を推定できる。

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

gamma.shape(fit1)

##                

## Alpha: 9.918626

## SE:    1.379759

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

Alphaが形状母数の推定値である。ネタばらしをすれば、今回のデータは形状母数k = 10としてデータを生成しており、gamma.shape関数のほうがよりよい推定値を与えていることがわかる。

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

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

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

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

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

さらに狭い範囲で図示すると次のようになる。

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

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

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

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

対数尤度関数の頂点付近は分散exp(b0)/(k×(yの総和))の正規分布で近似できるので

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

1/sqrt(shape_e*sum(y)*exp(-b0e))

## [1] 0.03206063

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

Wald統計量とp値は以下の通り。

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

b0e*sqrt(shape_e*sum(y)*exp(-b0e))

## [1] 32.69964

sign2 <- function(x){

  if(x < 0){

    return(0)

  } else {

    return(1)

  }

}

pt(b0e*sqrt(shape_e*sum(y)*exp(-b0e)), df = fit1$df.residual, lower.tail = !sign2(b0e)) * 2

## [1] 7.285818e-55

coef(fit2)[1,4]

## [1] 7.285818e-55

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

　逸脱度残差のガンマ分布および対数リンクの場合の計算は下記のように行う。

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

d <- 2*(log(dgamma(x = y, shape = shape_e, rate = shape_e/y)) - log(dgamma(x = y, shape = shape_e, rate = shape_e/exp(b0e))))/shape_e

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

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

## -0.82294587 -0.22943459 -0.02225094  0.18535799  0.75500763

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

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

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

sum(d)

## [1] 10.25128

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

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

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

-2*(f(y, shape_e, b0e) -  sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/y)))) /shape_e

## [1] 10.25128

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

さらにf(y, shape_e, b0e)は（本当なら）最大の対数尤度であるが、求まる値からずれている。一応、この値を-2倍し、推定したパラメータの数×2だけ、値を足したものがAICである。なお、ガンマ分布の場合、b0以外に、形状母数も推定されるパラメータとして計上されているので、2×2 = 4だけ足す。

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

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

## [1] 261.1468

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

　続いて、説明変数xを含む解析を考えよう。次のようなデータを解析する。

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

x <- c(3.2, 24.3, 9.1, 16.8, 7.2, 6.7, 10.8, 18.8, 29.6, 1.9,

       23.2, 17.1, 5.6, 23.2, 4.6, 10.3, 12.9, 4.4, 26.7, 20.3,

       10.6, 14.9, 17.8, 2.6, 28.4, 29.5, 12.6, 18.5, 2.8, 4.8, 

       22.8, 29.7, 21.6, 26.1, 14.6, 22.9, 24.8, 9.1, 26.6, 18.6, 

       3.8, 6.2, 26.2, 7.7, 8.5, 13.5, 20.6, 21.8, 4, 0.9)

y <- c(1.73, 6.79, 0.33, 0.74, 2.7, 0.06, 1.63, 1.34, 6.23, 0.04, 

       0.49, 1.76, 0.25, 0.4, 0.52, 2.3, 1.64, 0.23, 1.52, 6.05,

       0.75, 0.81, 0.8, 0.38, 0.28, 5.17, 0.3, 2.64, 0.45, 2.74, 

       3.07, 10.25, 1.5, 8.18, 2.29, 5.49, 4.26, 0.1, 0.62, 2.77, 

       0.68, 0.84, 0.15, 1.35, 1.33, 1.71, 0.52, 1.15, 0.79, 1.66)

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

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

データを見て、線形回帰ではなくガンマ分布を使ったGLMを想定するべき理由がいくつかある。1. 応答変数yは正の連続値、2. 応答変数yはxが大きくなるにつれてばらつきが大きくなっている、これはガンマ分布の特徴が出ている。では解析をしてみよう。回帰を行い、回帰曲線を描く。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.0875  -0.9815  -0.1579   0.3071   1.4025  

## 

## Coefficients:

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

## (Intercept) -0.43749    0.23093  -1.894   0.0642 .  

## x            0.06396    0.01328   4.817  1.5e-05 ***

## ---

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

## 

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

## 

##     Null deviance: 63.392  on 49  degrees of freedom

## Residual deviance: 46.384  on 48  degrees of freedom

## AIC: 157.39

## 

## Number of Fisher Scoring iterations: 5

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

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

　fit2内に格納されているdispersion parameterの逆数が形状母数の推定値であり、以下のような値である。

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

shape_e <- 1/fit2$dispersion

shape_e

## [1] 1.454732

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

dispersion parameterは以下のように計算できて、一致することがわかるだろう。

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

fit2$dispersion

## [1] 0.6874119

sum(((y - exp(b0e + b1e*x))/(exp(b0e + b1e*x)))^2)/((length(y) - 2))

## [1] 0.6874119

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

MASSパッケージのgamma.shape関数を使えば、最尤推定で形状母数を推定できる。

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

gamma.shape(fit1)

##                 

## Alpha: 1.2176133

## SE:    0.2180724

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

Alphaが形状母数の推定値である。ネタばらしをすれば、今回のデータは形状母数k = 1としてデータを生成しており、gamma.shape関数のほうがよりよい推定値を与えていることがわかる。

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

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

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

  mapply(function(tmp1, tmp2) sum(log(dgamma(x = y, shape = s, rate = s/exp(tmp1+tmp2*x)))), b0t, b1t)

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

db0 <- 0.01

db1 <- 0.001

b0min <- b0e - db0

b0max <- b0e + db0

b1min <- b1e - db1

b1max <- b1e + db1

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

              seq(b1min, b1max, length = 200), 

              f, x = x, y = y)

#対数尤度の計算

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

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

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

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

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

d <- 2*(log(dgamma(x = y, shape = shape_e, rate = shape_e/y)) - log(dgamma(x = y, shape = shape_e, rate = shape_e/exp(b0e+b1e*x))))/shape_e

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

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

## -2.0874958 -0.9814865 -0.1579090  0.3071104  1.4024735

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

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

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

sum(d)

## [1] 46.38362

-2*(f(x, y, shape_e, b0e, b1e) -  sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/y))))/shape_e

## [1] 46.38362

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

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

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

-2*(sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/mean(y)))) -  sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/y))))/shape_e

## [1] 63.392

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

さらにf(x, y, shape_e, b0e, b1e)は（本当なら）最大の対数尤度であるが、求まる値からずれている。一応、この値を-2倍し、推定したパラメータの数×2だけ、値を足したものがAICである。なお、ガンマ分布の場合、b0、b1以外に、形状母数も推定されるパラメータとして計上されているので、2×3 = 6だけ足す。

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

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

## [1] 157.9786

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

最後に、次のようなデータを解析する。

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

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)

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

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

データを見て、直上と同じく、ガンマ分布を使ったGLMを想定するべきだ。では解析をしてみよう。回帰を行い、回帰曲線を描く。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

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

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

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

　fit2内に格納されているdispersion parameterの逆数が形状母数の推定値であり、以下のような値である。

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

shape_e <- 1/fit2$dispersion

shape_e

## [1] 8.289461

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

dispersion parameterは以下のように計算できて、一致することがわかるだろう。

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

fit2$dispersion

## [1] 0.1206351

sum(((y - exp(b0e + b1e*x))/(exp(b0e + b1e*x)))^2)/((length(y) - 2))

## [1] 0.1206351

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

MASSパッケージのgamma.shape関数を使えば、最尤推定で形状母数を推定できる。

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

gamma.shape(fit1)

##                 

## Alpha: 9.372556

## SE:    1.842108

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

Alphaが形状母数の推定値である。ネタばらしをすれば、今回のデータは形状母数k = 10としてデータを生成しており、gamma.shape関数のほうがよりよい推定値を与えていることがわかる。

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

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

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)

#対数尤度の計算

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

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

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

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

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

d <- 2*(log(dgamma(x = y, shape = shape_e, rate = shape_e/y)) - log(dgamma(x = y, shape = shape_e, rate = shape_e/exp(b0e+b1e*x))))/shape_e

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

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

## -0.6812663 -0.2806954 -0.1275530  0.1916897  0.7045343

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

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

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

sum(d)

## [1] 5.429481

-2*(f(x, y, shape_e, b0e, b1e) -  sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/y))))/shape_e

## [1] 5.429481

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

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

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

-2*(sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/mean(y)))) -  sum(log(dgamma(x = y, shape = shape_e, rate = shape_e/y))))/shape_e

## [1] 44.94576

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

さらにf(x, y, shape_e, b0e, b1e)は（本当なら）最大の対数尤度であるが、求まる値からずれている。一応、この値を-2倍し、推定したパラメータの数×2だけ、値を足したものがAICである。なお、ガンマ分布の場合、b0、b1以外に、形状母数も推定されるパラメータとして計上されているので、2×3 = 6だけ足す。

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

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

## [1] 56.3383

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

以上で、正の連続値のようなデータをガンマ分布を使って解析することができた。