PA2

Quanitifying Injuries and Economic impact of Weather Events in the United States Using Box-Plots Based on NWS Storm Data

Synopsis

This study analyzes the National Weather Service (NWS) Storm Data set to examine the relationships between weather events and their annual impact on health and economic damage across the United States. Annual fatalities, injury counts, and property damage are investigated via box-plots. A non-parametrics Kruskal-Walis test is applied to examine the relative impact of the worst weather event types. Our results show that Tornadoes pose a statistically significant higher risk to human health compared to the second leading cause of health-related impacts, namely heat (p-value < 0.05). Additionally, floods and thunderstorm winds are identified as the primary contributors to economic damage.

Data Processing

Raw Data Analysis

Raw data from the National Weather Service is used for the analysis (access date: 2025-06-02) and is stored in the “ProjectData” directory. Documentation of the data set is also available (access date: 2025-06-02).

# Create data directory
if (!dir.exists("ProjectData")) {
  dir.create("ProjectData")
}

# Load the data
download.file("https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2", 
              "storm_data.csv.bz2")
data_set_raw <- read.csv(bzfile("ProjectData/storm_data.csv.bz2"))

The possible weather event types according to the documentation are manually retrieved.

doc_event_types <- c(
  "ASTRONOMICAL LOW TIDE", "AVALANCHE", "BLIZZARD", "COASTAL FLOOD",
  "COLD/WIND CHILL", "DEBRIS FLOW", "DENSE FOG", "DENSE SMOKE", "DROUGHT",
  "DUST DEVIL", "DUST STORM", "EXCESSIVE HEAT", "EXTREME COLD/WIND CHILL",
  "FLASH FLOOD", "FLOOD", "FROST/FREEZE", "FUNNEL CLOUD", "FREEZING FOG",
  "HAIL", "HEAT", "HEAVY RAIN", "HEAVY SNOW", "HIGH SURF", "HIGH WIND",
  "HURRICANE (TYPHOON)", "ICE STORM", "LAKE-EFFECT SNOW", "LAKESHORE FLOOD",
  "LIGHTNING", "MARINE HAIL", "MARINE HIGH WIND", "MARINE STRONG WIND",
  "MARINE THUNDERSTORM WIND", "RIP CURRENT", "SEICHE", "SLEET",
  "STORM SURGE/TIDE", "STRONG WIND", "THUNDERSTORM WIND", "TORNADO",
  "TROPICAL DEPRESSION", "TROPICAL STORM", "TSUNAMI", "VOLCANIC ASH",
  "WATERSPOUT", "WILDFIRE", "WINTER STORM", "WINTER WEATHER"
)

doc_event_types_ordered <- doc_event_types[order(nchar(doc_event_types), decreasing = TRUE)]

A check is performed on the data set to confirm conformity to the documentation. According to the documentation, 48 types of weather events exists.

unique_event_types_raw <- unique(data_set_raw$EVTYPE)

However, 985 unique weather event types are detected in the data set. This is because the documentation allows event types to deviate from the standards. For example, measurements can be appended to the event type. Weather events may also consists of multiple events.

Data cleaning

Before we process the data events, we first look at the important variables of this study. For this, we focus on the date of the weather event, fatalities and injuries that occurred and property damage.

Variable event_date is constructed and represents the date when a weather event ends. If this information is missing from the data set, we use the start date of the event. If the start date is missing or invalid, the data sample is discarded.

# Create the EVENT_DATE column
data_set_raw <- data_set_raw %>%
  mutate(
    # Clean the input
    END_DATE_CLEAN = ifelse(END_DATE == "", NA, END_DATE),
    BGN_DATE_CLEAN = ifelse(BGN_DATE == "", NA, BGN_DATE),
    EVENT_DATE_RAW = ifelse(!is.na(END_DATE_CLEAN), 
                            END_DATE_CLEAN, BGN_DATE_CLEAN),
    # Convert to date
    EVENT_DATE = as.Date(sub(" .*", "", EVENT_DATE_RAW), 
                         format = "%m/%d/%Y")
  )

# Only keep valid rows
data_set_clean_1 <- data_set_raw %>%
  filter(!is.na(EVENT_DATE))

