0
$\begingroup$

I am trying to determine distribution if this is possible with my data set. After I analyze my data I find that there is possible zero inflated model since there are around 80% of zeros in data set (not missing values). I used data frame with one vector variable M1. I created several models trying to find best fit.

model1p<- glm(M1 ~ .,data= moj.data,family=poisson)
> model1nb<- glm.nb(M1  ~ . , data=moj.data)
> model1zip<-zeroinfl(M1~.|.,data = moj.data, link = "logit", dist = "poisson")
> model1zinb<-zeroinfl(M1~.|.,data = moj.data, link = "logit", dist = "negbin")
> model1plh<-hurdle(M1~.|.,data=moj.data,zero.dist="binomial",link="logit",dist="poisson")
> model1nblh<-hurdle(M1~.|.,data=moj.data, zero.dist="binomial",link="logit",dist="negbin")
> Bloc

kquote

I used even hurdle model (even I know it treats data set as 2 parts) just for comparison. After I ploted this model on hist of my data best alignment is with zero inflated negative binomial. After that i done

Voungov test

> vuong(model1nb,model1zinb)
> 

>Vuong Non-Nested Hypothesis Test-Statistic: 
>(test-statistic is asymptotically distributed N(0,1) under the null that the models are indistinguishible)
-------------------------------------------------------------
              Vuong z-statistic             H_A  p-value
>Raw                  -1.6783428 model2 > model1 0.046640

>AIC-corrected        -1.3532836 model2 > model1 0.087983

>BIC-corrected        -0.7739436 model2 > model1 0.219482

AIC paramters are:

> t(AIC(model1p,model1nb,model1zip,model1zinb,model1plh,model1nblh))

>     model1p model1nb model1zip model1zinb model1plh model1nblh
>df     1.000   2.0000    2.0000     3.0000    2.0000     3.0000

>AIC 1325.297 584.0711  657.4628   575.7447  657.4628   575.7447

Summary of my model is 
> summary(model1zinb)
>
>Call:
>zeroinfl(formula = M1 ~ . | ., data = moj.data, dist = "negbin", link = >"logit")
>
>Pearson residuals:
>
>    Min      1Q  Median      3Q     Max 

>-0.3765 -0.3765 -0.3765 -0.3765  7.8643 
>
>Count model coefficients (negbin with log link):

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

>(Intercept)   1.5215     0.1642   9.266   <2e-16 ***
>Log(theta)    0.1378     0.4122   0.334    0.738    

>Zero-inflation model coefficients (binomial with logit link):
>            Estimate Std. Error z value Pr(>|z|)    
>(Intercept)   1.0491     0.1994   5.262 1.42e-07 ***
>---

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

>Theta = 1.1477 

>Number of iterations in BFGS optimization: 30 

>Log-likelihood: -284.9 on 3 Df

enter image description here I am not expert in this field so I am doing my best :). Question is can I use ZINB model for my data? What should I do next for finding distribution for my dataset?

$\endgroup$

1 Answer 1

1
$\begingroup$

Your regression models just contain intercepts and no real regressors. Hence, the zero-inflated models (zeroinfl) and the hurdle models (hurdle) lead to exactly the same log-likelihoods. They just parameterize the probability for a zero outcome differently.

Judging from the information criteria there is clearly overdispersion (i.e., NB to be preferred over Poisson) and also some amount of excess zeros (ZINB/HNB to be preferred over NB). However, it is conceivable that the picture changes somewhat when including regressors.

For graphically comparing observed and fitted frequencies, the histogram in your figure does not work so well. It is dominated by the probability for zero and deviations at other counts are not visible very well. A better choice would be a rootogram on a square-root scale. This is easily available in the R package countreg (containing the successor functions to hurdle/zeroinfl from pscl). See also Christian Kleiber, Achim Zeileis (2016). "Visualizing Count Data Regressions Using Rootograms." The American Statistician, 70(3), 296--303. http://dx.doi.org/10.1080/00031305.2016.1173590

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.