Case of large sample sizes

If the sample size is large enough (n > 30), we can ignore the distribution of the data and use parametric tests.

The central limit theorem tells us that no matter what distribution things have, the sampling distribution tends to be normal if the sample is large enough (n > 30).

However, to be consistent, normality can be checked by visual inspection [normal plots (histogram), Q-Q plot (quantile-quantile plot)] or by significance tests].

Visual methods

Density plot and Q-Q plot can be used to check normality visually.

1. Density plot: the density plot provides a visual judgment about whether the distribution is bell shaped.

library("ggpubr")
ggdensity(my_data$len, 
          main = "Density plot of tooth length",
          xlab = "Tooth length")

2. Q-Q plot: Q-Q plot (or quantile-quantile plot) draws the correlation between a given sample and the normal distribution. A 45-degree reference line is also plotted.

library(ggpubr)
ggqqplot(my_data$len)

It’s also possible to use the function qqPlot() [in car package]:

library("car")
qqPlot(my_data$len)

As all the points fall approximately along this reference line, we can assume normality.

Normality test

Visual inspection, described in the previous section, is usually unreliable. It’s possible to use a significance test comparing the sample distribution to a normal one in order to ascertain whether data show or not a serious deviation from normality.

There are several methods for normality test such as Kolmogorov-Smirnov (K-S) normality test and Shapiro-Wilk’s test.

The null hypothesis of these tests is that “sample distribution is normal”. If the test is significant, the distribution is non-normal.

Shapiro-Wilk’s method is widely recommended for normality test and it provides better power than K-S. It is based on the correlation between the data and the corresponding normal scores.

Note that, normality test is sensitive to sample size. Small samples most often pass normality tests. Therefore, it’s important to combine visual inspection and significance test in order to take the right decision.

The R function shapiro.test() can be used to perform the Shapiro-Wilk test of normality for one variable (univariate):

shapiro.test(my_data$len)

    Shapiro-Wilk normality test
data:  my_data$len
W = 0.96743, p-value = 0.1091

From the output, the p-value > 0.05 implying that the distribution of the data are not significantly different from normal distribution. In other words, we can assume the normality.

How to assess the normality of the data?

1. With large enough sample sizes (n > 30) the violation of the normality assumption should not cause major problems (central limit theorem). This implies that we can ignore the distribution of the data and use parametric tests.

2. However, to be consistent, we can use Shapiro-Wilk’s significance test comparing the sample distribution to a normal one in order to ascertain whether data show or not a serious deviation from normality.

How to assess the equality of variances?

The standard Student’s t-test (comparing two independent samples) and the ANOVA test (comparing multiple samples) assume also that the samples to be compared have equal variances.

If the samples, being compared, follow normal distribution, then it’s possible to use:

• F-test to compare the variances of two samples
• Bartlett’s Test or Levene’s Test to compare the variances of multiple samples.

Measure of central tendency: mean, median, mode

Roughly speaking, the central tendency measures the “average” or the “middle” of your data. The most commonly used measures include:

• the mean: the average value. It’s sensitive to outliers.
• the median: the middle value. It’s a robust alternative to mean.
• and the mode: the most frequent value

In R,

• The function mean() and median() can be used to compute the mean and the median, respectively;
• The function mfv() [in the modeest R package] can be used to compute the mode of a variable.

The R code below computes the mean, median and the mode of the variable Sepal.Length [in my_data data set]:

# Compute the mean value
mean(my_data$Sepal.Length)

[1] 5.843333

# Compute the median value
median(my_data$Sepal.Length)

[1] 5.8

# Compute the mode
# install.packages("modeest")
require(modeest)
mfv(my_data$Sepal.Length)

[1] 5

Measure of variablity

Measures of variability gives how “spread out” the data are.

# Compute the minimum value
min(my_data$Sepal.Length)

[1] 4.3

# Compute the maximum value
max(my_data$Sepal.Length)

[1] 7.9

# Range
range(my_data$Sepal.Length)

[1] 4.3 7.9

quantile(x, probs = seq(0, 1, 0.25))

quantile(my_data$Sepal.Length)

  0%  25%  50%  75% 100% 
 4.3  5.1  5.8  6.4  7.9 

quantile(my_data$Sepal.Length, seq(0, 1, 0.1))

