class: center, middle, inverse, title-slide # Lecture 10 ## Introduction to linear regression ### Jung-Jin Lee ### Mar 31, 2020 --- ## Simple linear regression - In simple linear regression, we describe the relationship between two numerical variables, `\(x\)` and `\(y\)`, by determining the straight line that approximates the data points on a scatter diagram most closely. -- - We call `\(y\)` the **outcome**, **dependent**, or **response** variable and we represent this on the vertical axis of the scatter diagram. -- - We call `\(x\)` the **explanatory**, **independent**, or **predictor** variable, or **regressor**. --- ## Simple linear regression - We regard the `\(x\)` variable as one whose values can be measured without error. On the other hand, the `\(y\)` variable is a random variable which is subject to experimental variation. -- - To be more precise, we assume that `\(y_i=\beta_0+\beta_1 x_i+e_i\)`, where `\(e_i\)` is a random variable following `\(N(0,\sigma^2)\)`. -- - Roughly speaking, simple linear regression is a combination of parameter estimation and best fitting line. -- - We will estimate the unknown true parameters `\(\beta_0\)`, `\(\beta_1\)`, and `\(\sigma^2\)` based on our sample, just like we estimate the population mean `\(\mu\)` and population variance `\(\sigma^2\)` with the sample mean `\(\bar{X}\)` and the sample variance `\(S^2\)`, respectively. --- ## Least square line ```r x <- c(2, 4, 5, 9, 12); y <- c(3, 2, 4, 6, 5) ``` <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-3-1.png" width="50%" style="display: block; margin: auto;" /> Want to find the line that best describes the trend --- ## Least square line .pull-left[ <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-4-1.png" style="display: block; margin: auto;" /> ] .pull-right[ <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-5-1.png" style="display: block; margin: auto;" /> ] --- ## Least square line .pull-left[ <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-6-1.png" style="display: block; margin: auto;" /> ] .pull-right[ <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-7-1.png" style="display: block; margin: auto;" /> ] Determine the `\(y\)`-intercept and slope so that the sum of squares of the <font color="blue"> blue dashed line segment </font> is minimized. --- ## Fitting with `lm()` ```r fit <- lm(y ~ x) summary(fit) ``` ``` ## ## Call: ## lm(formula = y ~ x) ## ## Residuals: ## 1 2 3 4 5 ## 0.3497 -1.2638 0.4294 1.2025 -0.7178 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.0368 1.0330 1.972 0.143 ## x 0.3067 0.1406 2.182 0.117 ## ## Residual standard error: 1.135 on 3 degrees of freedom ## Multiple R-squared: 0.6135, Adjusted R-squared: 0.4847 ## F-statistic: 4.762 on 1 and 3 DF, p-value: 0.1171 ``` --- ## Fitting with `lm()` <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-9-1.png" width="80%" style="display: block; margin: auto;" /> --- ## Interpretation of `lm()` output `lm()` returns lots of information regarding the regression model: ```r names(fit) ``` ``` ## [1] "coefficients" "residuals" "effects" "rank" ## [5] "fitted.values" "assign" "qr" "df.residual" ## [9] "xlevels" "call" "terms" "model" ``` These _names_ can be extracted using `$` sign, e.g. `fit$residuals`. --- ## Interpretation of `lm()` output - `\(\hat{\beta}_0\)`: fitted intercept - `\(\hat{\beta}_1\)`: fitted slope ```r *fit$coefficients ``` ``` ## (Intercept) x ## 2.0368098 0.3067485 ``` ```r beta0 <- fit$coefficients[1]; beta1 <- fit$coefficients[2] beta0 ``` ``` ## (Intercept) ## 2.03681 ``` ```r beta1 ``` ``` ## x ## 0.3067485 ``` --- ## Interpretation of `lm()` output - `\(\hat{y}_i\)`: fitted value `$$\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i\quad\text{for each }i$$` ```r fitted_values <- beta0 + beta1*x fitted_values ``` ``` ## [1] 2.650307 3.263804 3.570552 4.797546 5.717791 ``` -- ```r *fit$fitted.values ``` ``` ## 1 2 3 4 5 ## 2.650307 3.263804 3.570552 4.797546 5.717791 ``` --- ## Interpretation of `lm()` output - Residual: the difference between observation and fitted value (blue dashed lines): `$$r_i = y_i- \hat{y}_i = y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i)\quad\text{for each }i$$` ```r residuals <- y - fitted_values residuals ``` ``` ## [1] 0.3496933 -1.2638037 0.4294479 1.2024540 -0.7177914 ``` -- ```r *fit$residuals ``` ``` ## 1 2 3 4 5 ## 0.3496933 -1.2638037 0.4294479 1.2024540 -0.7177914 ``` --- ## Interpretation of `lm()` output - RSS (residual sum of squares): sum of the squares of residuals: `$$RSS=\sum_i r_i^2 = \sum_i\left[y_i - \underbrace{(\hat{\beta}_0 + \hat{\beta}_1 x_i)}_{\text{fitted value}}\right]^2$$` ```r rss <- sum(fit$residuals^2) rss ``` ``` ## [1] 3.865031 ``` --- ## Interpretation of `lm()` output - SSreg (regression sum of squares): sum of squares of the difference between fitted value and `\(\bar{y}\)`: `$$SSreg = \sum_i \left[(\hat{\beta}_0 + \hat{\beta}_1 x_i)-\bar{y}\right]^2$$` ```r ssreg <- sum((fit$fitted.values-mean(y))^2) ssreg ``` ``` ## [1] 6.134969 ``` --- ## Interpretation of `lm()` output - TSS (total sum of squares): sum of the squares of the difference between each `\(y_i\)` value and the mean `\(\bar{y}\)`: `$$TSS=\sum_i\left(y_i-\bar{y}\right)^2$$` ```r tss <- sum((y-mean(y))^2) tss ``` ``` ## [1] 10 ``` --- ## Interpretation of `lm()` output - It is known that RSS + SSreg = TSS: https://en.wikipedia.org/wiki/Partition_of_sums_of_squares ```r rss + ssreg ``` ``` ## [1] 10 ``` ```r tss ``` ``` ## [1] 10 ``` -- - `\(R^2\)` (R-squared or coefficient of determination): `$$R^2 = \frac{SSreg}{TSS} = 1 - \frac{RSS}{TSS}$$` -- - `\(R^2\)` is a measure of goodness of model (model is good if `\(R^2\)` is close to 1, bad if `\(R^2\)` is close to 0) --- ## Interpretation of `lm()` output ```r r2 <- ssreg/tss r2 ``` ``` ## [1] 0.6134969 ``` -- ```r names(summary(fit)) ``` ``` ## [1] "call" "terms" "residuals" "coefficients" ## [5] "aliased" "sigma" "df" "r.squared" ## [9] "adj.r.squared" "fstatistic" "cov.unscaled" ``` -- ```r summary(fit)$r.squared ``` ``` ## [1] 0.6134969 ``` --- ## Interpretation of `lm()` output - Estimate of `\(\sigma\)`: it is well known that `$$\frac{RSS}{n-2}$$` is an _unbiased_ estimator of `\(\sigma^2\)`, where `\(n\)` is the sample size. -- ```r sqrt(rss/(length(x)-2)) ``` ``` ## [1] 1.135052 ``` -- ```r summary(fit)$sigma ``` ``` ## [1] 1.135052 ``` --- ## Example: weight/girth of sheep Download a tab-delimited file [sheep.tsv](sheep.tsv). It consists of 66 observations and 2 variables. ```r sheep <- read.table("sheep.tsv", header = T, sep = "\t") dim(sheep) ``` ``` ## [1] 66 2 ``` .pull-left[ ```r head(sheep) ``` ``` ## weight girth ## 1 30 76 ## 2 24 71 ## 3 20 63 ## 4 25 69 ## 5 25 67 ## 6 19 62 ``` ] .pull-right[ - `weight`: live weight of a sheep (in kg) - `girth`: chest girth (in cm) ] --- ## Example: weight/girth of sheep Strong linear trend is visible: ```r ggplot(sheep, aes(girth, weight)) + geom_point() ``` <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-27-1.png" width="50%" style="display: block; margin: auto;" /> --- ## Example: weight/girth of sheep - Find the best fitting line and overlay onto the scatter plot using `geom_abline()` -- ```r fit <- lm(sheep$weight ~ sheep$girth) # or fit <- lm(weight ~ girth, data = sheep) ``` -- ```r fit$coefficients ``` ``` ## (Intercept) sheep$girth ## -46.041373 1.043237 ``` -- ```r beta0 <- fit$coefficients[1]; beta1 <- fit$coefficients[2] beta0; beta1 ``` ``` ## (Intercept) ## -46.04137 ``` ``` ## sheep$girth ## 1.043237 ``` --- ## Example: weight/girth of sheep ```r ggplot(sheep, aes(girth, weight)) + geom_point() + geom_abline(intercept = beta0, slope = beta1, color = "red") ``` <img src="Drexel_2020_Lecture_10_files/figure-html/unnamed-chunk-31-1.png" width="50%" style="display: block; margin: auto;" /> --- ## Example: weight/girth of sheep - Find RSS, SSreg, and TSS -- ```r rss <- sum((fit$residuals)^2) rss ``` ``` ## [1] 452.3425 ``` -- ```r ssreg <- sum((fit$fitted.values - mean(sheep$weight))^2) ssreg ``` ``` ## [1] 3972.93 ``` -- ```r tss <- rss + ssreg # or tss <- sum((sheep$weight - mean(sheep$weight))^2) tss ``` ``` ## [1] 4425.273 ``` --- ## Example: weight/girth of sheep - Compute `\(R^2\)` -- ```r r2 <- ssreg/tss; print(r2) ``` ``` ## [1] 0.897782 ``` ```r # or alternatively summary(fit)$r.squared ``` ``` ## [1] 0.897782 ``` -- - Estimate `\(\sigma\)` -- ```r sqrt(rss/(length(sheep$girth)-2)) ``` ``` ## [1] 2.658543 ``` ```r summary(fit)$sigma # alternatively ``` ``` ## [1] 2.658543 ```