class: center, middle, inverse, title-slide # Lecture 11 ## Multiple linear regression ### Jung-Jin Lee ### Apr 7, 2020 --- ## 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_11_files/figure-html/unnamed-chunk-4-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_11_files/figure-html/unnamed-chunk-8-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 ``` --- ## Example: weight/girth of sheep - The best fitting line equation: `$$y = -46.041 + 1.043x$$` -- - Interpretation: the live weight of a sheep with `\(x\)` cm chest girth is approximately `\(-46.041 + 1.043x\)` kg. -- - For example, if there is a sheep of which chest girth is 85 cm, we can predict that its live weight would be approximately $$-46.041 + 1.043 \times 85 = 42.614 \text{ kg} $$ -- - What if there is another sheep with chest girth 86 cm? The corresponding live weight is approximately $$-46.041 + 1.043 \times 86 = 43.657 \text{ kg} $$ which is 1.043 more than previous. -- - In general `\(\hat{\beta}_1\)` is interpreted as: a unit increase in `\(x\)` results in the increase of `\(y\)` by `\(\hat{\beta}_1\)`. --- ## Making inference on `\(\beta_0\)` and `\(\beta_1\)` - Point estimates of `\(\beta_0\)` and `\(\beta_1\)` were already covered: ```r fit$coefficients ``` ``` ## (Intercept) sheep$girth ## -46.041373 1.043237 ``` -- - How about 95% confidence of `\(\beta_0\)` and `\(\beta_1\)`? ```r *confint(fit) ``` ``` ## 2.5 % 97.5 % ## (Intercept) -52.596253 -39.486493 ## sheep$girth 0.955333 1.131141 ``` -- We are 95% certain that `\(\beta_1\)` is between 0.955333 and 1.1311406. --- ## Making inference on `\(\beta_0\)` and `\(\beta_1\)` - Statistical test on `\(\beta_1\)` - `\(H_0\)`: `\(\beta_1=0\)` - `\(p\)`-value is obtained from `coefficients` of `summary` of `lm` output ```r summary(fit)$coefficients ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -46.041373 3.28116478 -14.03202 3.364110e-21 ## sheep$girth 1.043237 0.04400187 23.70892 2.112011e-33 ``` - With an extremely small `\(p\)`-value `\(2.1120107\times 10^{-33}\)`, we conclude that `\(\beta_1\neq 0\)`. - Interpretation of this conclusion will be explained later --- ## Multiple linear regression Read in data set [donkey.tsv](donkey.tsv). ```r donkey <- read.table("donkey.tsv", header = T, sep = "\t") dim(donkey) ``` ``` ## [1] 386 8 ``` ```r head(donkey) ``` ``` ## Donkey Sex Age Bodywt Heartgirth Umbgirth Length Height ## 1 s1 male 3.5 128 111 126 95 106 ## 2 s2 male 2.0 102 106 120 75 101 ## 3 s3 male 3.0 100 103 125 76 101 ## 4 s4 male 3.0 92 102 121 73 101 ## 5 s5 male 3.0 98 105 121 70 101 ## 6 s6 male 4.0 110 101 126 80 101 ``` --- ## Multiple linear regression - Donkey: donkey ID - Sex: sex - Age: age in years - Bodywt: body weight in kg - Heartgirth: girth at the level of the heart (cm) - Umbgirth: girth at the umbilicus (cm) - Length: length from the olecranon to the tuber ischii (cm) - Height: height at the withers (cm) -- Researchers were interested in estimating the **weight** of donkeys using other measurements that can be obtained more easily such as girth at the level of the heart (cm), girth at the umbilicus (cm), length from the olecranon to the tuber ischii (cm), height at the withers (cm), and the age (years) and sex. --- ## Multiple linear regression - Simple linear regression using Length ```r ggplot(donkey, aes(Length, Bodywt)) + geom_point() ``` <img src="Drexel_2020_Lecture_11_files/figure-html/unnamed-chunk-18-1.png" width="50%" style="display: block; margin: auto;" /> --- ## Multiple linear regression ```r fit_Length <- lm(Bodywt ~ Length, data = donkey) summary(fit_Length) ``` ``` ## ## Call: ## lm(formula = Bodywt ~ Length, data = donkey) ## ## Residuals: ## Min 1Q Median 3Q Max ## -45.243 -10.584 -1.357 10.625 62.075 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -115.0015 10.0375 -11.46 <2e-16 *** ## Length 2.8940 0.1219 23.75 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 15.82 on 384 degrees of freedom ## Multiple R-squared: 0.5949, Adjusted R-squared: 0.5939 ## F-statistic: 564 on 1 and 384 DF, p-value: < 2.2e-16 ``` --- ## Multiple linear regression ```r ggplot(donkey, aes(Length, Bodywt)) + geom_point() + geom_smooth(method = "lm", se = F) ``` <img src="Drexel_2020_Lecture_11_files/figure-html/unnamed-chunk-20-1.png" width="50%" style="display: block; margin: auto;" /> --- ## Multiple linear regression - Simple linear regression using Height ```r ggplot(donkey, aes(Height, Bodywt)) + geom_point() ``` <img src="Drexel_2020_Lecture_11_files/figure-html/unnamed-chunk-21-1.png" width="50%" style="display: block; margin: auto;" /> --- ## Multiple linear regression ```r fit_Height <- lm(Bodywt ~ Height, data = donkey) summary(fit_Height) ``` ``` ## ## Call: ## lm(formula = Bodywt ~ Height, data = donkey) ## ## Residuals: ## Min 1Q Median 3Q Max ## -105.241 -9.825 -0.825 8.505 75.131 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -137.2379 11.7708 -11.66 <2e-16 *** ## Height 2.5055 0.1132 22.13 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 16.48 on 384 degrees of freedom ## Multiple R-squared: 0.5605, Adjusted R-squared: 0.5594 ## F-statistic: 489.8 on 1 and 384 DF, p-value: < 2.2e-16 ``` --- ## Multiple linear regression ```r ggplot(donkey, aes(Height, Bodywt)) + geom_point() + geom_smooth(method = "lm", se = F) ``` <img src="Drexel_2020_Lecture_11_files/figure-html/unnamed-chunk-23-1.png" width="50%" style="display: block; margin: auto;" /> --- ## Multiple linear regression - In multiple linear regression, we use more than one explanatory variables. -- - For example, we assume that `$$Bodywt_i=\beta_0+\beta_1 Length_i +\beta_2 Heartgirth_i +\beta_3 Height_i + e_i,$$` where `\(e_i\)` is a random variable following `\(N(0,\sigma^2)\)`. -- - Again `lm()` is used for fitting --- ## Multiple linear regression ```r fit <- lm(Bodywt ~ Length + Heartgirth + Height, data = donkey) summary(fit)$coefficients ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -217.7061239 7.6606181 -28.418872 2.913804e-96 ## Length 0.9163415 0.1193591 7.677181 1.375154e-13 ## Heartgirth 2.1241803 0.1017086 20.884954 3.681853e-65 ## Height 0.2608338 0.1117453 2.334182 2.010392e-02 ``` -- Based on the output, we conclude that `\begin{align} Bodywt = &-217.706 \\ &+0.916 \times Length \\ &+2.124 \times Heartgirth \\ &+0.261 \times Height \end{align}` Interpretation: 0.916 represents the change in Bodywt for a unit change in Length, when **the other explanatory variables, Heartgirth and Height, are held constant** (i.e. after controlling for these variables) --- ## Model comparison - Consider a model `$$\textbf{Model 1: } \texttt{Bodywt}=\beta_0+\beta_1 \texttt{Length}.$$` How good is this model? -- - One approach to assess a model is to compute `\(R^2\)`. - A good model will yield an `\(R^2\)` close to 1. ```r fit1 <- lm(Bodywt ~ Length, data = donkey) ## Model 1 r.squared1 <- summary(fit1)$r.squared r.squared1 ``` ``` ## [1] 0.5949296 ``` --- ## Model comparison - Consider another model `$$\textbf{Model 2: } \texttt{Bodywt}=\beta_0+\beta_1 \texttt{Height}.$$` ```r fit2 <- lm(Bodywt ~ Height, data = donkey) ## Model 2 r.squared2 <- summary(fit2)$r.squared r.squared2 ``` ``` ## [1] 0.5605374 ``` - **Model 1** has bigger `\(R^2\)` than **Model 2**, so we can conclude that **Model 1** is a better model than **Model 2**. --- ## Model comparison .pull-left[ <img src="Drexel_2020_Lecture_11_files/figure-html/unnamed-chunk-27-1.png" style="display: block; margin: auto;" /> ] .pull-right[ <img src="Drexel_2020_Lecture_11_files/figure-html/unnamed-chunk-28-1.png" style="display: block; margin: auto;" /> ] --- ## Model comparison - Now consider the model `$$\textbf{Model 3: } \texttt{Bodywt}=\beta_0+\beta_1\texttt{Length}+\beta_2\texttt{Height}.$$` How good is this model? -- ```r fit3 <- lm(Bodywt ~ Length + Height, data = donkey) ## Model 3 r.squared3 <- summary(fit3)$r.squared r.squared3 ``` ``` ## [1] 0.6669952 ``` - This has the biggest `\(R^2\)` so far. Is this surprising? -- - **Model 1** and **Model 2** are _nested_ within **Model 3**, so `\(R^2\)` being the bigger in **Model 3** than **Model 1** or **Model 2** is expected. -- - In general, the more predictors are added, the bigger `\(R^2\)` is. --- ## Model comparison - Now consider yet another model `\begin{align} \textbf{Model 4: } \texttt{Bodywt} = & \beta_0 + \beta_1\texttt{Length} + \beta_2\texttt{Height} \\ & +\beta_3\texttt{Heartgirth} + \beta_4\texttt{Umbgirth} \end{align}` ```r fit4 <- lm(Bodywt ~ Length + Height + Heartgirth + Umbgirth, data = donkey) ## Model 4 r.squared4 <- summary(fit4)$r.squared r.squared4 ``` ``` ## [1] 0.8561892 ``` -- - Although this model the biggest `\(R^2\)` so far, it is more complex than other models. -- - Need to consider trade-off between simple model (less number of variables) vs good fit (higher `\(R^2\)`) -- - Idea: if a bigger model does not considerably increase `\(R^2\)`, then stick to simpler model.