Forecasting count data after fitting a Poisson regression model

Question

I have created a Poisson regression model in Python for predicting the number of orders a company will receive a day based on a dataset of order counts and a number (~5) factors that contribute to those order counts. After some fine tuning of the model, I am happy that it is fitting the data well.

I know that the steps I want to take from here are that of being able to forecast the counts in the future, and after some research I am unsure how to convert what I have, into something that can do this. Ideally I would like to produce a model that can forecast for a nominal amount of days into the future, basing the predictions from the linear data of days of the week and the expected values of the other factors, such as website clicks, current market size/growth etc.

Answer 1

What you are looking for is a prediction interval for a GLM. The prediction interval should take account of the inherent uncertainty in the quantity being predicted, but also the uncertainty in the coefficient estimator used to estimate the distribution of the former quantity. There are various ways of constructing this type of interval, depending on how sophisticated you want to be. The essence of these methods is to "integrate" the distribution of the quantity of interest (in this case a Poisson distribution for the count variable) over the uncertainty in the model coefficients. To do this, we would usually apply the delta method and take the normal approximation to the distribution of the coefficient estimator in the model.

---

I'll give you an example of a fairly crude method for getting a prediction interval. Given the standard WLS estimator $\hat{\boldsymbol{\beta}} \in \mathbb{R}^m$ with estimated covariance matrix $\hat{\mathbf{\Sigma}}_\boldsymbol{\beta} \in \mathbb{R}^{m \times m}$ we can take the approximate distribution $\hat{\boldsymbol{\beta}} - \boldsymbol{\beta} \sim \text{N}(\mathbf{0}, \hat{\mathbf{\Sigma}}_\boldsymbol{\beta})$ for the error in the coefficient estimator. We can use this to estimate the predictive distribution using the expression:

$$\hat{f} (Y_i=y | \mathbf{x}_i) = \int \limits_\mathbb{{R}^m} \text{Pois}(y | \exp(\mathbf{x}_i^\text{T} \boldsymbol{\beta})) \cdot \text{N}(\boldsymbol{\beta} | \hat{\boldsymbol{\beta}}, \hat{\mathbf{\Sigma}}_\boldsymbol{\beta}) \ d \boldsymbol{\beta}.$$

This estimmator is fairly crude, but it gives the basic idea of how you estimate the predictive distribution. It takes account of the uncertainty in the coefficient estimator for $\boldsymbol{\beta}$, but it does not take account of uncertainty in the estimator $\hat{\mathbf{\Sigma}}_\boldsymbol{\beta}$. (You can generalise to use a multivariate T-distribution to take some account of this additional uncertainty if you want. More sophisticated methods also exist for specific models, but it is a big field.) Since this integral expression does not have a closed-form solution, usually you would simulate this distribution by generating random values for the true parameter and using a Monte-Carlo estimate of the integral. With $K$ simulations you could use:

$$\hat{f} (Y_i=y | \mathbf{x}_i) \approx \frac{1}{K} \sum_{k = 1}^K \text{Pois}(y | \exp(\mathbf{x}_i^\text{T} \boldsymbol{\beta}_k)) \quad \quad \quad \boldsymbol{\beta}_1,...,\boldsymbol{\beta}_K \sim \text{IID N}(\hat{\boldsymbol{\beta}}, \hat{\mathbf{\Sigma}}_\boldsymbol{\beta}).$$

You can then form prediction intervals from the estimated predictive distribution for the response variable. You can find some more sophisticated analysis for the Poisson GLM in Wood (2005), though this also uses some conservative methods.

There are various packages in Python and R that implement various types of prediction intervals for GLMs including the Poisson regression model. I am not an expert in Python so my knowledge of particular packages and their underlying methods is thin. In any case, if you look up prediction intervals for GLMs in that program you should be able to find some packages to implement this. Alternatively, you can program crude methods from scratch fairly simply from the output of the Poisson regression model.

Answer 2

The standard Poisson regression models observations as Poisson distributed with a parameter $\lambda=e^{x\beta}$ for a predictor row vector $x$. Python has presumably given you parameter estimates $\hat{\beta}$.

So you can forecast a future realization's Poisson parameter $\hat{\lambda}=e^{x\hat{\beta}}$, where $x$ is of course the corresponding vector of predictors, as translated into a design matrix row (dummy coding etc.).

Now, $\hat{\lambda}$ is also the expected value of your future realization. If you want a median or other quantile forecast (e.g., for setting safety amounts), you can just extract the quantiles from the $\text{Pois}(\hat{\lambda})$ distribution.

(Note that this approach completely disregards the sampling variability in $\hat{\beta}$, but it is almost always done this way, and usually this makes no major difference.)