Data Visualization and Descriptive Statistics in R

Hands-on Training Using the Diarrhea Dataset

1 Introduction

This training module develops practical skills in data inspection, data cleaning, descriptive statistics, exploratory data analysis (EDA), and data visualization in R.

The examples use the diarrhea dataset available from:

https://raw.githubusercontent.com/bijayprad/data/refs/heads/main/dirrhea.csv

Dataset note: The training examples use the variables exactly as represented in the source CSV. Always verify definitions, coding, and units against the original data documentation before using the dataset for substantive research.

The purpose is not only to produce graphs and numerical summaries. Trainees should learn to:

1. understand the structure and meaning of variables;
2. identify missing, inconsistent, and unusual values;
3. choose an appropriate statistical summary;
4. choose an appropriate graph;
5. compare groups carefully;
6. describe relationships without overstating them;
7. communicate results in a research-oriented way.

2 Learning Objectives

By the end of this module, trainee should be able to:

• import a CSV file directly from a URL;
• inspect the dimensions and structure of a dataset;
• identify numerical and categorical variables;
• detect missing and inconsistent values;
• convert variables to appropriate R classes;
• calculate measures of central tendency and dispersion;
• create frequency and contingency tables;
• produce publication-ready exploratory graphs with ggplot2;
• compare numerical variables across groups;
• visualize relationships between numerical variables;
• calculate and interpret correlation appropriately;
• identify potential outliers;
• distinguish an unusual observation from a data-entry error;
• write concise research-style interpretations.

3 Packages and Setup

Install packages once if necessary:

install.packages(c(
  "tidyverse",
  "janitor",
  "skimr",
  "scales",
  "patchwork"
))

Load the packages:

library(tidyverse)
library(janitor)
library(skimr)
library(scales)
library(patchwork)

theme_set(theme_classic())

3.0.1 Importing the Diarrhea Dataset

3.1 Read the CSV directly from GitHub

url <- "https://raw.githubusercontent.com/bijayprad/data/refs/heads/main/dirrhea.csv"

dirdata <- read_csv(
  url,
  na = c("", "NA")
)

## Rows: 1296 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Department, Fever, DegreeDehy, OutComeDisc
## dbl (5): ID, Sex, Age, DiaMaxNo_24hrs, DiaDurInDays
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

dirdata

3.2 First inspection

head(dirdata)

tail(dirdata)

dim(dirdata)

## [1] 1296    9

names(dirdata)

## [1] "ID"             "Department"     "Sex"            "Age"           
## [5] "Fever"          "DiaMaxNo_24hrs" "DiaDurInDays"   "DegreeDehy"    
## [9] "OutComeDisc"

The current CSV contains the following variables:

Variable	Description for training	Suggested type
ID	Record identifier	Integer
Department	Department/category such as OPD or IPD	Factor
Sex	Sex coded as 1/2 in the source	Factor
Age	Age value recorded in the source	Numeric
Fever	Fever status	Factor
DiaMaxNo_24hrs	Maximum number of diarrhea episodes in 24 hours	Integer
DiaDurInDays	Duration of diarrhea	Numeric
DegreeDehy	Degree of dehydration	Factor
OutComeDisc	Recorded outcome/disposition	Factor

Data documentation note: The CSV records Age as a numerical variable and Sex as codes 1 and 2. The unit of Age is not encoded in the CSV itself, so the original data documentation should be checked before describing it as months, years, or another unit.

4 Data Audit Before Analysis

A good statistical analysis starts with an audit. Do not immediately make graphs from an unfamiliar dataset.

4.1 Structure

str(dirdata)

## spc_tbl_ [1,296 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ ID            : num [1:1296] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Department    : chr [1:1296] "OPD" "OPD" "IPD" "IPD" ...
##  $ Sex           : num [1:1296] 1 2 1 1 1 2 2 1 2 1 ...
##  $ Age           : num [1:1296] 6 8 12 2 12 12 10 21 6 37 ...
##  $ Fever         : chr [1:1296] "No" "No" "No" "Yes" ...
##  $ DiaMaxNo_24hrs: num [1:1296] 7 5 10 9 8 7 12 8 8 10 ...
##  $ DiaDurInDays  : num [1:1296] 2 1 3 1 2 5 2 2 2 2 ...
##  $ DegreeDehy    : chr [1:1296] "some" "some" "some" "some" ...
##  $ OutComeDisc   : chr [1:1296] NA NA "Improved" "Improved" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   ID = col_double(),
##   ..   Department = col_character(),
##   ..   Sex = col_double(),
##   ..   Age = col_double(),
##   ..   Fever = col_character(),
##   ..   DiaMaxNo_24hrs = col_double(),
##   ..   DiaDurInDays = col_double(),
##   ..   DegreeDehy = col_character(),
##   ..   OutComeDisc = col_character()
##   .. )
##  - attr(*, "problems")=<pointer: 0x000001db41464fb0>

glimpse(dirdata)

## Rows: 1,296
## Columns: 9
## $ ID             <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ Department     <chr> "OPD", "OPD", "IPD", "IPD", "IPD", "OPD", "IPD", "IPD",…
## $ Sex            <dbl> 1, 2, 1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 2, 1, 1, 1…
## $ Age            <dbl> 6, 8, 12, 2, 12, 12, 10, 21, 6, 37, 12, 3, 23, 7, 6, 4,…
## $ Fever          <chr> "No", "No", "No", "Yes", "Yes", "No", "No", "No", "No",…
## $ DiaMaxNo_24hrs <dbl> 7, 5, 10, 9, 8, 7, 12, 8, 8, 10, 7, 9, 3, 11, 10, 6, 12…
## $ DiaDurInDays   <dbl> 2, 1, 3, 1, 2, 5, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 3, 2, 2…
## $ DegreeDehy     <chr> "some", "some", "some", "some", "some", "some", "some",…
## $ OutComeDisc    <chr> NA, NA, "Improved", "Improved", "Improved", NA, "Improv…

4.2 Number of observations and variables

nrow(dirdata)

## [1] 1296

ncol(dirdata)

## [1] 9

4.3 Summary

summary(dirdata)

