10
$\begingroup$

I have been puzzling over a toy regression problem with simulated Poisson distributed data and hoping someone more educated in statistics can help me gain some insight about the following observation.

Libraries used

library(tidyverse)
library(cowplot)
library(broom)
library(modelbased)
library(parameters)
library(ggbeeswarm)

Data generation

I simulated count values using rpois for two scenarios:

  1. Counts of traffic accidents where lambda is linearly proportional to traffic volume.
  2. Counts of traffic accidents when lambda scales as the exponent of traffic volume.
# number of observations
n_obs = 10

# generate log dependent data
traffic_volume = log(c(1, 2, 4, 7, 10, 15))
log_data = tibble(
  volume=traffic_volume,
  lambda=exp(0.43*volume + 0.2)
) %>%
  rowwise %>%
  mutate(accident_counts = list(rpois(n_obs, lambda = lambda))) %>%
  mutate(observed_avg_accidents = mean(accident_counts)) 

# generate linear dependent data
linear_data = tibble(
  volume=traffic_volume,
  lambda = 0.43*volume + 0.2
) %>%
  rowwise %>%
  mutate(accident_counts = list(rpois(n_obs, lambda = lambda))) %>%
  mutate(observed_avg_accidents = mean(accident_counts)) 

Modelling

I fit two glms to each data set. One using the gaussian family and one using the poisson family. I used the "log" linker for log data and the "identity" linker for the linear data.

# fit
proc_list = list(
  log=list(data=log_data, linker="log"),
  linear=list(data=linear_data, linker="identity")
)
models = map(proc_list, function(proc) {
  poisson_model <- glm(
    accident_counts ~ volume,
    data = proc$data %>% unnest(accident_counts),
    family = poisson(link=proc$linker),
  )
  gaussian_model = glm(
    accident_counts ~ volume,
    data = proc$data %>% unnest(accident_counts),
    family = gaussian(link=proc$linker),
    start=c(1, 1)
  )
  return(list("poisson"=poisson_model, "gaussian"=gaussian_model))
})

Results

> compare_models(unlist(models, recursive=FALSE))

Parameter    |        log.poisson |        log.gaussian |    linear.poisson |    linear.gaussian
------------------------------------------------------------------------------------------------
(Intercept)  | 0.01 (-0.40, 0.43) | -0.03 (-0.59, 0.53) | 0.39 (0.06, 0.73) | 0.43 (-0.08, 0.94)
volume       | 0.52 ( 0.32, 0.72) |  0.54 ( 0.30, 0.79) | 0.36 (0.13, 0.59) | 0.34 ( 0.05, 0.62)
------------------------------------------------------------------------------------------------
Observations |                 60 |                  60 |                60 |                 60

Visualization

# create the visualization grid and predict values
viz_grid = modelbased::visualisation_matrix(tibble(volume=traffic_volume)) %>% as_tibble
augmented = map_df(unlist(models, recursive=FALSE), function(.x) {
  augment(.x, newdata=viz_grid, type.predict="response")
}, .id="model")

# separate model and data labels
augmented = augmented %>%
  separate("model", c("data", "regression"), sep="\\.")

p = map_df(proc_list, ~.x$data, .id="data") %>%
  unnest(accident_counts) %>%
  ggplot(aes(volume, accident_counts)) +
  # geom_violin(adjust=1.5) +
  geom_quasirandom() +
  geom_point(
    data=~.x %>% distinct(data, volume, observed_avg_accidents),
    aes(volume, observed_avg_accidents),
    color="red"
  ) +
  geom_line(data=augmented, aes(volume, .fitted, color=regression)) +
  facet_wrap(~data, labeller=label_both) +
  theme_gray(base_size=16)
p %>% ggsave(file="temp.pdf", w=8, h=4)

enter image description here

Question

Why is there essentially no difference in the fits regardless of family function? I purposely chose a small number of observations and relatively small lambda values hoping the fit using the gaussian family would break down. But this did not happen. If the data generation process here is truly Poisson, wouldn't the choice of family function affect the fit?

