R tutorial on linear regression

Introduction

Linear regression models (also known as “Ordinary Least Squares” model) allow us to determine if changing the values on a variable is associated with the values of another variable. In other words, if I make a 1-unit change in \(X\), how much does \(Y\) change? In fact, linear regression is similar to the algebraic equation for a simple line \((Y=mx+b)\), where \(m\) is the slope, \(X\) is the parameter that is changing, and \(b\) is the \(Y\)-intercept). In statistics, we use linear regression models to test the association between two or more variables where the outcome is a continuous data type.

In many text books, the linear regression model is called the “Ordinary Least Squares” or OLS model because it minimizes the squared errors (e.g., distance from the best-fit line).

The linear regression model is pretty robust when the assumptions don’t hold. Regardless, it’s good practice to test these assumptions.

There are several conditions that need to be satisfied in order for us to use the results from a linear regression model. These include:

• Outcome variable is normally distributed (parametric)
• Observations are independent
• Residuals are not correlated with the \(X\) variables (homoscedasticity)
• Association between \(X\) and \(Y\) is linear

However, the linear regression model is pretty robust to violations of these assumptions; hence, its popularity. Moreover, interpretting the output of a linear regression is pretty straighforward and communicable to those who may not have statsistical training.

Simple linear regression

The structural form of a linear regression model:

\(Y_i = \beta_0 + \beta_1 X_{1i} + \epsilon_i\)

Typical notations of the linear regression include:

• \(Y_i\) denotes the outcome (or dependent) variable for subject \(i\)
• \(X_{1i}\) denotes the predictor of interest \((X_1)\) for subject \(i\)
• \(\beta_0\) denotes the \(Y\)-intercept when \(X\) is zero
• \(\beta_1\) denotes the slope or the change in \(Y\) with a 1-unit change in \(X\)
• \(\epsilon_i\) denotes the error or residual for subject \(i\)

Data

This notebook will use the Pima Indians Diabetes Dataset, which you can access from the mlbench library. The data was originally posted on the UC Irvine Machine Learning Repository, and can be found on Kaggle. The Pima Indians Diabetes Dataset orginated from the National Institute of Diabetes and Digestive and Kidney Diseases as part of a study to predict diabetes among adult females (21 years old and older) of Pima Indian heritage.

The following variables are included in the data

• pregnant: Number of times pregnant
  
• glucose: Plasma glucose concentration (glucose tolerance test)
  
• pressure: Diastolic blood pressure (mm Hg)
  
• triceps: Triceps skin fold thickness (mm)
  
• insulin: 2-Hour serum insulin (mu U/ml)
  
• mass: Body mass index
  
• pedigree: Diabetes pedigree function
• age: Age (years)
• diabetes: Class variable (test for diabetes)

Packages

You will need to install and load the following packages:

#### Load the libraries
library("ggplot2")
# install.packages("devtools")
library("devtools")

Loading required package: usethis
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

# install.packages("predict3d")
## Note: 02-11-2023: This has been removed from CRAN, so you need to use the following code to install predict3d:
## devtools::install_github("cardiomoon/predict3d") ## Make sure to have the devtools installed
library("predict3d")
# install.packages("psych")
library("psych")

Attaching package: ‘psych’

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha

# install.packages("magrittr")
# library("magrittr") ## allows for rounding using the %>%
library("dplyr")

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

# install.packages("gtsummary") ## Allows for publication ready tables
library("gtsummary")
# install.packages("DescTools")
library("DescTools") ## Needed for normality testing

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: ‘DescTools’

The following objects are masked from ‘package:psych’:

    AUC, ICC, SD

# install.packages("nortest")
library("nortest") ## Needed for normality testing
# install.packages("lmtest")
library("lmtest") ## Need for heteroskedasticity testing

Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

# install.packages("sandwich")
library("sandwich")  ## Needed for estimating Huber-White sandwich standard errors

Load Data

Load the data from the mlbench package using the data() function. You can use the knitr::kable function to generate a table.

#### Load Data
data("PimaIndiansDiabetes", package = "mlbench")

knitr::kable(
  head(PimaIndiansDiabetes), caption = "Table 1. First six rows of the Pima Indians Diabetes Dataset"
)

