Computational error running regression model

Question

This problem has held me up for three days now, so I really hope somebody here has a solution for the problem.

I have a model with an excessive number of zeros, so I use a zero-inflated poisson regression model with the following code and summary.

cr_f1 = formula(cr ~ depth + habtype2 + month + year + lightregime + depth*month + depth*lightregime + depth*habtype2 + habtype2*year | depth + habtype2 + month + year + lightregime + depth*month + depth*lightregime + depth*habtype2)
summary(zeroinfl(cr_f1, dist = "poisson", link = "logit", data = allUVCdata))

Call:
zeroinfl(formula = cr_f1, data = allUVCdata, dist = "poisson", link = "logit")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.6430 -0.5680 -0.2893  0.1426 16.8090 

Count model coefficients (poisson with log link):
                            Estimate Std. Error z value Pr(>|z|)   
(Intercept)                -4.515522   2.182503  -2.069  0.03855 * 
depth                       0.108941   0.072278   1.507  0.13175   
habtype2Pinnacles           0.879765   0.791166   1.112  0.26614   
habtype2Unexposed          -0.604246   0.786129  -0.769  0.44211   
month2                      0.628468   0.380450   1.652  0.09855 . 
month3                      0.309282   0.367690   0.841  0.40026   
month4                      0.649411   0.371667   1.747  0.08059 . 
month5                      0.758717   0.364079   2.084  0.03717 * 
month6                      0.467611   0.341024   1.371  0.17031   
month7                      0.523043   0.343363   1.523  0.12768   
month8                      0.563272   0.356843   1.578  0.11445   
month9                      0.204509   0.400398   0.511  0.60952   
month10                     0.662415   0.341616   1.939  0.05249 . 
month11                     0.934844   0.335077   2.790  0.00527 **
month12                     0.252216   0.360512   0.700  0.48417   
year2013                   -1.271010   1.282158  -0.991  0.32154   
year2014                    1.221887   0.753644   1.621  0.10495   
year2015                   -0.463176   0.771131  -0.601  0.54808   
lightregimeLight            2.754925   1.948779   1.414  0.15746   
depth:month2               -0.019864   0.008906  -2.230  0.02572 * 
depth:month3               -0.014157   0.008106  -1.747  0.08071 . 
depth:month4               -0.020553   0.008332  -2.467  0.01364 * 
depth:month5               -0.021213   0.008373  -2.533  0.01129 * 
depth:month6               -0.013561   0.007393  -1.834  0.06663 . 
depth:month7               -0.015043   0.007544  -1.994  0.04615 * 
depth:month8               -0.017383   0.008011  -2.170  0.03003 * 
depth:month9               -0.012340   0.008990  -1.373  0.16988   
depth:month10              -0.019631   0.007629  -2.573  0.01008 * 
depth:month11              -0.024101   0.007611  -3.167  0.00154 **
depth:month12              -0.014319   0.007952  -1.801  0.07174 . 
depth:lightregimeLight     -0.079860   0.071024  -1.124  0.26084   
depth:habtype2Pinnacles    -0.006819   0.011178  -0.610  0.54182   
depth:habtype2Unexposed     0.014857   0.011103   1.338  0.18086   
habtype2Pinnacles:year2013  1.351509   1.277930   1.058  0.29025   
habtype2Unexposed:year2013  1.538282   1.256047   1.225  0.22069   
habtype2Pinnacles:year2014 -1.213233   0.754305  -1.608  0.10775   
habtype2Unexposed:year2014 -0.495275   0.726863  -0.681  0.49563   
habtype2Pinnacles:year2015  0.389117   0.775476   0.502  0.61582   
habtype2Unexposed:year2015  0.659117   0.750396   0.878  0.37975   

