Descriptive Statistics

Libraries

library(tidyverse)
library(ggpubr)
library(gridExtra)
library(forecast)
library(scales)
library(fmsb)
# library(sn)
# library(modeest)
# library(fpp)

Intro to Descriptive Statistics

# Create vector classA and classB
classA <- c(102,128,131,98,140,93,110,115,109,89,106,119,97)
classB <- c(127,131,96,80,93,120,109,162,103,111,109,87,105)

# Calculate the mean for each class
mean(classA)

[1] 110.5385

mean(classB)

[1] 110.2308

Types of Data

Summarizing Data

Mean

# See-saw illustration
see_saw <- c(93,106,131)
mean(see_saw)

[1] 110

# Density plot
set.seed(123)
income <- rnorm(50, mean = 50, sd = 5)
p1 <- ggplot(data = data.frame(income), aes(income)) + 
        geom_density(fill = "red", alpha = 0.1) +
        geom_vline(xintercept = mean(income), 
                   linetype="dotted", 
                   color = "red", size=1.5)
p1

mean(income)

[1] 50.17202

# Outlier affects mean
income_outlier <- c(income,100,100,100,100,100,100)
p2 <- ggplot(data = data.frame(income = income_outlier), aes(income)) + 
      geom_density(fill = "red", alpha = 0.1) +
        geom_vline(xintercept = mean(income_outlier), linetype="dotted", 
                   color = "red", size=1.5)
p2

mean(income_outlier)

[1] 55.51073

# Normal data
set.seed(123)
sample_normal <- rnorm(n = 10000, mean = 0, sd = 1)

# Skew Normal data
set.seed(123)
sample_sn <- sn::rsn(n=10000, xi=0, omega=1, alpha=10, tau=0)

# Create density plots
p3 <- ggplot(data = data.frame(x = sample_normal), aes(x=x)) + 
        geom_density(fill = "red", alpha = 0.1) +
        geom_vline(xintercept = mean(sample_normal), linetype="dotted", 
                   color = "red", size=1) +
        ggtitle("Normal Distribution")
      
p4 <- ggplot(data = data.frame(x = sample_sn), aes(x=x)) + 
        geom_density(fill = "red", alpha = 0.1) +
        geom_vline(xintercept = mean(sample_sn), linetype="dotted", 
                   color = "red", size=1) +
        ggtitle("Skew Normal Distribution")

grid.arrange(p3,p4, ncol=1)

Median

grid.arrange(p1 + geom_vline(xintercept = median(income), linetype="dotted", 
                             color = "blue", size=1.5),
             p2 + geom_vline(xintercept = median(income_outlier), 
                             linetype="dotted", color = "blue", size=1.5))

summary(income)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  40.17   47.20   49.64   50.17   53.49   60.84 

summary(income_outlier)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  40.17   47.53   50.60   55.51   54.41  100.00 

grid.arrange(p3 + geom_vline(xintercept = median(sample_normal), 
                             linetype="dotted", color = "blue", size=1),
             p4 + geom_vline(xintercept = median(sample_sn), 
                             linetype="dotted", color = "blue", size=1))

summary(sample_normal)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-3.845320 -0.667969 -0.011089 -0.002372  0.673347  3.847768 

# Quartile examples
sample_quart <- quantile(sample_normal,probs = c(0.25, 0.5, 0.75))
sample_quart

        25%         50%         75% 
-0.66796903 -0.01108922  0.67334733 

# Percentile examples
sample_percent <- quantile(sample_normal,probs = seq(0.1,0.9,by=0.1))
sample_percent

        10%         20%         30%         40%         50%         60%         70%         80%         90% 
-1.27943893 -0.83997128 -0.51872750 -0.24831290 -0.01108922  0.24768364  0.51856015  0.84214001  1.26280398 

# Density plot with quartiles
p5 <- ggplot(data = data.frame(x = sample_normal), aes(x=x)) + 
        geom_density(fill = "red", alpha = 0.1) +
        geom_vline(xintercept = sample_quart, linetype="dotted", 
                   color = "red", size=1) +
        ggtitle("Normal Distribution")
p5

Mode

classAB <- c(classA, classB)
modeest::mfv(classAB)

[1] 109

Range

# Calculate range
range(classA)

[1]  89 140

range(classB)

[1]  80 162

range(classA)[2] - range(classA)[1]

