StatsModels

Section 11.1 StatsModels

In the previous chapter we used the SciPy function linregress to compute least squares fits. This function performs simple regression, but not multiple regression. For that, we’ll use StatsModels, a package that provides several forms of regression and other analyses.

🔗

As a first example, we’ll continue exploring the penguin data.

🔗

The following cell downloads the data from a repository created by Allison Horst.

🔗

Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. R package version 0.1.0. https://allisonhorst.github.io/palmerpenguins/. doi: 10.5281/zenodo.3960218.

🔗

The data was collected as part of the research that led to this paper: Gorman KB, Williams TD, Fraser WR (2014). Ecological sexual dimorphism and environmental variability within a community of Antarctic penguins (genus Pygoscelis). PLoS ONE 9(3):e90081. https://doi.org/10.1371/journal.pone.0090081

🔗

Listing 11.1.1. Python Code

🔗

download(
    "https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
)

When we load the data, we’ll use the following dictionary to give the columns names that don’t contain spaces, which will make them easier to use with StatsModels.

🔗

Listing 11.1.2. Python Code

🔗

columns = {
    "Body Mass (g)": "mass",
    "Flipper Length (mm)": "flipper_length",
    "Culmen Length (mm)": "culmen_length",
    "Culmen Depth (mm)": "culmen_depth",
}

Now we can load the data, drop the rows with missing mass, and rename the columns.

🔗

Listing 11.1.3. Python Code

🔗

$ penguins = (
    pd.read_csv("penguins_raw.csv")
    .dropna(subset=["Body Mass (g)"])
    .rename(columns=columns)
)
penguins.shape
(342, 17)

The dataset contains three species of penguins. We’ll work with just the Adélie penguins.

🔗

Listing 11.1.4. Python Code

🔗

$ adelie = penguins.query('Species.str.startswith("Adelie")').copy()
len(adelie)
151

In the previous chapter, we computed a least squares fit between the penguins’ flipper lengths and weights.

🔗

Listing 11.1.5. Python Code

🔗

flipper_length = adelie["flipper_length"]
body_mass = adelie["mass"]

As a reminder, here’s how we did that with linregress.

🔗

Listing 11.1.6. Python Code

🔗

$ from scipy.stats import linregress

result_linregress = linregress(flipper_length, body_mass)
result_linregress.intercept, result_linregress.slope
(np.float64(-2535.8368022002514), np.float64(32.83168975115009))

StatsModels provides two interfaces (APIs) -- we’ll use the "formula" API, which uses the Patsy formula language to specify the response and explanatory variables. The following formula string specifies that the response variable, mass, is a linear function of one explanatory variable, flipper_length.

🔗

Listing 11.1.7. Python Code

🔗

formula = "mass ~ flipper_length"

We can pass this formula to the StatsModels function ols, along with the data.

🔗

Listing 11.1.8. Python Code

🔗

$ import statsmodels.formula.api as smf

model = smf.ols(formula, data=adelie)
type(model)
statsmodels.regression.linear_model.OLS

The name ols stands for "ordinary least squares", which indicates that this function computes a least squares fit under the most common, or "ordinary", set of assumptions.

🔗

The result is an OLS object that represents the model. In general, a model is a simplified description of the relationship between variables. In this example, it’s a linear model, which means that it assumes that the response variable is a linear function of the explanatory variables.

🔗

The fit method fits the model to the data and returns a RegressionResults object that contains the result.

🔗

Listing 11.1.9. Python Code

🔗

result_ols = model.fit()

Listing 11.1.10. Python Code

🔗

$ # Technically it's a RegressionResultsWrapper

type(result_ols)
statsmodels.regression.linear_model.RegressionResultsWrapper

The RegressionResults object contains a lot of information, so thinkstats provides a function that just displays the information we need for now.

🔗

Listing 11.1.11. Python Code

🔗

from thinkstats import display_summary

display_summary(result_ols)

Table 11.1.12.

🔗

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	-2535.8368	964.798	-2.628	0.009	-4442.291	-629.382
flipper_length	32.8317	5.076	6.468	0.000	22.801	42.862

Table 11.1.13.

🔗

R-squared:

0.2192

The first column contains the intercept and slope, which are the coefficients of the model. We can confirm that they are the same as the coefficients we got from linregress.

🔗

Listing 11.1.14. Python Code

🔗

$ result_linregress.intercept, result_linregress.slope
(np.float64(-2535.8368022002514), np.float64(32.83168975115009))

The second column contains the standard errors of the coefficients -- again, they are the same as the values we got from linregress.

🔗

Listing 11.1.15. Python Code

🔗

$ result_linregress.intercept_stderr, result_linregress.stderr,
(np.float64(964.7984274994059), np.float64(5.076138407990821))

The next column reports \(t\) statistics, which are used to compute p-values -- we can ignore them, because the p-values are in the next column, labeled P>|t|. The p-value for flipper_length is rounded down to 0, but we can display it like this.

🔗

Listing 11.1.16. Python Code

🔗

$ result_ols.pvalues["flipper_length"]
np.float64(1.3432645947789321e-09)

And confirm that linregress computed the same result.

🔗

Listing 11.1.17. Python Code

🔗

$ result_linregress.pvalue
np.float64(1.3432645947790051e-09)

The p-value is very small, which means that if there were actually no relationship between weight and flipper length, it is very unlikely we would see a slope as big as the estimated value by chance.

🔗

The last two columns, labeled [0.025 and 0.975], report 95% confidence intervals for the intercept and slope. So the 95% CI for the slope is [22.8, 42.9].

🔗

The last line reports the \(R^2\) value of the model, which is about 0.22 -- that means we can reduce MSE by about 22% if we use flipper length to predict weight, compared to using only the average weight.

🔗

The \(R^2\) value we get from simple correlation is the square of the correlation coefficient, \(r\text{.}\) So we can compare rsquared computed by ols with the square of the rvalue computed by linregress.

🔗

Listing 11.1.18. Python Code

🔗

$ result_ols.rsquared, result_linregress.rvalue**2
(np.float64(0.21921282646854878), np.float64(0.21921282646854875))

They are the same except for a small difference due to floating-point approximation.

🔗

Before we go on to multiple regression, let’s compute one more simple regression, with culmen_length as the explanatory variable (the culmen is the top ridge of the bill).

🔗

Listing 11.1.19. Python Code

🔗

formula = "mass ~ culmen_length"
result = smf.ols(formula, data=adelie).fit()
display_summary(result)

Table 11.1.20.

🔗

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	34.8830	458.439	0.076	0.939	-870.998	940.764
culmen_length	94.4998	11.790	8.015	0.000	71.202	117.798

Table 11.1.21.

🔗

R-squared:

0.3013

Again, the p-value of the slope is very small, which means that if there were actually no relationship between mass and culmen length, it is unlikely we would see a slope this big by chance. You might notice that the p-value associated with the intercept is large, but that’s not a problem because we are not concerned about whether the intercept might be zero. In this model, the intercept is close to zero, but that’s just a coincidence -- it doesn’t indicate a problem with the model.

🔗

The \(R^2\) value for this model is about 0.30, so the reduction in MSE is a little higher if we use culmen length rather than flipper length as an explanatory variable (the \(R^2\) value with flipper length is 0.22). Now, let’s see what happens if we combine them.

🔗

Prev Top Next