Exercise 6.1: Linear regression on the fish tank dataset - solution

1 Fish tank dataset

In this experiment, 96 fish (dojofish, goldfish and zebrafish) were placed separately in a tank with two liters of water and a certain dose (in mg) of the poison EI-43,064. The resistance of the fish against the poison was measured as the amount of minutes the fish survived after being exposed to the poison (Surv_time, in minutes). Additionally, the weight of each fish was measured.

2 Goal

The research goal is to study the association between the dose of the poison that was administered to the fish and their survival time.

Load the required libraries

library(tidyverse)

3 Import the data

poison <- read_csv("https://raw.githubusercontent.com/statOmics/PSLSData/main/poison.csv")

glimpse(poison)

## Rows: 96
## Columns: 4
## $ species   <dbl> 0, 2, 1, 2, 0, 0, 2, 1, 0, 1, 0, 1, 0, 0, 2, 0, 0, 0, 2, 0, …
## $ Weight    <dbl> 1.884359, 1.732261, 2.830895, 1.748562, 2.114166, 1.854602, …
## $ Dose      <dbl> 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, …
## $ Surv_time <dbl> 3.461493, 2.110417, 11.419131, 1.821191, 3.277619, 2.955556,…

head(poison)

4 Data tidying

We can see a couple of things in the data that can be improved:

1. Capitalise the fist column name. Hint: use the rename function.

2. Set the Species column as a factor. Hint: use the mutate and as.factor functions.

3. Change the species factor levels from 0, 1 and 2 to Dojofish, Goldfish and Zebrafish. Hint: use the fct_recode function inside the mutate function.

poison <- poison %>%
  rename("Species" = "species") %>%
  mutate(Species = as.factor(Species)) %>%
  mutate(
    Species = fct_recode(Species,
      Dojofish = "0", Goldfish = "1", Zebrafish = "2"
    )
  )

poison

5 Data exploration

How many fish do we have per species?

poison %>%
  count(Species)

Make a suitable visualization of the association between the dose and the survival time. Additionally, add the fish species as a color to the plot.

poison %>%
  ggplot(aes(x = Dose, y = Surv_time)) +
  geom_point() +
  stat_smooth(method = "loess") +
  geom_smooth(method = "lm", col = "black") +
  ylab("Survival time (min)") +
  xlab("Dose (mg)") +
  geom_point(aes(col = Species)) + ## not necessary, additional layer of information
  scale_color_manual(values = c("red", "darkgoldenrod", "black")) +
  theme_bw()

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

The linear regression line (black) is a good approximation of the best fitting smooth line (blue) through the data. Based on this figure, it seems realistic to suggest a linear relationship between dose and survival, where higher doses have lower survival times (as expected). However, before we start with the regression analysis we must check if all required assumptions are met.

6 Important note on the dataset

In this dataset, there are multiple variables that can have an effect on the survival time of the fish. The most obvious one is the dose of poison that was administered (as displayed above). However, we could also imagine that heavier fish are less prone to the poison than light fish. Additionally, one fish species may be more resistant to the poison than the other.

As such, not taking into account species into the analysis will introduce problems regarding the assumption of independence of the data. The assumption of independence entails that knowing the response value of one observation (fish) does not learn us anything about the response values for another variable. This will not be the case here; the response values (survival times) of fish of the same species will be much more alike compared to those of fish from another species.

Thus, in order to correctly analyze this data, fish species and weight should be taken into consideration. For now, we will “avoid” taking into account the effect of species into the analysis by only performing an analysis for the dojofish. As such, we will filter the dataset so that it only contains observations for the Dojofish.

# filter the data so it only contains dojofish
poison_dojo <- poison %>%
  filter(Species == "Dojofish")

Later in the course, we will come back to this dataset and perform a analysis that allows us to assess the different research hypotheses for all fish species at once!

The second predictor variable, Weight, still may influence our analysis; fish with higher weight are expected to be more resistant to the poison. If there would be a systematic bias in the experiment, e.g. if all lighter fish obtained a low dose of poison and heavier fish obtained a high dose of poison (so, if the dose was not correctly randomized accros the different weights), then we would not be able to correctly estimate the linear association between Dose and survival; indeed, this relationship would be confounded by the effect of weight.

