Shapiro-Wilk normality test:

It is used to test the normality in the distribution of the data for a variable. It can be used in groups with no more than 5000 observations (N < 5000).

Hypothesis

• H0 = The variable shows a normal distribution
  
• H1 = The variable does not show a normal distribution
  

Interpretation

• p-value > alfa: Do not reject H0 (normal)
  
• p-value < alfa: Reject H0 (not normal)
  

Note
hypothetical alpha 5% (0,05)

shapiro.test(penguins$bill_length_mm)

## 
##  Shapiro-Wilk normality test
## 
## data:  penguins$bill_length_mm
## W = 0.97485, p-value = 1.12e-05

Bartlett homoscedasticity test

It is used to test for homogeneity of variance (homoscedasticity) in k groups of samples, where k can be greater than two. It is adapted for normally distributed data.

Hypothesis

• H0 = Groups show homoscedasticity
  
• H1 = Groups (at least 2) do not show homoscedasticity.
  

Interpretación

• p-value > alfa: Do not reject H0 (homocedastics)
  
• p-value < alfa: Reject H0 (not homocedastics)
  

# Una variable independiente
bartlett.test(bill_length_mm ~ sex, data = penguins)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  bill_length_mm by sex
## Bartlett's K-squared = 1.3429, df = 1, p-value = 0.2465

# Múltiples variables independientes
bartlett.test(bill_length_mm ~ interaction(sex,species), data = penguins)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  bill_length_mm by interaction(sex, species)
## Bartlett's K-squared = 22.178, df = 5, p-value = 0.0004843

Fligner-Killeen homoscedasticity test

To test homoscedasticity in k groups of samples, where k can be greater than two. More robust against deviations from normality or when there are problems related to outliers.

Hypothesis

• H0 = Groups show homoscedasticity
  
• H1 = Groups (at least 2) do not show homoscedasticity.
  

Interpretación

• p-value > alfa: Do not reject H0 (homocedastics)
  
• p-value < alfa: Reject H0 (not homocedastics)
  

# One independent variable
fligner.test(bill_length_mm ~ sex, data = penguins)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  bill_length_mm by sex
## Fligner-Killeen:med chi-squared = 3.0297, df = 1, p-value = 0.08175

# Multiple independent variables
fligner.test(bill_length_mm ~ interaction(sex,species), data = penguins)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  bill_length_mm by interaction(sex, species)
## Fligner-Killeen:med chi-squared = 5.9733, df = 5, p-value = 0.3088

ANOVA (One-way)

ANOVA is a statistical technique used to compare the means of two or more groups. ANOVA requires that the groups meet the 3 assumptions.

Important!

It is worth noting that none of the continuous variables in our dataset penguins has a normal distribution, so in this example we will use the variable bill_length_mm only for demonstrative purposes, and considering that the result obtained is not reliable.

peng.aov <- aov(bill_length_mm ~ species, data = penguins)
summary(peng.aov)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## species       2   7194    3597   410.6 <2e-16 ***
## Residuals   339   2970       9                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness

Post-Hoc test for ANOVA

Tukey’s HSD (honestly significant difference) test

TukeyHSD(peng.aov)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = bill_length_mm ~ species, data = penguins)
## 
## $species
##                       diff       lwr        upr     p adj
## Chinstrap-Adelie 10.042433  9.024859 11.0600064 0.0000000
## Gentoo-Adelie     8.713487  7.867194  9.5597807 0.0000000
## Gentoo-Chinstrap -1.328945 -2.381868 -0.2760231 0.0088993

ANOVA (two-way)

Two-way ANOVA allows us to analyze the variance in the mean between groups when there is interaction between two different categorical and independent variables.

peng.aov2 <- aov(bill_length_mm ~ species*sex, data = penguins)
summary(peng.aov2)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## species       2   7015    3508 654.189 <2e-16 ***
## sex           1   1136    1136 211.807 <2e-16 ***
## species:sex   2     24      12   2.284  0.103    
## Residuals   327   1753       5                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 11 observations deleted due to missingness

TukeyHSD(peng.aov2)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = bill_length_mm ~ species * sex, data = penguins)
## 
## $species
##                       diff       lwr        upr     p adj
## Chinstrap-Adelie 10.009851  9.209424 10.8102782 0.0000000
## Gentoo-Adelie     8.744095  8.070779  9.4174106 0.0000000
## Gentoo-Chinstrap -1.265756 -2.094535 -0.4369773 0.0010847
## 
## $sex
##                 diff      lwr      upr p adj
## MALE-FEMALE 3.693368 3.194089 4.192646     0
## 
## $`species:sex`
##                                      diff       lwr        upr     p adj
## Chinstrap:FEMALE-Adelie:FEMALE   9.315995  7.937732 10.6942588 0.0000000
## Gentoo:FEMALE-Adelie:FEMALE      8.306259  7.138639  9.4738788 0.0000000
## Adelie:MALE-Adelie:FEMALE        3.132877  2.034137  4.2316165 0.0000000
## Chinstrap:MALE-Adelie:FEMALE    13.836583 12.458320 15.2148470 0.0000000
## Gentoo:MALE-Adelie:FEMALE       12.216236 11.064727 13.3677453 0.0000000
## Gentoo:FEMALE-Chinstrap:FEMALE  -1.009736 -2.443514  0.4240412 0.3338130
## Adelie:MALE-Chinstrap:FEMALE    -6.183118 -7.561382 -4.8048548 0.0000000
## Chinstrap:MALE-Chinstrap:FEMALE  4.520588  2.910622  6.1305547 0.0000000
## Gentoo:MALE-Chinstrap:FEMALE     2.900241  1.479553  4.3209291 0.0000002
## Adelie:MALE-Gentoo:FEMALE       -5.173382 -6.341002 -4.0057622 0.0000000
## Chinstrap:MALE-Gentoo:FEMALE     5.530325  4.096547  6.9641020 0.0000000
## Gentoo:MALE-Gentoo:FEMALE        3.909977  2.692570  5.1273846 0.0000000
## Chinstrap:MALE-Adelie:MALE      10.703707  9.325443 12.0819703 0.0000000
## Gentoo:MALE-Adelie:MALE          9.083360  7.931851 10.2348686 0.0000000
## Gentoo:MALE-Chinstrap:MALE      -1.620347 -3.041035 -0.1996591 0.0149963

