hurdle model

昨日紹介した遠洋水研外洋資源部鯨類管理研究室の岡村寛さんのありがたいＰＤＦにhurdle modelというのも紹介されています。

遠洋水産研究所といえば、資源選択検定のレビューとＲコードを書いてくれた人達ではないですか。
これもかなりありがたいです。

話をもとに戻して、
Hurdle modelもpsclパッケージに扱われています。

以前KuboLogにて0の切れたポアソン分布の話題が出ていました。

これもHurdle modelの一種のような気がします。

KuboLogではVGAMパッケージのvglm(,family=pospoisson)で解いているが、

psclパッケージならhurdle()で解けるのかな？
hurdle()は、0の部分とそれ以外の部分を別々に推定してくれるので。

というわけでKuboLogコードをパクって検証。

v.sample <- rpois(1000, exp(-1))
y <- v.sample[v.sample > 0]
cat("# length(y) =", length(y), ", mean(y) = ", mean(y), "\n")
# length(y) = 323 , mean(y) =  1.235294
fit1 <- vglm(y ~ 1, family = pospoisson())
fit2 <- hurdle(y ~ 1)
以下にエラー hurdle(y ~ 1) :
  invalid dependent variable, minimum count is not zero
fit3 <- hurdle(v.sample ~ 1)
summary(fit1)

Call:
vglm(formula = y ~ 1, family = pospoisson())

Pearson Residuals:
                 Min       1Q   Median       3Q    Max
log(lambda) -0.46948 -0.46948 -0.46948 -0.46948 7.5117

Coefficients:
               Value Std. Error t value
(Intercept) -0.82411    0.11102  -7.423

Number of linear predictors:  1

Name of linear predictor: log(lambda)

Dispersion Parameter for pospoisson family:   1

Log-likelihood: -136.0555 on 322 degrees of freedom

Number of Iterations: 5

summary(fit3)

Call:
hurdle(formula = v.sample ~ 1)

Count model coefficients (truncated poisson with log link):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.8241     0.1110  -7.423 1.15e-13 ***
Zero hurdle model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.74002    0.06762  -10.94   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 6
Log-likelihood:  -823 on 2 Df

ごめんなさい。まったくの別物でした。

で、ここまで書いて、統計学なんでもあり掲示板の５２７２番からその一連のレスを読むともうすでにこのことについて議論されたあとでしたね。

完全に恥さらしなので、全部消そうかと思ったけど、ここまで書いてしまったのでこのままアップすることにします。

で、Hurdle modelっていったい何なんだろう？

０とそれ以外を別個に計算するのはZero Inflated modelも同じなわけだし。

雰囲気としてHurdle modelとZero Inflated modelはお仲間で計算方法が違うものという認識でいいのかな。

overdispersionが0に引っ張られるときに、0がすごくたくさんある時は、Zero Inflated Modelで、そうでなさそなときはHurdle Modelということか？