Why does my Poisson regression fail? I'm using the R command glm(), x and y are non-negative integers, I'm using a linear link

I'm trying to regress one set of hurricane numbers per year onto another, as a way to estimate the proportion of the hurricanes that hit the US, and the uncertainty around that proportion, all in the context of using a Poisson distribution for the distribution of the number hurricanes hitting the US.

I'm using the R command:

reg1=glm(y~0+x,family=poisson(link="identity"),start=0.5)

where x is the number of hurricanes, and y is the number of hurricanes that hit the US.

This seems to me to be the obvious statistical model to use.

I've got various data-sets. It works for most of them, but fails for those data-sets in which the number of hurricanes in the y variable is quite low. For instance, here's one that fails:

x= c(1,1,3,0,1,3,2,1,2,1,1,0,5,6,1,3,5,3,4,2,3,6,7,2,2,5,2,5,4,2,0,2,2,4,6,2,3,7,4)
y= c(1,0,1,0,0,0,1,0,0,1,1,0,1,1,0,0,1,0,0,0,0,3,4,0,0,0,0,0,0,0,0,0,0,0,2,1,0,2,1)


The error I get is:

"Error: cannot find valid starting values: please specify some"

But whatever starting values I specify, it fails. From a mathematical/statistical point of view, it seems to me to be a reasonable question to ask, that should have a solution.

Any suggestions for: a) how I might get glm to work b) any alternative model I could use?

Things I've tried:

a) gaussian regression: that works, but doesn't really make sense since y is non-negative integers.

b) poisson regression with log link: that works, but that doesn't make sense either since it doesn't use a constant proportion.

c) I can estimate the proportion mean p and standard error s using

p=sum(y)/sum(x)
q=1-p
s=sqrt(pq/sum(x))


but that doesn't really take the Poisson context into account, so doesn't seem like the best solution (and it gives lower results for the standard error s than the glm does in the cases where it works).

thanks

Steve

• Wouldn't a binomial model make more sense here? If there are $X$ hurricanes and each one hits the US with some probability $p$, we have $(Y|X=x) \sim \mathrm{Bin}(x, p)$. Mar 29, 2023 at 15:55
• The "things I've tried c)" is based on the binomial. But for the cases where it works, the glm() gives larger standard errors, and better results (where by better I mean: I calculate out of sample log-likelihoods, and they are higher). Do you know of a more accurate formula for the variance of a binomial than pq/n?...there are certainly lots of formulas for confidence intervals for the binomial on wikipedia, and maybe I could use one of them to give a standard error. At the moment my conclusion is that I should use glm() where it works and the simple binomial estimators otherwise. Mar 29, 2023 at 18:25
• what do you mean by "it doesn't use a constant proportion." for your criticism of log-link poisson glm? also the binomial standard error/variance formulas you quote are for the normal approximation. There are lots of better options available, especially for the small counts as in your example Mar 29, 2023 at 19:04
• Sorry if I wasn't clear...with the log link the proportion (the slope of the link curve) varies with the x variable, and I'd rather assume it's constant. For this problem, a constant proportion is more reasonable, I think, unless we have strong evidence otherwise. Have you got a link to a better expression for binomial standard error? I haven't been able to find one yet. I can find lots of discussion about methods for confidence intervals, but not specifically about standard errors. Perhaps my use of glm() is effectively a way to get a better estimate for the standard error. thanks, Steve Mar 29, 2023 at 19:30
• why do you specifically want standard errors, as opposed to just a confidence interval? Mar 30, 2023 at 17:29

You could remove the blockade by hacking the definition of the Poisson family and allowing values of zero

fam = poisson(link="identity")
fam\$validmu = function(mu) {
# all(is.finite(mu)) && all(mu > 0)
all(is.finite(mu)) && all(mu >= 0)
}
glm(y~0+x, family=f, start = 0.3)


but this will just give a new error

Error in glm.fit(x = c(1, 1, 3, 0, 1, 3, 2, 1, 2, 1, 1, 0, 5, 6, 1, 3,  :
0s in V(mu)
Calls: glm -> eval -> eval -> glm.fit
Execution halted


The problem is that the iterative reweighted least squares procedure that is performed by the glm function is using a division by the current estimate of the variance.

The solution of the maximum likelihood equations is equivalent to an iterative weighted least-squares procedure with a weight function $$w = (d\mu/dY)^2/V$$ and a modified dependendent variable (the working probit of probit analysis) $$y = Y + (z-\mu)/(d\mu/dY)$$

For Poisson regression, $$V$$, the variance of the distribution at the current estimated mean, $$\mu$$, is equal to $$V = \mu$$. So you get a division by zero if your starting values estimate $$\mu = 0$$ which happens in the point $$x=0$$.

In your case you can exclude these points from the regression since for $$x=0$$ your model predicts $$y=0$$ anyway, no matter what the estimate is for the model coefficients.

So if asside from the above hack with the valid values for mu, we also add weights to get rid of the cases where $$x=0$$, then you can run the glm function

Call:  glm(formula = y ~ 0 + x, family = fam, weights = 1 * (x > 0),
start = 0.3)

Coefficients:
x
0.1858

Degrees of Freedom: 36 Total (i.e. Null);  35 Residual
Null Deviance:      Inf
Residual Deviance: 34.88    AIC: 68.36


This is probably because there are issues with estimating parameters that are on the boundaries, eg $$\lambda=0$$ for a Poisson($$\lambda$$). with a log link this doesn't become an issue since $$\log\lambda =x$$ is 0 only for $$x\to -\infty$$, but for identity link like you use this becomes an issue.

• It seems to me that, for the specific problem I'm looking at (regressing non-negative integers onto non-negative integers), it should be possible to write code that uses a linear link and never fails. But I guess the glm() code is a general algorithm. Perhaps I should stop typing on stack-exchange and see if I can write that code. Thanks for all the comments, btw. Hang on though...is there an analytic solution? Going look into that this morning Mar 31, 2023 at 7:42

Here's the solution to my specific problem.

(I never got glm() working, but I found an easier solution to the actual problem).

1. I don't need to use glm() at all, because there is an analytic solution to the problem I posed.

2. The analytic solution (written using R code) is:

p=sum(y)/sum(x)
s=sqrt(p/sum(x))

1. This solution can be derived analytically by maximising the log-likelihood, and looking at the curvature of the log-likelihood
2. Note how the standard error is not the usual binomial pq/n, because it's all in the context of a Poisson distribution
3. In my particular prediction context, this gives better out-of-sample predictions than the glm() solution (for those cases when the glm solution does work).