Module 8: Exploratory Data Analysis

Introduction

Before running a statistical model, calculating an index, or interpreting ecological patterns, we should first understand our data.

This process is called Exploratory Data Analysis, or EDA.

EDA helps us answer questions such as:

  • How much data do we have?
  • What does one row represent?
  • Which variables are available?
  • Are there missing or unusual values?
  • Which species were recorded?
  • How are records distributed among camera stations?
  • How are records distributed through time?
  • Are there patterns that deserve further investigation?

EDA is about understanding the dataset well enough to decide what questions we can reasonably ask next.

Importante

Exploration comes before interpretation.

Finding a pattern in a graph does not necessarily mean that we have found an ecological mechanism.

Let’s get started

Open your intro-r-course project and create a new script.

Save it inside:

scripts/

as:

module-08.R

Load the tidyverse:

library(tidyverse)

Then import our camera-trap dataset:

recordTable <- read_csv(
  "data/recordTable.csv"
)
Download the example data

We will continue using the same camera-trap record table from the previous modules.

Download recordTable.csv

Step 1: Understand the structure

Before calculating anything, ask:

What does one row represent?

In recordTable, each row represents one camera-trap record in our dataset.

Let’s inspect the data:

dim(recordTable)
[1] 3898   25
names(recordTable)
 [1] "RootFolder"        "File"              "RelativePath"     
 [4] "Station"           "Camera"            "DateTime"         
 [7] "Species"           "N_individuals"     "Sex"              
[10] "Age"               "Identifier"        "Comments"         
[13] "Favorite"          "DateTime_original" "year"             
[16] "date_shift"        "n_images"          "Date"             
[19] "Time"              "delta.time.secs"   "delta.time.mins"  
[22] "delta.time.hours"  "delta.time.days"   "DateTime2"        
[25] "Date2"            
glimpse(recordTable)
Rows: 3,898
Columns: 25
$ RootFolder        <chr> "Field Station 2025", "Field Station 2025", "Field S…
$ File              <chr> "DSCF0007.JPG", "DSCF0134.JPG", "DSCF0425.JPG", "DSC…
$ RelativePath      <chr> "Control 0\\Control 0 A", "Control 1\\Control 1 B", …
$ Station           <chr> "Control 0", "Control 1", "Control 1", "Control 1", …
$ Camera            <chr> "Control 0 A", "Control 1 B", "Control 1 B", "Contro…
$ DateTime          <dttm> 2025-05-15 14:53:20, 2025-06-05 04:25:19, 2025-06-2…
$ Species           <chr> "White-tailed deer", "Armadillo", "Armadillo", "Arma…
$ N_individuals     <dbl> 0, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ Sex               <chr> NA, "Unknown", "Unknown", "Unknown", "Unknown", "Unk…
$ Age               <chr> "Unknown", "Adult", "Unknown", "Unknown", "Unknown",…
$ Identifier        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ Comments          <chr> NA, NA, NA, NA, NA, NA, NA, NA, "ukn Lepidoptera", "…
$ Favorite          <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL…
$ DateTime_original <dttm> 2025-05-15 14:53:20, 2025-03-08 04:25:19, 2025-03-2…
$ year              <dbl> 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025…
$ date_shift        <dbl> NA, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, 89, …
$ n_images          <dbl> 1, 2, 2, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 1, 2…
$ Date              <date> 2025-05-15, 2025-06-05, 2025-06-26, 2025-06-26, 202…
$ Time              <time> 14:53:20, 04:25:19, 07:39:44, 08:07:20, 08:39:46, 0…
$ delta.time.secs   <dbl> 0, 0, 1826065, 1656, 1946, 2046685, 3601, 0, 2332886…
$ delta.time.mins   <dbl> 0.0, 0.0, 30434.4, 27.6, 32.4, 34111.4, 60.0, 0.0, 3…
$ delta.time.hours  <dbl> 0.0, 0.0, 507.2, 0.5, 0.5, 568.5, 1.0, 0.0, 648.0, 0…
$ delta.time.days   <dbl> 0.0, 0.0, 21.1, 0.0, 0.0, 23.7, 0.0, 0.0, 27.0, 0.0,…
$ DateTime2         <dttm> 2025-05-15 14:53:20, 2025-06-05 04:25:19, 2025-06-2…
$ Date2             <date> 2025-05-15, 2025-06-05, 2025-06-26, 2025-06-26, 202…

