Linear Regression with Python

Joaquín Amat RodrigoFebruary 2025

More about machine learning: Machine Learning with Python

Introduction¶

Linear regression is a statistical method that aims to model the relationship between a continuous variable and one or more independent variables by fitting a linear equation. It is called simple linear regression when there is only one independent variable and multiple linear regression when there is more than one. Depending on the context, the modeled variable is referred to as the dependent variable or response variable, while the independent variables are known as regressors, predictors, or features.

Throughout this document, the theoretical foundations of linear regression, key practical considerations, and examples of how to create these models in Python are progressively described.

Linear Regression Models in Python

Two of the most commonly used implementations of linear regression models in Python are: scikit-learn and statsmodels. Although both are highly optimized, Scikit-learn is mainly focused on prediction, which means it lacks functionalities to display many of the model characteristics necessary for inference. Statsmodels, on the other hand, is much more comprehensive in this regard.

Mathematical Definition¶

The linear regression model (Legendre, Gauss, Galton, and Pearson) assumes that, given a set of observations
$\{y_i, x_{i1},...,x_{ip}\}^{n}_{i=1}$, the mean $\mu$ of the response variable $y$ is linearly related to the regressors $x_1, ..., x_p$ according to the equation:

$$ \mu_y = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p} $$

The result of this equation is known as the population regression line, which captures the relationship between the predictors and the mean of the response variable.

Another common definition found in statistical literature is:

$$ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ... + \beta_p x_{ip} + \epsilon_i $$

In this case, the equation refers to the value of $y$ for a specific observation $i$. Since an individual observation will never be exactly equal to the mean, an error term $\epsilon$ is included.

In both cases, the interpretation of the model components is the same:

• $\beta_0$: The intercept, representing the average value of the response variable $y$ when all predictors are zero.

• $\beta_j$: The partial regression coefficient, representing the average effect on $y$ when predictor $x_j$ increases by one unit, keeping the other variables constant.

• $\epsilon$: The residual (error term), representing the difference between the observed and estimated values of $y$. It captures the influence of all factors affecting $y$ that are not included in the model.

In most cases, the population values $\beta_0$ and $\beta_j$ are unknown. Therefore, their estimates $\hat{\beta}_0$ and $\hat{\beta}_j$ are obtained from a sample. Fitting the model consists of estimating the regression coefficients from the available data in a way that maximizes the likelihood, meaning it identifies the model most likely to have generated the observed data.

The most commonly used method for this estimation is Ordinary Least Squares (OLS), which finds the best-fitting line (or plane, in the case of multiple regression) by minimizing the sum of squared vertical deviations between each training data point and the regression line.

Linear regression model and its errors: the gray line represents the regression line (the model) and the red segments show the error between the line and each observation.

Computationally, these calculations are more efficient when performed in matrix form:

$$ \mathbf{y} = \mathbf{X}^T \mathbf{\beta} + \epsilon $$