IQR(my_data$Sepal.Length)

[1] 1.3

# Compute the variance
var(my_data$Sepal.Length)
# Compute the standard deviation =
# square root of th variance
sd(my_data$Sepal.Length)

# Compute the median
median(my_data$Sepal.Length)
# Compute the median absolute deviation
mad(my_data$Sepal.Length)

Computing an overall summary of a variable and an entire data frame

summary(my_data$Sepal.Length)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   5.100   5.800   5.843   6.400   7.900 

summary(my_data, digits = 1)

  Sepal.Length  Sepal.Width  Petal.Length  Petal.Width        Species  
 Min.   :4     Min.   :2    Min.   :1     Min.   :0.1   setosa    :50  
 1st Qu.:5     1st Qu.:3    1st Qu.:2     1st Qu.:0.3   versicolor:50  
 Median :6     Median :3    Median :4     Median :1.3   virginica :50  
 Mean   :6     Mean   :3    Mean   :4     Mean   :1.2                  
 3rd Qu.:6     3rd Qu.:3    3rd Qu.:5     3rd Qu.:1.8                  
 Max.   :8     Max.   :4    Max.   :7     Max.   :2.5                  

# Compute the mean of each column
sapply(my_data[, -5], mean)

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
    5.843333     3.057333     3.758000     1.199333 

# Compute quartiles
sapply(my_data[, -5], quantile)

     Sepal.Length Sepal.Width Petal.Length Petal.Width
0%            4.3         2.0         1.00         0.1
25%           5.1         2.8         1.60         0.3
50%           5.8         3.0         4.35         1.3
75%           6.4         3.3         5.10         1.8
100%          7.9         4.4         6.90         2.5

install.packages("pastecs")

# Compute descriptive statistics
library(pastecs)
res <- stat.desc(my_data[, -5])
round(res, 2)

             Sepal.Length Sepal.Width Petal.Length Petal.Width
nbr.val            150.00      150.00       150.00      150.00
nbr.null             0.00        0.00         0.00        0.00
nbr.na               0.00        0.00         0.00        0.00
min                  4.30        2.00         1.00        0.10
max                  7.90        4.40         6.90        2.50
range                3.60        2.40         5.90        2.40
sum                876.50      458.60       563.70      179.90
median               5.80        3.00         4.35        1.30
mean                 5.84        3.06         3.76        1.20
SE.mean              0.07        0.04         0.14        0.06
CI.mean.0.95         0.13        0.07         0.28        0.12
var                  0.69        0.19         3.12        0.58
std.dev              0.83        0.44         1.77        0.76
coef.var             0.14        0.14         0.47        0.64

Case of missing values

Note that, when the data contains missing values, some R functions will return errors or NA even if just a single value is missing.

For example, the mean() function will return NA if even only one value is missing in a vector. This can be avoided using the argument na.rm = TRUE, which tells to the function to remove any NAs before calculations. An example using the mean function is as follow:

mean(my_data$Sepal.Length, na.rm = TRUE)

Graphical display of distributions

The R package ggpubr will be used to create graphs.

# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

install.packages("ggpubr")

library(ggpubr)

ggboxplot(my_data, y = "Sepal.Length", width = 0.5)

gghistogram(my_data, x = "Sepal.Length", bins = 9, 
             add = "mean")

ggecdf(my_data, x = "Sepal.Length")

ggqqplot(my_data, x = "Sepal.Length")

Create some data

Distribution of hair and eye color by sex of 592 students:

# Hair/eye color data
df <- as.data.frame(HairEyeColor)
hair_eye_col <- df[rep(row.names(df), df$Freq), 1:3]
rownames(hair_eye_col) <- 1:nrow(hair_eye_col)
head(hair_eye_col)

   Hair   Eye  Sex
1 Black Brown Male
2 Black Brown Male
3 Black Brown Male
4 Black Brown Male
5 Black Brown Male
6 Black Brown Male

# hair/eye variables
Hair <- hair_eye_col$Hair
Eye <- hair_eye_col$Eye

Simple frequency distribution: one categorical variable

• Table of counts

# Frequency distribution of hair color
table(Hair)

Hair
Black Brown   Red Blond 
  108   286    71   127 

# Frequency distribution of eye color
table(Eye)

Eye
Brown  Blue Hazel Green 
  220   215    93    64 

