Exercises for Section 4.1 Simple linear regression

Exercises for Section 4.1 Simple linear regression#

This notebook contains the exercises from Section 4.1 Simple linear regression in the No Bullshit Guide to Statistics.

Notebook setup#

# Ensure required Python modules are installed
%pip install --quiet numpy scipy seaborn pandas ministats

[notice] A new release of pip is available: 26.2.1 -> 26.2.2
[notice] To update, run: pip install --upgrade pip

Note: you may need to restart the kernel to use updated packages.

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Figures setup
plt.clf()  # needed otherwise `sns.set_theme` doesn't work
sns.set_theme(
    context="paper",
    style="whitegrid",
    palette="colorblind",
    rc={"figure.figsize": (6, 1.6)},
)
%config InlineBackend.figure_format = "retina"

<Figure size 640x480 with 0 Axes>

# Simplified int and float __repr__
np.set_printoptions(legacy='1.25')

# Check the `ministats` version to make sure all the datasets will be in place.
import ministats
from packaging.version import Version

if Version(ministats.__version__) <= Version("0.5.19"):
    print("You have version", ministats.__version__, "of the `ministats` library.")
    print("You must use ministats version 0.5.20 or later.")
    print("Run the following command in a new notebook cell:")
    print("   %pip install --upgrade ministats")

You have version 0.5.7 of the `ministats` library.
You must use ministats version 0.5.20 or later.
Run the following command in a new notebook cell:
   %pip install --upgrade ministats

# %pip install --upgrade ministats

# Download datasets/ directory if necessary
from ministats import ensure_datasets
ensure_datasets()

datasets/ directory already exists.

$\def\stderr#1{\mathbf{se}_{#1}}$ $\def\stderrhat#1{\hat{\mathbf{se}}_{#1}}$ $\newcommand{\Mean}{\textbf{Mean}}$ $\newcommand{\Var}{\textbf{Var}}$ $\newcommand{\Std}{\textbf{Std}}$ $\newcommand{\Freq}{\textbf{Freq}}$ $\newcommand{\RelFreq}{\textbf{RelFreq}}$ $\newcommand{\DMeans}{\textbf{DMeans}}$ $\newcommand{\Prop}{\textbf{Prop}}$ $\newcommand{\DProps}{\textbf{DProps}}$ $$ \newcommand{\CI}[1]{\textbf{CI}_{#1}} \newcommand{\CIL}[1]{\textbf{L}_{#1}} \newcommand{\CIU}[1]{\textbf{U}_{#1}} \newcommand{\ci}[1]{\textbf{ci}_{#1}} \newcommand{\cil}[1]{\textbf{l}_{#1}} \newcommand{\ciu}[1]{\textbf{u}_{#1}} $$ (this cell contains the macro definitions like $\stderr{\overline{\mathbf{x}}}$, $\stderrhat{}$, $\Mean$, …)

Exercises#

E4.1#

Use statsmodels to fit a linear model to the doctors dataset (datasets/doctors.csv) based on the formula "score 1 + alc". Interpret the estimated slope parameter in your own words.

import pandas as pd
doctors = pd.read_csv("datasets/doctors.csv")
doctors.head(3)

	permit	loc	work	hours	caf	alc	weed	exrc	score
0	93273	rur	hos	21	2	0	5.0	0.0	63
1	90852	urb	cli	74	26	20	0.0	4.5	16
2	92744	urb	hos	63	25	1	0.0	7.0	58

import statsmodels.formula.api as smf
...

Ellipsis

E4.2#

Consider the following dataset with six observations $[\mathbf{x},\mathbf{y}]=$ $(105, 45)$, $(113,63)$, $(125,86)$, $(137,118)$, $(141,112)$, $(153,169)$. We’re interested in the linear model $Y|x \sim \mathcal{N}(\beta_0 + \beta_1x, \; \sigma)$ for this data.

a) Use the formulas for $\widehat{\beta}_{0}$ and $\widehat{\beta}_{1}$ from page to estimate the intercept $\beta_0$ and slope $\beta_1$ model parameters.

b) Use the formula from page to estimate the model’s standard deviation parameter $\sigma$.

c) Use statsmodels to fit a linear model to this dataset. Verify your answers to a) and b) by accessing the appropriate attributes of the fitted model.

Hint: For part c), you can “package” the data $\mathbf{x}=\texttt{xs}$ and $\mathbf{y}=\texttt{ys}$ as a data frame object using df_xys=pd.DataFrame({"x":xs,"y":ys}).

# Data
xs = [105, 113, 125, 137, 141, 153]
ys = [45, 63, 86, 118, 112, 169]

E4.3#

