Case study: asphalt

Packages for this section

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(broom)
library(MASS, exclude = "select") # explained later

Metabolizing alcohol

It is believed that women have a lower tolerance for alcohol compared to men, and they develop alcohol-related liver disease more easily than men (on average).
Italian researchers developed a theory about why this is. To test their theory, they collected some data (variables on next page).

The variables

Subject: subject number (ignored by us)
Metabol: first-pass metabolism of alcohol in stomach (millimoles per litre-hour), response
Gastric: gastric alcohol dehydrogenase activity in stomach (micromoles per minute per gram of tissue)
Sex: whether subject was Female or Male
Alcohol: whether subject was an alcoholic or not (assessed by how much alcohol they drank and how frequently)

Comments

one response variable, Metabol
several explanatory variables
- hence, multiple regression
- some of them are categorical.
according to theory, more gastric activity should mean higher metabolism
metabolism is expected to be higher for males than females even at same level of gastric activity
may be an effect of alcoholism, but direction not predicted by theory.

The data xxx

my_url <- "http://datafiles.ritsokiguess.site/alcohol-metabolism.csv"
alc <- read_csv(my_url)
alc

# A tibble: 32 × 5
   Subject Metabol Gastric Sex    Alcohol      
     <dbl>   <dbl>   <dbl> <chr>  <chr>        
 1       1     0.6     1   Female Alcoholic    
 2       2     0.6     1.6 Female Alcoholic    
 3       3     1.5     1.5 Female Alcoholic    
 4       4     0.4     2.2 Female Non-alcoholic
 5       5     0.1     1.1 Female Non-alcoholic
 6       6     0.2     1.2 Female Non-alcoholic
 7       7     0.3     0.9 Female Non-alcoholic
 8       8     0.3     0.8 Female Non-alcoholic
 9       9     0.4     1.5 Female Non-alcoholic
10      10     1       0.9 Female Non-alcoholic
# ℹ 22 more rows

A graph xxx

Have two quantitative variables (suggests scatterplot), but also two categorical ones (suggests colours, and ???)
Can also use shape as well as colour to distinguish observations.
Hence:

ggplot(alc, aes(x = Gastric, y = Metabol, colour = Sex, 
                shape = Alcohol)) + geom_point()

The graph xxx

Comments

As gastric activity increases, metabolism also increases.
Trend is apparently linear and of moderate strength.
Two points unusually high on both variables, but seem to be on the trend.
Male points (blue) mostly above female ones (red) for similar Gastric
No apparent effect of Alcohol.

Fit regression xxx

alc.1 <- lm(Metabol ~ Gastric + Sex + Alcohol, data = alc)
summary(alc.1)


Call:
lm(formula = Metabol ~ Gastric + Sex + Alcohol, data = alc)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.3119 -0.6339 -0.0927  0.6070  4.5629 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -2.0159     0.6348  -3.175  0.00362 ** 
Gastric                1.9466     0.2884   6.749  2.5e-07 ***
SexMale                1.6535     0.5514   2.999  0.00564 ** 
AlcoholNon-alcoholic   0.1183     0.6008   0.197  0.84527    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.354 on 28 degrees of freedom
Multiple R-squared:  0.7657,    Adjusted R-squared:  0.7406 
F-statistic: 30.51 on 3 and 28 DF,  p-value: 5.748e-09

Comments

no effect of Alcohol (suggests removing)
significant effects of both Gastric and Sex
both slopes positive (as expected)
but, before we go further, check residuals:
- vs. fitted values
- normal quantile plot
- also, vs. each explanatory variable (if any trends, have form of relationship with that explanatory variable wrong).
talk about how categoricals dealt with. xxx

Create dataframe for plots xxx

alc.1 %>% augment(alc) -> alc.1a
glimpse(alc.1a)

Rows: 32
Columns: 11
$ Subject    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, …
$ Metabol    <dbl> 0.6, 0.6, 1.5, 0.4, 0.1, 0.2, 0.3, 0.3,…
$ Gastric    <dbl> 1.0, 1.6, 1.5, 2.2, 1.1, 1.2, 0.9, 0.8,…
$ Sex        <chr> "Female", "Female", "Female", "Female",…
$ Alcohol    <chr> "Alcoholic", "Alcoholic", "Alcoholic", …
$ .fitted    <dbl> -0.06926514, 1.09870918, 0.90404679, 2.…
$ .resid     <dbl> 0.66926514, -0.49870918, 0.59595321, -1…
$ .hat       <dbl> 0.17712194, 0.19504634, 0.18979002, 0.0…
$ .sigma     <dbl> 1.371461, 1.374635, 1.372891, 1.320444,…
$ .cooksd    <dbl> 1.597846e-02, 1.021005e-02, 1.400360e-0…
$ .std.resid <dbl> 0.544915310, -0.410544719, 0.489003445,…