Zero-inflation model coefficients (binomial with logit link):
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -4.61555    7.04621  -0.655 0.512442    
depth                    0.28728    0.28211   1.018 0.308524    
habtype2Pinnacles        9.41037    3.82210   2.462 0.013813 *  
habtype2Unexposed        2.11213    1.46465   1.442 0.149282    
month2                   8.67847    3.91193   2.218 0.026523 *  
month3                   7.12210    3.86428   1.843 0.065320 .  
month4                   4.10296    2.41285   1.700 0.089044 .  
month5                  12.76919    4.28035   2.983 0.002852 ** 
month6                   3.57695    2.49820   1.432 0.152198    
month7                   5.85534    3.27394   1.788 0.073700 .  
month8                   5.59503    3.33054   1.680 0.092974 .  
month9                   4.22953    3.76919   1.122 0.261807    
month10                  6.35022    3.59424   1.767 0.077265 .  
month11                  5.92079    3.36405   1.760 0.078404 .  
month12                  4.36214    3.17233   1.375 0.169113    
year2013                -0.18722    0.42651  -0.439 0.660688    
year2014                -1.50194    0.45263  -3.318 0.000906 ***
year2015                -9.79773    4.87536  -2.010 0.044469 *  
lightregimeLight         0.79826    5.62419   0.142 0.887133    
depth:month2            -0.39212    0.16795  -2.335 0.019557 *  
depth:month3            -0.36363    0.16695  -2.178 0.029397 *  
depth:month4            -0.21521    0.10211  -2.108 0.035059 *  
depth:month5            -0.57543    0.16933  -3.398 0.000678 ***
depth:month6            -0.24336    0.10398  -2.341 0.019256 *  
depth:month7            -0.33704    0.13975  -2.412 0.015874 *  
depth:month8            -0.35343    0.14683  -2.407 0.016082 *  
depth:month9            -0.31787    0.16903  -1.881 0.060026 .  
depth:month10           -0.37550    0.16021  -2.344 0.019087 *  
depth:month11           -0.34650    0.14821  -2.338 0.019397 *  
depth:month12           -0.29639    0.14221  -2.084 0.037142 *  
depth:lightregimeLight   0.08117    0.21795   0.372 0.709571    
depth:habtype2Pinnacles -0.57765    0.17049  -3.388 0.000704 ***
depth:habtype2Unexposed -0.17897    0.06252  -2.863 0.004200 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 146 
Log-likelihood: -3977 on 72 Df

So I included the interaction 'habtype2*year' in the count part of the formula, but now want to include it in the second model aswel (the binomial), but if I do I get the following error:

cr_f1 = formula(cr ~ depth + habtype2 + month + year + lightregime + depth*month + depth*lightregime + depth*habtype2 + habtype2*year | depth + habtype2 + month + year + lightregime + depth*month + depth*lightregime + depth*habtype2 + habtype2*year)

summary(zeroinfl(cr_f1, dist = "poisson", link = "logit", data = allUVCdata))

Error in solve.default(as.matrix(fit$hessian)) : 
  system is computationally singular: reciprocal condition number = 2.08629e-37

This also happens if I want to try to include any of the other interaction terms that I still want to put into the model ("monthyear", "monthlightregime" and "month*habtype2").

I searched here on the forum and on google, seems like more people have encountered this error (also in other functions that doing a zeroinfl), but I have not found any suitable solution.

Data: sightings of as species on 29 different locations, >5200 observations (including zeros).

What could possibly solve this, so that I can run the model with the interaction terms that I want?

EDIT: added some new output to give insight to the problem.

allUVCdata$year = as.numeric(as.character(allUVCdata$year))
cr_f1 = formula(cr ~ depth + lightregime + month + year + habtype2 + month*habtype2 + month*year + habtype2*year + depth*month)
summary(hurdle(cr_f1, dist = "poisson", link = "logit", data = allUVCdata))

Error in solve.default(as.matrix(fit_count$hessian)) : 
  system is computationally singular: reciprocal condition number = 6.23277e-26
> allUVCdata$year = as.factor(as.character(allUVCdata$year))
> table(allUVCdata$year, allUVCdata$cr)

          0    1    2    3    4    5    6    7
  2012  750  149   25   12    3    0    0    0
  2013 1133  209   69   16    4    1    1    0
  2014  844  387  142   42   11    7    0    1
  2015  833  401  125   31    5    3    1    2
> table(allUVCdata$month, allUVCdata$cr)

       0   1   2   3   4   5   6   7
  1  299  53  18   7   1   1   1   2
  10 346 104  40   9   4   4   0   0
  11 328 114  43  17   5   0   0   0
  12 350 112  29  10   2   0   0   0
  2  248  80  16   1   1   0   0   1
  3  303  82  24   6   0   0   1   0
  4  329  93  32   4   0   0   0   0
  5  277 105  28   9   1   1   0   0
  6  312 111  36  12   3   2   0   0
  7  362 113  46  14   5   2   0   0
  8  213 100  25   8   1   1   0   0
  9  193  79  24   4   0   0   0   0