Kruskal-Wallis

Kruskal-Wallis test is a statistical and nonparametric test (it does not assume normality in the data) used to compare the means of two or more groups. It is similar to ANOVA but with data by categories.

kruskal.test(bill_length_mm ~ species, data = penguins)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  bill_length_mm by species
## Kruskal-Wallis chi-squared = 244.14, df = 2, p-value < 2.2e-16

pairwise.wilcox.test(penguins$bill_length_mm, penguins$species, p.adjust.method="fdr")

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  penguins$bill_length_mm and penguins$species 
## 
##           Adelie Chinstrap
## Chinstrap <2e-16 -        
## Gentoo    <2e-16 0.0018   
## 
## P value adjustment method: fdr

Densiplot

ggplot(data = penguins, 
       mapping = aes(x = bill_length_mm)) + 
       geom_density()

The densiplot calculates and plots the density estimate of the data (it is a smoothed version of the histogram). A densiplot allows us to see the distribution of our data, and the area under the densiplot curve adds up to 1 which corresponds to the probability of finding 100% of our observations.

The {ggplot2} package, allows us to perform several types of visualizations due to the large variety of geom objects. As we have already reviewed our data set “penguins”, and we have identified continuous variables such as “bill_length_mm” (peak length in millimeters) and “bill_depth_mm” (peak depth in millimeters). Next we will make a point plot (selected through the geom_point object), where we will map two variables (i.e. x and y) and not only one (i.e. x) as in the previous example. To do this, we will map bill_length_mm to the x axis, and bill_depth_mm to the y axis.

Scatterplot

ggplot(data = penguins, 
       mapping = aes(x = bill_length_mm, 
                     y = bill_depth_mm)) + 
       geom_point()

Components of a plot in {ggplot2}

• ggplot() : It is the function that creates a coordinate system - in general form, i.e. for all later incorporated layers - that will be incorporated into the layers. The first argument of the function is the data set (data =). By itself this function does not generate a layer, but it provides all the information needed to add one.
  

• data : The dataset (in this case penguins) is a rectangular collection of data with the variables (columns) and their observations/values (rows) to be mapped.
  

• mapping (aes()) : This is where you specify the set of variables and observations that are “mapped” or assigned to the visual properties to be used in the graph and which axes to assign to these values (x = bill_length_mm, y = bill_depth_mm). If they are not specified -in general-, they must be indicated in each layer added to the chart.
  

• geom_point() : The plot layers are incorporated by Geom functions. In this case the function adds a layer of points to the plot (i.e., scatterplot). ggplot2 includes more than 30 geom functions, in addition to those developed by other authors.
  

Boxplot

ggplot(data = penguins, 
       mapping = aes(x = species, 
                     y = bill_length_mm,
                     color = species)) + 
       geom_boxplot()

The boxplot is a very useful type of visualization because it allows us to identify in a single graph the quantiles (0.25 (Q1), 0.5 (Q2) and 0.75 (Q3)), as well as the outliers (outliers) of a data set.

How to interpret a boxplot?

Violin Plot

ggplot(data = penguins, 
       mapping = aes(x = species, 
                     y = bill_length_mm,
                     color = species)) + 
       geom_violin()

The violin plot is another type of visualization, in which we can observe the distribution of the data, allowing us to observe the density more clearly.

Normal distribution

ggplot(data = penguins, 
       mapping = aes(x = bill_length_mm)) + 
       geom_density() +
stat_overlay_normal_density(color = "darkblue", linetype = "dashed")

With this function we can see that our distribution is not perfectly normal, which visually verifies the result obtained previously with the Shapiro-Wilk test.

Plot an analysis of variance

The following function stat_compare_means() will allow us to perform an analysis of variance with the Kruskal-Wallis test on our boxplot. The second comparison includes the Wilcoxon post-hoc rank sum test, in the same way as we performed the previous example with the Kruskal-Wallis test. By default, the stat_compare_means() function calculates the Kruskal-Wallis test when there are 3 or more groups. The ref.group = ".all." argument refers to a comparison of all groups against all groups. The label.y argument specifies where (with respect to the y axis) our legends will be placed on the graph.

ggplot(data = penguins, 
       mapping = aes(x = species, 
                     y = bill_length_mm,
                     color = species)) + 
       geom_boxplot() +
  stat_compare_means() +
  stat_compare_means(aes(label = ..p.signif..), 
                         method = "wilcox.test", 
                         ref.group = ".all.", 
                         label.y = 65)

If we want to perform the ANOVA test we must follow the following code, where the paired t-test will be performed. But let’s remember that in this case our data does not comply with normality, so this example is just a demonstration of the code.

ggplot(data = penguins, 
       mapping = aes(x = species, 
                     y = bill_length_mm,
                     color = species)) + 
       geom_boxplot() +
  stat_compare_means(method = "anova") +
  stat_compare_means(aes(label = ..p.signif..), 
                         method = "t.test", 
                         ref.group = ".all.", 
                         label.y = 65)