被説明変数が0以上の値をとるデータのモデリング

　一般化線形モデルGLMの項で紹介したようにこの世のすべてのデータが正規分布に従うと考えることには無理がある。その最たる例が、取れる値に制限があるデータや離散値データである。例えば、種子の数を計測した時に、その値は負になることはないし、（平均のようなパラメータ化した数値の場合を除いて）1.23というような整数の中間値をとることもない。このような被説明変数が-∞から+∞までの値をとりうる正規分布に従うと考えるのは無理があるだろう。そこで被説明変数が従う分布を拡張し、さらにリンク関数によって被説明変数と線形予測子の関係を拡張することで解析する方法がGLMなのであった。本項から数回は、正規分布以外の代表的な確率分布に従う変数をGLMで解析することを目標とする。

ポアソン分布：0以上の離散値をとる確率分布

　種子の数のようなカウントデータを解析するときによく用いられる確率分布がポアソン分布Poisson distributionである。この分布は0以上の離散値（0, 1, 2……）をとる分布である。このような確率分布が定義される数値空間を台supportと呼ぶ。正規分布が平均と分散という2つのパラメータで分布の形が規定されていたように、ポアソン分布にも分布の形を決定するパラメータが存在する。ただし、パラメータは平均（多くの場合、λで表される）、ただ一つである。ポアソン分布の特徴として、平均と分散が同じである特徴がある。それゆえ、平均が唯一のパラメータである。これは言い換えれば、平均が大きくなるにつれて分散も大きくなる分布ということである。この点は、平均と分散が独立に制御される正規分布とは異なる。ポアソン分布を表す確率質量関数Probability mass functionはλをパラメータとして、下記のとおりである。

では、この定義から分布の平均=分散=λであることを示そう。まず、期待値=平均は以下のように計算できる。

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

ゆえにポアソン分布は平均と分散がともにλであることがわかる。では続いて、具体的なデータを生成しながら、ポアソン分布の特徴をつかんでいこう。rpois関数を使うことで、ポアソン分布に従うデータを生成できる。引数nは生成する個数、lambdaは平均値を指定できる。

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

library(plotn)

library(viridis)

rpois(n = 30, lambda = 5)

##  [1]  2 11  5  4  6  8  5  6  4  1  7  7  5  5  4  5  6  3  6  5  4  9  7  6  5

## [26]  5  4  2  5  4

rpois(n = 30, lambda = 0.1)

##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0

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

このように乱数を生成してみるとわかるが、すべて0以上の整数であることがわかるだろう。続いて、平均のパラメータによって確率分布がどうなるかをいくつか確認してみる。試しに、λ=0.5, 1, 1.5を指定してみよう。

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

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

dd$lambda05 <- dpois(dd$X, lambda = 0.5)

dd$lambda10 <- dpois(dd$X, lambda = 1)

dd$lambda15 <- dpois(dd$X, lambda = 1.5)

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

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

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

次はもっと大きなλを指定してみる。

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

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

dd$lambda1 <- dpois(dd$X, lambda = 1)

dd$lambda3 <- dpois(dd$X, lambda = 3)

dd$lambda5 <- dpois(dd$X, lambda = 5)

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

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

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

さらに大きなλを指定してみる。

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

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

dd$lambda1 <- dpois(dd$X, lambda = 1)

dd$lambda10 <- dpois(dd$X, lambda = 10)

dd$lambda20 <- dpois(dd$X, lambda = 20)

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

      leg.lab = c(1, 10, 20), leg.sp = 3, ylab = "確率質量分布")#図3の描画

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

このようにポアソン分布はλの大きさによって分布の平均とともにばらつきも変化していることがわかるだろう。また、λが大きくなってくると分布の形がほぼ左右対称といってもいいような形になる。実際、大きいλに関して言えば、ポアソン分布は正規分布と近似が可能になることが知られている。

ポアソン分布を使ったGLM（ポアソン回帰）

　では、ポアソン分布を使ったGLMについて考察してゆこう。このような回帰分析はポアソン回帰Poisson regressionとも呼ばれる。ポアソン分布を使う場合、Rのプログラム上では、リンク関数はデフォルトで対数関数を使うことになっている。もちろん、対数関数以外をリンク関数に指定可能であるが、ポアソン分布の性質と最も相性がいいのは対数関数であろう。本項ではリンク関数を対数関数にした場合（対数リンク）で議論する。

　対数リンクであるということは下記のように、対数をとることで線形モデルの形に帰着できるということである。逆に元のモデル式はネイピア数の指数部分に線形予測子がのった状態になっている。線形予測子は適宜β・x = β0 + β1×x1 + β2×x2……と置き換えてもらえばよい。 

