Simple Linear Regression in R

The workhorse of statistical analysis is the linear model, particularly regression. Invented by Francis Galton to study relationship between parents and children described as regressing to the mean, it has become one of the most widely used modeling techniques. In this tutorial, we will focus on simple linear regression.

SIMPLE LINEAR REGRESSION

The simple linear regression determines the relationship between two variables. That is we are using one variable to tell us what we can expect from the other variable. The outcome variable, what we are trying to predict, is called the response variable, and the input variable, what we are using to predict, is called the predictor variable. The general idea of the simple linear regression is to use the predictor to come up with some average value of the response. The relationship is defined as: y = a + bx + E, where a is the intercept; and b is the slope; E is the error term; x is the predictor variable; and y is the outcome/response variable. The equation describes is essentially a straight line that goes through the data where a is the y-intercept and b is the slope. First, install the required packages for our analysis.

# install usingR and ggplot2 packages; packages already installed; loading them using library()
library(UsingR)

Loading required package: MASS Loading required package: HistData Loading required package: Hmisc Loading required package: lattice Loading required package: survival Loading required package: Formula Loading required package: ggplot2

Above is the list of items contained in UsingR package collection. To illustrate the simple linear regression, we use the fathers’ and sons’ data. We will also require ggplot2 for plotting.

# Require ggplot2 and UsingR
require(UsingR)
require(ggplot2)

EXPLORATORY ANALYSIS

Let’s check the first 10 and the last 10 observations for the fathers’ and sons’ data.

# The first 10 observation of our dataset using the print(head(data, n = 10)) function
print(head(father.son, n = 10))

# The last 10 observations of our dataset using print(tail(data, n = 10)) function
print(tail(father.son, n = 10))

The father.son data is Pearson’s data set on heights of fathers and their sons. The data set contains 1078 measurements of fathers’ heights and sons’ heights. Its class is a data.frame with 1078 observations on the following 2 variables: fheight for fathers’ heights in inches and sheight for sons’ heights in inches. This data set was used by Pearson to investigate regression. The above information can be found by just typing ?father.son in the R console. Alternatively, using the following commands can get you accustomed with the data set.

# Checking the class of the data set using class(dataset) function
class(father.son)

[1] "data.frame"

A data.frame is a collection of variables which share many of the properties of matrices and of lists, used as the fundamental data structure by most R’s modeling software.

# Checking the structure of the data set using str(dataset) function
str(father.son)

'data.frame':   1078 obs. of  2 variables:
 $ fheight: num  65 63.3 65 65.8 61.1 ...
 $ sheight: num  59.8 63.2 63.3 62.8 64.3 ...

The str(father.son) command tells us that our data set contains 1078 observations and 2 variables, which are numeric variables. Additionally, we can do a brief sommary of our data set using the summary() function.

#Summary statistics of father.son data set
summary(father.son)

    fheight         sheight     
 Min.   :59.01   Min.   :58.51  
 1st Qu.:65.79   1st Qu.:66.93  
 Median :67.77   Median :68.62  
 Mean   :67.69   Mean   :68.68  
 3rd Qu.:69.60   3rd Qu.:70.47  
 Max.   :75.43   Max.   :78.36  

We can visualize the distribution of fathers’ heights by dividing the x axis into bins and counting the number of observations in each bin. Histogram, geom_histogram in the ggplot2 package, displays the count with bars, while geom_freqpoly displays the counts with lines. Frequency polygons are more suitable when you want to compare the distribution across the levels of a categorical variable.

# Histogram of father's height distribution
ggplot(data = father.son, mapping = aes(x = fheight)) +
  geom_histogram(bins = 30, fill = "seagreen") +
  ggtitle("Histogram of Father's Height") +
  theme(plot.title = element_text(hjust = 0.5)) 

# Histogram of father's height using geom_freqpoly displaying the counts with lines
ggplot(data = father.son, mapping = aes(x = fheight)) +
  geom_freqpoly(bins = 30) +
  ggtitle("Histogram of Height Counts With Lines") +
  theme(plot.title = element_text(hjust = 0.5))

We can also plot the histogram of the son’s height. We can apply the same technique as the one used in the previous graphs.

# Histogram of the son's height using the ggplot2 package
ggplot(data = father.son, mapping = aes(x = sheight)) +
  geom_histogram(bins = 30, fill = "seagreen") +
  ggtitle(expression(atop("Histogram of Son's Height", atop(italic("Heights by Count", ""))))) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.subtitle = element_text(hjust = 0.5))

We have added new elements, which are title and subtitle, by passing these arguments to ggtitle() function. Let’s apply the same technique to the line graph.

# Histogram of son's height using geom_freqpoly() function displaying counts by line.
ggplot(data = father.son, mapping = aes(x = sheight)) +
  geom_freqpoly(bins = 30) +
  ggtitle(expression(atop("Histogram of Son's Height", atop(italic("Height With Line Count", ""))))) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.subtitle = element_text(hjust = 0.5))