##        ID             Department        Sex             Age      
##  Min.   :   1.0   Length   :1296   Min.   :1.000   Min.   : 1.0  
##  1st Qu.: 324.8   N.unique :   2   1st Qu.:1.000   1st Qu.: 5.0  
##  Median : 648.5   N.blank  :   0   Median :2.000   Median : 9.0  
##  Mean   : 648.5   Min.nchar:   3   Mean   :1.667   Mean   :11.2  
##  3rd Qu.: 972.2   Max.nchar:   3   3rd Qu.:2.000   3rd Qu.:12.0  
##  Max.   :1296.0                    Max.   :2.000   Max.   :59.0  
##                                                                  
##        Fever      DiaMaxNo_24hrs    DiaDurInDays        DegreeDehy  
##  Length   :1296   Min.   : 2.000   Min.   : 1.000   Length   :1296  
##  N.unique :   2   1st Qu.: 7.000   1st Qu.: 1.000   N.unique :   4  
##  N.blank  :   0   Median : 9.000   Median : 2.000   N.blank  :   0  
##  Min.nchar:   2   Mean   : 9.836   Mean   : 2.316   Min.nchar:   4  
##  Max.nchar:   3   3rd Qu.:12.000   3rd Qu.: 3.000   Max.nchar:   6  
##  NAs      :   1   Max.   :35.000   Max.   :52.000                   
##                                    NAs    :2                        
##     OutComeDisc  
##  Length   :1296  
##  N.unique :   1  
##  N.blank  :   0  
##  Min.nchar:   8  
##  Max.nchar:   8  
##  NAs      : 196  
## 

4.4 Variable classes

sapply(dirdata, class)

##             ID     Department            Sex            Age          Fever 
##      "numeric"    "character"      "numeric"      "numeric"    "character" 
## DiaMaxNo_24hrs   DiaDurInDays     DegreeDehy    OutComeDisc 
##      "numeric"      "numeric"    "character"    "character"

4.5 Unique values of categorical variables

lapply(dirdata[c("Department", "Sex", "Fever", "DegreeDehy", "OutComeDisc")], unique)

## $Department
## [1] "OPD" "IPD"
## 
## $Sex
## [1] 1 2
## 
## $Fever
## [1] "No"  "Yes" NA   
## 
## $DegreeDehy
## [1] "some"   "Severe" "severe" "Some"  
## 
## $OutComeDisc
## [1] NA         "Improved"

4.6 Missing values

colSums(is.na(dirdata))

##             ID     Department            Sex            Age          Fever 
##              0              0              0              0              1 
## DiaMaxNo_24hrs   DiaDurInDays     DegreeDehy    OutComeDisc 
##              0              2              0            196

5 Data Cleaning and Preparation

Cleaning should be transparent and reproducible. Do not delete observations simply because they look unusual.

5.1 Standardize text

The source contains inconsistent capitalization in DegreeDehy, such as some and Some, and severe and Severe. We can standardize these values without changing their meaning.

dirdata <- dirdata %>%
  mutate(
    Department = factor(str_to_upper(str_trim(Department))),
    Sex = factor(Sex, levels = c(1, 2), labels = c("Male", "Female")),
    Fever = factor(str_to_title(str_trim(Fever))),
    DegreeDehy = factor(str_to_title(str_trim(DegreeDehy))),
    OutComeDisc = factor(str_trim(OutComeDisc)),
    ID = as.integer(ID)
  )

glimpse(dirdata)

## Rows: 1,296
## Columns: 9
## $ ID             <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ Department     <fct> OPD, OPD, IPD, IPD, IPD, OPD, IPD, IPD, IPD, IPD, IPD, …
## $ Sex            <fct> Male, Female, Male, Male, Male, Female, Female, Male, F…
## $ Age            <dbl> 6, 8, 12, 2, 12, 12, 10, 21, 6, 37, 12, 3, 23, 7, 6, 4,…
## $ Fever          <fct> No, No, No, Yes, Yes, No, No, No, No, No, No, No, No, N…
## $ DiaMaxNo_24hrs <dbl> 7, 5, 10, 9, 8, 7, 12, 8, 8, 10, 7, 9, 3, 11, 10, 6, 12…
## $ DiaDurInDays   <dbl> 2, 1, 3, 1, 2, 5, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 3, 2, 2…
## $ DegreeDehy     <fct> Some, Some, Some, Some, Some, Some, Some, Some, Some, S…
## $ OutComeDisc    <fct> NA, NA, Improved, Improved, Improved, NA, Improved, Imp…

5.2 Check duplicates

sum(duplicated(dirdata))

## [1] 0

Check duplicated IDs:

dirdata %>%
  count(ID) %>%
  filter(n > 1)

An observation that looks duplicated should be investigated against the original records. Never remove duplicates automatically.

5.3 Check numerical ranges

dirdata %>%
  summarise(
    Age_min = min(Age, na.rm = TRUE),
    Age_max = max(Age, na.rm = TRUE),
    DiaMax_min = min(DiaMaxNo_24hrs, na.rm = TRUE),
    DiaMax_max = max(DiaMaxNo_24hrs, na.rm = TRUE),
    DiaDur_min = min(DiaDurInDays, na.rm = TRUE),
    DiaDur_max = max(DiaDurInDays, na.rm = TRUE)
  )

5.4 Investigate unusual values

For example:

dirdata %>%
  filter(DiaDurInDays > quantile(DiaDurInDays, 0.99, na.rm = TRUE)) %>%
  arrange(desc(DiaDurInDays))

The purpose of this step is investigation, not automatic deletion.

6 Understanding Variable Types

6.1 Categorical variables

Examples in this dataset include:

• Department
• Sex
• Fever
• DegreeDehy
• OutComeDisc

6.2 Numerical variables

Examples include:

• Age
• DiaMaxNo_24hrs
• DiaDurInDays

DiaMaxNo_24hrs is a count and is therefore discrete in its measurement scale, while Age and DiaDurInDays are numerical measurements as recorded in the dataset.

6.3 Choosing a graph

Question	Variable structure	Suitable graph
How common is each category?	One categorical	Bar chart
How is a numerical variable distributed?	One numerical	Histogram/density
Are there potential outliers?	One numerical	Boxplot
How does a numerical variable differ by group?	Categorical + numerical	Boxplot/violin plot
Are two numerical variables associated?	Two numerical	Scatterplot
How do two categorical variables compare?	Two categorical	Grouped/proportional bar chart
How does a distribution differ by group?	Numerical + categorical	Faceted histogram/density

