Please note: In the global options I’ve set ‘warning’, ‘message’, and ‘error’ to FALSE in this document. As a result, you will see all R code and selected results, but not the warnings or messages produced by running this code. I did this to enhance readability and save the reader’s time.

Synopsis

This analysis uses the NOAA Storm Database to consider the severe weather events across the United States that are 1) the most harmful for public health, and 2) the most economically damaging. The document outlines Data Processing steps including reading the data, consolidating closely related weather event types, transforming column values to calculate total values, and consolidating data by severe weather event type. This produces analyzable data, and the document describes an analysis approach to address both public health harms and economic damages. The report then conducts this analysis, providing associated figures and charts that show that Tornadoes cause the most harm to public health across the US and that Floods cause the most economic damage. In fact, both of these event types scored in the top 5 in both public health harms and economic damage, making them clearly the two most threatening severe weather events across the United States.

Data Processing

Ultimately, we want to answer two questions:

1. Across the United States, which types of events are most harmful with respect to population health?

2. Across the United States, which types of events have the greatest economic consequences?

Data processing will be tailored to exploring these two questions. The following code reads the data file into a tibble in R:

library(readr)
file <- "repdata_data_StormData.csv.bz2"
stormdata <- read_csv(bzfile(file))

Let’s take a look at what we’re working with here. What are the dimensions of the tibble?

## look at the dimensions of the data
dim(stormdata)
## [1] 902297     37

Ok so there are a lot of observations here. Over 900,000! Let’s see how many different severe weather event types are in the data set:

## get the count of unique sever weather event types
length(unique(stormdata$EVTYPE))
## [1] 977

And what kinds of event names are we looking at here? Let’s look at a sample:

head(unique(stormdata$EVTYPE), 20)
##  [1] "TORNADO"                   "TSTM WIND"                
##  [3] "HAIL"                      "FREEZING RAIN"            
##  [5] "SNOW"                      "ICE STORM/FLASH FLOOD"    
##  [7] "SNOW/ICE"                  "WINTER STORM"             
##  [9] "HURRICANE OPAL/HIGH WINDS" "THUNDERSTORM WINDS"       
## [11] "RECORD COLD"               "HURRICANE ERIN"           
## [13] "HURRICANE OPAL"            "HEAVY RAIN"               
## [15] "LIGHTNING"                 "THUNDERSTORM WIND"        
## [17] "DENSE FOG"                 "RIP CURRENT"              
## [19] "THUNDERSTORM WINS"         "FLASH FLOOD"

With over 900,000 observations and 977 different weather event types, we’re going to want to narrow this data set down in order to more effectively explore the research questions. The event names also suggest that there are some repeats that could be consolidated. Here’s a step by step plan for doing that cleaning:

  1. As the event names show, there are some event names that should probably be merged together, such as “Snow” and “Snow/Ice,” or “Hurricane Opal” and “Hurricane Erin.” We should consolidate and streamline the event type list.
  2. Because the PROPDMGEXP column specifies whether the PROPDMG column describes amounts in the thousands, millions, or billions, we’ll first mutate the PROPDMGEXP so that it represents the full dollar amount. The PROPDMGEXP column has other values than ‘k’, ‘m’, and ‘b’, but as these are not explained in the documentation and are not associated with any actual property damage values, we will treat these as irrelevant.
  3. Do the same for the CROPDMG column, based on CROPDMGEXP.
  4. Group the data by severe weather event type.
  5. Create columns that calculate the total Fatalities, Injuries, Property Damage, and Crop Damage per weather event type (that is, per group). At this step, we’ll also create a Frequency column showing the total number of observations for each weather event type.

The result will be a new tibble that has a consolidated list of event types, and columns for Event Type, Total Fatalities, Total Injuries, Total Property Damage, and Total Crop Damage. These are the values we’ll use in order to answer our research questions.

