18
$\begingroup$

I'm trying to fit a regression to explain the number of homicides in each district of a city. Although I know that my data follows a Poisson distribution, I tried to fit an OLS like this:

$log(y+1) = \alpha + \beta X + \epsilon $

Then, I also tried (of course!) a Poisson regression. The problem is that I have better results in the OLS regression: the pseudo-$R^2$ is higher (0.71 vs 0.57) and the RMSE as well (3.8 vs 8.88. Standardized to have the same unit).

Why? Is it normal? What's wrong on using the OLS no matter what the distribution of the data is?

edit Following the suggestions of kjetil b halvorsen and others, I fitted the data through two models: OLS and Negative Binomial GLM (NB). I started with all the features I have, then I recursively removed one by one the features which were not significant. OLS is

$\sqrt{\frac{crime}{area}} = \alpha + \beta X + \epsilon $

with weights = $area$.

summary(w <- lm(sqrt(num/area) ~  RNR_nres_non_daily + RNR_nres_daily + hType_mix_std + area_filtr + num_community_places+ num_intersect + pop_rat_num + employed + emp_rat_pop + nden_daily + nden_non_daily+ bld_rat_area + bor_rat_area + mdist_highways+ mdist_parks, data=p, weights=area))

error2 <- p$num - (predict(w, newdata=p[,-1:-2], type="response")**2)*p$area

rmse(error2)
[1] 80.64783

The NB predicts the number of crime, with the district's area as offset.

summary(m3 <- glm.nb(num ~  LUM5_single  + RNR_nres + mdist_daily + mdist_non_daily+ hType_mix_std + ratio_daily_nondaily_area + area_filtr + num_community_places  + employed  + nden_daily + nden_non_daily+ bld_rat_area + bor_rat_area + mdist_smallparks + mdist_highways+ mdist_parks + offset(log(area)), data=p, maxit = 1000))

error <- p$num - predict(m3, newdata=p[,-1:-2], type="response")

rmse(error)
[1] 121.8714

OLS residuals:

enter image description here

NB residuals

enter image description here

So the RMSE is lower in the OLS but it seems that the residuals are not so Normal....

Improve this question
$\endgroup$
13
  • $\begingroup$ Can you post some more details? What are the nature of the data? that is, what is the response variable counting? what is the explanatory variables? $\endgroup$
    – kjetil b halvorsen ♦
    Jun 1, 2016 at 8:24
  • $\begingroup$ @kjetilbhalvorsen the dependent variable is the number of homicides per district (112 districts). The indipendent ones are the structural characteristics of the city (street intersections, POIs etc) $\endgroup$
    – marcodena
    Jun 1, 2016 at 9:15
  • 2
    $\begingroup$ If I were fitting this model using a Poisson regression I would include log(districtsize) as an offset to account for districts not all being the ame size. Unless they are. $\endgroup$
    – mdewey
    Jun 1, 2016 at 13:46
  • 1
    $\begingroup$ What is your rationale behind thinking that comparing the OLS $R^2$ with the $pseudo-R^2$ from a ML estimation (and the $RMSE$), gives you an indication for how good a certain model is? OLS, by construction, maximizes $R^2$. Is the Poison regression constructed so as to maximize the $pseudo-R^2$? I don't think so, and I don't think that this comparison is useful. $\endgroup$
    – coffeinjunky
    Jun 1, 2016 at 14:51
  • 1
    $\begingroup$ Another thing to add - the $R^2$ from ols is telling the % of variance explained in $z=\log (y+1) $ whereas the poisson psuedo $R^2$ is trying to give an indication of % of the variance of $y $ that is explained. This could also explain the difference $\endgroup$
    – probabilityislogic
    Jun 1, 2016 at 15:03

4 Answers 4

Reset to default
16
$\begingroup$

I suspect that part of the problem may lie in your choice of performance metric. If you measure test performance using RMSE then training the model to minimise the MSE matches the test criterion, giving a hint as to what is considered important. You may find that if you measure the test performance using the negative log-likelihood of the test set using a Poisson likelihood that the Poisson model works better (as might be expected). This may be a minor issue compared with the other issues raised, but it might be a useful sanity check.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ +1. If the OPs objective was prediction, there might actually be a rationale for using an OLS model instead! Nonetheless, the classical error-based inference arising from OLS cannot/should not be applied in GLMs. One could inspect studentized residuals, or a better option would be comparing models with AIC. $\endgroup$
    – AdamO
    Jul 19, 2016 at 12:31