7 Descriptive Statistics

7.1 Mean

The arithmetic mean is

\[\bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_i.\]

mean(dirdata$Age, na.rm = TRUE)

## [1] 11.2037

7.2 Median

median(dirdata$Age, na.rm = TRUE)

## [1] 9

7.3 Minimum, maximum and range

min(dirdata$Age, na.rm = TRUE)

## [1] 1

max(dirdata$Age, na.rm = TRUE)

## [1] 59

range(dirdata$Age, na.rm = TRUE)

## [1]  1 59

diff(range(dirdata$Age, na.rm = TRUE))

## [1] 58

7.4 Variance and standard deviation

var(dirdata$Age, na.rm = TRUE)

## [1] 97.14921

sd(dirdata$Age, na.rm = TRUE)

## [1] 9.85643

7.5 Quartiles and IQR

quantile(
  dirdata$Age,
  probs = c(0.25, 0.50, 0.75),
  na.rm = TRUE
)

## 25% 50% 75% 
##   5   9  12

IQR(dirdata$Age, na.rm = TRUE)

## [1] 7

7.6 Five-number summary

fivenum(dirdata$Age)

## [1]  1  5  9 12 59

7.7 A reusable descriptive-statistics function (Generate your own function)

describe_numeric <- function(x) {
  tibble(
    n = sum(!is.na(x)),
    missing = sum(is.na(x)),
    mean = mean(x, na.rm = TRUE),
    median = median(x, na.rm = TRUE),
    sd = sd(x, na.rm = TRUE),
    min = min(x, na.rm = TRUE),
    q1 = quantile(x, 0.25, na.rm = TRUE),
    q3 = quantile(x, 0.75, na.rm = TRUE),
    max = max(x, na.rm = TRUE),
    IQR = IQR(x, na.rm = TRUE)
  )
}

describe_numeric(dirdata$Age)

7.8 Compare several numerical variables

dirdata %>%
  summarise(
    across(
      c(Age, DiaMaxNo_24hrs, DiaDurInDays),
      list(
        n = ~ sum(!is.na(.)),
        mean = ~ mean(.x, na.rm = TRUE),
        median = ~ median(.x, na.rm = TRUE),
        sd = ~ sd(.x, na.rm = TRUE),
        IQR = ~ IQR(.x, na.rm = TRUE)
      )
    )
  )

8 Frequency Tables for Categorical Variables

8.1 Frequency

table(dirdata$Department)

## 
##  IPD  OPD 
## 1017  279

table(dirdata$Sex)

## 
##   Male Female 
##    432    864

table(dirdata$Fever)

## 
##  No Yes 
## 869 426

table(dirdata$DegreeDehy)

## 
## Severe   Some 
##    238   1058

table(dirdata$OutComeDisc)

## 
## Improved 
##     1100

8.2 Percentage

prop.table(table(dirdata$Department)) * 100

## 
##      IPD      OPD 
## 78.47222 21.52778

A tidy version:

dirdata %>%
  count(Department) %>%
  mutate(Percent = 100 * n / sum(n))

8.3 Reusable frequency table

frequency_table <- function(data, variable) {
  data %>%
    count({{ variable }}) %>%
    mutate(Percent = 100 * n / sum(n))
}

frequency_table(dirdata, DegreeDehy)

9 Cross-Tabulation of Categorical Variables

9.1 Two-way table

table(dirdata$Department, dirdata$Sex)

##      
##       Male Female
##   IPD  341    676
##   OPD   91    188

9.2 Row percentages

prop.table(
  table(dirdata$Department, dirdata$Sex),
  margin = 1
) * 100

##      
##           Male   Female
##   IPD 33.52999 66.47001
##   OPD 32.61649 67.38351

9.3 Column percentages

prop.table(
  table(dirdata$Department, dirdata$Sex),
  margin = 2
) * 100

##      
##           Male   Female
##   IPD 78.93519 78.24074
##   OPD 21.06481 21.75926

10 Visualizing of Categorical Variables

10.1 Bar Diagram

barplot(table(dirdata$Department),
        main = "Records by Department",
        xlab = "Department",
        ylab = "Frequency",
        col = "skyblue")

10.2 Horizontal Bar Diagram

barplot(table(dirdata$Department),
        horiz = TRUE,
        main = "Records by Department",
        xlab = "Frequency",
        col = "lightgreen",
        las = 1)

10.3 Percentage Bar Diagram

dept <- prop.table(table(dirdata$Department)) * 100

barplot(dept,
        main = "Percentage by Department",
        ylab = "Percentage",
        col = "lightblue")

10.4 Grouped Bar Diagram

first we need a cross table

table(dirdata$Department, dirdata$Sex)

##      
##       Male Female
##   IPD  341    676
##   OPD   91    188

Bar diagram

tab <- table(dirdata$Department, dirdata$Sex)

barplot(tab,
        beside = TRUE,
        col = c("skyblue", "pink"),
        main = "Department by Sex",
        xlab = "Department",
        ylab = "Frequency",
        legend.text = colnames(tab))

10.5 Pie Diagram

pie(table(dirdata$Department),
    main = "Distribution by Department",
    col = rainbow(length(table(dirdata$Department))))

10.6 Pie Diagram with Percentage

dept <- table(dirdata$Department)

percent <- round(100 * dept / sum(dept), 1)

pie(dept,
    labels = paste(names(dept), percent, "%"),
    main = "Department Distribution",
    col = rainbow(length(dept)))

11 Visualizing of Numerical Variables

11.1 Scatter Diagram

plot(dirdata$Age,
     dirdata$DiaDurInDays,
     main = "Age and Diarrhea Duration",
     xlab = "Age",
     ylab = "Duration",
     pch = 16,
     col = "blue")

11.2 Scatter Diagram with Regression Line

plot(dirdata$Age,
     dirdata$DiaDurInDays,
     main = "Age and Diarrhea Duration",
     xlab = "Age",
     ylab = "Duration",
     pch = 16,
     col = "gray")

model <- lm(DiaDurInDays ~ Age, data = dirdata)

abline(model,
       col = "red",
       lwd = 2)

11.3 Histogram

hist(dirdata$Age,
     breaks = 20,
     main = "Distribution of Age",
     xlab = "Age",
     ylab = "Frequency",
     col = "lightblue")

