# Closest approximation of a Poisson GLM using weighted least squares analysis to take into account mean/variance relationship

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)


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 sqrt transforming x and 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 nnls or glmnet with lower.limits=0 (as 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 glmnet or 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.

• I am curious about why you want to do this. Seems like you should just be using poisson glm... – probabilityislogic Apr 25 '18 at 9:25
• Well mainly to be able to be able to use some penalized regression frameworks that don't support families other than gaussian but do support weights, and still take into account that my noise is poisson – Tom Wenseleers Apr 25 '18 at 9:27
• And also for performance reasons - GLMs are fitted iteratively, whereas linear models are not... – Tom Wenseleers Apr 25 '18 at 9:29
• And I also wanted a Poisson model with an identity link rather than a log link, which is not incorporated in glmnet – Tom Wenseleers Apr 25 '18 at 9:34
• You should look into the CRAN package h2o which support regularized poisson regression with identity link function. – kjetil b halvorsen Apr 25 '18 at 19:41

The square root transformation is a totally acceptable way to proceed.

I'd drop the weights = 1/(y+epsilon) idea. Good thinking but this is a hack that could depend pretty heavily on this choice of arbitrary epsilon.

Check out the gls function from nlme. It will likely have what you need.

> glmfit2 <- gls(y ~ 0 + x, weights = varPower(), data = data)
> summary(glmfit2)
Generalized least squares fit by REML
Model: y ~ 0 + x
Data: data
AIC      BIC   logLik
632.624 640.4093 -313.312

Variance function:
Structure: Power of variance covariate
Formula: ~fitted(.)
Parameter estimates:
power
0.5014736

Coefficients:
Value  Std.Error  t-value p-value
x 1.022571 0.01231209 83.05423       0

Standardized residuals:
Min         Q1        Med         Q3        Max
-2.6004374 -0.7538303  0.1200103  0.6402228  2.2430706

Residual standard error: 0.8599757
Degrees of freedom: 100 total; 99 residual


Look at that power estimate!

Need to ensure your coefficient is greater than 0? Then reparameterize $$\beta = \exp(\theta)$$ and use the nonlinear version of gls, gnls:

> glmfit3 <- gnls(y ~ exp(theta) * x, start = list(theta = 2), weights = varPower(),
data = data)
> summary(glmfit3)
Generalized nonlinear least squares fit
Model: y ~ exp(theta) * x
Data: data
AIC      BIC    logLik
625.6613 633.4768 -309.8306

Variance function:
Structure: Power of variance covariate
Formula: ~fitted(.)
Parameter estimates:
power
0.49894

Coefficients:
Value  Std.Error  t-value p-value
theta 0.02232537 0.01202577 1.856461  0.0664

Standardized residuals:
Min         Q1        Med         Q3        Max
-2.5979856 -0.7547549  0.1201196  0.6395667  2.2459128

Residual standard error: 0.8679793
Degrees of freedom: 100 total; 99 residual
> exp(.022325)
[1] 1.022576


You mention OLS is "more accurate" than GLS for your example. Let the seed go and run the simulation 1000 times. See which is closer to the truth on average.

• Thanks for this - but is gls with this correlation structure not just the same as iteratively reweighted least squares, where the weights are initialised with 1, the model is estimated using OLS, and the fitted model is used to reestimate the variance weights using whatever functional form is requested (power relationship here), which is then further iterated? (en.wikipedia.org/wiki/Iteratively_reweighted_least_squares) Thanks also for the pointer to gnls - that's using the port algo though, which in my experience is not that robust, especially not for larger models... – Tom Wenseleers Apr 12 '19 at 13:54
• The sqrt transform option I am afraid will also not work for me for models with >1 covariate, as the sqrt transform breaks the additivity of my predictors... You were right about my weighted LS fit being slightly more biased though than GLM... I had another application where my stability of the estimated using weighted least squares was better thought than with GLM, but I'll probably post that as a separate question... – Tom Wenseleers Apr 12 '19 at 13:57
• Tom are you sure gls and gnls are not using iteratively reweighted least squares? There must be some kind of loop going on because these functions slow down pretty quickly as the sample size gets larger. For the square root transformation, my response is that additivity was just an approximation in the first place, unless there was something special about the original scale so that you expected additivity to hold exactly. – Ben Ogorek Apr 12 '19 at 14:34
• No worries - and yes, this approximation indeed only works well as a correction if you have large counts, which is the case in my applications. In that case, it's just a simple fix for heteroscedasticity. But for any counts > say 100 it was really quite accurate, and the computational cost is ca 10 x lower than doing a GLM or a nonnegative Poisson fit. I was using all this in the context of nonnegative matrix factorization where I have to do this millions of times and then a speed difference of 10x is quite large... For a single fit it doesn't matter of course. – Tom Wenseleers Apr 26 '19 at 20:20
• Also added one more example with those 1/(1+y) weights here: stats.stackexchange.com/questions/203355/… – Tom Wenseleers Apr 28 '19 at 0:22

Just figured I was being stupid - with identity link the iteratively reweighted LS algo to fit GLMs just ends up regressing yhat on X in each iteration (cf irls algo http://bwlewis.github.io/GLM/, with identity link z=y and W=1/expected variance) and as initialization for yhat R uses in the glm.fit source code y+eps with eps=0.1. So what I was doing with my weighed least squares was in fact just a single iteration of the IRLS algo to fit GLMs [which boils down to a Fisher scoring variant of the Newton-Raphson algorithm], which explains why it was working so well for me. It also answers my problem - what I was doing then was in fact the closest noniterative (single step) approximation of a Poisson identity link GLM using weighted least squares... For nonnegativity constrained regression the same recipe applies and all I have to do is change lm.wfit by nnls (and iterate this if I would like to approach the true maximum likelihood objective more accurately). I posted a worked out example here.

Square root transforming I cannot do btw as that would break the additivity of my predictors, which I know a priori should hold in my application.

• So are you going to write the loop yourself? I'd be interested in seeing your code if you do. – Ben Ogorek May 15 '19 at 16:09
• @ben-ogorek Yes I coded it - I was planning to write it up as a CRAN package to fit nonnegative identity link GLMs, with additional functions also for ridge and adaptive ridge penalized nonnegative GLMs (which can be done by row augmenting X with diag(sqrt(lambda)) plus (building on this) L0-norm penalized nonnegative GLMs (aka best subset), as those can be fit using an iterative adaptive ridge procedure. So far I only did a crude version of the nnglm function with or without a ridge penalty. – Tom Wenseleers May 16 '19 at 10:33
• @ben-ogorek I was also still contemplating how to best do post-selection inference. I did a first attempt here stats.stackexchange.com/questions/373253/… but Jonathan Tayler tells me I should do it like this gist.github.com/jonathan-taylor/… as in the selectiveInference package... – Tom Wenseleers May 16 '19 at 10:34
• @ben-ogorek But if you send me an email (tom.wenseleers@kuleuven.be) I can send you what I already have... – Tom Wenseleers May 16 '19 at 11:21