We exclude samples where the event type is missing. For records where fatalities, injuries, or property damage are ivnalid or missing, we assume that the intended value was zero.

# Clean the variables or interest
data_set_clean_2 <- data_set_clean_1 %>%
  mutate(
    # Convert fatalities to numeric if possible
    # Set "" to NA, then set missing values to 0
    FATALITIES = as.numeric(ifelse(FATALITIES == "", NA, FATALITIES)),
    FATALITIES = ifelse(is.na(FATALITIES), 0, FATALITIES),
    
    # Convert injuries to numeric if possible
    # Set "" to NA, then set missing values to 0
    INJURIES = as.numeric(ifelse(INJURIES == "", NA, INJURIES)),
    INJURIES = ifelse(is.na(INJURIES), 0, INJURIES),
    
    # Convert property damage to numeric if possible
    # Set "" to NA, then set missing values to 0
    PROPDMG = as.numeric(ifelse(PROPDMG == "", NA, PROPDMG)),
    PROPDMG = ifelse(is.na(PROPDMG), 0, PROPDMG)
  ) %>%
  # Remove rows with missing EVTYPE (NA or "")
  filter(!is.na(EVTYPE) & EVTYPE != "")

Weather event grouping

Using the weather events defined in the documentation, it is possible to use pattern checking to reclassify the data. The data is prepared for analysis as follows:

1. Weather event types in the data set are reclassified according to weather event notations of the documentation.
2. For every data sample, the weather event type reclassification is performed using pattern comparison, i.e., a weather event is checked if it contains any of the weather events of the documentation.
3. Weather events that do not exists in the documentation will be reclassified as “Other”.
4. Data samples that contain multiple weather event types will have all event types reclassified.

# Pattern match or set to "Other" if no match is found.
matches_list <- lapply(data_set_clean_2$EVTYPE, function(x) {
  # Pattern check
  matches <- doc_event_types_ordered[
    str_detect(x, fixed(doc_event_types_ordered, ignore_case = TRUE))
  ]
  if (length(matches) > 0) {
    matches
  } else {
    "Other"
  }
})

# Store the results
data_set_clean_2$EVTYPE_MATCHES <- matches_list

The data samples are distributed as follows:

# Count how many times each doc_event_type occurs
long_df <- data_set_clean_2 %>%
  mutate(row_id = row_number()) %>%
  select(row_id, EVTYPE, EVTYPE_MATCHES) %>%
  unnest(EVTYPE_MATCHES)

evtype_group_counts <- long_df %>%
  group_by(EVTYPE_MATCHES) %>%
  summarise(count = n()) %>%
  arrange(desc(count))

print(evtype_group_counts)

## # A tibble: 47 × 2
##    EVTYPE_MATCHES     count
##    <chr>              <int>
##  1 HAIL              290401
##  2 Other             238109
##  3 THUNDERSTORM WIND 109446
##  4 FLOOD              82731
##  5 TORNADO            60700
##  6 FLASH FLOOD        55668
##  7 HIGH WIND          21953
##  8 HEAVY SNOW         15802
##  9 LIGHTNING          15776
## 10 HEAVY RAIN         11817
## # ℹ 37 more rows

Within the “Other” group, the distribution is as follows:

# Display the distribution
other_counts <- long_df %>%
  filter(EVTYPE_MATCHES == "Other") %>%
  count(EVTYPE) %>%
  arrange(desc(n))

head(other_counts)

## # A tibble: 6 × 2
##   EVTYPE                    n
##   <chr>                 <int>
## 1 TSTM WIND            219940
## 2 MARINE TSTM WIND       6175
## 3 URBAN/SML STREAM FLD   3392
## 4 WILD/FOREST FIRE       1457
## 5 EXTREME COLD            655
## 6 LANDSLIDE               600

TSTM is by far the largest occurrence and has many counts. For this reason, we make TSTM a separate category and keep the rest in “Other”.

# Copy the data
data_set_clean <- data_set_clean_2

# Add a new category
doc_event_types_mod <- c(doc_event_types, "TSTM WIND")
doc_event_types_ordered_mod <- doc_event_types_mod[order(nchar(doc_event_types_mod), decreasing = TRUE)]