11.4 Choosing the Number of Bins

par(mfrow = c(1, 3))

hist(dirdata$Age,
     breaks = 10,
     main = "10 Bins",
     col = "skyblue")

hist(dirdata$Age,
     breaks = 20,
     main = "20 Bins",
     col = "lightgreen")

hist(dirdata$Age,
     breaks = 30,
     main = "30 Bins",
     col = "pink")

par(mfrow = c(1, 1))

11.5 Histogram with Mean and Median

hist(dirdata$Age,
     breaks = 20,
     col = "lightblue",
     main = "Age with Mean and Median",
     xlab = "Age")

abline(v = mean(dirdata$Age, na.rm = TRUE),
       col = "red",
       lwd = 2,
       lty = 2)

abline(v = median(dirdata$Age, na.rm = TRUE),
       col = "blue",
       lwd = 2,
       lty = 3)

legend("topright",
       legend = c("Mean", "Median"),
       col = c("red", "blue"),
       lty = c(2, 3),
       lwd = 2)

11.6 Boxplot

boxplot(dirdata$Age,
        main = "Boxplot of Age",
        ylab = "Age",
        col = "lightgreen")

11.7 Boxplot by Department

boxplot(Age ~ Department,
        data = dirdata,
        main = "Age by Department",
        xlab = "Department",
        ylab = "Age",
        col = "lightblue")

11.8 Diarrhea Duration by Department

boxplot(DiaDurInDays ~ Department,
        data = dirdata,
        main = "Diarrhea Duration by Department",
        xlab = "Department",
        ylab = "Duration (Days)",
        col = "lightpink")

11.9 Potential Outliers

boxplot.stats(dirdata$Age)$out

##   [1] 37 23 23 32 24 48 24 48 27 24 49 24 24 48 24 36 32 24 28 24 24 59 48 36 24
##  [26] 24 48 24 48 36 24 48 24 48 24 24 36 31 25 28 56 24 42 39 36 24 24 24 24 24
##  [51] 36 48 23 24 36 23 26 36 25 36 36 48 36 45 24 24 24 24 24 24 24 48 48 24 48
##  [76] 24 24 24 48 24 56 48 24 36 48 48 36 36 24 24 24 24 24 36 24 24 36 24 36 27
## [101] 24 24 36 23 24 31 23 49 36 42 24 48 24 29 30 30 24 30 36 48 30 29 24 24 29
## [126] 24 24 56 59 24 37 59 48 36 36 36 23 29 36 48 36 30 24 42 36 42 52 24 30 30
## [151] 24 24 59 36 24 36 24 36 36

11.10 Density Plot

plot(density(dirdata$Age, na.rm = TRUE),
     main = "Density Distribution of Age",
     xlab = "Age",
     ylab = "Density",
     col = "blue",
     lwd = 2)

11.11 Histogram with Density Curve

hist(dirdata$Age,
     probability = TRUE,
     breaks = 20,
     col = "lightblue",
     main = "Histogram with Density Curve",
     xlab = "Age")

lines(density(dirdata$Age, na.rm = TRUE),
      col = "red",
      lwd = 2)

11.12 Line Diagram

plot(dirdata$Age,
     type = "l",
     main = "Age Pattern",
     xlab = "Observation",
     ylab = "Age",
     col = "blue",
     lwd = 2)

11.13 Line Diagram with Points

plot(dirdata$Age,
     type = "b",
     main = "Age Pattern",
     xlab = "Observation",
     ylab = "Age",
     col = "red",
     pch = 16)

12 Correlation Analysis

12.1 Pearson Correlation

cor(dirdata$Age,
    dirdata$DiaDurInDays,
    use = "complete.obs")

## [1] 0.1011935

12.2 Correlation Test