We can infer that both father’s height and son’s height seem to be normally distributed, according to the plotted histograms. We can extend our analysis by looking at how father’s height influences son’s height using simple linear regression.

SIMPLE LINEAR REGRESSION: USING FATHERS’HEIGHTS to DETERMINE SONS’ HEIGHTS

# Using fathers' heights to predict sons' heights with a geom_smooth () function
ggplot(data = father.son, mapping = aes(x = fheight, y = sheight)) + 
  geom_point() +
  geom_smooth(method = "lm") + 
  labs(x = "Fathers", y = "Sons") +
  ggtitle(expression(atop("Scatterplot of Fathers vs. Sons", atop(italic("With Fitted Regression Line", ""))))) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.subtitle = element_text(hjust = 0.5))

This plot is built using fathers’ heights to predict sons’ heights using simple linear regression. The fathers’ heights are the predictors and the sons’ heights are the responses. The blue line running through the points is the regression line and the grey band around it represents uncertainty in the fit.

To generate this nice plot, we pass the argument, (method = “lm”), to the geom_smooth() function. To make the results of the regression available to us, we need to calculate a regression using the function lm.

CALCULATE LINEAR REGRESSION USING THE lm()FUNCTION

# Calculate Linear regression using lm() function
(height.lm <- lm(sheight ~ fheight, data = father.son))

Call:
lm(formula = sheight ~ fheight, data = father.son)

Coefficients:
(Intercept)      fheight  
    33.8866       0.5141  

The above formula specifies to regress sheight, the response or outcome variable, on fheight, the predictor variable, using the father.son data, and adds the intercept term automatically. The results show coefficients for (intercept) and fheight which is the slope for fheight, the predictor variable. How can we then interpret these results? The interpretation of this is that, for every extra inch of father’s height, we expect an extra half inch increase in son’s height. The intercept doesn’t make sense in this case because it represents the height of a son whose father height is zero height, which is nonsensical. To get a full report of the regression model, just pass the height.lm model to the summary function.

# Complete regression results using summary() function
(summary(height.lm))

Call:
lm(formula = sheight ~ fheight, data = father.son)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8772 -1.5144 -0.0079  1.6285  8.9685 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 33.88660    1.83235   18.49   <2e-16 ***
fheight      0.51409    0.02705   19.01   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.437 on 1076 degrees of freedom
Multiple R-squared:  0.2513,    Adjusted R-squared:  0.2506 
F-statistic: 361.2 on 1 and 1076 DF,  p-value: < 2.2e-16

This prints out a lot of information about the model, including the standard error, the t-values and p-values for the coefficients, the degrees of freedom, residual summary statistics and the results of the F-test.

ANOVA AS AN ALTERNATIVE REGRESSION

An alternative to running an ANOVA test is to fit a regression with just one categorical variable and no intercept term. To run ANOVA, we will use the tips data in the reshape2 package on which we will fit a regression. To learn more about the reshape2 package, just type ?reshape2 in the R console.

# Pass tips and package reshape2 to the data() function
data(tips, package = "reshape2")

Let’s check the content of the dataset tips.

str(tips)

'data.frame':   244 obs. of  7 variables:
 $ total_bill: num  17 10.3 21 23.7 24.6 ...
 $ tip       : num  1.01 1.66 3.5 3.31 3.61 4.71 2 3.12 1.96 3.23 ...
 $ sex       : Factor w/ 2 levels "Female","Male": 1 2 2 2 1 2 2 2 2 2 ...
 $ smoker    : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
 $ day       : Factor w/ 4 levels "Fri","Sat","Sun",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ time      : Factor w/ 2 levels "Dinner","Lunch": 1 1 1 1 1 1 1 1 1 1 ...
 $ size      : int  2 3 3 2 4 4 2 4 2 2 ...

# Checking the first 10 observations
(print(head(tips, n = 10)))

# Checking the last 10 observations
(print(tail(tips, n = 10)))

We are going to fit our model with tip regresses against day minus the intercept. Tip is our outcome variable and day is our predictor variable. We will use aov() function to fit a model with a categorical variable without intercept.

# Pass arguments to aov() function for an ANOVA test
(tips.anova <- aov(tip ~ day - 1, data = tips))

Call:
   aov(formula = tip ~ day - 1, data = tips)

Terms:
                      day Residuals
Sum of Squares  2203.0066  455.6866
Deg. of Freedom         4       240

Residual standard error: 1.377931
Estimated effects are balanced

To see the full regression output, use summary() function.

# Pass the model tips.anova to summary()function
summary(tips.anova)

           Df Sum Sq Mean Sq F value Pr(>F)    
day         4 2203.0   550.8   290.1 <2e-16 ***
Residuals 240  455.7     1.9                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Writting - 1 in the formula indicates that the intercept should not be included in the model; the categorical variable day is automatically setup to have a coefficient for each level. You can apply the same technique to the linear model using the lm() function.

# Linear model using the lm() function
(tips.lm <- lm(tip ~ day - 1, data = tips))