# Pattern match or set to "Other" if no match is found.
matches_list_clean <- lapply(data_set_clean$EVTYPE, function(x) {
  # Pattern match
  matches <- doc_event_types_ordered_mod[
    str_detect(x, fixed(doc_event_types_ordered_mod, ignore_case = TRUE))
  ]
  if (length(matches) > 0) {
    matches
  } else {
    "Other"
  }
})

# Store the results
data_set_clean$EVTYPE_MATCHES <- matches_list_clean

# Group by event
long_df_clean <- data_set_clean %>%
  mutate(row_id = row_number()) %>%
  select(row_id, EVTYPE, EVTYPE_MATCHES) %>%
  unnest(EVTYPE_MATCHES)

# Count how many times each doc_event_type occurs
evtype_group_counts_clean <- long_df_clean %>%
  group_by(EVTYPE_MATCHES) %>%
  summarise(count = n()) %>%
  arrange(desc(count))

head(evtype_group_counts_clean)

## # A tibble: 6 × 2
##   EVTYPE_MATCHES     count
##   <chr>              <int>
## 1 HAIL              290401
## 2 TSTM WIND         227230
## 3 THUNDERSTORM WIND 109446
## 4 FLOOD              82731
## 5 TORNADO            60700
## 6 FLASH FLOOD        55668

Finally, we construct the yearly injury counts and economic damage.

# Create the sum of injuries
data_set_clean$SUM_INJURIES <- data_set_clean$FATALITIES + data_set_clean$INJURIES

# Group by year
annual_summary <- data_set_clean %>%
  mutate(YEAR = year(EVENT_DATE)) %>%
  group_by(YEAR) %>%
  summarise(
    TOTAL_SUM_INJURIES = sum(SUM_INJURIES, na.rm = TRUE),
    TOTAL_PROPDMG = sum(PROPDMG, na.rm = TRUE)
  )

Results

The figures below show box plots of the sum of the injuries and property damage for the top 10 weather events.

# Unpack samples with multiple events
df_long <- data_set_clean %>%
  unnest(EVTYPE_MATCHES)

# Determine the annual sums
annual_sums_by_type <- df_long %>%
  mutate(YEAR = year(EVENT_DATE)) %>%
  group_by(EVTYPE_MATCHES, YEAR) %>%
  summarise(
    ANNUAL_SUM_INJURIES = sum(SUM_INJURIES, na.rm = TRUE),
    .groups = "drop"
  )

# Determine the top n-events
top_events <- annual_sums_by_type %>%
  group_by(EVTYPE_MATCHES) %>%
  summarise(MEDIAN_INJURIES = median(ANNUAL_SUM_INJURIES, na.rm = TRUE)) %>%
  arrange(desc(MEDIAN_INJURIES)) %>%
  slice_head(n = 10) %>%
  pull(EVTYPE_MATCHES)

annual_sums_top_events <- annual_sums_by_type %>%
  filter(EVTYPE_MATCHES %in% top_events)

# Plot the top n-events
ggplot(
  annual_sums_top_events,
  aes(
    x = reorder(EVTYPE_MATCHES, ANNUAL_SUM_INJURIES, FUN = function(x) -median(x, na.rm = TRUE)),
    y = ANNUAL_SUM_INJURIES
  )
) +
  geom_boxplot(outlier.size = 1) +
  scale_y_log10() +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  ) +
  labs(
    x = "Event Type",
    y = "Annual Total Injuries",
    title = "Boxplot of Annual Injuries by Top 10 Weather Events"
  )

The boxplot above shows that Tornado, Heat, Exessive heat are the top 3 leading causes of health-related issues. A Kruskal-Wallis test is done below with alpha 0.05 and a Bonferroni correction of 2 to show statistical significance between Tornado and Heat:

top_two_events <- top_events[1:2]

# Filter for the top 2 events
top_two_data <- annual_sums_by_type %>%
  filter(EVTYPE_MATCHES %in% top_two_events)

# Apply the Kruskal-Wallis test
kw_result <- kruskal.test(
  ANNUAL_SUM_INJURIES ~ EVTYPE_MATCHES,
  data = top_two_data
)

# Apply Bonferroni correction
p_value_corrected_1 <- min(kw_result$p.value * 2, 1)

The Bonferroni-corrected p-value is 8.6069185^{-5} indicating statistical signficance.