New contributor
Mikhail is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
5
  • $\begingroup$ Welcome to Cross Validated! Shouldn’t the $y$-coordinates be the same in each side of your plot? $\endgroup$
    – Dave
    May 15 at 23:41
  • $\begingroup$ @Dave No, OP is fitting 4 models: (linear, log) x (gaussian, poisson). The question is why the blue line (Gaussian fit) is so similar to the red line (Poisson fit) in each example, which I think you touched on in your answer $\endgroup$
    – ischmidt20
    May 16 at 18:32
  • 3
    $\begingroup$ Although there is little difference in the $\hat\beta$ estimates here, their confidence intervals differ: they are a bit wider for the Gaussian GLM than for the Poisson GLM. By using a more-appropriate model for your data, you can get more-appropriate inferences (p-values and CI widths), even if the point estimates don't change much. $\endgroup$
    – civilstat
    May 16 at 19:37
  • 3
    $\begingroup$ Welcome to the site. Please consider using base R, & commenting it extensively, when illustrating posts here with R code. Not everyone who will come to this page will be familiar with R, & not all of those will be able to read tidy-code. This is a Q&A site for statistics, not R. $\endgroup$ May 16 at 19:47
  • 2
    $\begingroup$ @civilstat Although there is little difference in the β^ estimates here, their confidence intervals differ (+1): This question of mine should illustrate what you are saying: stats.stackexchange.com/questions/646657/… $\endgroup$
    – dariober
    May 17 at 5:34

2 Answers 2

12
$\begingroup$

You’ve got models to two different data sets.

For the Poisson regression, your true conditional expected values (lambda values) are given by $\exp(0.43x+0.2)$.

For the linear regression, your true conditional expected values are given by $0.43x+0.2$.

What you might be more interested in is fitting both models to the same data set, rather than using different $y$ variables.

As far as why the OLS linear model achieves a good fit despite the outcome being Poisson-distributed instead of Gaussian, OLS linear models are rather robust to deviations from normality. If you bootstrap your residuals and calculate the mean, you are likely to find that the distribution looks rather normal.

$\endgroup$
2
  • $\begingroup$ Thank you! The point about robustness of OLS is well taken. I will try your suggestion about bootstrapping residuals to understand that better. I did try fitting each model to each dataset. I think it is just hard to see on the plots and I didn't explain the visualization well. The lines are the predicted values from the fits and there are two lines per plot (one for gaussian and one for poisson), they are just almost perfectly overlapping. $\endgroup$
    – Mikhail
    May 16 at 3:28
  • 1
    $\begingroup$ @Mikhail After reading your comment and revisiting your post, I’m not so sure this answers your question, though if this was helpful, I’m happy to have helped. $\endgroup$
    – Dave
    May 16 at 10:50
6
$\begingroup$

@Dave's answer is a good one (and OP seems to be satisfied with it), but for anyone else finding this question later I'll offer a different perspective.

The TL;DR is that if you think of fitting a GLM and finding the MLE for the coefficients as an optimization problem (which it is), the only real differences that arise when using a different likelihood function are how the model responds to outliers and how overestimates are penalized relative to underestimates.

Any GLM (including OLS, which is just a GLM with a Gaussian likelihood) will try to minimize the sum of the negative log-likelihood across the data points. For Gaussian this will be:

$$n \log\left(\sigma \times \sqrt{2 \pi}\right) + \sum_{i = 1}^{n} \frac{1}{2} \frac{(y_i - \hat{y}_i)^2}{\sigma^2}$$

with $\hat{y}_i = X_i \beta$. For simplicity let's assume $\sigma$ to be known (of course, in practice it isn't, but our choice of $\sigma$ does not affect which $\beta$ minimizes the above expression). If we only look at the minimization with respect to $\beta$, really this is equivalent to minimizing $\sum_{i = 1}^{n} (y_i - \hat{y}_i)^2$, hence the term "least squares."

For Poisson, the sum negative log-likelihood is:

$$\sum_{i = 1}^{n} \left(\hat{y}_i -y_i \log{\hat{y}_i} + \log \left(y_i !\right) \right)$$

Both of these functions are minimized when $\hat{y}_i = y_i$ or $\hat{y}_i - y_i =0$, thus both models are incentivized to return similar coefficients.

The differences come in when $\hat{y}_i - y_i$ moves away from 0. For example, Gaussian likelihood is "symmetric" so it will punish overestimates and underestimates equally, whereas Poisson likelihood will penalize underestimates more severely. Poisson also tends to be more sensitive to outliers, in part because it doesn't have a spread parameter.

But in a situation where the data are roughly symmetric (especially in your log example) and there aren't many outliers, it makes sense that the two fits are very similar.

$\endgroup$

Your Answer

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

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