Call:
lm(formula = tip ~ day - 1, data = tips)

Coefficients:
 dayFri   daySat   daySun  dayThur  
  2.735    2.993    3.255    2.771  

# The complete regression output
summary(tips.lm)

Call:
lm(formula = tip ~ day - 1, data = tips)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2451 -0.9931 -0.2347  0.5382  7.0069 

Coefficients:
        Estimate Std. Error t value Pr(>|t|)    
dayFri    2.7347     0.3161   8.651 7.46e-16 ***
daySat    2.9931     0.1477  20.261  < 2e-16 ***
daySun    3.2551     0.1581  20.594  < 2e-16 ***
dayThur   2.7715     0.1750  15.837  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.378 on 240 degrees of freedom
Multiple R-squared:  0.8286,    Adjusted R-squared:  0.8257 
F-statistic: 290.1 on 4 and 240 DF,  p-value: < 2.2e-16

From the output, we notice that the F-value of an ANOVA test is 290.1 and the F-statistic of the linear model is also 290.1. They are basically the same, as are degrees of freedom. This shows that the ANOVA and regression were derived along the same lines and can accomplish the same analysis.

VISUALIZING THE COEFFICIENTS

Visualizing the coefficients and the standard errors should show the same results as computing them using the ANOVA formula. The POINT ESTIMATES for the MEAN are identical and the CONFIDENCE INTERVALS are similar, the difference due to slightly different calculations as shown in the following plots.

1. Calculating The Means And Confidence Interval Manually

# Install the plyr package : install it andd then load it
library(plyr)
require(plyr)

Having downloaded the required package, you can create an object called tips.by.day using the ddply() function.

# Create an object called tips.by.day using ddply() function from the plyr package
(tips.by.day <- ddply(tips, "day", summarize,
                     tip.mean = mean(tip),
                     tip.sd = sd(tip),
                     Length = NROW(tip),
                     tfrac = qt(p = .90, df = Length - 1),
                     Lower = tip.mean - tfrac * tip.sd / sqrt(Length),
                     Upper = tip.mean + tfrac * tip.sd / sqrt(Length)))

Now that we have the output of tips.by.day, we can extract them from summary for tips.lm. It is important to point out that nrow and ncol return the number of rows and columns present in a data set. NCOL() and NROW() do the same treating a vector as 1-column matrix. The tidyverse package would be helpful here as well.

# Extract them from the summary for tips.lm
(tips.Info <- summary(tips.lm))

Call:
lm(formula = tip ~ day - 1, data = tips)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2451 -0.9931 -0.2347  0.5382  7.0069 

Coefficients:
        Estimate Std. Error t value Pr(>|t|)    
dayFri    2.7347     0.3161   8.651 7.46e-16 ***
daySat    2.9931     0.1477  20.261  < 2e-16 ***
daySun    3.2551     0.1581  20.594  < 2e-16 ***
dayThur   2.7715     0.1750  15.837  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.378 on 240 degrees of freedom
Multiple R-squared:  0.8286,    Adjusted R-squared:  0.8257 
F-statistic: 290.1 on 4 and 240 DF,  p-value: < 2.2e-16

# Extract point estimates for the mean and standard errors of tips per day
(tips.Coef <- as.data.frame(tips.Info$coefficients[, 1:2]))

tips.Coef <- within(tips.Coef, {
  Lower <- Estimate - qt(p = 0.90, df = tips.Info$df[2]) * `Std. Error`
  Upper <- Estimate + qt(p = 0.90, df = tips.Info$df[2]) * `Std. Error`
  day <- row.names(tips.Coef)
  })

PLOT TIP.MEAN AGAINST DAY

library(ggplot2)
# PLOTTING tip.mean against day with horizontal error bars
ggplot(data = tips.by.day, mapping = aes(x = tip.mean, y = day)) +
  geom_point() +
  geom_errorbarh(aes(xmin = Lower, xmax = Upper), height = .4) +
  ggtitle("Tips by Day Calculated Manually") +
  theme(plot.title = element_text(hjust = 0.5))

# Plotting tip estimates against day with horizontal error bars displayed
ggplot(data = tips.Coef, mapping = aes(x = Estimate, y = day)) +
  geom_point() +
  geom_errorbarh(aes(xmin = Lower, xmax = Upper), height = .4) +
  ggtitle("Tips by Day Calculated From Regression Model") +
  theme(plot.title = element_text(hjust = 0.5))

Regression coefficients and confidence intervals are taken from a regression model and calculated manually. The point estimates for the mean are identical and the confidence intervals are very similar, the difference due to slightly different calculations. The y-axis labels are also different because when dealing with factors lm tacks on the name of the variable to the level value.

A new function and a new feature were used here. We have introduced within(), which is similar to with() in that it lets us refer to columns in a data.frame by name, but different in that we can create new columns within that data.frame, hence the name. Second, one of the columns was named Std.Error with a space. To refer to a variable with spaces in its name, even as a column in a data.frame, we must enclose the name in back tips.