cor.test(dirdata$Age,
         dirdata$DiaDurInDays,
         method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  dirdata$Age and dirdata$DiaDurInDays
## t = 3.6561, df = 1292, p-value = 0.0002664
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.04695768 0.15483435
## sample estimates:
##       cor 
## 0.1011935

12.3 Spearman Correlation

cor.test(dirdata$Age,
         dirdata$DiaDurInDays,
         method = "spearman",
         use = "complete.obs")

## Warning in cor.test.default(dirdata$Age, dirdata$DiaDurInDays, method =
## "spearman", : cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  dirdata$Age and dirdata$DiaDurInDays
## S = 313095263, p-value = 1.582e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1329879

12.4 Correlation Matrix

numeric_data <- dirdata[, c(
  "Age",
  "DiaMaxNo_24hrs",
  "DiaDurInDays"
)]

cor_matrix <- cor(numeric_data,
                  use = "pairwise.complete.obs")

cor_matrix

##                        Age DiaMaxNo_24hrs DiaDurInDays
## Age             1.00000000    -0.04324656   0.10119346
## DiaMaxNo_24hrs -0.04324656     1.00000000   0.06194965
## DiaDurInDays    0.10119346     0.06194965   1.00000000

12.5 Visualizing Correlation Matrix

image(1:ncol(cor_matrix),
      1:nrow(cor_matrix),
      cor_matrix,
      axes = FALSE,
      col = hcl.colors(20, "Blue-Red 3"),
      main = "Correlation Matrix")

axis(1,
     at = 1:ncol(cor_matrix),
     labels = colnames(cor_matrix))

axis(2,
     at = 1:nrow(cor_matrix),
     labels = rownames(cor_matrix))

13 Multiple Diagrams in One Window

par(mfrow = c(2, 2))

barplot(table(dirdata$Department),
        main = "Bar Diagram",
        col = "skyblue")

pie(table(dirdata$Sex),
    main = "Pie Diagram")

hist(dirdata$Age,
     main = "Histogram",
     col = "lightgreen")

boxplot(dirdata$Age,
        main = "Boxplot",
        col = "pink")

14 Diagrams using ggplot

ggplot is a function from the ggplot2 package in R used to create high-quality and attractive graphs and charts.

ggplot2 is one of the most widely used packages for data visualization in R.

The ggplot2 has not been designed in our course however the trainee can practice using this tutorials.

14.1 Histogram

ggplot(dirdata, aes(x = Age)) +
  geom_histogram(
    bins = 20,
    color = "white"
  ) +
  labs(
    title = "Distribution of Age",
    x = "Age",
    y = "Frequency"
  )

14.2 Choosing the number of bins

The number of bins affects the appearance of a histogram.

Try several values:

p10 <- ggplot(dirdata, aes(x = Age)) +
  geom_histogram(bins = 10, color = "white") +
  labs(title = "10 bins")

p20 <- ggplot(dirdata, aes(x = Age)) +
  geom_histogram(bins = 20, color = "white") +
  labs(title = "20 bins")

p30 <- ggplot(dirdata, aes(x = Age)) +
  geom_histogram(bins = 30, color = "white") +
  labs(title = "30 bins")

p10 + p20 + p30

14.3 Density plot

ggplot(dirdata, aes(x = Age)) +
  geom_density(na.rm = TRUE) +
  labs(
    title = "Density of Age",
    x = "Age",
    y = "Density"
  )

14.4 Histogram with density

ggplot(dirdata, aes(x = Age)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 20,
    color = "white"
  ) +
  geom_density(linewidth = 1) +
  labs(
    title = "Age Distribution with Density Curve",
    x = "Age",
    y = "Density"
  )

14.5 Mean and median

age_mean <- mean(dirdata$Age, na.rm = TRUE)
age_median <- median(dirdata$Age, na.rm = TRUE)

ggplot(dirdata, aes(x = Age)) +
  geom_histogram(bins = 20, color = "white") +
  geom_vline(
    xintercept = age_mean,
    linetype = "dashed",
    linewidth = 1
  ) +
  geom_vline(
    xintercept = age_median,
    linetype = "dotted",
    linewidth = 1
  ) +
  labs(
    title = "Age Distribution with Mean and Median",
    x = "Age",
    y = "Frequency"
  )

15 Boxplots and Potential Outliers

15.1 Boxplot

ggplot(dirdata, aes(y = Age)) +
  geom_boxplot() +
  labs(
    title = "Boxplot of Age",
    y = "Age",
    x = NULL
  )

15.2 Boxplot statistics

boxplot.stats(dirdata$Age)

## $stats
## [1]  1  5  9 12 22
## 
## $n
## [1] 1296
## 
## $conf
## [1] 8.692778 9.307222
## 
## $out
##   [1] 37 23 23 32 24 48 24 48 27 24 49 24 24 48 24 36 32 24 28 24 24 59 48 36 24
##  [26] 24 48 24 48 36 24 48 24 48 24 24 36 31 25 28 56 24 42 39 36 24 24 24 24 24
##  [51] 36 48 23 24 36 23 26 36 25 36 36 48 36 45 24 24 24 24 24 24 24 48 48 24 48
##  [76] 24 24 24 48 24 56 48 24 36 48 48 36 36 24 24 24 24 24 36 24 24 36 24 36 27
## [101] 24 24 36 23 24 31 23 49 36 42 24 48 24 29 30 30 24 30 36 48 30 29 24 24 29
## [126] 24 24 56 59 24 37 59 48 36 36 36 23 29 36 48 36 30 24 42 36 42 52 24 30 30
## [151] 24 24 59 36 24 36 24 36 36

15.3 Potential outliers

boxplot.stats(dirdata$Age)$out

##   [1] 37 23 23 32 24 48 24 48 27 24 49 24 24 48 24 36 32 24 28 24 24 59 48 36 24
##  [26] 24 48 24 48 36 24 48 24 48 24 24 36 31 25 28 56 24 42 39 36 24 24 24 24 24
##  [51] 36 48 23 24 36 23 26 36 25 36 36 48 36 45 24 24 24 24 24 24 24 48 48 24 48
##  [76] 24 24 24 48 24 56 48 24 36 48 48 36 36 24 24 24 24 24 36 24 24 36 24 36 27
## [101] 24 24 36 23 24 31 23 49 36 42 24 48 24 29 30 30 24 30 36 48 30 29 24 24 29
## [126] 24 24 56 59 24 37 59 48 36 36 36 23 29 36 48 36 30 24 42 36 42 52 24 30 30
## [151] 24 24 59 36 24 36 24 36 36

A boxplot outlier is an observation identified by a statistical rule. It is not automatically an error. Check the original record before deciding whether an observation should be corrected, retained, or excluded.

16 Comparing Groups

16.1 Age by Department

ggplot(dirdata, aes(x = Department, y = Age)) +
  geom_boxplot() +
  labs(
    title = "Age Distribution by Department",
    x = "Department",
    y = "Age"
  )

16.2 Add individual observations

ggplot(dirdata, aes(x = Department, y = Age)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(
    width = 0.12,
    alpha = 0.5
  ) +
  labs(
    title = "Age by Department with Individual Observations",
    x = "Department",
    y = "Age"
  )

16.3 Distribution of diarrhea duration by Department

ggplot(dirdata, aes(x = Department, y = DiaDurInDays)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.12, alpha = 0.4) +
  labs(
    title = "Diarrhea Duration by Department",
    x = "Department",
    y = "Duration"
  )

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

16.4 Dehydration by Department

ggplot(dirdata, aes(x = DegreeDehy, fill = Department)) +
  geom_bar(position = "dodge") +
  labs(
    title = "Degree of Dehydration by Department",
    x = "Degree of Dehydration",
    y = "Frequency",
    fill = "Department"
  )

17 Bar Charts for Categorical Variables

17.1 Basic bar chart

ggplot(dirdata, aes(x = Department)) +
  geom_bar() +
  labs(
    title = "Records by Department",
    x = "Department",
    y = "Frequency"
  )

17.2 Horizontal bar chart

ggplot(dirdata, aes(x = Department)) +
  geom_bar() +
  coord_flip() +
  labs(
    title = "Records by Department",
    x = NULL,
    y = "Frequency"
  )

18 Relationships Between Numerical Variables

18.1 Scatterplot

ggplot(
  dirdata,
  aes(x = Age, y = DiaDurInDays)
) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Age and Diarrhea Duration",
    x = "Age",
    y = "Duration"
  )

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

18.2 Add a fitted linear trend

ggplot(
  dirdata,
  aes(x = Age, y = DiaDurInDays)
) +
  geom_point(alpha = 0.5) +
  geom_smooth(
    method = "lm",
    se = TRUE
  ) +
  labs(
    title = "Age and Diarrhea Duration",
    x = "Age",
    y = "Duration"
  )

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

The regression line is a summary of the observed association. It does not establish causation.

18.3 Another relationship

ggplot(
  dirdata,
  aes(x = DiaMaxNo_24hrs, y = DiaDurInDays)
) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    title = "Maximum Episodes and Diarrhea Duration",
    x = "Maximum episodes in 24 hours",
    y = "Duration"
  )

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

