Statistical Inference

1 Preliminarities

• This document is inspired by the scripts of the course Data Science Academy: Master Data Science In R

• An interactive version of this document is available in R Studio Cloud

• The packages needed to execute this document are as follows:

library(tidyverse) # modern data science library
library(ggpubr) # for beautiful figures using ggplot in the background
library(car) # companion to applied regression analysis

2 What is Hypothesis Testing?

3 T-Test

4 Non-Parametric T-Tests

5 One-way ANOVA

• An extension of independent two-samples t-test
• Useful for comparing means when we have two or more groups
• Data is organized into two or more groups based on one single grouping variable (also called factor variable)

data("PlantGrowth")
PlantGrowth

str(PlantGrowth)

'data.frame':   30 obs. of  2 variables:
 $ weight: num  4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
 $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...

summary(PlantGrowth)

     weight       group   
 Min.   :3.590   ctrl:10  
 1st Qu.:4.550   trt1:10  
 Median :5.155   trt2:10  
 Mean   :5.073            
 3rd Qu.:5.530            
 Max.   :6.310            

levels(PlantGrowth$group)

[1] "ctrl" "trt1" "trt2"

ggboxplot(PlantGrowth, x = "group", y = "weight", 
          color = "group", palette = c("blue", "red", "green"),
          order = c("ctrl", "trt1", "trt2"),
          ylab = "Weight", xlab = "Treatment")

md <- aov(weight ~ group, data = PlantGrowth)
summary(md) 

            Df Sum Sq Mean Sq F value Pr(>F)  
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Since the p-value is less than 5 %, then there are significant differences between groups

5.1 Which of the groups are significantly different?

5.1.1 Post hoc test

t <- TukeyHSD(md)
t

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = weight ~ group, data = PlantGrowth)

$group
            diff        lwr       upr     p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

plot(t)

Because only one of the comparison (trt2-trt1) shows a p-value is less than 5 %, then there is a statistically difference only between the groups trt2 and trt1

5.1.2 Pairwise t-test

pairwise.t.test(PlantGrowth$weight, PlantGrowth$group,
                p.adjust.method = "BH")

    Pairwise comparisons using t tests with pooled SD 

data:  PlantGrowth$weight and PlantGrowth$group 

     ctrl  trt1 
trt1 0.194 -    
trt2 0.132 0.013

P value adjustment method: BH 

5.2 Conditions of one way ANOVA

5.2.1 Homogeneity of the variance of the residuals

plot(md, 1)

The variance of the residuals seems to be largely homogeneous

5.2.2 Normality of residuals.

plot(md, 2)

# Extract the residuals
aov_residuals <- residuals(object = md )
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals ) 

    Shapiro-Wilk normality test

data:  aov_residuals
W = 0.96607, p-value = 0.4379

Since the p-value is more than 5% in the Shapiro-Wilk normality test, then the errors are normally distributed

6 Non-parametric ANOVA

7 Two way ANOVA

• Useful to evaluate the effect of two factor variables on a response variable
• Useful to evaluate the interaction effect of two factor variables The nule hypothesis is H0: response mean for all factor levels are equal

data("ToothGrowth")
df <- ToothGrowth
df

summary(df)

      len        supp         dose      
 Min.   : 4.20   OJ:30   Min.   :0.500  
 1st Qu.:13.07   VC:30   1st Qu.:0.500  
 Median :19.25           Median :1.000  
 Mean   :18.81           Mean   :1.167  
 3rd Qu.:25.27           3rd Qu.:2.000  
 Max.   :33.90           Max.   :2.000  

str(df)

'data.frame':   60 obs. of  3 variables:
 $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
 $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
 $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

We need to convert the variable dose into a factor variable

df$dose <- factor(df$dose, 
                       levels = c(0.5, 1, 2),
                       labels = c("D0.5", "D1", "D2"))

boxplot(len ~ supp * dose, data=df, frame = FALSE, 
        col = c("red", "blue"), ylab="Tooth Length")

Let’s see if tooth length depends on supp and dose.

7.1 Assuming that two factor variables are independent

md1 <- aov(len ~ supp + dose, data = df)
summary(md1)

            Df Sum Sq Mean Sq F value   Pr(>F)    
supp         1  205.4   205.4   14.02 0.000429 ***
dose         2 2426.4  1213.2   82.81  < 2e-16 ***
Residuals   56  820.4    14.7                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Both supp and dose are statistically significant

Assuming Mean tooth length for both supp and dose are NOT equal

7.2 Testing the interaction effect

Let’s test wheater two variables might interact to create a interactive a synergistic effect

md2 <- aov(len ~ supp * dose, data = df)
summary(md2)

            Df Sum Sq Mean Sq F value   Pr(>F)    
supp         1  205.4   205.4  15.572 0.000231 ***
dose         2 2426.4  1213.2  92.000  < 2e-16 ***
supp:dose    2  108.3    54.2   4.107 0.021860 *  
Residuals   54  712.1    13.2                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The two main effects (supp and dose) are statistically significant, as well as their interaction

• relationships between dose and tooth length depends on the supp method
• influence the difference between mean tooth length

7.3 Post Hoc test

To identify which dosage group is different

TukeyHSD(md2, which = "dose")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = len ~ supp * dose, data = df)

$dose
          diff       lwr       upr   p adj
D1-D0.5  9.130  6.362488 11.897512 0.0e+00
D2-D0.5 15.495 12.727488 18.262512 0.0e+00
D2-D1    6.365  3.597488  9.132512 2.7e-06

7.4 Pairwise t-tests

pairwise.t.test(df$len, df$dose,
                p.adjust.method = "BH")

    Pairwise comparisons using t tests with pooled SD 

data:  df$len and df$dose 

   D0.5    D1     
D1 1.0e-08 -      
D2 4.4e-16 1.4e-05

P value adjustment method: BH 

7.5 Testing when we have unequal sample numbers

mya <- aov(len ~ supp * dose, data = df)
Anova(mya, type = "III") # Type-III sums of squares used when we have unequal numbers per group

Anova Table (Type III tests)

Response: len
             Sum Sq Df F value    Pr(>F)    
(Intercept) 1750.33  1 132.730 3.603e-16 ***
supp         137.81  1  10.450  0.002092 ** 
dose         885.26  2  33.565 3.363e-10 ***
supp:dose    108.32  2   4.107  0.021860 *  
Residuals    712.11 54                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note that his Anova() functions is from the car package

7.6 Testing the condtions for a two-way anova

7.6.1 Homogeneity of variance of residuals

plot(md2, 1)

7.6.2 Normality of residuals

plot(md2, 2)

aov_residuals <- residuals(object = md2 )
shapiro.test(x = aov_residuals ) #errors are normally distributed

    Shapiro-Wilk normality test

data:  aov_residuals
W = 0.98499, p-value = 0.6694

8 Power Test for Detecting Effect