Residuals vs fitted

ggplot(alc.1a, aes(x = .fitted, y = .resid)) + geom_point()

Normal quantile plot

ggplot(alc.1a, aes(sample = .resid)) +
  stat_qq() + stat_qq_line()

Residuals vs `Gastric`

ggplot(alc.1a, aes(x = Gastric, y = .resid)) + geom_point()

Residuals vs `Sex`

ggplot(alc.1a, aes(x = Sex, y = .resid)) + geom_point()

Residuals vs `Alcohol`

ggplot(alc.1a, aes(x = Alcohol, y = .resid)) + geom_point()

Comments

residuals vs fitted and vs Gastric: suggestion of curve?
residuals vs Sex: males more spread out
residuals vs Alcohol: non-alcoholics more spread out

What to do?

curves plus unequal spreads suggest transformation of response
problems with only one or two explanatory variables suggest transforming just those (seem to be bigger problems here)
if we have some theory about relationship, can use that to guide transformation (eg. take logs, but not here)
can use Box-Cox to suggest good transformation.
- aims to find transformation of response that promotes straightness plus equal spread (will also often reduce influence of large values).

Ladder of powers

from here:

Box-Cox

Box and Cox developed a method for choosing a transformation from the ladder of powers.
In package MASS as boxcox
Technical detail: MASS also has a select that we don’t want to have interfere with the tidyverse select, so load MASS as shown:

library(MASS, exclude = "select")

Running Box-Cox

use boxcox with a model formula (such as you would use in lm):

boxcox(Metabol ~ Gastric + Sex + Alcohol, data = alc)

Comments

Output is a graph.
Peak of curve is single best power transformation
Outer vertical lines mark 95% CI for transformation power
Goal: find value from ladder of powers that is within CI
Here that is 0.5, square root
Note that 1, “do nothing”, and 0, “take logs”, are not supported by data.

Re-do regression

with square root of Metabol as response:

alc.2 <- lm(sqrt(Metabol) ~ Gastric + Sex + Alcohol, data = alc)

Output

summary(alc.2)


Call:
lm(formula = sqrt(Metabol) ~ Gastric + Sex + Alcohol, data = alc)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.74769 -0.26586  0.01346  0.25088  0.75600 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           0.22929    0.17494   1.311 0.200607    
Gastric               0.50406    0.07948   6.342 7.33e-07 ***
SexMale               0.55856    0.15195   3.676 0.000995 ***
AlcoholNon-alcoholic -0.04690    0.16556  -0.283 0.779042    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3731 on 28 degrees of freedom
Multiple R-squared:  0.768, Adjusted R-squared:  0.7431 
F-statistic: 30.89 on 3 and 28 DF,  p-value: 5.03e-09

Remove `Alcohol`:

alc.4 <- lm(sqrt(Metabol) ~ Gastric + Sex, data = alc)
summary(alc.4)


Call:
lm(formula = sqrt(Metabol) ~ Gastric + Sex, data = alc)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.77319 -0.23494  0.02865  0.26615  0.74255 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.20185    0.14334   1.408 0.169709    
Gastric      0.49655    0.07373   6.735 2.17e-07 ***
SexMale      0.57286    0.14103   4.062 0.000338 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3671 on 29 degrees of freedom
Multiple R-squared:  0.7673,    Adjusted R-squared:  0.7513 
F-statistic: 47.81 on 2 and 29 DF,  p-value: 6.584e-10

Check residuals again

Comments

The asphalt data

31 asphalt pavements prepared under different conditions. How does quality of pavement depend on these?
Variables:
- pct.a.surf Percentage of asphalt in surface layer
- pct.a.base Percentage of asphalt in base layer
- fines Percentage of fines in surface layer
- voids Percentage of voids in surface layer
- rut.depth Change in rut depth per million vehicle passes
- viscosity Viscosity of asphalt
- run 2 data collection periods: 1 for run 1, 0 for run 2.
rut.depth response. Depends on other variables, how?

Packages for this section

library(MASS, exclude = "select")
library(tidyverse)
library(broom)
library(leaps)

Make sure to load MASS before tidyverse (for annoying technical reasons), or to load MASS excluding its select (as above).