We can also look at the first observations:

head(recordTable)
# A tibble: 6 × 25
  RootFolder       File  RelativePath Station Camera DateTime            Species
  <chr>            <chr> <chr>        <chr>   <chr>  <dttm>              <chr>  
1 Field Station 2… DSCF… "Control 0\… Contro… Contr… 2025-05-15 14:53:20 White-…
2 Field Station 2… DSCF… "Control 1\… Contro… Contr… 2025-06-05 04:25:19 Armadi…
3 Field Station 2… DSCF… "Control 1\… Contro… Contr… 2025-06-26 07:39:44 Armadi…
4 Field Station 2… DSCF… "Control 1\… Contro… Contr… 2025-06-26 08:07:20 Armadi…
5 Field Station 2… DSCF… "Control 1\… Contro… Contr… 2025-06-26 08:39:46 Armadi…
6 Field Station 2… DSCF… "Control 1\… Contro… Contr… 2025-07-20 01:11:11 Armadi…
# ℹ 18 more variables: N_individuals <dbl>, Sex <chr>, Age <chr>,
#   Identifier <chr>, Comments <chr>, Favorite <lgl>, DateTime_original <dttm>,
#   year <dbl>, date_shift <dbl>, n_images <dbl>, Date <date>, Time <time>,
#   delta.time.secs <dbl>, delta.time.mins <dbl>, delta.time.hours <dbl>,
#   delta.time.days <dbl>, DateTime2 <dttm>, Date2 <date>

These simple functions already tell us a lot about the dataset.

For example:

dim()       → How large is the dataset?

names()     → What variables are available?

glimpse()   → What type of information does each variable contain?

head()      → What do the actual observations look like?
Tip

Before analyzing any new dataset, try to answer:

What does one row represent?

What does one column represent?

Many analysis problems begin when we misunderstand the unit represented by the data.

Step 2: How much information do we have?

We can create a very simple overview of the dataset.

recordTable |>
  summarise(
    records = n(),
    stations = n_distinct(Station),
    species = n_distinct(Species)
  )
# A tibble: 1 × 3
  records stations species
    <int>    <int>   <int>
1    3898       17      22

This gives us three useful numbers:

number of records
number of camera stations
number of species/categories

We can also look at the temporal coverage:

recordTable |>
  summarise(
    first_record = min(Date),
    last_record = max(Date)
  )
# A tibble: 1 × 2
  first_record last_record
  <date>       <date>     
1 2025-02-12   2025-10-02

This gives us a quick description of the data available for exploration.

Step 3: Look for missing information

Real ecological datasets almost always contain missing values.

Missing information is not necessarily a problem.

But we need to know where it occurs.

For example:

recordTable |>
  summarise(
    missing_sex = sum(is.na(Sex)),
    missing_age = sum(is.na(Age)),
    missing_identifier = sum(is.na(Identifier))
  )
# A tibble: 1 × 3
  missing_sex missing_age missing_identifier
        <int>       <int>              <int>
1         100          87               3692

Remember:

is.na()

returns TRUE when a value is missing.

Because R treats TRUE as 1 and FALSE as 0 in this context:

sum(is.na(Age))

counts the number of missing values.

Missing values are not the whole story

Sometimes missing or uncertain information is represented using a category instead of NA.

For example:

recordTable |>
  count(
    Age,
    sort = TRUE
  )
# A tibble: 4 × 2
  Age          n
  <chr>    <int>
1 Adult     2563
2 Unknown   1177
3 <NA>        87
4 Juvenile    71

and:

recordTable |>
  count(
    Sex,
    sort = TRUE
  )
# A tibble: 4 × 2
  Sex         n
  <chr>   <int>
1 Unknown  2727
2 Famale    547
3 Male      524
4 <NA>      100

A value such as:

Unknown

is not technically missing to R.

It is a character value.

But biologically, it may still represent information that we do not know.

Importante

R only knows what is written in the dataset.

It does not know that:

Unknown

may represent missing biological information.

