You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: myst_nbs/generalized_linear_models/GLM-robust.myst.md
+108-43
Original file line number
Diff line number
Diff line change
@@ -6,47 +6,56 @@ jupytext:
6
6
format_version: 0.13
7
7
jupytext_version: 1.13.7
8
8
kernelspec:
9
-
display_name: Python 3
9
+
display_name: Python 3 (ipykernel)
10
10
language: python
11
11
name: python3
12
12
---
13
13
14
+
(GLM-robust)=
14
15
# GLM: Robust Linear Regression
15
16
16
-
This tutorial first appeard as a post in small series on Bayesian GLMs on:
17
+
:::{post} August, 2013
18
+
:tags: regression, linear model, robust
19
+
:category: beginner
20
+
:author: Thomas Wiecki, Chris Fonnesbeck, Abhipsha Das, Conor Hassan, Igor Kuvychko, Reshama Shaikh, Oriol Abril Pla
21
+
:::
17
22
18
-
1.[The Inference Button: Bayesian GLMs made easy with PyMC3](http://twiecki.github.com/blog/2013/08/12/bayesian-glms-1/)
19
-
2.[This world is far from Normal(ly distributed): Robust Regression in PyMC3](http://twiecki.github.io/blog/2013/08/27/bayesian-glms-2/)
20
-
3.[The Best Of Both Worlds: Hierarchical Linear Regression in PyMC3](http://twiecki.github.io/blog/2014/03/17/bayesian-glms-3/)
23
+
+++
24
+
25
+
# GLM: Robust Linear Regression
26
+
27
+
The tutorial is the second of a three-part series on Bayesian *generalized linear models (GLMs)*, that first appeared on [Thomas Wiecki's blog](https://twiecki.io/):
28
+
29
+
1. {ref}`Linear Regression <pymc:GLM_linear>`
30
+
2. {ref}`Robust Linear Regression <GLM-robust>`
31
+
3. {ref}`Hierarchical Linear Regression <GLM-hierarchical>`
21
32
22
33
In this blog post I will write about:
23
34
24
35
- How a few outliers can largely affect the fit of linear regression models.
25
36
- How replacing the normal likelihood with Student T distribution produces robust regression.
26
-
- How this can easily be done with the `Bambi` library by passing a `family` object or passing the family name as an argument.
27
-
28
-
This is the second part of a series on Bayesian GLMs (click [here for part I about linear regression](http://twiecki.github.io/blog/2013/08/12/bayesian-glms-1/)). In this prior post I described how minimizing the squared distance of the regression line is the same as maximizing the likelihood of a Normal distribution with the mean coming from the regression line. This latter probabilistic expression allows us to easily formulate a Bayesian linear regression model.
37
+
38
+
In the {ref}`linear regression tutorial <pymc:GLM_linear>` I described how minimizing the squared distance of the regression line is the same as maximizing the likelihood of a Normal distribution with the mean coming from the regression line. This latter probabilistic expression allows us to easily formulate a Bayesian linear regression model.
29
39
30
40
This worked splendidly on simulated data. The problem with simulated data though is that it's, well, simulated. In the real world things tend to get more messy and assumptions like normality are easily violated by a few outliers.
31
41
32
42
Lets see what happens if we add some outliers to our simulated data from the last post.
33
43
34
44
+++
35
45
36
-
Again, import our modules.
46
+
First, let's import our modules.
37
47
38
48
```{code-cell} ipython3
39
49
%matplotlib inline
40
50
51
+
import aesara
52
+
import aesara.tensor as at
41
53
import arviz as az
42
-
import bambi as bmb
43
54
import matplotlib.pyplot as plt
44
55
import numpy as np
45
56
import pandas as pd
46
-
import pymc3 as pm
47
-
import theano
48
-
49
-
print(f"Running on pymc3 v{pm.__version__}")
57
+
import pymc as pm
58
+
import xarray as xr
50
59
```
51
60
52
61
```{code-cell} ipython3
@@ -79,7 +88,7 @@ data = pd.DataFrame(dict(x=x_out, y=y_out))
79
88
Plot the data together with the true regression line (the three points in the upper left corner are the outliers we added).
80
89
81
90
```{code-cell} ipython3
82
-
fig = plt.figure(figsize=(7, 7))
91
+
fig = plt.figure(figsize=(7, 5))
83
92
ax = fig.add_subplot(111, xlabel="x", ylabel="y", title="Generated data and underlying model")
Lets see what happens if we estimate our Bayesian linear regression model using the `bambi`. This function takes a [`formulae`](https://bambinos.github.io/formulae/api_reference.html) string to describe the linear model and adds a Normal likelihood by default.
102
+
Lets see what happens if we estimate our Bayesian linear regression model by fitting a regression model with a normal likelihood.
103
+
Note that the bambi library provides an easy to use such that an equivalent model can be built using one line of code.
104
+
A version of this same notebook using Bambi is available at {doc}`bambi's docs <bambi:notebooks/t_regression>`
To evaluate the fit, I am plotting the posterior predictive regression lines by taking regression parameters from the posterior distribution and plotting a regression line for each (this is all done inside of `plot_posterior_predictive()`).
124
+
To evaluate the fit, the code below calculates the posterior predictive regression lines by taking regression parameters from the posterior distribution and plots a regression line for 20 of them.
plt.title("Posterior predictive for normal likelihood");
108
136
```
109
137
110
138
As you can see, the fit is quite skewed and we have a fair amount of uncertainty in our estimate as indicated by the wide range of different posterior predictive regression lines. Why is this? The reason is that the normal distribution does not have a lot of mass in the tails and consequently, an outlier will affect the fit strongly.
111
139
112
140
A Frequentist would estimate a [Robust Regression](http://en.wikipedia.org/wiki/Robust_regression) and use a non-quadratic distance measure to evaluate the fit.
113
141
114
-
But what's a Bayesian to do? Since the problem is the light tails of the Normal distribution we can instead assume that our data is not normally distributed but instead distributed according to the [Student T distribution](http://en.wikipedia.org/wiki/Student%27s_t-distribution) which has heavier tails as shown next (I read about this trick in ["The Kruschke"](http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/), aka the puppy-book; but I think [Gelman](http://www.stat.columbia.edu/~gelman/book/) was the first to formulate this).
142
+
But what's a Bayesian to do? Since the problem is the light tails of the Normal distribution we can instead assume that our data is not normally distributed but instead distributed according to the [Student T distribution](http://en.wikipedia.org/wiki/Student%27s_t-distribution) which has heavier tails as shown next {cite:p}`gelman2013bayesian,kruschke2014doing`.
115
143
116
144
Lets look at those two distributions to get a feel for them.
As you can see, the probability of values far away from the mean (0 in this case) are much more likely under the `T` distribution than under the Normal distribution.
130
158
131
-
To define the usage of a T distribution in `Bambi` we can pass the distribution name to the `family` argument -- `t` -- that specifies that our data is Student T-distributed. Note that this is the same syntax as `R` and `statsmodels` use.
159
+
Below is a PyMC model, with the `likelihood` term following a `StudentT` distribution with $\nu=3$ degrees of freedom, opposed to the `Normal` distribution.
plt.title("Posterior predictive for Student-T likelihood")
145
189
```
146
190
147
191
There, much better! The outliers are barely influencing our estimation at all because our likelihood function assumes that outliers are much more probable than under the Normal distribution.
@@ -150,20 +194,41 @@ There, much better! The outliers are barely influencing our estimation at all be
150
194
151
195
## Summary
152
196
153
-
-`Bambi` allows you to pass in a `family` argument that contains information about the likelihood.
154
197
- By changing the likelihood from a Normal distribution to a Student T distribution -- which has more mass in the tails -- we can perform *Robust Regression*.
155
198
156
-
The next post will be about logistic regression in PyMC3 and what the posterior and oatmeal have in common.
157
-
158
199
*Extensions*:
159
200
160
201
- The Student-T distribution has, besides the mean and variance, a third parameter called *degrees of freedom* that describes how much mass should be put into the tails. Here it is set to 1 which gives maximum mass to the tails (setting this to infinity results in a Normal distribution!). One could easily place a prior on this rather than fixing it which I leave as an exercise for the reader ;).
161
-
- T distributions can be used as priors as well. I will show this in a future post on hierarchical GLMs.
162
-
- How do we test if our data is normal or violates that assumption in an important way? Check out this [great blog post](http://allendowney.blogspot.com/2013/08/are-my-data-normal.html) by Allen Downey.
202
+
- T distributions can be used as priors as well. See {ref}`GLM-hierarchical`
203
+
- How do we test if our data is normal or violates that assumption in an important way? Check out this [great blog post](http://allendowney.blogspot.com/2013/08/are-my-data-normal.html) by Allen Downey.
0 commit comments