4.9. Summary of steps to build and investigate a linear model¶

Plot the data to assess model structure and degree of correlation between the \(\mathrm{x}\) and \(\mathrm{y}\) variable.

S, 2 lines

plot(x, y)           # plot the raw data
lines(lowess(x,y))   # superimpose non-parametric smoother to see correlation

Fit the model and examine the printed output.
S, 3 lines
```
model <- lm(y ~ x)   # fit the model: "y as described by variable x"
summary(model)
confint(model)
```
- Investigate the model’s standard error, how does it compare to the range of the \(\mathrm{y}\) variable?
- Calculate confidence intervals for the model parameters and interpret them.

Visualize the model’s predictions in the context of the model building data.

S, 3 lines

plot(x, y)
lines(lowess(x,y))        # show the smoother
abline(model, col="red")  # and show the least squares model

Plot a normal probability plot, or a q-q plot, of the residuals. Are they normally distributed? If not, investigate if a transformation of the \(\mathrm{y}\) variable might improve them. But also see the additional plots on checking for non-linearity and consider adding extra explanatory variables.
S, 2 lines
```
library(car)
qqPlot(resid(model))
```
Plot the residuals against the \(\mathrm{x}\)-values. We expect to see no particular structure. If you see trends in the data, it indicates that a transformation of the \(\mathrm{x}\) variable might be appropriate, or that there are unmodelled phenomena in the \(\mathrm{y}\) variable - we might need an additional \(\mathrm{x}\) variable.
S, 2 lines
```
plot(x, resid(model))
abline(h=0, col="red")
```
Plot the residuals in time (sequence) order. We expect to see no particular trends in the data. If there are patterns in the plot, assess whether autocorrelation is present in the \(\mathrm{y}\) variable (use the acf(y) function in R). If so, you might have to sub-sample the data, or resort to proper time-series analysis tools to fit your model.
S, 3 lines
```
plot(resid(model))
abline(h=0, col="red")
lines(lowess(resid(model), f=0.2))   # use a shorter smoothing span
```
Plot the residuals against the fitted-values. By definition of the least-squares model, the covariance between the residuals and the fitted values is zero. You can verify that \(e^T\hat{y} = \sum_i^n{e_i\hat{y}_i} = 0\). A fan-shape to the residuals indicates the residual variance is not constant over the range of data: you will have to use weighted least squares to counteract that. It is better to use studentized residuals, rather than the actual residuals, since the actual residuals can show non-constant variance even though the errors have constant variance.
S, 3 lines
```
plot(predict(model), rstudent(model))
lines(lowess(predict(model), rstudent(model)))
abline(h=0, col="red")
```

Plot the predictions of \(\mathrm{y}\) against the actual values of \(\mathrm{y}\). We expect the data to fall around a 45 degree line.

S, 3 lines

plot(y, predict(model))
lines(lowess(y, predict(model), f=0.5))     # a smoother
abline(a=0, b=1, col="red")                 # a 45 degree line

4.9.1. A worked example of the workflow in Python¶

The steps above describe a general workflow that applies to any data analysis project, not just least squares. We illustrate the early steps of that workflow on the blender efficiency data set, which records the result of designed experiments where four factors were varied to study blending efficiency: ParticleSize, MixerDiameter, MixerRotation, and BlendingTime.

The six-step workflow on these data:

Define the objective: understand which factors drive BlendingEfficiency.

Get the data:

Python, 5 lines

import pandas as pd

blender = pd.read_csv(
    "https://openmv.net/file/blender-efficiency.csv"
)

Explore: look at a few rows, the column types, and a numeric summary.
Python, 4 lines
```
blender.head()
blender.tail()
blender.describe()
blender.info()
```
Clean: in this case the data are pre-cleaned, so we move on.

Calculate: a correlation matrix and scatter-plot matrix highlight which factors are most related to the outcome variable.

Python, 14 lines

from pandas.plotting import scatter_matrix

# Numeric correlation matrix:
blender.corr()

# Visual version: scatter plot of every
# pair of variables, with a kde on the
# diagonal.
scatter_matrix(
    blender,
    alpha=0.2,
    figsize=(10, 8),
    diagonal="kde",
)

Filtering and grouping are part of the daily work of anyone working with data, and Pandas makes both very compact:

Python, 11 lines

# Boolean indexing returns only the rows
# where the condition is true:
blender[blender["ParticleSize"] == 2]
blender[blender["ParticleSize"] <= 5]
blender[blender["ParticleSize"] > 5]

# groupby applies the same calculation
# to each value of the grouping variable:
blender.groupby("ParticleSize").mean()
blender.groupby("ParticleSize").std()
blender.groupby("ParticleSize").max()

Communicate: create a separate plot per particle size, so each subgroup can be inspected on its own axes:
Python, 6 lines
```
for psize, subset in blender.groupby("ParticleSize"):
    subset.plot.scatter(
        x="BlendingTime",
        y="BlendingEfficiency",
        title=f"When particle size = {psize}",
    )
```
Once the patterns are clear, you can fit a least squares model using process_improve or statsmodels and continue with the residual diagnostics outlined in the steps above.

Summary of steps to build and investigate a linear model

Contents

4.9. Summary of steps to build and investigate a linear model¶

4.9.1. A worked example of the workflow in Python¶