First, we’ll try to clean and streamline the data a bit, especially the EVTYPE variable. With 977 different entries, we won’t be able to be comprehensive here. But we can use grepl() and regex coding to at least identify major themes and convert over-detailed event names to uniform categories. For example, “Hurricane Opal” and “Hurricane Erin” should both become HURRICANE; and “Heat” and “Excessive Heat” can both become “Heat”. In many cases, the EVTYPE value combines multiple severe whether events, such as “thunderstorm winds” or “snow/rain/hail.” For these cases, we’ll develop a hierarchy of damaging/harmful weather events that allows us to choose which event type prevails. That hierarchy is:

  1. Blizzards
  2. Snow
  3. Heat
  4. Flood
  5. Tornado
  6. Lightning
  7. Hail
  8. Rain
  9. Cold
  10. Fire
  11. Freeze
  12. Windchill
  13. Hurricane
  14. High Surf
  15. Thunderstorm

This is not to say that this hierarchy is true in all cases, and there are certainly assumptions made here. Nonetheless, it’s an effective way of ordering weather events that are likely to overlap or be combined in EVTYPE names in order to consolidate meaningful categories of events. Importantly, this hierarchy is not designed to rank events by most harmful/damaging. Rather it is designed to ensure that in cases of likely overlap or combination, the more harfmul/damaging of the two will be ranked higher.

Here’s the code that applies these transformations to the data:

library(dplyr)
stormdata_cleaner <- stormdata %>%
     ## change entries to upper case to avoid inconsistencies
     mutate(EVTYPE = toupper(EVTYPE)) %>%
     ## use regex to consolidate related weather event types
     mutate(EVTYPE = ifelse(grepl("BLIZZARD", EVTYPE), "BLIZZARD", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("SNOW", EVTYPE), "SNOW", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("HEAT", EVTYPE), "HEAT", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("FLOOD", EVTYPE), "FLOOD", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("TORNADO", EVTYPE), "TORNADO", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("LIGHT", EVTYPE), "LIGHTNING", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("HAIL", EVTYPE), "HAIL", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("RAIN", EVTYPE), "RAIN", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("COLD", EVTYPE), "COLD", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("FIRE", EVTYPE), "FIRE", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("FREEZE", EVTYPE), "FREEZE", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("CHILL", EVTYPE), "WCHILL", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("HURRICANE", EVTYPE), "HURRICANE", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("WIND", EVTYPE), "HIGH WINDS", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("SURF", EVTYPE), "HIGH SURF", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("THUNDER", EVTYPE), "THUNDERSTORM", EVTYPE)) %>%
     mutate(EVTYPE = ifelse(grepl("TSTM", EVTYPE), "THUNDERSTORM", EVTYPE))

Now that we’ve got a consolidated list of event types, here’s the code that calculates total property damages by applying the information in the PROPDMGEXP column to the value in the PRODMG column:

storms <- stormdata_cleaner %>%
     ## standardize PROPDMGEXP values to lower case
     mutate(PROPDMGEXP = tolower(PROPDMGEXP)) %>%
     ##apply calculations to only the cases where PROPDMGEXP is k, m, or b
     mutate(prop_damage = case_when(
          grepl("[kmb]", PROPDMGEXP) & PROPDMGEXP == 'k' ~ PROPDMG*1000,
          grepl("[kmb]", PROPDMGEXP) & PROPDMGEXP == 'm' ~ PROPDMG*1000000,
          grepl("[kmb]", PROPDMGEXP) & PROPDMGEXP == 'b' ~ PROPDMG*1000000000,
          TRUE ~ PROPDMG))

Now we’ll do the same with the CROPDMG and CROPDMGEXP columns. Here’s the code:

storms<- storms %>%
     ## standardize CROPDMGEXP values to lower case
     mutate(CROPDMGEXP = tolower(CROPDMGEXP)) %>%
     ##apply calculations to only the cases where CROPCMGEXP is k, m, or b
     mutate(crop_damage = case_when(
          grepl("[kmb]", CROPDMGEXP) & CROPDMGEXP == 'k' ~ CROPDMG*1000,
          grepl("[kmb]", CROPDMGEXP) & CROPDMGEXP == 'm' ~ CROPDMG*1000000,
          grepl("[kmb]", CROPDMGEXP) & CROPDMGEXP == 'b' ~ CROPDMG*1000000000,
          TRUE ~ CROPDMG))

