1
$\begingroup$

I am working on analyzing panel data of countries with several independent variables. I am aware from previous literature on panel data that an OLS model could be performed. However, because the estimation results for count data will suffer bias in OLS where dependent values are allowed to take both negative and positive values, I chose the Poisson model. I have dependent variables that are of count data, over dispersed, and having excess of zeroes.

I performed a count data hurdle regression with Poisson and to prevent multicollinearity, variance inflation factors are checked to arrive at my optimal model. Every single output that I got (be it univariate or multivariate analysis) had narrow confidence intervals, low coefficents and extremely low p values (< 2.2e-16). Here is my final model -

    Call:
hurdle(formula = Y ~ A+B+C+ 
    offset(log(Pop.Density)) | 1, data = Data, dist = "poisson")

Pearson residuals:
     Min       1Q   Median       3Q      Max 
 -0.7931  -0.7807  -0.7401   0.2798 301.3965 

Count model coefficients (truncated poisson with log link):

                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -0.694339   0.004247 -163.47   <2e-16 ***
A                        -0.396599   0.005468  -72.53   <2e-16 ***
B                        -0.328605   0.004792  -68.58   <2e-16 ***
C                        -0.240072   0.004170  -57.58   <2e-16 ***

Zero hurdle model coefficients (binomial with logit link):

            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4635     0.0393  -11.79   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 18 
Log-likelihood: -1.381e+05 on 5 Df
> vif(Mv7)
                   A                         B                       C
                3.564376                 2.419123                 3.663317

My primary question is, what could be causing these p values to be so low in my output? My dataset is not extremely huge, but I would say fairly large (3083 obs of 15 variables). A, B,C have values that range from negative to positive.

Could there be a time effect that effect my results that I did not consider? Thanks!

enter image description here