Understanding what values mean requires knowledge of how the data were collected.

Step 4: Look for unusual values

EDA can also help us identify values that deserve closer inspection.

For example:

summary(recordTable$N_individuals)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.000   1.000   1.074   1.000   7.000

and:

summary(recordTable$n_images)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   1.000   2.000   2.491   2.000  61.000

We can visualize them too:

ggplot(
  recordTable,
  aes(x = N_individuals)
) +
  geom_histogram(
    binwidth = 1,
    fill = "steelblue",
    color = "white"
  ) +
  labs(
    x = "Number of individuals",
    y = "Number of records"
  ) +
  theme_bw()

Now ask:

  • Are there zeros?
  • Are there unusually large values?
  • Do those values make biological sense?
  • Could they represent data-entry problems?
  • Could they represent legitimate observations?

For example, we can inspect records with zero individuals:

recordTable |>
  filter(
    N_individuals == 0
  ) |>
  select(
    Station,
    Species,
    N_individuals,
    Date,
    File
  )
# A tibble: 71 × 5
   Station   Species           N_individuals Date       File            
   <chr>     <chr>                     <dbl> <date>     <chr>           
 1 Control 0 White-tailed deer             0 2025-05-15 DSCF0007.JPG    
 2 Control 2 Bug                           0 2025-09-01 DSCF0322 (7).JPG
 3 Control 2 Bug                           0 2025-09-01 DSCF0325 (3).AVI
 4 Control 2 Bug                           0 2025-09-11 DSCF0738 (7).JPG
 5 Control 2 Bug                           0 2025-09-11 DSCF0741 (3).AVI
 6 Control 2 Human                         0 2025-05-13 DSCF0451 (2).JPG
 7 Control 4 Human                         0 2025-06-03 DSCF0001 (8).JPG
 8 Control 4 Human                         0 2025-06-03 DSCF0045.JPG    
 9 Control 4 Human                         0 2025-06-08 DSCF0133.JPG    
10 Control 4 Wild turkey                   0 2025-06-18 DSCF0248.AVI    
# ℹ 61 more rows
Tip

An unusual value is not automatically an error.

EDA helps us identify observations that should be reviewed, not automatically deleted.

Step 5: Explore categorical variables

Categorical variables should also be inspected.

For example:

recordTable |>
  count(
    Species,
    sort = TRUE
  )
# A tibble: 22 × 2
   Species                   n
   <chr>                 <int>
 1 White-tailed deer      1266
 2 Armadillo               554
 3 Racoon                  522
 4 Hog                     454
 5 Virginia opossum        295
 6 Eastern gray squirrel   244
 7 Bug                      99
 8 Human                    80
 9 Bird                     78
10 Coyote                   75
# ℹ 12 more rows

This shows which species or categories occur most frequently in the record table.

We can visualize the result:

recordTable |>
  count(Species,sort = TRUE ) |>
  ggplot( aes(  x = reorder(Species, n),  y = n  )) +
  geom_col( fill = "steelblue") +
  coord_flip() +
  labs( x = "Species",
    y = "Number of records",
    title = "Camera-trap records by species"  ) +
  theme_bw()

Notice that we are combining:

data manipulation
        ↓
visualization

inside the same pipeline.

Be careful with record counts

A species with more records is not necessarily more abundant.

For example:

more records

could result from:

  • more activity,
  • longer visits in front of the camera,
  • differences in detection probability,
  • differences in camera effort,
  • or actual differences in abundance.
Records are not the same as abundance

Throughout this exploratory analysis, we are describing the records contained in the dataset.

We should not automatically interpret:

more records = more animals

Camera-trap detections are influenced by both ecological and observation processes.

Step 6: Explore camera stations

We can also ask:

How many records came from each camera station?

station_records <- recordTable |>
  count(
    Station,
    sort = TRUE
  )

station_records
# A tibble: 17 × 2
   Station               n
   <chr>             <int>
 1 Control 2          1029
 2 Fish 0              759
 3 Perfume 0           430
 4 Fish 2              402
 5 Perfume 2           252
 6 Control 4           232
 7 Control 3           214
 8 Hardwood Canopy 1   158
 9 Fish 1              156