Registered S3 method overwritten by 'rmarkdown':
  method         from
  print.paged_df     

Table 1. First six rows of the Pima Indians Diabetes Dataset

pregnant	glucose	pressure	triceps	insulin	mass	pedigree	age	diabetes
6	148	72	35	0	33.6	0.627	50	pos
1	85	66	29	0	26.6	0.351	31	neg
8	183	64	0	0	23.3	0.672	32	pos
1	89	66	23	94	28.1	0.167	21	neg
0	137	40	35	168	43.1	2.288	33	pos
5	116	74	0	0	25.6	0.201	30	neg

Visualize the data

Once the data have been loaded, take a look at the summary statistics. You can do this using the describeBy function, which is part of the psych package.

knitr::kable(
  describeBy(PimaIndiansDiabetes) %>% round(2) 
)

Warning: no grouping variable requested

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
pregnant	1	768	3.85	3.37	3.00	3.46	2.97	0.00	17.00	17.00	0.90	0.14	0.12
glucose	2	768	120.89	31.97	117.00	119.38	29.65	0.00	199.00	199.00	0.17	0.62	1.15
pressure	3	768	69.11	19.36	72.00	71.36	11.86	0.00	122.00	122.00	-1.84	5.12	0.70
triceps	4	768	20.54	15.95	23.00	19.94	17.79	0.00	99.00	99.00	0.11	-0.53	0.58
insulin	5	768	79.80	115.24	30.50	56.75	45.22	0.00	846.00	846.00	2.26	7.13	4.16
mass	6	768	31.99	7.88	32.00	31.96	6.82	0.00	67.10	67.10	-0.43	3.24	0.28
pedigree	7	768	0.47	0.33	0.37	0.42	0.25	0.08	2.42	2.34	1.91	5.53	0.01
age	8	768	33.24	11.76	29.00	31.54	10.38	21.00	81.00	60.00	1.13	0.62	0.42
diabetes*	9	768	1.35	0.48	1.00	1.31	0.00	1.00	2.00	1.00	0.63	-1.60	0.02

There are 768 subjects with nine variables. We can see that the average number of pregnancies is 3.85 with a standard deviation (SD) of 3.37) We also see that the average age of the sample was 33.24 years (SD, 11.76), average BMI was 31.99 (SD, 7.8), and the average glucose level was 120.89 (SD, 31.97).

Motivating example: Evalaute the association between Age and Glucose level

Let’s suppose our main research question is to determine whether age was associated with glucose level. We will set the glucose level as our dependent variable and age as our independent variable (or predictor of interest). Then, we update our linear regression model’s structural form

\(E[Glucose_i|Age_i] = \beta_0 + \beta_1Age_i + \epsilon_i,\)
where \(E[Glucose_i|Age_i]\) denotes the expected Glucose level for subject \(i\) given the Age of subject \(i\).

We call this the “expected” value because we are predicting this using a model. All regression models take existing data and attempt to make predictions. However, if your assumptions are violated, then these predictions are erroneous. As George Box once wrote, “All models are wrong, but some are useful.”

We can also represent this relationship in a directed acyclic graph (DAG) diagram.

I used DAGitty to generate the DAG diagram.

# Install packages if not already installed
#install.packages(c("dagitty", "ggdag"))

# Load libraries
library(dagitty)
library(ggdag)

Attaching package: ‘ggdag’

The following object is masked from ‘package:stats’:

    filter

# Define the causal DAG
dag <- dagitty("dag {
  Age -> Glucose
}")

# Adjust DAG layout
coordinates(dag) <- list(
  x = c(Age = 1, Glucose = 2),
  y = c(Age = 1, Glucose = 1)
)

ggdag(dag, text = TRUE) +
  theme_minimal() +
  ggtitle("Causal Relationship Between Age and Glucose Level")

DAG diagram illustrating the causal relationship between Age and Glucose level.

Visualize the association between Age and Glucose level

Let’s look at how Age is related to Glucose level by plotting their relationship. We’ll do this using ggplot. As Age increases, the Glucose level also increases. There appears to be a positive relationship between Age and Glucose level.

### Plot the association between the subject's age and glucose level
ggplot(PimaIndiansDiabetes, aes(x = age, y = glucose)) +
  geom_point() +
  stat_smooth()

