I have an application where I would like to approximate a Poisson GLM with identity link, i.e. a
glm of the form
fit = glm(y~x1+x2+x3+..., family=poisson(link=identity), data=data)
using a (noniterative) weighted least squares model fit instead, ie using
fit = lm(y~x1+x2+x3+...,, weights=..., data=data)
which comes down to using a weighted nonnegative least squars model (
nnls) where both y and the covariate matrix are multiplied by
I was wondering what would be the correct/best weights to use though to take into account the mean/variance relationship?
For Poisson the variance=mean, so would that imply that I should use weights=1/variance, i.e.
weights = 1/(y+epsilon) with some small
epsilon (e.g. 1), given that
y should be a good estimator of the mean & of the variance? I would like to use weighted OLS instead of a GLM mainly for computational reasons (a
lm.fit typically being 5-10x faster than a
glm.fit) and also because the identity link Poisson GLM in R doesn't seem to be implemented in a very robust way (it often complains about it not finding valid starting values; probably this is also related to it not allowing for nonnegativity constraints on the fitted coefficients, so that predictions can go negative in some cases).
When working with a single covariate
y to stabilize the variance and doing a regular OLS fit on those also seems to give a nearly identical fit as the Poisson GLM. But this only works with a single covariate though.
I was wondering if both of these procedures are generally accepted approximations of Poisson GLMs? And how they would typically be referred to, and what a good citation would be for both recipes?
To enforce nonnegativity coefficients on the fitted coefficients I would then also like to use
glmnet does not support identity link Poisson;
h20 does but does not support nonnegativity constraints;
nnlm with Kullback-Leibler divergence & the
nnpois function from the
addreg package are other options and do have nonnegativity constraints, but both are much slower than
nnls). So would using
1/variance weights within a regular (nonnegative) weighted least squares framework be an OK thing to do to approach a Poisson error structure?
In this small example the expected slope of 1 is in fact estimated more accurately by a weighted OLS regression with 1/variance weights than by a Poisson regression, which surprised me - what would be the intuitive reason for this ? :
x=1:100 n=length(x) set.seed(1) y=rpois(n=n, lambda=x) data=data.frame(x=x,y=y) glmfit = glm(y~0+x, family=poisson(link=identity), data=data) # Poisson GLM coef(glmfit) # 1.022574 olsfit = lm(y~0+x, data=data) # regular unweighted OLS fit coef(olsfit) # 1.024758 sqrlmfit = lm(sqrt(y)~0+sqrt(x), data=data) # OLS with x & y sqrt transformed to stabilize variance coef(sqrlmfit) # 1.009303 data$weights = 1/(y+1) # approx 1/(y+eps)=1/variance for Poisson with eps here = 1 wlmfit = lm(y ~ 0+x, weights=weights, data=data) # weighted OLS with 1/variance weights coef(wlmfit) # 1.008069 - closest to expected simulated slope of 1!
Also relevant perhaps is this link.