library(ggplot2)
library(dplyr)
library(tidyverse)
library(knitr)
library(ggthemes)
library(ISLR2)
library(gglm)
theme_set(theme_bw())

Outline of Today

  1. Announcements
  2. Linear Regression (3.1)
  3. Other considerations (3.3)
  4. Lab

Introduction to Linear Regression

Questions that we can answer

  1. Is there a relationship between X and Y
  2. How strong is the relationship between X and Y
  3. Which X are associated with Y
  4. How strong are X associated with Y
  5. How accurately can we predict Y
  6. Is the relationship linear
  7. Are there interactions between X

Simple Linear Regression

\[Y\approx \beta_0+\beta_1X\]

  • \(Y\) is our output
  • X is our Input
  • \(\beta\) are our coefficients or parameters
  • Y is approximately modeled as (intercept) + (slope) times X

\[\hat{y}=\hat{\beta_0}+\hat{\beta_1}x\]

Estimating Coefficients

\[e_i=y_i-\hat{y_i}\]

We define the residual sum of squares as \(\sum e_i^2\). We can show that the optimal estimators are:

\[\hat{\beta_1}=\frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}\] \[\hat{\beta_0}=\bar{y}-\hat{\beta_1}\bar{x}\]

These are the least squares coefficient estimates for simple linear regression.

Math if you are curious:

\[\frac{\delta RSS(\beta_0,\beta_1)}{\delta\beta_0}= \frac{\delta \sum (y_i-\beta_0x_i-\beta_1x_i)^2}{\delta\beta_0} =-2\sum(y_i-\beta_0-\beta_1x_i)\implies n\beta_0=\sum(y_i)-\sum(\beta_1x_i)\] \[\implies \beta_0=\bar{y}-\beta_1\bar{x}\]

\[\frac{\delta RSS(\beta_0,\beta_1)}{\delta\beta_0}= -2\sum(x_iy_i-\beta_0x_i-\beta_1x_i^2)= -2\sum x_iy_i-\bar{y}x_i+\beta_1\bar{x}x_i-\beta_1x_i^2 \implies \beta_1\sum(x_i^2-\bar{x}x_i)=\sum(x_iy_i-x_i\bar{y})\]

\[\implies \beta_1=\frac{\sum(x_iy_i-x_i\bar{y})}{\sum(x_i^2-\bar{x}x_i)} =\frac{\sum(xiyi)-n\bar{x}\bar{y}}{\sum(x_i^2)-n\bar{x}^2} =\frac{\sum(xiyi)-n\bar{x}\bar{y}+n\bar{x}\bar{y}-n\bar{x}\bar{y}}{\sum(x_i^2)-n\bar{x}^2-n\bar{x}^2+n\bar{x}^2} =\frac{\sum(xiyi)-\bar{x}\sum y_i+\sum\bar{x}\bar{y}-\bar{y}\sum x_i}{\sum(x_i^2)-2\bar{x}\sum(x_i)+\sum(\bar{x}^2)}\]

\[=\frac{\sum (xiyi-\bar{x} y_i+\bar{x}\bar{y}-\bar{y} x_i)}{\sum(x_i^2-2\bar{x}x_i+\bar{x}^2)} \implies \beta_1=\frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}\]

Accuracy of Coefficient Estimates

set.seed(5)
X <- runif(50, 0,10)
Y <- 5+X+ rnorm(50, mean=0, sd=1)
my_df <- data.frame(X,Y)
g_train_pts <- my_df %>% ggplot(aes(x = X, y = Y))+
  geom_point()+
  xlim(c(0,10))+
  ylim(c(4,15))
g_train_pts
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

lines<-data.frame(matrix(lm(Y~X,my_df)$coefficients,ncol=2)) # full set
for(i in 1:10){
  tempDF<-my_df[sample(nrow(my_df)),][1:20,]#Randomly sort, take first 20
  lines<-rbind(lines,data.frame(matrix(lm(Y~X,tempDF)$coefficients,ncol=2)))#Find slope
}

withlines<-g_train_pts+
  geom_abline(slope=1,intercept=5,color="red",size=2)+ #trueline
  geom_abline(slope=lines[1,2],intercept=lines[1,1],color="blue",size=2) #full set
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
withlines
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

withlines2<-withlines+
  geom_abline(data=lines,aes(slope=X2,intercept=X1))#random sample

withlines2
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

The OLS regression estimates will not be exact, but they are generally reasonable (unbiased).

Aside: Unbiased Estimators

An estimator is unbiased if on average we expect the estimate \(\hat{\mu}=\mu\). If we could average a huge number of estimates of \(\hat{\mu}\) from a huge number of sets then it would exactly equal \(\mu\) (Law of Large Numbers). It does not systematically over or underestimate

The OLS estimates are unbiased.

How good is \(\hat{\mu}\)? Standard error:

\[\text{Var}(\mu)=\text{SE}(\hat{\mu})^2=\frac{\sigma^2}{n}\]

Back to OLS