The figure below show a box plot of the sum of property damage of the top 10 weather events.

# Unpack samples with multiple events
df_long <- data_set_clean %>%
  unnest(EVTYPE_MATCHES)

# Summarize annual property damage totals by event type
annual_sums_by_type <- df_long %>%
  mutate(YEAR = year(EVENT_DATE)) %>%
  group_by(EVTYPE_MATCHES, YEAR) %>%
  summarise(
    ANNUAL_SUM_PROPDMG = sum(PROPDMG, na.rm = TRUE),
    .groups = "drop"
  )

# Determine the top n-events
top_events <- annual_sums_by_type %>%
  group_by(EVTYPE_MATCHES) %>%
  summarise(MEDIAN_PROPDMG = median(ANNUAL_SUM_PROPDMG, na.rm = TRUE)) %>%
  arrange(desc(MEDIAN_PROPDMG)) %>%
  slice_head(n = 10) %>%
  pull(EVTYPE_MATCHES)

annual_sums_top_events <- annual_sums_by_type %>%
  filter(EVTYPE_MATCHES %in% top_events)

# Plot the top n-events
ggplot(
  annual_sums_top_events,
  aes(
    x = reorder(EVTYPE_MATCHES, ANNUAL_SUM_PROPDMG, FUN = function(x) -median(x, na.rm = TRUE)),
    y = ANNUAL_SUM_PROPDMG
  )
) +
  geom_boxplot(outlier.size = 1) +
  scale_y_log10() +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  ) +
  labs(
    x = "Weather Event Type",
    y = "Annual Total Property Damage",
    title = "Boxplot of Annual Property Damage by Top 10 Weather Events"
  )

The boxplot above shows that floods, thunderstorm wind and flash floods are the top 3 leading causes to economic damage. A Kruskal-Wallis test is done below with alpha 0.05 and a Bonferroni correction of 2 to show statistical significance between floods and thunderstorm winds:

top_two_events <- top_events[1:2]

# Filter annual_sums_by_type for the top 2 events
top_two_data <- annual_sums_by_type %>%
  filter(EVTYPE_MATCHES %in% top_two_events)

# Kruskal-Wallis test
kw_result <- kruskal.test(
  ANNUAL_SUM_PROPDMG ~ EVTYPE_MATCHES,
  data = top_two_data
)

# Apply Bonferroni correction
p_value_corrected_2 <- min(kw_result$p.value * 2, 1)

The Bonferroni-corrected p-value is 1, which does not indicate a statistical difference between floods and thunderstorm winds.

Discussion

In this document, a data processing method was described that uses pattern matching to reclassify event types in the data set to events defined in the Documentation. Note that some weather event types listed in the documentation are substrings of other event types, e.g., “COLD/WIND CHILL” is a substring of “EXTREME COLD/WIND CHILL”. Other overlapping weather event types include:

# Create all valid pair combinations
pairs <- expand.grid(i = seq_along(doc_event_types),
                     j = seq_along(doc_event_types)) %>%
  filter(i != j)

# Check if doc_event_types[i] is a substring of doc_event_types[j]
overlaps <- pairs %>%
  rowwise() %>%
  mutate(
    pattern = doc_event_types[i],
    text    = doc_event_types[j],
    is_match = str_detect(text, fixed(pattern))
  ) %>%
  ungroup() %>%
  filter(is_match) %>%
  transmute(
    shorter = pattern,
    longer  = text
  ) %>%
  distinct() %>%
  arrange(shorter, longer)

# Display all overlapping pairs
print(overlaps)

## # A tibble: 9 × 2
##   shorter           longer                  
##   <chr>             <chr>                   
## 1 COLD/WIND CHILL   EXTREME COLD/WIND CHILL 
## 2 FLOOD             COASTAL FLOOD           
## 3 FLOOD             FLASH FLOOD             
## 4 FLOOD             LAKESHORE FLOOD         
## 5 HAIL              MARINE HAIL             
## 6 HEAT              EXCESSIVE HEAT          
## 7 HIGH WIND         MARINE HIGH WIND        
## 8 STRONG WIND       MARINE STRONG WIND      
## 9 THUNDERSTORM WIND MARINE THUNDERSTORM WIND