This project aims to address below 2 questions:
Which types of events are most harmful to population health?
Which types of events have the greatest economic consequences?
As such, I intend to put together the event types with highest amount of fatalities and/or injuries as the most harmful contributors to population health. Then I use the event types with highest total damage to property & crop as those with greatest economic consequences.
The data is based on Storm Data measured from 1950 to November 2011. I tried testing the hypothesis that more recent data, specifically from 1993 till November 2011), is more complete & more reflective for the purpose of this project, but my findings suggest that such hypothesis is not reliable.
Throughout the data processing phase I found out that event types recorded in Storm Data are not named unanimously, which leads to hundreds of event types being report separately while they should be grouped up together. I used visual check & manually grouped them up together on best-effort basic.
library(dplyr)
library(lubridate)
url <- "https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2"
download.file(url,destfile = "Storm.data")
download.time <- Sys.time()
raw.data <- read.csv("storm.data")
Below steps added a “YEAR” column for future use:
transformed.data <- raw.data
transformed.data$BGN_DATE <- as.Date(transformed.data$BGN_DATE,format = "%m/%d/%Y")
transformed.data$YEAR <- year(transformed.data$BGN_DATE)
year.period <- max(year(transformed.data$BGN_DATE)) - min(year(transformed.data$BGN_DATE)) + 1
Below I tried transforming data into summary based on event type (EVTYPE column), together with fatalities, injuries, prop damage (PROPDMG), crop damage (CROPDMG column), & added number of occurrence throughout the record period (count), plus the yearly average occurrence (year.count):
event.data <- transformed.data %>% group_by(EVTYPE) %>% summarise(fatalities = sum(FATALITIES), injuries = sum(INJURIES), propdmg = sum(PROPDMG), cropdmg = sum(CROPDMG),count = n(),year.count = n()/year.period)
dim(event.data)
## [1] 985 7
View(event.data)
Here I notice there are too many event types (985 of them) & there are blanks before & after some EVTYPE values, so I try to trim the EVTYPE first
transformed.data$EVTYPE <- trimws(transformed.data$EVTYPE)
event.data <- transformed.data %>% group_by(EVTYPE) %>% summarise(fatalities = sum(FATALITIES), injuries = sum(INJURIES), propdmg = sum(PROPDMG), cropdmg = sum(CROPDMG),count = n(),year.count = n()/year.period)
dim(event.data)
## [1] 977 7
Even types got reduced to 977. Then since a lot of event types got inproperly separated. For example: “SNOW/ BITTER COLD”, “SNOW/BLOWING SNOW”, “SNOW/FREEZING RAIN”. So I try to spot the similar types & group them up a bit:
event.data$reduced.EVTYPE <- substr(event.data$EVTYPE,1,6)
event.data <- event.data %>% group_by(reduced.EVTYPE) %>% mutate (duplicate.index = n())
duplicate.check <- event.data %>% select(EVTYPE,reduced.EVTYPE,duplicate.index) %>% arrange(desc(duplicate.index))
print(duplicate.check)
## # A tibble: 977 × 3
## # Groups: reduced.EVTYPE [348]
## EVTYPE reduced.EVTYPE duplicate.index
## <chr> <chr> <int>
## 1 THUNDEERSTORM WINDS THUNDE 79
## 2 THUNDERESTORM WINDS THUNDE 79
## 3 THUNDERSNOW THUNDE 79
## 4 THUNDERSTORM THUNDE 79
## 5 THUNDERSTORM WINDS THUNDE 79
## 6 THUNDERSTORM DAMAGE THUNDE 79
## 7 THUNDERSTORM DAMAGE TO THUNDE 79
## 8 THUNDERSTORM HAIL THUNDE 79
## 9 THUNDERSTORM W INDS THUNDE 79
## 10 THUNDERSTORM WIND THUNDE 79
## # ℹ 967 more rows
View(duplicate.check)
The “duplicate.check” table is used to spot event types that should be grouped up. Due to the messiness of raw even types, I found no alternative but to group them up manually one by one, with priority given to event types with more occurrence & damage. Group up event types manually:
transformed.data$reduced.EVTYPE <- substr(transformed.data$EVTYPE,1,6)
# THUNDERSTORM
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(reduced.EVTYPE == "THUNDE","THUNDERSTORM",EVTYPE))
# RAIN
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("RAIN",transformed.data$newEVTYPE),"RAIN",newEVTYPE))
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("SHOWER",transformed.data$newEVTYPE),"RAIN",newEVTYPE))
# SNOW
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("SNOW",transformed.data$newEVTYPE),"SNOW",newEVTYPE))
# SURF
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("SURF",transformed.data$newEVTYPE),"SURF",newEVTYPE))
#COLD
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("COLD",transformed.data$newEVTYPE),"COLD",newEVTYPE))
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("LOW TEMP",transformed.data$newEVTYPE),"COLD",newEVTYPE))
#HEAT
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("HEAT",transformed.data$newEVTYPE),"HEAT",newEVTYPE))
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("HIGH TEMP",transformed.data$newEVTYPE),"HEAT",newEVTYPE))
#WIND
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("WIND",transformed.data$newEVTYPE),"WIND",newEVTYPE))
#FLOOD
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("FLOOD",transformed.data$newEVTYPE),"FLOOD",newEVTYPE))
#UNSEASONABLE WEATHER
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("^UNSEASONABLY",transformed.data$newEVTYPE),"UNSEASONABLE WEATHER",newEVTYPE))
#LIGHTNING
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("LIGHTNING",transformed.data$newEVTYPE),"LIGHTNING",newEVTYPE))
#TORNADO
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("TORNADO",transformed.data$newEVTYPE),"TORNADO",newEVTYPE))
#HURRICANE
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("HURRICANE",transformed.data$newEVTYPE),"HURRICANE",newEVTYPE))
#BLIZZARD
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("BLIZZARD",transformed.data$newEVTYPE),"BLIZZARD",newEVTYPE))
#RIP CURRENT
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("RIP CURRENT",transformed.data$newEVTYPE),"RIP CURRENT",newEVTYPE))
#FOG
transformed.data <- transformed.data %>% mutate(newEVTYPE = ifelse(grepl("FOG",transformed.data$newEVTYPE),"FOG",newEVTYPE))
Now I transform data into (new) event type again. The number of event type should have been reduced significantly:
event.data <- transformed.data %>% group_by(newEVTYPE) %>% summarise(fatalities = sum(FATALITIES), injuries = sum(INJURIES), propdmg = sum(PROPDMG), cropdmg = sum(CROPDMG),count = n(),year.count = n()/year.period)
dim(event.data)
## [1] 464 7
So event types got reduced to 464. Still insane imo, but (I assume) grouping event types is not core analysis here.
event.data <- mutate(event.data,totaldmg = cropdmg + propdmg)
event.data <- mutate(event.data,meanfatalities = fatalities/count)
event.data <- mutate(event.data,meaninjuries = injuries/count)
event.data <- mutate(event.data,meandmg = totaldmg/count)
Above steps added “totaldmg” column for total damage done to both crop & property (totaldmg), then the average damage per event (meandmg),average fatalities per event (meanfatalities), average injuries per event (meaninjuries)
top.10.fatalities <- head(event.data %>% arrange(desc(fatalities)),n=10)
top.10.injuries <- head(event.data %>% arrange(desc(injuries)),n=10)
top.10.dmg <- head(event.data %>% arrange(desc(totaldmg)),n=10)
In order to reveal the impact of the top 10 event types in term of percentage of total, I tried drawing below figure of 3 pie charts
par(mfrow = c(1,3),mar = c(4,1,4,1))
pie(
c(top.10.fatalities$fatalities,
sum(event.data[!event.data$newEVTYPE %in%
top.10.fatalities$newEVTYPE,"fatalities"])),
labels = c(top.10.fatalities$newEVTYPE,"OTHERS"),
main = "Fatalities by Event Type",
cex = 0.7
)
pie(
c(top.10.injuries$injuries,
sum(event.data[!event.data$newEVTYPE %in%
top.10.injuries$newEVTYPE,"injuries"])),
labels = c(top.10.injuries$newEVTYPE,"OTHERS"),
main = "Injuries by Event Type",
cex = 0.7
)
pie(
c(top.10.dmg$totaldmg,
sum(event.data[!event.data$newEVTYPE %in%
top.10.dmg$newEVTYPE,"totaldmg"])),
labels = c(top.10.dmg$newEVTYPE,"OTHERS"),
main = "Economic damage by Event Type",
cex = 0.7
)
Then for below sections of top 10 event types with most human harm &
greatest economic consequences per event I rule out
those event types with rare occurrence (count <10 for the record
period of 62 years) since those are either too rare to be significant
for this analysis or the name of event type is not matching with above
grouping exercise:
top.10.avg.fatalities <- head(event.data %>% filter(count > 9) %>% arrange(desc(meanfatalities)),n=10)
top.10.avg.injuries <- head(event.data %>% filter(count > 9) %>% arrange(desc(meaninjuries)),n=10)
top.10.avg.dmg <- head(event.data %>% filter(count > 9) %>% arrange(desc(meandmg)),n=10)
Creating a summary table of the findings in this section
summary <- data.frame(top.10.fatalities = top.10.fatalities$newEVTYPE)
summary$top.10.injuries = top.10.injuries$newEVTYPE
summary$top.10.dmg = top.10.dmg$newEVTYPE
summary$top.10.avg.fatalities = top.10.avg.fatalities$newEVTYPE
summary$top.10.avg.injuries = top.10.avg.injuries$newEVTYPE
summary$top.10.avg.dmg = top.10.avg.dmg$newEVTYPE
As mentioned in the project intro: In the earlier years of the database there are generally fewer events recorded, most likely due to a lack of good records. More recent years should be considered more complete. Here I’m just trying to see if I can find a cut-off time between low-quality data & high-quality data, by checking the amount of events recorded each year for the recording period from 1950 to 2011
year.count <- transformed.data %>% group_by(YEAR) %>% summarise(count = n())
plot(year.count$YEAR,year.count$count,type = "l")
The line chart suggests that somewhere in early 1990s may be a good cut-off time. Thus I check the “year.count” table & see a spike from 1993 (12,607 events recorded) to 1994 (20,631 events recorded). Therefore I try extracting the event summary with data after 1993 only:
recent.event.data <- transformed.data %>% filter(YEAR > 1993) %>% group_by(EVTYPE) %>% summarise(fatalities = sum(FATALITIES), injuries = sum(INJURIES), propdmg = sum(PROPDMG), cropdmg = sum(CROPDMG),count = n(),year.count = n()/year.period)
There are still 922 event types (compared to 977 in the table with all data) in this table, which rules out my hypothesis that recent data records better event types, so I give up on this hypothesis.
Here I think better to share the summary chart & table first, before moving on to provide comments on the findings. So here is the chart:
par(mfrow = c(1,3),mar = c(4,1,4,1))
pie(
c(top.10.fatalities$fatalities,
sum(event.data[!event.data$newEVTYPE %in%
top.10.fatalities$newEVTYPE,"fatalities"])),
labels = c(top.10.fatalities$newEVTYPE,"OTHERS"),
main = "Fatalities by Event Type",
cex = 0.7
)
pie(
c(top.10.injuries$injuries,
sum(event.data[!event.data$newEVTYPE %in%
top.10.injuries$newEVTYPE,"injuries"])),
labels = c(top.10.injuries$newEVTYPE,"OTHERS"),
main = "Injuries by Event Type",
cex = 0.7
)
pie(
c(top.10.dmg$totaldmg,
sum(event.data[!event.data$newEVTYPE %in%
top.10.dmg$newEVTYPE,"totaldmg"])),
labels = c(top.10.dmg$newEVTYPE,"OTHERS"),
main = "Economic damage by Event Type",
cex = 0.7
)
Then the summary table:
print(summary)
## top.10.fatalities top.10.injuries top.10.dmg top.10.avg.fatalities
## 1 TORNADO TORNADO TORNADO TSUNAMI
## 2 HEAT HEAT FLOOD HEAT
## 3 FLOOD WIND WIND RIP CURRENT
## 4 WIND FLOOD THUNDERSTORM AVALANCHE
## 5 LIGHTNING LIGHTNING HAIL HURRICANE
## 6 RIP CURRENT THUNDERSTORM LIGHTNING Cold
## 7 COLD ICE STORM SNOW GLAZE
## 8 AVALANCHE HAIL WINTER STORM MIXED PRECIP
## 9 WINTER STORM HURRICANE WILDFIRE COLD
## 10 THUNDERSTORM WINTER STORM RAIN ICY ROADS
## top.10.avg.injuries top.10.avg.dmg
## 1 GLAZE TYPHOON
## 2 TSUNAMI ICE
## 3 HURRICANE HURRICANE
## 4 HEAT Gusty Winds
## 5 MIXED PRECIP TROPICAL STORM
## 6 ICE STORM SURGE
## 7 BLACK ICE TORNADO
## 8 TORNADO STORM SURGE/TIDE
## 9 ICY ROADS SEICHE
## 10 DUST STORM TSUNAMI
# based on fatalities:
print(top.10.fatalities)
## # A tibble: 10 × 11
## newEVTYPE fatalities injuries propdmg cropdmg count year.count totaldmg
## <chr> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 TORNADO 5636 91407 3215748. 100027. 60698 979 3315775.
## 2 HEAT 3138 9154 3233. 1473. 2652 42.8 4706.
## 3 FLOOD 1522 8603 2432422. 362784. 82622 1333. 2795206.
## 4 WIND 1019 8999 1800038. 135463. 259592 4187. 1935501.
## 5 LIGHTNING 817 5232 603387. 3581. 15761 254. 606967.
## 6 RIP CURRENT 572 529 163 0 774 12.5 163
## 7 COLD 435 320 13967. 8910. 2432 39.2 22877.
## 8 AVALANCHE 224 170 1624. 0 386 6.23 1624.
## 9 WINTER STORM 206 1321 132721. 1979. 11433 184. 134700.
## 10 THUNDERSTORM 202 2453 1328911. 85735. 103719 1673. 1414645.
## # ℹ 3 more variables: meanfatalities <dbl>, meaninjuries <dbl>, meandmg <dbl>
# based on injuries:
print(top.10.injuries)
## # A tibble: 10 × 11
## newEVTYPE fatalities injuries propdmg cropdmg count year.count totaldmg
## <chr> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 TORNADO 5636 91407 3215748. 100027. 60698 979 3315775.
## 2 HEAT 3138 9154 3233. 1473. 2652 42.8 4706.
## 3 WIND 1019 8999 1800038. 135463. 259592 4187. 1935501.
## 4 FLOOD 1522 8603 2432422. 362784. 82622 1333. 2795206.
## 5 LIGHTNING 817 5232 603387. 3581. 15761 254. 606967.
## 6 THUNDERSTORM 202 2453 1328911. 85735. 103719 1673. 1414645.
## 7 ICE STORM 89 1975 66001. 1689. 2006 32.4 67690.
## 8 HAIL 15 1361 688693. 579596. 288661 4656. 1268290.
## 9 HURRICANE 133 1326 23757. 10803. 285 4.60 34560.
## 10 WINTER STORM 206 1321 132721. 1979. 11433 184. 134700.
## # ℹ 3 more variables: meanfatalities <dbl>, meaninjuries <dbl>, meandmg <dbl>
Here the top 5 event types are the same re. either fatalities or injuries as a metric for impact on population health. Down below we RIP CURRENT, AVALANCHE & COLD appearance in top 10 contributors to fatalities, not to injuries, & HURRICANE, ICE STORM and HAIL appearing in top 10 contributors to injuries, & not fatalities. Percentage-wise TORNADO has absolute impact on fatalities & injuries. HEAT also has absolute impact on fatalities.
print(top.10.dmg)
## # A tibble: 10 × 11
## newEVTYPE fatalities injuries propdmg cropdmg count year.count totaldmg
## <chr> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 TORNADO 5636 91407 3215748. 100027. 60698 979 3315775.
## 2 FLOOD 1522 8603 2432422. 362784. 82622 1333. 2795206.
## 3 WIND 1019 8999 1800038. 135463. 259592 4187. 1935501.
## 4 THUNDERSTORM 202 2453 1328911. 85735. 103719 1673. 1414645.
## 5 HAIL 15 1361 688693. 579596. 288661 4656. 1268290.
## 6 LIGHTNING 817 5232 603387. 3581. 15761 254. 606967.
## 7 SNOW 161 1158 150710. 2196. 17569 283. 152906.
## 8 WINTER STORM 206 1321 132721. 1979. 11433 184. 134700.
## 9 WILDFIRE 75 911 84459. 4364. 2761 44.5 88824.
## 10 RAIN 114 301 59376. 12922. 12209 197. 72298.
## # ℹ 3 more variables: meanfatalities <dbl>, meaninjuries <dbl>, meandmg <dbl>
The list here sees similar appearance in top 10 contributors to fatalities/injuries, with the new appearance of SNOW, WILDFIRE and RAIN. Percentage-wise the top 5 contributors compose a very big portion of total economic impact.
print(top.10.avg.fatalities)
## # A tibble: 10 × 11
## newEVTYPE fatalities injuries propdmg cropdmg count year.count totaldmg
## <chr> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 TSUNAMI 33 129 905. 20 20 0.323 925.
## 2 HEAT 3138 9154 3233. 1473. 2652 42.8 4706.
## 3 RIP CURRENT 572 529 163 0 774 12.5 163
## 4 AVALANCHE 224 170 1624. 0 386 6.23 1624.
## 5 HURRICANE 133 1326 23757. 10803. 285 4.60 34560.
## 6 Cold 3 0 54 0 10 0.161 54
## 7 GLAZE 7 216 311. 0 32 0.516 311.
## 8 MIXED PRECIP 2 26 0 0 10 0.161 0
## 9 COLD 435 320 13967. 8910. 2432 39.2 22877.
## 10 ICY ROADS 5 31 341. 0 28 0.452 341.
## # ℹ 3 more variables: meanfatalities <dbl>, meaninjuries <dbl>, meandmg <dbl>
print(top.10.avg.injuries)
## # A tibble: 10 × 11
## newEVTYPE fatalities injuries propdmg cropdmg count year.count totaldmg
## <chr> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 GLAZE 7 216 311. 0 32 0.516 311.
## 2 TSUNAMI 33 129 905. 20 20 0.323 925.
## 3 HURRICANE 133 1326 23757. 10803. 285 4.60 34560.
## 4 HEAT 3138 9154 3233. 1473. 2652 42.8 4706.
## 5 MIXED PRECIP 2 26 0 0 10 0.161 0
## 6 ICE 6 137 7660 0 61 0.984 7660
## 7 BLACK ICE 1 24 0 0 14 0.226 0
## 8 TORNADO 5636 91407 3215748. 100027. 60698 979 3315775.
## 9 ICY ROADS 5 31 341. 0 28 0.452 341.
## 10 DUST STORM 22 440 5050. 1602. 427 6.89 6651
## # ℹ 3 more variables: meanfatalities <dbl>, meaninjuries <dbl>, meandmg <dbl>
print(top.10.avg.dmg)
## # A tibble: 10 × 11
## newEVTYPE fatalities injuries propdmg cropdmg count year.count totaldmg
## <chr> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 TYPHOON 0 5 1.43e3 825 11 0.177 2254.
## 2 ICE 6 137 7.66e3 0 61 0.984 7660
## 3 HURRICANE 133 1326 2.38e4 10803. 285 4.60 34560.
## 4 Gusty Winds 0 1 9.38e2 0 10 0.161 938
## 5 TROPICAL STORM 58 340 4.84e4 5899. 690 11.1 54323.
## 6 STORM SURGE 13 38 1.94e4 5 261 4.21 19398.
## 7 TORNADO 5636 91407 3.22e6 100027. 60698 979 3315775.
## 8 STORM SURGE/TI… 11 5 6.78e3 850 148 2.39 7627.
## 9 SEICHE 0 0 9.8 e2 0 21 0.339 980
## 10 TSUNAMI 33 129 9.05e2 20 20 0.323 925.
## # ℹ 3 more variables: meanfatalities <dbl>, meaninjuries <dbl>, meandmg <dbl>
Based on the 3 lists here, & the summary table, top 10 damaging event types on per-event basic include extreme event types like TSUNAMI, RIP CURRENT, AVALANCHE, HURRICANE, etc. Then still TORNADO appears in both top 10 average injuries & top 10 average economic damage, presenting itself as among the most damaging event type not only at total level but at average level.