Getting set up

my_url <- "http://datafiles.ritsokiguess.site/asphalt.txt"
asphalt <- read_delim(my_url, " ")

Rows: 31 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: " "
dbl (7): pct.a.surf, pct.a.base, fines, voids, rut.depth, viscosity, run

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Quantitative variables with one response: multiple regression.
Some issues here that don’t come up in “simple” regression; handle as we go. (STAB27/STAC67 ideas.)

The data (some)

asphalt

# A tibble: 31 × 7
   pct.a.surf pct.a.base fines voids rut.depth viscosity   run
        <dbl>      <dbl> <dbl> <dbl>     <dbl>     <dbl> <dbl>
 1       4.68       4.87   8.4  4.92      6.75      2.8      1
 2       5.19       4.5    6.5  4.56     13         1.4      1
 3       4.82       4.73   7.9  5.32     14.8       1.4      1
 4       4.85       4.76   8.3  4.86     12.6       3.3      1
 5       4.86       4.95   8.4  3.78      8.25      1.7      1
 6       5.16       4.45   7.4  4.40     10.7       2.9      1
 7       4.82       5.05   6.8  4.87      7.28      3.7      1
 8       4.86       4.7    8.6  4.83     12.7       1.7      1
 9       4.78       4.84   6.7  4.86     12.6       0.92     1
10       5.16       4.76   7.7  4.03     20.6       0.68     1
# ℹ 21 more rows

Plotting response “rut depth” against everything else

Same idea as for plotting separate predictions on one plot:

asphalt %>%
  pivot_longer(
    -rut.depth,
    names_to="xname", values_to="x"
  ) %>%
  ggplot(aes(x = x, y = rut.depth)) + geom_point() +
  facet_wrap(~xname, scales = "free") -> g

“collect all the x-variables together into one column called x, with another column xname saying which x they were, then plot these x’s against rut.depth, a separate facet for each x-variable.”

I saved this graph to plot later (on the next page).

The plot

Interpreting the plots

One plot of rut depth against each of the six other variables.
Get rough idea of what’s going on.
Trends mostly weak.
viscosity has strong but non-linear trend.
run has effect but variability bigger when run is 1.
Weak but downward trend for voids.
Non-linearity of rut.depth-viscosity relationship should concern us.

Log of `viscosity`: more nearly linear?

Take this back to asphalt engineer: suggests log of viscosity:

ggplot(asphalt, aes(y = rut.depth, x = log(viscosity))) +
  geom_point() + geom_smooth(se = FALSE) -> g

(plot overleaf)

Rut depth against log-viscosity

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Comments and next steps

Not very linear, but better than before.
In multiple regression, hard to guess which x’s affect response. So typically start by predicting from everything else.
Model formula has response on left, squiggle, explanatories on right joined by plusses:

rut.1 <- lm(rut.depth ~ pct.a.surf + pct.a.base + fines +
  voids + log(viscosity) + run, data = asphalt)

Regression output:

summary(rut.1)


Call:
lm(formula = rut.depth ~ pct.a.surf + pct.a.base + fines + voids + 
    log(viscosity) + run, data = asphalt)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1211 -1.9075 -0.7175  1.6382  9.5947 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)    -12.9937    26.2188  -0.496   0.6247   
pct.a.surf       3.9706     2.4966   1.590   0.1248   
pct.a.base       1.2631     3.9703   0.318   0.7531   
fines            0.1164     1.0124   0.115   0.9094   
voids            0.5893     1.3244   0.445   0.6604   
log(viscosity)  -3.1515     0.9194  -3.428   0.0022 **
run             -1.9655     3.6472  -0.539   0.5949   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.324 on 24 degrees of freedom
Multiple R-squared:  0.806, Adjusted R-squared:  0.7575 
F-statistic: 16.62 on 6 and 24 DF,  p-value: 1.743e-07

Comments

R-squared 81%, not so bad.
P-value in glance asserts that something helping to predict rut.depth.
Table of coefficients says log(viscosity).
But confused by clearly non-significant variables: remove those to get clearer picture of what is helpful.

Before we do anything, look at residual plots:

1. of residuals against fitted values (as usual)
1. of residuals against each explanatory.
Problem fixes:
- with (a): fix response variable;
- with some plots in (b): fix those explanatory variables.

Plot fitted values against residuals

ggplot(rut.1, aes(x = .fitted, y = .resid)) + geom_point()

Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
ℹ Please use `broom::augment(<lm>)` instead.
ℹ The deprecated feature was likely used in the ggplot2 package.
  Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.

