TMBによるポアソン回帰モデルとゼロ過剰ポアソンモデルのあてはめ

Applied Statistical Modelling for Ecologistsで使われているということもありましたので、TMBをためしてみました。

準備

TMBライブラリと、データ可視化用にggplot2を読み込みます。擬似乱数も固定します。

library(TMB)
library(ggplot2)
set.seed(123)

データ

ゼロ過剰ポアソン分布にしたがう模擬データを生成します。

N <- 100
p <- 0.3
alpha <- -1
beta <- 0.8
x <- runif(N, 0, 5)
y <- rbinom(100, 1, 1 - p) * rpois(N, exp(alpha + beta * x))
example <- data.frame(Y = y, X = x)

このようなデータになりました。

ggplot(example, aes(x = X, y = Y)) +
  geom_point()

ポアソン回帰

まずはポアソン回帰モデルをあてはめます。

モデル

例を参考にして、以下のようなモデルを書きました。対数尤度×-1を返り値にするみたいです。このあたりはStanでtarget構文をつかう場合にちょっと似ている気がしました。あとで、Rで対数尤度×-1を最小化します。

TMBのC++のコードです。"models/poisson.cpp"に保存しておきます。

poisson.cpp

// Poisson regression
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_VECTOR(Y);
  DATA_VECTOR(X);
  PARAMETER(alpha);
  PARAMETER(beta);
  Type lp = 0;
  vector<Type> lambda = exp(alpha + beta * X);
  lp += -sum(dpois(Y, lambda, true));
  return lp;
}

コンパイルと最適化

モデルをコンパイルして、できたライブラリをロードします。

モデルをコンパイルしなおしても挙動がかわらなくて困ったことがあったのですが、そういうときはいったんライブラリをアンロードするとよいのかもしれません（Rを再起動してました）。

model_name <- "poisson"
file.path("models", paste(model_name, "cpp", sep = ".")) |>
  compile()
file.path("models", dynlib(model_name)) |>
  dyn.load()

MakeADFun関数で、最適化関数に渡すオブジェクトを作成して、最適化します。最適化関数にはnlminbを使用しました。Introductionではoptim関数を使用していましたので、そちらもつかえるみたいです。

data <- list(Y = example$Y, X = example$X)
parameters <- list(alpha = 1, beta = 1)
obj_pois <- MakeADFun(data, parameters, DLL = model_name)
## Constructing atomic D_lgamma
opt_pois <- nlminb(obj_pois$par, obj_pois$fn, obj_pois$gr)
## outer mgc:  29879.53 
## outer mgc:  771.6844 
## outer mgc:  796.7531 
## outer mgc:  183.2513 
## outer mgc:  39.24369 
## outer mgc:  22.22143 
## outer mgc:  3.306709 
## outer mgc:  0.8717667 
## outer mgc:  0.09010732 
## outer mgc:  0.01217738 
## outer mgc:  9.899208e-08

結果

結果です。データにあわないモデルなので、推定値はまあそれなりです。

print(opt_pois)
## $par
##     alpha      beta 
## -1.367507  0.815415 
## 
## $objective
## [1] 247.2812
## 
## $convergence
## [1] 0
## 
## $iterations
## [1] 10
## 
## $evaluations
## function gradient 
##       14       11 
## 
## $message
## [1] "relative convergence (4)"

glm関数による結果と比較してみます。

fit_glm <- glm(Y ~ X, family = poisson, data = example)
summary(fit_glm)
## 
## Call:
## glm(formula = Y ~ X, family = poisson, data = example)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.36751    0.20423  -6.696 2.14e-11 ***
## X            0.81542    0.05152  15.827  < 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: 646.12  on 99  degrees of freedom
## Residual deviance: 300.58  on 98  degrees of freedom
## AIC: 498.56
## 
## Number of Fisher Scoring iterations: 5
logLik(fit_glm)
## 'log Lik.' -247.2812 (df=2)

パラメータ・対数尤度とも、同じ値になっていました。

ただし、残差逸脱度が自由度の3倍くらいあるので、やはり過分散になっているようです。

ゼロ過剰ポアソンモデル

つづいて、同じデータをゼロ過剰ポアソンモデルにあてはめます。

モデル

dzipoisという分布関数が用意されていますので、これを使用します。

zipois.cpp

// Zero-Inflated Poisson
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_VECTOR(Y);
  DATA_VECTOR(X);
  PARAMETER(p);
  PARAMETER(alpha);
  PARAMETER(beta);
  Type lp = 0;
  vector<Type> lambda = exp(alpha + beta * X);
  lp += -sum(dzipois(Y, lambda, p, true));
  return lp;
}

コンパイルと最適化

モデルをコンパイルして、できたライブラリをロードします。

model_name <- "zipois"
file.path("models", paste(model_name, "cpp", sep = ".")) |>
  compile()
file.path("models", dynlib(model_name)) |>
  dyn.load()

最適化します。

data <- list(Y = example$Y, X = example$X)
parameters <- list(p = 0.5, alpha = 1, beta = 1)
obj_zip <- MakeADFun(data, parameters, DLL = model_name)
## Constructing atomic D_lgamma
opt_zip <- nlminb(obj_zip$par, obj_zip$fn, obj_zip$gr)
## outer mgc:  21261.9 
## outer mgc:  939.6188
## Warning in nlminb(obj_zip$par, obj_zip$fn, obj_zip$gr): NA/NaN function
## evaluation
## outer mgc:  873.141 
## outer mgc:  364.9433 
## outer mgc:  365.5757 
## outer mgc:  70.04271 
## outer mgc:  69.5802
## Warning in nlminb(obj_zip$par, obj_zip$fn, obj_zip$gr): NA/NaN function
## evaluation
## outer mgc:  80.84804 
## outer mgc:  23.68681 
## outer mgc:  25.78416 
## outer mgc:  92.0502 
## outer mgc:  17.66233 
## outer mgc:  56.65431 
## outer mgc:  10.246 
## outer mgc:  2.013517 
## outer mgc:  0.09686809 
## outer mgc:  0.01954563 
## outer mgc:  7.131502e-05

結果

結果です。最適化で警告が出たりしましたが、convergence=0なので収束したようです。

推定値も、データを生成した元の値に近いものになっていました。

print(opt_zip)
## $par
##          p      alpha       beta 
##  0.2509950 -1.1137171  0.8300888 
## 
## $objective
## [1] 174.6459
## 
## $convergence
## [1] 0
## 
## $iterations
## [1] 18
## 
## $evaluations
## function gradient 
##       23       18 
## 
## $message
## [1] "relative convergence (4)"