Constructing the linear regression model

Now, we can construct a linear regression model with Glucose level as the dependent variable and Age as the independent variable (or predictor of interest). We will use the lm() function with Glucose level as the \(Y\) variable and Age as the \(X\) variable. The formula for a regression model in R uses the ~ symbol. For example, if was want to regress Age on Glucose level, we use the notation Glucose ~ Age.

By using the lm() function, we can construct the linear regression model: lm(Glucose ~ Age, data = PimaIndiansDiabetes). We can create an object that will contain this linear model; I called this object linear_model1.

We generate the 95% confidence interval (CI) by using the confint() function.

\(\beta_1\) coefficient is in the linear regression output as Age, which is 0.71642.

Here is how we put all of this together in R:

### Linear regression model (Y = Glucose, X = Age)
linear_model1 <- lm(glucose ~ age, data = PimaIndiansDiabetes)
summary(linear_model1)

Call:
lm(formula = glucose ~ age, data = PimaIndiansDiabetes)

Residuals:
     Min       1Q   Median       3Q      Max 
-126.453  -20.849   -3.058   18.304   86.159 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 97.08016    3.34095   29.06  < 2e-16 ***
age          0.71642    0.09476    7.56 1.15e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30.86 on 766 degrees of freedom
Multiple R-squared:  0.06944,   Adjusted R-squared:  0.06822 
F-statistic: 57.16 on 1 and 766 DF,  p-value: 1.15e-13

The lm() function does not generate 95% CI, so you will need to use the confint() function.

### Generate the 95% CI
confint(linear_model1)

                 2.5 %      97.5 %
(Intercept) 90.5216601 103.6386585
age          0.5304001   0.9024361

Interpret the linear regression output

Make sure to have the library gtsummary loaded to create tables from the regression output.

We are interested in the coefficients. To make interpreting the output easier, we can create a table to visualize the critical elements.

#### Present the output in a table
model1 <- tbl_regression(linear_model1, intercept = TRUE)
as_gt(model1) |> 
              gt::tab_header("Table 2. Linear regression model output (Glucose ~ Age)") |> 
              gt::tab_options(table.align='left')

Table 2. Linear regression model output (Glucose ~ Age)
Characteristic	Beta	95% CI1	p-value
(Intercept)	97	91, 104	<0.001
age	0.72	0.53, 0.90	<0.001
1 CI = Confidence Interval

The (Intercept) denotes the \(Y\) intercept when \(X\) is equal to zero. In this case, it would be where Glucose level would be on the linear plot when Age is equal to zero, which is 97.08-units of Glucose.

When possible, present the coefficient’s point estimate and the 95% CI. For example, a 1-year increase in Age is associated with a 0.72-unit increase in Glucose level (95% CI: 0.53, 0.90).

The Age coefficient denotes the change in Glucose level for a one-unit increase in Age. In other words, a 1-year increase in Age is associated with a 0.72-unit increase in Glucose level. Since the 95% CI is between 0.53-units and 0.90-units of Glucose, it does not include zero, so this association is statistically significant. We can also look a the p-value of the Age coefficient to determine whether this is statistically significant (<0.0001). However, it is preferable to present the 95% CI when describing the association between \(X\) and \(Y\).

While p-values indicate statistical significance, confidence intervals are more informative because they show the effect size, direction, and precision of the estimate. That’s why scientific reports often emphasize 95% CIs over p-values alone.

The Adjusted R squared denotes the amount of variance in our outcome variable (glucose) that is explained by the linear regression model. In other words, the current linear regression model explains 6.8% of the data.

Visualize predicted model

We can plot the linear form of the model against the actual data using the ggPredict() function.

ggPredict(linear_model1, digits = 1, show.point = TRUE, se = TRUE, xpos = 0.5)

Here is a version without the scatterplot.

ggPredict(linear_model1, digits = 1, show.point = FALSE, se = TRUE, xpos = 0.5)

Adding a confounder

We generate a new variable called pregnancy_history that is a dichotomous variable (0 = no history of pregnancy and 1 = history of pregnancy).