Normal quantile plot of residuals

ggplot(rut.1, aes(sample = .resid)) + stat_qq() + 
  stat_qq_line()

Plotting residuals against \(x\) variables

Problem here is that residuals are in the fitted model, and the observed \(x\)-values are in the original data frame asphalt.
Package broom contains a function augment that combines these two together so that they can later be plotted: start with a model first, and then augment with a data frame:

rut.1 %>% augment(asphalt) -> rut.1a
rut.1a

# A tibble: 31 × 13
   pct.a.surf pct.a.base fines voids rut.depth viscosity   run .fitted .resid
        <dbl>      <dbl> <dbl> <dbl>     <dbl>     <dbl> <dbl>   <dbl>  <dbl>
 1       4.68       4.87   8.4  4.92      6.75      2.8      1    10.4 -3.65 
 2       5.19       4.5    6.5  4.56     13         1.4      1    13.7 -0.718
 3       4.82       4.73   7.9  5.32     14.8       1.4      1    13.1  1.60 
 4       4.85       4.76   8.3  4.86     12.6       3.3      1    10.4  2.22 
 5       4.86       4.95   8.4  3.78      8.25      1.7      1    12.1 -3.87 
 6       5.16       4.45   7.4  4.40     10.7       2.9      1    11.2 -0.577
 7       4.82       5.05   6.8  4.87      7.28      3.7      1    10.1 -2.81 
 8       4.86       4.7    8.6  4.83     12.7       1.7      1    12.4  0.221
 9       4.78       4.84   6.7  4.86     12.6       0.92     1    14.0 -1.46 
10       5.16       4.76   7.7  4.03     20.6       0.68     1    16.0  4.57 
# ℹ 21 more rows
# ℹ 4 more variables: .hat <dbl>, .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>

What does rut.1a contain?

names(rut.1a)

 [1] "pct.a.surf" "pct.a.base" "fines"      "voids"      "rut.depth" 
 [6] "viscosity"  "run"        ".fitted"    ".resid"     ".hat"      
[11] ".sigma"     ".cooksd"    ".std.resid"

all the stuff in original data frame, plus:
quantities from regression (starting with a dot)

Plotting residuals against \(x\)-variables

rut.1a %>%
  mutate(log_vis=log(viscosity)) %>% 
  pivot_longer(
    c(pct.a.surf:voids, run, log_vis),
    names_to="xname", values_to="x"
  ) %>%
  ggplot(aes(x = x, y = .resid)) +
  geom_point() + facet_wrap(~xname, scales = "free") -> g

The plot

Comments

There is serious curve in plot of residuals vs. fitted values. Suggests a transformation of \(y\).
The residuals-vs-\(x\)’s plots don’t show any serious trends. Worst probably that potential curve against log-viscosity.
Also, large positive residual, 10, that shows up on all plots. Perhaps transformation of \(y\) will help with this too.
If residual-fitted plot OK, but some residual-\(x\) plots not, try transforming those \(x\)’s, eg. by adding \(x^2\) to help with curve.

Which transformation?

Best way: consult with person who brought you the data.
Can’t do that here!
No idea what transformation would be good.
Let data choose: “Box-Cox transformation”.
Scale is that of “ladder of powers”: power transformation, but 0 is log.

Running Box-Cox

From package MASS:

boxcox(rut.depth ~ pct.a.surf + pct.a.base + fines + voids +
  log(viscosity) + run, data = asphalt)

Comments on Box-Cox plot

\(\lambda\) represents power to transform \(y\) with.
Best single choice of transformation parameter \(\lambda\) is peak of curve, close to 0.
Vertical dotted lines give CI for \(\lambda\), about (−0.05, 0.2).
\(\lambda = 0\) means “log”.
Narrowness of confidence interval mean that these not supported by data:
- No transformation (\(\lambda = 1\))
- Square root (\(\lambda = 0.5\))
- Reciprocal (\(\lambda = −1\)).

Relationships with explanatories

As before: plot response (now log(rut.depth)) against other explanatory variables, all in one shot:

asphalt %>%
  mutate(log_vis=log(viscosity)) %>% 
  pivot_longer(
    c(pct.a.surf:voids, run, log_vis),
    names_to="xname", values_to="x"
  ) %>%
  ggplot(aes(y = log(rut.depth), x = x)) + geom_point() +
  facet_wrap(~xname, scales = "free") -> g3

The new plots

g3

Modelling with transformed response

