# Computational error running regression model

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

## 1 Answer

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).

• 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? – Guido167 Jul 30 '16 at 9: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. – Guido167 Jul 30 '16 at 9:17
• 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 9:17
• How do I make month a quadratic term, simply by adding month^2? – Guido167 Jul 30 '16 at 9:19
• 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 2 '16 at 20:36