19 Correlation

19.1 Pearson correlation

cor(
  dirdata$Age,
  dirdata$DiaDurInDays,
  use = "complete.obs"
)

## [1] 0.1011935

19.2 Correlation test

cor.test(
  dirdata$Age,
  dirdata$DiaDurInDays,
  use = "complete.obs",
  method = "pearson"
)

## 
##  Pearson's product-moment correlation
## 
## data:  dirdata$Age and dirdata$DiaDurInDays
## t = 3.6561, df = 1292, p-value = 0.0002664
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.04695768 0.15483435
## sample estimates:
##       cor 
## 0.1011935

19.3 Spearman correlation

For ordinal, non-normal, or strongly skewed variables, Spearman correlation may be useful:

cor.test(
  dirdata$Age,
  dirdata$DiaDurInDays,
  use = "complete.obs",
  method = "spearman"
)

## Warning in cor.test.default(dirdata$Age, dirdata$DiaDurInDays, use =
## "complete.obs", : cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  dirdata$Age and dirdata$DiaDurInDays
## S = 313095263, p-value = 1.582e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1329879

Correlation measures association, not causation. Also inspect a scatterplot before relying on a single correlation coefficient.

20 Correlation Matrix

numeric_data <- dirdata %>%
  select(
    Age,
    DiaMaxNo_24hrs,
    DiaDurInDays
  )

cor(
  numeric_data,
  use = "pairwise.complete.obs"
)

##                        Age DiaMaxNo_24hrs DiaDurInDays
## Age             1.00000000    -0.04324656   0.10119346
## DiaMaxNo_24hrs -0.04324656     1.00000000   0.06194965
## DiaDurInDays    0.10119346     0.06194965   1.00000000

A simple visual representation:

cor_matrix <- cor(
  numeric_data,
  use = "pairwise.complete.obs"
)

image(
  1:ncol(cor_matrix),
  1:nrow(cor_matrix),
  cor_matrix,
  axes = FALSE,
  xlab = "",
  ylab = "",
  col = hcl.colors(20, "Blue-Red 3")
)

axis(
  1,
  at = 1:ncol(cor_matrix),
  labels = colnames(cor_matrix)
)

axis(
  2,
  at = 1:nrow(cor_matrix),
  labels = rownames(cor_matrix)
)

21 Faceting

Faceting is useful when the same graph should be repeated across groups.

21.1 Age distributions by Department

ggplot(dirdata, aes(x = Age)) +
  geom_histogram(bins = 20, color = "white") +
  facet_wrap(~ Department) +
  labs(
    title = "Age Distribution by Department",
    x = "Age",
    y = "Frequency"
  )

21.2 Age distribution by Sex and Department

ggplot(dirdata, aes(x = Age)) +
  geom_histogram(bins = 15, color = "white") +
  facet_grid(Department ~ Sex) +
  labs(
    title = "Age Distribution by Department and Sex",
    x = "Age",
    y = "Frequency"
  )

22 Research-Oriented Summary Tables

22.1 Numerical variables

numeric_summary <- dirdata %>%
  summarise(
    across(
      c(Age, DiaMaxNo_24hrs, DiaDurInDays),
      list(
        N = ~ sum(!is.na(.)),
        Mean = ~ mean(.x, na.rm = TRUE),
        SD = ~ sd(.x, na.rm = TRUE),
        Median = ~ median(.x, na.rm = TRUE),
        IQR = ~ IQR(.x, na.rm = TRUE),
        Min = ~ min(.x, na.rm = TRUE),
        Max = ~ max(.x, na.rm = TRUE)
      )
    )
  )

numeric_summary

22.2 Summary by Department

dirdata %>%
  group_by(Department) %>%
  summarise(
    n = n(),
    Mean_Age = mean(Age, na.rm = TRUE),
    SD_Age = sd(Age, na.rm = TRUE),
    Median_Age = median(Age, na.rm = TRUE),
    IQR_Age = IQR(Age, na.rm = TRUE),
    Mean_Duration = mean(DiaDurInDays, na.rm = TRUE),
    Median_Duration = median(DiaDurInDays, na.rm = TRUE),
    .groups = "drop"
  )

23 Confidence Interval for a Mean

A confidence interval describes uncertainty around an estimated mean.

mean_ci <- function(x, conf = 0.95) {
  x <- x[!is.na(x)]
  n <- length(x)
  m <- mean(x)
  se <- sd(x) / sqrt(n)
  tcrit <- qt((1 + conf) / 2, df = n - 1)

  tibble(
    n = n,
    mean = m,
    lower = m - tcrit * se,
    upper = m + tcrit * se
  )
}

mean_ci(dirdata$Age)

By Department:

ci_department <- dirdata %>%
  group_by(Department) %>%
  summarise(
    mean_ci(Age),
    .groups = "drop"
  )

ci_department

Plot:

ggplot(
  ci_department,
  aes(x = Department, y = mean)
) +
  geom_point(size = 3) +
  geom_errorbar(
    aes(ymin = lower, ymax = upper),
    width = 0.15
  ) +
  labs(
    title = "Mean Age with 95% Confidence Intervals",
    x = "Department",
    y = "Mean Age"
  )

24 Optional Inferential Extensions

The main purpose of this module is descriptive analysis. If students continue to inferential statistics, the dataset can be used for introductory examples.

24.1 Chi-square test

Example research question:

Is Department associated with Sex?

tab <- table(dirdata$Department, dirdata$Sex)

chisq.test(tab)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab
## X-squared = 0.046246, df = 1, p-value = 0.8297

Before interpreting the result, inspect expected counts:

chisq.test(tab)$expected

##      
##       Male Female
##   IPD  339    678
##   OPD   93    186

24.2 Comparing a numerical variable between two departments

A t-test may be considered when its assumptions are reasonably supported:

t.test(
  Age ~ Department,
  data = dirdata
)