These trends look pretty straight, especially with log.viscosity.
Values of log.rut.depth for each run have same spread.
Other trends weak, but are straight if they exist.
Start modelling from the beginning again.
Model log.rut.depth in terms of everything else, see what can be removed:

rut.2 <- lm(log(rut.depth) ~ pct.a.surf + pct.a.base +
  fines + voids + log(viscosity) + run, data = asphalt)

use tidy from broom to display just the coefficients.

Output

tidy(rut.2)

# A tibble: 7 × 5
  term           estimate std.error statistic     p.value
  <chr>             <dbl>     <dbl>     <dbl>       <dbl>
1 (Intercept)     -1.57      2.44      -0.646 0.525      
2 pct.a.surf       0.584     0.232      2.52  0.0190     
3 pct.a.base      -0.103     0.369     -0.280 0.782      
4 fines            0.0978    0.0941     1.04  0.309      
5 voids            0.199     0.123      1.62  0.119      
6 log(viscosity)  -0.558     0.0854    -6.53  0.000000945
7 run              0.340     0.339      1.00  0.326

Taking out everything non-significant

Try: remove everything but pct.a.surf and log.viscosity:

rut.3 <- lm(log(rut.depth) ~ pct.a.surf + log(viscosity), data = asphalt)
tidy(rut.3)

# A tibble: 3 × 5
  term           estimate std.error statistic  p.value
  <chr>             <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)       0.900    1.08       0.833 4.12e- 1
2 pct.a.surf        0.391    0.219      1.79  8.46e- 2
3 log(viscosity)   -0.619    0.0271   -22.8   1.27e-19

summary(rut.3)


Call:
lm(formula = log(rut.depth) ~ pct.a.surf + log(viscosity), data = asphalt)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.61938 -0.21361  0.06635  0.14932  0.63012 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.90014    1.08059   0.833   0.4119    
pct.a.surf      0.39115    0.21879   1.788   0.0846 .  
log(viscosity) -0.61856    0.02713 -22.797   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3208 on 28 degrees of freedom
Multiple R-squared:  0.9509,    Adjusted R-squared:  0.9474 
F-statistic: 270.9 on 2 and 28 DF,  p-value: < 2.2e-16

Check that removing all those variables wasn’t too much

anova(rut.3, rut.2)

Analysis of Variance Table

Model 1: log(rut.depth) ~ pct.a.surf + log(viscosity)
Model 2: log(rut.depth) ~ pct.a.surf + pct.a.base + fines + voids + log(viscosity) + 
    run
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     28 2.8809                           
2     24 2.2888  4   0.59216 1.5523 0.2191

\(H_0\) : two models equally good; \(H_a\) : bigger model better.
Null not rejected here; small model as good as the big one, so prefer simpler smaller model rut.3.

Find the largest P-value by eye:

tidy(rut.2)

# A tibble: 7 × 5
  term           estimate std.error statistic     p.value
  <chr>             <dbl>     <dbl>     <dbl>       <dbl>
1 (Intercept)     -1.57      2.44      -0.646 0.525      
2 pct.a.surf       0.584     0.232      2.52  0.0190     
3 pct.a.base      -0.103     0.369     -0.280 0.782      
4 fines            0.0978    0.0941     1.04  0.309      
5 voids            0.199     0.123      1.62  0.119      
6 log(viscosity)  -0.558     0.0854    -6.53  0.000000945
7 run              0.340     0.339      1.00  0.326

Largest P-value is 0.78 for pct.a.base, not significant.
So remove this first, re-fit and re-assess.
Or, as over.

Get the computer to find the largest P-value for you

Output from tidy is itself a data frame, thus:

tidy(rut.2) %>% arrange(p.value)

# A tibble: 7 × 5
  term           estimate std.error statistic     p.value
  <chr>             <dbl>     <dbl>     <dbl>       <dbl>
1 log(viscosity)  -0.558     0.0854    -6.53  0.000000945
2 pct.a.surf       0.584     0.232      2.52  0.0190     
3 voids            0.199     0.123      1.62  0.119      
4 fines            0.0978    0.0941     1.04  0.309      
5 run              0.340     0.339      1.00  0.326      
6 (Intercept)     -1.57      2.44      -0.646 0.525      
7 pct.a.base      -0.103     0.369     -0.280 0.782

Largest P-value at the bottom.

Take out `pct.a.base`

Copy and paste the lm code and remove what you’re removing:

rut.4 <- lm(log(rut.depth) ~ pct.a.surf + fines + voids + 
              log(viscosity) + run, data = asphalt)
