Exercise 6.1: Linear regression on puromycin data

library(tidyverse)

1 Puromycin data

Data on the velocity of an enzymatic reaction were obtained by Treloar (1974). The number of counts per minute of radioactive product from the reaction was measured as a function of substrate concentration in parts per million (ppm) and from these counts the initial rate (or velocity) of the reaction was calculated (counts/min/min). The experiment was conducted once with the enzyme treated with Puromycin, and once with the enzyme untreated.

Assess if there is an association between the substrate concentration and rate for the treated enzyme.

data(Puromycin)
Puromycin <- Puromycin %>%
  filter(state == "treated")

1.1 Data Exploration

Puromycin %>%
  ggplot(aes(conc, rate)) +
  geom_point() +
  stat_smooth(method = "lm") +
  stat_smooth(color = "firebrick", se = FALSE) +
  labs(
    x = "Concentration (ppm)",
    y = "Reaction rate (counts/min/min)"
  ) +
  ggtitle("Reaction rates for puromycin-treated enzymes")

## `geom_smooth()` using formula 'y ~ x'

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

The plot shows that there is a relation between the velocity and the concentration, however, the relation does not seem to be linear.

We will assess the impact of log-transforming the concentration. Because the concentration is measured in ppm we will log\(_{10}\) transform the data.

Puromycin %>%
  ggplot(aes(log10(conc), rate)) +
  geom_point() +
  stat_smooth(method = "lm") +
  stat_smooth(color = "firebrick", se = FALSE) +
  labs(
    x = "Concentration (log10(ppm))",
    y = "Reaction rate (counts/min/min)"
  ) +
  ggtitle("Reaction rates for puromycin-treated enzymes")

## `geom_smooth()` using formula 'y ~ x'

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

The relation between the velocity and the log\(_{10}\) transformed concentration seems to be linear.

2 Linear regression

We will fit the following model to the data

\(Y_i = \beta_0 + \beta_1 log_{10}(x_i) + \epsilon_i\)

with \(Y_i\) the reaction rate, \(x_i\) the adopted concentration, \(\beta_0\) the intercept, \(\beta_1\) the slope or the effect of the log\(_{10}\)(concentratie) on the average reaction rate, and \(\epsilon_i\) the error term that is \(\epsilon_i\) i.i.d. normally distributed with mean 0 and variance \(\sigma^2\).

Note, that we write the substrate concentratio with a small letter because the predictor is not random. The researchers have chosen the substrate concentrations in the design phase and it is thus no random variable.

mod <- lm(rate ~ log10(conc), Puromycin)
summary(mod)

## 
## Call:
## lm(formula = rate ~ log10(conc), data = Puromycin)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.0176  -6.2455   0.6039   7.4262  13.3228 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  209.194      5.045   41.47 1.59e-12 ***
## log10(conc)   85.450      5.133   16.65 1.28e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.37 on 10 degrees of freedom
## Multiple R-squared:  0.9652, Adjusted R-squared:  0.9617 
## F-statistic: 277.1 on 1 and 10 DF,  p-value: 1.28e-08

Before we perform inference we will first assess the assumptions

2.1 Assumptions

1. Independence
2. Linearity
3. Normal distribution of the residuals
4. Homoscedasticity

We assume that the experiment was well designed and that the different reactions that were use in the experiment are independent.

2.2 Linearity

We assess linearity in a residual analysis

plot(mod, which = 1)

There seem to be small deviations for the residuals at high concentrations. However, we do not have many observation in the study so this might be caused by random varibility.

We will use simulation to assess if we can observe similar residual plots if all assumptions for the linear model hold.

set.seed(1031)
sigma <- sigma(mod)
dataHlp <- Puromycin

simModels <- list()

par(mfrow = c(3, 3))
for (i in 1:9) {
  nobs <- nrow(Puromycin)
  dataHlp$ySim <- mod$fit + rnorm(nobs, sd = sigma)
  simModels[[i]] <- lm(ySim ~ log10(conc), dataHlp)
  plot(simModels[[i]], which = 1)
}

We observe that similar deviations might be generated for a model where all assumptions are valid.

2.2.1 Normality

plot(mod, which = 2)

The QQ-plot does not show large deviations from normality.

2.2.2 Homoscedasticity: equality of the variance

We can again use the residual plot for assessing this assumption or the plot were we plot the square root of the standardized residuals in function of the fit.

plot(mod, which = 3)

We see that the spread of the majority of the residuals is more or less similar. However, the two residuals at low concentrations are slightly higher. Again we can use our simulations to show that similar deviations might occur for models where all assumptions hold.

par(mfrow = c(3, 3))
for (i in 1:9) {
  plot(simModels[[i]], which = 3)
}

2.3 Inference

2.3.1 Hypotheses

We want to assess if there is a linear association between the reaction rate and the log\(_{10}\) transformed concentration.

So we will assess

\(H_0\): \(\beta_1 = 0\)

against

\(H_A\): \(\beta_1 \ne 0\)

using a t-test on the slope parameter.

summary(mod)

## 
## Call:
## lm(formula = rate ~ log10(conc), data = Puromycin)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.0176  -6.2455   0.6039   7.4262  13.3228 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  209.194      5.045   41.47 1.59e-12 ***
## log10(conc)   85.450      5.133   16.65 1.28e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.37 on 10 degrees of freedom
## Multiple R-squared:  0.9652, Adjusted R-squared:  0.9617 
## F-statistic: 277.1 on 1 and 10 DF,  p-value: 1.28e-08

confint(mod)

##                 2.5 %    97.5 %
## (Intercept) 197.95335 220.43564
## log10(conc)  74.01255  96.88732

2.3.2 Conclusion

There is an extremely significant linear association between the substrate concentration on the log-scale and the reaction rate (p<<0.001). A reaction at a substrate concentration that is 10 times higher will have a reaction speed that is on average 85.4 counts/min higher (95% CI [74, 96.9] counts/min).

---

This work is licensed under the CC BY-NC-SA 4.0 licence.