## 
##  Welch Two Sample t-test
## 
## data:  Age by Department
## t = 4.0155, df = 692.95, p-value = 6.581e-05
## alternative hypothesis: true difference in means between group IPD and group OPD is not equal to 0
## 95 percent confidence interval:
##  1.066338 3.106847
## sample estimates:
## mean in group IPD mean in group OPD 
##         11.652901          9.566308

A non-parametric alternative is:

wilcox.test(
  Age ~ Department,
  data = dirdata
)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Age by Department
## W = 151322, p-value = 0.08727
## alternative hypothesis: true location shift is not equal to 0

These tests should not be presented as automatic steps. Students should first examine the study design, measurement scale, distribution, independence, missing data, and assumptions.

25 Interpreting Graphs Correctly

A graph should lead to a description supported by the graph, not an unsupported conclusion.

25.1 Example: histogram

A useful interpretation discusses:

• centre;
• spread;
• shape;
• skewness;
• unusual observations.

Avoid:

“Age causes longer diarrhea duration.”

A scatterplot cannot support a causal statement.

Prefer:

“The scatterplot shows the observed association between age and diarrhea duration. The direction and strength should be assessed from the pattern of points and an appropriate correlation or model.”

25.2 Example: boxplot

A good interpretation may discuss:

• median;
• interquartile range;
• overall spread;
• potential outliers;
• differences between groups.

26 Publication-Quality Graphs

A research graph should have:

• a meaningful title or caption;
• clear axis labels;
• correct units where known;
• readable text;
• an appropriate scale;
• a useful legend;
• minimal unnecessary decoration;
• an appropriate graph type.

Example:

p <- ggplot(
  dirdata,
  aes(x = Department, y = DiaDurInDays)
) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(
    width = 0.1,
    alpha = 0.35
  ) +
  labs(
    title = "Distribution of Diarrhea Duration by Department",
    x = "Department",
    y = "Diarrhea duration"
  ) +
  theme_classic()

p

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

27 Saving Figures

ggsave(
  "diarrhea_duration_by_department.png",
  plot = p,
  width = 7,
  height = 5,
  units = "in",
  dpi = 300
)

For PDF:

ggsave(
  "diarrhea_duration_by_department.pdf",
  plot = p,
  width = 7,
  height = 5
)

28 Mini Analysis Workflow

The following workflow is recommended for a real analysis.

### Import
data <- read_csv(
  url,
  na = c("", "NA", "N/A", "NULL")
)

## Rows: 1296 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Department, Fever, DegreeDehy, OutComeDisc
## dbl (5): ID, Sex, Age, DiaMaxNo_24hrs, DiaDurInDays
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

### Clean variable types
data <- data %>%
  mutate(
    Department = factor(str_to_upper(str_trim(Department))),
    Sex = factor(Sex, levels = c(1, 2),
                 labels = c("Male", "Female")),
    Fever = factor(str_to_title(str_trim(Fever)),
                   levels = c("No", "Yes")),
    DegreeDehy = factor(
      str_to_lower(str_trim(DegreeDehy)),
      levels = c("some", "severe"),
      labels = c("Some", "Severe")
    ),
    OutComeDisc = factor(str_trim(OutComeDisc))
  )

### Audit
dim(data)

## [1] 1296    9

str(data)

## tibble [1,296 × 9] (S3: tbl_df/tbl/data.frame)
##  $ ID            : num [1:1296] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Department    : Factor w/ 2 levels "IPD","OPD": 2 2 1 1 1 2 1 1 1 1 ...
##  $ Sex           : Factor w/ 2 levels "Male","Female": 1 2 1 1 1 2 2 1 2 1 ...
##  $ Age           : num [1:1296] 6 8 12 2 12 12 10 21 6 37 ...
##  $ Fever         : Factor w/ 2 levels "No","Yes": 1 1 1 2 2 1 1 1 1 1 ...
##  $ DiaMaxNo_24hrs: num [1:1296] 7 5 10 9 8 7 12 8 8 10 ...
##  $ DiaDurInDays  : num [1:1296] 2 1 3 1 2 5 2 2 2 2 ...
##  $ DegreeDehy    : Factor w/ 2 levels "Some","Severe": 1 1 1 1 1 1 1 1 1 1 ...
##  $ OutComeDisc   : Factor w/ 1 level "Improved": NA NA 1 1 1 NA 1 1 1 1 ...

colSums(is.na(data))

##             ID     Department            Sex            Age          Fever 
##              0              0              0              0              1 
## DiaMaxNo_24hrs   DiaDurInDays     DegreeDehy    OutComeDisc 
##              0              2              0            196

### Describe numerical variables
data %>%
  summarise(
    across(
      c(Age, DiaMaxNo_24hrs, DiaDurInDays),
      list(
        mean = ~ mean(.x, na.rm = TRUE),
        median = ~ median(.x, na.rm = TRUE),
        sd = ~ sd(.x, na.rm = TRUE),
        IQR = ~ IQR(.x, na.rm = TRUE)
      )
    )
  )

### Describe categorical variables
data %>%
  count(Department) %>%
  mutate(Percent = 100 * n / sum(n))

### Visualize
ggplot(data, aes(x = Age)) +
  geom_histogram(bins = 20, color = "white") +
  labs(
    title = "Distribution of Age",
    x = "Age",
    y = "Frequency"
  )

### Compare groups
ggplot(data, aes(x = Department, y = Age)) +
  geom_boxplot() +
  labs(
    title = "Age by Department",
    x = "Department",
    y = "Age"
  )

