# Linear Regression ## Libraries This is a little modified version of the code from the book `Introduction to Statistical Learning` by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani. The code is written in R and I have used the `tidyverse` package to make it more readable and easy to understand instead of using the base R functions. The first thing we need to do is install and then load the `tidyverse` set of R packages to provide us with lots of extra functionality. You only need to install this once: once it’s installed we can simply load it into the workspace using the `library(packagename)` function each time we open a new R session. Understanding data sets requires many hours/days or in some cases weeks.There are many commercially available software but open source community based software have now dominated and R is one of these. Here I also load the `MASS` package, which is a very large collection of data sets and functions. We also load the `ISLR2` package, which includes the data sets associated with the book `Introduction to Statistical Learning`. ```{r chunk1} #| warning: false library(MASS) library(ISLR2) library(tidyverse) ``` ## Simple Linear Regression Here's a tabular format representing the meta data for the `Boston` data set from the `ISLR2` library. This table provides the column name, a brief description, and the type of data in each column. | Column Name | Description | Data Type | |-------------|-------------|-----------| | `crim` | Per capita crime rate by town | Numeric | | `zn` | Proportion of residential land zoned for lots over 25,000 sq. ft. | Numeric | | `indus` | Proportion of non-retail business acres per town | Numeric | | `chas` | Charles River dummy variable (1 if tract bounds river; 0 otherwise) | Categorical | | `nox` | Nitrogen oxides concentration (parts per 10 million) | Numeric | | `rm` | Average number of rooms per dwelling | Numeric | | `age` | Proportion of owner-occupied units built prior to 1940 | Numeric | | `dis` | Weighted distances to five Boston employment centers | Numeric | | `rad` | Index of accessibility to radial highways | Categorical | | `tax` | Full-value property-tax rate per $10,000 | Numeric | | `ptratio` | Pupil-teacher ratio by town | Numeric | | `black` | 1000(Bk - 0.63)^2 where Bk is the proportion of Black people by town | Numeric | | `lstat` | Percentage of lower status of the population | Numeric | | `medv` | Median value of owner-occupied homes in $1000's | Numeric | This meta data outlines the features available in the `Boston` dataset that can be used to predict the `medv` variable. ```{r chunk2} head(Boston) ``` To find out more about the data set, we can type `?Boston`. We will start by using the `lm()` function to fit a simple linear regression model, with `medv` as the response and `lstat` as the predictor. The basic syntax is `lm(y ~ x, data)`, where `y` is the response, `x` is the predictor, and `data` is the data set in which these two variables are kept. ```{r chunk3, error=TRUE} lm_fit <- lm(medv ~ lstat) ``` The command causes an error because `R` does not know where to find the variables `medv` and `lstat`. The next line tells `R` that the variables are in `Boston`. If we attach `Boston`, the first line works fine because `R` now recognizes the variables. ```{r chunk4} lm_fit <- lm(medv ~ lstat, data = Boston) attach(Boston) # attach the Boston data set lm_fit <- lm(medv ~ lstat) # fit the model with data already attached ``` If we type `lm_fit`, some basic information about the model is output. For more detailed information, one may use `summary(lm_fit)` and I am using here `modelsummary` package to present summary in a more readable format. The `summary()` function outputs the coefficients of the model as well as their standard errors, $t$-statistics, and $p$-values. The `summary()` function also outputs the $R^2$ statistic and an analysis of variance (ANOVA) table, which breaks down the variance associated with the regression model and the residuals. This gives us $p$-values and standard errors for the coefficients, as well as the $R^2$ statistic and $F$-statistic for the model. ```{r chunk5} library(modelsummary) msummary(lm_fit) ``` One may extract the coefficients of the model using the `coef()` function. The `names()` function can be used to extract the names of the coefficients. ```{r chunk6} names(lm_fit) coef(lm_fit) ``` For confidence intervals, we can use the `confint()` function. By default, `confint()` provides 95 \% confidence intervals; however, this can be changed using the `level` argument. ```{r chunk7} confint(lm_fit) ``` To predict the median house value for a given percentage of lower status of the population, we can use the `predict()` function. The `predict()` function can be used to produce confidence intervals and prediction intervals for the prediction. ```{r chunk8} predict(lm_fit, data.frame(lstat = (c(5, 10, 15))), interval = "confidence") predict(lm_fit, data.frame(lstat = (c(5, 10, 15))), interval = "prediction") ``` ```{r} library(modelsummary) library(broom) # Fit the linear model lm_fit <- lm(medv ~ lstat, data = Boston) # Create predictions with confidence intervals predictions <- predict(lm_fit, newdata = data.frame(lstat = c(5, 10, 15)), interval = "confidence") # Convert predictions to a data frame pred_df <- as.data.frame(predictions) pred_df$lstat <- c(5, 10, 15) library(gt) pred_df |> gt() %>% tab_header(title = "Predicted Values with Confidence Intervals") ``` For instance, the 95 \% confidence interval associated with a `lstat` value of 10 is $(24.47, 25.63)$, and the 95 \% prediction interval is $(12.828, 37.28)$. As expected, the confidence and prediction intervals are centered around the same point (a predicted value of $25.05$ for `medv` when `lstat` equals 10), but the latter are substantially wider. We will now plot `medv` and `lstat` along with the least squares regression line using ggplot2. ```{r chunk9} ggplot(Boston, aes(x = lstat, y = medv)) + geom_point() + geom_smooth(method = "lm", se = FALSE) ``` There is some evidence for non-linearity in the relationship between `lstat` and `medv`. We will explore this issue later in this lab. ```{r chunk10} ggplot(Boston, aes(x = lstat, y = medv)) + geom_point() + geom_smooth(method = "lm", se = FALSE) + geom_abline(intercept = 34, slope = -1, color = "red", size = 3) + geom_point(color = "red") + geom_point(shape = 20) ``` Next we examine some diagnostic plots, several of which were discussed in Section 3.3.3. Four diagnostic plots are automatically produced by applying the `plot()` function directly to the output from `lm()`. In general, this command will produce one plot at a time, and hitting *Enter* will generate the next plot. However, it is often convenient to view all four plots together. We can achieve this by using the `par()` and `mfrow()` functions, which tell `R` to split the display screen into separate panels so that multiple plots can be viewed simultaneously. For example, `par(mfrow = c(2, 2))` divides the plotting region into a $2 \times 2$ grid of panels. ```{r chunk11} # Load necessary libraries library(broom) # Fit the linear model lm_fit <- lm(medv ~ lstat, data = Boston) # Extract augmented data for diagnostics augmented_data <- augment(lm_fit) ``` ```{r chunk12} # Residuals vs Fitted Plot ggplot(augmented_data, aes(.fitted, .resid)) + geom_point() + geom_smooth(method = 'loess', col = 'blue', se = FALSE) + labs(title = "Residuals vs Fitted", x = "Fitted values", y = "Residuals") + theme_minimal() ``` ```{r chunk13} # Normal Q-Q Plot ggplot(augmented_data, aes(sample = .std.resid)) + geom_qq() + geom_qq_line() + labs(title = "Normal Q-Q", x = "Theoretical Quantiles", y = "Standardized Residuals") + theme_minimal() ``` ```{r chunk14} # Scale-Location Plot (Spread of residuals) ggplot(augmented_data, aes(.fitted, sqrt(abs(.std.resid)))) + geom_point() + geom_smooth(method = 'loess', col = 'blue', se = FALSE) + labs(title = "Scale-Location", x = "Fitted values", y = "Sqrt(|Standardized Residuals|)") + theme_minimal() ``` ```{r chunk15} # Residuals vs Leverage Plot ggplot(augmented_data, aes(.hat, .std.resid)) + geom_point() + geom_smooth(method = 'loess', col = 'blue', se = FALSE) + geom_hline(yintercept = c(-3, 3), linetype = "dashed", color = "red") + labs(title = "Residuals vs Leverage", x = "Leverage", y = "Standardized Residuals") + theme_minimal() ``` ```{r chunk-infpts} # Highlight influential points # Replace .rownames with row_number if needed # Add a row number column to identify rows augmented_data <- augmented_data %>% mutate(row_id = row_number()) # Identify influential points using the row_id influential_points <- augmented_data %>% filter(.hat > 2 * mean(.hat)) %>% pull(row_id) # Print the identified influential points print(influential_points) print(influential_points) ``` On the basis of the residual plots, there is some evidence of non-linearity. Leverage statistics can be computed for any number of predictors using the `hatvalues()` function. ```{r hatvalues} leverage <- hatvalues(lm_fit) ``` The `which.max()` function identifies the index of the largest element of a vector. In this case, it tells us which observation has the largest leverage statistic. ```{r leverage} which.max(leverage) ``` ## Multiple Linear Regression Mutiple linear regression is an extension of simple linear regression to the case of two or more predictors. Each predictor has a regression coefficient (i.e. a slope), and there is one intercept. The model is given by $$ Y_i = \beta_0 + \beta_1 X_1i + \beta_2 X_2i + \ldots + \beta_p X_pi + \epsilon_i. $$ $$ Y_i = \beta_0 + \beta_1 X_1i + \beta_2 X_2i + \ldots + \beta_p X_pi + \epsilon_i. $$ The `summary()` function produces a detailed summary of the regression fit. I will use `modelsummary` package to produce a more readable output. ```{r chunk16} lm_fit <- lm(medv ~ lstat + age, data = Boston) library(modelsummary) modelsummary(lm_fit, estimate = "std.error", stars = TRUE) ``` The `Boston` data set contains 12 variables, and so it would be cumbersome to have to type all of these in order to perform a regression using all of the predictors. Instead, we can use the following short-hand: ```{r chunk17} lm_fit <- lm(medv ~ ., data = Boston) msummary(lm_fit, estimate = "std.error", stars = TRUE) ``` We can access the individual components of a summary object by name. The `names()` function can be used to obtain the names of the components. The `summary()` function returns a list with components such as `call`, `terms`, `residuals`, `coefficients`, `aliased`, `sigma`, `df`, `r.squared`, `adj.r.squared`, `fstatistic`, `cov.unscaled`, `na.action`. The `vif()` function, part of the `car` package, can be used to compute variance inflation factors. Most VIF's are low to moderate for this data. The `car` package is not part of the base `R` installation so it must be downloaded the first time you use it via the `install.packages()` function in `R`. ```{r chunk18} library(car) ``` What if we would like to perform a regression using all of the variables but one? For example, in the above regression output, `age` has a high $p$-value. So we may wish to run a regression excluding this predictor. The following syntax results in a regression using all predictors except `age`. ```{r chunk19} lm_fit1 <- lm(medv ~ . - age, data = Boston) msummary(lm_fit1) ``` Alternatively, the `update()` function can be used. ```{r chunk20} lm_fit1 <- update(lm_fit, ~ . - age) ``` ## Interaction Terms It is easy to include interaction terms in a linear model using the `lm()` function. The syntax `lstat:age` tells `R` to include an interaction term between `lstat` and `age`. The syntax `lstat * age` simultaneously includes `lstat`, `age`, and the interaction term `lstat`$\times$`age` as predictors; it is a shorthand for `lstat + age + lstat:age`. %We can also pass in transformed versions of the predictors. ```{r chunk21} msummary(lm(medv ~ lstat * age, data = Boston), estimate = "std.error", stars = TRUE) ``` ## Non-linear Transformations of the Predictors The `lm()` function can also accommodate non-linear transformations of the predictors. For instance, given a predictor $X$, we can create a predictor $X^2$ using `I(X^2)`. The function `I()` is needed since the `^` has a special meaning in a formula object; wrapping as we do allows the standard usage in `R`, which is to raise `X` to the power `2`. We now perform a regression of `medv` onto `lstat` and `lstat^2`. ```{r chunk22} lm_fit2 <- lm(medv ~ lstat + I(lstat^2)) msummary(lm_fit2, estimate = "std.error", stars = TRUE) ``` The near-zero $p$-value associated with the quadratic term suggests that it leads to an improved model. We use the `anova()` function to further quantify the extent to which the quadratic fit is superior to the linear fit. ```{r chunk23} summary(lm(medv ~ ., data = Boston)) anova(lm_fit, lm_fit2) ``` Here Model 1 represents the linear submodel containing only one predictor, `lstat` $$ medv_i = \beta_0 + \beta_1 \times lstat_i + \epsilon_i, $$ while Model 2 corresponds to the larger quadratic model that has two predictors, `lstat` and `lstat^2` $$ medv_i = \beta_0 + \beta_1 \times lstat_i + \beta_2 \times lstat_i^2 + \epsilon_i, $$ The `anova()` function performs a hypothesis test comparing the two models. The null hypothesis is that the two models fit the data equally well, and the alternative hypothesis is that the full model is superior. Here the $F$-statistic is $135$ and the associated $p$-value is virtually zero. This provides very clear evidence that the model containing the predictors `lstat` and `lstat^2` is far superior to the model that only contains the predictor `lstat`. This is not surprising, since earlier we saw evidence for non-linearity in the relationship between `medv` and `lstat`. Diagnostic plots can be used to further investigate the quality of the model fit using residuals. ```{r chunk24} par(mfrow = c(2, 2)) plot(lm_fit2) ``` ```{r diagnostic} par(mfrow = c(2, 2)) # Base R diagnostic plots plot(lm_fit2) # Enhanced diagnostic plots using the `car` package library(car) residualPlots(lm_fit2) qqPlot(lm_fit2) spreadLevelPlot(lm_fit2) influencePlot(lm_fit2) ``` ```{r} library(olsrr) ols_plot_resid_fit(lm_fit2) ``` ```{r} library(olsrr) ols_plot_resid_hist(lm_fit2) ``` ```{r} library(olsrr) ols_plot_resid_qq(lm_fit2) ``` ```{r} ols_test_correlation(lm_fit2) ``` then we see that when the `lstat^2` term is included in the model, there is little discernible pattern in the residuals. In order to create a cubic fit, we can include a predictor of the form `I(X^3)`. However, this approach can start to get cumbersome for higher-order polynomials. A better approach involves using the `poly()` function to create the polynomial within `lm()`. For example, the following command produces a fifth-order polynomial fit: ```{r chunk25} lm_fit5 <- lm(medv ~ poly(lstat, 5)) msummary(lm_fit5) ``` This suggests that including additional polynomial terms, up to fifth order, leads to an improvement in the model fit! However, further investigation of the data reveals that no polynomial terms beyond fifth order have significant $p$-values in a regression fit. By default, the `poly()` function orthogonalizes the predictors: this means that the features output by this function are not simply a sequence of powers of the argument. However, a linear model applied to the output of the `poly()` function will have the same fitted values as a linear model applied to the raw polynomials (although the coefficient estimates, standard errors, and p-values will differ). In order to obtain the raw polynomials from the `poly()` function, the argument `raw = TRUE` must be used. Of course, we are in no way restricted to using polynomial transformations of the predictors. Here we try a log transformation. ```{r chunk26} msummary(lm(medv ~ log(rm), data = Boston)) ``` ## Qualitative Predictors We will now examine the `Carseats` data, which is part of the `ISLR2` library. We will attempt to predict `Sales` (child car seat sales) in $400$ locations based on a number of predictors. Carseats have `r nrow(Carseats)` rows and `r ncol(Carseats)` columns. ```{r chunk27} head(Carseats) Carseats |> janitor::clean_names() |> head() ``` The `Carseats` data includes qualitative predictors such as `shelveloc`, an indicator of the quality of the shelving location---that is, the space within a store in which the car seat is displayed---at each location. The predictor `shelveloc` takes on three possible values: *Bad*, *Medium*, and *Good*. Given a qualitative variable such as `shelveloc`, `R` generates dummy variables automatically. Below we fit a multiple regression model that includes some interaction terms. ```{r chunk28} lm_fit <- lm(Sales ~ . + Income:Advertising + Price:Age, data = Carseats) msummary(lm_fit) ``` The `contrasts()` function returns the coding that `R` uses for the dummy variables. ```{r chunk29} attach(Carseats) contrasts(ShelveLoc) ``` Use `?contrasts` to learn about other contrasts, and how to set them. `R` has created a `ShelveLocGood` dummy variable that takes on a value of 1 if the shelving location is good, and 0 otherwise. It has also created a `ShelveLocMedium` dummy variable that equals 1 if the shelving location is medium, and 0 otherwise. A bad shelving location corresponds to a zero for each of the two dummy variables. The fact that the coefficient for `ShelveLocGood` in the regression output is positive indicates that a good shelving location is associated with high sales (relative to a bad location). And `ShelveLocMedium` has a smaller positive coefficient, indicating that a medium shelving location is associated with higher sales than a bad shelving location but lower sales than a good shelving location. ## Writing Functions R comes with numerous built-in functions, and many more are accessible through R libraries. However, there are times when a specific operation is needed for which no existing function is available. In such cases, we can create our own custom functions. Below, we illustrate this by creating a function called `LoadLibraries()` that loads the `tidyverse`, `broom`, and `car` packages. Attempting to call the function before its definition results in an error, as shown: Next, we define the function. The `+` symbols are printed by R automatically and should not be typed by the user. The `{` symbol indicates that a block of commands is about to follow. Pressing *Enter* after `{` causes R to display the `+` prompt, allowing for the input of multiple commands. The block is closed with `}` to signal the end of the function. ```{r} LoadLibraries <- function() { library(tidyverse) library(broom) library(car) library(MASS) library(ISLR2) print("The libraries have been loaded.") } ``` Typing `LoadLibraries` in the console after defining the function will display its content: ```{r} LoadLibraries ``` Calling the function will load the specified libraries and display the print message: ```{r} LoadLibraries() ```