> Call:
hurdle(formula = Y ~ A+B+C+ offset(log(Pop.Density)) | 1, data = Data, dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.6743 -0.3575 -0.3489 -0.2177 37.8530 

Count model coefficients (truncated negbin with log link):
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               0.67061    0.06041  11.100  < 2e-16 ***
A                         0.16976    0.10964   1.548   0.1215    
B                        -0.43295    0.09676  -4.474 7.66e-06 ***
C                        -0.13336    0.07475  -1.784   0.0744 .  
Log(theta)               -1.06590    0.06579 -16.201  < 2e-16 ***
Zero hurdle model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4635     0.0393  -11.79   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 0.3444
Number of iterations in BFGS optimization: 25 
Log-likelihood: -8014 on 6 Df
> vif(Mv7)
             A                    B                         C
       5.387540                 3.981045                 2.848343 

When I use population total as an offset, I get an error in R - Error: no valid set of coefficients has been found: please supply starting values

When I use log(Population Total) as an offset this is what I get -

> 
Call:
hurdle(formula = Y ~ A+B+C+
    offset(log(Pop.Total)), data = Data, dist = "negbin")

Pearson residuals:
      Min        1Q    Median        3Q       Max 
 -0.64145  -0.40354  -0.23862  -0.07396 250.78934 

Count model coefficients (truncated negbin with log link):
                          Estimate Std. Error  z value Pr(>|z|)    
(Intercept)              -11.98678    0.05473 -219.016  < 2e-16 ***
A                        -0.33166    0.09544   -3.475 0.000511 ***
B                        -0.22369    0.09157   -2.443 0.014567 *  
C                         0.55690    0.07117    7.825 5.07e-15 ***
Log(theta)               -0.84623    0.06062  -13.960  < 2e-16 ***
Zero hurdle model coefficients (binomial with logit link):
                          Estimate Std. Error  z value Pr(>|z|)    
(Intercept)              -16.44459    0.05300 -310.294  < 2e-16 ***
A                         -0.38690    0.09239   -4.188 2.82e-05 ***
B                          0.12640    0.08411    1.503   0.1329    
C                         -0.13685    0.07855   -1.742   0.0815 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 0.429
Number of iterations in BFGS optimization: 43 
Log-likelihood: -7389 on 9 Df
> AIC(Mv8)
[1] 14795.51
> vif(Mv8)
          A                        B                        C
       3.815796                 3.371319                 3.406296 

This is my outcome without using an offset.

Call:
hurdle(formula = Y~A+B+C, 
    data = Data, dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.8101 -0.4350 -0.3287 -0.1220 21.8952 

Count model coefficients (truncated negbin with log link):
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               4.27904    0.04229 101.191  < 2e-16 ***
 A                       -0.01475    0.07926  -0.186   0.8524    
 B                       -0.11324    0.05358  -2.113   0.0346 *  
 C                       -0.49852    0.06001  -8.308  < 2e-16 ***
Log(theta)               -0.28185    0.04783  -5.893  3.8e-09 ***
Zero hurdle model coefficients (binomial with logit link):
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -0.55545    0.04305 -12.904   <2e-16 ***
A                         0.09582    0.07769   1.233    0.217    
B                         0.06450    0.06931   0.931    0.352    
C                        -0.97937    0.06872 -14.251   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 0.7544
Number of iterations in BFGS optimization: 18 
Log-likelihood: -7558 on 9 Df
> AIC(Mv7)
[1] 15133.16
> vif(Mv7)
                     A                        B                        C
                3.683669                 2.178420                 2.897415 

Offset with population total scaled to 1. ( I divided each country's population total with largest population) -

Call:
hurdle(formula = Y~A+B+C+ 
    offset(Pop.Total.Test), data = POP, dist = "negbin", link = "logit")

Pearson residuals:
     Min       1Q   Median       3Q      Max 
-0.88522 -0.45546 -0.34142 -0.06619 18.90660 

Count model coefficients (truncated negbin with log link):
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               4.15278    0.03762 110.386   <2e-16 ***
A                        -0.11584    0.06897  -1.680    0.093 .  
B                         0.02039    0.05119   0.398    0.690    
C                        -0.45935    0.05182  -8.864   <2e-16 ***
Log(theta)               -0.07232    0.04557  -1.587    0.113    
Zero hurdle model coefficients (binomial with logit link):
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -0.58007    0.04310 -13.458   <2e-16 ***
A                         0.08308    0.07774   1.069    0.285    
B                         0.06679    0.06962   0.959    0.337    
C                        -0.95615    0.06867 -13.923   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 0.9302
Number of iterations in BFGS optimization: 13 
Log-likelihood: -7433 on 9 Df

> AIC(Mv7)
[1] 14884.03

For my second dependent variable (occurrence) - My optimal model was A+B+log(pop.density). To calculate occurrence rate, I offset it with log(pop.total). The code follows below.

Call:
hurdle(formula = Occurence ~ A +B+ log(Pop.Density) + 
    offset(log(Pop.Total)), data = Data, dist = "negbin", link = "logit")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.1211 -0.5798 -0.2758  0.0748 20.7457 

Count model coefficients (truncated negbin with log link):
                       Estimate Std. Error  z value Pr(>|z|)    
(Intercept)           -16.08087    0.15542 -103.467  < 2e-16 ***
A                     -0.35991    0.07428   -4.846 1.26e-06 ***
B                     -0.13195    0.06224   -2.120  0.03399 *  
log(Pop.Density)       -0.21688    0.03435   -6.314 2.73e-10 ***
Log(theta)              0.38368    0.12342    3.109  0.00188 ** 
Zero hurdle model coefficients (binomial with logit link):
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -16.02224    0.16485 -97.194  < 2e-16 ***
A                     -0.28540    0.07408  -3.852 0.000117 ***
B                     -0.09963    0.07744  -1.287 0.198255    
log(Pop.Density)       -0.08668    0.03829  -2.264 0.023590 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 1.4677
Number of iterations in BFGS optimization: 12 
Log-likelihood: -3090 on 9 Df
> AIC(Mvfrequency)
[1] 6198.895
> vif(Mvfrequency)
                 A                     B      log(Pop.Density) 
             2.776789              3.469304             11.629956 

The output when I remove population density as a factor.

Call:
hurdle(formula = Occurence ~ A+B + offset(log(Pop.Total)), 
    data = Data, dist = "negbin", link = "logit")

Pearson residuals:
     Min       1Q   Median       3Q      Max 
-1.04306 -0.57164 -0.26873  0.05331 18.46993 

Count model coefficients (truncated negbin with log link):
                       Estimate Std. Error  z value Pr(>|z|)    
(Intercept)           -17.02648    0.06424 -265.048  < 2e-16 ***
A                     -0.38514    0.07957   -4.840  1.3e-06 ***
B                     -0.13069    0.06732   -1.941   0.0522 .  
Log(theta)              0.14126    0.12476    1.132   0.2575    
Zero hurdle model coefficients (binomial with logit link):
                       Estimate Std. Error  z value Pr(>|z|)    
(Intercept)           -16.38060    0.05235 -312.893  < 2e-16 ***
A                     -0.31557    0.07265   -4.344  1.4e-05 ***
B                     -0.07864    0.07668   -1.026    0.305    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 1.1517
Number of iterations in BFGS optimization: 10 
Log-likelihood: -3115 on 7 Df
> AIC(Mvfrequency)
[1] 6243.296
> vif(Mvfrequency)
               A                     B 
             2.694534              3.328758 
$\endgroup$
  • $\begingroup$ Are you referring to both parts of the model or is this just for the count part or just for the hurdle part? It might help to post output from a simple model with just one predictor. $\endgroup$ – mdewey Jul 30 '16 at 15:32
  • $\begingroup$ You seem to have (at least) one outlying observation. The largest residual looks suspicious to me. Perhaps you should check that? $\endgroup$ – mdewey Jul 30 '16 at 17:09
  • $\begingroup$ You are right @mdewey. There are outliers in my data. Here is a screenshot of my plot of my Y (Y axis) with A (X axis) - I added it to my post. I see these outliers. I don't know how to define these outliers so I can filter them to fit in my model. $\endgroup$ – Sofia Jul 30 '16 at 19:26
  • $\begingroup$ First thing is to see whether this is just over dispersion. Try using negative binomial rathe than Poisson. $\endgroup$ – mdewey Jul 30 '16 at 20:59
  • $\begingroup$ @mdwey, I performed a negative binomial and the results are included in my post. 1) I did not remove the outliers 2) what does negative log(theta) mean? I have already read the post on theta on this forum (stats.stackexchange.com/questions/10419/…) - I have read Mr. Joesph Hilbe's detailed comment but I still don't understand what a negative theta means. Could you please succinctly explain the significance of a negative log of theta? $\endgroup$ – Sofia Jul 30 '16 at 21:55
1
$\begingroup$

The problem here is that the distribution of the counts does not follow a Poisson. This was evidenced by the very large residuals. Fitting a negative binomial improves that aspect but as is shown by the very small value of $\theta$ the counts are still very concentrated in the low values but with a few large ones. Investigating these large counts would prove instructive. They may be (a) errors (b) values from another process which should have been excluded (c) really, really interesting values which cast light on aspects of the process hitherto unknown. Only careful examination will reveal which is which.

To get a feel for how $\theta$ affects the shape try plotting using dnbinomp for a few example values of $\theta$ and see what happens.

The role of the offset may need some attention too. Recall that an offset is simply a covariate whose coefficient is known to be 1. The most common example is in Poisson regression where a count is being modelled via a log link and the denominator is included as its log. So if we model number of traffic accidents in cities we use log(citysize) as an offset since we believe that if he city is twice as large it will have twice as many accidents. This is a testable assumption, we can include log(citysize) as a covariate in addition to as an offset and if its coefficent as a covariate is different from zero it means our simple assumption about the relationship is wrong.

$\endgroup$
  • $\begingroup$ Thank you @mdewey! This is quite helpful. The large counts in my dataset are (c ) interesting values which originate from highly populated countries. What would you recommend to improve my model in this case? The options I see are to (1) add another variable to explain the large values but if I use population total, I have a correlation issue (2) To transform my dependent variable by applying a log function. However, because of the zeroes, this doesn't seem feasible. (3) Remove outliers by inserting a threshold. Is there a way (AIC?) to measure which is the best distribution - poisson or NB? $\endgroup$ – Sofia Aug 1 '16 at 10:47
  • $\begingroup$ @Sofia is there any reason why you are using density as an offset rather than total? $\endgroup$ – mdewey Aug 2 '16 at 7:56
  • $\begingroup$ I thought population density would mean more dense the area, the more people would be affected. However, I am probably wrong - and population total is the right choice to model for affectedness/occurrence rates. $\endgroup$ – Sofia Aug 2 '16 at 15:21
  • $\begingroup$ @mdeweyI performed a negative binomial hurdle with and without population total as an offset (I added the results to the initial post). I am very surprised to see that the results with an offset show much higher residuals, a smaller theta, and a very slight AIC improvement. Could the log function applied to the population be the reason? I am not able to offset directly by the population total - R error $\endgroup$ – Sofia Aug 2 '16 at 15:31
  • $\begingroup$ @Sofia I am at a loss now, I think you have tried everything which I can think of. It sees that the model does not fit the data very well so you may have to report that and explain about the outliers. $\endgroup$ – mdewey Aug 2 '16 at 15:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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