summary(rut.4)


Call:
lm(formula = log(rut.depth) ~ pct.a.surf + fines + voids + log(viscosity) + 
    run, data = asphalt)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.51610 -0.18785 -0.02248  0.18364  0.57160 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -2.07850    1.60665  -1.294   0.2076    
pct.a.surf      0.59299    0.22526   2.632   0.0143 *  
fines           0.08895    0.08701   1.022   0.3165    
voids           0.20047    0.12064   1.662   0.1091    
log(viscosity) -0.55249    0.08184  -6.751 4.48e-07 ***
run             0.35977    0.32533   1.106   0.2793    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3031 on 25 degrees of freedom
Multiple R-squared:  0.9608,    Adjusted R-squared:  0.953 
F-statistic: 122.7 on 5 and 25 DF,  p-value: < 2.2e-16

tidy(rut.4) %>% arrange(p.value) %>% select(term, p.value)

# A tibble: 6 × 2
  term               p.value
  <chr>                <dbl>
1 log(viscosity) 0.000000448
2 pct.a.surf     0.0143     
3 voids          0.109      
4 (Intercept)    0.208      
5 run            0.279      
6 fines          0.316

fines is next to go, P-value 0.32.

“Update”

Another way to do the same thing:

rut.4 <- update(rut.2, . ~ . - pct.a.base)
tidy(rut.4) %>% arrange(p.value)

# A tibble: 6 × 5
  term           estimate std.error statistic     p.value
  <chr>             <dbl>     <dbl>     <dbl>       <dbl>
1 log(viscosity)  -0.552     0.0818     -6.75 0.000000448
2 pct.a.surf       0.593     0.225       2.63 0.0143     
3 voids            0.200     0.121       1.66 0.109      
4 (Intercept)     -2.08      1.61       -1.29 0.208      
5 run              0.360     0.325       1.11 0.279      
6 fines            0.0889    0.0870      1.02 0.316

Again, fines is the one to go. (Output identical as it should be.)

Take out fines:

rut.5 <- update(rut.4, . ~ . - fines)
summary(rut.5)


Call:
lm(formula = log(rut.depth) ~ pct.a.surf + voids + log(viscosity) + 
    run, data = asphalt)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.57275 -0.20080  0.01061  0.17711  0.59774 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -1.25533    1.39147  -0.902   0.3753    
pct.a.surf      0.54837    0.22118   2.479   0.0200 *  
voids           0.23188    0.11676   1.986   0.0577 .  
log(viscosity) -0.58039    0.07723  -7.516 5.59e-08 ***
run             0.29468    0.31931   0.923   0.3646    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3033 on 26 degrees of freedom
Multiple R-squared:  0.9592,    Adjusted R-squared:  0.9529 
F-statistic: 152.8 on 4 and 26 DF,  p-value: < 2.2e-16

tidy(rut.5) %>% arrange(p.value) %>% select(term, p.value)

# A tibble: 5 × 2
  term                p.value
  <chr>                 <dbl>
1 log(viscosity) 0.0000000559
2 pct.a.surf     0.0200      
3 voids          0.0577      
4 run            0.365       
5 (Intercept)    0.375

Can’t take out intercept, so run, with P-value 0.36, goes next.

Take out run:

rut.6 <- update(rut.5, . ~ . - run)
summary(rut.6)


Call:
lm(formula = log(rut.depth) ~ pct.a.surf + voids + log(viscosity), 
    data = asphalt)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53548 -0.20181 -0.01702  0.16748  0.54707 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -1.02079    1.36430  -0.748   0.4608    
pct.a.surf      0.55547    0.22044   2.520   0.0180 *  
voids           0.24479    0.11560   2.118   0.0436 *  
log(viscosity) -0.64649    0.02879 -22.458   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3025 on 27 degrees of freedom
Multiple R-squared:  0.9579,    Adjusted R-squared:  0.9532 
F-statistic: 204.6 on 3 and 27 DF,  p-value: < 2.2e-16

tidy(rut.6) %>% arrange(p.value) %>% select(term, p.value)

# A tibble: 4 × 2
  term            p.value
  <chr>             <dbl>
1 log(viscosity) 5.29e-19
2 pct.a.surf     1.80e- 2
3 voids          4.36e- 2
4 (Intercept)    4.61e- 1

Again, can’t take out intercept, so largest P-value is for voids, 0.044. But this is significant, so we shouldn’t remove voids.

Comments

