Seeing if a Poisson Regression follows the necessary assumptions

Question

I am trying to understand how simulation plays a role in checking model assumptions when the residuals do not have exact distributions.

I took this Poisson Regression Model:

$$\lambda_i = e^{(\beta_0 + \beta_1X_i)}$$

$$P(Y_i | X_i) \sim Poisson(\lambda_i = e^{(\beta_0 + \beta_1X_i)})$$

How can I see if the model assumptions are being met?

Some ideas come to mind:

1. Approach 1: I can bootstrap the data, fit the model on each bootstrap sample. If the model assumptions are met - at each $x_i$: the distribution of residuals over all bootstrap samples will have a mean of 0 and a variance of $\lambda_i$. (this approach seems like the most work computationally since multiple models are being fit)

2. Approach 2: I fit the model on the original data, then simulate new datasets from the model ... and compare all real $Y$ vs the simulated $Y$. If the model assumptions are met, the model will be able to very closely "reproduce" the observed data. (I don't know how to measure "closeness" over all $x_i$) . This approach seems like less work since only one model is being fit.

3. Approach 3 (DHARMA approach https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html): I fit the model on the original data. After fitting the model, I generate multiple values of $y_i$ at each $x_i$. For a given $x_i$ , I make an Empirical CDF of the simulated $y_i$. For the actual observed $y_i$ at that $x_i$, I find out where this $y_i$ is situated on the CDF. If the model fits the data well, this $y_i$ should be at the 50% level (I am not sure why - is this because that residuals have a uniform distribution?). I then repeat this empirical CDF for each $x_i$ and see on average if the true $y_i$ are close to the 50% level in each CDF on average. (this approach also only involves fitting one model)

All 3 approaches correct? Is this the correct use of simulation to validate model assumptions when the residuals do not have exact distributions?

Answer 1

In addition to the resampling based methods you suggested, let me offer you the "traditional" old-fashined methods. For Poisson regression we write $y_i\sim\text{Poisson}(\lambda_i)$ with $E[y_i|x_i]=\lambda_i=\exp(\beta^Tx_i+\epsilon_i)$. We have the following model assumptions:

1. Error independence ($\epsilon_i\perp\epsilon_j$)

2. Error normality ($\hat\beta$ has multivariate normal distribution)

3. Unbiasedness ($E[\epsilon]=0$)

4. Homoskedasticity ($Var(\epsilon)=const.$)

5. Poisson property ($E[\hat y]=Var(\hat y)$)

Assumptions 1-4 are standard for all GLMs. Unlike linear regression, we have to define here standardized residuals as (see Section 6.2 of Gelman & Hill)

$$z_i=\frac{y_i-\lambda_i}{\sqrt{\lambda_i}}$$

With these at hand, you can go for Q-Q plots and all traditional diagnostics to validate the error assumptions (see here). The only special thing here is assumption 5. The case where $Var(\hat y)$ is too large is called overdispersion, and we diagnose it using a $\chi^2$ deviance statistic:

$$\text{Est. Overdispersion}=\frac{1}{n-p}\sum_{i=1}^{n}z_i^2$$

with $n$ the number of samples and $p$ the length of the vector $\beta$. The above mean is assumed to follow the $\chi^2_{n-p}$ distribution, and so a simple $\chi^2$ test will help you here. Of course, you have to keep in mind that like all other goodness-of-fit tests, the test tends to reject to null when $n$ is large.

As per additional resampling based methods, I recommend taking a look at Wickham's Lineup method, which might become handy here.