Case study: Alcohol metabolism

Packages for this section

library(tidyverse)
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 (categorical)
Alcohol: whether subject was an alcoholic or not (assessed by how much alcohol they drank and how frequently). Values Alcoholic, Non-alcoholic (categorical).

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

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

A graph

Have two quantitative variables (suggests scatterplot), but also two categorical ones (suggests colours, but also something else)
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

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

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

Interpretation of Estimates

There is one intercept, and each explanatory variable has a slope that expresses effect of that variable, all else equal.
Slope for Gastric (quantitative) says that if Gastric increases by 1, and everything else same, Metabol predicted to increase by 1.95.
Each categorical variable has a “baseline” category, the first one alphabetically: Female, Alcoholic.
For categorical explanatory variables, the slope says how Metabol compares for the named category vs. the baseline:
- for males, Metabol is 1.65 higher than for females, all else equal
- for non-alcoholics, Metabol is 0.12 higher than for alcoholics, all else equal
Intercept is the predicted value of Metabol when all the quantitative variables are 0 and all the categorical variables are at baseline (not usually very interesting).

Create dataframe for plots

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,…
$ Metabol    <dbl> 0.6, 0.6, 1.5, 0.4, 0.1, 0.2, 0.3,…
$ Gastric    <dbl> 1.0, 1.6, 1.5, 2.2, 1.1, 1.2, 0.9,…
$ Sex        <chr> "Female", "Female", "Female", "Fem…
$ Alcohol    <chr> "Alcoholic", "Alcoholic", "Alcohol…
$ .fitted    <dbl> -0.06926514, 1.09870918, 0.9040467…
$ .resid     <dbl> 0.66926514, -0.49870918, 0.5959532…
$ .hat       <dbl> 0.17712194, 0.19504634, 0.18979002…
$ .sigma     <dbl> 1.371461, 1.374635, 1.372891, 1.32…
$ .cooksd    <dbl> 1.597846e-02, 1.021005e-02, 1.4003…
$ .std.resid <dbl> 0.544915310, -0.410544719, 0.48900…

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 next?

curves plus unequal spreads suggest transformation of response
problems with only one or two explanatory variables suggest transforming just those (our problems here are bigger than this)
if we have some theory about relationship, can use that to guide transformation (eg. take logs, but don’t have that 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 (link):

Box-Cox

Box and Cox developed a method for choosing a transformation from the ladder of powers.
- They worked out how to estimate the transformation as a power, but making zero be log (since power zero makes no sense otherwise)
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.3 <- lm(sqrt(Metabol) ~ Gastric + Sex, data = alc)
summary(alc.3)


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

Comments; set up to check residuals again

Alcohol was not significant in the regression, so removed it. If you have more than one non-significant explanatory variable, remove only the one least significant one (highest P-value), refit, re-evaluate.
Everything else significantly adds to the regression, so cannot remove anything else.

alc.3 %>% augment(alc) -> alc.3a

Expect the residual plots to be a lot better.
In principle, check again:
- residuals vs. fitted
- normal quantile plot of residuals
- residuals against each explanatory

Residuals against fitted

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

Normal quantile plot of residuals

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

Residuals against `Sex`

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

Comments

Fitted vs residual looks less curved
Residuals look much more normal (high values brought down)
Residuals vs Sex have no outliers and look equally spread

Back to regression output

glance(alc.3)

tidy(alc.3)

Comments

R-squared is reasonably high
Slope for Gastric is significantly positive: as gastric activity increases, metabolism increases (for everybody)
Slope for SexMale positive: compared to baseline Sex (Female), metabolism is higher for Males even allowing for the effect of Gastric.
We removed Alcohol: being alcoholic does not affect metabolism over and above the other explanatory variables.
This all seems to support the researchers’ theory.

More details of researchers’ theory

from Ramsey and Schafer, “The Statistical Sleuth”, 3rd edition, case 11.1.1, page 311:

Women exhibit a lower tolerance for alcohol and develop alcohol-related liver disease more readily than men. When men and women of the same size and drinking history consume equal amounts of alcohol, the women on average carry a higher concentration of alcohol in their bloodstream. According to a team of Italian researchers, this occurs because alcohol-degrading enzymes in the stomach (where alcohol is partially metabolized before it enters the bloodstream and is eventually metabolized by the liver) are more active in men than in women. The researchers studied the extent to which the activity of the enzyme explained the first-pass alcohol metabolism and the extent to which it explained the diﬀerences in first-pass metabolism between women and men.

Study details (also from Ramsey and Schafer)

The subjects were 18 women and 14 men, all volunteers living in Trieste.
Three of the women and five of the men were categorized as alcoholic.
All subjects received ethanol, at a dose of 0.3 grams per kilogram of body weight, orally one day and intravenously another, in randomly determined order.
Since the intravenous administration bypasses the stomach, the diﬀerence in blood alcohol concentration — the concentration after intravenous administration minus the concentration after oral administration — provides a measure of the “first-pass metabolism” in the stomach.
In addition, gastric alcohol dehydrogenase (AD) activity (activity of the key enzyme) was measured in mucus samples taken from the stomach linings.

Case study: Alcohol metabolism

Packages for this section

Metabolizing alcohol

The variables

Comments

The data

A graph

The graph

Comments

Fit regression

Comments

Interpretation of Estimates

Create dataframe for plots

Residuals vs fitted

Normal quantile plot

Residuals vs Gastric

Residuals vs Sex

Residuals vs Alcohol

Comments

What to do next?

Ladder of powers

Box-Cox

Running Box-Cox

Comments

Re-do regression

Output

Remove Alcohol:

Comments; set up to check residuals again

Residuals against fitted

Normal quantile plot of residuals

Residuals against Sex

Comments

Back to regression output

Comments

More details of researchers’ theory

Study details (also from Ramsey and Schafer)

Residuals vs `Gastric`

Residuals vs `Sex`

Residuals vs `Alcohol`

Remove `Alcohol`:

Residuals against `Sex`