\[\text{SE}(\hat\beta_0)^2=\sigma^2\left[\frac{1}{n}+\frac{\bar{x}^2}{\sum(x_i-\bar{x})^2}\right]\] \[\text{SE}(\hat\beta_1)^2=\frac{\sigma^2}{\sum(x_i-\bar{x})^2}\]

where \(\sigma^2\)=Var(\(\epsilon\)). We assume that the variance is uncorrelated with x and the same for all points (homoskedastic). This is not always true, but tends to be a reasonable approximation. We obviously do not know \(\sigma^2\) but we can estimate it from the data using the residual standard error (strictly speaking, should have a hat, but it is normally dropped.) RSE=\(\sqrt{RSS/(n-2)}\).

Confidence Intervals:

95% confidence interval:\(\hat{\beta_1}\pm2\cdot\text{SE}(\hat{\beta_1})\). This is approximate (guassian errors, other assumptions, but it is close enough for this course right now. The 2 should be the 97.5% quantile of a t distribution for n-2 degrees of freedom)

Confidence tests:

normal test: \(H_0:\) there is no relationship between X and Y \(H_a\): there is a relationship between X and Y

\(H_0:\beta_1=0\)

\(H_a:\beta_1\neq0\)

t statistic:

\[t=\frac{\hat{\beta_1}-0}{SE(\hat{\beta_1})}\]

Basically, how many standard deviations is \(\hat{\beta_1}\) away from 0. If there is no relationship, we expect to follow a t distribution with n-2 degrees of freedom, which is approximately standard normal when n\(\geq 30\).

P value: the probability that \(|t|\) or larger would have been observed assuming \(\beta_1=0\).

Assessing the Accuracy of the model

  • Residual Standard Error (RSE):
    • \(\sqrt{\frac{1}{n-2}\sum(y_i-\hat{y})^2}\)
    • a measure of the lack of fit.
    • Units of Y
  • \(R^2\):
    • Alternative measure of fit
    • proportion of variance explained [0,1]
    • \(R^2=\frac{TSS-RSS}{TSS}=1-\frac{RSS}{TSS}=1-\frac{\sum(y_i-\hat{y})^2}{\sum(y_i-\bar{y})^2}\)
    • 1- Variance of fit (RSS)/Total variance
    • a “good” \(R^2\) depends on context
    • in simple Linear Regression \(R^2=r^2=Cor(X,Y)^2\)
    • tends to be biased, so we use the adjusted\(R^2=1-\frac{RSS}{TSS}\frac{n-1}{n-p-1}\) in practice. (although it is still overly optimistic)
  • F-Statistic: Thursday

Other Considerations (3.3)

Qualitative Predictors

Normally we talk about quantitative predictors, but qualitative ones can be used as well.

  • 2 levels:
    • Use an indicator dummy variable
    • 0/1 indicator: (e.g. M/F would be 0/1)
      • Intercept is then the average among those in class 0
      • \(\beta_1\) is the average difference between class 0 and 1
      • 1/0 is arbitrary, no effect on fit, but will effect interpretation of the fits
    • -1/1 indicator: (e.g. M/F would be -1/1)
      • intercept is average
      • \(\beta_1\) is the amount above/below average for the groups.
    • Final predictions are the same either way
    • We will (mostly) use 0/1, it’s easier.
  • More than 2 variables:
    • Use classes-1 dummy variables (0/1)
    • the rest is approximately the same as above

Extensions of the Linear Model

We have so far assumed that the relationship between predictors and response are additive and linear. We will cover this later in the course, although examples are below:

  • Interaction terms: \(\beta_3X_1X_2\)
    • hierarchical principal: if you include an interaction, you should include the main effects regardless of p values of their coefficients.
  • Non-linear relationships: polynomial regression \(\beta_2X_1^2\)
    • This is still a linear model!
    • So we can still do OLS

Potential Problems

Identification and fixes are more of an art than a science. Tend to identify with residual plots.

  1. Non-linearity
    • identify: U shape resid plot
    • fix with: nonlinear transformations of input variables
  2. Correlation of error terms
    • extreme example: doubled points. narrows confidence interval by \(\sqrt{2}\)
    • often found in time series
    • identify: knowing an error gives idea of next error
    • fix: not in this class
  3. Non-constant error terms
    • heteroscedasticity
    • identify: funnel shape
    • fix: transform response variable
  4. Outliers
    • identify: Significantly different from other values (studentized residuals: divide by standard error, expected between [-3,3])
    • fix: can potentially remove if you can justify. Can also indicate missing predictor
  5. High-leverage points
    • identify: predictor value outside normal range of observations
      • leverage statistic: \(h_i=\frac{1}{n}+\frac{(x_i-\bar{x})^2}{\sum(x_j-\bar{x})^2}\)
      • always between 1/n and 1
      • if significantly higher than (p+1)/n it has high leverage
  6. Collinearity
    • two+ predictor variables are closely related
    • issue: hard to separate out effects, much more uncertainty
    • power of the hypothesis test is reduced
    • identify: correlation matrix (only works for 2)
    • multicollinearity: 3+ variables related
    • Variance inflation factor (VIF): \(VIF(\hat{\beta_j})=\frac{1}{1-R^2_{X_j|X_{-j}}}\)
      • Variance of \(\beta_j\) when fitting the full model divided by variance when fit on its own
      • Smallest: 1
      • Rule of thumb: \(>5\) or 10 is a problem