Let’s add to our current model by including a confounder. We have a variable called Pregnancies but this provides the number of past pregnancies. We want to create a new variable that has a dichotomous outcome: History of Pregnancy or No history of pregnancy. To do this, we need to to subset the data and create rules. Anyone with 1 or more pregnancies will be coded as 1. Anyone who does not have a history of pregnancy will be coded as 0.

#### Generate groups based on pregnancy history (Group 0 = 0, Group 1 = 1 or more pregnancies)
PimaIndiansDiabetes <- PimaIndiansDiabetes |>
  mutate(pregnancy_history = as.factor(ifelse(pregnant == 0, 0, 1)))

knitr::kable(table(PimaIndiansDiabetes$pregnancy_history),
             col.names = c("Pregnancy History Group", "Count"),
             caption = "Frequency of Pregnancy History Groups")

Frequency of Pregnancy History Groups

Pregnancy History Group	Count
0	111
1	657

We see that there are 657 women who a history of pregnancy and 111 women with no history of pregnancy.

A DAG diagram illustrating the relationship between Pregnancy History as a confounder on the Age to Glucose relationship can be drawn.

dag <- dagitty("dag {
  Age -> Glucose
  Pregnancy -> Age
  Pregnancy -> Glucose
}")

# Define coordinates for layout
coordinates(dag) <- list(
  x = c(Age = 2, Pregnancy = 2, Glucose = 3),
  y = c(Age = 2, Pregnancy = 1, Glucose = 2)
)

# Convert DAG to a tidy format for ggdag
dag_tidy <- tidy_dagitty(dag) %>%
  mutate(node_fill = ifelse(name == "Pregnancy", "lightblue", "white"))  # Assign colors


# Plot DAG with custom node colors
ggdag(dag_tidy) +
  geom_dag_point(aes(fill = node_fill), shape = 21, size = 21, color="lightgrey") +  # Colored nodes
  geom_dag_text(color = "black") +  # Node labels
  scale_fill_identity() +  # Use manually assigned colors
  theme_minimal() +
  ggtitle("Pregnancy History as a Confounder in Age-Glucose Relationship")

DAG diagram illustrating the causal relationship between Age and Glucose level and Pregnancy History as a confounder.

In our linear regression model, we have the following structural form:

\(E[Glucose_i|Age_i] = \beta_0 + \beta_1Age_i + \epsilon_i,\)
where \(E[Glucose_i|Age_i]\) denotes the expected Glucose level for subject \(i\) given the Age of subject \(i\).

But we can include a confounder pregnancy_history:

\(E[Glucose_i|Age_i, PregnancyHistory_i] = \beta_0 + \beta_1Age_i + \beta_2PregnancyHistory_i + \epsilon_i,\)
where \(E[Glucose_i|Age_i, PregnancyHistory_i]\) denotes the expected Glucose level for subject \(i\) given the Age of subject \(i\) controlling for Pregnancy History of subject \(i\).

All you need to do to add another variable in the regression model is to use the + symbol. For example lm(glucose ~ age + pregnancy_history, data = PimaIndiansDiabetes).

Using the lm() function, we can add pregnancy_history to the linear regression model. Here is how we put all of this together in R:

### Linear regression model (Y = Glucose, X1 = Age, X2 = Pregnancy History)
linear_model2 <- lm(glucose ~ age + pregnancy_history, data = PimaIndiansDiabetes)
summary(linear_model2)

Call:
lm(formula = glucose ~ age + pregnancy_history, data = PimaIndiansDiabetes)

Residuals:
     Min       1Q   Median       3Q      Max 
-125.715  -20.546   -2.991   17.316   87.734 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        102.00752    3.95100  25.818  < 2e-16 ***
age                  0.76050    0.09638   7.891 1.04e-14 ***
pregnancy_history1  -7.47264    3.22137  -2.320   0.0206 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30.77 on 765 degrees of freedom
Multiple R-squared:  0.07594,   Adjusted R-squared:  0.07352 
F-statistic: 31.43 on 2 and 765 DF,  p-value: 7.592e-14

confint(linear_model2)

                         2.5 %      97.5 %
(Intercept)         94.2514334 109.7636012
age                  0.5712954   0.9497004
pregnancy_history1 -13.7964201  -1.1488593

We can present the model output into a table.