Here we stop: pct.a.surf, voids and log.viscosity would all make fit significantly worse if removed. So they stay.
Different final result from taking things out one at a time (top), than by taking out 4 at once (bottom):

coef(rut.6)

   (Intercept)     pct.a.surf          voids log(viscosity) 
    -1.0207945      0.5554686      0.2447934     -0.6464911

coef(rut.3)

   (Intercept)     pct.a.surf log(viscosity) 
     0.9001389      0.3911481     -0.6185628

Point: Can make difference which way we go.

Comments on variable selection

Best way to decide which \(x\)’s belong: expert knowledge: which of them should be important.
Best automatic method: what we did, “backward selection”.
Do not learn about “stepwise regression”! eg. here
R has function step that does backward selection, like this:

step(rut.2, direction = "backward", test = "F")

Gets same answer as we did (by removing least significant x).

Removing non-significant \(x\)’s may remove interesting ones whose P-values happened not to reach 0.05. Consider using less stringent cutoff like 0.20 or even bigger.
Can also fit all possible regressions, as over (may need to do install.packages("leaps") first).

All possible regressions (output over)

Uses package leaps:

leaps <- regsubsets(log(rut.depth) ~ pct.a.surf + 
                      pct.a.base + fines + voids + 
                      log(viscosity) + run, 
                    data = asphalt, nbest = 2)
s <- summary(leaps)
with(s, data.frame(rsq, outmat)) -> d

The output

d %>% rownames_to_column("model") %>% arrange(desc(rsq))

      model       rsq pct.a.surf pct.a.base fines voids log.viscosity. run
1  6  ( 1 ) 0.9609642          *          *     *     *              *   *
2  5  ( 1 ) 0.9608365          *                *     *              *   *
3  5  ( 2 ) 0.9593265          *          *     *     *              *    
4  4  ( 1 ) 0.9591996          *                      *              *   *
5  4  ( 2 ) 0.9589206          *                *     *              *    
6  3  ( 1 ) 0.9578631          *                      *              *    
7  3  ( 2 ) 0.9534561          *                *                    *    
8  2  ( 1 ) 0.9508647          *                                     *    
9  2  ( 2 ) 0.9479541                                 *              *    
10 1  ( 1 ) 0.9452562                                                *    
11 1  ( 2 ) 0.8624107                                                    *

Comments

Problem: even adding a worthless x increases R-squared. So try for line where R-squared stops increasing “too much”, eg. top line (just log.viscosity), first 3-variable line (backwards-elimination model). Hard to judge.
One solution (STAC67): adjusted R-squared, where adding worthless variable makes it go down.
data.frame rather than tibble because there are several columns in outmat.

All possible regressions, adjusted R-squared

with(s, data.frame(adjr2, outmat)) %>% 
  rownames_to_column("model") %>% 
  arrange(desc(adjr2))

      model     adjr2 pct.a.surf pct.a.base fines voids log.viscosity. run
1  3  ( 1 ) 0.9531812          *                      *              *    
2  5  ( 1 ) 0.9530038          *                *     *              *   *
3  4  ( 1 ) 0.9529226          *                      *              *   *
4  4  ( 2 ) 0.9526007          *                *     *              *    
5  6  ( 1 ) 0.9512052          *          *     *     *              *   *
6  5  ( 2 ) 0.9511918          *          *     *     *              *    
7  3  ( 2 ) 0.9482845          *                *                    *    
8  2  ( 1 ) 0.9473550          *                                     *    
9  2  ( 2 ) 0.9442365                                 *              *    
10 1  ( 1 ) 0.9433685                                                *    
11 1  ( 2 ) 0.8576662                                                    *

Revisiting the best model

Best model was our rut.6:

tidy(rut.6)

# A tibble: 4 × 5
  term           estimate std.error statistic  p.value
  <chr>             <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)      -1.02     1.36      -0.748 4.61e- 1
2 pct.a.surf        0.555    0.220      2.52  1.80e- 2
3 voids             0.245    0.116      2.12  4.36e- 2
4 log(viscosity)   -0.646    0.0288   -22.5   5.29e-19

Revisiting (2)

Regression slopes say that rut depth increases as log-viscosity decreases, pct.a.surf increases and voids increases. This more or less checks out with out scatterplots against log.viscosity.
We should check residual plots again, though previous scatterplots say it’s unlikely that there will be a problem:

g <- ggplot(rut.6, aes(y = .resid, x = .fitted)) + 
geom_point()

Residuals against fitted values

Normal quantile plot of residuals

