2

Suppose I have some problem typical linear inverse problem, where I have a relation $Ax = \vec{b}$, and want to solve for $x$ in terms of $A$ and $\vec{b}$. Furthermore, assume that $\vec{b}$ is a composed of independent heteroskedastic Poisson distributed terms. Given this information, OLS shouldn't be applicable. Yet when I run simulations, I find that OLS and Poisson regression produce closely comparable results.

The experiments in detail:

  1. Start with some known (say 30x2) matrix $A$ and some true (2x1) vector $\vec{x}$, and multiply them out to get a no noise true value for $\vec{b}$.
  2. Then, replace each $b_i$ with a random number drawn from a Poisson distribution with mean $b_i$ (simulating Poisson noise).
  3. Finally, using the noisy $\vec{b}$ value and the matrix $A$, attempt to recover as closely as possible the vector $\vec{x}$.
  4. Repeat steps 2 and 3 many times, and evaluate the mean of the error and the standard deviation of the error.

When I do this experiment, the average standard deviation for the ols model errors is in some cases lower than the standard deviation for the Poisson. Is there an analytical reason for this? Can I predict when this will happen based on the matrix $A$?

danton
  • 23
  • 1
    What if $A x$ has negative entries? How do you generate $b_i$? How do you do Poisson regression? The standard one assumes that $\log E(y|x)$ is linear in $x$ not $E(y|x)$. – passerby51 Oct 09 '22 at 17:21
  • This doesn't describe a statistical model. Maybe you mean $E[\vec b] = Ax$ and you have independent observations of $\vec b$? If most of these expectations are largish (above $10$ ought to work fine) or if you have loads of data, OLS ought to yield a comparable estimate. You can't really compare the OLS errors to the Poisson ones unless perhaps you are using an appropriate weighted version of OLS, because OLS will not produce good estimates of heteroscedastic errors. – whuber Oct 09 '22 at 17:26
  • @passerby51 Neither A nor x are allowed to have negative entries. And the $b_i$ are the entries in the vector $\vec{b}$, generated by multiplying Ax. Then those are replaced with noisy equivalents. – danton Oct 09 '22 at 17:45
  • @whuber I am essentially generating independent observations of $\vec{b}$ by analytically multiplying Ax, then noising the $b$, and then trying to recover $x$ from the noised b and then A matrix. Can I ask why the OLS should yield comparable estimates if the expectations are above 10 in your model? And the interesting thing which I have observed is that the errors are similar between the models, despite the heteroscdastic nature of them. That phenomena is what I want to understand. – danton Oct 09 '22 at 17:48
  • @SextusEmpiricus Using standard libraries, linear link Poisson using R's GLM package. – danton Oct 09 '22 at 20:19
  • @SextusEmpiricus the average standard deviation for the ols model errors is in some cases lower than the average standard deviation for the Poisson – danton Oct 09 '22 at 21:01

1 Answers1

2

I can reproduce this situation in a simple way by creating a situation where all $E[b_i]$ are equal or unequal by changing the slope parameter in a linear model. This will influence the homo/heteroscedasticity and relative accuracy of OLS and GLM.

I use $x_i = 1, 2, 3, \dots, 30$ and $$b_i|x_i \sim Poisson(intercept + slope \times x_i)$$

Then I compute different simulations of fitting and compute the variances of the intercept coefficient when the true slope coefficient has different values.

As I expected, when the slope is around zero, which is when the the values of $b_i$ are very similar, then the OLS estimated intercept has a smaller variance than with GLM.

The reasoning behind it can be seen intuitively as following

  1. For simple cases, the ordinary least squares regression (OLS) is the best linear unbiased estimator. More precisely, when the conditional distributions of the observations are all independent and have unequal variances (which is when the slope is zero).

  2. For more complicated cases, generalized least squares regression (GLS) is the best linear unbiased estimator. More precisely, when the variances are not equal (and when there are correlations).

  3. However GLS only performs best when the (unequal) variances are known. When some method is used that estimate the variances then the performance will be reduced.

  4. In some way we could see the general linear model (GLM), with the iterated least squares regression, as some method to do GLS with estimating the variances based on the estimates of the means. And it will perform less well than GLS (if we know the variance instead of estimate the variance).

  5. So for cases when the variance of the data points are nearly equal, then OLS will perform better than GLM/Poisson regression.

Whether one can make an exact computation or express the conditions, I am not sure. But typically the Poisson regression might perform worse when the outcome variable $b_i$ has relatively little variation such that the sample is more or less homoscedastic.

comparison for different heteroscedasticity

set.seed(1)
a = 5
n = 30
x = 1:n

sim = function(b = 0){


generate observations


y = rpois(n, a + b * x)


OLS and GLM models


mod1 = lm(y ~ x)
  mod2 = glm(y ~ x, family = poisson(link='identity'), start = c(a,b))


return intercept coefficients


return(c(coef(mod1)[1], coef(mod2)[1]))
}


bs = c()
var1 = c()
var2 = c()


compute for different slopes 'btest'


for (i in 1:10) {
  btest = (i-1)/50
  sims = replicate(10^4, sim(btest))


store values in vector bs, var1 and var2


bs = c(bs, btest)
  var1 = c(var1, var(sims[1,]))
  var2 = c(var2, var(sims[2,]))
}


plot results


plot(bs, var1, xlab = "value of true slope parameter", ylab = "variance of estimated coefficients from simulations")
points(bs,var2, col = 2 )


legend(0, 0.9, c("OLS", "GLM/Poisson"), col = c(1,2), pch = 1)
Sextus Empiricus
  • 77,915