Linear Regression and Poisson Regression¶

We shall review linear regression in the context of an applied dataset. We will also use this example to motivate Poisson Regression.

Review of Linear Regression with an Economics Dataset¶

The dataset (called MROZ) comes from an Econometrica paper by Mroz in 1987 and contains data on a bunch of variables for married women in the year 1975. This dataset is analyzed in the Econometrics book by Wooldridge (this book contains many other interesting economic and political datasets). The specific variables included in the Mroz dataset are:

1. inlf: binary variable equaling 1 if the individual worked (i.e., they were ‘in the labor force’) in the year 1975 and 0 otherwise

2. hours: number of hours worked in 1975

3. kidslt6: number of kids < 6 years of age

4. kidsge6: number of kids 6-18 years of age

5. age: age in years

6. educ: years of schooling

7. wage: hourly wage in 1975

8. repwage: reported wage at interview in 1976

9. hushrs: hours worked by husband in 1975

10. husage: husband’s age

11. huseduc: husband’s years of schooling

12. huswage: husband’s hourly wage in 1975

13. faminc: family income in 1975

14. mtr: federal marginal tax rate facing woman

15. motheduc: mother’s years of schooling

16. fatheduc: father’s years of schooling

17. unem: unemployment rate in county of residence

18. city: =1 if live in Standard metropolitan statistical area

19. exper: actual labor market experience

20. nwifeinc: (faminc - wage*hours)/1000

21. lwage: log(wage)

22. expersq: $\text{exper}^2$ (the square of the experience variable)

Several regressions can be fit on this dataset. Let us start by fitting a linear regression model with “hours” as the response variable, and “kidslt6”, “kidsge6”, “age”, “educ”, “exper”, “expersq”, “huswage”, “huseduc”, “hushrs”, “motheduc” and “fatheduc” as the covariates (or predictor variables). Note that we are using both “exper” and “expersq” (which is the square of the variable “exper”) as covariates in the regression. This is because we would expect “hours” to increase with “exper” for small values of “exper” (for women new to the workforce) but to perhaps decrease with “exper” for large values of “exper” (for example, people with very large “exper” might be closer to retirement and work less leading to smaller “hours”). Such a relationship cannot be captured by linear functions but can be captured by quadratic functions which is why both “exper” and “expersq” are included in the model.

Using the statsmodels package in Python, this regression can be carried out in the following way.

The table above gives estimates of the regression coefficients and their standard errors for each of the covariates. These estimates and standard errors are computed in the following way. The data is denoted by $(y_i, x_{i1}, x_{i2}, \dots, x_{im})$ for $i = 1, \dots, n$ (here $m$ is the number of covariates and $n$ is the number of observations). One often represents the response data ($y_1, \dots, y_n$) in a $n \times 1$ vector called $Y$, and the covariate data $(x_{i1}, \dots, x_{im})$ for $i = 1, \dots, n$ in a $n \times (m+1)$ matrix $X$ as follows:

Note the additional column of ones in the $X$ matrix. In this notation, the estimates of the regression coefficients are computed by solving equation:

The standard errors are computed in the following way. First one computes the matrix:

The standard errors are given by the square roots of the diagonal entries of the $(m+1) \times (m+1)$ matrix $\hat{\sigma}^2 (X^T X)^{-1}$.

Looking back at the statsmodels regression summary, it is common practice to judge the importance of each variable by the value of the estimated regression coefficient as well as the corresponding standard error. The third column (the column titled t) of the regression output gives the ratio of the estimate of each regression coefficient to the corresponding standard error. Variables for which the magnitude of this ratio is large are considered important. For example, in this example, the variable “kidslt6” has a t value of -7.479 i.e., 7.479 in magnitude which is the largest of all other variables so this can be considered the most important variable among all the covariates. On the other hand, the variable “fatheduc” has t value equal to -0.412 (i.e., 0.412 in absolute value) which is small so this can be considered the least important variable. The right thing to compare these t-values are the quantiles of the Student t probability distribution with $n-m-1$ degrees of freedom. This comparison is done by the OLS function and its values are given in the fourth column of the regression table. This column (titled $P>|t|$) goes by the same ‘p-values’ and values close to zero indicate that the variable is important. For example, for ‘kidslt6’, this $p$-value is rounded to zero while for ‘fathedu’, this $p$-value is quit large equalling 0.681.