#### Present the output in a table
model2 <- tbl_regression(linear_model2, intercept = TRUE)
as_gt(model2) %>% 
              gt::tab_header("Table 3. Linear regression model output with confounder (Glucose ~ Age + Pregnancy History)") %>% 
              gt::tab_options(table.align='left')

Table 3. Linear regression model output with confounder (Glucose ~ Age + Pregnancy History)
Characteristic	Beta	95% CI1	p-value
(Intercept)	102	94, 110	<0.001
age	0.76	0.57, 0.95	<0.001
pregnancy_history
0	—	—
1	-7.5	-14, -1.1	0.021
1 CI = Confidence Interval

You can see that the Age coefficient is slightly different from our first model. It is 0.76 with a 95% CI of 0.57, 0.95. Compare this to the previous model’s result, which was 0.72; 95% CI: 0.53, 0.90.

#### Merge the two linear regression model's outputs
model1 <- tbl_regression(linear_model1, intercept = TRUE)
model2 <- tbl_regression(linear_model2, intercept = TRUE)
table1 <- tbl_merge(tbls = list(model1, model2),
          tab_spanner = c("**Model 1**", "**Model 2**"))
as_gt(table1) %>% 
              gt::tab_header("Table 4. Comparison between linear regression models [Model 1 (crude) v. Model 2 (adjusted)]") %>% 
              gt::tab_options(table.align='left')

Table 4. Comparison between linear regression models [Model 1 (crude) v. Model 2 (adjusted)]
Characteristic	Model 1			Model 2
	Beta	95% CI1	p-value	Beta	95% CI1	p-value
(Intercept)	97	91, 104	<0.001	102	94, 110	<0.001
age	0.72	0.53, 0.90	<0.001	0.76	0.57, 0.95	<0.001
pregnancy_history
0				—	—
1				-7.5	-14, -1.1	0.021
1 CI = Confidence Interval

Model 1 is considered the crude model or the unadjusted model. Model 2 is the adjusted model because it is adjusting based on the Pregnancy History confounder. Notice that the \(\beta_{1, unadjusted}\) is 0.72 which is lower than the \(\beta_{1, adjusted}\) result which is 0.76.

Additionally, the Adjusted R squared is higher in model 2 (7.35%) compared to model 1, which was 6.82%. This means that Model 2 does a better job of explaining the data than Model 1.

Let’s plot Model 2’s results.

ggPredict(linear_model2, digits = 1, show.point = TRUE, se = TRUE, xpos = 0.5)

You can derive the difference between the groups with a history of pregnancy and without a history of pregnancy by substracting the intercept coefficients 102 - 94.5, which is 7.5, the \(\beta_2\) estimate for pregnancy_history.

We can see that the group that had a history of pregnancy is lower than the group that did not have a history of pregnancy. This makes sense when you look at the pregnancy_history coefficient. It is -7.5, which means that a subject with a history of pregnancy is associated with a 7.5 decrease in Glucose level (95% CI: -13.80, -1.15) compared to a subject without a history of pregnancy controlling for age. Therefore, for all ranges of Age, the group with a history of pregnancy will have Glucose levels that are 7.5 units lower than a group without a history of pregnancy. You can visualize this on the plot; the linear lines do not cross and remain constant across all ranges of Age. But there is a positive correlation between Age and Glucose level.

Evaluate residual plots

It is good practice to look at the residuals of the regression model to make sure that the assumptions hold.

Recall the assumptions for the linear regression model:

• Association between \(X\) and \(Y\) is linear
  
• Residuals are not correlated with the \(X\) variables (homoscedasticity)
  
• Observations are independent
  
• Outcome variable is normally distributed (parametric)

Let’s use the crude linear regression model from our previous example:

\(E[Glucose_i|Age_i] = \beta_0 + \beta_1Age_i + \epsilon_i,\)
where \(E[Glucose_i|Age_i]\) denotes the expected Glucose level for subject \(i\) given the Age of subject \(i\).

linear_model1 <- lm(glucose ~ age, data = PimaIndiansDiabetes)
summary(linear_model1)

Call:
lm(formula = glucose ~ age, data = PimaIndiansDiabetes)

Residuals:
     Min       1Q   Median       3Q      Max 