[1] 51

range(classB)[2] - range(classB)[1]

[1] 82

Interquartile Range

# Calculate quartiles
income_outlier_quart <- quantile(income_outlier,probs = c(0.25, 0.5, 0.75),
                                 names = FALSE)

# Density plot with quartiles
p6 <- ggplot(data = data.frame(income = income_outlier), aes(income)) + 
      geom_density(fill = "red", alpha = 0.1) +
        geom_vline(xintercept = income_outlier_quart, linetype="dotted", 
                   color = "red", size=1)
p6

# Range of income_outlier
range(income_outlier)[2] - range(income_outlier)[1]

[1] 59.83309

# Interquartile range of income_outlier
income_outlier_quart[3] - income_outlier_quart[1]

[1] 6.879677

Variance and Standard Deviation

set.seed(123)
sample_multinormal <- data.frame(type = as.factor(c(rep(4:6,each=10000))),
                                 data = c(rnorm(10000,0,4),
                                          rnorm(10000,0,5),
                                          rnorm(10000,0,6)))

p7 <- ggplot(data = sample_multinormal, aes(x = data, fill = type, colour = type)) +
        geom_density(alpha = 0.2)
p7

# Comparing MLE variance (denominator n) and and unbiased variance (denominator n-1)

# n = 10
set.seed(123)
sample1 <- rnorm(10, mean = 0, sd = 2)

# Unbiased variance
sd(sample1)^2

[1] 3.638816

# MLE variance
sd(sample1)^2 * (length(sample1)-1)/length(sample1)

[1] 3.274934

# n = 30
set.seed(123)
sample2 <- rnorm(30, mean = 0, sd = 2)

# Unbiased variance
sd(sample2)^2

[1] 3.849685

# MLE variance
sd(sample2)^2 * (length(sample2)-1)/length(sample2)

[1] 3.721362

# n = 100
set.seed(123)
sample3 <- rnorm(100, mean = 0, sd = 2)

# Unbiased variance
sd(sample3)^2

[1] 3.332931

# MLE variance
sd(sample3)^2 * (length(sample3)-1)/length(sample3)

[1] 3.299602

# n = 1000
set.seed(123)
sample4 <- rnorm(1000, mean = 0, sd = 2)

# Unbiased variance
sd(sample4)^2

[1] 3.933836

# MLE variance
sd(sample4)^2 * (length(sample4)-1)/length(sample4)

[1] 3.929902

Graphs

Pie chart

df_pie <- data.frame(
            group = c("Male", "Female", "Child"),
            value = c(20, 30, 50)
            )
df_pie

pie <- ggplot(df_pie, aes(x="", y=value, fill=group)) +
        geom_bar(width = 1, stat = "identity") + 
        coord_polar("y", start=0) +
        # scale_fill_brewer(palette="Blues") + 
        theme_minimal() +
        theme(
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.border = element_blank(),
        panel.grid=element_blank(),
        axis.ticks = element_blank(),
        plot.title=element_text(size=14, face="bold")
        ) +
        theme(axis.text.x=element_blank()) +
        geom_text(aes(y = value/3 + c(0, cumsum(value)[-length(value)]), 
                      label = scales::percent(value/100)), size=5)
pie

Bar chart

df_vehicles <- data.frame(kendaraan=c("Sepeda Motor", "Mobil", "Bus"),
                          jumlah=c(30, 10, 5))
df_vehicles

bar <- ggplot(data=df_vehicles, aes(x=kendaraan, y=jumlah)) + 
        geom_bar(stat="identity", fill="steelblue") +
        geom_text(aes(label=jumlah), vjust=1.6, color="white", size=5)
bar

Line chart

tourist <- fpp::austourists
autoplot(tourist, color = 'black', linetype = 'solid') + geom_point() +
      ggtitle("Quarterly visitor nights spent by International Tourists to Australia")

Histogram

histogram <- ggplot(data = data.frame(x = sample_normal), aes(x=x)) + 
              geom_histogram(bins=30, color="darkblue", fill="lightblue") +
              ggtitle("Normal Distribution")
histogram

# Bimodal distribution
sample_bimodal <- c(rnorm(n = 10000, mean = 0, sd = 1), 
                    rnorm(n = 10000, mean = 4, sd = 1.2))
