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 errors
library(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, aka causes of death. That way we can better review the options.

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 rows
EventDiagnosisUnique <- 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 rows

By 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 diagnoses, and assume that all individuals listed in the observation were afflicted by the diagnosis. This is because it is unrealistic labor-wise 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.

Metadata Source: U.S. Geological Survey National Wildlife Health Center, [2026], 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. We can anticipate that Cholera and Trauma will be our leading causes of death, but let’s confirm with further analysis.

Statistical Analysis

H0 : μ1=μ2=μ3=…μ162  All causes of death will show an equal average of mortalities.

Ha: μ1μ2≠μ3≠…μ162  The average of mortalities between causes of death will not be equal.

Now, we are going to calculate the mean mortality of 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. We can create a visualization to show our leading causes of death. Let’s see if we can get the top 5 leading causes of death on average, and we can do this by adjusting the subset of Average_Mortality.

newdata <- subset(WHISPersinfo, Average_Mortality >= 1000,
                  select=c(Event_Diagnosis, Average_Mortality))
                  
newdata %>% ggplot(aes(Event_Diagnosis, Average_Mortality))+
  geom_col(aes(fill=Event_Diagnosis))+
  labs(title = "Leading Causes of Death on Average in Maryland, DC, Virginia Bird Mortality Events",
       subtitle = "Avian Cholera claims most lives across 50 years")+
  scale_fill_manual(values = c("coral", "lightgreen", "yellow", "orange", "red"))+
  coord_flip()+
  theme(legend.position="none")

It appears that the leading causes of death on average are 1. Avian Cholera 2. Bacterial Infection (Erysipelothrix rhusiopathae) 3. Trauma 4. Toxicosis (fenthion) and 5. Steatitis.

For our hypothesis test, we are going to use an ANOVA which compares our means across groups. There are over 150 unique groups in our Event_Diagnosis column, so other tests to compare means would not be effective. We are looking for a low p-value (<0.05) to determine significance.

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

These results with a p-value of 2e-16 is extremely close to zero. That shows we can interpret our means as significant in respect to one another, and that the diagnosis has a meaningful effect on the value of the mean. We can then assume that our mean leading causes of death is interpretable data.

Conclusion

Previously, we saw that one event with 35,000 deaths had both Trauma and Cholera listed as the diagnosis. But, because we’re using the mean, we adjusted for this and showed that based on number of events, Trauma is not as impactful. In general, a single large mortality event is more likely to be caused by disease, and while monitoring traumas and preventing them is certainly necessary, when considering time and resources it would be more effective for Wildlife monitoring to target diseases such as Cholera, Bacterial Infection, Toxicosis, and Steatitis. Avian Cholera is responsible for more deaths on average than Trauma and the other 3 causes of death combined.

Some limitations with this data are that it’s reported by multiple different agencies, and the data is taken opportunistically. There is no guarantee that this is a compilation of all significant death events in the last 50 years, since wildlife research is often disjointed and underfunded. The data does represent a best-effort approach and solid collaboration from multiple agencies and organizations. Additionally, in some events small samples of total deaths were actually analyzed and then the results were projected onto all of the birds (this is due to making the most of time, money, and effort). This analysis would be best used to indicate where future research should take place–such as investigating Avian Cholera, Bacterial Infections, Toxicosis, in the Maryland, DC, Virginia area.