ggplot(rut.6, aes(sample = .resid)) + stat_qq() + stat_qq_line()

Plotting residuals against x’s

Do our trick again to put them all on one plot:

augment(rut.6)

# A tibble: 31 × 10
   `log(rut.depth)` pct.a.surf voids `log(viscosity)` .fitted  .resid   .hat
              <dbl>      <dbl> <dbl>            <dbl>   <dbl>   <dbl>  <dbl>
 1             1.91       4.68  4.92           1.03      2.12 -0.207  0.0772
 2             2.56       5.19  4.56           0.336     2.76 -0.197  0.111 
 3             2.69       4.82  5.32           0.336     2.74 -0.0503 0.142 
 4             2.53       4.85  4.86           1.19      2.09  0.441  0.0528
 5             2.11       4.86  3.78           0.531     2.26 -0.150  0.178 
 6             2.37       5.16  4.40           1.06      2.23  0.134  0.0839
 7             1.99       4.82  4.87           1.31      2.00 -0.0170 0.0510
 8             2.54       4.86  4.83           0.531     2.52  0.0216 0.0723
 9             2.53       4.78  4.86          -0.0834    2.88 -0.347  0.105 
10             3.03       5.16  4.03          -0.386     3.08 -0.0570 0.139 
# ℹ 21 more rows
# ℹ 3 more variables: .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>

augment(rut.6, data = asphalt) %>%
  mutate(log_vis=log(viscosity)) %>% 
  pivot_longer(
    c(pct.a.surf:voids, run, log_vis),
    names_to="xname", values_to="x",
  ) %>%
  ggplot(aes(y = .resid, x = x)) + geom_point() +
  facet_wrap(~xname, scales = "free") -> g2

Residuals against the x’s

g2

Comments

None of the plots show any sort of pattern. The points all look random on each plot.
On the plot of fitted values (and on the one of log.viscosity), the points seem to form a “left half” and a “right half” with a gap in the middle. This is not a concern.
One of the pct.a.surf values is low outlier (4), shows up top left of that plot.
Only two possible values of run; the points in each group look randomly scattered around 0, with equal spreads.
Residuals seem to go above zero further than below, suggesting a mild non-normality, but not enough to be a problem.

Variable-selection strategies

Expert knowledge.
Backward elimination.
All possible regressions.
Taking a variety of models to experts and asking their opinion.
Use a looser cutoff to eliminate variables in backward elimination (eg. only if P-value greater than 0.20).
If goal is prediction, eliminating worthless variables less important.
If goal is understanding, want to eliminate worthless variables where possible.
Results of variable selection not always reproducible, so caution advised.

Case study: asphalt

Packages for this section

Metabolizing alcohol

The variables

Comments

The data xxx

A graph xxx

The graph xxx

Comments

Fit regression xxx

Comments

Create dataframe for plots xxx

Residuals vs fitted

Normal quantile plot

Residuals vs Gastric

Residuals vs Sex

Residuals vs Alcohol

Comments

What to do?

Ladder of powers

Box-Cox

Running Box-Cox

Comments

Re-do regression

Output

Remove Alcohol:

Check residuals again

Comments

The asphalt data

Packages for this section

Getting set up

The data (some)

Plotting response “rut depth” against everything else

The plot

Interpreting the plots

Log of viscosity: more nearly linear?

Rut depth against log-viscosity

Comments and next steps

Regression output:

Comments

Before we do anything, look at residual plots:

Plot fitted values against residuals

Normal quantile plot of residuals

Plotting residuals against \(x\) variables

What does rut.1a contain?

Plotting residuals against \(x\)-variables

The plot

Comments

Which transformation?

Running Box-Cox

Comments on Box-Cox plot

Relationships with explanatories

The new plots

Modelling with transformed response

Output

Taking out everything non-significant

Check that removing all those variables wasn’t too much

Find the largest P-value by eye:

Get the computer to find the largest P-value for you

Take out pct.a.base

“Update”

Take out fines:

Take out run:

Comments

Comments on variable selection

All possible regressions (output over)

The output

Comments

All possible regressions, adjusted R-squared

Revisiting the best model

Revisiting (2)

Residuals against fitted values

Normal quantile plot of residuals

Plotting residuals against x’s

Residuals against the x’s

Comments

Variable-selection strategies

Residuals vs `Gastric`

Residuals vs `Sex`

Residuals vs `Alcohol`

Remove `Alcohol`:

Log of `viscosity`: more nearly linear?

Take out `pct.a.base`