• Graphics: to create the graphics, we start by converting the table as a data frame.

# Compute table and convert as data frame
df <- as.data.frame(table(Hair))
df

   Hair Freq
1 Black  108
2 Brown  286
3   Red   71
4 Blond  127

# Visualize using bar plot
library(ggpubr)
ggbarplot(df, x = "Hair", y = "Freq")

Two-way contingency table: Two categorical variables

tbl2 <- table(Hair , Eye)
tbl2

       Eye
Hair    Brown Blue Hazel Green
  Black    68   20    15     5
  Brown   119   84    54    29
  Red      26   17    14    14
  Blond     7   94    10    16

It’s also possible to use the function xtabs(), which will create cross tabulation of data frames with a formula interface.

xtabs(~ Hair + Eye, data = hair_eye_col)

• Graphics: to create the graphics, we start by converting the table as a data frame.

df <- as.data.frame(tbl2)
head(df)

   Hair   Eye Freq
1 Black Brown   68
2 Brown Brown  119
3   Red Brown   26
4 Blond Brown    7
5 Black  Blue   20
6 Brown  Blue   84

# Visualize using bar plot
library(ggpubr)
ggbarplot(df, x = "Hair", y = "Freq",
          color = "Eye", 
          palette = c("brown", "blue", "gold", "green"))

# position dodge
ggbarplot(df, x = "Hair", y = "Freq",
          color = "Eye", position = position_dodge(),
          palette = c("brown", "blue", "gold", "green"))

Multiway tables: More than two categorical variables

• Hair and Eye color distributions by sex using xtabs():

xtabs(~Hair + Eye + Sex, data = hair_eye_col)

, , Sex = Male

       Eye
Hair    Brown Blue Hazel Green
  Black    32   11    10     3
  Brown    53   50    25    15
  Red      10   10     7     7
  Blond     3   30     5     8

, , Sex = Female

       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

• You can also use the function ftable() [for flat contingency tables]. It returns a nice output compared to xtabs() when you have more than two variables:

ftable(Sex + Hair ~ Eye, data = hair_eye_col)

      Sex   Male                 Female                
      Hair Black Brown Red Blond  Black Brown Red Blond
Eye                                                    
Brown         32    53  10     3     36    66  16     4
Blue          11    50  10    30      9    34   7    64
Hazel         10    25   7     5      5    29   7     5
Green          3    15   7     8      2    14   7     8

Compute table margins and relative frequency

Table margins correspond to the sums of counts along rows or columns of the table. Relative frequencies express table entries as proportions of table margins (i.e., row or column totals).

The function margin.table() and prop.table() can be used to compute table margins and relative frequencies, respectively.

1. Format of the functions:

margin.table(x, margin = NULL)

prop.table(x, margin = NULL)

• x: table
• margin: index number (1 for rows and 2 for columns)

2. compute table margins:

Hair <- hair_eye_col$Hair
Eye <- hair_eye_col$Eye
# Hair/Eye color table
he.tbl <- table(Hair, Eye)
he.tbl

       Eye
Hair    Brown Blue Hazel Green
  Black    68   20    15     5
  Brown   119   84    54    29
  Red      26   17    14    14
  Blond     7   94    10    16

# Margin of rows
margin.table(he.tbl, 1)

Hair
Black Brown   Red Blond 
  108   286    71   127 

# Margin of columns
margin.table(he.tbl, 2)

Eye
Brown  Blue Hazel Green 
  220   215    93    64 

3. Compute relative frequencies:

# Frequencies relative to row total
prop.table(he.tbl, 1)

       Eye
Hair         Brown       Blue      Hazel      Green
  Black 0.62962963 0.18518519 0.13888889 0.04629630
  Brown 0.41608392 0.29370629 0.18881119 0.10139860
  Red   0.36619718 0.23943662 0.19718310 0.19718310
  Blond 0.05511811 0.74015748 0.07874016 0.12598425

# Table of percentages
round(prop.table(he.tbl, 1), 2)*100

       Eye
Hair    Brown Blue Hazel Green
  Black    63   19    14     5
  Brown    42   29    19    10
  Red      37   24    20    20
  Blond     6   74     8    13

To express the frequencies relative to the grand total, use this:

he.tbl/sum(he.tbl)