histogram2 <- ggplot(data = data.frame(x = sample_bimodal), aes(x=x)) + 
              geom_histogram(bins=30, color="black", fill="white") +
              ggtitle("Bimodal Distribution")
histogram2

Radar chart

# Demo data
exam_scores <- data.frame(
    row.names = c("Student 1", "Student 2", "Student 3"),
      Biology = c(7.9, 3.9, 9.4),
      Physics = c(10, 20, 0),
        Maths = c(3.7, 11.5, 2.5),
        Sport = c(8.7, 20, 4),
      English = c(7.9, 7.2, 12.4),
    Geography = c(6.4, 10.5, 6.5),
          Art = c(2.4, 0.2, 9.8),
  Programming = c(0, 0, 20),
        Music = c(20, 20, 20)
)

# Define the variable ranges: maximum and minimum
max_min <- data.frame(
  Biology = c(20, 0), Physics = c(20, 0), Maths = c(20, 0),
  Sport = c(20, 0), English = c(20, 0), Geography = c(20, 0),
  Art = c(20, 0), Programming = c(20, 0), Music = c(20, 0)
)
rownames(max_min) <- c("Max", "Min")

# Bind the variable ranges to the data
df_exam <- rbind(max_min, exam_scores)
df_exam

student1_data <- df_exam[c("Max", "Min", "Student 1"), ]
radarchart(student1_data)

create_beautiful_radarchart <- function(data, color = "#00AFBB", 
                                        vlabels = colnames(data), vlcex = 0.7,
                                        caxislabels = NULL, title = NULL, ...){
  radarchart(
    data, axistype = 1,
    # Customize the polygon
    pcol = color, pfcol = scales::alpha(color, 0.5), plwd = 2, plty = 1,
    # Customize the grid
    cglcol = "grey", cglty = 1, cglwd = 0.8,
    # Customize the axis
    axislabcol = "grey", 
    # Variable labels
    vlcex = vlcex, vlabels = vlabels,
    caxislabels = caxislabels, title = title, ...
  )
}

# Define colors and titles
colors <- c("#00AFBB", "#E7B800", "#FC4E07")
titles <- c("Student 1", "Student 2", "Student 3")

# Reduce plot margin using par()
# Split the screen in 3 parts
par(mfrow = c(1,3))

# Create the radar chart
for(i in 1:3){
  create_beautiful_radarchart(
    data = df_exam[c(1, 2, i+2), ], caxislabels = c(0, 5, 10, 15, 20),
    color = colors[i], title = titles[i]
    )
}

Bubble chart

See https://www.gapminder.org/tools/#$chart-type=bubbles&url=v1

# Load data
df_cars <- mtcars

# Convert cyl as a grouping variable
df_cars$cyl <- as.factor(df_cars$cyl)

# Bubble cart
plot_cars <- ggplot(df_cars, aes(x = wt, y = mpg)) + 
              geom_point(aes(color = cyl, size = qsec), alpha = 0.5) +
              scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07")) +
              scale_size(range = c(0.5, 12))  # Adjust the range of points size

plot_cars

Scatter plot

# Scatter plot
plot_cars2 <- ggplot(df_cars, aes(x = hp, y = mpg)) + 
              geom_point(size = 2)

plot_cars2

Dot plot

homework <- data.frame(hour = seq(0,10),
                       freq = c(1,5,4,3,4,2,3,0,2,1,0))

hours_spent <- rep(homework$hour[1],homework$freq[1])
for (i in 2:11){
  hours_spent <- c(hours_spent, rep(homework$hour[i],homework$freq[i]))
}

hours_spent

 [1] 0 1 1 1 1 1 2 2 2 2 3 3 3 4 4 4 4 5 5 6 6 6 8 8 9

plot_hour <- ggplot(data = data.frame(hour = hours_spent), aes(x=hour)) +
              geom_dotplot(binwidth = 1, dotsize=0.4) +
              scale_x_continuous(breaks=seq(0,10))
plot_hour

Boxplot

df_plant <- PlantGrowth

box1 <- ggboxplot(df_plant, x = "group", y = "weight", color = "group",
          ylab = "Weight", xlab = "Treatment")

box2 <- ggline(df_plant, x = "group", y = "weight", 
          add = c("mean_se", "jitter"), 
          order = c("ctrl", "trt1", "trt2"),
          ylab = "Weight", xlab = "Treatment")

grid.arrange(box1, box2, ncol=2)