Seasonal Prediction Model
James B. Elsner and Thomas H. Jagger
Regular grids across multiple states. The framework uses spatial rasters.
Get the tornado data. Set the domain and grid resolution.
if(!file.exists("torn")){
download.file(url = "http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
destfile = "tornado.zip")
unzip("tornado.zip")
}
library("rgdal")## Loading required package: sp## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3TornL = readOGR(dsn = "torn", layer = "torn",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "torn", layer: "torn"
## with 60114 features
## It has 22 fieldslibrary("raster")
xmn = -104; xmx = -80; ymn = 25; ymx = 45
r = raster(xmn = xmn, xmx = xmx,
ymn = ymn, ymx = ymx,
resolution = 2)
dim(r)## [1] 10 12 1nCells = ncell(r)
sp = as(r, 'SpatialPolygons')
spT = spTransform(sp, CRS(proj4string(TornL)))Create a domain area map.
library("ggmap")## Loading required package: ggplot2Map = get_map(location = c(mean(c(xmn, xmx)), mean(c(ymn, ymx))),
source = "google",
zoom = 4,
color = "bw",
maptype = "terrain")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=35,-92&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=falsep1 = ggmap(Map, dev = "extent") +
geom_segment(aes(x = xmn, xend = xmx, y = ymn, yend = ymn),
color = "red") +
geom_segment(aes(x = xmn, xend = xmx, y = ymx, yend = ymx),
color = "red") +
geom_segment(aes(x = xmn, xend = xmn, y = ymn, yend = ymx),
color = "red") +
geom_segment(aes(x = xmx, xend = xmx, y = ymn, yend = ymx),
color = "red") +
labs(x = expression(paste("Longitude (", degree, "W)")),
y = expression(paste("Latitude (", degree, "N)")))
library("maps")##
## # maps v3.1: updated 'world': all lakes moved to separate new #
## # 'lakes' database. Type '?world' or 'news(package="maps")'. #library("maptools")## Checking rgeos availability: TRUElibrary("grid")
p1 + theme(panel.grid.minor = element_line(colour = NA),
panel.grid.minor = element_line(colour = NA),
panel.background = element_rect(fill = NA, colour = NA),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.title = element_blank(),
rect = element_blank())
Figure: Prediction domain
Subset the tornado data. Overlay tornado paths on the raster. Each element of the list is a cell subset of the tornado track file. The data frame nTor.df contains the per-cell tornado count, area, and rate. The PathsAll.df data frame contains the tornado attributes listed by cell. Information is repeated for each cell in cases where the path intersects more than one cell. Area is square meter. The unadjusted climatological rate (climoU) has units of tornadoes per 100 x 100 km square per year.
library("rgeos")## rgeos version: 0.3-19, (SVN revision 524)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-3
## Polygon checking: TRUEbeginYear = 1954
endYear = 2015
TornL2 = subset(TornL, yr >= beginYear & mag >= 0)
TornL2ll = spTransform(TornL2, proj4string(r))
TnT = dim(TornL2)[1]
TornL2$Sid = 1:TnT
ct = over(spT, TornL2, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
Nyears = endYear - beginYear + 1
area = rgeos::gArea(spT, byid = TRUE)
nTor.df = data.frame(ID = 1:nCells,
Nyears = Nyears,
nT = nT,
area = area,
climoU = nT/area * 10^6 * 100 * 100 / Nyears
)
PathsAll.df = plyr::ldply(ct, data.frame)
ID = rep(nTor.df$ID, nT)
PathsAll.df$ID = ID
TnTdom = length(unique(PathsAll.df$Sid))Percentage of all tornadoes within the domain by EF rating.
TnTdom/TnT * 100## [1] 81.70599length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 1]))/dim(TornL2[TornL2$mag >= 1, ])[1] * 100## [1] 83.90558length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 2]))/dim(TornL2[TornL2$mag >= 2, ])[1] * 100## [1] 87.53938length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 3]))/dim(TornL2[TornL2$mag >= 3, ])[1] * 100## [1] 89.89864length(unique(PathsAll.df$Sid[PathsAll.df$mag >= 4]))/dim(TornL2[TornL2$mag >= 4, ])[1] * 100## [1] 91.60839Get annual counts by domain. The data frame PathsAll.df contains the list of tornadoes by grid cell.
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionnTor.year = PathsAll.df %>%
group_by(yr) %>%
summarize(nT = length(unique(Sid)),
nTd = length(unique(date)))
nTor.year = as.data.frame(nTor.year)
nTor.year$yr2 = nTor.year$yr
sum(nTor.year$nT)## [1] 48200library("ggplot2")
ggplot(nTor.year, aes(x = yr, y = nT)) +
geom_line() +
geom_text(aes(label = nT), vjust = -.5, size = 4, color = "darkblue") +
xlab("Year") + ylab("Annual Number of EF0+ Tornadoes") +
scale_y_continuous(limits = c(0, 2000)) +
# scale_y_continuous(limits = c(0, 500)) +
theme(axis.text.x = element_text(size = 9),
legend.position = "none")
Figure: Time series tornado counts over the domain
Create a space-time data frame with 0s in years/cells without tornadoes.
dat = PathsAll.df %>%
# filter(mo > 10) %>%
# filter(mag >= 1) %>%
group_by(ID, yr) %>%
summarize(nT = n(),
nTd = length(unique(date)))
zz = expand.grid(ID = 1:nCells, yr = beginYear:endYear)
all = left_join(zz, dat)## Joining, by = c("ID", "yr")all0 = all
all0[is.na(all)] = 0Create a space-time raster stack.
STrs = stack()
for(i in beginYear:endYear){
df = all0[all0$yr == i, ]
ra = r
values(ra) = df$nT
STrs = stack(STrs, ra)
}
names(STrs) = beginYear:endYearMap boundaries
library("mapproj")
library("maptools")
ext = as.vector(extent(r))
bndryC = map("county", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(bndryC$names, ":"), function(x) x[1])
bndryCP = map2SpatialPolygons(bndryC, IDs = IDs,
proj4string = CRS(projection(r)))
bndryS = map("state", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(bndryS$names, ":"), function(x) x[1])
bndrySP = map2SpatialPolygons(bndryS, IDs = IDs,
proj4string = CRS(projection(r)))
bndryUSA = map("usa", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(bndryUSA$names, ":"), function(x) x[1])
bndryUSAP = map2SpatialPolygons(bndryUSA, IDs = IDs,
proj4string = CRS(projection(r)))Put NAs in grid cells that lie outside U.S. land borders.
nTr = STrs
values(nTr) = all0$nT
nTrm = mask(nTr, bndryUSAP)
all0$nT = as.vector(values(nTrm))Plots by year. Only last 16 years.
s = names(STrs)
s = s[(length(s) - 16 + 1):length(s)]
STrsP = STrs[[s]]
library('rasterVis')## Loading required package: lattice## Loading required package: latticeExtra## Loading required package: RColorBrewer##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layerlibrary('RColorBrewer')
library("wesanderson")
range(log2(values(STrsP)), na.rm = TRUE)## [1] -Inf 6.189825rng = seq(0, 7, 1)
crq = wes_palette(7, name = "Zissou", type = "continuous")
labs = as.character(round(2^rng, 1))levelplot(log2(STrsP), margin = FALSE,
sub = " Number of Tornadoes",
xlab = NULL, ylab = NULL,
names.attr = 2000:endYear,
layout = c(4, 4),
col.regions = crq, at = rng,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP,
col = gray(.35), lwd = 1)) Figure: Space-time plot
Merge area and climotology to the all0 data frame.
df = data.frame(ID = 1:nCells, area, climo = nTor.df$climo)
allp1 = merge(all0, df, by = "ID")Add climate covariates. http://www.esrl.noaa.gov/psd/data/correlation/censo.data http://www1.ncdc.noaa.gov/pub/data/cmb/ersst/v3b/pdo/ Things tried: Global angular momentum, solar flux, QBO, EA, WP, TPI.
ENSO = read.table("ENSO.txt", header = TRUE) #Bivariate ENSO (BEST) + -> El Nino Mar-Jul + -> El Nino
ENSO = subset(ENSO, Year >= beginYear & Year <= endYear)
NAO = read.table("NAO.txt", header = TRUE)
NAO = subset(NAO, Year >= beginYear & Year <= endYear)
CI = read.table("PNA.txt", header = TRUE) # Mar PNA +
CI = subset(CI, Year >= beginYear & Year <= endYear)
WCA = read.table("WCA.txt", header = TRUE)
GAK = read.table("GAK.txt", header = TRUE)
WCA = subset(WCA, Year >= beginYear & Year <= endYear)
GAK = subset(GAK, Year >= beginYear & Year <= endYear)
Cov.df0 = data.frame(
enso = (ENSO$Mar + ENSO$Apr + ENSO$May)/3 ,
nao = (NAO$Apr + NAO$May + NAO$Jun)/3,
wca = WCA$Feb,
gak = GAK$Apr,
td = WCA$Feb - GAK$Apr)
scaled = scale(Cov.df0)
scaling = data.frame(t(data.frame(attributes(scaled)[c(3:4)])))
Cov.df = cbind(yr = ENSO$Year, scaled)
allp3 = merge(allp1, Cov.df, by = "yr", all = TRUE)
allp3 = arrange(allp3, yr, ID)
allp3$yr2 = allp3$yr
allp3$ID2 = allp3$ID
allp3$ID3 = allp3$ID
allp3$ID4 = allp3$IDSome global controls for INLA. Use as needed.
library("INLA")## Loading required package: Matrix## Loading required package: splines## This is INLA 0.0-1468872408, dated 2016-07-18 (14:43:05+0100).
## See www.r-inla.org/contact-us for how to get help.control = list(
predictor = list(compute = TRUE),
inla = list(strategy = "laplace",
fast = FALSE,
stencil = 7,
npoints = 198,
int.strategy = "grid",
dz = .5),
results = list(return.marginals.random = TRUE),
compute = list(config = TRUE, mlik = TRUE, cpo = TRUE, dic = TRUE, po = TRUE),
family = list(variant = 1, hyper = list(theta = list(prior = "loggamma",
param = c(1, 1)))))Spatial neighborhood definition as an inla graph.
library('spdep')
nb = poly2nb(sp)
nb2INLA("g", nb)
g = inla.read.graph("g")The area times the Nyears indicates the square meter years exposed to tornadoes (E). Scale the exposure to have a mean of unity. The model has a spatial id.
http://www.ospo.noaa.gov/Products/ocean/sst/anomaly/index.html
The model.
allp3$E = allp3$area/10^10
allp3$ID1 = factor(allp3$ID)
predcov = subset(allp3, yr == 1991)
predcov$nT = NA
#values for prediction
predcov$yr = mean(allp3$yr) #No population bias inflation
# 2016
# predcov$nao = 0
# predcov$enso = 0
# predcov$wca = 1.5
# predcov$gak = .25
# 2011
predcov$nao = -.5
predcov$enso = -2
predcov$wca = 0
predcov$gak = -1.5
predcov$td = with(predcov, (wca * scaling$wca[2] - gak * scaling$gak[2])/scaling$td[2])
# values for the null model... using means, gives mean observed rates...
allp3p = rbind(predcov, allp3)
#Use allp3p in the model for prediction...
formula = nT ~ -1 + ID1 + I(yr - 1984) +
f(ID2, enso, model = "besag", graph = g, constr = FALSE) +
f(ID3, nao, model = "besag", graph = g, constr = FALSE) +
f(ID4, td, model = "besag", graph = g, constr = FALSE)
model = inla(formula = formula, family = "nbinomial",
E = E ,
data = allp3p,
quantiles = c(.05, .5, .95),
control.compute = control$compute,
control.fixed = list(expand.factor.strategy = "inla"),
control.predictor = list(link = 1),
control.results = control$results,
control.family = control$family
)
summary(model)##
## Call:
## c("inla(formula = formula, family = \"nbinomial\", data = allp3p, ", " quantiles = c(0.05, 0.5, 0.95), E = E, control.compute = control$compute, ", " control.predictor = list(link = 1), control.family = control$family, ", " control.results = control$results, control.fixed = list(expand.factor.strategy = \"inla\"))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 3.8345 19.5325 1.7642 25.1312
##
## Fixed effects:
## mean sd 0.05quant 0.5quant 0.95quant mode kld
## ID11 -0.2167 0.1184 -0.4111 -0.2172 -0.0209 -0.2181 0
## ID12 -0.0649 0.1145 -0.2525 -0.0655 0.1247 -0.0668 0
## ID13 0.4833 0.1051 0.3119 0.4823 0.6582 0.4803 0
## ID14 0.9719 0.0997 0.8098 0.9707 1.1382 0.9683 0
## ID15 0.7678 0.1014 0.6027 0.7667 0.9368 0.7644 0
## ID16 0.8380 0.1010 0.6736 0.8368 1.0063 0.8345 0
## ID17 0.2329 0.1093 0.0543 0.2321 0.4144 0.2304 0
## ID18 0.8146 0.1012 0.6498 0.8135 0.9833 0.8111 0
## ID19 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID110 0.0007 0.1131 -0.1844 0.0001 0.1882 -0.0013 0
## ID111 -0.4036 0.1221 -0.6043 -0.4039 -0.2019 -0.4045 0
## ID112 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID113 0.5920 0.1021 0.4256 0.5910 0.7620 0.5889 0
## ID114 0.2591 0.1064 0.0853 0.2583 0.4360 0.2565 0
## ID115 0.7123 0.1009 0.5478 0.7112 0.8804 0.7090 0
## ID116 1.0285 0.0980 0.8692 1.0273 1.1920 1.0250 0
## ID117 0.8797 0.0989 0.7188 0.8785 1.0445 0.8762 0
## ID118 1.0780 0.0972 0.9199 1.0768 1.2402 1.0743 0
## ID119 0.9659 0.0982 0.8061 0.9647 1.1298 0.9624 0
## ID120 0.8637 0.0993 0.7021 0.8625 1.0291 0.8602 0
## ID121 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID122 0.7383 0.1010 0.5738 0.7372 0.9065 0.7350 0
## ID123 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID124 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID125 1.0157 0.0963 0.8590 1.0145 1.1763 1.0122 0
## ID126 0.8713 0.0977 0.7123 0.8702 1.0342 0.8680 0
## ID127 1.1621 0.0955 1.0068 1.1609 1.3214 1.1585 0
## ID128 1.1012 0.0959 0.9452 1.1000 1.2612 1.0976 0
## ID129 1.0377 0.0967 0.8804 1.0365 1.1990 1.0342 0
## ID130 0.4113 0.1034 0.2424 0.4104 0.5832 0.4086 0
## ID131 0.6931 0.0995 0.5310 0.6920 0.8588 0.6899 0
## ID132 1.1812 0.0952 1.0265 1.1799 1.3400 1.1775 0
## ID133 0.9450 0.0974 0.7864 0.9439 1.1074 0.9416 0
## ID134 0.8466 0.0985 0.6862 0.8455 1.0107 0.8433 0
## ID135 0.4196 0.1033 0.2511 0.4187 0.5913 0.4169 0
## ID136 -0.1652 0.1137 -0.3517 -0.1657 0.0230 -0.1666 0
## ID137 0.4961 0.1000 0.3330 0.4952 0.6625 0.4933 0
## ID138 0.9177 0.0958 0.7618 0.9166 1.0774 0.9144 0
## ID139 1.2880 0.0933 1.1364 1.2868 1.4437 1.2843 0
## ID140 1.1193 0.0946 0.9655 1.1182 1.2771 1.1159 0
## ID141 0.9704 0.0959 0.8143 0.9693 1.1303 0.9671 0
## ID142 0.6661 0.0984 0.5058 0.6651 0.8299 0.6631 0
## ID143 0.6442 0.0987 0.4832 0.6432 0.8085 0.6412 0
## ID144 0.8488 0.0964 0.6917 0.8477 1.0095 0.8455 0
## ID145 0.5788 0.0995 0.4166 0.5778 0.7443 0.5759 0
## ID146 0.5571 0.0994 0.3950 0.5562 0.7226 0.5542 0
## ID147 -0.4117 0.1171 -0.6041 -0.4120 -0.2184 -0.4125 0
## ID148 -1.5493 0.1648 -1.8242 -1.5472 -1.2819 -1.5427 0
## ID149 0.0050 0.1072 -0.1706 0.0044 0.1827 0.0032 0
## ID150 1.1052 0.0932 0.9536 1.1040 1.2606 1.1018 0
## ID151 1.0955 0.0936 0.9432 1.0944 1.2517 1.0922 0
## ID152 1.4723 0.0912 1.3242 1.4711 1.6246 1.4687 0
## ID153 1.2101 0.0925 1.0596 1.2089 1.3645 1.2066 0
## ID154 0.6532 0.0973 0.4945 0.6522 0.8151 0.6502 0
## ID155 0.7476 0.0962 0.5909 0.7466 0.9078 0.7446 0
## ID156 0.7335 0.0962 0.5767 0.7325 0.8937 0.7304 0
## ID157 0.6369 0.0968 0.4791 0.6359 0.7981 0.6340 0
## ID158 0.2156 0.1025 0.0482 0.2148 0.3857 0.2133 0
## ID159 -0.5260 0.1183 -0.7206 -0.5261 -0.3309 -0.5264 0
## ID160 0.1213 0.1043 -0.0494 0.1206 0.2944 0.1192 0
## ID161 0.8275 0.0948 0.6729 0.8265 0.9854 0.8245 0
## ID162 1.1067 0.0920 0.9570 1.1056 1.2601 1.1034 0
## ID163 1.0745 0.0924 0.9241 1.0735 1.2286 1.0713 0
## ID164 1.2181 0.0915 1.0693 1.2170 1.3706 1.2147 0
## ID165 0.7400 0.0950 0.5853 0.7391 0.8981 0.7371 0
## ID166 0.6651 0.0955 0.5094 0.6641 0.8240 0.6622 0
## ID167 0.7009 0.0952 0.5457 0.7000 0.8595 0.6980 0
## ID168 0.6654 0.0956 0.5096 0.6645 0.8246 0.6626 0
## ID169 1.0625 0.0919 0.9129 1.0614 1.2158 1.0593 0
## ID170 0.6784 0.0954 0.5228 0.6774 0.8373 0.6755 0
## ID171 0.3912 0.0988 0.2299 0.3903 0.5553 0.3886 0
## ID172 0.6320 0.0958 0.4758 0.6311 0.7914 0.6292 0
## ID173 0.3105 0.0993 0.1481 0.3098 0.4754 0.3083 0
## ID174 0.6002 0.0950 0.4453 0.5993 0.7582 0.5975 0
## ID175 0.7968 0.0931 0.6451 0.7959 0.9519 0.7939 0
## ID176 1.1431 0.0905 0.9959 1.1421 1.2940 1.1399 0
## ID177 0.8965 0.0924 0.7461 0.8955 1.0504 0.8935 0
## ID178 0.7369 0.0934 0.5846 0.7359 0.8924 0.7340 0
## ID179 0.9565 0.0915 0.8074 0.9555 1.1090 0.9534 0
## ID180 0.9486 0.0914 0.7997 0.9475 1.1009 0.9455 0
## ID181 0.5339 0.0956 0.3780 0.5330 0.6928 0.5313 0
## ID182 0.6691 0.0943 0.5153 0.6682 0.8260 0.6663 0
## ID183 0.4933 0.0966 0.3356 0.4924 0.6538 0.4907 0
## ID184 -0.2153 0.1079 -0.3923 -0.2157 -0.0368 -0.2165 0
## ID185 -1.0775 0.1323 -1.2965 -1.0769 -0.8609 -1.0755 0
## ID186 -1.1341 0.1339 -1.3558 -1.1333 -0.9150 -1.1317 0
## ID187 0.0097 0.1021 -0.1574 0.0092 0.1788 0.0081 0
## ID188 0.7425 0.0927 0.5914 0.7416 0.8967 0.7398 0
## ID189 1.0485 0.0901 0.9018 1.0475 1.1986 1.0454 0
## ID190 0.7194 0.0923 0.5689 0.7185 0.8731 0.7167 0
## ID191 0.5778 0.0937 0.4249 0.5769 0.7336 0.5752 0
## ID192 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID193 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID194 0.2197 0.0985 0.0587 0.2190 0.3831 0.2176 0
## ID195 0.1543 0.0996 -0.0087 0.1537 0.3196 0.1524 0
## ID196 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID197 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID198 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID199 -0.6121 0.1154 -0.8022 -0.6121 -0.4220 -0.6120 0
## ID1100 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1101 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1102 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1103 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1104 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1105 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1106 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1107 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1108 0.8692 0.0898 0.7228 0.8682 1.0186 0.8664 0
## ID1109 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1110 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1111 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1112 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1113 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1114 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1115 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1116 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1117 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1118 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1119 0.0000 31.6228 -52.0820 -0.0009 52.0918 0.0000 0
## ID1120 0.7480 0.0893 0.6024 0.7472 0.8966 0.7454 0
## I(yr - 1984) 0.0123 0.0006 0.0112 0.0123 0.0133 0.0123 0
##
## Random effects:
## Name Model
## ID2 Besags ICAR model
## ID3 Besags ICAR model
## ID4 Besags ICAR model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 5.565e-01 1.430e-02
## Precision for ID2 5.479e+01 2.567e+01
## Precision for ID3 3.451e+02 4.927e+02
## Precision for ID4 1.891e+04 1.865e+04
## 0.05quant 0.5quant
## size for the nbinomial observations (overdispersion) 0.5329 5.565e-01
## Precision for ID2 26.4047 4.860e+01
## Precision for ID3 66.4163 2.010e+02
## Precision for ID4 1653.2694 1.332e+04
## 0.95quant mode
## size for the nbinomial observations (overdispersion) 5.799e-01 0.5568
## Precision for ID2 1.039e+02 38.9929
## Precision for ID3 1.071e+03 101.4554
## Precision for ID4 5.468e+04 2750.6498
##
## Expected number of effective parameters(std dev): 132.48(9.213)
## Number of equivalent replicates : 42.59
##
## Deviance Information Criterion (DIC) ...: 33578.40
## Effective number of parameters .........: 135.65
##
## Watanabe-Akaike information criterion (WAIC) ...: 33594.75
## Effective number of parameters .................: 145.27
##
## Marginal log-Likelihood: -17566.07
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedPredictions relative rates:
nID = nrow(predcov)
climatology = with(model$summary.fixed, exp(mean + .5 * sd^2)[1:nID])
rates = model$summary.fitted$mean[1:nID]/climatology
Rp = r
covnam = c("enso", "nao", "td")
covval = predcov[1, covnam]
values(Rp) = rates * 100
Rpm = mask(Rp, bndrySP)
range(values(Rpm), na.rm = TRUE)## [1] 79.80353 149.63269#rngp = seq(80, 120, 5) # 2015
#rngp = seq(96, 104, 1) # 2016
rngp = seq(50, 150, 10) # 2011
cr = brewer.pal(10, "BrBG")
cr = cr[-(1:3)]
rngp = rngp[-(1:3)]
#cr = cr[-c(1, 8)] # 2015
#rngp = rngp[-c(1, 9)] # 2015
labs = paste(as.character(rngp), "%", sep ="")
levelplot(Rpm, margin = FALSE,
# sub = paste("2011 Forecast (% of normal)\n",
# paste(covnam,format(covval, digits = 2), collapse = " ", sep = " = "), sep = " "),
sub = paste(" Hindcast of Tornado Activity for 2011 [% of Long-Term Rate]"),
xlab = NULL, ylab = NULL,
col.regions = cr, at = rngp,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP,
col = gray(.15), lwd = 1)) +
latticeExtra::layer(sp.lines(TornL2ll[TornL2ll$yr == 2011, ], col = "white"))
as.data.frame(TornL2) %>%
filter(st == "KS") %>%
summarize(nY = length(unique(yr)),
avg = length(st)/nY)## nY avg
## 1 62 62.45161