Diagnostic plots

Code taken from Professor Nate Wells.

set.seed(700)
X <- runif(80, 0, 1)
e <- rnorm(80, 0, .25)
Y <- 1 + 2*X + e
my_data <- data.frame(X,Y)
ggplot(my_data, aes(x = X , y = Y)) + geom_point()

Question 1 What is the underlying model? error term?

my_mod<-lm(Y ~ X, data = my_data)
gglm(data = my_mod)

Question 2 For each diagnostic image, what is represented along each axis? Why? What should we look for?

Lab: Simple Linear Regression

Lab: Partner work

Based on the lab section of ISLR

Import the Boston Dataset from ISLR2

Question 3 Using the ?Boston command, tell me about the dataset

Fitting a Linear Model in R

We can fit a linear model in R using the lm function. Syntax: lm(y~x,data=DataFrame)

lm.fit <- lm(lstat ~ age , data = bostonData)

We can then call it directly to get basic information

lm.fit
## 
## Call:
## lm(formula = lstat ~ age, data = bostonData)
## 
## Coefficients:
## (Intercept)          age  
##      2.1744       0.1528

Additional functions that are useful:

  • summary(): gives significantly more information about the lm
  • names(): gives you access to the underlying values, e.g. lm.fit$coefficients gives the coefficients
  • confint(): confidence interval of coefficient intervals
  • predict(): confidence and prediction intervals for the predicted value
  • residuals(): pull the residuals from the linear model
  • rstudent: get the studentized residuals
  • hatvalues(): leverage statistics
  • gglm(): diagnostic plots

Question 4 Using the Boston Dataset

  1. Fit a linear model where x=lstat and y=medv

  2. Report the estimates, error, and p values for the two coeficients.

  3. Using geom_point() and geom_abline(), plot the datapoints for lstat and medv, and plot the OLS line.

  4. Plot the residuals using ggplot2

  5. Use gglm to create the diagnostic plots

  6. Based on the plots, do you believe that a linear model is appropriate? Why or why not?

  7. Did you use X or Y for the residual value plots? Why? what happens if you use the other variable? When would you want to use one over the other?

Non-linear Transformations

Syntax: I(), allows for transforming variables within lm

Syntax: poly(x,power), creates a polynomial fit. By default, it orthogonalizes the predictors. use raw=TRUE to get the raw higher order terms.

For example:

lm.fit <- lm(lstat ~ age+I(age^2) , data = bostonData)

lm.fit
## 
## Call:
## lm(formula = lstat ~ age + I(age^2), data = bostonData)
## 
## Coefficients:
## (Intercept)          age     I(age^2)  
##    8.177692    -0.108972     0.002175
lm.fit <- lm(lstat ~ poly(age,2,raw=TRUE) , data = bostonData)

lm.fit
## 
## Call:
## lm(formula = lstat ~ poly(age, 2, raw = TRUE), data = bostonData)
## 
## Coefficients:
##               (Intercept)  poly(age, 2, raw = TRUE)1  
##                  8.177692                  -0.108972  
## poly(age, 2, raw = TRUE)2  
##                  0.002175

Question 5

  1. Add a quadratic term to your model from Question 4. Does this improve the model? Why?
  2. Use anova(Question2Model, Question3a model) to verify your answer to part a.
  3. use poly() to fit up to a 5th degree polynomial. What degree polynomial do you think is most appropriate? Why?
  4. use poly() to fit that polynomial, one higher, two higher, and one lower. What do you notice about the significance (p value)? Which would you use?

Fitting Linear Models

Question 6 Using the code from earlier as a start, explore how linear regressions change (get more/less accurate):

set.seed(5)
X <- runif(50, 0,10)
Y <- 5+X+ rnorm(50, mean=0, sd=1)
my_df <- data.frame(X,Y)
  1. As the number of points used to fit the line increases (from 2 to 50+)
  2. As the error \(\epsilon\) changes
  3. As we allow for duplicated points (you may need to figure out how to make this happen)
  4. If the error (rnorm) changes to a different distribution (e.g. runif)
  5. Are there any interactions between these changes?
  6. Do any of these changes cause the estimates to be biased in a particular way? Or are the estimated parameters still unbiased?

Wrap-up

Lab goals:

In this lab we practiced the following:

  • Simple Linear Regression
  • Other Considerations to Linear Regression
  • Diagnostic Plots
  • helpful functions when running regressions in R

Course Schedule:

  1. Today: Simple Linear Regression
  2. Tomorrow: Multiple Linear Regression

Reminders for next class:

  • Homework 2 is due Friday at 10pm
  • Labs 5,6 are dur Friday at 10pm
  • The reading assignment for Chapter 3.2 of the text is due Wednesday at 10pm