Exercise 8.1: Additive linear model on the poison dataset - solution

1 The poison 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

In this tutorial session we will focus on Dojofish, and we will model the survival time in function of the poison dose while correcting for the weight of the fish.

1. We will first analyse the survival data by only considering the dose as an explanatory variable for survival time

2. Next we will model the survival data with and additive model for dose and weight

Load libraries

library(tidyverse)
library(ggplot2)
#install.packages("GGally")
library(GGally)
library(car)
library(multcomp)

3 Import the data

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

4 Data tidying

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

1. Capitalise the fist column name

2. Set the Species column as a factor

3. Change the species factor levels from “0” to Dojofish. Hint: use the fct_recode function.

4. In the previous analysis on this dataset (Simple linear regression session), we performed a log-transformation on the response variable Surv_time to meet the normality and homoscedasticity assumptions of the linear model. Here, we will immediately work with log-transformed survival times; store these in the new variable log2Surv_time and remove the non-transformed values.

5. Subset the data to only retain Dojofish.

poison <- poison %>%
  rename("Species" = "species") %>%
  mutate(Species = as.factor(Species)) %>%
  mutate(Species = fct_recode(Species, Dojofish = "0", Goldfish = "1", Zebrafish = "2")) %>%
  mutate(log2Surv_time = Surv_time %>% log2) %>%
  dplyr::select(-Surv_time) %>%
  filter(Species=="Dojofish")

poison

5 Data exploration

Prior to the analysis, we should explore our data. To start our data exploration, we will make use of the ggpairs function of the GGally R package. This function will generate a visualization containing multiple panels, which display (1) univariate plots of the different variables in our dataset, (2) bivariate plots and (3) correlation coefficients between the different variables.

poison %>%
  ggpairs + theme_bw()

Based on these plots, we observe that:

• The survival time seems to be associated with dose and fish weight.

From the tutorial of H6 we have seen that the fish weights were not nicely uniform across the different poison dosages due to the randomisation.

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

6 Simple linear regression

This is the same regression model that we have already fit in the exercise session on simple linear regression, with Dose as the only explanatory variable for log2Surv_time.

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

## display the diagnostic plots of the model
par(mfrow=c(2,2))
plot(lm_simple)

1. The independence assumption was met because the fish were randomized to the dose.
2. The linearity assumption is met.
3. The normality assumption is met.
4. The homoscedasticity assumption is met.

Finally, we look at the output of the model.

summary(lm_simple)

## 
## Call:
## lm(formula = log2Surv_time ~ Dose, data = poison)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6891 -0.3804 -0.1076  0.3452  1.0856 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.0469     0.3462   8.800 1.33e-10 ***
## Dose         -0.9063     0.2208  -4.104 0.000215 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4828 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

Or for an interpretation at the original scale (minutes in stead of log2 minutes):

2^(lm_simple$coef)["Dose"]

##      Dose 
## 0.5335551

2^(confint(lm_simple))[2,]

##     2.5 %    97.5 % 
## 0.3912681 0.7275855

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 0.53 per gram of poison that is added (95% CI [0.39, 0.73]).

7 Analysis with additive effect for weight

7.1 Model specification

Here, we will estimate the effect of the poison while correcting for weight and we will add it as an additional covariate to our linear regression model, such that

\[ y_i=\beta_0+\beta_d x_d + \beta_g x_g + \epsilon_i, \]

with \(\epsilon_i \text{ i.i.d. } N(0,\sigma^2)\).

7.2 Assumptions

The model will again be fit to allow for assessing the model assumptions

lm_additive <- lm(log2Surv_time ~ Dose + Weight, data = poison)
par(mfrow=c(2,2))
plot(lm_additive)

The assumption of independence, linearity and homoscedasticity are met.

The QQ-plot suggest that 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. 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 <- lm_additive %>% sigma 

dataHlp <- poison
simModels <- list()

par(mfrow=c(3,3))
for (i in 1:9)
{
  nobs <- poison %>% nrow
  dataHlp$ySim <- lm_additive$fit + rnorm(nobs, sd = sigma)
  simModels[[i]] <- lm(ySim~Dose + Weight, dataHlp)
  plot(simModels[[i]], which = 2)
}

It seems that deviations of the size that we see in our real data may be expected even when normality is met. As such, all assumptions for linear regression seem to be valid.

7.3 Inference

We then inspect the results.

summary(lm_additive)

## 
## Call:
## lm(formula = log2Surv_time ~ Dose + Weight, data = poison)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.59629 -0.33110 -0.06836  0.32507  0.83315 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.8294     0.6457   1.285 0.207119    
## Dose         -0.9590     0.1888  -5.081 1.17e-05 ***
## Weight        1.0783     0.2792   3.862 0.000451 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4116 on 36 degrees of freedom
## Multiple R-squared:  0.5141, Adjusted R-squared:  0.4871 
## F-statistic: 19.04 on 2 and 36 DF,  p-value: 2.282e-06

7.4 Interpretation of model parameters

We see that the effect of dose on survival time remains similar, however, it has become more significant after we have incorporated weight in our model. Indeed, from the data exploration, we learned that weight is associated with survival.

1. As such, by incorporating weight in our model, we are able to explain a larger part of the variability in the response variable survival time. As a consequence, the variability in the residuals of the model will decrease, which in turn will lead to smaller standard error estimates for the different parameter estimates in the model.

2. From the data exploration, we additionally found that the dojo-fish weights were not uniform across the different poison dosages due to the randomisation. Therefore, we can estimate the effect between dose and survival time better while accounting for the weight.

In this model, the effect of Dose can be interpreted as the average change in the log2 survival time between two groups of dojofish with the same weight that are exposed to a poison dosage that differs 1 mg/L. In symbols:

\[\begin{eqnarray} \hat{\mu_1}&=& \beta_0 + \beta_d x_{1d} + \beta_g x_g \text{ (average log2-survival time for dose 1 for a certain weight)}\\ \hat{\mu_2}&=& \beta_0 + \beta_d x_{2d} + \beta_g x_g \text{ (average log2-survival time for dose 2 for that same weight)}\\ \hat \mu_2- \hat \mu_1&=&\beta_0 + \beta_d x_{2d} + \beta_g x_g - (\beta_0 + \beta_d x_{1d} + \beta_g x_g) \text{ (difference in average log2-survival time between dose 2 en dose 1)}\\ \hat \mu_2-\hat \mu_1&=&\beta_d (x_{2d}-x_{1d}) \end{eqnarray}\]

7.5 Conclusion

2^(lm_additive$coef)

## (Intercept)        Dose      Weight 
##    1.777004    0.514402    2.111484

2^(confint(lm_additive))

##                 2.5 %    97.5 %
## (Intercept) 0.7169728 4.4042734
## Dose        0.3945101 0.6707292
## Weight      1.4259918 3.1265018

The dose of the poison has an extremely significant effect on the log2 transformed survival time of dojofish (p-value = 1e-05). The geometric average of the survival time for dojofish that are exposed to a poison dose that is 1mg/L larger is approximately halved (factor = \(2^{\beta_d}=\) 0.51) .

The effect of weight on the survival time of dojofish is also extremely significant (p-value = 5e-04). The geometric average of the survival time of a dojofish that weighs 1 gram more than another dojofish is approximately twice as long (factor = \(2^{ \beta_g}\)= 2.1, 95% BI [1.4,3.1]).

7.6 Remarks

In the lm_additive model, we included only a main effect for weight. However, there could also be an interaction effect between weight and dose. An interaction between weight and dose implies that the dose effect on the survival time changes according to the weight of the fish.

---

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