廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

optim による glm:最尤推定 + 信頼区間, Wald 検定

あらまし

  1. 自分で尤度を書いてみて, R の glm 関数がやってることを再現する.
  2. 尤度さえ書ければパラメータの点推定, 区間推定ができるし検定もできる.
  3. それができるようになれば, パッケージなどが用意されていない新しいモデルでも計算できるようになる(と思う).

最尤推定

glm そのものについては、久保『データ解析のための統計モデリング入門』などを参照して欲しい.


今回はポアソン回帰を例にとる.

目的変数 y_i が平均 \log \lambda_i= \beta_0 +\beta_1 x_iポアソン分布に従うとすると, 確率関数は,

\displaystyle f(y_i) = \frac{\lambda_i^{y_i} e^{-\lambda_i}}{y_i !}

である.

データは ポアソン回帰(カウントデータ) | 一般化線形モデル より.

このモデルのパラメータ \beta_0, \beta_1 を推定するには以下のようにすればよい.

x<- 1:12
y <- c(6, 5, 14, 14, 21, 31, 38, 51, 86, 136, 169, 234)
fit1 <- glm(y ~ x, family = poisson(link = "log"))

これがなにをやっているのか理解するために, 自分で対数尤度関数を書いてみる.

対数尤度関数は

\displaystyle \log \prod f(y_i) = \sum \log f(y_i)

だった.

\log \lambda_i= \beta_0 x_i+\beta_1\lambda_i= \exp(\beta_0 x_i+\beta_1) なので, こんな感じに書けばよい.

LL <- function(par){
  beta0 <- par[1]
  beta1 <- par[2]
  sum(dpois(y,exp(beta0+beta1*x),log = TRUE))
}

尤度関数を最大にする \beta_0, \beta_1最尤推定量だった.

R には optim という関数を最小化してくれる関数がある.

デフォルトだと最小化なので control = list(fnscale=-1) とすると最大化してくれる.

opt1 <- optim(c(1,1),LL,control = list(fnscale=-1))
opt2 <- optim(opt1$par,LL,control = list(fnscale=-1),method = "BFGS",hessian = TRUE)

ここでは一度ネルダーミード法で求めた値を初期値にして, BFGS でパラメータを求めることした.

ぼくの経験ではこうするとうまくいくことが多い.

最適化の手法の特徴については, 【 R言語 】最適化問題 でよく使う optim( )関数・optimize( )関数 まとめ - Qiita などを参照して欲しい.

> coef(fit1)
(Intercept)           x 
  1.2989403   0.3487019 
> opt2$par
[1] 1.2992184 0.3486741

だいたい glm の結果と近い値が求まった.

せっかくなのでプロット.

f:id:abrahamcow:20160409085408p:plain

補遺

  • なぜ \lambda_i= \beta_0+\beta_1x_i としないで, \log \lambda_i= \beta_0 +\beta_1 x_i とするの?
    1. \lambda_i= \beta_0+\beta_1x_i としてもいい.
    2. \lambda_i は正の値しかとれない. そういう制約をつけるよりも \log \lambda_i として負の値もとれるようにしたほうが最適化がかんたん.
    3. 指数型分布族のナチュラルパラメータは一般化線形モデルの自然リンク関数である.

ああああ: リンク:GLM

信頼区間

glm によるパラメータの推定結果を改めて見てみる.

> coef(summary(fit1))
             Estimate Std. Error   z value      Pr(>|z|)
(Intercept) 1.2989403 0.14658401  8.861405  7.901234e-19
x           0.3487019 0.01453532 23.989967 3.539181e-127

Estimate, Std. Error, z value, Pr(>|z|) と四つの値が出力されている.

順番に見ていこう.

Estimate はパラメータの点推定値でさっき求めた.

Std. Error は標準誤差(推定量の標準偏差)だ.

ヘシアン(対数尤度関数の二回微分)を計算することで、その負の逆行列で近似的に分散共分散行列を求めることができる.

optim のオプションで hessian = TRUE を選ぶと、近似ヘシアンを求めてくれる.

対角成分が分散なので,

> (my.Std <-sqrt(-diag(solve(opt2$hessian))))
[1] 0.14657254 0.01453409

glm の出力と近い値が求まった.

最尤推定量は漸近的に正規分布に従うので, 標準偏差がわかれば信頼区間を得ることができる.

95%信頼区間を出すには以下のようにする.

> up <- qnorm(0.975,opt2$par,my.Std)
> low <- qnorm(0.025,opt2$par,my.Std)
> round(matrix(c(low,up),2,2),2)
     [,1] [,2]
[1,] 1.01 1.59
[2,] 0.32 0.38
> round(confint(fit1),2) #glm で確かめ算
Waiting for profiling to be done...
            2.5 % 97.5 %
(Intercept)  1.01   1.58
x            0.32   0.38

Wald(ワルド)検定

最後に係数の推定値の表にあらわれた z value と Pr(>|z|) について.

z はパラメータの推定値を標準誤差で割った値で、ワルド検定の検定統計量になる.

最尤推定量は漸近的に正規分布に従うので、標準誤差で割った値は分散 1 の正規分布に従う.

帰無仮説を “回帰係数は 0” とすると、帰無仮説の下で、z は平均 0、分散 1 の正規分布に従う.

その二乗は自由度 1 のカイ二乗分布に従う.

このことを利用して検定ができる.

> (z <- opt2$par/my.Std)
[1]  8.863996 23.990082
#Wald test
pchisq((opt2$par/my.Std)^2,df=1,lower.tail=FALSE)
[1]  7.719623e-19 3.529427e-127

これも大体 glm と同じ値になった.