Code to reproduce Figure8.4
grid <- unique(euso3[, 1:3])
head(grid)
## s.index Longitude Latitude
## 1 1 -85.802 33.281
## 1531 2 -86.137 32.498
## 3061 3 -86.915 33.486
## 4591 4 -87.004 33.331
## 6121 5 -86.817 33.386
## 7651 6 -86.669 33.705
useastmap <- map_data(database="state",regions=eaststates)
# head(useastmap)
head(knots.coords.151)
## lon lat
## 23 -96.59 48.225
## 24 -96.59 46.950
## 25 -96.59 45.675
## 44 -94.98 48.225
## 45 -94.98 46.950
## 46 -94.98 45.675
pknots <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_point(data=knots.coords.151, aes(x=lon, y =lat), inherit.aes = FALSE, col="Red", size=1.5) +
labs( title= "151 knot sites in the Eastern US", x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank()) +
ggsn::north(data=useastmap, location="bottomright", symbol=12) +
ggsn::scalebar(data =useastmap, dist =500, location ="bottomleft", transform=T, dist_unit = "km",
st.dist = .05, st.size = 5, height = .06, st.bottom=T, model="WGS84")
pknots
ggsave(file=paste0(figurepath, "usagrid_knots.png"))
## Saving 7 x 5 in image
euso3$xsqcmax <- scale(euso3$sqcmax)
mean(euso3$xsqcmax)
sd(euso3$xsqcmax)
time.data <- spT.time(t.series =153, segments=10 )
f2 <- o8hrmax ~ xsqcmax
M2 <- Bsptime(formula=f2, data=euso3, package="spTimer", model="GPP",
coordtype = "lonlat", coords=2:3, time.data = time.data,
prior.phi="Gamm", scale.transform = "SQRT", # validrows=vs,
knots.coords = knots.coords.151, n.report=20, N=3000, burn.in = 1000)
plot(M2$fit)
summary(M2)
table8.2 <- M2$params
dput(table8.2, file=paste0(filepath, "table8.2.txt"))
#outputpath <- "~/Dropbox/sks/bookproject/rbook/mcmc_output_files/uspollution/"
# save(M2, euso3, f2, time.data, knots.coords.151, file=paste0(outputpath, "a1.RData"))
load(file=paste0(outputpath, "a1.RData"))
table8.2 <- M2$params
dput(table8.2, file=paste0(filepath, "table8.2.txt"))
summary(M2)
##
## The GPP model has been fitted using the spTimer package.
## Call:
## Bsptime(formula = f1, data = euso3, package = "spTimer", model = "GPP",
## coordtype = "lonlat", coords = 2:3, validrows = vs, scale.transform = "SQRT",
## prior.phi = "Gamm", time.data = time.data, knots.coords = knots.coords.151,
## N = 3000, burn.in = 1000, n.report = 20)
## 15 - Hour/s. 15 - Mins. 59.49 - Sec.
## Model formula
## o8hrmax ~ xsqcmax
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) 6.737 0.025 6.709 6.811
## xsqcmax 0.311 0.004 0.304 0.318
## rho 0.355 0.003 0.350 0.361
## sig2eps 0.264 0.000 0.263 0.265
## sig2eta 0.747 0.004 0.739 0.756
## phi 0.003 0.000 0.003 0.003
##
## Validation Statistics:
## rmse mae crps cvg
## 10.454 7.693 15.733 99.555
plot(M2$fit)
modfit <- M2$fit
names(modfit)
## [1] "phip" "accept" "nup"
## [4] "sig2eps" "sig2etap" "betap"
## [7] "rhop" "mu_lp" "sig2lp"
## [10] "w0p" "wp" "fitted"
## [13] "X" "Y" "call"
## [16] "tol.dist" "distance.method" "cov.fnc"
## [19] "scale.transform" "sampling.sp.decay" "covariate.names"
## [22] "Distance.matrix.knots" "knots.coords" "coords"
## [25] "KM" "n" "r"
## [28] "T" "p" "knots"
## [31] "initials" "priors" "PMCC"
## [34] "iterations" "nBurn" "computation.time"
## [37] "combined.fit.pred" "model" "parameter"
## [40] "data"
dim(modfit$fitted)
## [1] 1057230 2
head(modfit$fitted)
## Mean SD
## [1,] 7.227619 0.5446847
## [2,] 7.126412 0.5475193
## [3,] 6.919424 0.5496207
## [4,] 7.528879 0.5401945
## [5,] 7.887672 0.5525673
## [6,] 7.928116 0.5330198
summary(modfit$fitted)
## Mean SD
## Min. : 2.723 Min. :0.4925
## 1st Qu.: 6.303 1st Qu.:1.1923
## Median : 7.029 Median :1.7194
## Mean : 7.011 Mean :2.0661
## 3rd Qu.: 7.731 3rd Qu.:2.5893
## Max. :11.491 Max. :7.3601
dim(modfit$w0p)
## [1] 1510 2000
dim(modfit$wp)
## [1] 231030 2000
dim(euso3)
## [1] 1057230 11
dim(modfit$betap)
## [1] 2 2000
dim(modfit$wp)
## [1] 231030 2000
dim(modfit$X)
## [1] 1057230 2
head(modfit$X)
## (Intercept) xsqcmax
## 1 1 -0.34839411
## 2 1 -0.05992703
## 3 1 -0.64090492
## 4 1 -1.45357816
## 5 1 -0.88620073
## 6 1 -0.30624939
ofitted <- modfit$fitted
dim(ofitted)
## [1] 1057230 2
# a <- 1:3
# matrix(rep(a, each=4), byrow=F, ncol=3)
n <- nrow(modfit$X)
sn <- modfit$n
head(euso3)
## s.index Longitude Latitude Year Month Day o8hrmax cMAX torder sqcmax
## 1 1 -85.802 33.281 1997 5 1 59.75 25.59053 1 5.058708
## 2 1 -85.802 33.281 1997 5 2 52.75 27.03593 2 5.199608
## 3 1 -85.802 33.281 1997 5 3 51.75 24.16542 3 4.915833
## 4 1 -85.802 33.281 1997 5 4 59.38 20.42035 4 4.518888
## 5 1 -85.802 33.281 1997 5 5 66.17 23.00181 5 4.796020
## 6 1 -85.802 33.281 1997 5 6 75.50 25.79922 6 5.079294
## xsqcmax
## 1 -0.34839411
## 2 -0.05992703
## 3 -0.64090492
## 4 -1.45357816
## 5 -0.88620073
## 6 -0.30624939
r <- modfit$r
tn <- modfit$T[1]
itmax <- modfit$iterations - modfit$nBurn
n
## [1] 1057230
sn
## [1] 691
tn
## [1] 153
library(ggplot2)
oz4th<-function(x){
y<-sort(x,na.last=F)
an4th<-y[length(y)-3]
}
#
## Provides 3 year rolling averages
three_year_rolling_av <- function(x) {
n <- length(x)
y <- rep(NA, n-2)
for (k in 3:n) {
y[k-2] <- mean(x[(k-2):k])
}
y
}
head(ofitted)
## Mean SD
## [1,] 7.227619 0.5446847
## [2,] 7.126412 0.5475193
## [3,] 6.919424 0.5496207
## [4,] 7.528879 0.5401945
## [5,] 7.887672 0.5525673
## [6,] 7.928116 0.5330198
dim(ofitted)
## [1] 1057230 2
fit4th <- aggregate.data.frame(ofitted[, 1],
by=list(s.index=euso3$s.index, Year=euso3$Year), oz4th)
dim(fit4th)
## [1] 6910 3
head(fit4th)
## s.index Year x
## 1 1 1997 8.542603
## 2 2 1997 8.352123
## 3 3 1997 8.297536
## 4 4 1997 8.285875
## 5 5 1997 8.338782
## 6 6 1997 8.305429
fit4th$x <- (fit4th$x)^2
summary(fit4th$x)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 49.93 71.22 77.18 77.63 83.85 109.54
pfit1 <- ggplot(data=fit4th, aes(x=Year, y=x, color=factor(s.index))) +
geom_line() +
labs(x="Year", y="Fitted annual 4th highest ozone") +
geom_abline(slope=0, intercept=85, col="black", size=1)+
scale_colour_discrete(guide="none")
pfit1
ggsave(filename = paste0(figurepath, "fitted_an4thmax_97_06_base85.png"))
## Saving 7 x 5 in image
library(doBy)
a <- summaryBy(x ~ s.index, data = fit4th, FUN = three_year_rolling_av)
dim(a)
## [1] 691 9
b <- as.vector(t(a[,-1]))
rolling <- data.frame(s.index=rep(a$s.index, each=8), Year=rep(1999:2006, 691), ave3=b)
head(rolling)
## s.index Year ave3
## 1 1 1999 82.78449
## 2 1 2000 85.23552
## 3 1 2001 81.43925
## 4 1 2002 80.91200
## 5 1 2003 77.36652
## 6 1 2004 75.06333
pfit2 <- ggplot(data=rolling, aes(x=Year, y=ave3, color=factor(s.index))) +
geom_line() +
labs(x="Year", y="Fitted three year rolling average") +
geom_abline(slope=0, intercept=85, col="black", size=1)+
scale_colour_discrete(guide="none")
pfit2
ggsave(filename = paste0(figurepath, "fitted_3yave_99_06_base85.png"))
## Saving 7 x 5 in image
## met adjusted
library(bmstdr)
hatbeta <- get_parameter_estimates(t(modfit$betap))
hatbeta <- hatbeta$mean
u <- modfit$X
head(u)
## (Intercept) xsqcmax
## 1 1 -0.34839411
## 2 1 -0.05992703
## 3 1 -0.64090492
## 4 1 -1.45357816
## 5 1 -0.88620073
## 6 1 -0.30624939
hatbeta
## [1] 6.7371527 0.3113843
hatxbeta <- modfit$X %*% hatbeta
udat <- ofitted[, 1] - hatxbeta
dim(udat)
## [1] 1057230 1
udat <- udat^2
summary(udat)
## V1
## Min. : 0.0000
## 1st Qu.: 0.1106
## Median : 0.4754
## Mean : 0.9707
## 3rd Qu.: 1.3163
## Max. :17.9384
library(ggplot2)
met4th <- aggregate.data.frame(udat[, 1],
by=list(s.index=euso3$s.index, Year=euso3$Year), oz4th)
dim(met4th)
## [1] 6910 3
head(met4th)
## s.index Year x
## 1 1 1997 2.815393
## 2 2 1997 2.565974
## 3 3 1997 2.358718
## 4 4 1997 2.410436
## 5 5 1997 2.410536
## 6 6 1997 2.160183
summary(met4th$x)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.133 3.211 4.188 4.334 5.262 12.228
p11 <- ggplot(data=met4th, aes(x=Year, y=x, color=factor(s.index))) +
geom_line() +
labs(x="Year", y="Met adjusted annual 4th highest ozone") +
scale_colour_discrete(guide="none")
p11
ggsave(filename = paste0(figurepath, "met_adjusted_an4thmax_97_06_base85.png"))
## Saving 7 x 5 in image
library(doBy)
a <- summaryBy(x ~ s.index, data = met4th, FUN = three_year_rolling_av)
dim(a)
## [1] 691 9
b <- as.vector(t(a[,-1]))
rolling <- data.frame(s.index=rep(a$s.index, each=8), Year=rep(1999:2006, 691), ave3=b)
head(rolling)
## s.index Year ave3
## 1 1 1999 4.457021
## 2 1 2000 4.912820
## 3 1 2001 4.567063
## 4 1 2002 4.612157
## 5 1 2003 5.005674
## 6 1 2004 4.988578
p21 <- ggplot(data=rolling, aes(x=Year, y=ave3, color=factor(s.index))) +
geom_line() +
labs(x="Year", y="Met adjusted three year rolling average") +
scale_colour_discrete(guide="none")
p21
ggsave(filename = paste0(figurepath, "met_adjusted_3yave_99_06_base85.png"))
## Saving 7 x 5 in image
head(met4th)
## s.index Year x
## 1 1 1997 2.815393
## 2 2 1997 2.565974
## 3 3 1997 2.358718
## 4 4 1997 2.410436
## 5 5 1997 2.410536
## 6 6 1997 2.160183
gyear <- unique(euso3[, 1:4]) # includes year
dim(gyear)
## [1] 6910 4
head(gyear)
## s.index Longitude Latitude Year
## 1 1 -85.802 33.281 1997
## 154 1 -85.802 33.281 1998
## 307 1 -85.802 33.281 1999
## 460 1 -85.802 33.281 2000
## 613 1 -85.802 33.281 2001
## 766 1 -85.802 33.281 2002
dim(met4th)
## [1] 6910 3
mdat <- met4th[order(met4th$s.index, met4th$Year), ]
head(mdat)
## s.index Year x
## 1 1 1997 2.815393
## 692 1 1998 5.735940
## 1383 1 1999 4.819730
## 2074 1 2000 4.182790
## 2765 1 2001 4.698668
## 3456 1 2002 4.955011
adat <- merge(mdat, gyear)
head(adat)
## s.index Year x Longitude Latitude
## 1 1 1997 2.815393 -85.802 33.281
## 2 1 1998 5.735940 -85.802 33.281
## 3 1 1999 4.819730 -85.802 33.281
## 4 1 2000 4.182790 -85.802 33.281
## 5 1 2001 4.698668 -85.802 33.281
## 6 1 2002 4.955011 -85.802 33.281
adat <- adat[order(adat$s.index, adat$Year), ]
head(adat)
## s.index Year x Longitude Latitude
## 1 1 1997 2.815393 -85.802 33.281
## 2 1 1998 5.735940 -85.802 33.281
## 3 1 1999 4.819730 -85.802 33.281
## 4 1 2000 4.182790 -85.802 33.281
## 5 1 2001 4.698668 -85.802 33.281
## 6 1 2002 4.955011 -85.802 33.281
summary(adat)
## s.index Year x Longitude Latitude
## Min. : 1 Min. :1997 Min. : 1.133 Min. :-95.95 Min. :25.39
## 1st Qu.:173 1st Qu.:1999 1st Qu.: 3.211 1st Qu.:-87.42 1st Qu.:35.55
## Median :346 Median :2002 Median : 4.188 Median :-83.06 Median :39.30
## Mean :346 Mean :2002 Mean : 4.334 Mean :-82.78 Mean :38.26
## 3rd Qu.:519 3rd Qu.:2004 3rd Qu.: 5.262 3rd Qu.:-78.44 3rd Qu.:41.61
## Max. :691 Max. :2006 Max. :12.228 Max. :-68.41 Max. :48.41
adat$s.index <- as.factor(adat$s.index)
head(adat)
## s.index Year x Longitude Latitude
## 1 1 1997 2.815393 -85.802 33.281
## 2 1 1998 5.735940 -85.802 33.281
## 3 1 1999 4.819730 -85.802 33.281
## 4 1 2000 4.182790 -85.802 33.281
## 5 1 2001 4.698668 -85.802 33.281
## 6 1 2002 4.955011 -85.802 33.281
dim(adat)
## [1] 6910 5
y1 <- 1999
y2 <- 2005
idy1 <- which(adat$Year==y1)
idy2 <- which(adat$Year==y2)
length(idy1)
## [1] 691
length(idy2)
## [1] 691
useastmap <- map_data(database="state",regions=eaststates)
head(adat)
## s.index Year x Longitude Latitude
## 1 1 1997 2.815393 -85.802 33.281
## 2 1 1998 5.735940 -85.802 33.281
## 3 1 1999 4.819730 -85.802 33.281
## 4 1 2000 4.182790 -85.802 33.281
## 5 1 2001 4.698668 -85.802 33.281
## 6 1 2002 4.955011 -85.802 33.281
ptrend <- 100*(adat[idy2, "x"] - adat[idy1, "x"])/adat[idy1, "x"]
summary(ptrend)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -69.02 -38.96 -27.19 -13.29 -12.46 237.31
v <- cbind(adat[idy2, "x"], adat[idy1, "x"], ptrend)
head(v)
## ptrend
## [1,] 5.371857 4.819730 11.45554
## [2,] 5.417788 4.129024 31.21231
## [3,] 4.661045 3.625279 28.57064
## [4,] 5.026614 3.448682 45.75463
## [5,] 4.985810 3.625865 37.50679
## [6,] 4.190312 3.400243 23.23566
ptrend1 <- 0.5 * (ptrend - min(ptrend))/(max(ptrend) - min(ptrend))
head(gyear)
## s.index Longitude Latitude Year
## 1 1 -85.802 33.281 1997
## 154 1 -85.802 33.281 1998
## 307 1 -85.802 33.281 1999
## 460 1 -85.802 33.281 2000
## 613 1 -85.802 33.281 2001
## 766 1 -85.802 33.281 2002
ggrid <- unique(gyear[, 1:3])
ptsum <- data.frame(ggrid, ptrend =ptrend, rtrend=ptrend1)
ptsum$rtrend[ptsum$ptrend<0] <- ptsum$rtrend[ptsum$ptrend<0] * 2
ptsum$col <- sign(ptsum$ptrend)
table(ptsum$col)
##
## -1 1
## 571 120
ptsum$trend <- factor(ptsum$col, levels=c("-1", "1"), labels=c("-ve", "+ve"))
trendp <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r = rtrend, fill=trend), inherit.aes = FALSE,data = ptsum) +
scale_fill_manual(values = c("#006600", "#FF0000")) +
scale_color_manual(values = c("#006600", "#FF0000")) +
labs( title= paste("Met adjusted RPCT for ", y2, " (base ",y1, ")"), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank()) +
ggsn::north(data=useastmap, location="bottomright", symbol=12) +
ggsn::scalebar(data =useastmap, dist =500, location ="bottomleft", transform=T, dist_unit = "km",
st.dist = .05, st.size = 5, height = .06, st.bottom=T, model="WGS84")
trendp
ggsave(filename = paste0(figurepath, "met_adjusted_4thmax_99_05.png"))
## Saving 7 x 5 in image
## Above plot for rolling 3 year averages
head(rolling)
## s.index Year ave3
## 1 1 1999 4.457021
## 2 1 2000 4.912820
## 3 1 2001 4.567063
## 4 1 2002 4.612157
## 5 1 2003 5.005674
## 6 1 2004 4.988578
dim(rolling)
## [1] 5528 3
rgyear <- unique(euso3[, 1:4]) # includes year
dim(rgyear)
## [1] 6910 4
rgyear <- rgyear[rgyear$Year>1998, ]
dim(rgyear)
## [1] 5528 4
head(rgyear)
## s.index Longitude Latitude Year
## 307 1 -85.802 33.281 1999
## 460 1 -85.802 33.281 2000
## 613 1 -85.802 33.281 2001
## 766 1 -85.802 33.281 2002
## 919 1 -85.802 33.281 2003
## 1072 1 -85.802 33.281 2004
rdat <- rolling[order(rolling$s.index, rolling$Year), ]
head(rdat)
## s.index Year ave3
## 1 1 1999 4.457021
## 2 1 2000 4.912820
## 3 1 2001 4.567063
## 4 1 2002 4.612157
## 5 1 2003 5.005674
## 6 1 2004 4.988578
adat <- merge(rdat, rgyear)
head(adat)
## s.index Year ave3 Longitude Latitude
## 1 1 1999 4.457021 -85.802 33.281
## 2 1 2000 4.912820 -85.802 33.281
## 3 1 2001 4.567063 -85.802 33.281
## 4 1 2002 4.612157 -85.802 33.281
## 5 1 2003 5.005674 -85.802 33.281
## 6 1 2004 4.988578 -85.802 33.281
adat <- adat[order(adat$s.index, adat$Year), ]
head(adat)
## s.index Year ave3 Longitude Latitude
## 1 1 1999 4.457021 -85.802 33.281
## 2 1 2000 4.912820 -85.802 33.281
## 3 1 2001 4.567063 -85.802 33.281
## 4 1 2002 4.612157 -85.802 33.281
## 5 1 2003 5.005674 -85.802 33.281
## 6 1 2004 4.988578 -85.802 33.281
summary(adat)
## s.index Year ave3 Longitude Latitude
## Min. : 1 Min. :1999 Min. :1.477 Min. :-95.95 Min. :25.39
## 1st Qu.:173 1st Qu.:2001 1st Qu.:3.647 1st Qu.:-87.42 1st Qu.:35.55
## Median :346 Median :2002 Median :4.445 Median :-83.06 Median :39.30
## Mean :346 Mean :2002 Mean :4.478 Mean :-82.78 Mean :38.26
## 3rd Qu.:519 3rd Qu.:2004 3rd Qu.:5.271 3rd Qu.:-78.44 3rd Qu.:41.61
## Max. :691 Max. :2006 Max. :9.277 Max. :-68.41 Max. :48.41
adat$s.index <- as.factor(adat$s.index)
head(adat)
## s.index Year ave3 Longitude Latitude
## 1 1 1999 4.457021 -85.802 33.281
## 2 1 2000 4.912820 -85.802 33.281
## 3 1 2001 4.567063 -85.802 33.281
## 4 1 2002 4.612157 -85.802 33.281
## 5 1 2003 5.005674 -85.802 33.281
## 6 1 2004 4.988578 -85.802 33.281
dim(adat)
## [1] 5528 5
y1 <- 1999
y2 <- 2005
idy1 <- which(adat$Year==y1)
idy2 <- which(adat$Year==y2)
length(idy1)
## [1] 691
length(idy2)
## [1] 691
head(v)
## ptrend
## [1,] 5.371857 4.819730 11.45554
## [2,] 5.417788 4.129024 31.21231
## [3,] 4.661045 3.625279 28.57064
## [4,] 5.026614 3.448682 45.75463
## [5,] 4.985810 3.625865 37.50679
## [6,] 4.190312 3.400243 23.23566
head(adat)
## s.index Year ave3 Longitude Latitude
## 1 1 1999 4.457021 -85.802 33.281
## 2 1 2000 4.912820 -85.802 33.281
## 3 1 2001 4.567063 -85.802 33.281
## 4 1 2002 4.612157 -85.802 33.281
## 5 1 2003 5.005674 -85.802 33.281
## 6 1 2004 4.988578 -85.802 33.281
ptrend <- 100*(adat[idy2, "ave3"] - adat[idy1, "ave3"])/adat[idy1, "ave3"]
summary(ptrend)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -44.974 -28.395 -20.552 -6.301 -3.646 175.630
v <- cbind(adat[idy2, "ave3"], adat[idy1, "ave3"], ptrend)
head(v)
## ptrend
## [1,] 5.127526 4.457021 15.04378
## [2,] 5.866974 3.958567 48.20954
## [3,] 4.946557 3.264785 51.51249
## [4,] 5.404193 3.252318 66.16435
## [5,] 5.246021 3.332351 57.42703
## [6,] 4.173957 3.249867 28.43472
ptrend1 <- 0.5 * (ptrend - min(ptrend))/(max(ptrend) - min(ptrend))
head(rgyear)
## s.index Longitude Latitude Year
## 307 1 -85.802 33.281 1999
## 460 1 -85.802 33.281 2000
## 613 1 -85.802 33.281 2001
## 766 1 -85.802 33.281 2002
## 919 1 -85.802 33.281 2003
## 1072 1 -85.802 33.281 2004
ggrid <- unique(rgyear[, 1:3])
ptsum <- data.frame(ggrid, ptrend =ptrend, rtrend=ptrend1)
ptsum$rtrend[ptsum$ptrend<0] <- ptsum$rtrend[ptsum$ptrend<0] *2
ptsum$col <- sign(ptsum$ptrend)
table(ptsum$col)
##
## -1 1
## 544 147
ptsum$trend <- factor(ptsum$col, levels=c("-1", "1"), labels=c("-ve", "+ve"))
rtrendp <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
scale_fill_manual(values = c("#006600", "#FF0000")) +
scale_color_manual(values = c("#006600", "#FF0000")) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r = rtrend, fill=trend), inherit.aes = FALSE,data = ptsum) +
labs( title= paste("Met adjusted RPCT in rolling averages for ", y2, " (base ",y1, ")"), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank()) +
ggsn::north(data=useastmap, location="bottomright", symbol=12) +
ggsn::scalebar(data =useastmap, dist =500, location ="bottomleft", transform=T, dist_unit = "km",
st.dist = .05, st.size = 5, height = .06, st.bottom=T, model="WGS84")
rtrendp
ggsave(filename = paste0(figurepath, "met_adjusted_3ave_99_05.png"))
## Saving 7 x 5 in image
## Probability of compliance
head(euso3)
## s.index Longitude Latitude Year Month Day o8hrmax cMAX torder sqcmax
## 1 1 -85.802 33.281 1997 5 1 59.75 25.59053 1 5.058708
## 2 1 -85.802 33.281 1997 5 2 52.75 27.03593 2 5.199608
## 3 1 -85.802 33.281 1997 5 3 51.75 24.16542 3 4.915833
## 4 1 -85.802 33.281 1997 5 4 59.38 20.42035 4 4.518888
## 5 1 -85.802 33.281 1997 5 5 66.17 23.00181 5 4.796020
## 6 1 -85.802 33.281 1997 5 6 75.50 25.79922 6 5.079294
## xsqcmax
## 1 -0.34839411
## 2 -0.05992703
## 3 -0.64090492
## 4 -1.45357816
## 5 -0.88620073
## 6 -0.30624939
a <- (sqrt(euso3$o8hrmax) - ofitted[,1])^2
sqrt(mean(a, na.rm=T))
## [1] 0.4570944
head(ofitted)
## Mean SD
## [1,] 7.227619 0.5446847
## [2,] 7.126412 0.5475193
## [3,] 6.919424 0.5496207
## [4,] 7.528879 0.5401945
## [5,] 7.887672 0.5525673
## [6,] 7.928116 0.5330198
class(ofitted)
## [1] "matrix" "array"
ofitted <- data.frame(ofitted)
ofitted$rownumber <- 1:nrow(ofitted)
ofitted$s.index <- euso3$s.index
ofitted$Year <- euso3$Year
ofitted$time <- rep(1:153, 6910)
head(ofitted)
## Mean SD rownumber s.index Year time
## 1 7.227619 0.5446847 1 1 1997 1
## 2 7.126412 0.5475193 2 1 1997 2
## 3 6.919424 0.5496207 3 1 1997 3
## 4 7.528879 0.5401945 4 1 1997 4
## 5 7.887672 0.5525673 5 1 1997 5
## 6 7.928116 0.5330198 6 1 1997 6
head(euso3)
## s.index Longitude Latitude Year Month Day o8hrmax cMAX torder sqcmax
## 1 1 -85.802 33.281 1997 5 1 59.75 25.59053 1 5.058708
## 2 1 -85.802 33.281 1997 5 2 52.75 27.03593 2 5.199608
## 3 1 -85.802 33.281 1997 5 3 51.75 24.16542 3 4.915833
## 4 1 -85.802 33.281 1997 5 4 59.38 20.42035 4 4.518888
## 5 1 -85.802 33.281 1997 5 5 66.17 23.00181 5 4.796020
## 6 1 -85.802 33.281 1997 5 6 75.50 25.79922 6 5.079294
## xsqcmax
## 1 -0.34839411
## 2 -0.05992703
## 3 -0.64090492
## 4 -1.45357816
## 5 -0.88620073
## 6 -0.30624939
oz4throw <- function(x){
k <- length(x)
y <- data.frame(x=x, number=1:k)
y <- y[order(y$x), ]
y$number[k-3]
}
library(doBy)
ids <- summaryBy(Mean ~ s.index+Year, data = ofitted, FUN = oz4throw)
head(ids)
## s.index Year Mean.oz4throw
## 1 1 1997 142
## 2 1 1998 121
## 3 1 1999 127
## 4 1 2000 121
## 5 1 2001 49
## 6 1 2002 133
colnames(ids)[3] <- "time"
head(ids)
## s.index Year time
## 1 1 1997 142
## 2 1 1998 121
## 3 1 1999 127
## 4 1 2000 121
## 5 1 2001 49
## 6 1 2002 133
u <- merge(ids, ofitted)
dim(u) # 6910 rows
## [1] 6910 6
head(u)
## s.index Year time Mean SD rownumber
## 1 1 1997 142 8.542603 0.5432913 142
## 2 1 1998 121 9.441528 0.5462202 274
## 3 1 1999 127 9.286278 0.5464004 433
## 4 1 2000 121 8.962653 0.5515905 580
## 5 1 2001 49 8.817803 0.5377157 661
## 6 1 2002 133 9.200718 0.5519402 898
highsum <- cbind(u$Mean, u$SD)
class(highsum)
## [1] "matrix" "array"
dim(highsum)
## [1] 6910 2
head(highsum)
## [,1] [,2]
## [1,] 8.542603 0.5432913
## [2,] 9.441528 0.5462202
## [3,] 9.286278 0.5464004
## [4,] 8.962653 0.5515905
## [5,] 8.817803 0.5377157
## [6,] 9.200718 0.5519402
gyear <- unique(euso3[, 1:4])
head(gyear)
## s.index Longitude Latitude Year
## 1 1 -85.802 33.281 1997
## 154 1 -85.802 33.281 1998
## 307 1 -85.802 33.281 1999
## 460 1 -85.802 33.281 2000
## 613 1 -85.802 33.281 2001
## 766 1 -85.802 33.281 2002
dim(gyear)
## [1] 6910 4
head(euso3)
## s.index Longitude Latitude Year Month Day o8hrmax cMAX torder sqcmax
## 1 1 -85.802 33.281 1997 5 1 59.75 25.59053 1 5.058708
## 2 1 -85.802 33.281 1997 5 2 52.75 27.03593 2 5.199608
## 3 1 -85.802 33.281 1997 5 3 51.75 24.16542 3 4.915833
## 4 1 -85.802 33.281 1997 5 4 59.38 20.42035 4 4.518888
## 5 1 -85.802 33.281 1997 5 5 66.17 23.00181 5 4.796020
## 6 1 -85.802 33.281 1997 5 6 75.50 25.79922 6 5.079294
## xsqcmax
## 1 -0.34839411
## 2 -0.05992703
## 3 -0.64090492
## 4 -1.45357816
## 5 -0.88620073
## 6 -0.30624939
euso3$time <- rep(1:153, 6910)
v <- merge(ids, euso3)
dim(v)
## [1] 6910 12
head(v)
## s.index Year time Longitude Latitude Month Day o8hrmax cMAX torder
## 1 1 1997 142 -85.802 33.281 9 19 69.25 35.07221 142
## 2 1 1998 121 -85.802 33.281 8 29 99.00 36.06224 274
## 3 1 1999 127 -85.802 33.281 9 4 84.00 33.93251 433
## 4 1 2000 121 -85.802 33.281 8 29 80.50 35.32145 580
## 5 1 2001 49 -85.802 33.281 6 18 68.25 32.42145 661
## 6 1 2002 133 -85.802 33.281 9 10 83.00 35.85556 898
## sqcmax xsqcmax
## 1 5.922180 1.4194059
## 2 6.005184 1.5893424
## 3 5.825162 1.2207798
## 4 5.943185 1.4624108
## 5 5.693983 0.9522157
## 6 5.987952 1.5540616
v <- v[order(v$s.index, v$Year, v$torder), ]
head(v)
## s.index Year time Longitude Latitude Month Day o8hrmax cMAX torder
## 1 1 1997 142 -85.802 33.281 9 19 69.25 35.07221 142
## 2 1 1998 121 -85.802 33.281 8 29 99.00 36.06224 274
## 3 1 1999 127 -85.802 33.281 9 4 84.00 33.93251 433
## 4 1 2000 121 -85.802 33.281 8 29 80.50 35.32145 580
## 5 1 2001 49 -85.802 33.281 6 18 68.25 32.42145 661
## 6 1 2002 133 -85.802 33.281 9 10 83.00 35.85556 898
## sqcmax xsqcmax
## 1 5.922180 1.4194059
## 2 6.005184 1.5893424
## 3 5.825162 1.2207798
## 4 5.943185 1.4624108
## 5 5.693983 0.9522157
## 6 5.987952 1.5540616
sum4thstat <- cbind.data.frame(v, fitmean=highsum[,1], fitsd=highsum[,2])
head(sum4thstat)
## s.index Year time Longitude Latitude Month Day o8hrmax cMAX torder
## 1 1 1997 142 -85.802 33.281 9 19 69.25 35.07221 142
## 2 1 1998 121 -85.802 33.281 8 29 99.00 36.06224 274
## 3 1 1999 127 -85.802 33.281 9 4 84.00 33.93251 433
## 4 1 2000 121 -85.802 33.281 8 29 80.50 35.32145 580
## 5 1 2001 49 -85.802 33.281 6 18 68.25 32.42145 661
## 6 1 2002 133 -85.802 33.281 9 10 83.00 35.85556 898
## sqcmax xsqcmax fitmean fitsd
## 1 5.922180 1.4194059 8.542603 0.5432913
## 2 6.005184 1.5893424 9.441528 0.5462202
## 3 5.825162 1.2207798 9.286278 0.5464004
## 4 5.943185 1.4624108 8.962653 0.5515905
## 5 5.693983 0.9522157 8.817803 0.5377157
## 6 5.987952 1.5540616 9.200718 0.5519402
# Not used any more
Rcpp::sourceCpp("metadjusted.cpp")
a <- three_year_rolling(high4summary=highsum, sn=sn, r=r, itmax=1000)
dim(a)
## [1] 5528 1000
# head(a)
pgreatsqrty85 <- function(x) {
length(x[x>sqrt(85)])/length(x)
}
v <- apply(a, 1, pgreatsqrty85)
summary(v)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0050 0.3205 0.2841 0.4620 0.9860
gyear <- unique(euso3[, 1:4])
head(gyear)
## s.index Longitude Latitude Year
## 1 1 -85.802 33.281 1997
## 154 1 -85.802 33.281 1998
## 307 1 -85.802 33.281 1999
## 460 1 -85.802 33.281 2000
## 613 1 -85.802 33.281 2001
## 766 1 -85.802 33.281 2002
dim(gyear)
## [1] 6910 4
rgyear <- gyear[gyear$Year>1998, ]
dim(rgyear)
## [1] 5528 4
head(rgyear)
## s.index Longitude Latitude Year
## 307 1 -85.802 33.281 1999
## 460 1 -85.802 33.281 2000
## 613 1 -85.802 33.281 2001
## 766 1 -85.802 33.281 2002
## 919 1 -85.802 33.281 2003
## 1072 1 -85.802 33.281 2004
rgyear$pg85 <- v
head(rgyear)
## s.index Longitude Latitude Year pg85
## 307 1 -85.802 33.281 1999 0.367
## 460 1 -85.802 33.281 2000 0.514
## 613 1 -85.802 33.281 2001 0.275
## 766 1 -85.802 33.281 2002 0.243
## 919 1 -85.802 33.281 2003 0.083
## 1072 1 -85.802 33.281 2004 0.039
hist(v)
b <- c(-Inf, 0.25, 0.5, 0.75, 1.01)
labss <- c("0-0.25", "0.25-0.50", "0.50-0.75", "0.75-1")
# b <- c(-Inf, 0.5, 0.90, 1.01)
# labss <- c("0-0.5", "0.5-0.9", "0.9-1")
rgyear$compliant <- cut(rgyear$pg85, breaks = b, labels =labss)
table(rgyear$compliant)
##
## 0-0.25 0.25-0.50 0.50-0.75 0.75-1
## 2333 2209 916 70
com <- c("green4", "dodgerblue4", "orange", "red")
rgyear$pnew85 <- rgyear$pg85
# rgyear$pnew85[rgyear$pnew85<0.75] <- rgyear$pg85[rgyear$pnew85<0.75] * 2
# rgyear$pnew85[rgyear$pnew85<0.5] <- rgyear$pg85[rgyear$pnew85<0.5] * 2
rgyear$pnew85 <- rgyear$pnew85
y1 <- 1999
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 307 1 -85.802 33.281 1999 0.367 0.25-0.50 0.367
## 1837 2 -86.137 32.498 1999 0.022 0-0.25 0.022
## 3367 3 -86.915 33.486 1999 0.337 0.25-0.50 0.337
## 4897 4 -87.004 33.331 1999 0.353 0.25-0.50 0.353
## 6427 5 -86.817 33.386 1999 0.304 0.25-0.50 0.304
## 7957 6 -86.669 33.705 1999 0.289 0.25-0.50 0.289
complot99 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot99
y1 <- 2000
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 460 1 -85.802 33.281 2000 0.514 0.50-0.75 0.514
## 1990 2 -86.137 32.498 2000 0.119 0-0.25 0.119
## 3520 3 -86.915 33.486 2000 0.382 0.25-0.50 0.382
## 5050 4 -87.004 33.331 2000 0.349 0.25-0.50 0.349
## 6580 5 -86.817 33.386 2000 0.369 0.25-0.50 0.369
## 8110 6 -86.669 33.705 2000 0.197 0-0.25 0.197
complot00 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot00
y1 <- 2001
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 613 1 -85.802 33.281 2001 0.275 0.25-0.50 0.275
## 2143 2 -86.137 32.498 2001 0.061 0-0.25 0.061
## 3673 3 -86.915 33.486 2001 0.313 0.25-0.50 0.313
## 5203 4 -87.004 33.331 2001 0.333 0.25-0.50 0.333
## 6733 5 -86.817 33.386 2001 0.365 0.25-0.50 0.365
## 8263 6 -86.669 33.705 2001 0.178 0-0.25 0.178
complot01 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot01
y1 <- 2002
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 766 1 -85.802 33.281 2002 0.243 0-0.25 0.243
## 2296 2 -86.137 32.498 2002 0.049 0-0.25 0.049
## 3826 3 -86.915 33.486 2002 0.507 0.50-0.75 0.507
## 5356 4 -87.004 33.331 2002 0.537 0.50-0.75 0.537
## 6886 5 -86.817 33.386 2002 0.576 0.50-0.75 0.576
## 8416 6 -86.669 33.705 2002 0.184 0-0.25 0.184
complot02 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot02
y1 <- 2003
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 919 1 -85.802 33.281 2003 0.083 0-0.25 0.083
## 2449 2 -86.137 32.498 2003 0.012 0-0.25 0.012
## 3979 3 -86.915 33.486 2003 0.520 0.50-0.75 0.520
## 5509 4 -87.004 33.331 2003 0.562 0.50-0.75 0.562
## 7039 5 -86.817 33.386 2003 0.555 0.50-0.75 0.555
## 8569 6 -86.669 33.705 2003 0.291 0.25-0.50 0.291
complot03 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot03
y1 <- 2004
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 1072 1 -85.802 33.281 2004 0.039 0-0.25 0.039
## 2602 2 -86.137 32.498 2004 0.012 0-0.25 0.012
## 4132 3 -86.915 33.486 2004 0.360 0.25-0.50 0.360
## 5662 4 -87.004 33.331 2004 0.427 0.25-0.50 0.427
## 7192 5 -86.817 33.386 2004 0.343 0.25-0.50 0.343
## 8722 6 -86.669 33.705 2004 0.189 0-0.25 0.189
complot04 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot04
y1 <- 2005
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 1225 1 -85.802 33.281 2005 0 0-0.25 0
## 2755 2 -86.137 32.498 2005 0 0-0.25 0
## 4285 3 -86.915 33.486 2005 0 0-0.25 0
## 5815 4 -87.004 33.331 2005 0 0-0.25 0
## 7345 5 -86.817 33.386 2005 0 0-0.25 0
## 8875 6 -86.669 33.705 2005 0 0-0.25 0
complot05 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot05
y1 <- 2006
idy <- which(rgyear$Year==y1)
probmat <- rgyear[idy, ]
head(probmat)
## s.index Longitude Latitude Year pg85 compliant pnew85
## 1378 1 -85.802 33.281 2006 0 0-0.25 0
## 2908 2 -86.137 32.498 2006 0 0-0.25 0
## 4438 3 -86.915 33.486 2006 0 0-0.25 0
## 5968 4 -87.004 33.331 2006 0 0-0.25 0
## 7498 5 -86.817 33.386 2006 0 0-0.25 0
## 9028 6 -86.669 33.705 2006 0 0-0.25 0
complot06 <- ggplot(data=useastmap, aes(x=long, y=lat, group = group)) +
geom_polygon(colour='black',size=0.8, fill=NA) +
geom_circle(aes(x0 = Longitude, y0 = Latitude, r =pnew85, col=compliant, fill=compliant), inherit.aes = FALSE,data = probmat) +
scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE))+
scale_color_manual(values =com) +
guides(colour=FALSE) +
labs( title= paste("Probability of non-compliance: ", y1), x="", y = "") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
complot06
png(filename = paste0(figurepath, "compliance_plot_05.png"))
complot05
dev.off()
## png
## 2
png(filename = paste0(figurepath, "compliance_plot_06.png"))
complot06
dev.off()
## png
## 2
library(ggpubr)
a <- get_legend(complot99, position = "bottom")
ggarrange(
complot99, complot00, complot01, complot02,
common.legend = TRUE, legend = "top", nrow = 2, ncol = 2)
ggsave(filename=paste0(figurepath, "prob99-02.png"), device = "png")
## Saving 7 x 5 in image
ggarrange(
complot03, complot04, complot05, complot06,
common.legend = TRUE, legend = "top", nrow = 2, ncol = 2)
ggsave(filename=paste0(figurepath, "prob03-06.png"), device = "png")
## Saving 7 x 5 in image
##
# All done
end.time <- proc.time()[3]
comp.time<- (end.time-start.time)/60
# comp.time<-fancy.time(comp.time)
print(comp.time)
## elapsed
## 1.965983