Now we’ll group the data by severe weather event:

## group data by severe weather event type
storms_grouped <- storms %>%
     group_by(EVTYPE)

And here’s the code that calculates total values for fatalities, injuries, and economic damage (defined as property damage + crop damage):

## create calculated columns for total values per group (that is, per weather event type)
## this also creates a 'frequency' column that shows the number of observations per group
storms_tidy <- storms_grouped %>%
     summarize(
          frequency = n(),
          total_fatalities = sum(FATALITIES),
          total_injuries = sum(INJURIES),
          total_damage = sum(prop_damage+crop_damage),
          )

## and we'll rename the EVTYPE column to event_type for reading ease
storms_tidy <- storms_tidy %>%
     rename(event_type = EVTYPE)

Now we have a cleaned, processed, and tidy data set ‘storms_tidy’ that gives us all the information we need for our analysis. Let’s take a quick look:

dim(storms_tidy)
## [1] 317   5

We see that the tibble storms_tidy has 317 rows- one for each consolidated weather event type in the data set- and just 5 columns. We can look at what those are and how they’ll guide our analysis:

head(storms_tidy)
## # A tibble: 6 × 5
##   event_type             frequency total_fatalities total_injuries total_damage
##   <chr>                      <int>            <dbl>          <dbl>        <dbl>
## 1 ?                              1                0              0         5000
## 2 ABNORMAL WARMTH                4                0              0            0
## 3 ABNORMALLY DRY                 2                0              0            0
## 4 ABNORMALLY WET                 1                0              0            0
## 5 APACHE COUNTY                  1                0              0         5000
## 6 ASTRONOMICAL HIGH TIDE       103                0              0      9425000

And that’s it for data processing. On to analysis!

Analysis Approach

So now we have to define how we will analyze the following two questions:

1. Across the United States, which types of events are most harmful with respect to population health?

2. Across the United States, which types of events have the greatest economic consequences?

We’ve already made a couple decisions here:

We could use mean values instead of totals here, but that seems less relevant than totals for answering the question. If our questions were seeking to know the events that on average cause the most public health and economic risk, respectively, then means would make more sense. But here our questions imply an interest in overall or total harms.

Defining Public Health Harms

There are still more decisions to make, though. How do we define public health harms when we know both fatalities and injuries? We could combine them into a summed ‘casualties’ metric, but that seems inappropriate as fatalities are such a more drastic harm than injury. Though it’s a bit coarse and arbitrary, we’ll create a weighting system that defines fatalities as twice as harmful as injuries. We’ll use min-max scaling to set a “Health Harms Score” on a 0-1 scale. The most harmful event will therefore have a score of 1, and the least harmful will have a score of 0. All events in between will be scaled proportionally, relative to the min and max values.

Here’s the code that accomplishes this:

## Define the weighting scale, with fatalities twice as harmful as injuries
weights <- c(
     total_fatalities = 0.666667,
     total_injuries = 0.333333
)

## run the normalize function on the fatalities and injuries values
health_events <- storms_tidy %>%
     mutate(
          scaled_fatalities = (total_fatalities - min(total_fatalities)) / (max(total_fatalities) - min(total_fatalities)),
          scaled_injuries = (total_injuries - min(total_injuries)) / (max(total_injuries) - min(total_injuries))
     ) %>%
## calculate composite score using scaled values
     mutate(health_score = (scaled_fatalities * weights["total_fatalities"]) +
                 (scaled_injuries * weights["total_injuries"])
            ) %>%
# order event types by composite score and select relevant columns
     arrange(desc(health_score)) %>%
     select(event_type, frequency, health_score, total_fatalities, total_injuries)

Defining Economic Consequences

