Regression model for integer response

Question

Let the response be $Y_i \in \mathbb{Z}$ and the covariate $X_i \in \mathbb{R}^p$. For counting data where $Y_i$ are restricted to be nonnegative, we have Poisson regression or negative binomial regression. Now, I consider the situation where the response is a kind of integral score which can range from negative to positive integers, and that we don't know the bounds for the response. For example, the number of games that a person won in the past ten years, which can be negative if he/she lost more than won.

Is it still reasonable to use linear regression with normal error? If not, any alternatives should I be looking for? Thanks in advance.

Answer 1

As an alternative we would try using a Skellam distribution as our response distribution. (i.e. our outcome variable is the difference between two Poisson-distributed random variables). It is sometimes used to predict football/rugby goal difference. Do note though that maybe more straightforward to just have to Poisson models and compute their difference. SE.SO has a relevant thread on: How can I fit a Skellam regression?.

Answer 2

Ideally, you would use a regression model with a discrete distribution with support over the integers for your response variable. This could be approximated by a continuous distribution, so long as the standard error in the model (i.e., the standard deviation of the error term) is not too small compared to the unit interval between integers (so that the continuous density is not changing much between integers).

If you want to, you can take a regression model that uses a continuous response distribution $f$ (e.g., from the Gaussian linear model) and then "discretise" the response distribution by taking:

$$\mathbb{P}(Y_i = y | \mathbf{x}_i, \boldsymbol{\beta}) = \int \limits_{y-1/2}^{y+1/2} f(r | \mathbf{x}_i, \boldsymbol{\beta}) \ dr.$$

For example, in the Gaussian linear model, instead of having the likelihood function:

$$\begin{align} L_{\mathbf{y}, \mathbf{x}}(\boldsymbol{\beta}) &= \prod_{i=1}^n \phi \bigg( \frac{y_i - \mathbf{x}_i \cdot \boldsymbol{\beta}}{\sigma} \bigg), \end{align}$$

you would instead have the "discretised" version:

$$\begin{align} L_{\mathbf{y}, \mathbf{x}}(\boldsymbol{\beta}) &= \prod_{i=1}^n \bigg[ \Phi \bigg( \frac{y_i + \tfrac{1}{2} - \mathbf{x}_i \cdot \boldsymbol{\beta}}{\sigma} \bigg) - \Phi \bigg( \frac{y_i - \tfrac{1}{2} - \mathbf{x}_i \cdot \boldsymbol{\beta}}{\sigma} \bigg) \bigg] \end{align}$$

The MLE in both cases is going to be quite similar, so long as $\sigma$ is substantially bigger than one. The main drawback of the discretised version is that the MLE is no longer the OLS estimator and does not have a closed form, so the resulting theoretical properties are a bit trickier to deal with (but certainly not impossible).