Usually one which looks at the table above and drops variables for which the $p$-value $P > |t|$ is large. In this problem, it might be reasonable to drop the variables “motheduc”, “fatheduc”, “hushrs”, “huseduc” and “kidsge6” deeming these to be unimportant for explaining the response “hours”. Below we carry out linear regression on “hours” with the other six variables “kidslt6”, “age”, “educ”, “huswage”, “exper”, “expersq”.

Instead of using the package StatsModels (specifically the function sm.OLS), we can also fit this linear regression model via the Bayesian approach (in PyMC) with flat priors on the regression coefficients (as well as $\log \sigma$). This is done as follows. Note that with these priors, it can be mathematically proved (by calculating the posterior explicitly) that the Bayesian approach gives identical results as the frequentist approach. This means that implementation in PyMC should give essentially the same results as found in the sm.OLS summary (with perhaps slight differences due to Monte Carlo fluctuations).

The above PyMC code specifies the linear regression model in the following way:

Let us check that the Bayesian analysis leads to basically the same answers as statsmodels OLS.

The estimated regression coefficent $\hat{\beta}_1$ for “kidslt6” is about -430 which can be interpreted as follows: having a small child reduces the mean hours worked by about 430. This result is clearly meaningless for women with small working hours (e.g., if someone only works for about 300 hours, what does a reduction of 430 mean?). It is also not very meaningful for women with lots of working hours (e.g., if someone works for more than 3000 hours, a reduction of 430 due to a small kid sounds too small). A much more meaningful and useful interpretation would report reductions in terms of percentages and not absolute hours. For example, a statement such as “having a small child reduces the mean hours worked by about 40%” would apply to all working women and be much more interpretable. Such an intepretation cannot be drawn based on this linear regression.

A percentage interpretation can be realized if we change our model by taking the response variable to be $\log(\text{Hours worked})$ as opposed to “hours”. Indeed, in the linear regression model for $\log(\text{Hours})$:

the coefficient $\beta_1$ gives the reduction in $\log(\text{Hours})$ due to a small child. This implies a percentage increase of $100 \times(\exp(\beta_1) - 1)$ in hours worked because:

It is much more preferable to use $\log(\text{Hours})$ as the response variable in these problems compared to $\text{Hours}$. However, one big problem is that, in this dataset, there are quite a few observations for which the Hours variable equals 0. So we cannot really work with the variable $\log \text{Hours}$. In such case, Poisson regression provides a great alternative model.

Poisson Regression¶

In Poisson regression, one modifies the linear regression model:

(which cannot be used as the variable $\log(Y_i)$ is meaningless when $Y_i = 0$) in the following way:

The main change is that the response variable $Y_i$ is now modelled as $\text{Poisson}(\mu_i)$. The Poisson distribution is often used to model counts. In this case, $Y_i$ is the number of hours worked in 1975 which is a count variable. The linear relationship between the response and covariates is specified through $\log \mu_i$ which enables percentage interpretation for the regression coefficients.

As in the case of linear regression, there are two ways of implementing Poisson regression:

1. Through an in-built function in the library statsmodels. This function is called sm.GLM. GLM stands for Generalized Linear Models which is a general class of models which includes linear regression and Poisson regression (and also logistic regression). While using sm.GLM, we need to specify which GLM we are fitting using the family argument.

2. Using PyMC to do Bayesian inference with flat priors: $\beta_0, \dots, \beta_m \overset{\text{i.i.d}}{\sim} \text{Flat}$.

Most practitioners use the first option above. Let us first see how it works before turning to the Bayesian PyMC approach.

The output above looks very similar to the linear regression output. This model is much more interpretable compared to the linear regression model that we fit earlier. To see this, consider the interpretation of the coefficient -0.8075 for the “kidslt6” variable: having a small kid reduces mean hours worked by 55%. This is clearly a much more interpretable result compared to before.

Maximum Likelihood Estimation in Poisson Regression¶

Let us now understand how sm.GLM is obtaining the estimates and standard errors for the regression coefficients. The estimates are obtained by Maximum Likelihood Estimation. The likelihood in this Poisson Regression Model is given by