GLMではこのように確率分布とリンク関数を決めたのち、尤度を最大化するようなパラメータβを求める。なので、次に尤度関数を求めてみよう。データの生成過程を図示して様子をつかむことにする。 

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

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

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

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

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

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

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

b0 <- 0.5

b1 <- 0

n <- 30

x <- runif(n, 0, 30)

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

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

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

データの生成過程通り、xに対して特に関係のなさそうなデータが生成された。では、このデータをGLMで解析する。関数はglm関数を使い、方法はlm関数を用いた線形回帰の時とほとんど同じである。ただし、lmと異なるのは、family引数が存在することである。family引数にはデータが従う確率分布とリンク関数を指定する。ここではポアソン分布と対数リンクを指定するので、poisson(link="log")と指定する。ちなみにlink="identity"とすれば恒等リンク関数（無変換）などのように指定できる。併せて、回帰曲線も描いてみよう。

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

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

fit2 <- summary(fit1)

## 

## Call:

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

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.7634  -0.5868  -0.1026   0.4471   1.9593  

## 

## Coefficients:

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

## (Intercept)  0.56927    0.27750   2.051   0.0402 *

## x           -0.01774    0.01846  -0.961   0.3367  

## ---

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

## 

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

## 

##     Null deviance: 28.602  on 29  degrees of freedom

## Residual deviance: 27.659  on 28  degrees of freedom

## AIC: 89.188

## 

## Number of Fisher Scoring iterations: 5

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

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

plotn(y ~ x)

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

yy <- exp(b0e + b1e * xx)

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

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

GLMの結果をsummary関数を通して得られる結果の解釈も線形回帰の時とそれほど変わらない。今回の場合、線形予測子の切片の推定値は0.57、xの係数の推定値は-0.018であり、初めの設定と大きく変わらないことがわかる。見慣れない変数はDevianceやAICである。Devianceは逸脱度、AICは赤池情報量基準と呼ばれる値である。これら指標値はモデル間比較をするときに重要であるが、詳細は該当ページに譲る。

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

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

b0 <- 0.5