-126.453  -20.849   -3.058   18.304   86.159 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 97.08016    3.34095   29.06  < 2e-16 ***
age          0.71642    0.09476    7.56 1.15e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30.86 on 766 degrees of freedom
Multiple R-squared:  0.06944,   Adjusted R-squared:  0.06822 
F-statistic: 57.16 on 1 and 766 DF,  p-value: 1.15e-13

confint(linear_model1)

                 2.5 %      97.5 %
(Intercept) 90.5216601 103.6386585
age          0.5304001   0.9024361

Check Linearity

To check linearity, we look for a straight-line relationship between Age and Glucose using a scatter plot with the regression line.

# Scatter plot with regression line
ggplot(PimaIndiansDiabetes, aes(x = age, y = glucose)) +
  geom_point() +
  geom_smooth(method = "lm", color = "blue") +
  theme_minimal() +
  ggtitle("Linearity: Scatter Plot with Regression Line")

What to look for: A roughly straight line suggests that the relationship between Age and Glucose is linear.

Check Homoscedasticity (Constant Variance)

Homoscedasticity v. Heteroscedasticity.

We can test to see if the residuals are uncorrelated to with the \(X\) variables (homoscedasticity).

Since the model generates predictions, we check to see if the residual are correlated with fitted or predicted values of Glucose. If there is an association, then we have heteroscedasticity, which is a violation of the linear regression model assumption.

To check homoscedasticity, we can plot the residuals and see if they are associated with increasing values of the expected or predicted Glucose level. If there is no association, we should expect to see a uniform distribution across all ranges of the expected values of Glucose. That is, If the spread of residuals is roughly constant across the fitted values, the assumption holds.

#### Plot the residuals against the predicted model (Is it homoscedastic?)
# Residuals vs Fitted plot
residuals <- resid(linear_model1)
fitted_values <- fitted(linear_model1)