10 Control 1            93
11 Hardwood Canopy 2    64
12 Fish 3               42
13 Perfume 1            42
14 Pine Canopy 1        11
15 Pine Canopy 3         9
16 Pine Canopy 2         4
17 Control 0             1

And visualize them:

ggplot(
  station_records,
  aes(
    x = reorder(Station, n),
    y = n
  )
) +
  geom_col(
    fill = "steelblue"
  ) +
  coord_flip() +
  labs(
    x = "Camera station",
    y = "Number of records",
    title = "Records by camera station"
  ) +
  theme_bw()

Some stations may have many more records than others.

Why?

Possible explanations include:

animal activity
habitat
camera placement
camera operation time
technical failures
sampling duration

The graph alone cannot tell us which explanation is correct.

Records are not sampling effort

This distinction is especially important for camera-trap data.

Suppose:

Camera A = 500 records

Camera B = 50 records

That does not tell us how long each camera was operating.

Camera A may have operated for 100 days and Camera B for only 10 days.

Our recordTable contains detection records.

It does not by itself describe the complete sampling effort of cameras with no detections or periods when cameras were inactive.

Importante

Do not use the number of records at a station as a measure of camera effort.

To calculate effort, we need information about when each camera was operating.

We will return to this idea when we begin working with camera-trap-specific tools.

Step 7: Species richness by station

We can ask another descriptive question:

How many different species were recorded at each station?

station_summary <- recordTable |>
  group_by(Station) |>
  summarise(
    records = n(),
    species_richness = n_distinct(Species),
    .groups = "drop"
  )

station_summary
# A tibble: 17 × 3
   Station           records species_richness
   <chr>               <int>            <int>
 1 Control 0               1                1
 2 Control 1              93                9
 3 Control 2            1029               13
 4 Control 3             214               14
 5 Control 4             232               14
 6 Fish 0                759               16
 7 Fish 1                156               11
 8 Fish 2                402               14
 9 Fish 3                 42                8
10 Hardwood Canopy 1     158                6
11 Hardwood Canopy 2      64                7
12 Perfume 0             430               11
13 Perfume 1              42                6
14 Perfume 2             252               13
15 Pine Canopy 1          11                5
16 Pine Canopy 2           4                2
17 Pine Canopy 3           9                4

Now each row represents:

one camera station

We can visualize species richness:

ggplot(
  station_summary,
  aes(
    x = reorder(Station, species_richness),
    y = species_richness
  )
) +
  geom_col(
    fill = "steelblue"
  ) +
  coord_flip() +
  labs(
    x = "Camera station",
    y = "Number of species",
    title = "Observed species richness by station"
  ) +
  theme_bw()

Notice the wording:

observed species richness

We are describing the number of species recorded in our dataset.

We are not claiming that this is the true number of species using each site.

Records and observed richness

We can also explore whether stations with more records tend to contain more observed species.

ggplot(
  station_summary,
  aes(
    x = records,
    y = species_richness
  )
) +
  geom_point(
    size = 3
  ) +
  labs(
    x = "Number of records",
    y = "Observed species richness"
  ) +
  theme_bw()

You may notice that stations with more records often also contain more species.

Does that mean that those stations have better habitat?

Not necessarily.

Nota

Stations with more sampling effort or more detections have more opportunities to record additional species.

This is why exploratory relationships should not automatically be interpreted as ecological effects.

Step 8: Explore species among stations

Sometimes we want to know where particular species were recorded.

Let’s focus on a few common species:

common_species <- c(
  "White-tailed deer",
  "Armadillo",
  "Racoon",
  "Hog",
  "Virginia opossum"
)

Now summarize the number of records for each station and species:

station_species <- recordTable |>
  filter(
    Species %in% common_species
  ) |>
  count(
    Station,
    Species
  )

station_species
# A tibble: 62 × 3
   Station   Species               n
   <chr>     <chr>             <int>
 1 Control 0 White-tailed deer     1
 2 Control 1 Armadillo             6
 3 Control 1 Hog                   7
 4 Control 1 Racoon               12
 5 Control 1 Virginia opossum      2
 6 Control 1 White-tailed deer    29
 7 Control 2 Armadillo           182
 8 Control 2 Hog                 210
 9 Control 2 Racoon               99