11
$\begingroup$

First, with such data I would expect overdispersion (if you don't know what that is, see https://stats.stackexchange.com/search?q=what+is+overdispersion%3F ).

That would have to be addressed with a Poisson glm, but is not an issue with usual linear regression. As said in a comment, with a poisson glm you want to include $\log(\text{DistrictSize})$ as an offset, with a linear regression you will need to use as response variable $\frac{\text{Nr. homicides}}{\text{District Size}}$. One possible reason for the discrepancy of results is that you have treated this problem differently in the two cases. You could post here some plots of results, like residual plots, so we can see what is happening. Or you could post your data as a table in the original post .... could be interesting to have a look.

Another issue is the transformation you used with the linear regression. The usual variance stabilizing transformation used with count data is the square root, not the logarithm.

Another issue is the choice of transformation used with linear regression. When using as response $Y_i/x_i$, you wil need weighted linear regression. Assuming as an approximation that $Y_i \sim \text{Poisson}(\lambda x_i)$, we have $$ \DeclareMathOperator{\E}{\mathbb{E}} \DeclareMathOperator{\V}{\mathbb{V}} \E \frac{Y_i}{x_i} \propto \lambda \\ \V \frac{Y_i}{x_i} \propto x_i^{-1} $$ So you should use weighted linear regression with $x_i$ as weight. A simple analysis shows that, as an approximation, the same weights are appropriate with $\sqrt{Y_i/x_i}$ or $\log (Y_i/x_i +1)$ as responses. 

    EDIT

As for your additional analysis in the post, note that rmse's cannot be compared directly between the two models, as different responses are used! To make a direct comparison, you will need to backtransform the predicted values to the original scale. Then you can calculate rmse's yourself, and see. But note that predictions obtained after backtransformation can be biased, because of nonlinearities. So some adjustment to the backtransformed predictions could make them more useful. In some cases, such could be calculated theoretically, ot you could just use a bootstrap.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I fitted the models as you suggested, although I did not really understand the reson behind the weighted OLS. What do you think? $\endgroup$
    – marcodena
    Jun 2, 2016 at 14:09
6
$\begingroup$

There are lots of choices of pseudo $R^2$'s. Lots of them are very flawed. Generally speaking, there's usually no reason that the $R^2$ produced from OLS will be a comparable value to a given pseudo $R^2$; rather, pseudo $R^2$'s are typically used for comparing models of the same distributional family.

Improve this answer
$\endgroup$
2
$\begingroup$

It is true that your data is not Normally distributed (which I presume is why you also ran a Poisson regression) but your data is likely not a Poisson distribution either. The Poisson distribution assumes that the mean and the variance are the same, which is likely not the case (as mentioned in other answers - you can capture this discrepancy and incorporate it into the model). Since your data isn't really a perfect fit for either model, it makes sense that OLS may perform better.

Another thing to note is that the ordinary least squares estimates are robust to non-Normality which may be why you're getting a reasonable model. The Gauss-Markov Theorem tells us that OLS coefficients estimates are the best (in terms of mean squared error) linear unbiased estimators (BLUE) under the following assumptions,

  • The errors have a mean of zero
  • The observations are uncorrelated
  • The errors have have constant variance

There is no assumption of Normality here so your data can very well be reasonable for this model! With that being said, I would look into a Poisson model with an over-dispersion parameter baked in there and you should get better results.

Improve this answer
$\endgroup$
3
  • $\begingroup$ @TynnaDoStat thanks! I fitted two models now, one with dispersion parameter. What do you think? $\endgroup$
    – marcodena
    Jun 2, 2016 at 14:09
  • 2
    $\begingroup$ Variance = mean for a Poisson distribution is often invoked as a problematic assumption for Poisson regression, but the point is not as difficult as is implied here. Despite its name the main idea of Poisson regression is that of a log link function; assumptions about conditional distribution are not nearly so important. What's likely if the assumptions don't all hold is mainly that the standard errors are off unless you adjust, but the fit will often make sense. $\endgroup$
    – Nick Cox
    Jun 2, 2016 at 16:04
  • 2
    $\begingroup$ Indeed Poisson regression can make sense for non-negative measured responses where variance and mean don't even have the same dimensions. See e.g. blog.stata.com/2011/08/22/… $\endgroup$
    – Nick Cox
    Jun 2, 2016 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.