ggplot(PimaIndiansDiabetes, aes(x = fitted_values, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red") +
  theme_minimal() +
  ggtitle("Homoscedasticity: Residuals vs Fitted")

What to look for: No pattern (e.g., no funnel shape). Residuals should be evenly spread around the horizontal line at zero.

Reviewing the residual plot along the fitted values, there doesn’t appear to be any evidence of heteroskedasticity. Upon visual inspection, we can see that the residuals are uniform across the expected Glucose level range (also called the “fitted” values). We can verify this visual inspection by performing the Breusch-Pagan test of heteroskedasticity. We will need to install and load the lmtest package and use the bptest() function.

bptest(linear_model1)

    studentized Breusch-Pagan test

data:  linear_model1
BP = 2.4585, df = 1, p-value = 0.1169

According to the BP-test results, the p-value is 0.1169, which means that we fail to reject the null that the variance of the residuals are constant. In other words, there are no associations between the residuals of the model and the predicted values generated from the model.

If, however, there was heteroskedasticity, we could address this by estimating robust standard errors. For linear regression models, the most common method is to use the Huber-White sandwich estimation method. To do this, we need to use the sandwich package and the coeftest() function. (Note: In practice, I default to the robust standard errors rather than let the lm() function estimate these for me.)

### Huber-White sandwich estimation
robust1 <- coeftest(linear_model1, vcov = vcovHC(linear_model1, type = "HC1"))
robust1

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 97.080159   3.265127  29.732 < 2.2e-16 ***
age          0.716418   0.095573   7.496  1.82e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

confint(robust1)

                 2.5 %      97.5 %
(Intercept) 90.6704993 103.4898193
age          0.5288023   0.9040339

Recall that the 95% CI is calculated using the standard error (SE).

Notice that the standard errors are slightly difference from the ones estimated in the previous model. In the previous model, the 95% CI was between 0.530 and 0.902. In the model with the robust standard errors, the 95% CI was between 0.529 and 0.904. The differences are trivial in this example, but could be important when the 95% CI is close to the null value.

Check Independence

Independence of residuals is difficult to assess with a plot in the same way as other assumptions, but in practice, this is usually verified based on the study design. For time series or clustered data, the Durbin-Watson test or plots of residuals against time can be used to check for autocorrelation.

# Durbin-Watson test for autocorrelation in residuals
dwtest(linear_model1)

    Durbin-Watson test

data:  linear_model1
DW = 1.8627, p-value = 0.02839
alternative hypothesis: true autocorrelation is greater than 0

What to look for:

• DW ≈ 2 → No autocorrelation (ideal).
  
• DW < 1.5 → Positive autocorrelation (common in time-series data).
  
• DW > 2.5 → Negative autocorrelation.
  
• p-value < 0.05 → Suggests significant autocorrelation (violating independence assumption).

# Extract residuals
residuals <- resid(linear_model1)

# Create lagged residuals (shift residuals by one step)
dw_data <- data.frame(
  residuals = residuals[-1],       # Current residuals
  lagged_residuals = residuals[-length(residuals)]  # Previous residuals
)

# Durbin-Watson plot
ggplot(dw_data, aes(x = lagged_residuals, y = residuals)) +
  geom_point(color = "blue", size = 3, alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", linetype = "dashed", se = FALSE) +
  theme_minimal() +
  ggtitle("Durbin-Watson Plot: Residuals vs Lagged Residuals") +
  xlab("Lagged Residuals") +
  ylab("Residuals")

What to look for:

• Random scatter: No autocorrelation (good).
• Upward or downward trend: Suggests positive or negative autocorrelation (bad).
• Strong clustering: Indicates possible model misspecification

Check Normality of Residuals

We can also evaluate if the residuals are normally distributed. We can generate a histogram and a Q-Q plot. The histogram has a slight left skew and the Q-Q plot has its tails deviate from the neutral line.

# Extract residuals from the model
residuals <- resid(linear_model1)

# Plot histogram of residuals
ggplot(data.frame(residuals), aes(x = residuals)) +
  geom_histogram(binwidth = 5, color = "black", fill = "lightblue", alpha = 0.7) +
  theme_minimal() +
  ggtitle("Histogram of Residuals") +
  xlab("Residuals") +
  ylab("Frequency")

# Q-Q plot for normality of residuals
ggplot(PimaIndiansDiabetes, aes(sample = residuals)) +
  geom_qq() +
  geom_qq_line() +
  theme_minimal() +
  ggtitle("Normality: Q-Q Plot of Residuals")

What to look for:
- Histogram should appear roughly bell-shaped centered arount zero.
- Points should follow a straight line if the residuals are normally distributed.

You only need to use one of these tests. They will generally give the same results.

We can also test for the normality of the residuals. Common tests of normality include the Shapiro-Wilk’s test, the Jarque Bera test, and the Kolmogorov-Smirnov (Lilliforms) test. I provided their codes below. Despite the differences in their output, the conclusions are all the same: the residuals are not normally distributed. Despite not being normally distributed, the linear regression model is pretty robust to violations of this assumption. You can make a concluding statement that there is an association between Age and Glucose level based on these findings.

#### Test normality using Shapiro-Wilk's test
shapiro.test(linear_model1$res)

    Shapiro-Wilk normality test

data:  linear_model1$res
W = 0.97467, p-value = 2.838e-10

 
#### Test normality using Jarque Bera test
JarqueBeraTest(linear_model1$res, robust = FALSE) ### Does not use robust method

    Jarque Bera Test

data:  linear_model1$res
X-squared = 25.565, df = 2, p-value = 2.81e-06

 
#### Test normality using the Kolmogorov-Smirnov test
lillie.test(linear_model1$res)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  linear_model1$res
D = 0.065397, p-value = 2.941e-08

Conclusions

Linear regression models are useful for understanding the relationship between a predictor variable and outcome varaible if the outcome variable is continuous. Additionally, you can add confounders into the regression model to control for their effects. Once you control for confounders, you should compare the results with the crude model to see how the relationship between the predictor of interest and outcome changes. Finally, after reviewing the results of the linear regression model, it is good practice to look at the residuals and verify that the assumptions of homoscedasticity and normality continue to hold.

References

• The gtsummary package is great at merging outcomes from regression models into publication quality tables. Daniel D. Sjoberg authored the gtsummary package with instructions on his website. You can also use external functions by converting the gtsummary table into an object using as_gt(); instructions can be found here.

• The sandwich package was used to estimate the Huber-White sandwich standard errors.

• Bruno Rodrigues has a wonderful article on dealing with heteroskedasticity.

• Czar Yobero wrote a great article on how to test for heteroskedasticity.