### Examine a relationship
ggplot(data, aes(x = Age, y = DiaDurInDays)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

29 Exercises

29.1 Exercise 1: Data Audit

Using dirdata:

1. How many observations are present?
2. How many variables are present?
3. Which variables are categorical?
4. Which variables are numerical?
5. Which variables contain missing values?
6. Are there duplicated IDs?
7. What unique values occur in DegreeDehy before cleaning?
8. What unique values occur after cleaning?

29.2 Exercise 2: Descriptive Statistics

Calculate for Age, DiaMaxNo_24hrs, and DiaDurInDays:

1. mean;
2. median;
3. standard deviation;
4. minimum;
5. maximum;
6. quartiles;
7. IQR.

Then explain why the mean and median may differ.

29.3 Exercise 3: Categorical Variables

Calculate frequency and percentage tables for:

1. Department;
2. Sex;
3. Fever;
4. Degree of dehydration;
5. Outcome.

Then create a bar chart for each variable.

29.4 Exercise 4: Group Comparison

Compare Age between Departments using:

1. a boxplot;
2. a boxplot with jittered observations;
3. group-specific descriptive statistics.

Write a short paragraph describing what the graphs show.

29.5 Exercise 5: Diarrhea Duration

Study DiaDurInDays:

1. create a histogram;
2. create a density plot;
3. create a boxplot;
4. identify potential outliers;
5. investigate unusually large values;
6. calculate mean, median, SD and IQR.

Do not remove observations simply because they are outliers.

29.6 Exercise 6: Relationships

Investigate:

1. Age versus DiaDurInDays;
2. DiaMaxNo_24hrs versus DiaDurInDays.

For each:

• make a scatterplot;
• add a fitted trend line;
• calculate Pearson correlation;
• calculate Spearman correlation;
• write a short interpretation.

29.7 Exercise 7: Two Categorical Variables

Investigate the relationship between:

• Department and Sex;
• Department and Fever;
• Department and DegreeDehy.

For each:

1. create a contingency table;
2. calculate appropriate percentages;
3. create a grouped or proportional bar chart;
4. describe the observed pattern.

29.8 Exercise 8: Research-Style Figure

Create one final figure suitable for inclusion in a report or manuscript. It must contain:

• a meaningful title/caption;
• clear axis labels;
• appropriate units where known;
• readable text;
• no unnecessary decoration.

29.9 Exercise 9: Short Research Summary

Write 150–250 words describing the dataset using:

• one numerical summary table;
• one categorical summary;
• one distribution graph;
• one group-comparison graph;
• one relationship graph.

Do not make causal claims.

30 Suggested Final Student Assignment

Conduct a complete exploratory analysis of the diarrhea dataset.

Your report should contain:

30.0.1 Data description

• number of observations;
• number of variables;
• variable types;
• missing-data summary.

30.0.2 Descriptive statistics

Report appropriate statistics for numerical variables.

30.0.3 Categorical summaries

Report frequencies and percentages.

30.0.4 Visualizations

Include at least:

1. one histogram;
2. one boxplot;
3. one categorical bar chart;
4. one group-comparison graph;
5. one scatterplot.

30.0.5 Data-quality discussion

Discuss:

• missing values;
• inconsistent labels;
• potential outliers;
• any values requiring verification.

30.0.6 Interpretation

Explain what the analysis shows without making unsupported causal claims.

31 Common Mistakes to Avoid

31.1 Calling every numerical variable continuous

Counts such as DiaMaxNo_24hrs are discrete counts.

31.2 Assuming units that are not documented

Do not label Age as months, years, or another unit unless the original documentation confirms it.

31.3 Treating missing values as zero

A missing value is not automatically zero.

31.4 Deleting outliers automatically

First investigate whether the value is a valid observation.

31.5 Confusing correlation with causation

A correlation does not demonstrate a causal relationship.

31.6 Using the wrong graph

• categorical → usually bar chart;
• numerical distribution → histogram/density;
• numerical by group → boxplot/violin;
• two numerical variables → scatterplot.

31.7 Overusing pie charts

Bar charts generally make category comparisons easier.

31.8 Using too many colors

Use colour when it communicates information, not simply for decoration.

31.9 Reporting only p-values

Statistical analysis should also communicate effect size, uncertainty, direction, and practical meaning where appropriate.

32 Quick Reference

32.1 Descriptive statistics

Purpose	R function
Mean	mean()
Median	median()
Minimum	min()
Maximum	max()
Range	range()
Variance	var()
Standard deviation	sd()
Quartiles	quantile()
IQR	IQR()
Five-number summary	fivenum()

32.2 Frequency analysis

Purpose	R function
Frequency	table() / count()
Proportion	prop.table()
Cross-tabulation	table()
Tidy cross-tabulation	count()

32.3 ggplot2

Purpose	Function
Start graph	ggplot()
Histogram	geom_histogram()
Density	geom_density()
Boxplot	geom_boxplot()
Violin plot	geom_violin()
Points	geom_point()
Jitter	geom_jitter()
Bar chart	geom_bar()
Trend line	geom_smooth()
Labels	labs()
Faceting	facet_wrap() / facet_grid()
Flip coordinates	coord_flip()
Save figure	ggsave()

33 Graph Selection Guide

START
  |
  +-- What type of variable?
        |
        +-- Categorical
        |      |
        |      +-- One variable
        |      |      -> Bar chart
        |      |
        |      +-- Two variables
        |             -> Grouped / proportional bar chart
        |
        +-- Numerical
               |
               +-- One variable
               |      -> Histogram / Density / Boxplot
               |
               +-- Numerical + categorical
               |      -> Boxplot / Violin / Faceted distribution
               |
               +-- Two numerical variables
                      -> Scatterplot
                      -> Add trend line if appropriate
                      -> Consider correlation

34 Final Checklist

Before presenting an analysis, ask:

• Did I inspect the raw data first?
• Did I verify the number of observations and variables?
• Did I check missing values?
• Did I check inconsistent category labels?
• Did I check duplicate IDs?
• Did I verify variable types?
• Did I avoid assuming undocumented units?
• Did I choose an appropriate statistic?
• Did I choose an appropriate graph?
• Are axes and labels clear?
• Are units reported when known?
• Did I investigate unusual observations?
• Did I avoid automatically deleting outliers?
• Did I distinguish association from causation?
• Can another researcher reproduce my analysis from the code?

35 Final Remarks

Good statistical analysis is more than producing numbers and attractive graphs.

A strong exploratory analysis follows a logical sequence:

1. understand the data;
2. audit data quality;
3. clean and document variables;
4. summarize the data numerically;
5. visualize distributions;
6. compare relevant groups;
7. examine relationships;
8. investigate unusual observations;
9. interpret results cautiously;
10. communicate the evidence clearly.

The diarrhea dataset provides a useful teaching example because it contains numerical and categorical variables, missing values, inconsistent text labels, group structure, and observations that warrant investigation. These features make it suitable not only for learning ggplot2 and descriptive statistics, but also for learning an essential research skill: checking the data before trusting the analysis.