You’re looking to rent an apartment. To get an idea of how rent depends on apartment size, you check the classifieds and find the following offers: a 200 sq ft studio for $\$500$, a 300 sq ft studio for $\$620$, a 400 sq ft small apartment for $\$750$, a 500 sq ft one-bedroom condo for $\$890$, a 900 sq ft two-bedroom apartment for $\$1250$, and a 1000 sq ft apartment for $\$1300$. Find the best-fitting linear model $p(a) = \beta_0 + \beta_1 a$ for the above data, where $p(a)$ represents the rent price for an apartment with $a$ square feet of livable area. What is the predicted price for a 700 sq ft apartment?

areas = [200, 300, 400, 500, 900, 1000]    # in square feet
prices = [500, 620, 750, 890, 1250, 1300]  # in dollars

# Package `areas` and `prices` as a data frame
df_apts = pd.DataFrame({ ... })

# Fit the linear model
import statsmodels.formula.api as smf
...

Ellipsis

E4.4#

Repeat the calculations to obtain the regression result object lm1 for the linear model $S|e \sim \mathcal{N}(\beta_0 + \beta_1 e, \sigma)$, which describes the students’ scores $S$ as a function of their effort $e$.

a) Use the method lm1.get_prediction to predict the mean score for a student that invests $10$ hours of $\texttt{effort}$, $\mu_{\widehat{S}}(10)$.

b) Construct a 90% confidence interval $\ci{\mu_{S},0.9}(10)$ for the mean score $\mu_{S}(10)$ by calling the .conf_int method.

c) Use lm1.get_prediction to predict the score of a student that invests $10$ hours of $\texttt{effort}$, $\hat{s}(10)$.

d) Report a 90% prediction interval $\ci{S,0.9}(10)$ for the score.

Hint: Use the code examples from Section 4.1 as a reference.

# Load the students dataset
import pandas as pd
students = pd.read_csv("datasets/students.csv")

# Fit the linear model
import statsmodels.formula.api as smf
lm1 = smf.ols("score ~ 1 + effort", data=students).fit()
lm1.params

Intercept    32.465809
effort        4.504850
dtype: float64

E4.5#

The OECD Better Life Index (BLI) dataset (datasets/oecd_bli.csv) contains information about life satisfaction collected from different countries around the world. We want to know if being rich, as measured by the gross domestic product per capita, as measured by the variable gdp, is linked to higher life satisfaction as measured by the variable lifesat. Use simple linear regression to model this relationship. Generate a scatter plot of the data and the regression line.

Hint: You can use the ministats function plot_reg for the plot.

oecd_bli = pd.read_csv("datasets/oecd_bli.csv")
oecd_bli.head()

	Country	educ	income	lifesat	salary	gdp
0	Mexico	42.0	16269.0	6.0	16230.0	23659.0
1	Greece	76.0	20791.0	5.8	27207.0	38397.0
2	Slovak Republic	92.0	21149.0	6.5	23619.0	40586.0
3	Hungary	86.0	21026.0	6.0	25409.0	43476.0
4	Portugal	55.0	24877.0	5.8	28410.0	44963.0

import statsmodels.formula.api as smf
...

Ellipsis

E4.6#

The marketing dataset (datasets/marketing.csv) contains information about advertising spending on different platforms and the resulting sales figures of a business.

a) Fit a model for the relationship between youtube ad spend and sales.

b) Use the ministats function plot_resid to plot the residuals.

c) Comment on the modelling assumptions based on your observations from the residuals plot you generated in part b).

marketing = pd.read_csv("datasets/marketing.csv")
marketing.head(3)

	youtube	facebook	newspaper	sales
0	276.12	45.36	83.04	26.52
1	53.40	47.16	54.12	12.48
2	20.64	55.08	83.16	11.16

# a)
import statsmodels.formula.api as smf
...

Ellipsis

# b)
from ministats import plot_resid
...

Ellipsis

E4.7#

The function get_rdataset defined in statsmodels.datasets allows us to downlaod datasets provided with R packages as pandas data frames. For example, if we want to download the Duncan occupational prestige dataset from the carData package, we can use the code get_rdataset("Duncan","carData").data. Fit the linear model for an occupation’s prestige as a function of income variable based on the data in the Duncan dataset.

Hint: Use smf.ols based on the formula "prestige 1 + income".

from statsmodels.datasets import get_rdataset
Duncan = get_rdataset("Duncan", "carData").data
Duncan.head()

	type	income	education	prestige
rownames
accountant	prof	62	86	82
pilot	prof	72	76	83
architect	prof	75	92	90
author	prof	55	90	76
chemist	prof	64	86	90

import statsmodels.formula.api as smf
...

Ellipsis