Regression revisited

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

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:
my_url <- "http://datafiles.ritsokiguess.site/sleep.txt"
sleep <- read_delim(my_url, " ")

Check data

sleep

Make scatter plot of ATST (response) vs. age (explanatory).

The scatterplot

ggplot(sleep, aes(x = age, y = atst)) + geom_point()

Plot with smooth trend

The regression

  • Scatterplot shows no obvious curve, and a pretty clear downward trend. So we can run the regression:
sleep.1 <- lm(atst ~ age, data = sleep)
  • and look at the output via
summary(sleep.1)

(over)

The output


Call:
lm(formula = atst ~ age, data = sleep)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.011  -9.365   2.372   6.770  20.411 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  646.483     12.918   50.05 2.49e-14 ***
age          -14.041      1.368  -10.26 5.70e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.15 on 11 degrees of freedom
Multiple R-squared:  0.9054,    Adjusted R-squared:  0.8968 
F-statistic: 105.3 on 1 and 11 DF,  p-value: 5.7e-07

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):

tidy(sleep.1)

also one-line summary of model

glance(sleep.1)

augment

Useful for plotting residuals:

sleep.1 %>% augment(sleep) -> sleep.1a
sleep.1a

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

sleep.1a %>% 
  select(atst, age, .hat) %>% 
  mutate(cutoff = 2 * 2 / 13) %>% 
  slice_max(.hat, n = 5)
  • 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

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:

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:

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:

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\):

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.

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

my_url <- "http://datafiles.ritsokiguess.site/ice_crystals.csv"
ice_crystal <- read_csv(my_url)
ice_crystal

A graph

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

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:
boxcox(m ~ log(t), data = ice_crystal)

Regression with transformed response

ice.2 <- lm(log(m) ~ log(t), data = ice_crystal)
summary(ice.2)

Output on next slide.

Regression output


Call:
lm(formula = log(m) ~ log(t), data = ice_crystal)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6662 -0.2347 -0.0412  0.2735  0.4703 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -5.7278     0.7381   -7.76 1.42e-09 ***
log(t)        2.0315     0.1558   13.04 3.59e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3091 on 41 degrees of freedom
Multiple R-squared:  0.8056,    Adjusted R-squared:  0.8009 
F-statistic: 169.9 on 1 and 41 DF,  p-value: 3.589e-16

Are the residuals better now?

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:
my_url <- "http://datafiles.ritsokiguess.site/punting.txt"
punting <- read_table(my_url)

The data

punting

Regression and output

punting.1 <- lm(punt ~ left + right + avg, data = punting)
summary(punting.1)

Call:
lm(formula = punt ~ left + right + avg, data = punting)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.9325 -11.5618  -0.0315   9.0415  20.0886 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -4.6855    29.1172  -0.161    0.876
left          0.2679     2.1111   0.127    0.902
right         1.0524     2.1477   0.490    0.636
avg          -0.2672     4.2266  -0.063    0.951

Residual standard error: 14.68 on 9 degrees of freedom
Multiple R-squared:  0.7781,    Adjusted R-squared:  0.7042 
F-statistic: 10.52 on 3 and 9 DF,  p-value: 0.00267

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

cor(punting)
           left     right      punt       avg
left  1.0000000 0.8957224 0.8117368 0.9722632
right 0.8957224 1.0000000 0.8805469 0.9728784
punt  0.8117368 0.8805469 1.0000000 0.8679507
avg   0.9722632 0.9728784 0.8679507 1.0000000

  • 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

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:
summary(punting.1)$r.squared
[1] 0.7781401
  • Only right:
summary(punting.2)$r.squared
[1] 0.7753629
  • Basically no difference. In regression (over), right significant:

Regression results

summary(punting.2)

Call:
lm(formula = punt ~ right, data = punting)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.7576 -11.0611   0.3656   7.8890  19.0423 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.6930    25.2649  -0.146    0.886    
right         1.0427     0.1692   6.162 7.09e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.36 on 11 degrees of freedom
Multiple R-squared:  0.7754,    Adjusted R-squared:  0.7549 
F-statistic: 37.97 on 1 and 11 DF,  p-value: 7.088e-05

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

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):

punting.3 <- lm(punt ~ left + I(left^2) + right,
  data = punting
)

Regression with left-squared

summary(punting.3)

Call:
lm(formula = punt ~ left + I(left^2) + right, data = punting)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.3777  -5.3599   0.0459   4.5088  13.2669 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -4.623e+02  9.902e+01  -4.669  0.00117 **
left         6.888e+00  1.462e+00   4.710  0.00110 **
I(left^2)   -2.302e-02  4.927e-03  -4.672  0.00117 **
right        7.396e-01  2.292e-01   3.227  0.01038 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.931 on 9 degrees of freedom
Multiple R-squared:  0.9352,    Adjusted R-squared:  0.9136 
F-statistic:  43.3 on 3 and 9 DF,  p-value: 1.13e-05

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.