This one’s easier, as we already have a ‘total_damages’ column that combines property damages with economic damages. So that’s how we’ll define economic damages. As we did with the public health analysis, we’ll create a normalized 0-1 score for the dollar amount damages, because that’s helpful for seeing the distance between damage values. Both measures are in USD, so we don’t need to set a weighting scale. We’ll create a separate table for economic damages with those values and other relevant variables from the data set.

# Define Min-max normalization function
normalize <- function(x) {
     (x - min(x)) / (max(x) - min(x))
}
#create object with prop_damage and crop_damage values by group
econ_grouped <- storms_grouped %>%
                    summarize(prop_damage = sum(prop_damage),
                              crop_damage = sum(crop_damage))

## create economic damage table ordered from highest to lowest damage
econ_events <- storms_tidy %>%
     mutate(econ_score = normalize(total_damage))%>%
     arrange(desc(econ_score)) %>%
     select(event_type, 
            frequency, 
            econ_score, 
            total_damage) %>%
## join the prop_damage and crop_damage values from the econ_grouped table
     right_join(
          select(econ_grouped, EVTYPE, prop_damage, crop_damage),
          by = c("event_type" = "EVTYPE"))

Results

Public Health Harms

Let’s take a look at the results of the health analysis. These are the top five events in the table ranked by our Health Harms Score:

# View the ranked event types
head(health_events, 5)
## # A tibble: 5 × 5
##   event_type frequency health_score total_fatalities total_injuries
##   <chr>          <int>        <dbl>            <dbl>          <dbl>
## 1 TORNADO        60700        1                 5661          91407
## 2 HEAT            2648        0.403             3138           9224
## 3 FLOOD          82725        0.211             1525           8604
## 4 HIGH WINDS    361884        0.180             1175          11314
## 5 LIGHTNING      15804        0.115              817           5232

We see that based on our scoring system that considers the total fatalities and injuries while supposing fatalities are twice as harmful as injuries, Tornadoes are by far the most harmful severe weather event for public health across the United States. We can learn more about these events:

freq <- health_events[1,2]
fatalities <- health_events[1,4]
injuries <- health_events[1,5]
casualties <- sum(health_events[1,4]+health_events[1,5])

print(casualties)
## [1] 97068
print(fatalities)
## # A tibble: 1 × 1
##   total_fatalities
##              <dbl>
## 1             5661
print(injuries)
## # A tibble: 1 × 1
##   total_injuries
##            <dbl>
## 1          91407
print(freq)
## # A tibble: 1 × 1
##   frequency
##       <int>
## 1     60700

Tornadoes have caused 97068 total casualties- 5661 fatalities and 91407 injuries. And they have occurred quite frequently, with 60700 observations of Tornadoes in this data set.

Let’s take a look at these data points plotted out to see if the Health Harms Score really is correlating with higher public health harms. For these plots, we’ll define public health harms as Total Casualties, which is the sum of fatalities and injuries. And we’ll separate the 5 event types with the highest Health Harms Score (in red) from all the other events (in black):

library(ggplot2)

## create variable to isolate the top 5 elements
health_events <- health_events %>%
  mutate(top_5 = ifelse(rank(-health_score) <= 5, "Top 5 Most Harmful Events", "All Others"))

# Create two-panel plot comparing Casualties with Health Score, with Top5 events disaggregated
g <- ggplot(health_events, aes(x = total_fatalities + total_injuries, y = health_score, color = top_5, shape = top_5)) +
  geom_point(size = 3) +
  geom_smooth(method = "lm") +
  labs(
    x = "Total Casualties (Fatalities + Injuries)",
    y = "Health Harms Score",
    title = "Relationships Between Health Risk Score and Total Casualties",
    shape = "Events",
    color = "Top 5 Health Scores"
  ) +
  scale_color_manual(
    values = c("Top 5 Most Harmful Events" = "red", "All Others" = "black"),
    labels = c("Top 5 Most Harmful Events", "All Others")
  ) +
  facet_wrap(~top_5, ncol = 1, scales = "free_x") +
  guides(color = FALSE, shape = FALSE) +
  coord_cartesian(ylim = c(0, 1.3))