> table(allUVCdata$month, allUVCdata$year)

     2012 2013 2014 2015
  1    41  129   72  140
  10   91  149  152  115
  11  112  121  150  124
  12  112  124  154  113
  2    35  108   79  125
  3    33  149  101  133
  4    88  105  150  115
  5   101  108   95  117
  6    94  115  142  125
  7   118  133  153  138
  8    61  114  100   73
  9    53   78   86   83

table(allUVCdata$habtype2, allUVCdata$year)

            2012 2013 2014 2015
  Exposed     93  138  120  144
  Pinnacles  274  386  339  338
  Unexposed  572  909  975  919

What is happening in month 5 and year 2015? Do you have any non-zero counts there, or all zero? The values of the coefficients for those two are very large in magnitude in the zero part of your model. — mdewey, Jul 28 '16 at 12:52

score 3 · Answer 1 · edited Apr 13 '17 at 12:44

3

The problem is likely due to the model matrix for either the count or zero inflated part (though it looks like they are the same from your code) being close to singular.

One way forward is to inspect the model matrices and see if you can see if there is linear dependence between any of the rows or columns, or something else causing singularity.

That could prove tricky to do, so before you try that, try simplifying the model. Does month really need to be a factor ? It doesn't seem so to me. In the count part of the model, the estimates for the main effects and the interactions with depth are all of a similar magnitude so month could be numeric. With the zero inflated part, the estimates main effects of month indicate some nonlinearity, but the those for the interaction are similar, so I would again use month as numeric, and include a quadratic term (perhaps after centering to avoid collinearity between the linear and quadratic terms).

edited Apr 13 '17 at 12:44

Community

1

answered Jul 28 '16 at 10:54

Robert Long

53,316
10
84
148

Thank you @Robert Long. Month has to be part of the model as it is the only time unit that I have (except year). So if I understand you correctly, I have to make 'month' numeric? But it is an ordinal variable, so then it would differentiate between the value of the month 1/January would considered to be lower than 2/February.. That doesn't make much sense? – Shark167 Jul 30 '16 at 09:10
I tried making month numeric now, but still get the same error. The reciprocal condition number is now closer to zero `reciprocal condition number = 1.54988e-17`. – Shark167 Jul 30 '16 at 09:17
1

If you change month to numeric then you will get one estimate for it (the linear trend) which is normally what you want. Keeping it as a factor allows nonlinear change but you can also do that more cleanly by introducing a quadratic (and higher order if necessary) terms. Using month as a factor makes no sense to me and complicated the model and it's interpretation unnecessarily. – Robert Long Jul 30 '16 at 09:17
How do I make month a quadratic term, simply by adding `month^2`? – Shark167 Jul 30 '16 at 09:19
You can put `I(month^2)` but I would try to keep the model simple until the numerical error is solved. The problem seems to be connected to year. What happens if you make year numeric ? Also, what is the output from `table(year, cr)` – Robert Long Jul 30 '16 at 09:22
I added the output of both in the opening post under edit. Changing year to numeric still gives me the error, I also added the results of three tables. – Shark167 Jul 30 '16 at 09:49
I doubt it will fix the problem but can you change `habtype2*year` to `habtype2:year` and also post the result of `table(habtype2,year)` – Robert Long Jul 30 '16 at 10:22
changing the `*` to `:` didn't make a difference. I added the contingency table `table(habtype2, year) now to the results of the last edit I made. Can't seem to find a reason why the model doesn't run.. it had been bugging me for three days now.. – Shark167 Jul 30 '16 at 10:56
strange this is that replacing month with quarter does not solve the problem. But replacing it by hurseason (hurricane season y/n) does work. But the problem with this that I lose the overall temporal data.. and now I only have the hurseasons as temporal predictor (and year). – Shark167 Jul 30 '16 at 13:46
Let us [continue this discussion in chat](http://chat.stackexchange.com/rooms/43245/discussion-between-robert-long-and-guido167). – Robert Long Jul 30 '16 at 13:59
1

Just a quick comment about the solution to the problem: @Guido167 sent me his data and it turned out that indeed there was a problem with the conditioning of the regressor matrix as suspected by @RobertLong. For certain combinations of `month:lightregime` and `month:habtype2` there was quasi-separation. Thus, there were interaction terms in which no non-zero responses were observed. This leads to infinite coefficients in the zero-inflation part and lack of variation in the response for the count part. Omitting these (non-significant) interactions seems the best solution for this data set. – Achim Zeileis Aug 02 '16 at 20:36

Computational error running regression model

1 Answers1