Tutorial 6.2: Linear regression on the fish tank dataset

As a short exercise on linear regression, we will analyse the fish tank dataset.

1 The 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 poison EI-43,064. The resistance of the fish against the poison was measured as the amount of minutes the fish survived upon adding the poison (Surv_time, in minutes). Additionally, the weight of each fish was measured.

2 Goal

In this tutorial, we will study the association between dose and survival time by using a linear regression model.

Read the required libraries

library(tidyverse)

3 Import data

poison <- read_csv("poison.csv")

4 Data Exploration

head(poison)

## # A tibble: 6 x 4
##   species Weight  Dose Surv_time
##     <dbl>  <dbl> <dbl>     <dbl>
## 1       0   1.88   1        3.46
## 2       2   1.73   1.1      2.11
## 3       1   2.83   1.2     11.4 
## 4       2   1.75   1.3      1.82
## 5       0   2.11   1.4      3.28
## 6       0   1.85   1.5      2.96

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

1. Capitalize the fist column name
2. Set the Species column as a factor
3. Change the speciec factor levels from 0, 1 and 2 to Dojofish, Goldfish and Zebrafish. Hint: use the fct_recode function.

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

poison

## # A tibble: 96 x 4
##    Species   Weight  Dose Surv_time
##    <fct>      <dbl> <dbl>     <dbl>
##  1 Dojofish    1.88   1        3.46
##  2 Zebrafish   1.73   1.1      2.11
##  3 Goldfish    2.83   1.2     11.4 
##  4 Zebrafish   1.75   1.3      1.82
##  5 Dojofish    2.11   1.4      3.28
##  6 Dojofish    1.85   1.5      2.96
##  7 Zebrafish   1.22   1.6      1.85
##  8 Goldfish    2.05   1.7      8.05
##  9 Dojofish    1.80   1.8      1.65
## 10 Goldfish    2.50   1.9      7.95
## # … with 86 more rows

How many fish do we have per species?

poison %>%
  count(Species)

## # A tibble: 3 x 2
##   Species       n
##   <fct>     <int>
## 1 Dojofish     39
## 2 Goldfish     38
## 3 Zebrafish    19

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

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 analyis we must check if all required assumptions are met.

5 Analysis

5.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 not met. There are patterns of relatedness in the data. Indeed, we expect fish of the same species to have a more similar response to the poison as compared to fish of another species. Equivalently, we expect that fish of the same weight have a more similar response to the poison as compared to fish of different weight.

** For now, we will ignore this problem . In a later tutorial, we will see how we can account for this dependence by including the species data and weigth data in the linear model (class on multiple regression). As such, you may assume the assumption is met. **

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) 

## 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 lie on the y=0 line 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 way not continue by performing a linear regression.

One “trick” is to log-transform the data. Indeed, data with heavy right tails can often be normalized by applying this strategy

Generate a new column in the poison dataset, containing the log-transformed survival times.

poison <- poison %>%
  mutate(log.Surv_time = log(Surv_time))

Refit the (log-)linear model. Generate diagnostic plots. Assess the assumptions.

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

1. The independence assumption is nott met, but for this exercise, we will assume it is met (see above).
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 = log.Surv_time ~ Dose, data = poison)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.14507 -0.38506 -0.06518  0.37175  1.09357 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.1046     0.2396   8.782 7.13e-14 ***
## Dose         -0.5112     0.1509  -3.387  0.00103 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5104 on 94 degrees of freedom
## Multiple R-squared:  0.1088, Adjusted R-squared:  0.0993 
## F-statistic: 11.47 on 1 and 94 DF,  p-value: 0.001032

Compute the 95% confidence interval on the model parameters:

confint(log.model)

##                 2.5 %     97.5 %
## (Intercept)  1.628792  2.5804330
## Dose        -0.810830 -0.2115505

5.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 of survival time, we don’t have direct inference on survival time.

As such, we may interpret the output as follows:

• The slope:

If the dose at which the poison is administered increases with 1 mg, then the mean of the natural logarithm of survival time will decrease by 0.5112 units (95% CI [-0.8109,-0.2116]). This decrease is significant on the 5% significance level (p=0.00103).

• The intercept:

The average of the natural logarithm of survival time for fish that were given a dose of 0 mg is 2.1046 (95% CI [1.6288, 2.5804]). This value is significantly different from 0 on the 5% significance level (7.13e-14). This conclusion will be discussed in more detail below.

5.3 Interpretation on the original scale

The interpretation the log-scale is quite difficult: the natural logarithm survival time is not exactly a casual measure. To ease the interpetation, we may backtransform the results to the original scale (time in minutes). This we can do by taking the exponent of the outcomes:

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

## (Intercept)        Dose 
##   8.2039216   0.5997813

and of their confidence intervals:

exp(confint(log.model))

##                2.5 %     97.5 %
## (Intercept) 5.097711 13.2028535
## Dose        0.444489  0.8093284

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

• The slope:

If the dose at which the poison is administered increases with 1 mg, then the geometric mean of survival time will decrease by a factor of 0.5998 (95% CI [0.4445,0.8093]). This decrease is significant on the 5% significance level (p=0.00103).

• The intercept:

The geometric mean of the survival time for fish that were given a dose of 0 mg is 8.2039 minutes (95% CI [5.0977, 13.2029]). This value is significantly different from 0 on the 5% significance level (7.13e-14).

To go a little further on the latter interpretation, our analysis seems to suggest that fish in uncontaminated water would only suvive for about 8 minutes! This seems a little bit fishy…

Here, we make a diagnostic plot to explain this behavior:

plot(x=poison$Dose, y=poison$Surv_time, xlab="Dose (mg)", ylab="Survival time (min)", xlim=c(0,2.2))

abline(lsfit(poison$Dose,poison$Surv_time), lty=2) ## lty=2 draws a dashed line

## add the intercept as a point
summary_model <- summary(model)
points(0,summary_model$coefficients["(Intercept)","Estimate"],col="red", pch=16)

We can see that the dose value of 0 mg does not fall within the same range as our data. As such, it is unlikely that we can assume this model to have correct inference on the intercept: we expect that the linearity assumption will not be met for these low values. The model will only be trustworthy on the range of data that we use to build the model on!

This is not really a problem: we are mainly interested in the slope parameter in the given range of datapoints.

6 Conclusion

• On the log scale:

If the dose at which the poison is administered increases with 1 mg, then the mean of the natural logarithm of survival time will decrease by 0.5112 units (95% CI [-0.8109,-0.2116]). This decrease is highly significant on the 5% significance level (p=0.00103).

• On the original scale:

If the dose at which the poison is administered increases with 1 mg, then the geometric mean of survival time will decrease by a factor of 0.5998 (95% CI [0.4445,0.8093]). This decrease is highly significant on the 5% significance level (p=0.00103).

Again, note that in this dataset the assumption of independence is actually not met. In the tutorial of multiple regression, we will revisit this exercise. More specifically, we will there study the association between dose and survival while accounting for differences between species and fish of different weight.

---

Home