Lesson 07 - OLS Regression in R

Parameter Estimates

Recall that in ordinary least squares regression, we assume that the relationship between our variables \(X\) and \(Y\) is determined by a hypothetical model of the form: \[Y = \beta_0 + \beta_1 X + e\]

We collect a sample of paired observations \((x_i, y_i)\) which we use to create a fitted model:

\[\hat Y = \hat \beta_0 + \hat \beta_1 X\] \[Y = \hat \beta_0 + \hat \beta_1 X + \hat e\]

The parameter estimates that result in \(SSE = \sum\limits_{i=1}^n \hat e_i^2\) being minimized are given by the following formulas:

\[SXX = \sum\limits_{i=1}^n \left(x_i - \bar x \right)^2\] \[SXY = \sum\limits_{i=1}^n \left(x_i - \bar x \right)\left(y_i - \bar y \right)\]

\[\hat\beta_1= \frac{SXY}{SXX} \hspace{20px} \mathrm{and} \hspace{20px}\hat\beta_0 = \bar y - \hat\beta_1 \bar x\]

We will now walk through an example in which we use R to calculate these parameter estimates. We will begin the example by loading a real-world dataset from a text file.

Pearson Dataset

In this notebook, we will study the relationship between the heights of fathers and their adult sons. For this study, we will use the Pearson data set, collected in England in 1903. This data set consists of 1078 observations, each with two features:

• fheight gives the height of the father, in inches.
• sheight gives the height of the son, in inches.

Loading and Exploring the Data

We will begin by loading the data and generating a summary of it.

In the code chunk below, the function read.table() is used to read the date file “father_son.txt”. The parameter sep is set to “\t”, which indicates that the values in the data file are separated by tabs. The parameter header is set to TRUE, which tells R that the first row of the data file contains the names of the variables whose values are stored in the columns. The data is stored in a data frame named myData. We use the head() function to view the first few rows of this data frame.

myData <- read.table("father_son.txt", sep="\t", header=TRUE)
head(myData, n=8)

We will use the nrow() function to confirm that there are 1078 observations in the data frame.

nrow(myData)

[1] 1078

The summary() function can be used to provide information regarding the distribution of each of the variables in the data frame.

summary(myData)

    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  

The distributions of the heights of the fathers appears to be similar to that of the sons. The sons do seem to be slightly taller on average.

We will further explore the distributions of these variables by creating histograms.

par(mfrow=c(1,2))
hist(myData$fheight, ylim=c(0,200), xlim=c(55,80), breaks=20, col = "cyan",
     xlab="Father Height in Inches", main="Histogram of\n Father Height")
hist(myData$sheight, ylim=c(0,200), xlim=c(55,80), breaks=20, col = "orange",
     xlab="Son Height in Inches", main="Histogram of \n Son Height")

par(mfrow=c(1,1))

The histograms provide evidence that the heights of both populations are approximately normal.

Relationship Between Father and Son Height

We begin our study of the relationship between the heights of fathers and sons by calculating the correlation between the two variables.

cor(myData$fheight, myData$sheight)

[1] 0.5013383

A correlation of 0.5013383 indicates a moderate positive linear relationship.

We will further explore the relationship between the two variables by creating a scatterplot.

plot(myData$sheight ~ myData$fheight, pch=19, col=rgb(0.1, 0.2, 1.0, 0.6),
     xlab="Father Height (inches)", ylab="Son Height (inches)",
     main="Heights of Fathers and Sons")

This scatterplot does show a clear linear relationship, although it is one with a lot of noise. Given any particular father height, the range of possible son heights covers several inches.

Calulating Parameter Estimates Using the Formulas

Let’s use R to calculate the parameter estimates for the least squares regression lines using the formulas mentioned at the beginning of this lesson.

x <- myData$fheight
y <- myData$sheight
xbar <- mean(x)
ybar <- mean(y)
SXX <- sum((x - xbar)^2)
SXY <- sum((x - xbar)*(y - ybar))
beta1 <- SXY / SXX
beta0 <- ybar - beta1 * xbar
par_est <- c(beta0, beta1)
names(par_est) <- c("beta0", "beta1")
par_est

    beta0     beta1 
33.886604  0.514093 

plot(myData$sheight ~ myData$fheight, pch=19, col=rgb(0.1, 0.2, 1.0, 0.6),
     xlab="Father Height (inches)", ylab="Son Height (inches)",
     main="Heights of Fathers and Sons")
abline(beta0, beta1, col="red", lwd=2)

Calulating Parameter Estimates Using lm() Function

Conveniently, R has a built-in function for calculating parameter estimates (as well as many other quantities that are useful when performing regression analysis). The function is denoted lm(), which stands for “linear model”. We will demonstrate its usage below.

# mod <- lm(myData$sheight ~ myData$fheight)
mod <- lm(sheight ~ fheight, myData)
summary(mod)

Call:
lm(formula = sheight ~ fheight, data = myData)

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

Note that the estimates for the slope and the intercept in our model are exactly the same as what we calculated above. We can access these objects directly from the mod variable:

mod$coefficients

(Intercept)     fheight 
  33.886604    0.514093 

plot(myData$sheight ~ myData$fheight, pch=19, col=rgb(0.1, 0.2, 1.0, 0.6),
     xlab="Father Height (inches)", ylab="Son Height (inches)",
     main="Heights of Fathers and Sons")
abline(mod$coefficients, col="red", lwd=2)