Code in support of our paper submitted to PLoS ONE.
date()
## [1] "Thu May 15 16:16:08 2014"
Set working directory and load packages.
setwd("~/Dropbox/Tornadoes")
require(MASS)
require(RColorBrewer)
require(ggplot2)
require(rgdal)
require(quantreg)
require(plyr)
require(reshape)
require(spatstat)
require(maptools)
require(ordinal)
library(rgeos)
multiplot <- function(..., plotlist = NULL, file, cols = 1, layout = NULL) {
require(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel ncol: Number of columns of plots nrow: Number of rows
# needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols,
nrow = ceiling(numPlots/cols), byrow = TRUE)
}
if (numPlots == 1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col))
}
}
}
From the Storm Prediction Center (SPC) http://www.spc.noaa.gov/gis/svrgis/
download.file("http://myweb.fsu.edu/jelsner/data/tornado2011.zip", "tornado2011.zip",
mode = "wb")
unzip("tornado2011.zip")
Torn = readOGR(dsn = ".", layer = "tornado")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "tornado"
## with 56221 features and 28 fields
## Feature type: wkbPoint with 2 dimensions
Read the city data. City center locations.
download.file("http://myweb.fsu.edu/jelsner/data/ci08au12.zip", "ci08au12.zip",
mode = "wb")
unzip("ci08au12.zip")
ll = "+proj=longlat +ellps=GRS80 +datum=NAD83"
US = readOGR(dsn = ".", layer = "ci08au12", p4s = ll)
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "ci08au12"
## with 43603 features and 23 fields
## Feature type: wkbPoint with 2 dimensions
Cities = spTransform(US, CRS(proj4string(Torn)))
Subset for F-scale greater than 0.
Torn2 = subset(Torn, Torn$FSCALE > 0)
Torn2.df = as.data.frame(Torn2)
Add a distance-from-nearest city center column to the data frame. First create the function using a ppp object on the city centers and then apply the function to the set tornado touchdown locations.
Cities.ppp = as.ppp(Cities)
f = distfun(Cities.ppp)
xy = coordinates(Torn2)
dd = f(xy[, 1], xy[, 2])
## Warning: data contain duplicated points
Torn2.df$dd = dd
Occurrence rates
trpyWide = ddply(Torn2.df[Torn2.df$YEAR >= 1994, ], .(YEAR), summarize, F1 = sum(FSCALE ==
1), F2 = sum(FSCALE == 2), F3 = sum(FSCALE == 3), F4 = sum(FSCALE == 4),
F5 = sum(FSCALE == 5), F1F5 = sum(FSCALE >= 1), F2F5 = sum(FSCALE >= 2),
F3F5 = sum(FSCALE >= 3), F4F5 = sum(FSCALE >= 4))
trpyLong = melt(trpyWide, id = "YEAR", measure.vars = c("F1F5", "F2F5", "F3F5",
"F4F5"))
trpyLong$"F1-F5" = trpyLong$F1F5
trpyLong$"F2-F5" = trpyLong$F2F5
trpyLong$"F3-F5" = trpyLong$F3F5
trpyLong$"F4-F5" = trpyLong$F4F5
p1 = ggplot(trpyLong, aes(x = YEAR, y = value)) + geom_point() + facet_grid(variable ~
., scale = "free_y") + geom_smooth(method = "lm", span = 1, color = "gray") +
theme_bw() + ylab("Number of Tornado Reports (1994-2011)") + xlab("")
p1
FIGURE 1
Fraction of all tornadoes since 1994 that are violent.
nrow(Torn2.df[Torn2.df$YEAR >= 1994 & Torn2.df$FSCALE >= 4, ])/nrow(Torn2.df[Torn2.df$YEAR >=
1994, ])
## [1] 0.01529
Create an ordered factor for the EF category. Convert length in miles and the width in yards to meters.
Torn2.df$EF = factor(Torn2.df$FSCALE, ordered = TRUE)
Torn2.df$Length = Torn2.df$LENGTH * 1609.34
Torn2.df$Width = Torn2.df$WIDTH * 0.9144
Torn2.df$Area = 0.5 * Torn2.df$Length * 0.5 * Torn2.df$Width * pi
Plot univariate and bivariate relationships between EF category and path length/width
p2a = ggplot(Torn2.df[Torn2.df$YEAR >= 1994, ], aes(x = factor(FSCALE), y = Length/1000)) +
geom_boxplot() + theme_bw() + coord_flip() + ylab("Path Length (km)") +
xlab("EF Category")
p2b = ggplot(Torn2.df[Torn2.df$YEAR >= 1994, ], aes(x = factor(FSCALE), y = Width/1000)) +
geom_boxplot() + theme_bw() + coord_flip() + ylab("Path Width (km)") + xlab("EF Category")
Subset on year and compute table of frequency by EF category.
Torn3.df = subset(Torn2.df, YEAR >= 1994)
ddply(Torn3.df, ~FSCALE, summarize, avgL = mean(Length/1000), q25L = quantile(Length/1000,
prob = 0.25), q50L = quantile(Length/1000, prob = 0.5), q75L = quantile(Length/1000,
prob = 0.75), q25W = quantile(Width, prob = 0.25), q50W = quantile(Width,
prob = 0.5), q75W = quantile(Width, prob = 0.75), avgA = mean(Area/1e+06),
q25A = quantile(Area/1e+06, prob = 0.25), q50A = quantile(Area/1e+06, prob = 0.5),
q75A = quantile(Area/1e+06, prob = 0.75), nr = length(YEAR))
## FSCALE avgL q25L q50L q75L q25W q50W q75W avgA
## 1 1 6.238 1.609 3.524 8.047 45.72 91.44 137.2 0.8437
## 2 2 13.047 3.541 8.851 17.027 91.44 182.88 365.8 4.0250
## 3 3 26.891 9.897 20.921 35.405 228.60 402.34 804.7 13.8257
## 4 4 42.579 16.093 27.286 56.729 443.48 754.38 1063.0 36.9600
## 5 5 67.296 45.512 58.950 65.935 1207.01 1307.59 1609.3 84.9030
## q25A q50A q75A nr
## 1 0.05779 0.2312 0.7108 6179
## 2 0.34673 1.2714 3.5471 1843
## 3 2.17795 6.9347 16.2733 543
## 4 5.12125 16.5447 37.7516 120
## 5 39.98062 56.6007 99.8938 13
Bubble scatter plot
TornO.df = Torn3.df[order(Torn3.df$FSCALE), ]
p2c = ggplot(TornO.df, aes(x = Length, y = Width, color = EF)) + geom_point(alpha = 1) +
scale_x_log10() + scale_y_log10() + scale_color_brewer(palette = "YlOrRd") +
theme_bw() + xlab("Path Length (m)") + ylab("Path Width (m)") + guides(color = guide_legend(reverse = TRUE))
p2c = ggplot(TornO.df, aes(x = Length, y = Width, size = EF)) + geom_point(alpha = 0.25) +
scale_x_log10() + scale_y_log10() + theme_bw() + xlab("Path Length (m)") +
ylab("Path Width (m)") + guides(size = guide_legend(reverse = TRUE))
source("multiplot.txt")
## Warning: cannot open file 'multiplot.txt': No such file or directory
## Error: cannot open the connection
mat = matrix(c(1, 2, 3, 3), nrow = 2, byrow = TRUE)
p2a = p2a + ggtitle("a") + theme(plot.title = element_text(hjust = 0))
p2b = p2b + ggtitle("b") + theme(plot.title = element_text(hjust = 0))
p2c = p2c + ggtitle("c") + theme(plot.title = element_text(hjust = 0))
multiplot(p2a, p2b, p2c, layout = mat)
## Loading required package: grid
FIGURE 2
Measurement precision
head(Torn3.df$LENGTH[Torn3.df$YEAR == 2007])
## [1] 15.07 1.83 1.20 5.30 3.00 1.40
head(Torn3.df$LENGTH[Torn3.df$YEAR == 2006])
## [1] 1.0 1.0 0.1 8.4 6.4 0.6
Create new covariates
Torn3.df$SLength = sqrt(Torn3.df$Length)
Torn3.df$SWidth = sqrt(Torn3.df$Width)
Torn3.df$RWidth = Torn3.df$Width^(3/4)
Torn3.df$RLength = Torn3.df$Length^(0.6)
Torn3.df$RwLength = Torn3.df$Length^(1/3)
Torn3.df$FAT = Torn3.df$FATALITIES > 0
Torn3.df$F2EF = Torn3.df$YEAR >= 2007
Scale the variables
Torn3.df$SLengthS = as.numeric(scale(Torn3.df$SLength))
Torn3.df$SWidthS = as.numeric(scale(Torn3.df$SWidth))
Torn3.df$RWidthS = as.numeric(scale(Torn3.df$RWidth))
Torn3.df$RLengthS = as.numeric(scale(Torn3.df$RLength))
Torn3.df$RwLengthS = as.numeric(scale(Torn3.df$RwLength))
Torn3.df$YEARS = as.numeric(scale(Torn3.df$YEAR))
Torn3.df$ddS = as.numeric(scale(Torn3.df$dd))
Create weights
EF = Torn3.df$EF
Torn3.df$Weights = (length(EF)/table(EF))[match(EF, levels(EF))]/5
Determine optimal model
model1 = clm(EF ~ RLengthS * RWidthS + ddS, scale = ~RWidthS + RLengthS, data = Torn3.df,
link = "logit")
model2 = clm(EF ~ RwLengthS * SWidthS, weights = Weights, scale = ~SWidthS +
RwLengthS, data = Torn3.df, link = "logit")
model3 = clm(EF ~ RwLengthS * SWidthS, nominal = ~SWidthS, weights = Weights,
scale = ~RwLengthS, data = Torn3.df, link = "logit")
model4 = clm(EF ~ RwLengthS * SWidthS, nominal = ~RwLengthS, weights = Weights,
scale = ~SWidthS, data = Torn3.df, link = "logit")
AIC(model1, model2, model3, model4)
## df AIC
## model1 10 11290
## model2 9 20143
## model3 11 20149
## model4 11 20048
Sensitivity of coefficients to distance from nearest city, fatalities, and year.
m0 = summary(model2)
AIC(model2)
## [1] 20143
L0 = m0$coefficients[5, 1]
W0 = m0$coefficients[6, 1]
LW0 = m0$coefficients[7, 1]
sepL = 2 * m0$coefficients[5, 2]/L0 * 100
sepW = 2 * m0$coefficients[6, 2]/W0 * 100
sepLW = 2 * m0$coefficients[7, 2]/abs(LW0) * 100
sepL
## [1] 7.299
sepW
## [1] 5.12
sepLW
## [1] 21.41
AIC(model2)
## [1] 20143
mdd.fit = clm(EF ~ RwLengthS * SWidthS + ddS, weights = Weights, scale = ~SWidthS +
RwLengthS, data = Torn3.df, link = "logit")
mdd = summary(mdd.fit)
(1 - exp(mdd$coefficients[7, 1])) * 100
## [1] 12.46
xx = confint(mdd.fit, type = "Wald")
(1 - exp(xx[7, 1])) * 100
## [1] 16.14
(1 - exp(xx[7, 2])) * 100
## [1] 8.629
pcL = abs(L0 - mdd$coefficients[5, 1])/abs(L0) * 100
pcW = abs(W0 - mdd$coefficients[6, 1])/abs(W0) * 100
pcLW = abs(LW0 - mdd$coefficients[8, 1])/abs(LW0) * 100
pcL
## [1] 0.1618
pcW
## [1] 1.415
pcLW
## [1] 0.4294
AIC(model2) - AIC(mdd.fit)
## [1] 35.86
mYR.fit = clm(EF ~ RwLengthS * SWidthS + F2EF, weights = Weights, scale = ~SWidthS +
RwLengthS, data = Torn3.df, link = "logit")
mYR = summary(mYR.fit)
(1 - exp(mYR$coefficients[7, 1])) * 100
## [1] 32.99
xx = confint(mYR.fit, type = "Wald")
(1 - exp(xx[7, 1])) * 100
## [1] 39.02
(1 - exp(xx[7, 2])) * 100
## [1] 26.36
pcL = abs(L0 - mYR$coefficients[5, 1])/abs(L0) * 100
pcW = abs(W0 - mYR$coefficients[6, 1])/abs(W0) * 100
pcLW = abs(LW0 - mYR$coefficients[8, 1])/abs(LW0) * 100
pcL
## [1] 1.81
pcW
## [1] 3.908
pcLW
## [1] 14.2
AIC(model2) - AIC(mYR.fit)
## [1] 71.99
mFAT.fit = clm(EF ~ RwLengthS * SWidthS + FAT, weights = Weights, scale = ~SWidthS +
RwLengthS, data = Torn3.df, link = "logit")
mFAT = summary(mFAT.fit)
pcL = abs(L0 - mFAT$coefficients[5, 1])/abs(L0) * 100
pcW = abs(W0 - mFAT$coefficients[6, 1])/abs(W0) * 100
pcLW = abs(LW0 - mFAT$coefficients[8, 1])/abs(LW0) * 100
pcL
## [1] 7.58
pcW
## [1] 17.47
pcLW
## [1] 34.78
AIC(model2) - AIC(mFAT.fit)
## [1] 943.4
Width and length conditional on fatality
summary(Torn3.df$Width)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 46 91 197 206 4020
summary(subset(Torn3.df, I(FATALITIES > 0))$Width)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 183 393 549 805 4020
summary(Torn3.df$Length)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 1610 4830 9560 11800 240000
summary(subset(Torn3.df, I(FATALITIES > 0))$Length)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 290 9620 20900 31000 40200 240000
Trends in annual mean path dimensions
std = function(x) sd(x)/sqrt(length(x))
dat0 = ddply(Torn3.df, .(YEAR), summarize, AvgLength = mean(Length), AvgWidth = mean(Width),
seLength = std(Length), seWidth = std(Width))
dat0$loLength = dat0$AvgLength - 1.96 * dat0$seLength
dat0$hiLength = dat0$AvgLength + 1.96 * dat0$seLength
dat0$loWidth = dat0$AvgWidth - 1.96 * dat0$seWidth
dat0$hiWidth = dat0$AvgWidth + 1.96 * dat0$seWidth
mat = matrix(c(1, 2), nrow = 2, byrow = TRUE)
pAvgL = ggplot(dat0, aes(x = YEAR, y = AvgLength/1000)) + geom_point() + geom_smooth(method = lm) +
geom_linerange(aes(ymax = hiLength/1000, ymin = loLength/1000)) + ylab("Annual Mean Path Length (km)") +
xlab("Year") + theme_bw()
pAvgW = ggplot(dat0, aes(x = YEAR, y = AvgWidth)) + geom_point() + geom_smooth(method = lm) +
geom_linerange(aes(ymax = hiWidth, ymin = loWidth)) + ylab("Annual Mean Path Width (m)") +
xlab("Year") + theme_bw()
pAvgL = pAvgL + ggtitle("a") + theme(plot.title = element_text(hjust = 0))
pAvgW = pAvgW + ggtitle("b") + theme(plot.title = element_text(hjust = 0))
multiplot(pAvgL, pAvgW, layout = mat)
FIGURE 3
summary(lm(AvgLength ~ YEAR, data = dat0))
##
## Call:
## lm(formula = AvgLength ~ YEAR, data = dat0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2382 -960 294 984 2188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -360465.2 130735.9 -2.76 0.014 *
## YEAR 184.7 65.3 2.83 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1440 on 16 degrees of freedom
## Multiple R-squared: 0.333, Adjusted R-squared: 0.292
## F-statistic: 8 on 1 and 16 DF, p-value: 0.0121
summary(lm(AvgWidth ~ YEAR, data = dat0))
##
## Call:
## lm(formula = AvgWidth ~ YEAR, data = dat0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.65 -19.61 -0.31 28.99 41.38
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -12055.73 2797.52 -4.31 0.00054 ***
## YEAR 6.11 1.40 4.38 0.00047 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.8 on 16 degrees of freedom
## Multiple R-squared: 0.545, Adjusted R-squared: 0.516
## F-statistic: 19.2 on 1 and 16 DF, p-value: 0.000469
Add YEAR as a covariate to account for the increasing survey efforts over time.
model.Final = clm(EF ~ RwLengthS * SWidthS + FAT + YEARS, weights = Weights,
scale = ~SWidthS + RwLengthS, data = Torn3.df, link = "logit")
Final = summary(model.Final)
clmpred = list(model.Final = NULL)
clmpred[[1]] = as.character(predict(model.Final, type = "class")[[1]])
stats = lapply(clmpred, function(x) {
W <- sapply(-2:2, function(j) sapply(1:5, function(i) sum(x == i + j & Torn3.df$EF ==
i)/sum(Torn3.df$EF == i)))
dimnames(W) = list(Observed = c("EF1", "EF2", "EF3", "EF4", "EF5"), predicted = paste("OBS",
-2:2, sep = ""))
100 * round(W, 3)
})
stats
## $model.Final
## predicted
## Observed OBS-2 OBS-1 OBS0 OBS1 OBS2
## EF1 0.0 0.0 69.8 25.6 4.2
## EF2 0.0 29.5 43.9 20.3 5.9
## EF3 5.2 24.9 36.8 24.7 8.5
## EF4 13.3 25.0 32.5 28.3 0.0
## EF5 7.7 15.4 76.9 0.0 0.0
Model adequacy
pred = as.numeric(predict(model.Final, type = "class")$fit)
Torn4.df = Torn3.df
Torn4.df$obs = Torn4.df$FSCALE
Torn4.df$pred = pred
Torn4.df$dif = Torn4.df$obs - Torn4.df$pred
p4 = ggplot(Torn4.df, aes(x = factor(dif))) + geom_histogram(binwidth = 0.5) +
geom_histogram(data = subset(Torn4.df, dif == 0), binwidth = 1, fill = "grey") +
facet_wrap(~YEAR) + xlab("Observed - Predicted EF Category") + ylab("Number of Tornadoes") +
theme_bw()
p4
FIGURE 4
http://www.spc.noaa.gov/faq/tornado/ef-scale.html
EF scale implemented in the U.S. on 1 February 2007
Values in parentheses are m/s. Operational EF Scale EF0 = 65, 85 (29.06, 38.45) EF1 = 86, 110 (38.45, 49.62) EF2 = 111, 135 (49.62, 60.80) EF3 = 136, 165 (60.80, 74.21) EF4 = 166, 200 (74.21, 89.41) EF5 > 200 (> 89.41)
Derived EF Scale EF0 = 65, 85 (29.06, 38.45) EF1 = 86, 109 (38.45, 49.17) EF2 = 110, 137 (49.17, 61.69) EF3 = 138, 167 (61.69, 75.10) EF4 = 168, 199 (75.10, 89.41) EF5 = 200, 234 (89.41, 104.61)
Fujita Scale F0 = 45, 78 (20.12, 35.32) F1 = 79, 117 (35.32, 52.75) F2 = 118, 161 (52.75, 72.42) F3 = 162, 209 (72.42, 93.88) F4 = 210, 261 (93.88, 117.12) F5 = 262, 317 (117.12, 141.71)
Operational EF scale
windsLo = c(86, 111, 136, 166, 200) * 0.44704
windsHi = c(111, 136, 166, 200, 234) * 0.44704
windsMid = (windsLo + windsHi)/2
diff(windsMid)
## [1] 11.18 12.29 14.31 15.20
mean(diff(windsMid))
## [1] 13.24
Determine wind speeds
model.Final = clm(EF ~ RwLengthS * SWidthS + FAT + YEARS, weights = Weights,
scale = ~SWidthS + RwLengthS, data = Torn3.df, link = "logit")
newData = data.frame(RwLengthS = Torn4.df$RwLengthS, SWidthS = Torn4.df$SWidthS,
FAT = Torn4.df$FAT, YEARS = Torn4.df$YEARS)
probsF = predict(model.Final, newdata = newData, type = "prob", interval = FALSE)
WLo = probsF$fit[, 1:5] %*% windsLo
WHi = probsF$fit[, 1:5] %*% windsHi
WMid = probsF$fit[, 1:5] %*% windsMid
WRan = numeric()
for (j in 1:dim(Torn4.df)[1]) {
windsRan = numeric()
for (i in 1:length(windsLo)) windsRan[i] = runif(1, min = windsLo[i], max = windsHi[i])
WRan[j] = probsF$fit[j, 1:5] %*% windsRan
}
Torn4.df$WHi = WHi
Torn4.df$WLo = WLo
Torn4.df$WMid = WMid
Torn4.df$WRan = WRan
Uncertain.df = data.frame(EF = Torn4.df$EF, WMid, WLo, WHi)
Intensity distribution
Torn4.df$cut = Torn4.df$obs - Torn4.df$pred == 0
Torn4.df$cut2 = cut(Torn4.df$WMid, c(38.45, 49.62, 60.8, 74.21, 89.41))
p5a = ggplot(Torn4.df, aes(x = WMid)) + geom_histogram(binwidth = 1, aes(fill = cut2)) +
facet_grid(FSCALE ~ ., scales = "free_y") + guides(fill = FALSE) + ylab("Frequency") +
xlab(expression(paste("Estimated Tornado Intensity (m ", s^{
-1
}, ")"))) + theme_bw() + geom_vline(xintercept = windsLo, color = "gray") +
geom_histogram(binwidth = 1, aes(fill = cut2)) + scale_x_continuous(limits = c(35,
100))
p5b = ggplot(Torn4.df, aes(x = WRan)) + geom_histogram(binwidth = 1, aes(fill = cut2)) +
facet_grid(FSCALE ~ ., scales = "free_y") + guides(fill = FALSE) + ylab("Frequency") +
xlab(expression(paste("Estimated Tornado Intensity (m ", s^{
-1
}, ")"))) + theme_bw() + geom_vline(xintercept = windsLo, color = "gray") +
geom_histogram(binwidth = 1, aes(fill = cut2)) + scale_x_continuous(limits = c(35,
100))
mat = matrix(c(1, 2), nrow = 2, byrow = TRUE)
p5a = p5a + ggtitle("a") + theme(plot.title = element_text(hjust = 0))
p5b = p5b + ggtitle("b") + theme(plot.title = element_text(hjust = 0))
multiplot(p5a, p5b, layout = mat)
FIGURE 5
model.Final2 = clm(EF ~ RwLength * SWidth + FAT, weights = Weights, scale = ~SWidth +
RwLength, data = Torn3.df, link = "logit")
newData2 = data.frame(RwLength = Torn4.df$RwLength, SWidth = Torn4.df$SWidth,
FAT = Torn4.df$FAT)
probsF2 = predict(model.Final2, newdata = newData2, type = "prob", interval = FALSE)
WBest2 = probsF2$fit[, 1:5] %*% windsBest
Torn4.df$WBest2 = WBest2
IntByState.df = ddply(Torn4.df, .(STATE, YEAR), summarize, avgInt = mean(WBest2),
q95Int = quantile(WBest2, probs = 0.95), nT = length(YEAR))
TrendByState.df = ddply(IntByState.df, .(STATE), summarize, trend = lm(avgInt ~
YEAR)$coef[2], nY = length(STATE))
toofew = TrendByState.df$nY < 8
TrendByState.df$trend[toofew] = NA
TrendByState.df$stN = TrendByState.df$STATE
states.df = map_data("state")
states.df = subset(states.df, region != "alaska" & region != "district of columbia")
states.df$stN = state.abb[match(states.df$region, tolower(state.name))] # attach state abbreviations
df2 = merge(states.df, TrendByState.df, by = "stN")
ggplot(df2, aes(x = long, y = lat, group = group, fill = trend * 10)) + geom_polygon() +
geom_path(color = "gray75") + coord_map(project = "polyconic") + theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(), panel.background = element_blank(),
axis.ticks = element_blank(), axis.text = element_blank(), legend.position = "bottom") +
# labs(title='(E)F1+ Tornadoes\n Annual Intensity Trends (1994-2011)\n
# Data Source: U.S. Storm Prediction Center') +
labs(title = "Intensity Trends (1994-2011)") + xlab("") + ylab("") + scale_fill_gradient2("Trend\n[m/s/decade]",
low = "blue", mid = "white", high = "red", midpoint = 0, space = "rgb",
na.value = "grey50", guide = "colourbar")
model.Final = clm(EF ~ RwLengthS * SWidthS + FAT + YEARS, weights = Weights,
scale = ~SWidthS + RwLengthS, data = Torn3.df, link = "logit")
x = scale(Torn3.df$RwLength)
cRwLength = attr(x, "scaled:center")
sRwLength = attr(x, "scaled:scale")
x = scale(Torn3.df$SWidth)
cSWidth = attr(x, "scaled:center")
sSWidth = attr(x, "scaled:scale")
x = scale(Torn3.df$YEAR)
cYEAR = attr(x, "scaled:center")
sYEAR = attr(x, "scaled:scale")
Wayne, NE October 4, 2013 EF SCALE RATING: EF-4,ESTIMATED PEAK WIND: 170 mph, PATH LENGTH/STATUTE/: 19 PATH WIDTH/MAXIMUM/: 1.38 MILES FATALITIES: 0 INJURIES: 15
170 * 0.44704
## [1] 76
Length = 19 * 1609.34
Width = 1.38 * 1609.34
FAT = FALSE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
pWayne = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT,
YEARS), type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
pWayne[1:5]$lwr %*% windsMid
## [,1]
## 1 74.95
pWayne[1:5]$fit %*% windsMid
## [,1]
## 1 83.3
pWayne[1:5]$upr %*% windsMid
## [,1]
## 1 92.18
Wayne, NE tornado 10/04/13. Path length at the 94th percentile of all tornadoes since 1994 and path width above the 99th percentile.
Wayne, NE tornado 10/04/13. Our stat model predicts a peak wind speed of 83 m/s (75, 92) 95% uncertainty interval.
Newcastle/Moore OK tornado, May 20, 2013 Rating EF-5 Estimated peak wind: Path length: 17 mi Path Width: 1.3 mi
Length = 17 * 1609.34
Width = 1.3 * 1609.34
FAT = TRUE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 86.58
p[1:5]$fit %*% windsMid
## [,1]
## 1 91.63
p[1:5]$upr %*% windsMid
## [,1]
## 1 96.79
Estimated wind speeds from damage survey: 200-210 mph (89-94 m/s)
El Reno, OK
Length = 16.2 * 1609.34
Width = 2.6 * 1609.34
FAT = TRUE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 92.19
p[1:5]$fit %*% windsMid
## [,1]
## 1 95.83
p[1:5]$upr %*% windsMid
## [,1]
## 1 99.41
Bluestein radar estimate: 296 mph (132 m/s), Wurman radar estimate: 246-258 mph (110-115 m/s). The main funnel is believed to have been an EF4, with winds around 185 mph (298 km/h) [83 m/s]. EF5 winds were only found in the smaller vortices that rotated around this funnel at 110 mph (180 km/h). Jeff Masters (June 4, 2013). “Largest Tornado on Record: the May 31 El Reno, OK EF-5 Tornado”. Weather Underground. Retrieved June 4, 2013.
Adairsville GA tornado, Jan 30, 2013 Rating: High-end EF-3 Estimated peak wind: 160 mph (71.5 m/s) Path length: 21.8 miles Path's maximum width: 900 yards Fatalities: One Injuries: 17/9 in Bartow County and eight in Gordon County.
Length = 21.8 * 1609.34
Width = 900 * 0.9144
FAT = TRUE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 77.18
p[1:5]$fit %*% windsMid
## [,1]
## 1 82.21
p[1:5]$upr %*% windsMid
## [,1]
## 1 87.47
# quantile(Torn4.df$Length, probs=seq(.9, .99, .01))
# quantile(Torn4.df$Width, probs=seq(.95, .99, .005))
Our model estimates a wind speed of 85 (81, 89) m/s (95% CI). The 35 km path length puts it above the 95th percentile of all tornadoes since 1985 and the 822 m path width puts it above the 97th percentile of all tornadoes over the same period.
170 * 0.44704
## [1] 76
Length = 21 * 1609.34
Width = 0.75 * 1609.34
Length
## [1] 33796
Width
## [1] 1207
FAT = FALSE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 68.72
p[1:5]$fit %*% windsMid
## [,1]
## 1 74.42
p[1:5]$upr %*% windsMid
## [,1]
## 1 80.5
125 * 0.44704
## [1] 55.88
Length = 4.3 * 1609.34
Width = 220 * 0.9144
Length
## [1] 6920
Width
## [1] 201.2
FAT = FALSE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 50.14
p[1:5]$fit %*% windsMid
## [,1]
## 1 53.38
p[1:5]$upr %*% windsMid
## [,1]
## 1 56.78
145 * 0.44704
## [1] 64.82
Length = 68 * 1609.34
Width = 0.75 * 1609.34
Length
## [1] 109435
Width
## [1] 1207
FAT = TRUE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 83.14
p[1:5]$fit %*% windsMid
## [,1]
## 1 88.96
p[1:5]$upr %*% windsMid
## [,1]
## 1 94.94
180 * 0.44704
## [1] 80.47
Length = 2.75 * 1609.34
Width = 800 * 0.9144
Length
## [1] 4426
Width
## [1] 731.5
FAT = TRUE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 69.05
p[1:5]$fit %*% windsMid
## [,1]
## 1 75.69
p[1:5]$upr %*% windsMid
## [,1]
## 1 82.86
165 * 0.44704
## [1] 73.76
185 * 0.44704
## [1] 82.7
Length = 7 * 1609.34
Width = 1100 * 0.9144
Length
## [1] 11265
Width
## [1] 1006
FAT = FALSE
YEAR = 2013
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 63.97
p[1:5]$fit %*% windsMid
## [,1]
## 1 69.7
p[1:5]$upr %*% windsMid
## [,1]
## 1 75.9
180 * 0.44704
## [1] 80.47
Length = 26.5 * 1609.34
Width = 275 * 0.9144
Length
## [1] 42648
Width
## [1] 251.5
FAT = TRUE
YEAR = 2012
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 67.41
p[1:5]$fit %*% windsMid
## [,1]
## 1 74.26
p[1:5]$upr %*% windsMid
## [,1]
## 1 81.69
140 * 0.44704
## [1] 62.59
Length = 86 * 1609.34
Width = 1 * 1609.34
Length
## [1] 138403
Width
## [1] 1609
FAT = TRUE
YEAR = 2012
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 85.36
p[1:5]$fit %*% windsMid
## [,1]
## 1 91.34
p[1:5]$upr %*% windsMid
## [,1]
## 1 97.44
130 * 0.44704
## [1] 58.12
Length = 13.7 * 1609.34
Width = 200 * 0.9144
Length
## [1] 22048
Width
## [1] 182.9
FAT = FALSE
YEAR = 2012
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 53.02
p[1:5]$fit %*% windsMid
## [,1]
## 1 56.95
p[1:5]$upr %*% windsMid
## [,1]
## 1 61.12
111 * 0.44704
## [1] 49.62
135 * 0.44704
## [1] 60.35
Length = 5.7 * 1609.34
Width = 200 * 0.9144
Length
## [1] 9173
Width
## [1] 182.9
FAT = FALSE
YEAR = 2012
RwLength = Length^(1/3)
SWidth = sqrt(Width)
YEARS = (YEAR - cYEAR)/sYEAR
RwLengthS = (RwLength - cRwLength)/sRwLength
SWidthS = (SWidth - cSWidth)/sSWidth
p = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
p[1:5]$lwr %*% windsMid
## [,1]
## 1 50.69
p[1:5]$fit %*% windsMid
## [,1]
## 1 53.88
p[1:5]$upr %*% windsMid
## [,1]
## 1 57.22
Summarize
model = c(74, 91, 57, 54, 82, 74, 53, 89, 76, 70, 91, 96, 83)
survey = c(80, 63, 57, 55, 72, 76, 56, 64, 80, 78, 91, 112, 76)
cor.test(model, survey)
##
## Pearson's product-moment correlation
##
## data: model and survey
## t = 2.953, df = 11, p-value = 0.01315
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1798 0.8899
## sample estimates:
## cor
## 0.6649
sqrt(mean((model - survey)^2))
## [1] 12.23
cor.test(model[-c(2, 8)], survey[-c(2, 8)])
##
## Pearson's product-moment correlation
##
## data: model[-c(2, 8)] and survey[-c(2, 8)]
## t = 6.55, df = 9, p-value = 0.0001051
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6804 0.9765
## sample estimates:
## cor
## 0.9092
sqrt(mean((model[-c(2, 8)] - survey[-c(2, 8)])^2))
## [1] 6.974
# May 16, 1995 Hanston, KS
Ism = numeric()
xx = Torn3.df[Torn3.df$DATE == "1995/05/15", ]
RwLengthS = xx[5, ]$RwLengthS
SWidthS = xx[5, ]$SWidthS
FAT = xx[5, ]$FAT
YEARS = xx[5, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 25, 1996 Friona, TX
xx = Torn3.df[Torn3.df$DATE == "1996/05/25", ]
RwLengthS = xx[3, ]$RwLengthS
SWidthS = xx[3, ]$SWidthS
FAT = xx[3, ]$FAT
YEARS = xx[3, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# April 10, 1997 Tulia, TX
xx = Torn3.df[Torn3.df$DATE == "1997/04/10", ]
RwLengthS = xx[4, ]$RwLengthS
SWidthS = xx[4, ]$SWidthS
FAT = xx[4, ]$FAT
YEARS = xx[4, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 26, 1997 Kiefer, OK
xx = Torn3.df[Torn3.df$DATE == "1997/05/26", ]
RwLengthS = xx[2, ]$RwLengthS
SWidthS = xx[2, ]$SWidthS
FAT = xx[2, ]$FAT
YEARS = xx[2, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 30, 1998 Spencer, SD
xx = Torn3.df[Torn3.df$DATE == "1998/05/30", ]
RwLengthS = xx[4, ]$RwLengthS
SWidthS = xx[4, ]$SWidthS
FAT = xx[4, ]$FAT
YEARS = xx[4, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 1, 1999 Tarzan, TX
xx = Torn3.df[Torn3.df$DATE == "1999/05/01", ]
RwLengthS = xx[1, ]$RwLengthS
SWidthS = xx[1, ]$SWidthS
FAT = xx[1, ]$FAT
YEARS = xx[1, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 3, 1999 Moore, OK
xx = Torn3.df[Torn3.df$DATE == "1999/05/03", ]
RwLengthS = xx[4, ]$RwLengthS
SWidthS = xx[4, ]$SWidthS
FAT = xx[4, ]$FAT
YEARS = xx[4, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 31, 1999 Sitka, KS
xx = Torn3.df[Torn3.df$DATE == "1999/05/31", ]
RwLengthS = xx[3, ]$RwLengthS
SWidthS = xx[3, ]$SWidthS
FAT = xx[3, ]$FAT
YEARS = xx[3, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# June 4 1999 Thedford, NE
xx = Torn3.df[Torn3.df$DATE == "1999/06/04", ]
RwLengthS = xx[6, ]$RwLengthS
SWidthS = xx[6, ]$SWidthS
FAT = xx[6, ]$FAT
YEARS = xx[6, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 24 2004 Hebron, NE
xx = Torn3.df[Torn3.df$DATE == "2004/05/24", ]
RwLengthS = xx[3, ]$RwLengthS
SWidthS = xx[3, ]$SWidthS
FAT = xx[3, ]$FAT
YEARS = xx[3, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# May 10 2008 Stuttgart, AR
xx = Torn3.df[Torn3.df$DATE == "2008/05/10", ]
RwLengthS = xx[13, ]$RwLengthS
SWidthS = xx[13, ]$SWidthS
FAT = xx[13, ]$FAT
YEARS = xx[13, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
# June 5 2009 Goshen County, WY
xx = Torn3.df[Torn3.df$DATE == "2009/06/05", ]
RwLengthS = xx[1, ]$RwLengthS
SWidthS = xx[1, ]$SWidthS
FAT = xx[1, ]$FAT
YEARS = xx[1, ]$YEARS
pp = predict(model.Final, newdata = data.frame(RwLengthS, SWidthS, FAT, YEARS),
type = "prob", se.fit = TRUE, interval = TRUE, level = 0.95)
Ism = c(Ism, pp[1:5]$fit %*% windsMid)
Comparisons
Idow = c(82.6, 114.2, 54.1, 90.2, 205.5, 65.8, 117, 65, 96.4, 80.5, 113.4, 108.6)
Iwsr = c(65, 64, 46.5, 48.6, 111.9, 42.5, 92.9, 48, 74.8, 62.6, 81, 53)
Mdow = c(62.88, 64.38, 53.64, 37.69, 98.7, 38.97, 83.46, 57.88, 84.3, 59.44,
91.25, 52.67)
cor.test(Ism, Idow)
##
## Pearson's product-moment correlation
##
## data: Ism and Idow
## t = 3.184, df = 10, p-value = 0.009747
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2289 0.9121
## sample estimates:
## cor
## 0.7096
cor.test(Ism, Iwsr)
##
## Pearson's product-moment correlation
##
## data: Ism and Iwsr
## t = 7.498, df = 10, p-value = 2.067e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7375 0.9781
## sample estimates:
## cor
## 0.9214
cor.test(Ism, Mdow)
##
## Pearson's product-moment correlation
##
## data: Ism and Mdow
## t = 4.306, df = 10, p-value = 0.001546
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4319 0.9435
## sample estimates:
## cor
## 0.806
sqrt(mean((Ism - Mdow)^2))
## [1] 11.41