It is much more easier to maximize the logarithm of the likelihood (as opposed to the likelihood directly). We can also ignore the $y_i!$ terms in the denominator as they do not involve the parameters $\beta_0, \dots, \beta_m$. Let us also replace the sum $\beta_0 + \beta_1 x_{i1} + \dots + \beta_m x_{im}$ by the simpler expression $x_i^T \beta$ where $\beta$ denotes the column vector with components $\beta_0, \dots, \beta_m$ and $x_i$ denotes the column vector with components $1, x_{i1}, x_{i2}, \dots, x_{im}$. With these changes, we obtain the log-likelihood:

The goal is to maximize this log-likelihood over $\beta_0, \beta_1, \dots, \beta_m$. We can do this by taking the gradient (i.e., derivatives of $\ell(\beta)$ with respect to $\beta_0, \dots, \beta_m$ collected in a column vector), setting the gradient to zero and solving the resulting set of equations. The gradient of $\ell(\beta)$ is

where $\mu$ is the $n \times 1$ vector with entries $\mu_1, \dots, \mu_n$. Therefore setting the gradient equal to zero leads us to the equation:

This gives $m+1$ equations for the $m+1$ parameters $\beta_0, \dots, \beta_m$. Note that the parameters appear in the above equation through $\mu$. It is interesting to observe that in the case of linear regression (where we have the normal distribution instead of Poisson), if we attempt to calculate the MLE in the same way as above, we would obtain exactly the same equation: $X^T \mu = X^T Y$. However, in that case, $\mu = X \beta$ leading to the equation $X^T X \beta = X^T Y$ which can be solved as $\hat{\beta} = (X^T X)^{-1} X^T Y$. Unfortunately, in the present Poisson regression case, $\mu$, which has components $e^{x_1^T \beta}, \dots, e^{x_n^T \beta}$, depends on $\beta$ in a nonlinear fashion and, consequently, it is not possible to solve $X^T \mu = X^T Y$ for $\beta$ in closed form. One needs to resort to solving this equation using numerical methods.

MLE via Newton’s Method¶

A very classical method for solving nonlinear equations is the Newton method (also known as the Newton-Raphson method): see https://en.wikipedia.org/wiki/Newton’s_method. We apply this method to solve $\nabla \ell(\beta) = X^T Y - X^T \mu = 0$. Newton’s method is iterative (starts with a rough guess $\beta^{(0)}$ for the solution and then successively modifies into better approximations for the solution eventually converging to the correct solution). Newton’s interation for changing the current solution approximation $\beta^{(m)}$ into a better approximation $\beta^{(m+1)}$ of the solution to $\nabla \ell(\beta) = 0$ is given by:

where $H \ell(\beta^{(m)}$ is the Hessian Matrix of the function $\ell(\beta)$ evaluated at the point $\beta = \beta^{(m)}$. It may be recalled the Hessian matrix is the matrix of second derivatives (see Hessian matrix). To actually, implement Newton’s method, we would need to calculate the Hessian. We can do this by further differentiating the gradient $\nabla \ell(\beta) = \sum_{i=1}^n (y_i - e^{x_i^T \beta}) x_i = X^T(Y - \mu)$ with respect to $\beta$. It can be checked that this leads to:

We can rewrite this in matrix form as

where $M(\beta)$ is the diagonal matrix with diagonal entries $\mu_1 = e^{x_1^T \beta}, \dots, \mu_n = e^{x_n^T \beta}$. Thus Newton’s method for Poisson regression is:

where we wrote $\mu(\beta^{(m)})$ to emphasize that this vector $\mu$ is calculated as $e^{x_1^T \beta}, \dots, e^{x_n^T \beta}$ with $\beta = \beta^{(m)}$. This method can be implemented as follows.

Because of the exponential calculation in $\mu_i = e^{x_i^T \beta}$, one needs to be careful with initialization. For example, go back and change the initialization to $\beta^{(0)} = (0, \dots, 0)$. In this case, $\beta^{(1)}$ has some large entries leading to $e^{x_i^T \beta^{(1)}}$ blowing up. If this can be avoided, the method converges to the correct solution.

Bayesian Analysis with Flat Prior¶

Before seeing how sm.GLM is calculating the standard errors, let us first look at the Bayesian approach for fitting the Poisson regression model with flat priors. For the Bayesian approach with flat priors, we have

The posterior is therefore given by

This posterior in $\beta_0, \dots, \beta_m$ cannot be written in terms of standard distributions. One can use PyMC to generate samples from this posterior. However, PyMC can be somewhat unstable and might lead to (significantly different) answers for different runs. For example, in the following code, setting randomseed = 0 gives results similar to the frequentist output while setting randomseed = 4 gives quite different results.

Posterior Normal Approximation and Standard Errors¶

When the posterior cannot be written in terms of standard distributions, and when MCMC seems unreliable, an alternative option is to approximate the posterior using simple standard distributions. This idea is at the heart of the method of Variational Inference which is quite popular in Bayesian inference. In this problem, one can approximate the posterior via a (multivariate) normal distribution. This works in the following way.

Our posterior is proportional to $\exp(\ell(\beta))$. The normal distribution, on the other hand, is proportional to $\exp \left(\text{quadratic function of } \beta \right)$. Thus, in order to obtain a normal approximation of the posterior, it makes sense to approximate $\ell(\beta)$ by a quadratic function of $\beta$. This can be done by Taylor approximation but we need to figure out around which point does the Taylor approximation need to be done. The function $\ell(\beta)$ is the log-likelihood which is maximized at the MLE $\hat{\beta}$. It makes sense, therefore, to do the Taylor approximation around $\hat{\beta}$. The idea is that the approximation will be quite accurate for points near the MLE $\hat{\beta}$. For points far from the MLE, the posterior density is probably small anyway so bad approximation for these points might have an insignificant effect. Motivated by these considerations, we Taylor expand the log-likelihood $\ell(\beta)$ around $\hat{\beta}$ up to second order as follows:

Thus the posterior is approximated as:

In the last step above, we ignored $\exp \left(\ell(\hat{\beta}) \right)$ in proportionality because it does not involve the parameters $\beta$.

Now we can compare the posterior approximation to the density of a multivariate normal (with mean $\mu$ and covariance $\Sigma$) which is proportional to $\exp \left(- (\beta - \mu)^T \Sigma^{-1} (\beta - \mu) \right)$. It is then clear that $\mu = \hat{\beta}$ and $\Sigma^{-1} = -H\ell(\hat{\beta})$. Thus the posterior normal approximation is given by:

This gives a very simple and practical approximation to the complicated posterior, especially because we already know how to calculate the Hessian (it was used in the Newton algorithm for computing the MLE). If one wants, posterior samples can be easily generated from $N \left(\hat{\beta}, \left( - H\ell(\hat{\beta}) \right)^{-1}\right)$. This will be much faster than using PyMC. However, PyMC tries to obtain samples from the actual posterior while this normal distribution is an approximation to the actual posterior.

Standard Errors for Poisson Regression¶

The standard errors of Poisson Regression that were reported by sm.GLM are actually computed from the covariance of the posterior normal approximation. In other words, the standard errors are simply the square-roots of the diagonal entries of the matrix:

where $\hat{\beta}$ is the MLE and $M(\beta)$ is the diagonal matrix with diagonal entries $\mu_1 = e^{x_1^T \beta}, \dots, \mu_n = e^{x_n^T \beta}$. This can be readily verified in the context of our dataset as follows.

These standard errors can also be justified using frequentist arguments but that justification is much more complicated (it involves asymptotics and notions such as Fisher information). The Bayesian justification is much simpler and intuitive (based on normal approximation via Taylor expansion).

The important takeaways are:

1. Poisson Regression is an important alternative to linear regression which is applicable when the response variable involves counts. It gives more interpretable results than linear regression by allowing one to do log-linear modeling in situations where one cannot directly take logarithms of the response variable due to existence of zero counts.

2. In the case of linear regression, frequentist inference and Bayesian inference with flat priors coincide exactly as can be rigorously mathematically proved.

3. For Poisson regression, frequentist inference uses the MLE $\hat{\beta}$ and provides standard errors by taking square roots of the diagonal entries of the matrix $\left(-H\ell(\hat{\beta}) \right)^{-1}$. Bayesian inference with flat priors can give different answers to frequentist inference in Poisson regression. However, if one implements Bayesian inference by not working with the full posterior and instead approximating the posterior by a multivariate normal distribution centered at the MLE, then Bayesian inference gives identical answers to the frequentist inference.