print(g)

These plots evidence the utility of the Health Harms Score for identifying the most harmful weather events for public health across the US. At low scores, there’s no particular relationship between Health Risk Score and Total Casualties. That’s because casualties are so low that these events simply aren’t very harmful. Hence the low score.

But the bottom chart, in red, tells a different story. This isolates the 5 highest-scoring weather events. Notice how much higher the scale of casualties is along the x-axis for this plot than for the top chart. That’s how much more harmful these 5 weather events are than the others. And even within this group, we see a wide variation, with Tornadoes far and away the highest cause of fatalities.

And that’s our answer to the first question: based on this data, Tornadoes are the most harmful severe weather event for public health across the U.S.

Economic Consequences

Ok, so let’s take a look at the most economically damaging events. These are the top five rows of the econ_events table we created, which ranks the events by their overal economic damages in USD, combining property and crop damages:

# View the ranked event types
head(econ_events, 5)
## # A tibble: 5 × 6
##   event_type  frequency econ_score  total_damage   prop_damage crop_damage
##   <chr>           <int>      <dbl>         <dbl>         <dbl>       <dbl>
## 1 FLOOD           82725      1     179903280032. 167523190932. 12380089100
## 2 HURRICANE         288      0.502  90271472810   84756180010   5515292800
## 3 TORNADO         60700      0.328  59010559549.  58593098029.   417461520
## 4 STORM SURGE       261      0.241  43323541000   43323536000         5000
## 5 HAIL           290399      0.106  19131703945.  16019991072.  3111712873

We see that based on our scoring system that considers both the property and crop damages and weights them equally, Floods are by far the most economically damaging severe weather event in the United States. Let’s learn a bit more about these floods:

freq <- econ_events[1,2]
damage <- econ_events[1,4]

print(freq)
## # A tibble: 1 × 1
##   frequency
##       <int>
## 1     82725
print (damage)
## # A tibble: 1 × 1
##    total_damage
##           <dbl>
## 1 179903280032.

There have been 82725 floods across the U.S., collectively resulting in $179,724,664,032 in economic damages. That’s over 179 billion dollars!

Let’s take a look at these totals side by side, and see how much of these damages are coming from property damages vs. crop damages. We’ll look at the top 5 most economically damaging event types:

library(tidyr)
## isolate top 5 economically damaging events ordered by score and put in long format
top5_econ <- econ_events[1:5,] %>% 
     arrange(desc(econ_score)) %>%
     mutate(event_type = reorder(event_type, -econ_score)) %>%
     pivot_longer(cols = c(prop_damage, crop_damage),
                  names_to = "Damage_Type",
                  values_to = "Damages")

## create stacked bar chart showing prop_damage and crop_damage
ggplot(top5_econ, aes(x = event_type, y = Damages/1000000000, fill = Damage_Type)) +
     geom_bar(stat = "identity") +
     labs(title= "Total Economic Damages for the Top 5 Most Damaging Severe Weather Events",
          x = "Weather Event Type",
          y = "Damages (Billions of USD)",
          fill = "Damage Type") +
     theme(axis.text.x = element_text(angle = 90, hjust = 1))+
     scale_fill_manual(values = c("prop_damage" = "blue", "crop_damage" = "turquoise"),
                         labels = c("Crop Damage", "Property Damage"))

So we see that 1) Floods are far and away the most economically damaging severe weather event; and 2) across all of the top 5 most economically damaging events, property damage is the major source of economic damages.

So there we have our answers: Across the United States, Tornadoes are easily the most harmful severe weather event for public health, and Floods are easily the most economically damaging. Remember, Floods was also one of the top 5 most harmful event types for public health; and Tornadoes turned out to be one of the top 5 most economically damaging event types. So these two severe weather events are definitely the ones to pay attention to across the U.S.