10 Control 2 Virginia opossum     73
# ℹ 52 more rows

We can visualize this using a heatmap:

ggplot(
  station_species,
  aes(
    x = Species,
    y = Station,
    fill = n
  )
) +
  geom_tile() +
  labs(
    x = "Species",
    y = "Camera station",
    fill = "Records",
    title = "Camera-trap records by station and species"
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

This type of graph can help us quickly identify:

  • species recorded at many stations,
  • species concentrated at particular stations,
  • stations with few records,
  • possible spatial patterns worth investigating.

Again, these patterns are exploratory.

Step 9: Explore patterns through time

Camera-trap datasets also contain a temporal component.

We can summarize the number of records per day:

daily_records <- recordTable |>
  count(Date)

daily_records
# A tibble: 180 × 2
   Date           n
   <date>     <int>
 1 2025-02-12     2
 2 2025-03-05     3
 3 2025-03-09     2
 4 2025-03-10     1
 5 2025-03-15     2
 6 2025-04-03     1
 7 2025-04-06     2
 8 2025-04-07     3
 9 2025-04-08    10
10 2025-04-09     1
# ℹ 170 more rows

Then visualize the pattern:

ggplot(
  daily_records,
  aes(
    x = Date,
    y = n
  )
) +
  geom_line() +
  labs(
    x = "Date",
    y = "Number of records",
    title = "Camera-trap records through time"
  ) +
  theme_bw()

Large changes through time could reflect:

  • animal activity,
  • seasonal changes,
  • differences in sampling effort,
  • cameras being deployed or removed,
  • camera failures,
  • or other changes in the study.

The plot helps us identify the pattern.

It does not tell us its cause.

Explore one species through time

We can also focus on one species.

For example:

deer_daily <- recordTable |>
  filter(
    Species == "White-tailed deer"
  ) |>
  count(Date)

Then:

ggplot(
  deer_daily,
  aes(
    x = Date,
    y = n
  )
) +
  geom_line() +
  labs(
    x = "Date",
    y = "White-tailed deer records",
    title = "White-tailed deer records through time"
  ) +
  theme_bw()

This could help us generate new questions.

For example:

Why are records higher during some portions of the sampling period?

That question might later require information about:

  • sampling effort,
  • weather,
  • season,
  • human activity,
  • or other ecological covariates.

That is exactly what EDA is useful for:

finding questions worth investigating.

Step 10: Explore demographic information

Our record table also contains variables such as:

Sex
Age

For example, we can explore age classifications for White-tailed deer:

deer_demographics <- recordTable |>
  filter(
    Species == "White-tailed deer"
  ) |>
  count(
    Age,
    sort = TRUE
  )

deer_demographics
# A tibble: 4 × 2
  Age          n
  <chr>    <int>
1 Adult      986
2 Unknown    223
3 Juvenile    40
4 <NA>        17

And visualize them:

ggplot(
  deer_demographics,
  aes(
    x = Age,
    y = n
  )
) +
  geom_col(
    fill = "steelblue"
  ) +
  labs(
    x = "Age class",
    y = "Number of records",
    title = "White-tailed deer records by age class"
  ) +
  theme_bw()

But again, be careful.

A graph of:

records by age class

is not necessarily the same as:

population age structure

Differences in behavior and detectability can influence camera-trap records.

Tables are part of EDA too

EDA does not always require a graph.

Sometimes a summary table communicates the pattern more clearly.

For example:

species_summary <- recordTable |>
  group_by(Species) |>
  summarise(
    records = n(),
    stations = n_distinct(Station),
    total_individuals = sum(
      N_individuals,
      na.rm = TRUE
    ),
    .groups = "drop"
  ) |>
  arrange(
    desc(records)
  )

species_summary
# A tibble: 22 × 4
   Species               records stations total_individuals
   <chr>                   <int>    <int>             <dbl>
 1 White-tailed deer        1266       16              1366
 2 Armadillo                 554       15               576
 3 Racoon                    522       11               538
 4 Hog                       454       11               574
 5 Virginia opossum          295        9               299
 6 Eastern gray squirrel     244        9               252
 7 Bug                        99        8               113
 8 Human                      80       15                63
 9 Bird                       78        9                83
10 Coyote                     75       11                79
# ℹ 12 more rows

Now we can quickly compare several characteristics among species.

A useful exploratory workflow often combines:

tables
+
graphs
+
knowledge of the study

EDA is an iterative process

Exploration rarely follows a perfectly linear sequence.

You may begin with:

How many species do I have?

Then notice:

One station has very few records.

That may lead you to ask:

Was that camera operating throughout the study?

Then:

Do all species show the same pattern?

Then:

Does the pattern change through time?

EDA often looks like:

question
   ↓
summary
   ↓
visualization
   ↓
something interesting or strange
   ↓
new question
   ↓
check the data
   ↓
repeat

What EDA cannot tell us

Exploratory analysis is extremely useful.

But it has limits.

From these graphs alone, we cannot conclude that:

a habitat causes more detections

one station has greater abundance

one species prefers a treatment

a temporal change reflects population change

Those conclusions require appropriate study design, effort information, and statistical analyses.

Importante

EDA helps us describe patterns and generate hypotheses.

It does not automatically establish causation or account for detectability.

Practice at home

Tip

Open your intro-r-course project and create:

practice-08.R

inside the scripts folder.

Use:

recordTable.csv

to conduct a small exploratory analysis.

Part 1: Describe the dataset

Find:

  1. The number of records.
  2. The number of camera stations.
  3. The number of species/categories.
  4. The first and last dates in the dataset.

Part 2: Data quality

Explore:

Sex
Age
N_individuals

Look for:

  • missing values,
  • unknown categories,
  • unusual values,
  • possible inconsistencies.

Do not automatically modify anything.

Simply identify observations that may deserve further review.

Part 3: Choose one species

Choose any species in the dataset.

Determine:

  1. How many records it has.
  2. At how many stations it was recorded.
  3. Which station contains the most records.
  4. How its records are distributed through time.

Part 4: Create two figures

Create at least two different exploratory figures.

For example:

records among stations
records through time
age or sex categories
number of individuals
station × species

Add appropriate:

titles
axis labels
theme

Part 5: Interpret carefully

Write three observations about the patterns you found.

Then, for each observation, write at least one alternative explanation or limitation.

For example:

Observation:
Station X had more deer records than Station Y.

Possible explanation:
Deer may use the area around Station X more frequently.

Important limitation:
We have not yet accounted for differences in camera effort between stations.

There is no single correct solution because exploratory analysis depends on the questions you ask.

One possible workflow could be:

library(tidyverse)

recordTable <- read_csv(
  "data/recordTable.csv"
)


# ---------------------------
# Dataset overview
# ---------------------------

recordTable |>
  summarise(
    records = n(),
    stations = n_distinct(Station),
    species = n_distinct(Species),
    first_record = min(Date),
    last_record = max(Date)
  )


# ---------------------------
# Inspect data quality
# ---------------------------

recordTable |>
  count(
    Sex,
    sort = TRUE
  )

recordTable |>
  count(
    Age,
    sort = TRUE
  )

summary(
  recordTable$N_individuals
)


# ---------------------------
# Explore one species
# ---------------------------

deer_records <- recordTable |>
  filter(
    Species == "White-tailed deer"
  )

deer_records |>
  summarise(
    records = n(),
    stations = n_distinct(Station)
  )


# Records by station

deer_station <- deer_records |>
  count(
    Station,
    sort = TRUE
  )

ggplot(
  deer_station,
  aes(
    x = reorder(Station, n),
    y = n
  )
) +
  geom_col(
    fill = "steelblue"
  ) +
  coord_flip() +
  labs(
    x = "Camera station",
    y = "Number of records"
  ) +
  theme_bw()


# Records through time

deer_daily <- deer_records |>
  count(Date)

ggplot(
  deer_daily,
  aes(
    x = Date,
    y = n
  )
) +
  geom_line() +
  labs(
    x = "Date",
    y = "Number of records"
  ) +
  theme_bw()

The purpose is not to produce a specific answer.

The purpose is to practice asking:

What do I see?

Does it make sense?

What could explain it?

What information am I missing?

What should I investigate next?