Toothlength R dataset analysis

An application of an ANOVA model

Dataset

This is a classic R dataset and it concerns the vitamin C effect on tooth growth in Guinea pigs: there are sixty pigs and each one received one dose of vitamin C with three dosages (0.5, 1, and 2 mg/day) and each dosage can be delivered in two ways: by orange juice or by ascorbic acid (a form of vitamin C and coded as VC).

The sample stratification is the following one

sampleTable = table(myData$supp, myData$dose)
rownames(sampleTable) = c("Orange Juice", "Ascorbic acid")

sampleKable = kable(sampleTable) %>% 
  add_header_above(c("", "Dose" = 3))

sampleKable %>% kable_material(full_width = F) %>% 
      kableExtra::landscape()

First of all we proceed by a graphical inspection of the data: here we have the scatterplot and the boxplot:

myDataPlot = myData %>% group_by(supp) %>% plot_ly(x = ~dose, y = ~len, name = ~supp, type = 'scatter', mode = 'markers')


myDataBox =  myData %>% group_by(supp) %>% plot_ly(x = ~dose, y = ~len, name = ~supp, type = 'box')

plotS = subplot(myDataPlot, myDataBox, nrows = 1) 

config(plotS, displayModeBar = FALSE)

From the plots it is clear that the tooth growth differs by method of administration and dosage, also if it appears that at the highest dosage the main difference is in the variability of the length, more than the size.

Studying dependencies

Now, we need a confirmation of these qualitative observations, so we begin comparing the mean between the groups: applying directly apply a T test to the tooth length by administration method getting

> t.test(len ~ supp, data = myData)

	Welch Two Sample t-test

data:  len by supp
t = 1.9153, df = 55.309, p-value = 0.06063
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
 -0.1710156  7.5710156
sample estimates:
mean in group OJ mean in group VC 
        20.66333         16.96333 

Here we notice that the test result says that we can not reject the null hypothesis at level 5%. Furthermore, we point out that the value of the T statics is equal to the square of the F statistic in the one way ANOVA model of the length on the method of supply, but the p-value can be different for differents methods applied to calculate the degree of freedom in the two models.

From this result we would behave as if the supply method does not influence the length, but from our graphical exploratory analysis we know that this conclusion is suspect, indeed it seems that we have some kind of dependence. Why our statistical instruments does not confirm this dependence? They does not work well? No, this is not the case. One possible answer is that we are neglecting some aspect of the phenomenon, in this case we have a confounding effect due to the dosage variable, so we have to investigate to find out if this is the case. To this aim, now we test administration method effect for each given dosage level

tmpData = myData %>% group_split(dose, .keep = FALSE) %>% setNames(unique(myData$dose))


lapply(tmpData, function(l){t.test(len ~ supp, data = l)})

$`0.5`

	Welch Two Sample t-test

data:  len by supp
t = 3.1697, df = 14.969, p-value = 0.006359
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
 1.719057 8.780943
sample estimates:
mean in group OJ mean in group VC 
           13.23             7.98 


$`1`

	Welch Two Sample t-test

data:  len by supp
t = 4.0328, df = 15.358, p-value = 0.001038
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
 2.802148 9.057852
sample estimates:
mean in group OJ mean in group VC 
           22.70            16.77 


$`2`

	Welch Two Sample t-test

data:  len by supp
t = -0.046136, df = 14.04, p-value = 0.9639
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
 -3.79807  3.63807
sample estimates:
mean in group OJ mean in group VC 
           26.06            26.14 

From these results we can see that, given one of the first two dosage levels , the administration method has a big influence on the tooth length, but not the same for the third. Now, this is what we have expected from the above exploratory data inspection.

Remark

We have applied a Welch two sample T test to each of the groups defined above, this is an adaptation of Student’s T test to the case of two populations with unknown variances and not necessarily equals, but assuming normal distribution for the length.

One can also test if the dosage condition length inside administration classes and this can be done in the same way. But as we will see soon there are also alternative methods.

ANOVA model

We can analyze this dataset also by ANOVA models. First of all, as remarked above we need a two way ANOVA model, indeed if for istance we apply a one way ANOVA model of the length on the supply method, we get

anova(aov(data = myData, len ~ supp))

and in this case as with the Welch test we can not reject the null hypothesis, so we need a two-way ANOVA.

Now, to apply the ANOVA we can proceed in two way: first, we can consider the dose variable as a categorical variable with three levels and in this case we have a standard two way ANOVA model; second, we consider the dose variable as a quantitative variable and in this case we have an ANCOVA model, that is an analysis of covariance model.

In the first case the results are the following

myData %<>% mutate(doseFactor = as.factor(dose))
myDataLm = lm(data = myData, 
                len ~ supp * doseFactor,
                contrasts = list(supp = "contr.sum", 
                                  doseFactor = "contr.sum"))
anova(myDataLm)

From the above table we can see that all the variables are significant, furthermore the interaction too is significant.

In the second approach we consider the dose as a quantitative variable, in this case we point out that this can be suitable, indeed the dose levels are not equally spaced.

In this case we estimate the linear model with and without interaction to see the difference from a prediction point of view. The model result now is

myDataLm = lm(data = myData, len ~ supp + dose)
myDataLmInterAct = lm(data = myData, len ~ supp * dose)
summary(myDataLmInterAct)

Call:
lm(formula = len ~ supp * dose, data = myData)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.2264 -2.8463  0.0504  2.2893  7.9386 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   11.550      1.581   7.304 1.09e-09 ***
suppVC        -8.255      2.236  -3.691 0.000507 ***
dose           7.811      1.195   6.534 2.03e-08 ***
suppVC:dose    3.904      1.691   2.309 0.024631 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.083 on 56 degrees of freedom
Multiple R-squared:  0.7296,	Adjusted R-squared:  0.7151 
F-statistic: 50.36 on 3 and 56 DF,  p-value: 6.521e-16

aovKable(aov(myDataLmInterAct))

From the analysis of variance table we can conclude that the administration method and dosage are factors influencing the tooth growth, furthermore the phenomenon react to the dosage differently with the two administration method, indeed the estimated parameter ‘r suppVC:dose’ is significantly different from zero and this represent the difference in the line slope in the VC group. This aspect can be easily seen in the following plot contains data and prediction with or without interaction

myDataPlot =  myDataPlot %>% add_trace(x = myData$dose, y = myDataLm$fitted.values, mode = 'lines', line = list(width = 2.0))

myDataPlot %>% add_lines(x = myData$dose, y = myDataLmInterAct$fitted.values, size = 2, name = ~paste('inter', supp, sep = '-'), line = list(width = 2, dash = 'dash')) %>%
  config(displayModeBar = FALSE)