Project2
Mia Ramirez
2026-07-28
Which causes of wild bird mortality events on average, are the deadliest in the Maryland, DC, Virginia area?
The WHISPers (Wildlife Health Information Sharing Partnership) project is a resource where USGS compiles ongoing information about wildlife death and illness events. The information collected is the result of either USGS researchers or partner programs submitting the data. The data isn’t collected systematically, meaning it’s not recorded on a regular basis. Instead, the data is often “see something, say something” style, where the reporting agencies may come across a wildlife event and record the data as it’s happening, rather than trying to seek out the events.
For our analysis, we are downloading a subset of the data. The data includes all events in Maryland, District of Columbia, and Virginia between the dates January 1st, 1976 and December 31st, 2025. The only wildlife selected are Birds.
The data can be viewed here: https://whispers.usgs.gov/home, and downloaded after applying the correct filters.
We will use the columns “EventID” “Number Affected” “Event Start Date” “Event End Date” “States (or equivalent)” “Counties (or equivalent)” “Species” and “Event Diagnosis” in our research.
library(readr)
library(readxl)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ purrr 1.2.2
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(dplyr)Cleaning the Data
First, we load the data into the project using the readr function read_csv.
Clean column titles by replacing spaces with _. Also adjust the column titles to reflect that Washington, D.C. is included in this dataset.
WHISPersdata <- read_csv("WHISPersData.csv")## Rows: 507 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): Event Type, Countries, States (or equivalent), Counties (or equiva...
## dbl (2): Event ID, Number Affected
## lgl (1): Public
## date (2): Event Start Date, Event End Date
##
## ℹ 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.WHISPersdata <- rename_with(WHISPersdata, ~ gsub(" ", "_", .))
WHISPersdata <- rename(WHISPersdata, State_or_District = "States_(or_equivalent)", County_or_District = "Counties_(or_equivalent)")We’re interested in causes of death, so let’s create an object that only contains a list of the different Event Diagnoses, or causes of death.
EventDiagnosisList <- as_tibble(WHISPersdata$Event_Diagnosis)
EventDiagnosisList <- unique(EventDiagnosisList)
EventDiagnosisList <- arrange(EventDiagnosisList, value)
print(EventDiagnosisList)## # A tibble: 176 × 1
## value
## <chr>
## 1 Aflatoxicosis
## 2 Aspergillosis
## 3 Aspergillosis suspect
## 4 Aspergillosis suspect; Avian Pox suspect; Bacterial Infection (multisystemic…
## 5 Aspergillosis; Bacterial Infection (gastrointestinal/hepatic); Nutritional D…
## 6 Aspergillosis; Botulism Type C; Bacterial Infection (Escherichia coli) suspe…
## 7 Aspergillosis; Circovirus; Nutritional Deficiency; Undetermined
## 8 Aspergillosis; Egg Yolk Peritonitis; Emaciation (starvation)
## 9 Aspergillosis; Emaciation (not otherwise specified)
## 10 Aspergillosis; Highly Pathogenic Avian Influenza (AI virus H5N1)
## # ℹ 166 more rowsEventDiagnosisUnique <- EventDiagnosisList %>% separate_rows(value, sep = ";")
EventDiagnosisUnique <- unique(EventDiagnosisUnique)
print(EventDiagnosisUnique)## # A tibble: 162 × 1
## value
## <chr>
## 1 "Aflatoxicosis"
## 2 "Aspergillosis"
## 3 "Aspergillosis suspect"
## 4 " Avian Pox suspect"
## 5 " Bacterial Infection (multisystemic)"
## 6 " Bacterial Infection (musculoskeletal)"
## 7 " Emaciation (starvation)"
## 8 " Nutritional Deficiency"
## 9 " West Nile Virus suspect"
## 10 " Bacterial Infection (gastrointestinal/hepatic)"
## # ℹ 152 more rowsBy looking at this list, and selecting for unique values, we see that some of these events actually have multiple causes listed. For example, Row 7 (Event_ID 204681 when referencing original dataset) has Aspergillosis (a fungal disease), Circovirus (a broad range of viruses that present differently depending on species), Nutritional Deficiency, and Undetermined. The USGS website does offer further insight into this by looking up the event https://whispers.usgs.gov/event/204681. There, we find that of 3 dead individuals assessed by the researchers, 1 were diagnosed with Aspergillosis, 2 with Circovirus, and 3 with undetermined. By skimming through EventDiagnosisUnique, there are dozens of instances where an event has multiple diagnoses listed.
For the purposes of this analysis, we are going to split these observations and isolate the diagnosis, and assume that all individuals listed in the observation were afflicted by the diagnosis. This is because it is unrealistic to cross-reference the event data with all the extended information available on the USGS website.
We will use the separate_rows function from tidyverse, and luckily the USGS has helpfully used ; to divide between different diagnoses in one observation. We will use that value in the function. There are also some instances where the diagnosis is listed as “[diagnosis]” or “[diagnosis] suspect”. An event only has the diagnosis listed as suspect if there is no alternative non-suspect diagnosis. This was detailed in the metadata.
Because of this, we are also going to remove the word “suspect” from these diagnosis and treat them as if they are also confirmed.
Source:U.S. Geological Survey National Wildlife Health Center, [year accessed], Wildlife Health Information Sharing Partnership - event reporting system (WHISPers): U.S. Geological Survey database, accessed [July 29, 2026], at https://doi.org/10.5066/P9B0A3IM
WHISPersdataclean <- WHISPersdata %>% separate_rows(Event_Diagnosis, sep = ";")
WHISPersdataclean$Event_Diagnosis <- WHISPersdataclean$Event_Diagnosis %>% gsub("suspect", "", .)Next, we need to select the columns of importance. We also need to adjust the Event_ID from a number to a character, since the data is categorical (the ID points to a specific event, not a numerical value).
WHISPersdataclean <- WHISPersdataclean %>% select(Event_ID,
Number_Affected,
Event_Start_Date,
Event_End_Date,
State_or_District,
County_or_District,
Species,
Event_Diagnosis)
WHISPersdataclean$Event_ID <- as.character(WHISPersdataclean$Event_ID)
summary(WHISPersdataclean)## Event_ID Number_Affected Event_Start_Date Event_End_Date
## Length :671 Min. : 1.0 Min. :1976-08-10 Min. :1976-08-17
## N.unique :507 1st Qu.: 5.0 1st Qu.:1992-09-18 1st Qu.:1992-09-21
## N.blank : 0 Median : 18.0 Median :2005-09-14 Median :2005-12-01
## Min.nchar: 5 Mean : 248.3 Mean :2006-06-04 Mean :2006-06-12
## Max.nchar: 6 3rd Qu.: 60.0 3rd Qu.:2022-05-25 3rd Qu.:2022-05-27
## Max. :35000.0 Max. :2025-12-30 Max. :2025-12-30
## NAs :2
## State_or_District County_or_District Species Event_Diagnosis
## Length :671 Length : 671 Length :671 Length :671
## N.unique : 20 N.unique : 135 N.unique :247 N.unique :162
## N.blank : 0 N.blank : 0 N.blank : 0 N.blank : 0
## Min.nchar: 8 Min.nchar: 14 Min.nchar: 6 Min.nchar: 4
## Max.nchar:133 Max.nchar:1154 Max.nchar:498 Max.nchar: 65
## NAs : 3
## By looking at a summary of the data, we see some interesting things. For example, in Number_Affected aka the number of individuals killed in a mortality event, the 3rd quartile shows 60 individuals but the max is 35,000. That is a huge outlier, and must have been a truly devastating event. Whatever diagnoses are tied to that event are likely to be the leading cause of death since that number is going to skew the mean heavily. We can use the filter function to find this exact occurrence.
outlier35000 <- WHISPersdataclean %>% filter(Number_Affected == "35000")
print(outlier35000)## # A tibble: 2 × 8
## Event_ID Number_Affected Event_Start_Date Event_End_Date State_or_District
## <chr> <dbl> <date> <date> <chr>
## 1 12944 35000 1994-02-21 1994-04-07 Maryland; Virginia; …
## 2 12944 35000 1994-02-21 1994-04-07 Maryland; Virginia; …
## # ℹ 3 more variables: County_or_District <chr>, Species <chr>,
## # Event_Diagnosis <chr>Based on this, we can find the Event_ID, and the WHISPers database actually has a way to look into the specific details using an Event_ID! https://whispers.usgs.gov/event/12944
According to this detail there is not any more specific data on which birds were killed with Cholera vs. Trauma. We will assume that both are the potential cause of death.
Statistical Analysis
H0 : μ=μ0 All causes of death will show an equal average of mortalities.
Ha: μ≠μ0 The average of mortalities between causes of death will not be equal.
Now, we are going to calculate the mean mortality by each diagnosis. We will take the total number of affected / total number of events.
WHISPersinfo <- WHISPersdataclean %>%
group_by(Event_Diagnosis) %>%
summarise(
Average_Mortality = mean(Number_Affected),
Total_Mortality = sum(Number_Affected),
StDev_Mortality = sd(Number_Affected),
Number_of_Events = n())
summary(WHISPersinfo)## Event_Diagnosis Average_Mortality Total_Mortality StDev_Mortality
## Length :162 Min. : 1.00 Min. : 1.00 Min. : 0.00
## N.unique :162 1st Qu.: 10.00 1st Qu.: 13.25 1st Qu.: 14.11
## N.blank : 0 Median : 32.44 Median : 60.50 Median : 56.57
## Min.nchar: 4 Mean : 339.75 Mean : 1028.62 Mean : 280.13
## Max.nchar: 65 3rd Qu.: 100.00 3rd Qu.: 323.25 3rd Qu.: 137.89
## Max. :33147.50 Max. :66295.00 Max. :9009.63
## NAs :89
## Number_of_Events
## Min. : 1.000
## 1st Qu.: 1.000
## Median : 1.000
## Mean : 4.142
## 3rd Qu.: 3.000
## Max. :56.000
## Across 162 diagnoses, our highest Average Mortality is 33,147 and the lowest is 1. Quite a large range.
newdata <- subset(WHISPersinfo, Total_Mortality >= 5000,
select=c(Event_Diagnosis, Total_Mortality))
newdata %>% ggplot(aes(Event_Diagnosis,Total_Mortality)) +
geom_col()+
coord_flip()+
theme_grey()
anovadata <- WHISPersdataclean %>% select(Number_Affected,Event_Diagnosis)
anova_result <- aov(Number_Affected~Event_Diagnosis, anovadata)
summary(anova_result)## Df Sum Sq Mean Sq F value Pr(>F)
## Event_Diagnosis 161 2.301e+09 14288950 6.051 <2e-16 ***
## Residuals 509 1.202e+09 2361410
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1