--- title: "Regression revisited" editor: markdown: wrap: 72 --- ## Regression - Use regression when one variable is an outcome (*response*, $y$). - See if/how response depends on other variable(s), *explanatory*, $x_1, x_2,\ldots$. - Can have *one* or *more than one* explanatory variable, but always one response. - Assumes a *straight-line* relationship between response and explanatory. - Ask: - *is there* a relationship between $y$ and $x$'s, and if so, which ones? - what does the relationship look like? ## Packages ```{r bRegression-1, results='hide'} library(MASS, exclude = "select") # for Box-Cox, later library(tidyverse) library(broom) library(marginaleffects) ``` ## A regression with one $x$ - 13 children, measure average total sleep time (ATST, mins) and age (years) for each. - See if ATST depends on age. - Data in `sleep.txt`, ATST then age, values separated by a space. - Read in data: ```{r bRegression-2 } my_url <- "http://datafiles.ritsokiguess.site/sleep.txt" sleep <- read_delim(my_url, " ") ``` ## Check data \footnotesize ```{r} sleep ``` \normalsize Make scatter plot of ATST (response) vs. age (explanatory). ## The scatterplot ```{r suggo, fig.height=6} ggplot(sleep, aes(x = age, y = atst)) + geom_point() ``` ## Plot with smooth trend ```{r icko,fig.height=6} #| echo: false ggplot(sleep, aes(x = age, y = atst)) + geom_point() + geom_smooth() ``` ## The regression - Scatterplot shows no obvious curve, and a pretty clear downward trend. So we can run the regression: ```{r bRegression-6, echo=-1} # options(width = 60) sleep.1 <- lm(atst ~ age, data = sleep) ``` - and look at the output via ```{r} #| eval: false summary(sleep.1) ``` (over) ## The output \small ```{r bRegression-7} #| echo: false summary(sleep.1) ``` \normalsize ## Conclusions - The relationship appears to be a straight line, with a downward trend. - $F$-tests for model as a whole and $t$-test for slope (same) both confirm this (P-value $5.7\times 10^{-7}=0.00000057$). - Slope is $-14$, so a 1-year increase in age goes with a 14-minute decrease in ATST on average. - Here R-squared is 0.9054, pleasantly high. ## Doing things with the regression output - Output from regression (and eg. $t$-test) is all right to look at, but hard to extract and re-use information from. - Package `broom` extracts info from model output in way that can be used with pipe (later): ```{r bRegression-8, size="footnotesize"} tidy(sleep.1) ``` ## also one-line summary of model ```{r} #| echo: false w <- getOption("width") options(width = 60) ``` ```{r bRegression-9} glance(sleep.1) ``` ```{r} #| echo: false options(width = w) ``` ## `augment` Useful for plotting residuals: \footnotesize ```{r bRegression-10,warning=F} sleep.1 %>% augment(sleep) -> sleep.1a sleep.1a ``` \normalsize ## Leverage - The values in `.hat` say how much *leverage* an observation has, based on its explanatory variable values: that is, how much potential it has to affect the fit. - A `.hat` value bigger than $2p/n$ is considered unusually large, where: - $p$ is number of parameters estimated (here 2, intercept and slope) - $n$ is number of observations (here 13). ## The top leverages here \footnotesize ```{r} sleep.1a %>% select(atst, age, .hat) %>% mutate(cutoff = 2 * 2 / 13) %>% slice_max(.hat, n = 5) ``` \normalsize - The top two leverages are unusually large - These are the highest and lowest ages. ## Cook's distance - in `.cooksd` - measures the *influence* an observation has on the fit - combination of leverage and residual: if both large in size, Cook's distance large ## The top Cook's distances ```{r} sleep.1a %>% select(-.sigma, -.std.resid) %>% slice_max(.cooksd, n = 5) ``` - The largest Cook's distances are: - the oldest child (largest age) - an oldish child that has a large (negative) residual ## Comments - Large residual (in size): unusual *response* value - Large leverage: unusual *explanatory* value - Large Cook's distance: influential observation on fit - might be large residual or large leverage or both. - In multiple regression, a large leverage means that the explanatory variable values *in combination* are unusual. ## CI for mean response and prediction intervals Once useful regression exists, use it for prediction: - To get a single number for prediction at a given $x$, substitute into regression equation, eg. age 10: predicted ATST is $646.48-14.04(10)=506$ minutes. - To express uncertainty of this prediction, *CI for mean response* expresses uncertainty about mean ATST for all children aged 10, based on data. - Also do above for a child aged 5. ## The `marginaleffects` package 1/2 To get predictions for specific values, set up a dataframe with those values first: ```{r} new <- datagrid(model = sleep.1, age = c(10, 5)) new ``` Any variables in the dataframe that you don't specify are set to their mean values (quantitative) or most common category (categorical). ## The `marginaleffects` package 2/2 Then feed into `newdata` in `predictions`. This contains a lot of columns, so you probably want only to display the ones you care about: ```{r} cbind(predictions(sleep.1, newdata = new)) %>% select(estimate, conf.low, conf.high, age) %>% mutate(length = conf.high - conf.low) ``` The confidence limits are a 95% confidence interval for the mean response at that `age`. ## Plotting the confidence intervals for mean response again: ```{r} #| fig-height: 5 plot_predictions(sleep.1, condition = "age") ``` ## Comments - Age 10 closer to centre of data, so intervals are both narrower than those for age 5. - That is, a high-leverage observation will have a less accurate prediction. ## That grey envelope Marks confidence interval for mean for all $x$: ```{r bRegression-15, fig.height=4.5} ggplot(sleep, aes(x = age, y = atst)) + geom_point() + geom_smooth(method = "lm") + scale_y_continuous(breaks = seq(420, 600, 20)) ``` ## Diagnostics How to tell whether a straight-line regression is appropriate? - Before: check scatterplot for straight trend. - After: plot *residuals* (observed minus predicted response) against predicted values. Aim: a plot with no pattern. ## Residual plot Not much pattern here --- regression appropriate. ```{r akjhkadjfhjahnkkk,fig.height=4.5} ggplot(sleep.1, aes(x = .fitted, y = .resid)) + geom_point() ``` ## Ways to handle curves - Change $x$ (adding $x^2$) - Another way: change $y$ (transformation). - Can guess how to change $y$, or might be theory: - example: relationship $y=ae^{bx}$ (exponential growth): - take logs to get $\ln y=\ln a + bx$. - Taking logs has made relationship linear ($\ln y$ as response). - Or, *estimate* transformation, using Box-Cox method. - in package `MASS`, load with `library(MASS, exclude = "select")` ## Ice crystals - Ice crystals are introduced into a chamber that is kept at a fixed temperature of $-5$°C. The growth of the crystals over time is observed. Our data are measurements of mass `m` of an ice crystal in nanograms at time `t` in seconds after being added to the chamber. Each ice crystal was measured once, at a "random" time. - Model relationship between mass and time. - It is suggested by the scientist that `m` is actually related to `log(t)` in some way. ## Read in the data \small ```{r} my_url <- "http://datafiles.ritsokiguess.site/ice_crystals.csv" ice_crystal <- read_csv(my_url) ice_crystal ``` \normalsize ## A graph ```{r} #| fig-height: 4 ggplot(ice_crystal, aes(x = log(t), y = m)) + geom_point() + geom_smooth() ``` Suggestion of curve and maybe fanning out. ## Fit a regression and look at residuals ```{r} #| fig-height: 4.5 ice.1 <- lm(m ~ log(t), data = ice_crystal) ice.1 %>% augment(ice_crystal) -> ice.1a ggplot(ice.1a, aes(x = .fitted, y = .resid)) + geom_point() ``` ## Comments - The residual plot shows evidence of both a curve and fanning out - Suggests transformation of response, in this case log: ```{r} #| fig-height: 5 boxcox(m ~ log(t), data = ice_crystal) ``` ## Regression with transformed response ```{r} #| eval: false ice.2 <- lm(log(m) ~ log(t), data = ice_crystal) summary(ice.2) ``` Output on next slide. ## Regression output ```{r} #| echo: false ice.2 <- lm(log(m) ~ log(t), data = ice_crystal) summary(ice.2) ``` ## Are the residuals better now? ```{r} #| fig-height: 5 ice.2 %>% augment(ice_crystal) -> ice.2a ggplot(ice.2a, aes(x = .fitted, y = .resid)) + geom_point() ``` ## Comments - The residuals look much more random; no longer evidence of fanning out. - What does having a linear relationship between `log(m)` and `log(t)` mean? $$\ln(m) = a + b \ln(t)$$ Exp both sides: $$ m = e^{a + b \ln(t)}$$ Simplify the right: $$ m = e^a e^{b \ln(t)} $$ $$ m = A t^b $$ where $A = \exp(a)$. ## Multiple regression - What if more than one $x$? Extra issues: - Now one intercept and a slope for each $x$: how to interpret? - Which $x$-variables actually help to predict $y$? - Different interpretations of "global" $F$-test and individual $t$-tests. - R-squared no longer correlation squared, but still interpreted as "higher better". - In `lm` line, add extra $x$s after `~`. - Interpretation not so easy (and other problems that can occur). ## The punting data Data set `punting.txt` contains 4 variables for 13 right-footed football kickers (punters): left leg and right leg strength (lbs), distance punted (ft), average leg strength `avg`. Predict punting distance from other variables. ## Reading in - Separated by *multiple spaces* with *columns lined up*: ```{r bRegression-46 } my_url <- "http://datafiles.ritsokiguess.site/punting.txt" punting <- read_table(my_url) ``` ## The data \small ```{r bRegression-47} punting ``` \normalsize ## Regression and output \footnotesize ```{r bRegression-48} punting.1 <- lm(punt ~ left + right + avg, data = punting) summary(punting.1) ``` \normalsize ## Comments - Overall regression strongly significant, R-sq high. - None of the $x$'s significant! Why? - $t$-tests only say that you could take any one of the $x$'s out without damaging the fit; doesn't matter which one. - Explanation: look at *correlations*. ## The correlations ```{r bRegression-49 } cor(punting) ``` - *All* correlations are high: $x$'s with `punt` (good) and with each other (bad, at least confusing). - What to do? Probably do just as well to pick one variable, say `right` since kickers are right-footed. - If we take out just the least significant variable, the two explanatory variables that remain are still going to be correlated with each other. ## Just `right` ```{r bRegression-50 } punting.2 <- lm(punt ~ right, data = punting) anova(punting.2, punting.1) ``` No significant loss by dropping other two variables. ## Comparing R-squareds - All three $x$-variables: ```{r bRegression-51 } summary(punting.1)$r.squared ``` - Only `right`: ```{r} summary(punting.2)$r.squared ``` - Basically no difference. In regression (over), `right` significant: ## Regression results \footnotesize ```{r bRegression-52 } summary(punting.2) ``` \normalsize ## But... - Maybe we got the *form* of the relationship with `left` wrong. - Check: plot *residuals* from previous regression (without `left`) against `left`. - Residuals here are "punting distance adjusted for right leg strength". - If there is some kind of relationship with `left`, we should include in model. - Plot of residuals against original variable: `augment` from `broom`. ## Residuals against `left` ```{r basingstoke,fig.height=4.5} punting.2 %>% augment(punting) -> punting.2a ggplot(punting.2a, aes(x = left, y = .resid)) + geom_point() ``` ## Comments - There is a *curved* relationship with `left`. - We should add `left`-squared to the regression (and therefore put `left` back in when we do that): ```{r bRegression-54 } punting.3 <- lm(punt ~ left + I(left^2) + right, data = punting ) ``` ## Regression with `left-squared` \scriptsize ```{r bRegression-55} summary(punting.3) ``` ## Comments - This was definitely a good idea (R-squared has clearly increased). - We would never have seen it without plotting residuals from `punting.2` (without `left`) against `left`. - Negative slope for `left`-squared means that increased left-leg strength only increases punting distance up to a point: beyond that, it decreases again.