poison_dojo %>%
  ggplot(aes(x = Dose, y = Weight)) +
  geom_point() +
  ggtitle("Association between dose and weight") +
  theme_bw() +
  stat_summary(
    fun = mean, geom = "point",
    col = "black", size = 4, shape = 24,
    fill = "red"
  )

In our dataset, there seems to be no systematic bias, i.e., no clear trend in the association between Dose and Weight. However, we do notice that the fish weights are not perfectly uniformly distributed across the different doses; we don’t see a clear trend, but we do see some fluctuations. To check if fluctuations of the size that we observe here could happen by random change if we carefully randomise, we may use simulations where we simulate the allocation of the fish to the different doses.

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

set.seed(1351)

plots <- lapply(
  1:9,
  function(x) {
    poison_dojo %>%
      mutate(Dose = sample(Dose)) %>% # randomization
      ggplot(aes(x = Dose, y = Weight)) +
      geom_point() +
      theme_bw() +
      stat_summary(
        geom = "point",
        fun = "mean",
        col = "black",
        size = 2,
        shape = 24,
        fill = "red"
      )
  }
)
do.call(gridExtra::grid.arrange, plots)

In the above simulation, we randomly assigned the different fish to different poison doses. Even though we did this randomly, we still see that by random chance there can be fluctuations as large as the ones observed in our real dataset. As such, we may conclude that there is indeed no systematic bias, i.e., no clear trend in the association between Dose and Weight.

7 Linear regression analysis

7.1 Assumptions of linear regression

List assumptions:

1. The observations are independent of each other
2. Linearity between the response and predictor variable
3. The residues of the model must be normally distributed
4. Homoscedasticity of the data

The first assumption is met based on the experimental design. In the data exploration phase, we saw that there is no clear trend in the association between Dose and Weight. With respect to species, we have “avoided” taking into account the effect of species into the analysis by only performing an analysis for the dojofish.

We will check the other assumptions by first fitting the linear model and plotting the output. As such, we will get all the required diagnostic plots.

# fit a linear regression model with 'Surv_time' as response variable and
# 'Dose' as predictor variabele
model <- lm(Surv_time ~ Dose, data = poison_dojo)

## display the diagnostic plots of the model
plot(model)

We have four diagnostic plots:

2. Linearity with the Residuals vs fitted plot

• predictor of predictions \(\hat\beta_0+\hat\beta_1 x\) on \(X\)-axis
• residuals on \(Y\)-as \[e_i=y_i-\hat{g}(x_i)=y_i-\hat\beta_0-\hat\beta_1\times x_i,\]

If there would be a linear relationship in the data, the residuals are expected to be scattered at random around y=0 for the entire range of predicted values. This is clearly the case: the assumption of linearity is met.

3. Normal Q-Q

• QQ-plot of the residuals \(e_i\).

The residuals of the linear regression model should be normally distributed. Given the second diagnostic plot, this seems not to be the case. In stead, we observe a short left tail and a long right tail. The assumption of normal residues is not met.

4. Homoscedasticity

• Square-root of the absolute value of standardized residuals in function of the fitted values

To meet the third assumption of linear regression, the variance on the Square-root of the absolute value of standardized residuals must be similar over the entire range of fitted values. The smoother in the plot helps us with looking at this; it should be nicely horizontal over the entire range of fitted values. This seems not to be the case: for larger fitted values, the variance increase slightly. The assumption of homoscedasticity is not met.

Multiple assumptions are not met, so we may not continue by performing a linear regression.

We can try to remediate these deviations by log-transforming the data. Indeed, data with heavy right tails can often be normalized by applying this strategy.

Fit a (log-)linear model by log-transforming the response.Generate diagnostic plots and assess the assumptions.

log.model <- lm(log10(Surv_time) ~ Dose, data = poison_dojo)
plot(log.model)

Arguably, there might be some deviation from normality in the left tail of the distribution. However, when we would simulate data under the normality assumption, it seems that deviations of this size may be expected when normality is met:

set.seed(1406)
nobs <- nrow(poison_dojo)

data.frame(
  y = c(log.model$res, rnorm(nobs * 8, sd = sigma(log.model))),
  label = rep(c("orig", paste0("sim", 1:8)), each = nobs)
) %>%
  ggplot(aes(sample = y)) +
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~label)

1. The independence assumption is met.
2. Upon log-transformation, the linearity assumption is still met.
3. Upon log-transformation, the normality assumption is also met.
4. Upon log-transformation, the homoscedasticity assumption is also met.

Look at the output of the log-linear model:

summary(log.model)

## 
## Call:
## lm(formula = log10(Surv_time) ~ Dose, data = poison_dojo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2074 -0.1145 -0.0324  0.1039  0.3268 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.91722    0.10423   8.800 1.33e-10 ***
## Dose        -0.27282    0.06648  -4.104 0.000215 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1453 on 37 degrees of freedom
## Multiple R-squared:  0.3128, Adjusted R-squared:  0.2942 
## F-statistic: 16.84 on 1 and 37 DF,  p-value: 0.0002146

Compute the 95% confidence interval on the model parameters:

confint(log.model)

##                  2.5 %     97.5 %
## (Intercept)  0.7060342  1.1284104
## Dose        -0.4075255 -0.1381159

7.2 Interpretation on the log scale

Currently, all the outcomes should be interpreted on the log-scale. Indeed, since we are now modelling the log\(_{10}\) of survival time, we don’t have direct inference on survival time.

As such, we may interpret the output as follows:

• The slope:

Dojofish that are exposed to a higher dose of the poison will have a \(log_{10}\) survival time that is on average 0.27 lower per gram of poison that was administered additionally (95% CI [0.71, 0.41, 1.13, 0.14]). This decrease in survival time is significant on the 5% significance level (p<0.001).

Note, that the study is an experimental study so we can conclude that there is a causal effect of the poison dose and the survival time.

• The intercept:

The average of the log10 of survival time for fish that were given a dose of 0 mg is 0.92.

7.3 Interpretation on the original scale

The interpretation the log-scale is quite difficult. To ease the interpretation, we will backtransform the results to the original scale (time in minutes). This we can do by

10^(summary(log.model)$coefficients[, "Estimate"])

## (Intercept)        Dose 
##   8.2646094   0.5335551

and of their confidence intervals:

10^(confint(log.model))

##                 2.5 %     97.5 %
## (Intercept) 5.0819949 13.4403459
## Dose        0.3912681  0.7275855

Now, we can interpret the results in terms of the geometric mean:

• The slope:

Dojo fish that are exposed to a higher dose of the poison will have a geometric mean of the survival time that is factor 0.53 (95% CI [0.39, 0.73]) per gram of poison that they received more than fish exposed to a lower dose. This decrease is very significant on the 5% significance level (p<0.001).

• The intercept:

The geometric mean of the survival time for fish that were given a dose of 0 mg is 8.26 minutes.

In other words, our analysis seems to suggest that fish in uncontaminated water would only survive for about 8 minutes! This seems a little bit fishy…

Note, that this interpretation is not biologically relevant and is induced because we extrapolated considerably because all fish in the experiment were exposed to poison.

## Extract slope and intercept from the model
dojo_intercept <- coef(log.model)[1]
dojo_slope <- coef(log.model)[2]

ggplot(poison_dojo, aes(Dose, log10(Surv_time))) +
  geom_point() +
  geom_abline(
    slope = dojo_slope, intercept = dojo_intercept,
    color = "steelblue", size = 1.2
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

We can see that the dose value of 0 mg does not fall within the same range as our data. The model is only valid in the range of data that we use to build the model on!

However, in most applications we will never test for the intercept term.

8 Conclusion

There is a very significant effect of the poison on survival of Dojofish (p< 0.001). Dojofish that are exposed to a higher dose of the poison will have a survival time that decrease on average with a factor 1.87 per gram of poison that is added (95% CI [2.56, 1.37]).

In the tutorial of multiple regression, we will revisit this exercise. More specifically, we will there study the association between dose and survival for all the species in the dataset.

---

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