b1 <- 0

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

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

  fit1 <- glm(y ~ x, family = poisson(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] 4201

sum(p2 < 0.05)

## [1] 470

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

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

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

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

b0 <- 0.5

b1 <- 0.05

n <- 30

x <- runif(n, 0, 30)

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

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

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

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

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

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

fit2 <- summary(fit1)

## 

## Call:

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

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -1.85522  -0.52751  -0.09805   0.62741   1.27265  

## 

## Coefficients:

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

## (Intercept)   0.2578     0.2293   1.124    0.261    

## x             0.0591     0.0120   4.926 8.38e-07 ***

## ---

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

## 

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

## 

##     Null deviance: 47.451  on 29  degrees of freedom

## Residual deviance: 22.566  on 28  degrees of freedom

## AIC: 107.01

## 

## Number of Fisher Scoring iterations: 5

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

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

plotn(y ~ x)

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

yy <- exp(b0e + b1e * xx)

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

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

線形予測子の切片の推定値は0.26、xの係数の推定値は0.059であり、これも初めの設定と大きく変わらないことがわかる。回帰曲線は線形回帰の時と異なり、指数関数となっている。

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

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

b0 <- 0.5

b1 <- 0.05

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:10000){

  x <- runif(n, 0, 30)

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

  fit1 <- glm(y ~ x, family = poisson(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] 5567

sum(p2 < 0.05)

## [1] 9862

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

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

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

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

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

y <- c(3, 3, 3, 2, 3, 3, 3, 4, 1, 1, 1, 5, 3, 1, 3, 5, 4, 5,

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

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

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

       2, 3, 3, 1, 1, 1, 0, 1, 4, 1, 4, 5, 0, 4, 2, 0, 8, 4,

       3, 2, 1, 5, 2, 2, 3, 3, 3, 2)

histn(y)#図12の描画

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

データは0以上の離散値データであり、左右対称の分布にも見えないので、ポアソン分布を疑うべきである。一旦は練習なので、説明変数がないデータの解析をしてみることにする。これはデータが正規分布の時であれば、平均が有意に0から異なるかを調べる1標本のt検定に相当する。今回の場合は、切片の係数が0から有意に異なるか、あるいはポアソン分布の平均がexp(0) = 1から有意に異なるかを調べていると考えてもよいだろう。GLMで解析してみると下記のようになる。

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##    Min      1Q  Median      3Q     Max  

## -2.249  -1.097  -0.346   0.287   2.735  

## 

## Coefficients:

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

## (Intercept)  0.92822    0.06287   14.76   <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: 115.76  on 99  degrees of freedom

## Residual deviance: 115.76  on 99  degrees of freedom

## AIC: 374.02

## 

## Number of Fisher Scoring iterations: 5

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

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

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

mean(y)

## [1] 2.53

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

exp(b0e)

## [1] 2.53

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

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

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

f <- function(y, b0t){

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

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

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

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

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

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

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

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

overdraw(abline(v = b0e),

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

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

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

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

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

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

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

overdraw(abline(v = b0e),

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

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

このように対数尤度関数の頂点付近は近寄ればほとんど二次関数のように見える。これはつまり、対数尤度の頂点周りでは二次近似が有効に働くことを示唆しており、二次近似を橋渡しとして、対数をとった正規分布にも近似が可能である。正規分布に近似することでb0の標準偏差を求めることができる。具体的な計算としては、対数をとった正規分布の確率密度関数に関して、二次導関数を求め、平均値における値を求める。以下のように求めると、実は対数をとった正規分布の二次導関数は確率変数の値によらず一定であることがわかる。

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

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

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

1/sqrt(length(y)*exp(b0e))

## [1] 0.06286946

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

さらに、線形回帰と同様にb0の推定値を近似で求めたb0の標準偏差で割った値がsummary内で登場するz valueの値である。この値は計算の方法から検定統計量t（や検定統計量z）と変わらないが、尤度関数を正規分布で近似することで得られる値であり、区別してWald統計量と呼ばれる。Wald統計量は検定統計量tやzの拡張ととらえてもよいだろう。プログラムによって使われる帰無分布が異なるようでありglm関数では、ポアソン分布や二項分布では標準正規分布（ゆえにWald統計量がzと表示）、ガンマ分布では残差自由度を持つt分布（ゆえにWald統計量がtと表示）だったりする。他にはglmmTMB関数ではすべて標準正規分布（ゆえにz）で統一されている。今回のパラメータ推定値のp値はこのWald統計量が標準正規分布に従うことを利用した検定、Wald検定である。したがって、Wald統計量とそれを利用したp値は以下のように求まる。summaryでは省略されているが、具体的に中身を取り出すとp値がよく一致していることがわかる。

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

b0e*sqrt(length(y)*exp(b0e))

## [1] 14.76423

sign2 <- function(x){

  if(x < 0){

    return(0)

  } else {

    return(1)

  }

}

pnorm(b0e*sqrt(length(y)*exp(b0e)), mean = 0, sd = 1, lower.tail = !sign2(b0e)) * 2

## [1] 2.491447e-49

coef(fit2)[1,4]

## [1] 2.491442e-49

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

上記のWald検定の他に、説明変数が有意な効果を持つか判定する検定として説明変数を加えることによって尤度が有意に改善したかを判断する尤度比検定があげられる。こちらは別の項で解説しよう。

　またここまで説明してなかったが、今まで残差Residualsとよばれていた項目はGLMでは逸脱度残差Deviance residualsというものに名前が変わっている。これは逸脱度の計算に基づく「残差のようなもの」であり、かなり込み入った話になるので詳細はそのページおよび追記に譲るとして、ポアソン分布および対数リンクの場合の計算は下記のように行う。線形回帰で求めてきた残差と計算が違うことさえ覚えておけば、とりあえずは混乱を避けられるだろう。

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

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

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

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

## -2.2494444 -1.0970695 -0.3459935  0.2869807  2.7348777

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

ここで求めたdの総和は残差逸脱度Residual devianceと呼ばれる。今回のように説明変数を持たないモデルの場合、このモデルはヌルモデルNull modelと呼ばれ、残差逸脱度はヌル逸脱度Null devianceと一致する。

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

sum(d)

## [1] 115.7593

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

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

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

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

## [1] 115.7593

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

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

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

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

## [1] 374.0246

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

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

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

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)

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

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

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

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

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

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

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

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

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

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

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

  mapply(function(tmp1, tmp2) {

    sum(log(dpois(x = y, lambda = 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)#図17の図示

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

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

abline(v = b0e)

abline(h = b1e)

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

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

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

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

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

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

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

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

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

　説明変数が追加されると、もはや容易には係数の標準偏差を求めることはできない。図17を見ればわかるようにb0とb1の相関を考慮しながら多変量正規分布に近似をする必要があるだろう（たぶん、図17の最尤推定値周りの二次導関数を表しているのがヘッセ行列とかいうものだと思う。ヘッセ行列の逆行列が多変量正規分布の分散共分散行列になるらしい）。逸脱度残差は計算できて、次のようにする。

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

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

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

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

## -1.4696020 -0.9610454 -0.4634509  0.4594230  2.3777741

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

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

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

sum(d)

## [1] 53.30222

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

## [1] 53.30222

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

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

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

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

sum(dn)

## [1] 124.4836

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

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

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

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

## [1] 157.6728

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

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

追記：Full model、Null modelと逸脱度

　上記本文中では、逸脱度、残差逸脱度、逸脱度残差について、適当にごまかしてきた。というのも、この話をするとかなり込み入った話になってきて、推定値の話から離れてしまうからである。なので、セクションを分けてこちらで解説することにする。といっても詳しい解説はAICのページに譲り、ここでは最低限のことだけを解説する。

　まず、あるデータセット{xi, yi}が与えられ、それをモデル化するときのフルモデルFull modelとヌルモデルNull modelを理解する必要がある。もし、データペアxi, yiがn個あるとして、このデータを完全に再現するために必要なパラメータはいくつだろうか。例えば、(xi, yi) = (1,2), (3,3), (6,8)みたいなデータならどうだろう。この時、x1 = 1のときy1 = 2、x2 = 3のときy2 = 3、x3 = 6のときy3 = 8を返すようなモデルを組めば、データを完全に再現できる。パラメータをbiとでも置いて、yi = biとなるbi、つまりb1 = 2、b2 = 3、b3 = 8というように3つのパラメータを使えば元のデータを完全復元できるだろう。

　一般に、サンプルサイズnのデータに対してn個のパラメータを用意すればデータを完全再現できる。このようなデータを完全再現してしまうモデルをフルモデルとよぶ。この議論あるいは過去の議論（逆問題とメカニズムの想起）を見ればわかるように、フルモデルは元のデータセットを復唱しているだけに過ぎず、データをパラメータによって要約していないし、予測性もない。あるいは、こちらで議論したように、パラメータが過剰にあるゆえに過学習しているともいえるだろう。一方で、xとyには何の関係もないと考え、つまりどんなxに対しても平均的なyを返すようなモデルを考える。この時、yの平均 = 7.5という1つのパラメータだけでモデルを記述する。このようなモデルのことをヌルモデルと呼ぶ。この場合は、こちらでいうならば未学習な状態、と言えるだろう。

　さて、フルモデルとヌルモデルにおける最大対数尤度を求めてみよう。例えば、ポアソン分布を仮定した時どうなるかを考える。最大対数尤度はあるyiを最も高い確率で返すような確率分布を決定することである。yiの1つにつき、一つのパラメータbiを与えるのがフルモデルだったから、yiを平均とするようなポアソン分布を仮定すればよい。つまり、以下の計算を行う。

一方ヌルモデルはどんなxiにも平均E(y)を返すのだから、その最大対数尤度は以下のようになる。

これらの計算を終えることができれば、残差逸脱度の計算は簡単に計算できる。今、自分の手元で組んだ統計モデルの最大対数尤度をlogLmodelとすれば、残差逸脱度Dは、

すなわち、残差逸脱度とはフルモデルからどれくらい尤度が離れているかを示した指標である。残差逸脱度が小さければ、よりよくデータセットを説明できているモデルであると言えるだろう。つまりは、この残差逸脱度は最小二乗法における残差平方和に相当するものであると考えてもよい。残差二乗和のより広いモデルに適応したものが、残差逸脱度である。ヌル逸脱度は、logLmodelをlogLnullに置き換えたものである。

　この残差逸脱度にも、以前議論したような問題が同様に浮上する。つまり、残差逸脱度が小さいモデルほど、予測性の良いモデルとは限らないということである。実際、フルモデルはただの数値の復唱に過ぎず、予測とは縁遠い。しかし、説明変数を多くすれば多くするほど、残差逸脱度は小さくなる傾向がある。そこで、罰則項最小二乗法のように、残差逸脱度にもパラメータの数に対する罰則項をつけることで変数選択を目指す方法が、AICなのである。

　試しに、以下の問題を改めて考えて、残差逸脱度を計算してみる。

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

y <- c(3, 3, 3, 2, 3, 3, 3, 4, 1, 1, 1, 5, 3, 1, 3, 5, 4, 5,

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

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

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

       2, 3, 3, 1, 1, 1, 0, 1, 4, 1, 4, 5, 0, 4, 2, 0, 8, 4,

       3, 2, 1, 5, 2, 2, 3, 3, 3, 2)

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

fit2 <- summary(fit1)

fit2

## 

## Call:

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

## 

## Deviance Residuals: 

##    Min      1Q  Median      3Q     Max  

## -2.249  -1.097  -0.346   0.287   2.735  

## 

## Coefficients:

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

## (Intercept)  0.92822    0.06287   14.76   <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: 115.76  on 99  degrees of freedom

## Residual deviance: 115.76  on 99  degrees of freedom

## AIC: 374.02

## 

## Number of Fisher Scoring iterations: 5

f <- function(b0t){

  sapply(b0t, function(tmp){

    tmp * sum(y) - length(y)*exp(tmp) - sum(log(factorial(y)))

    }

    )

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

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

## [1] 115.7593

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

f(b0e)はこのモデルにおける最大対数尤度を計算しており、後半のsum(log(dpois(y, lambda = y)))によってフルモデルの最大対数尤度を計算している。