$$ \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, \quad \mathbf{X} = \begin{bmatrix} 1 & x_{11} & \dots & x_{1p} \\ 1 & x_{21} & \dots & x_{2p} \\ 1 & \dots & \dots & \dots \\ 1 & x_{n1} & \dots & x_{np} \end{bmatrix}, \quad \mathbf{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{bmatrix}, \quad \mathbf{\epsilon} = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{bmatrix} $$

$$ \hat{\beta} = \underset{\beta}{\arg\min} \left( (\mathbf{y} - \mathbf{X}^T \mathbf{\beta})^2 \right) $$

Once the coefficients are estimated, the estimates for each observation ($\hat{y}_i$) can be obtained as:

$$ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_{i1} + \hat{\beta}_2 x_{i2} + \dots + \hat{\beta}_p x_{ip} $$

Finally, the estimate of the model's variance ($\hat{\sigma}^2$) is obtained as:

$$ \hat{\sigma}^2 = \frac{\sum^n_{i=1} \hat{\epsilon}_i^2}{n-p} = \frac{\sum^n_{i=1} (y_i - \hat{y}_i)^2}{n-p} $$

where $n$ is the number of observations and $p$ is the number of predictors.

Model interpretation¶

The main elements to interpret in a linear regression model are the predictor coefficients:

• $\beta_{0}$ is the intercept, which corresponds to the expected value of the response variable $y$ when all predictors are zero.

• $\beta_{j}$, the partial regression coefficients of each predictor, indicate the expected average change in the response variable $y$ when the predictor variable $x_{j}$ increases by one unit, while keeping the other variables constant.

The magnitude of each partial regression coefficient depends on the units in which the corresponding predictor variable is measured, so its magnitude is not associated with the importance of each predictor.

To determine the impact of each variable in the model, standardized partial coefficients are used. These are obtained by standardizing (subtracting the mean and dividing by the standard deviation) the predictor variables before fitting the model. In this case, $\beta_{0}$ corresponds to the expected value of the response variable when all predictors are at their mean value, and $\beta_{j}$ represents the expected average change in the response variable when the predictor variable $x_{j}$ increases by one standard deviation, keeping the other variables constant.

While regression coefficients are usually the primary focus when interpreting a linear model, there are many other aspects to consider (overall model significance, predictor significance, normality condition, etc.). These aspects are often given little attention when the sole objective of the model is to make predictions. However, they become highly relevant when performing inference, that is, explaining the relationships between predictors and the response variable. Each of these aspects will be introduced throughout this document.

Meaning of "Linear"

The term "linear" in regression models refers to the fact that the parameters are incorporated into the equation in a linear manner, not necessarily that the relationship between each predictor and the response variable must follow a linear pattern.

The following equation shows a linear model in which the predictor $x_1$ is not linearly related to y:

$$y = \beta_0 + \beta_1x_1 + \beta_2log(x_1) + \epsilon$$

In contrast, the following equation does not represent a linear model:

$$y = \beta_0 + \beta_1x_1^{\beta_2} + \epsilon$$

However, some nonlinear relationships can be transformed into a linear form:

$$y = \beta_0x_1^{\beta_1}\epsilon $$

Taking the logarithm on both sides:

$$\log(y) = \log(\beta_0) + \beta_1 \log(x_1) + \log(\epsilon)$$

This transformation allows the model to be expressed in a linear form, making it suitable for linear regression techniques.

Goodness of Fit of the Model¶

Once the model has been fitted, it is necessary to verify its usefulness, as even though it represents the best-fitting line among all possible ones, it may still have a large error. The most commonly used metrics to assess the quality of fit are the Residual Standard Error (RSE) and the Coefficient of Determination ($R^2$).

Residual Standard Error (RSE)

The Residual Standard Error (RSE) measures the average deviation of any estimated point from the regression line. It has the same units as the response variable. A common way to determine whether the RSE is high is by dividing it by the mean value of the response variable, thereby obtaining a percentage deviation.

$$RSE = \sqrt{\frac{1}{n-p-1}RSS}$$

In simple linear regression models, since there is only one predictor, $(n-p-1) = (n-2)$.

Coefficient of Determination ($R^2$)

$R^2$ describes the proportion of variance in the response variable explained by the model relative to the total variance. Its value ranges between 0 and 1. Since it is dimensionless, it has the advantage of being easier to interpret compared to RSE.

$$R^{2}=\frac{\text{Total Sum of Squares - Residual Sum of Squares}}{\text{Total Sum of Squares}}=$$

$$1-\frac{\text{Residual Sum of Squares}}{\text{Total Sum of Squares}} = $$

$$1-\frac{\sum(\hat{y_i}-y_i)^2}{\sum(y_i-\overline{y})^2}$$

In simple linear regression models, the $R^2$ value corresponds to the square of the Pearson correlation coefficient (r) between $x$ and $y$. However, this is not the case in multiple regression.

In multiple linear regression models, the more predictors included in the model, the higher the $R^2$ value. This happens because each predictor, even if only slightly, helps explain part of the variability observed in $y$. For this reason, $R^2$ cannot be used to compare models with different numbers of predictors.

$R^2_{adjusted}$ introduces a penalty to the $R^2$ value for each predictor added to the model. The value of the penalty depends on the number of predictors used and the sample size, i.e., the degrees of freedom. The larger the sample size, the more predictors can be included in the model. Adjusted $R^2$ helps find the best model, the one that explains the variability in $y$ most effectively with the fewest predictors.

$$R^{2}_{adj} = 1 - \frac{SSE}{SST} \times \frac{n-1}{n-p-1} = R^{2} - (1-R^{2}) \frac{n-1}{n-p-1} = 1 - \frac{SSE/df_{e}}{SST/df_{t}}$$

where $SSE$ is the variability explained by the model (Sum of Squares Explained), $SST$ is the total variability of $y$ (Sum of Squares Total), $n$ is the sample size, and $p$ is the number of predictors included in the model.

Model Significance F-test¶

One of the first results to evaluate when fitting a linear model is the outcome of the $F$-test. This test answers the question of whether the model as a whole can predict the response variable better than expected by chance, or equivalently, if at least one of the predictors in the model contributes significantly. To perform this test, the sum of squared residuals of the model of interest is compared to that of the model without predictors, which consists only of the mean (also known as the total sum of squares, $TSS$).

$$F = \frac{(TSS - RSS)/(p-1)}{RSS/(n-p)}$$

Frequently, the null and alternative hypotheses for this test are described as:

• $H_0$: $\beta_1$ = ... = $\beta_{p-1} = 0$
• $H_a$: at least one $\beta_i \neq 0$

If the $F$ test is significant, it implies that the model is useful, but it does not necessarily mean it is the best model. It could happen that one of its predictors is unnecessary.

Significance of the Predictors¶

In most cases, although regression analysis is applied to a sample, the ultimate goal is to obtain a linear model that explains the relationship between the variables in the entire population. This means that the generated model is an estimate of the population relationship based on the relationship observed in the sample and, therefore, is subject to variations. For each of the coefficients in the linear regression equation ($\beta_j$), its significance (p-value) and confidence interval can be calculated. The most commonly used statistical test is the $T$-test$, although there are non-parametric alternatives.

The significance test for the coefficients ($\beta_j$) in the linear model considers the following hypotheses:

• $H_0$: the predictor $x_j$ does not contribute to the model ($\beta_j = 0$), in the presence of the other predictors. In the case of simple linear regression, it can also be interpreted as there being no linear relationship between the two variables, meaning the slope of the model is zero $\beta_j = 0$.

• $H_a$: the predictor $x_j$ does contribute to the model ($\beta_j \neq 0$), in the presence of the other predictors. In the case of simple linear regression, it can also be interpreted as there being a linear relationship between the two variables, meaning the slope of the model is not zero $\beta_j \neq 0$.

Calculation of the T-statistic and the p-value:

$$t = \frac{\hat{\beta}_j}{se(\hat{\beta}_j)}$$

where

$$SE(\hat{\beta}_j)^2 = \frac{\sigma^2}{\sum^n_{i=1}(x_{ji} -\overline{x})^2}$$

The error variance $\sigma^2$ is estimated from the Residual Standard Error (RSE), which can be understood as the average difference by which the response variable deviates from the true regression line. In the case of simple linear regression, RSE is equivalent to:

$$RSE = \sqrt{\frac{1}{n-2}RSS} = \sqrt{\frac{1}{n-2}\sum^n_{i=1}(y_i -\hat{y}_i)}$$

• Degrees of freedom (df) = number of observations - 2 = number of observations - number of predictors - 1

• p-value = P(|t| > calculated t value)

Confidence Intervals

$$\hat{\beta}_j \pm t_{df}^{\alpha/2} SE(\hat{\beta}_1)$$

The fewer the number of observations $n$, the less capacity there is to calculate the model's standard error. As a consequence, the accuracy of the estimated regression coefficients decreases. This is especially important in multiple regression.

In models generated with statsmodels, along with the regression coefficient values, the t-statistic value for each is returned, along with the p-values and corresponding confidence intervals. This allows to know, in addition to the estimates, whether they are significantly different from 0, i.e., if they contribute to the model.

For the calculations described above to be valid, it must be assumed that the residuals are independent and that they are normally distributed with a mean of 0 and variance $\sigma^2$. When the normality condition is not met, there is the possibility of using permutation tests to calculate significance (p-value) and bootstrapping to calculate confidence intervals.

Categorical Variables as Predictors¶

When a categorical variable is introduced as a predictor, one level is taken as the reference level (usually coded as 0) against which the other levels are compared. If the categorical predictor has more than two levels, so-called dummy variables or one-hot encoding are generated. These are variables created for each level of the categorical predictor that can take the value of 0 or 1. Each time the model is used to predict a value, only one dummy variable per predictor takes the value of 1 (the one corresponding to the level the predictor takes in that case), while the rest are considered 0. The partial regression coefficient $\beta_j$ for each dummy variable indicates the average percentage influence of that level on the dependent variable $y$ in comparison to the reference level of that predictor.

The concept of dummy variables is better understood with an example. Suppose a model in which the response variable weight is predicted based on height and nationality of the subject. The nationality variable is qualitative with 3 levels (American, European, and Asian). Even though the initial predictor is nationality, a new variable is created for each level, which can be either 0 or 1. Thus, the equation of the complete model is:

$$weight = \beta_0 + \beta_1height + \beta_2american + \beta_3european + \beta_4asian$$

If the subject is European, the dummy variables for American and Asian are considered 0, so the model for this case is:

$$weight = \beta_0 + \beta_1height + \beta_3european$$

Prediction¶

Once a valid model has been generated, it is possible to predict the value of the response variable $y$ for new values of the predictor variables $x$. It is important to note that predictions should a priori be restricted to the range of values within which the observations used to train the model fall. This is important because, although linear models can be extrapolated, it is only within this range that we are confident that the conditions for the model to be valid are met. The equation generated by the regression is used to calculate the predictions:

$$\hat{y}_i= \hat{\beta}_0 + \hat{\beta}_1 x_{i1} + \hat{\beta}_2 x_{i2} + ... + \hat{\beta}_p x_{ip}$$

Since the model is obtained from a sample, the estimates of the regression coefficients have an associated error, and so do the predicted values generated from them. There are two ways of measuring the uncertainty associated with a prediction:

• Confidence intervals: the interval for the average expected value of the response variable $y$ for a given value of $x$.

$$\hat{y} \pm t_{\alpha / 2, n - 2} \sqrt{MSE \left(\frac{1}{n} + \frac{(x_k - \bar{x})^2}{\sum(x_i - \bar{x}^2)} \right)}$$

• Prediction intervals: the interval for the expected value of the response variable $y$ for a given value of $x$.

$$\hat{y} \pm t_{\alpha / 2, n - 2} \sqrt{MSE \left(1 + \frac{1}{n} + \frac{(x_k - \bar{x})^2}{\sum (x_i - \bar{x})^2} \right)}$$

Although both seem similar, the difference lies in the fact that confidence intervals are applied to the average expected value of $y$ for a given value of $x$, whereas prediction intervals are not applied to the average. For this reason, the latter are always wider than the former.

A characteristic that stems from how the margin of error is calculated in both confidence and prediction intervals is that the interval widens as the values of $x$ approach the extremes of the observed range.

Why does this happen? By paying attention to the equation for the standard error of the confidence interval, the numerator contains the term $(x_k - \bar{x})^2$ (the same happens for the prediction interval).

$$\sqrt{MSE \left(\frac{1}{n} + \frac{(x_k - \bar{x})^2}{\sum^n_{i=1}(x_i - \bar{x}^2)} \right)}$$

This term corresponds to the squared difference between the value $x_k$ for which the prediction is being made and the mean $\bar{x}$ of the observed values of the predictor $x$. The further $x_k$ is from $\bar{x}$, the larger the numerator becomes, and therefore the larger the standard error. This is why the intervals widen as the values of $x$ move away from the mean of the observed values.

Model Validation¶

Once the best model that can be created with the available data is selected, its ability to predict new observations that were not used for training must be verified. This way, it can be checked if the model can generalize. A commonly used strategy is to randomly split the data into two groups, fit the model with the first group, and estimate the prediction accuracy using the second.

The appropriate size of the partitions largely depends on the amount of available data and the level of confidence needed in the error estimation. An 80%-20% split generally yields good results.

A more robust strategy is to use cross-validation, which consists of dividing the data into $k$ groups (folds) and using $k-1$ of them to train the model and the remaining one to estimate the prediction error. This process is repeated $k$ times, each time using a different fold as the test set. The final error is the average of the errors obtained in each iteration.

Conditions for Linear Regression¶

For a linear regression model using least squares, and the conclusions drawn from it, to be completely valid, the assumptions underlying its mathematical development must be verified. In practice, it is rare that all assumptions are met or can be proven to be met. However, this does not mean the model is not useful. What is important is to be aware of these assumptions and the impact they may have on the conclusions drawn from the model.

No Collinearity or Multicollinearity:

In multiple linear models, the predictors must be independent; there should be no collinearity between them. Collinearity occurs when a predictor is linearly related to one or more of the other predictors in the model. As a result of collinearity, the individual effect of each predictor on the response variable cannot be accurately identified, which leads to an increase in the variance of the estimated regression coefficients to the point where it becomes impossible to establish their statistical significance. Additionally, small changes in the data cause large changes in the coefficient estimates. While true collinearity only exists when the simple or multiple correlation coefficient between predictors is 1, which rarely happens in reality, it is common to encounter what is known as near-collinearity or imperfect multicollinearity.

There is no specific statistical method to determine the existence of collinearity or multicollinearity between the predictors of a regression model. However, many practical rules have been developed to assess the impact on the model. The recommended steps are:

• If the coefficient of determination $R^2$ is high, but none of the predictors are significant, there are indications of collinearity.

• Create a correlation matrix that calculates the linear relationship between each pair of predictors. It is important to note that even if no high correlation coefficients are obtained, multicollinearity cannot be ruled out. There may be a nearly perfect linear relationship between three or more variables, and the simple correlations between pairs of these variables may not exceed 0.5.

• Generate a simple linear regression model for each predictor against the others. If any of these models has a high coefficient of determination $R^2$, it may indicate potential collinearity.

• Tolerance (TOL) and Variance Inflation Factor (VIF). These are two parameters that essentially quantify the same thing (one is the inverse of the other). The VIF of each predictor is calculated using the following formula:

$$VIF_{\hat{\beta}_j} = \frac{1}{1 - R^{2}}$$

$$Tolerancia_{\hat{\beta}_j} = \frac{1}{VIF_{\hat{\beta}_j}}$$

Where $R^2$ is obtained from the regression of predictor $X_j$ on the other predictors. This is the most recommended approach, and the reference limits commonly used are:

• VIF = 1: complete absence of collinearity.

• 1 < VIF < 5: the regression may be affected by some collinearity.

• 5 < VIF < 10: the regression may be highly affected by some collinearity.

• The tolerance term is $\frac{1}{VIF}$, so the recommended limits are between 1 and 0.1.
  

If collinearity between predictors is found, there are two possible solutions. The first is to exclude one of the problematic predictors while attempting to retain the one that, in the researcher's judgment, is truly influencing the response variable. This measure generally has little impact on the model's predictive capacity, as, in the presence of collinearity, the information provided by one predictor is redundant in the presence of the other. The second option is to combine the collinear variables into a single predictor, though this carries the risk of losing its interpretation.

When trying to establish cause-and-effect relationships, collinearity can lead to very erroneous conclusions, making it seem that one variable is the cause when, in reality, another is influencing that predictor.

Linear Relationship Between Numeric Predictors and the Response Variable

Each numeric predictor must be linearly related to the response variable $y$, while keeping the other predictors constant. If this condition is not met, the predictor should not be included in the model. The most recommended way to check this is by plotting the residuals of the model against each of the predictors. If the relationship is linear, the residuals will be randomly distributed around zero. These analyses are only approximate, as it is not possible to definitively determine if the relationship is truly linear while keeping the other predictors constant.

Normal Distribution of the Response Variable

The response variable must follow a normal distribution. To verify this, histograms, normal quantile plots, or hypothesis tests for normality are commonly used.

Constant Variance of the Response Variable (Homoscedasticity)

The variance of the response variable must remain constant across the range of predictors. To verify this, it is common to plot the model's residuals against each predictor. If the variance is constant, the residuals will scatter randomly with consistent dispersion and no specific pattern. A conical distribution is a clear sign of a lack of homoscedasticity. Homoscedasticity tests, such as the Breusch-Pagan test, can also be used.

The reason for this condition is as follows: linear models assume that the response variable ($y$) follows a normal distribution, with a mean ($\mu$) that can be modeled as a function of other variables (predictors), and whose variance ($\sigma$) is calculated using a dispersion constant and a function $\nu(\mu)$. This means that the variance is not directly modeled as a function of the predictors, but indirectly through its relationship with the mean, and it is a unique value.

$$\hat{\sigma}^2 = \frac{\sum^n_{i=1} \hat{\epsilon}_i^2}{n-p} = \frac{\sum^n_{i=1} (y_i - \hat{\mu})^2}{n-p}$$

What impact can this have in practice? While the estimates of the mean obtained by these models are good, the uncertainties associated with them, and consequently the confidence and prediction intervals that can be calculated, may not be as accurate.

No Autocorrelation (Independence)

The values of each observation must be independent of the others. It is recommended to plot the residuals in the order of the time at which the observations were recorded. If any pattern is present, this may indicate autocorrelation. The Durbin-Watson hypothesis test can also be used to check for autocorrelation.

Outliers, High Leverage, or Influential Points

It is important to identify observations that are outliers or may be influencing the model. The easiest way to detect them is through the residuals (more details to follow).

Sample Size

This is not a condition in itself, but if there are not enough observations, predictors that are not truly influential might appear to be. A common recommendation is that the number of observations should be at least 10 to 20 times the number of predictors in the model.

Parsimony (Occam's razor)

This term refers to the idea that the best model is the one that can explain the observed variability in the response variable with the fewest predictors, and therefore, with fewer assumptions.

The vast majority of conditions are verified using the residuals, so the model is typically generated first, and then the conditions are validated. In fact, fitting a model should be seen as an iterative process where the model is adjusted, its residuals are evaluated, and improvements are made. This continues until an optimal model is reached.

Outliers¶

Regardless of whether the model is accepted, it is always advisable to identify any potential outliers, high leverage observations, or influential points, as they could significantly influence the model. The decision to remove such observations should be carefully considered, depending on the purpose of the model. If the goal is predictive, a model without outliers or highly influential observations is usually better at predicting most cases. However, it is very important to pay attention to these values because, if they are not measurement errors, they could represent the most interesting cases. The proper course of action when suspecting an outlier or influential point is to fit the regression model both including and excluding the value in question.

Example of Simple Linear Regression¶

Suppose a sports analyst wants to know if there is a relationship between the number of times a baseball team’s players hit and the number of runs they score. If such a relationship exists and a model can be established, this model could be used to predict the outcome of a game.

Libraries¶

Libraries used in this example are:

Data¶

Exploratory Data Analysis¶

The first step before generating a simple regression model is to visualize the data to intuitively assess if a relationship exists and quantify this relationship using a correlation coefficient.

The graph and the correlation test show a linear relationship, with considerable strength (r = 0.61) and significance (p-value = 0.000339). It makes sense to attempt generating a linear regression model with the goal of predicting the number of runs based on the team's number of hits.

Model Fitting¶

A model is fitted using runs as the response variable and hits as the predictor. As in any predictive study, it is not only important to fit the model but also to quantify its ability to predict new observations. To evaluate this, the data is split into two groups: a training set and a test set.

Scikit-learn¶

Once the model is trained, its predictive ability is evaluated using the test set.

Statsmodels¶

The Statsmodels linear regression implementation is more complete than that of Scikit-learn, as it not only fits the model but also allows the calculation of statistical tests and necessary analyses to verify that the conditions upon which this type of model is based are met. Statsmodels has two ways to train the model:

• Indicating the model formula and passing the training data as a dataframe that includes the response variable and predictors. This method is similar to the one used in R.

• Passing two matrices, one for the predictors and another for the response variable. This is the same as the one used by Scikit-learn, with the difference that the predictors matrix must include an additional column of 1s.

Confidence Intervals of the Coefficients¶

Predicctions¶

Once the model is trained, predictions can be obtained for new data. The Statsmodels models allow you to calculate predictions in two ways:

• .predict(): returns only the prediction values.

• .get_prediction().summary_frame(): returns the predictions along with the associated confidence intervals.

Graphical Representation of the Model¶

In addition to the least squares line, it is advisable to include the upper and lower bounds of the confidence interval. This helps identify the region where, according to the generated model and for a given confidence level, the average value of the response variable is expected to lie.

Test error¶

Interpretation¶

The column (coef) returns the estimated value for the two parameters of the linear model equation ($\hat{\beta}_0$ and $\hat{\beta}_1$), which correspond to the intercept (intercept or const) and the slope.

The standard errors, the t-statistic value, and the two-tailed p-value for each of the two parameters are also shown. This allows us to determine whether the predictors are significantly different from 0, meaning they have importance in the model. For the generated model, both the intercept and the slope are significant (p-values < 0.05).

The R-squared value indicates that the model can explain 27.1% of the observed variability in the response variable (runs). Additionally, the p-value obtained in the F-test (Prob (F-statistic) = 0.00906) suggests that there is evidence that the variance explained by the model is greater than what would be expected by chance (total variance).

The generated linear model follows the equation:

$$\text{runs = -2367.7028 + 0.5529 hits}$$

For each unit increase in the number of hits, the number of runs increases on average by 0.5529 units.

The test error of the model is 59.34. The predictions of the final model deviate on average by 59.34 units from the actual value.

Multiple Linear Regression Example¶

Suppose the sales department of a company wants to study the influence of advertising through different channels on the number of product sales. A dataset is available containing the revenue (in millions) obtained from sales in 200 regions, as well as the budget allocated (also in millions) for advertisements on radio, TV, and newspapers in each region.

Libraríes¶

The libraries used in this example are:

Data¶

Relationship Between Variables¶

The first step in establishing a multiple linear model is to study the relationship between variables. This information is critical for identifying the best predictors for the model and detecting collinearity among predictors. Additionally, it is recommended to visualize the distribution of each variable.

Model Fitting¶

A multiple linear model is fitted with the goal of predicting sales based on the investment in the three advertising channels.

The model with all variables introduced as predictors has a high $R^2$ (0.894), meaning it is able to explain 89.4% of the observed variability in sales. The p-value of the model is significant (1.01e-75), so we can accept that the model is better than expected by chance, and at least one of the partial regression coefficients is different from 0.

According to the p-value obtained for the partial regression coefficient of newspaper (0.723), this variable does not contribute significantly to the model in the presence of the other predictors. The model is retrained, but this time excluding the predictor newspaper.

Confidence Intervals of the Coefficients¶

Diagnostic of the reesiduals¶

Inspección visual¶

The residuals do not appear to be randomly distributed around zero, and the variability does not seem to remain consistent across the X-axis. This pattern suggests a lack of homoscedasticity and normal distribution.

Normality Test¶

The residuals are checked to see if they follow a normal distribution using two statistical tests: the Shapiro-Wilk test and D'Agostino's K-squared test. The latter is included in the statsmodels summary under the name Omnibus.

In both tests, the null hypothesis assumes that the data follows a normal distribution. Therefore, if the p-value is not lower than the selected significance level alpha, there is no evidence to reject the assumption that the data follows a normal distribution.

Both tests provide strong evidence to reject the hypothesis that the data follows a normal distribution (p-value << 0.01).

Predictions¶

Once the model is trained, predictions can be made for new data. Statsmodels models allow for the calculation of confidence intervals associated with each prediction.

Error on the test set¶

Interpretation¶

The multiple linear regression model:

$$\text{sells = 2.9004 + 0.0456tv + 0.1904radio}$$

is able to explain 89.4% of the observed variance in sales (R-squared: 0.894, Adj. R-squared: 0.892). The $F$ test is significant (p-value: 3.69e-77), providing clear evidence that the model explains the variance in sales better than expected by chance. The statistical tests for each variable confirm that tv and radio are related to the number of sales and contribute to the model.

The normality conditions are not satisfied, so the confidence intervals estimated for the coefficients and the predictions are not reliable.

The test error (rmse) is 1.696. The predictions of the final model deviate, on average, by 1.696 units from the actual value.

Interaction Between Predictors¶

The linear model from which the conclusions have been drawn assumes that the effect on sales due to an increase in the budget of one of the media channels is independent of the budget spent on the others. For example, the linear model considers that the average effect on sales from increasing the TV ad budget by one unit is always 0.0456, regardless of the amount spent on radio ads. However, this does not necessarily have to be the case—there may be interaction between predictors such that the effect of each predictor on the response variable depends, to some extent, on the value of the other predictor.

As previously defined, a linear model with two predictors follows the equation:

$$y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\epsilon$$

According to this definition, an increase of one unit in predictor $x_1$ results in an average increase of the response variable y by $\beta_{1}$. Modifications to predictor $x_2$ do not alter this, and the same applies to $x_2$ with respect to $x_1$. In order for the model to account for the interaction between both, a third predictor, called the interaction term, is introduced, which is the product of predictors $x_1$ and $x_2$.

$$y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{1}x_{2} + e$$

Reorganizing the terms gives:

$$y=\beta_{0}+(\beta_{1} + \beta_{3}x_2)x_{1}+\beta_{2}x_{2}+ e$$

Now, the effect of $x_{1}$ on $y$ is no longer constant but depends on the value of $x_{2}$.

In multiple linear regression models that incorporate interactions between predictors, the hierarchical principle must be taken into account. This principle states that if an interaction between predictors is included in the model, the individual predictors involved in the interaction must always be included, regardless of whether their p-value is significant or not.

In statsmodels, interaction between predictors can be introduced in two ways:

• By specifying it in the formula with the * term: sells ~ tv * radio

• By adding a new column with the multiplication of the predictors whose interaction is to be included in the model.

The results show clear evidence that the interaction tv radio is significant (p-value* very close to 0).

The model that incorporates the interaction has an Adjusted R-squared of 0.967, a value higher than the 0.892 of the model that only considered the individual effect of the predictors. However, is this difference enough to claim that the model with interaction is superior? One way to answer this question is by using the F-test (ANOVA).

Comparison of Models Using F-test (ANOVA)¶

Consider a model $M$ and another smaller model $m$, formed by a subset of the predictors contained in $M$. If the difference in fit is very small, according to the principle of parsimony, model $m$ is more appropriate. It is possible to test if the difference in fit is significant by comparing the residuals. Specifically, the statistic used is:

$$\frac{RSS_m - RSS_M}{RSS_M}$$

To avoid the model size influencing the test, the residual sum of squares RSS of each model is divided by its degrees of freedom. The resulting statistic follows an $F$ distribution:

$$\frac{(RSS_m - RSS_M)/(df_m - df_M)}{RSS_M/(df_M)} \sim F_{df_m-df_M, \ df_M}$$

where $df$ are the degrees of freedom of the model, which correspond to the number of observations minus the number of predictors.

In the previous sections, a model that included interactions and one that did not were trained. A hypothesis test is now performed to evaluate the null hypothesis that both models fit the data equally well. In other words, to determine whether one model is better than the other in predicting the response variable. This can be done using the anova_lm() function from statsmodels.stats.anova.

The ANOVA test provides clear evidence (p-value virtually zero) that the model including the interaction between predictors is able to model the response variable better.

This method of comparing models is well established in the statistical community. In the machine learning community, it is much more common to compare models using validation techniques such as cross-validation.

Miscellaneous Notes¶

In this section, I gather comments, definitions, and clarifications I have found in different sources and that I have considered better to keep as supplementary information.

Origin of the Least Squares Method and Regression¶

Linear Models with R, by Julian J. Faraway

The least squares method was published in 1805 by Adrien Marie Legendre. The term regression comes from a publication by Francis Galton in 1885 called Regression to Mediocrity. In it, Galton used the least squares method to demonstrate that children of tall parents tend to be tall as well, but not as much as their parents, and that children of shorter parents tend to be short, but not as much as their parents.

Estimation of the Variance of a Linear Model by Least Squares¶

Linear Models with R, by Julian J. Faraway

The estimation of the variance ($\sigma^2$) of a linear model obtained by least squares is:

$$\sigma^2 = \frac{RSS}{n - p}$$

The term Residual Standard Error (RSE) is the square root of $\sigma^2$:

$$RSE = \sqrt{\frac{RSS}{n - p}}$$

Advantages of the Least Squares Method for Estimating the Coefficients of a Linear Model¶

Linear Models with R, by Julian J. Faraway

Although there are alternatives to the least squares method for estimating the correlation coefficients $\hat{\beta}_j$ of a linear model, it offers several advantages:

• It has a geometric interpretation, which makes it easier to understand.

• If the errors are independent and normally distributed, its solution is equivalent to the maximum likelihood estimation.

• The $\hat{\beta}_j$ are unbiased estimations.

Identifiability or Collinearity¶

Linear Models with R, by Julian J. Faraway

The least squares method has a unique solution only if the matrix formed by the predictors has full rank, meaning all its columns (predictors) are linearly independent. In practice, this identifiability condition is often violated. The following are some scenarios where this occurs:

• When one of the predictors introduced in the model is a linear transformation or combination of other predictors already in the model. For example, if the variable weight is introduced in both kilograms and pounds, or if predictors like the number of years of basic education, number of years of university education, and total years of education are included. This problem can be avoided by studying the nature of the available variables and their relationships.

• Model oversaturation, when there are more predictors than observations.

Caution When Evaluating Residual Normality Through Hypothesis Testing¶

Linear Models with R, by Julian J. Faraway

For the residuals of a linear regression model to be normally distributed is a necessary condition for the significance (p-value) and the confidence intervals associated with the predictors (calculated from theoretical models) to be accurate. This condition is often evaluated with hypothesis tests such as the Shapiro-Wilk test. When this happens, it is important to understand the relationship between p-value and sample size. With more residuals, the test has greater power, so small deviations from normality become significant. Additionally, the central limit theorem suggests that, as the sample size increases, the results become more robust to deviations from normality. Given these properties, it is usually more advisable to assess the normality of the residuals graphically, using a theoretical quantile-quantile (qq-plot) plot.

Transformation of Variables¶

Linear Models with R, by Julian J. Faraway and Wikipedia

Transforming the response variable or the predictors can be a way to improve the fit of a model or correct violations of certain regression assumptions.

When the transformation is applied to the response variable, it is important to remember that the predictions generated by the model are on the transformed scale. Once the prediction is made, it can be transformed back to the original scale. The same applies to confidence intervals: first, the upper and lower limits are calculated, and then they are transformed back. In contrast, the regression coefficients obtained for the predictors must be interpreted in the context of the transformation. It is not possible to apply the inverse of the transformation and interpret them on the original scale.

Box-Cox Transformation

The exponential transformation of the response variable ($y^{\lambda}$) is often a useful solution to improve models. However, identifying the optimal exponent value is not always straightforward. The Box-Cox method allows determining if there is a value to which the response variable can be raised to improve the model fit, and if so, identifying which value it is. According to this method, a transformation $y \rightarrow g_{\lambda}(y)$ is sought such that:

$$ g_{\lambda}(y) = \begin{cases} \frac{y^{\lambda}-1}{\lambda} & \text{if } \lambda \neq 0 \\ \log(y) & \text{if } \lambda = 0 \end{cases} $$

When the goal of the model is predictive, $\frac{y^{\lambda}-1}{y}$ can be simplified to $y^{\lambda}$.

The optimal value of $\lambda$ is identified by maximum likelihood: a range of $\lambda$ values is selected, and for each value, the transformation is applied, and the residual sum of squares is calculated. The value that results in the smallest residual sum is the best option.

The transformation of the response variable can make the model more difficult to interpret, so it is not advisable unless truly necessary. To ensure it is, besides calculating the optimal $\lambda$ value, its 95% confidence interval can be estimated. If this interval includes the value $\lambda=1$, there is insufficient evidence to support the transformation. When the transformation is necessary, it is recommended to choose the value within the confidence interval that facilitates model interpretation. For example, if the interval is [-0.7, 0.7], the value 0.5 is the most suitable, as it corresponds to $\sqrt(y)$.

The Box-Cox transformation is only applicable when the response variable always takes positive values. If any $y_i < 0$, a constant can be added to $y$ to make it positive, but this is a somewhat inelegant solution. Another feature of this transformation is that it is highly influenced by outliers, causing the absolute values of $\lambda$ to be suspiciously high.

Logarithmic Transformation

There are multiple alternatives to the Box-Cox exponential transformation. One of them consists of adding a constant to the response variable and then applying the logarithm: $g_{\alpha} = \log(y + \alpha)$.

Interpretation of Regression Coefficients in a Linear Model¶

Linear Models with R, by Julian J. Faraway

When using a linear regression model to extract information about the relationship between predictors and the response variable, one must be cautious and avoid drawing conclusions that go beyond what the model tells us.

The gala dataset from the faraway package contains information about the islands that make up the Galápagos archipelago. Among the stored variables are: the number of species inhabiting the island (Species), the island's area (Area), elevation (Elevation), distance to the nearest island (Nearest), distance to the Santa Cruz island (Scruz), and the area of the adjacent island (Adjacent). The goal is to use a linear regression model to study the relationship between the predictor Elevation and the number of species on the island.

If the regression coefficient is interpreted as saying that for each unit increase in Elevation, the number of species increases on average by 0.3195 units, that would be an incorrect conclusion. In fact, if a simpler model is fitted using the same data, the coefficient value changes.

This happens because the regression coefficient of a given predictor must always be interpreted in the context of the other predictors included in the model, as they influence it. The correct interpretation is: an increase of one unit in predictor $x_j$ leads to an average change of $\hat{\beta}_j$ units in the response variable, while keeping the other predictors constant.

Frequently, when the scale on which each predictor is measured differs greatly (order of magnitude), the predictors are standardized before fitting the model. In this case, the interpretation of the regression coefficients changes to: an increase of one standard unit in predictor $x_j$ leads to an average change of $\hat{\beta}_j$ standard units in the response variable, while keeping the other predictors constant.

In observational studies, it is rarely possible to keep the predictors constant at the will of the researcher, which makes it difficult to establish causal relationships.

Permutation Test (Monte Carlo Simulation) for Model and Predictor Significance¶

Linear Models with R, by Julian J. Faraway

The hypothesis test based on the F-test relies on the assumption that the residuals are normally distributed. While the central limit theorem shows that, even when the distribution deviates from normality, inference results are good if the sample size is large, there is no exact size beyond which validity is guaranteed. Permutation tests offer an alternative that does not require assuming normality of the residuals.

Suppose we have a model with two predictors, Species ~ Nearest + Scruz, for which the F statistic is calculated. The associated p-value indicates how likely it is to obtain an F value as extreme or more extreme, assuming no relationship between the response variable and the predictors. When the condition $\epsilon \sim N(0, \sigma^2)$ holds, this probability can be derived from the F-distribution.

When the normality condition is not met, these values are unreliable. A better approach is to resort to a permutation test. For this, the null hypothesis of "no association between the response variable and all predictors" is simulated by randomly swapping the response variable between observations.

If this process is repeated sufficiently many times (1000 - 10000), adjusting the model each time and storing the F statistic, the percentage of cases in which F is equal to or greater than the value obtained in the initial model can be calculated. This percentage is the empirical approximation of the p-value.

The p-value obtained through permutation is very similar to the one obtained by the theoretical method based on normality (0.555).

The permutation method has the advantage of being just as good as the theoretical method (if enough permutations are performed) when the normality condition is met, and it can also be applied when it is not. The downside is that it requires much more computational effort.

The same permutation approach can be used to determine if a particular predictor contributes significantly to the model, or in other words, if it is associated with the response variable. In this case, only the value of the predictor of interest is permuted, thus simulating the null hypothesis that there is no association with the response variable. In each permutation, the model is fitted, and the t-statistic for the predictor is stored. After performing enough permutations, the percentage of cases where the absolute value of t is equal to or greater than the absolute value observed in the initial model is calculated.

Suppose a model with two predictors, Species ~ Nearest + Scruz, for which we want to calculate the significance of the predictor Scruz.

The p-value obtained through permutations is very similar to the one obtained through the theoretical method based on normality (0.2833).

Bootstrapping to Calculate Confidence Intervals¶

Linear Models with R, by Julian J. Faraway

The method used by statsmodels and scikitlearn to calculate confidence intervals for the regression coefficients of a linear model is based on the t-student distribution and follows the structure:

$$\hat{\beta}_j \pm t_{n-p}^{\alpha/2} se(\hat{\beta})$$

For the intervals based on this theoretical distribution to be valid, it is necessary for the model's residuals to follow a normal distribution. If this is not the case, the bootstrapping method offers an alternative that does not require this condition.

Below is an example of how to apply this method to estimate the 95% CIs for the model Species ~ Area + Elevation + Nearest + Scruz + Adjacent.

In the book "Linear Models with R" by Julian J. Faraway, resampling is done on the residuals of the model. In my opinion, it is more intuitive to do the resampling on the pairs (response, predictors). The confidence interval estimates obtained are not exactly the same, so I cannot guarantee the validity of this approach.

Robust Regression¶

The ordinary least squares (OLS) regression method described in the previous sections assumes that the residuals/errors are independent, with constant variance, and normally distributed. In practice, these conditions are often not met, making it necessary to have alternative fitting methods. Some of these methods include Robust Regression, Generalized Least Squares, and Weighted Least Squares.

Robust Regression

When errors do not follow a normal distribution, the results obtained by least squares are affected, with the impact being greater the longer the tails are. A simple solution is to eliminate the outliers that form these tails. However, if it is confirmed that these are not reading errors, the model should include them as they are part of the phenomenon being studied. Robust regression reduces the influence of outliers on the model fitting. The two most commonly used types of robust regression are: M-Estimation and Least Trimmed Squares.

• M-Estimation: This has different variants depending on the function used to reduce the weight of extreme observations. The Huber method decreases the influence of outliers, as long as they are few and not excessively large.

• Least Trimmed Squares: This type of robust regression fits the model by using least squares, but only employs the q smallest residuals, completely ignoring the rest of the observations. Therefore, it is more independent of outliers than M-estimation, although it depends on the specified value of q.

In most cases, the results obtained using robust regression are similar to those obtained by using least squares and analyzing the outlying or influential observations. When this latter analysis cannot be done, either because the researcher doesn't know or because it is an automated process, robust regression is safer.

Given the speed at which these adjustment methods can be applied, robust regression can be used as a confirmation method. First, the model is fitted using least squares, and if no significant differences are found when compared to robust regression, it indicates that the model does not contain influential observations.

Generalized Least Squares

This is an alternative to least squares regression when there is correlation between the residuals. This frequently occurs in time series data.

Weighted Least Squares

This is an alternative to least squares regression that is not affected by the lack of constant variance in the residuals (heteroscedasticity). Some scenarios where its use is recommended are:

• When the residuals are proportional to one of the predictors, i.e., as the value of the predictor increases, so do the residuals.

• When the variable $Y_i$ is obtained from $n_i$ observations (for example, the average of $n$ repetitions), it can occur that the variance is proportional to the group size, i.e., $Var(y_i) = \sigma^2/n_i$.

• When the response variable comes from different sources, for example, different instruments, each with a different precision, a weight is assigned to each source $i$ such that $weight_i = 1/sd(y_i)$.

Session Information¶

Bibliography¶

OpenIntro Statistics: Third Edition, David M Diez, Christopher D Barr, Mine Çetinkaya-Rundel

An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)

Linear Models with R, Julian J.Faraway

Points of Significance: Simple linear regression by Naomi Altman & Martin Krzywinski

Points of Significance: Multiple linear regression Martin Krzywinski & Naomi Altman

Métodos estadísticos en ingenieria Rafael Romero Villafranca, Luisa Rosa Zúnica Ramajo

Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.

Citation¶

How to cite this document

If you use this document or any part of it, please acknowledge the source, thank you!

Linear Regression with Python by Joaquín Amat Rodrigo, , available under Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) at https://cienciadedata.net/documentos/py10-linear-regression-python.html

Did you like the article? Your support is important

Your contribution will help me to continue generating free educational content. Many thanks! 😊

This work by Joaquín Amat Rodrigo and Javier Escobar Ortiz is licensed under a Attribution-NonCommercial-ShareAlike 4.0 International.

Allowed:

• Share: copy and redistribute the material in any medium or format.

• Adapt: remix, transform, and build upon the material.

Under the following terms:

• Attribution: You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.

• NonCommercial: You may not use the material for commercial purposes.

• ShareAlike: If you remix, transform, or build upon the material, you must distribute your contributions under the same license as the original.