Code to reproduce Table 7.1
head(nysptime)
## s.index Longitude Latitude utmx utmy Year Month Day y8hrmax xmaxtemp
## 1 1 -73.757 42.681 601838.1 4726140 2006 7 1 53.88 27.85772
## 2 1 -73.757 42.681 601838.1 4726140 2006 7 2 57.13 30.11563
## 3 1 -73.757 42.681 601838.1 4726140 2006 7 3 72.00 30.00001
## 4 1 -73.757 42.681 601838.1 4726140 2006 7 4 36.63 27.89656
## 5 1 -73.757 42.681 601838.1 4726140 2006 7 5 42.63 25.65698
## 6 1 -73.757 42.681 601838.1 4726140 2006 7 6 30.88 24.61968
## xwdsp xrh
## 1 5.459953 2.766221
## 2 8.211767 3.197750
## 3 4.459581 3.225186
## 4 3.692225 4.362334
## 5 4.374314 3.950320
## 6 4.178086 3.420533
M1 <- Bsptime(model="lm", formula=f2, data=nysptime, scale.transform = "SQRT",
N=5000)
## ##
## # Total time taken:: 0.01 - Sec.
summary(M1)
##
## The lm model has been fitted using bmstdr code in R.
## Call:
## Bsptime(formula = f2, data = nysptime, model = "lm", scale.transform = "SQRT",
## N = 5000)
## ##
## # Total time taken:: 0.01 - Sec.
## Model formula
## y8hrmax ~ xmaxtemp + xwdsp + xrh
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) 2.223 0.240 1.754 2.693
## xmaxtemp 0.175 0.007 0.162 0.188
## xwdsp 0.089 0.012 0.065 0.113
## xrh -0.173 0.028 -0.228 -0.119
## sigma2 0.555 0.019 0.519 0.593
M1$tn
## [1] 0
M1$sn
## [1] 1736
##
plot(M1)
##
## No other plots implemented for this model fitting method.
##
## Note that the residuals are provided on the transformed scale. Please see the scale.transform argument.
##
## Summary of the residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -2.364104 -0.527171 0.001326 -0.002186 0.501164 2.514540 24
## tn has not been supplied in residuals and it is
## not possible to figure this out. Hence a time series plot of the residuals
## has not been provided here.
##
## Note that the residuals are provided on the transformed scale.
## Please see the scale.transform argument.
a <- residuals(M1, numbers=list(sn=28, tn=62))
##
## Note that the residuals are provided on the transformed scale. Please see the scale.transform argument.
##
## Summary of the residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -2.364104 -0.527171 0.001326 -0.002186 0.501164 2.514540 24
a <- residuals(M1, numbers=list(sn=28))
##
## Note that the residuals are provided on the transformed scale. Please see the scale.transform argument.
##
##
## Summary of the residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -2.364104 -0.527171 0.001326 -0.002186 0.501164 2.514540 24
## tn has not been supplied in residuals and it is
## not possible to figure this out. Hence a time series plot of the residuals
## has not been provided here.
## Spatio-temporal model fitting
M2 <- Bsptime(model="separable", formula=f2, data=nysptime, scale.transform = "SQRT",
coordtype="utm", coords=4:5, N=5000)
## ##
## # Total time taken:: 0.01 - Sec.
names(M2)
## [1] "params" "fit" "max.d" "fitteds"
## [5] "phi.s" "phi.t" "residuals" "sn"
## [9] "tn" "formula" "scale.transform" "package"
## [13] "model" "call" "computation.time"
summary(M2)
##
## The separable model has been fitted using bmstdr code in R.
## Call:
## Bsptime(formula = f2, data = nysptime, model = "separable", coordtype = "utm",
## coords = 4:5, scale.transform = "SQRT", N = 5000)
## ##
## # Total time taken:: 0.01 - Sec.
## Model formula
## y8hrmax ~ xmaxtemp + xwdsp + xrh
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) 0.191 1.740 -3.219 3.602
## xmaxtemp 0.259 0.037 0.187 0.331
## xwdsp 0.010 0.035 -0.059 0.079
## xrh -0.011 0.109 -0.226 0.203
## sigma2 12.189 0.414 11.405 13.027
M2$phi.s
## [1] 0.005072598
M2$phi.t
## [1] 0.0483871
table7.1 <- cbind(M1$params, M2$params)
dput(table7.1, file=paste0(filepath, "table7.1.txt"))
print(round(table7.1, 2))
## mean sd 2.5% 97.5% mean sd 2.5% 97.5%
## (Intercept) 2.22 0.24 1.75 2.69 0.19 1.74 -3.22 3.60
## xmaxtemp 0.17 0.01 0.16 0.19 0.26 0.04 0.19 0.33
## xwdsp 0.09 0.01 0.06 0.11 0.01 0.04 -0.06 0.08
## xrh -0.17 0.03 -0.23 -0.12 -0.01 0.11 -0.23 0.20
## sigma2 0.55 0.02 0.52 0.59 12.19 0.41 11.40 13.03
# xtable(table7.1, digits=2)
# a <- phichoicep()
valids <- c(8,11,12,14,18,21,24,28)
vrows <- which(nysptime$s.index%in% valids)
M2.1 <- Bsptime(model="separable", formula=f2, data=nysptime,
validrows=vrows, coordtype="utm", coords=4:5,
phi.s=0.005, phi.t=0.05, scale.transform = "SQRT")
## ##
## # Total time taken:: 7.61 - Sec.
names(M2.1)
## [1] "params" "fit" "max.d" "fitteds"
## [5] "stats" "yobs_preds" "valpreds" "phi.s"
## [9] "phi.t" "residuals" "sn" "tn"
## [13] "formula" "scale.transform" "package" "model"
## [17] "call" "computation.time"
dim(M2.1$valpreds)
## [1] 2000 496
u <- get_validation_summaries(M2.1$valpreds)
dim(u)
## [1] 496 5
head(u)
## meanpred sdpred medianpred low up
## 1 67.86876 11.203631 66.89214 48.59173 91.42743
## 2 66.02045 11.012413 65.57377 45.88358 89.08762
## 3 58.56413 10.636774 58.43256 39.25781 80.54063
## 4 58.23229 10.559617 57.67607 39.85306 80.69852
## 5 42.09230 8.741677 41.86654 26.63715 59.88194
## 6 31.87050 7.653939 31.41477 18.79720 48.34924
M1 <- Bsptime(model="lm", formula=f2, data=nysptime,
validrows=vrows, scale.transform = "SQRT")
## ##
## # Total time taken:: 7.14 - Sec.
summary(M1)
##
## The lm model has been fitted using bmstdr code in R.
## Call:
## Bsptime(formula = f2, data = nysptime, model = "lm", validrows = vrows,
## scale.transform = "SQRT")
## ##
## # Total time taken:: 7.14 - Sec.
## Model formula
## y8hrmax ~ xmaxtemp + xwdsp + xrh
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) 2.303 0.293 1.728 2.878
## xmaxtemp 0.174 0.008 0.159 0.190
## xwdsp 0.089 0.015 0.061 0.118
## xrh -0.190 0.034 -0.256 -0.125
## sigma2 0.588 0.024 0.543 0.636
##
## Validation Statistics:
## rmse mae crps cvg
## 9.367 7.549 5.627 97.951
summary(M2.1)
##
## The separable model has been fitted using bmstdr code in R.
## Call:
## Bsptime(formula = f2, data = nysptime, model = "separable", coordtype = "utm",
## coords = 4:5, validrows = vrows, scale.transform = "SQRT",
## phi.s = 0.005, phi.t = 0.05)
## ##
## # Total time taken:: 7.61 - Sec.
## Model formula
## y8hrmax ~ xmaxtemp + xwdsp + xrh
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) 0.164 1.652 -3.074 3.402
## xmaxtemp 0.258 0.036 0.188 0.328
## xwdsp 0.002 0.034 -0.065 0.068
## xrh 0.008 0.107 -0.201 0.217
## sigma2 10.334 0.415 9.552 11.179
##
## Validation Statistics:
## rmse mae crps cvg
## 6.497 5.013 10.670 99.590
plot(M2.1)
##
## No other plots implemented for this model fitting method.
##
## Note that the residuals are provided on the transformed scale. Please see the scale.transform argument.
## Validation has been performed. The residuals include the validation observations as well.
## Expect the return value to be of the same length as the supplied data frame.
##
## Summary of the residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -3.11069 -0.97055 -0.47227 -0.47438 0.04332 2.13322 24
##
## Note that the residuals are provided on the transformed scale.
## Please see the scale.transform argument.
## Validation has been performed
## The fitted values include the validation observations as well.
## Expect the return value to be of the same length as the supplied data frame.
##
## Validation Statistics:
## rmse mae crps cvg
## 6.497 5.013 10.670 99.590
## $pwithseg
##
## $pwithoutseg
##
## $pordinary
plot(M2.1, segments = F)
##
## No other plots implemented for this model fitting method.
##
## Note that the residuals are provided on the transformed scale. Please see the scale.transform argument.
## Validation has been performed. The residuals include the validation observations as well.
## Expect the return value to be of the same length as the supplied data frame.
##
## Summary of the residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -3.11069 -0.97055 -0.47227 -0.47438 0.04332 2.13322 24
##
## Note that the residuals are provided on the transformed scale.
## Please see the scale.transform argument.
## Validation has been performed
## The fitted values include the validation observations as well.
## Expect the return value to be of the same length as the supplied data frame.
##
## Validation Statistics:
## rmse mae crps cvg
## 6.497 5.013 10.670 99.590
## $pwithseg
##
## $pwithoutseg
##
## $pordinary
f2 <- y8hrmax~xmaxtemp+xwdsp+xrh
M3 <- Bsptime(package="spTimer", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
N=5000, n.report=1)
##
## Output: GP models
## ---------------------------------------------------------------
## Sampled: 5000 of 5000, 100.00%.
## Batch Acceptance Rate (phi): 37.45%
## Checking Parameters:
## phi: 0.0163, sig2eps: 0.0165, sig2eta: 0.6449
## beta[1]: 1.1776 beta[2]: 0.1870 beta[3]: 0.1530 beta[4]: -0.0405
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 , Iterations = 5000 .
## # Overall Acceptance Rate (phi) = 37.44 %
## ##
## ##
## # Elapsed time: 6.16 Sec.
## ##
##
## # Model: GP
table7.2 <- M3$params
round(table7.2, 3)
## mean sd 2.5% 97.5%
## (Intercept) 2.140 0.452 1.259 3.031
## xmaxtemp 0.171 0.013 0.145 0.196
## xwdsp 0.114 0.021 0.073 0.155
## xrh -0.126 0.052 -0.225 -0.019
## sig2eps 0.018 0.004 0.012 0.026
## sig2eta 0.586 0.071 0.484 0.758
## phi 0.014 0.002 0.011 0.017
dput(table7.2, file=paste0(filepath, "table7.2.txt"))
nymap <- map_data(database="state",regions="new york")
head(nymap)
## long lat group order region subregion
## 1 -73.92874 40.80605 1 1 new york manhattan
## 2 -73.93448 40.78886 1 2 new york manhattan
## 3 -73.95166 40.77741 1 3 new york manhattan
## 4 -73.96312 40.75449 1 4 new york manhattan
## 5 -73.96885 40.73730 1 5 new york manhattan
## 6 -73.97458 40.72584 1 6 new york manhattan
s <- c(1, 5, 10)
fcoords <- nyspatial[-s, c("Longitude", "Latitude")]
vcoords <- nyspatial[s, c("Longitude", "Latitude")]
library(tidyr)
label <- tibble(
long = -76.1,
lat = 41.5,
label = "25 fitted (circles) & 3 \n validation (numbered) sites"
)
vsites3 <- ggplot() +
geom_polygon(data=nymap, aes(x=long, y=lat, group=group),
color="black", size = 0.6, fill=NA) +
geom_point(data =fcoords, aes(x=Longitude,y=Latitude)) +
geom_text(data=vcoords, aes(x=Longitude,y=Latitude, label=s), col=4) +
labs( title= "28 air pollution monitoring sites in New York", x="Longitude", y = "Latitude") +
geom_text(aes(label=label, x=long, y=lat), data = label, vjust = "top", hjust = "right") +
# geom_rect(mapping=aes(xmin=-80.2, xmax=-77.3, ymin=41, ymax=41.6), color="black", fill=NA) +
geom_rect(mapping=aes(xmin=-78.7, xmax=-75.8, ymin=41, ymax=41.6), color="black", fill=NA) +
ggsn::scalebar(data =nymap, dist = 100, location = "bottomleft", transform=T, dist_unit = "km",
st.dist = .05, st.size = 5, height = .06, st.bottom=T, model="WGS84") +
ggsn::north(data=nymap, location="topleft", symbol=12)
vsites3
ggsave(filename=paste0(figurepath, "3_validation_sites_nydata.png"))
## Saving 7 x 5 in image
f2 <- y8hrmax~xmaxtemp+xwdsp+xrh
set.seed(44)
tn <- 62
valids <- c(1, 5, 10)
validt <- sort(sample(1:62, size=31))
vrows <- getvalidrows(sn=28, tn=62, valids=valids, validt=validt)
ymat <- matrix(nysptime$y8hrmax, byrow=T, ncol=tn)
yholdout <- ymat[valids, validt]
M31 <- Bsptime(package="spTimer",formula=f2, data=nysptime,
coordtype="utm", coords=4:5,
validrows=vrows, model="GP",
mchoice=F, n.report =1, scale.transform = "NONE", N=5000)
##
## Output: GP models
## ---------------------------------------------------------------
## Sampled: 5000 of 5000, 100.00%.
## Batch Acceptance Rate (phi): 37.71%
## Checking Parameters:
## phi: 0.0139, sig2eps: 9.9769, sig2eta: 107.8830
## beta[1]: -10.1861 beta[2]: 2.0792 beta[3]: 1.8684 beta[4]: -1.8902
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 , Iterations = 5000 .
## # Overall Acceptance Rate (phi) = 37.7 %
## ##
## ##
## # Elapsed time: 6.07 Sec.
## ##
##
## # Model: GP
## ##
## # Total time taken:: 11.04 - Sec.
summary(M31)
##
## The GP model has been fitted using the spTimer package.
## Call:
## Bsptime(formula = f2, data = nysptime, package = "spTimer", model = "GP",
## coordtype = "utm", coords = 4:5, validrows = vrows, scale.transform = "NONE",
## N = 5000, n.report = 1, mchoice = F)
## ##
## # Total time taken:: 11.04 - Sec.
## Model formula
## y8hrmax ~ xmaxtemp + xwdsp + xrh
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) -17.146 6.379 -29.858 -4.307
## xmaxtemp 2.349 0.176 1.992 2.698
## xwdsp 1.545 0.294 0.964 2.112
## xrh -1.557 0.751 -3.010 -0.035
## sig2eps 8.623 1.390 6.178 11.511
## sig2eta 112.043 13.713 92.320 144.628
## phi 0.013 0.002 0.010 0.016
##
## Validation Statistics:
## rmse mae crps cvg
## 5.510 4.262 5.339 100.000
modfit <- M31$fit
## Extract the fits for the validation sites
fitall <- data.frame(modfit$fitted)
head(fitall)
## Mean SD
## 1 61.23963 10.103878
## 2 58.96459 9.862792
## 3 59.68148 9.992403
## 4 44.86823 10.154895
## 5 37.93188 10.459307
## 6 30.52798 10.290717
fitall$s.index <- rep(1:28, each=tn)
library(spTimer)
##
## ## spTimer version: 3.3.1
vdat <- spT.subset(data=nysptime, var.name=c("s.index"), s=valids)
fitvalid <- spT.subset(data=fitall, var.name=c("s.index"), s=valids)
head(fitvalid)
## s.index Mean SD
## 1 1 61.23963 10.103878
## 2 1 58.96459 9.862792
## 3 1 59.68148 9.992403
## 4 1 44.86823 10.154895
## 5 1 37.93188 10.459307
## 6 1 30.52798 10.290717
fitvalid$low <- fitvalid$Mean - 1.96 * fitvalid$SD
fitvalid$up <- fitvalid$Mean + 1.96 * fitvalid$SD
# fitvalid$yobs <- sqrt(vdat$y8hrmax)
fitvalid$yobs <- vdat$y8hrmax
yobs <- matrix(fitvalid$yobs, byrow=T, ncol=tn)
y.valids.low <- matrix(fitvalid$low, byrow=T, ncol=tn)
y.valids.med <- matrix(fitvalid$Mean, byrow=T, ncol=tn)
y.valids.up <- matrix(fitvalid$up, byrow=T, ncol=tn)
p1 <- fig11.13.plot(yobs[1, ], y.valids.low[1, ], y.valids.med[1, ], y.valids.up[1, ], misst=validt)
p1 <- p1 + ggtitle("Validation for Site 1")
p1
p2 <- fig11.13.plot(yobs[2, ], y.valids.low[2, ], y.valids.med[2, ], y.valids.up[2, ], misst=validt)
p2 <- p2 + ggtitle("Validation for Site 5")
p2
p3 <- fig11.13.plot(yobs[3, ], y.valids.low[3, ], y.valids.med[3, ], y.valids.up[3, ], misst=validt)
p3 <- p3 + ggtitle("Validation for Site 10")
p3
ggarrange(p1, p2, p3, common.legend = TRUE, legend = "top", nrow = 3, ncol = 1)
ggsave(filename=paste0(figurepath, "validation_plot_BCG_spTimer.png"), device = "png")
## Saving 7 x 5 in image
sitemeans <- function(a, sn, tn=62) {
u <- matrix(a, nrow=sn, ncol=tn, byrow=T)
b <- apply(u, 1, mean)
as.vector(b)
}
post <- M3$fit
library(spTimer)
gpred <- predict(post, newdata=gridnysptime, newcoords=~Longitude+Latitude)
u <- gpred$pred.samples
v <- apply(u, 2, sitemeans, sn=100)
a <- get_parameter_estimates(t(v))
b <- data.frame(gridnyspatial[, 1:5], a)
meanmat <- post$op
sig2eps <- post$sig2ep
sige <- sqrt(sig2eps)
itmax <- ncol(meanmat)
nT <- nrow(nysptime)
sigemat <- matrix(rep(sige, each=nT), byrow=F, ncol=itmax)
a <- matrix(rnorm(nT*itmax), nrow=nT, ncol=itmax)
ypreds <- meanmat + a * sigemat
ypreds <- (ypreds)^2
v <- apply(ypreds, 2, sitemeans, sn=28)
a <- get_parameter_estimates(t(v))
fits <- data.frame(nyspatial[, 1:5], a)
b <- rbind(b, fits)
coord <- nyspatial[, c("Longitude","Latitude")]
library(akima)
xo <- seq(from=min(coord$Longitude)-0.5, to = max(coord$Longitude)+0.8, length=200)
yo <- seq(from=min(coord$Latitude)-0.25, to = max(coord$Latitude)+0.8, length=200)
surf <- interp(b$Longitude, b$Latitude, b$mean, xo=xo, yo=yo)
v <- fnc.delete.map.XYZ(xyz=surf)
interp1 <- data.frame(long = v$x, v$z )
names(interp1)[1:length(v$y)+1] <- v$y
library(tidyr)
interp1 <- gather(interp1,key = lat,value =Predicted,-long,convert = TRUE)
library(ggplot2)
nymap <- map_data(database="state",regions="new york")
zr <- range(interp1$Predicted, na.rm=T)
P <- ggplot() +
geom_raster(data=interp1, aes(x = long, y = lat,fill = Predicted)) +
geom_polygon(data=nymap, aes(x=long, y=lat, group=group), color="black", size = 0.6, fill=NA) +
geom_point(data=coord, aes(x=Longitude,y=Latitude)) +
stat_contour(data=na.omit(interp1), aes(x = long, y = lat,z = Predicted), colour = "black", binwidth =2) +
scale_fill_gradientn(colours=colpalette, na.value="gray95", limits=zr) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(data =interp1, dist = 100, location = "bottomleft", transform=T, dist_unit = "km", st.dist = .05, st.size = 5, height = .06, st.bottom=T, model="WGS84") +
ggsn::north(data=interp1, location="topleft", symbol=12) +
labs(x="Longitude", y = "Latitude", size=2.5)
P
surf <- interp(b$Longitude, b$Latitude, b$sd, xo=xo, yo=yo)
v <- fnc.delete.map.XYZ(xyz=surf)
interp1 <- data.frame(long = v$x, v$z )
names(interp1)[1:length(v$y)+1] <- v$y
interp1 <- gather(interp1,key = lat,value =sd,-long,convert = TRUE)
nymap <- map_data(database="state",regions="new york")
zr <- range(interp1$sd, na.rm=T)
Psd <- ggplot() +
geom_raster(data=interp1, aes(x = long, y = lat,fill = sd)) +
geom_polygon(data=nymap, aes(x=long, y=lat, group=group), color="black", size = 0.6, fill=NA) +
geom_point(data=coord, aes(x=Longitude,y=Latitude)) +
stat_contour(data=na.omit(interp1), aes(x = long, y = lat,z = sd), colour = "black", binwidth =0.1) +
scale_fill_gradientn(colours=colpalette, na.value="gray95", limits=zr) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggsn::scalebar(data =interp1, dist = 100, location = "bottomleft", transform=T, dist_unit = "km",
st.dist = .05, st.size = 5, height = .06, st.bottom=T, model="WGS84") +
ggsn::north(data=interp1, location="topleft", symbol=12) +
labs(x="Longitude", y = "Latitude", size=2.5)
Psd
if (longrun) {
M1.c <- Bsptime(model="lm", formula=f2, data=nysptime,
scale.transform = "SQRT", mchoice=T)
M2.c <- Bsptime(model="separable", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, phi.s=0.005, phi.t=0.05,
scale.transform = "SQRT", mchoice=T)
M3.c <- Bsptime(package="spTimer", model="GP",
formula=f2, data=nysptime, coordtype="utm",
coords=4:5, scale.transform = "SQRT",
mchoice=T, N=5000, n.report=1)
M4.c <- Bsptime(package="stan",formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
N=1500, burn.in=500, mchoice=T, verbose = F)
table7.3 <- M4.c$params
dput(table7.3, file=paste0(filepath, "table7.3.txt"))
table7.4 <- cbind(M1.c$mchoice, M2.c$mchoice, M3.c$mchoice, M4.c$mchoice)
dput(table7.4, file=paste0(filepath, "table7.4.txt"))
aresid <- residuals(M4.c)
ggsave(filename=paste0(figurepath, "M4_residuals.png"))
} else {
table7.4 <- dget(file=paste0(filepath, "table7.4.txt"))
}
round(table7.4, 2)
## [,1] [,2] [,3] [,4]
## pdic 5.06 5.07 78.65 30.70
## pdicalt 5.38 5.38 841.96 31.37
## dicorig 3912.25 3214.02 3132.10 2695.77
## dicalt 3912.88 3214.64 4658.72 2697.11
## pwaic1 4.95 14.11 48.53 9.32
## pwaic2 4.97 14.18 132.90 10.31
## waic1 3912.14 2449.15 2603.86 2088.81
## waic2 3912.18 2449.30 2772.60 2090.79
## gof 963.68 285.75 216.75 328.83
## penalty 967.10 240.63 873.84 361.70
## pmcc 1930.78 526.38 1090.59 690.54
# xtable(table7.4, digits=2)
if (longrun) {
valids <- c(8,11,12,14,18,21,24,28)
vrows <- getvalidrows(sn=28, tn=62, valids=valids, allt=T)
M1.v <- Bsptime(model="lm", formula=f2, data=nysptime,
scale.transform = "SQRT", validrows=vrows)
M2.v <- Bsptime(model="separable", formula=f2, data=nysptime,
coordtype="utm", coords=4:5,
phi.s=0.005, phi.t=0.05, scale.transform = "SQRT",
validrows=vrows)
M3.v <- Bsptime(package="spTimer", model="GP", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
validrows=vrows, N=5000, n.report=1)
M4.v <- Bsptime(package="stan",formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
N=1500, burn.in=500, validrows=vrows,
verbose = F)
psums <- get_validation_summaries(t(M4.v$valpreds))
a <- obs_v_pred_plot(yobs=M4.v$yobs_preds$y8hrmax, predsums=psums, segments = FALSE)
ggsave(filename = paste0(figurepath, "obs_v_pred_stan.png"))
table7.5 <- cbind(unlist(M1.v$stats), unlist(M2.v$stats), unlist(M3.v$stats),
unlist(M4.v$stats))
dput(table7.5, file=paste0(filepath, "table7.5.txt"))
#xtable(table7.5, digits=2)
} else {
table7.5 <- dget(file=paste0(filepath, "table7.5.txt"))
}
table7.5
## [,1] [,2] [,3] [,4]
## rmse 9.367171 6.497137 6.401791 6.405868
## mae 7.549492 5.0126 4.938585 4.836409
## crps 5.627457 10.66969 6.791781 3.334424
## cvg 97.95082 99.59016 99.59016 93.03279
M5 <- Bsptime(package="spTimer", model="AR", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
mchoice=T, N=5000, n.report=1)
##
## Output: AR models
## ---------------------------------------------------------------
## Sampled: 5000 of 5000, 100.00%.
## Batch Acceptance Rate (phi): 27.93%
## Checking Parameters:
## phi: 0.0079, rho: 0.5568, sig2eps: 0.0160, sig2eta: 0.6058
## beta[1]: 2.3222 beta[2]: 0.0671 beta[3]: 0.0327 beta[4]: -0.3421
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 , Iterations = 5000 .
## # Overall Acceptance Rate (phi) = 27.92 %
## ##
## ##
## # Elapsed time: 6.51 Sec.
## ##
##
## # Model: AR
## ##
## # Total time taken:: 6.52 - Sec.
a <- residuals(M5)
##
## Note that the residuals are provided on the transformed scale. Please see the scale.transform argument.
##
## Summary of the residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1.64400 -0.25597 0.01703 0.01725 0.27419 1.76014 24
ggsave(filename=paste0(figurepath, "M5_residuals.png"))
## Saving 7 x 5 in image
summary(M5)
##
## The AR model has been fitted using the spTimer package.
## Call:
## Bsptime(formula = f2, data = nysptime, package = "spTimer", model = "AR",
## coordtype = "utm", coords = 4:5, scale.transform = "SQRT",
## N = 5000, n.report = 1, mchoice = T)
## ##
## # Total time taken:: 6.52 - Sec.
## Model formula
## y8hrmax ~ xmaxtemp + xwdsp + xrh
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) 1.437 0.540 0.383 2.513
## xmaxtemp 0.091 0.015 0.060 0.122
## xwdsp 0.031 0.023 -0.014 0.078
## xrh -0.191 0.060 -0.304 -0.074
## rho 0.512 0.021 0.471 0.554
## sig2eps 0.014 0.002 0.010 0.019
## sig2eta 0.570 0.032 0.511 0.638
## phi 0.009 0.001 0.008 0.010
##
## Model Choice Statistics:
## gof penalty pmcc
## 321.33 607.31 928.64
table7.6 <- M5$params
# xtable(table7.6, digits=3)
dput(table7.6, file=paste0(filepath, "table7.6.txt"))
if (longrun) {
M6 <- Bsptime(package="inla", model="AR", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
mchoice=T)
# 5 minutes
summary(M6)
table7.7 <- M6$params
dput(table7.7, file=paste0(filepath, "table7.7.txt"))
} else {
table7.7 <- dget(file=paste0(filepath, "table7.7.txt"))
}
round(table7.7, 3)
## mean sd 2.5% 97.5%
## xmaxtemp 0.241 0.003 0.235 0.247
## xwdsp 0.052 0.012 0.028 0.077
## xrh -0.044 0.018 -0.077 -0.009
## rho 0.889 0.049 0.771 0.960
## sigma2eps 0.442 0.016 0.412 0.474
## sig2eta 0.764 0.518 0.197 2.086
## phi 0.348 0.144 0.155 0.715
# xtable(table7.7, digits=3)
if (longrun) {
f2 <- y8hrmax~xmaxtemp+xwdsp+xrh
valids <- c(8,11,12,14,18,21,24,28)
vrows <- getvalidrows(sn=28, tn=62, valids=valids, allt=T)
M5.v <- Bsptime(package="spTimer", model="AR", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
validrows=vrows, mchoice=T, N=5000, n.report = 1)
M6.v <- Bsptime(package="inla", model="AR", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
validrows=vrows, mchoice=T)
a1 <- c(M5$mchoice, unlist(M5.v$stats))
b <- as.list(M6$mchoice)
b
b1 <- c(b$gof, b$penalty, b$pmcc, unlist(M6.v$stats))
b1 <- c(M6$mchoice$gof, M6$mchoice$penalty, M6$mchoice$pmcc, unlist(M6.v$stats))
table7.8 <- rbind(a1, b1)
dput(table7.8, file=paste0(filepath, "table7.8.txt"))
} else {
table7.8 <- dget(file=paste0(filepath, "table7.8.txt"))
}
table7.8
## gof penalty pmcc rmse mae crps cvg
## a1 321.3300 607.31000 928.6400 6.461426 4.991846 5.969121 99.38525
## b1 736.7078 21.35412 758.0619 9.719906 7.641264 2.644965 65.16393
# xtable(table7.8, digits=2)
library(spTDyn)
## Registered S3 method overwritten by 'spTDyn':
## method from
## print.spTpred spTimer
## ## spTDyn version: 2.0.1
f3 <- y8hrmax~ xmaxtemp + sp(xmaxtemp)+ tp(xwdsp) + xrh
M7 <- Bsptime(package="sptDyn", model="GP", formula=f3, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
mchoice=T, N=5000, n.report=1)
##
## Output: GP models
## ---------------------------------------------------------------
## Sampled: 5000 of 5000, 100.00%.
## Batch Acceptance Rate (phi): 38.07%
## Checking Parameters:
## phi: 0.0171, sig2eps: 0.0133, sig2eta: 0.2656, sig2beta: 0.0591,
## sig2delta: 0.0466, sig2op: 0.5043,
## rho[1]: 0.0164
## beta[1]: 2.4459 beta[2]: 0.1842 beta[3]: -0.2190
## ---------------------------------------------------------------
## ## Model used spatially and temporally varying dynamic parameters
## ## Spatial and dynamic beta parameters are omitted in the display
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 , Iterations = 5000 .
## # Overall Acceptance Rate (phi) = 38.06 %
## ##
## ##
## # Elapsed time: 22.26 Sec.
## ##
##
## # Model: GP
## ##
## # Total time taken:: 22.28 - Sec.
summary(M7)
##
## The GP model has been fitted using the spTDyn package.
## Call:
## Bsptime(formula = f3, data = nysptime, package = "sptDyn", model = "GP",
## coordtype = "utm", coords = 4:5, scale.transform = "SQRT",
## N = 5000, n.report = 1, mchoice = T)
## ##
## # Total time taken:: 22.28 - Sec.
## Model formula
## y8hrmax ~ xmaxtemp + sp(xmaxtemp) + tp(xwdsp) + xrh
##
##
## Parameter Estimates:
## mean sd 2.5% 97.5%
## (Intercept) 1.887 0.597 0.727 3.042
## xmaxtemp 0.192 0.015 0.163 0.221
## xrh -0.133 0.064 -0.258 -0.006
## rho1 0.263 0.231 -0.182 0.715
## sig2eps 0.015 0.003 0.011 0.021
## sig2eta 0.272 0.028 0.230 0.340
## sig2beta 0.067 0.018 0.040 0.110
## sig2delta 0.048 0.009 0.033 0.068
## sig2op 0.853 1.051 0.190 3.271
## phi 0.018 0.002 0.014 0.022
##
## Model Choice Statistics:
## gof penalty pmcc
## 166.93 434.90 601.83
table7.9 <- M7$params
round(table7.9, 3)
## mean sd 2.5% 97.5%
## (Intercept) 1.887 0.597 0.727 3.042
## xmaxtemp 0.192 0.015 0.163 0.221
## xrh -0.133 0.064 -0.258 -0.006
## rho1 0.263 0.231 -0.182 0.715
## sig2eps 0.015 0.003 0.011 0.021
## sig2eta 0.272 0.028 0.230 0.340
## sig2beta 0.067 0.018 0.040 0.110
## sig2delta 0.048 0.009 0.033 0.068
## sig2op 0.853 1.051 0.190 3.271
## phi 0.018 0.002 0.014 0.022
# xtable(table7.9, digits=3)
dput(table7.9, file=paste0(filepath, "table7.9.txt"))
out <- M7$fit
dim(out$betasp)
## [1] 28 4000
a <- out$betasp
u <- c(t(out$betasp))
sn <- nrow(a)
itmax <- ncol(a)
v <- rep(1:sn, each=itmax)
d <- data.frame(site=as.factor(v), sp = u)
p <- ggplot(data=d, aes(x=site, y=sp)) +
geom_boxplot(outlier.colour="black", outlier.shape=1,
outlier.size=0.5) +
geom_abline(intercept=0, slope=0, color="blue") +
labs(title= "Spatial effects of maximum temperature", x="Site", y = "Effects", size=2.5)
p
ggsave(filename=paste0(figurepath, "sptDyn_spat_effects.png"))
## Saving 7 x 5 in image
dim(out$betatp)
## [1] 62 4000
b <- out$betatp
dim(b)
## [1] 62 4000
tn <- nrow(b)
itmax <- ncol(b)
tids <- 1:tn
stat <- apply(b[tids,], 1, quantile, prob=c(0.025,0.5,0.975))
tstat <- data.frame(tids, t(stat))
dimnames(tstat)[[2]] <- c("Days", "low", "median", "up")
head(tstat)
## Days low median up
## 1 1 0.13888038 0.19322979 0.2488582012
## 2 2 0.04432397 0.07818191 0.1125185634
## 3 3 0.01682599 0.07758568 0.1360425933
## 4 4 0.05000672 0.11019927 0.1695704982
## 5 5 -0.12461666 -0.06321995 -0.0005297627
## 6 6 -0.14142194 -0.06967983 0.0032164089
yr <- c(min(c(stat)),max(c(stat)))
p1 <- ggplot(data=tstat, aes(x=Days, y=median)) +
geom_point(size=3) +
ylim(yr) +
geom_segment(data=tstat, aes(x=Days, y=median, xend=Days, yend=low), linetype=1) +
geom_segment(data=tstat, aes(x=Days, y=median, xend=Days, yend=up), linetype=1) +
geom_abline(intercept=0, slope=0, col="blue") +
labs(title="Temporal effects of wind speed", x="Days", y="Temporal effects")
p1
ggsave(filename=paste0(figurepath, "sptDyn_dynamic_effects.png"))
## Saving 7 x 5 in image
set.seed(44)
validt <- sort(sample(1:62, size=31))
valids <- c(1, 5, 10)
vrows <- getvalidrows(sn=28, tn=62, valids=valids, validt=validt)
ymat <- matrix(nysptime$y8hrmax, byrow=T, ncol=62)
yholdout <- ymat[valids, validt]
scale.transform <- "SQRT"
library(spTDyn)
M7.v <- Bsptime(package="sptDyn", model="GP", formula=f3, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
mchoice=T, N=5000, validrows=vrows, n.report=1)
M8 <- Bsptime(package="spBayes", formula=f2, data=nysptime,
prior.sigma2=c(2, 25),
prior.tau2 =c(2, 25),
prior.sigma.eta =c(2, 0.001),
coordtype="utm",
coords=4:5, scale.transform = "SQRT",
mchoice=T, N=5000, n.report=2)
M8.v <- Bsptime(package="spBayes", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
prior.sigma2=c(2, 25),
prior.tau2 =c(2, 25),
prior.sigma.eta =c(2, 0.001),
mchoice=T, N=5000, validrows=vrows, n.report=1)
# summary(M8.v)
a1 <- c(M7$mchoice, unlist(M7.v$stats))
M8$mchoice <- as.list(M8$mchoice)
b1 <- c(M8$mchoice$gof, M8$mchoice$penalty, M8$mchoice$pmcc, unlist(M8.v$stats))
table7.10 <- rbind(a1, b1)
rownames(table7.10) <- c("spTDyn", "spBayes")
table7.10
dput(table7.10, file=paste0(filepath, "table7.10.txt"))
# library(xtable)
# xtable(table7.10, digits=2)
modfit <- M8$fit
N <- 5000
burn.in <- 1000
tn <- 62
ysamples <- modfit$p.y.samples[,burn.in:N]
# if (scale.transform=="SQRT")
ysamples <- (ysamples)^2
quant95 <- function(x){
quantile(x, prob=c(0.025, 0.5, 0.975))
}
y.hat <- apply(ysamples, 1, quant95) ## 3 by 1736
dim(y.hat)
## [1] 3 1736
## Extract the ones for hold out
y.hat.med <- matrix(y.hat[1, ], ncol=tn)
y.hat.up <- matrix(y.hat[3,], ncol=tn)
y.hat.low <- matrix(y.hat[2,], ncol=tn)
y.valids.med <- y.hat.med[valids,]
y.valids.low <- y.hat.low[valids,]
y.valids.up <- y.hat.up[valids,]
dim(y.valids.up)
## [1] 3 62
yobs <- ymat[valids,]
## change here
par(mfrow=c(3,1))
p1 <- fig11.13.plot(yobs[1, ], y.valids.low[1, ], y.valids.med[1, ], y.valids.up[1, ], misst=validt)
p2 <- fig11.13.plot(yobs[2, ], y.valids.low[2, ], y.valids.med[2, ], y.valids.up[2, ], misst=validt)
p3 <- fig11.13.plot(yobs[3, ], y.valids.low[3, ], y.valids.med[3, ], y.valids.up[3, ], misst=validt)
# Not included
# p1
# p2
# p3
# Not included
## spBayes parameter plots
beta <- apply(modfit$p.beta.samples[burn.in:N,], 2, quant95)
theta <- apply(modfit$p.theta.samples[burn.in:N,], 2, quant95)
sigma.sq <- theta[,grep("sigma.sq", colnames(theta))]
tau.sq <- theta[,grep("tau.sq", colnames(theta))]
phi <- theta[,grep("phi", colnames(theta))]
adat <- data.frame(x=1:tn, low=sigma.sq[1, ], med=sigma.sq[2, ], up=sigma.sq[3, ])
head(adat)
## x low med up
## sigma.sq.t1 1 1.685153 2.957488 5.667349
## sigma.sq.t2 2 1.952880 3.406915 6.642714
## sigma.sq.t3 3 2.039829 3.674009 7.078543
## sigma.sq.t4 4 2.096181 3.770592 7.638405
## sigma.sq.t5 5 2.108906 3.898746 7.970373
## sigma.sq.t6 6 2.152113 3.897676 7.911681
psigma <- ggplot() +
geom_point(data = adat, aes(x =x, y = med, shape=19), shape=19, col="blue", size = 2) +
geom_ribbon(data = adat, aes(x =x, ymin =low, ymax = up), alpha = 0.2, color = "grey50") +
theme(legend.position ="none") +
labs(y = "sigma2", x = "Days")
psigma
adat <- data.frame(x=1:tn, low=tau.sq[1, ], med=tau.sq[2, ], up=tau.sq[3, ])
head(adat)
## x low med up
## tau.sq.t1 1 1.466560 2.445941 4.404950
## tau.sq.t2 2 1.549077 2.616873 4.773666
## tau.sq.t3 3 1.604033 2.694501 5.088641
## tau.sq.t4 4 1.599751 2.707378 4.940954
## tau.sq.t5 5 1.676406 2.808222 5.170507
## tau.sq.t6 6 1.664785 2.803294 5.248475
ptau <- ggplot() +
geom_point(data = adat, aes(x =x, y = med, shape=19), shape=19, col="blue", size = 2) +
geom_ribbon(data = adat, aes(x =x, ymin =low, ymax = up), alpha = 0.2, color = "grey50") +
theme(legend.position ="none") +
labs(y = "tau2", x = "Days")
ptau
adat <- data.frame(x=1:tn, low=3/phi[3, ], med=3/phi[2, ], up=3/phi[1, ])
head(adat)
## x low med up
## phi.t1 1 239.6378 253.9526 277.0626
## phi.t2 2 241.7177 263.4341 281.0002
## phi.t3 3 240.2357 272.0259 311.4584
## phi.t4 4 263.0467 286.4816 333.5377
## phi.t5 5 289.1202 300.1259 313.8120
## phi.t6 6 279.2400 300.7139 332.7002
prange <- ggplot() +
geom_point(data = adat, aes(x =x, y = med, shape=19), shape=19, col="blue", size = 2) +
geom_ribbon(data = adat, aes(x =x, ymin =low, ymax = up), alpha = 0.2, color = "grey50") +
theme(legend.position ="none") +
labs(y = "Range", x = "Days")
prange
library(ggpubr)
ggarrange(psigma, ptau, prange, common.legend = TRUE, legend = "none", nrow = 3, ncol = 1)
ggsave(filename=paste0(figurepath, "theta_spBayes_dynamic.png"))
## Saving 7 x 5 in image
vnames <- all.vars(f2)
xnames <- vnames[-1]
k <- 4
cat("\nOnly first 4 beta parameters are plotted\n")
##
## Only first 4 beta parameters are plotted
beta.0 <- beta[,grep("Intercept", colnames(beta))]
adat <- data.frame(x=1:tn, low=beta.0[1,], med=beta.0[2,], up=beta.0[3,])
head(adat)
## x low med up
## (Intercept).t1 1 -17.71215 1.670596 19.90322
## (Intercept).t2 2 -22.03606 3.055726 27.54663
## (Intercept).t3 3 -20.24276 4.442070 29.51893
## (Intercept).t4 4 -20.85011 6.097570 32.58254
## (Intercept).t5 5 -20.79564 7.474298 34.15032
## (Intercept).t6 6 -21.49947 6.509345 35.06609
pint <- ggplot() +
geom_point(data = adat, aes(x =x, y = med, shape=19), shape=19, col="blue", size = 2) +
geom_ribbon(data = adat, aes(x =x, ymin =low, ymax = up), alpha = 0.2, color = "grey50") +
geom_hline(yintercept=0, linetype="dashed", color = "red", size=1)+
theme(legend.position ="none") +
labs(y = "Intercept", x = "Days")
pint
j <- 2
betaj <- beta[,grep(xnames[j-1], colnames(beta))]
adat <- data.frame(x=1:tn, low=betaj[1,], med=betaj[2,], up=betaj[3,])
ptmp <- ggplot() +
geom_point(data = adat, aes(x =x, y = med, shape=19), shape=19, col="blue", size = 2) +
geom_ribbon(data = adat, aes(x =x, ymin =low, ymax = up), alpha = 0.2, color = "grey50") +
geom_hline(yintercept=0, linetype="dashed", color = "red", size=1)+
theme(legend.position ="none") +
labs(y = "Max temp", x = "Days")
ptmp
j <- 3
betaj <- beta[,grep(xnames[j-1], colnames(beta))]
adat <- data.frame(x=1:tn, low=betaj[1,], med=betaj[2,], up=betaj[3,])
head(adat)
## x low med up
## xwdsp.1.t1 1 -1.4985475 0.56598079 2.653102
## xwdsp.1.t2 2 -1.2169513 0.40998021 1.991235
## xwdsp.1.t3 3 -1.6480721 0.54715825 2.509891
## xwdsp.1.t4 4 -0.8275222 0.51282102 1.730727
## xwdsp.1.t5 5 -2.0267017 -0.17835602 1.522054
## xwdsp.1.t6 6 -1.8557356 0.07192384 1.917112
pwdsp <- ggplot() +
geom_point(data = adat, aes(x =x, y = med, shape=19), shape=19, col="blue", size = 2) +
geom_ribbon(data = adat, aes(x =x, ymin =low, ymax = up), alpha = 0.2, color = "grey50") +
theme(legend.position ="none") +
geom_hline(yintercept=0, linetype="dashed", color = "red", size=1) +
labs(y = "Wind speed", x = "Days")
pwdsp
j <- 4
betaj <- beta[,grep(xnames[j-1], colnames(beta))]
adat <- data.frame(x=1:tn, low=betaj[1,], med=betaj[2,], up=betaj[3,])
head(adat)
## x low med up
## xrh.1.t1 1 -4.673486 0.08340790 4.710658
## xrh.1.t2 2 -5.191933 -0.59409925 4.403957
## xrh.1.t3 3 -3.916520 -0.57157413 2.968728
## xrh.1.t4 4 -3.898812 -0.97632013 1.888477
## xrh.1.t5 5 -4.040512 -0.27236432 3.192834
## xrh.1.t6 6 -3.805560 -0.07178007 3.576767
prh <- ggplot() +
geom_point(data = adat, aes(x =x, y = med, shape=19), shape=19, col="blue", size = 2) +
geom_ribbon(data = adat, aes(x =x, ymin =low, ymax = up), alpha = 0.2, color = "grey50") +
theme(legend.position ="none") +
geom_hline(yintercept=0, linetype="dashed", color = "red", size=1)+
labs(y = "Rel humidity", x = "Days")
prh
library(ggpubr)
ggarrange(pint, ptmp, pwdsp, prh, common.legend = TRUE, legend = "none", nrow = 2, ncol = 2)
ggsave(filename=paste0(figurepath, "beta_spBayes_dynamic.png"), device = "png")
## Saving 7 x 5 in image
M9 <- Bsptime(package="spTimer", model="GPP", g_size=5,
formula=f2, data=nysptime, n.report=2,
coordtype="utm", coords=4:5, scale.transform = "SQRT")
##
## Output: GPP approximation models
## ---------------------------------------------------------------
## Sampled: 1000 of 2000, 50.00%.
## Batch Acceptance Rate (phi): 35.94%
## Checking Parameters:
## phi: 0.0071, rho: 0.1538, sig2eps: 0.1473, sig2eta: 0.8049
## beta[1]: 4.1054 beta[2]: 0.1055 beta[3]: 0.1185 beta[4]: -0.1933
## ---------------------------------------------------------------
## ---------------------------------------------------------------
## Sampled: 2000 of 2000, 100.00%.
## Batch Acceptance Rate (phi): 35.97%
## Checking Parameters:
## phi: 0.0078, rho: 0.2477, sig2eps: 0.1344, sig2eta: 0.7907
## beta[1]: 2.2545 beta[2]: 0.1592 beta[3]: 0.0656 beta[4]: -0.0172
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 . Iterations = 2000 .
## # Acceptance rate: (phi) = 35.95 %
## ##
## ##
## # Elapsed time: 1.57 Sec.
## ##
##
## # Model: GPP
## ##
## # Total time taken:: 1.58 - Sec.
table7.12 <- M9$params
dput(table7.12, file=paste0(filepath, "table7.12.txt"))
aresid <- residuals(M9)
##
## Note that the residuals are provided on the transformed scale. Please see the scale.transform argument.
##
## Summary of the residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -1.759101 -0.177575 0.015745 0.000779 0.185278 1.537252 24
ggsave(filename=paste0(figurepath, "M9_residuals.png"))
## Saving 7 x 5 in image
valids <- c(8,11,12,14,18,21,24,28)
vrows <- getvalidrows(sn=28, tn=62, valids=valids, allt=T)
M9.5 <- Bsptime(package="spTimer", model="GPP", g_size=5,
formula=f2, data=nysptime, validrow=vrows, n.report=2,
coordtype="utm", coords=4:5, scale.transform = "SQRT", mchoice=T)
##
## Output: GPP approximation models
## ---------------------------------------------------------------
## Sampled: 1000 of 2000, 50.00%.
## Batch Acceptance Rate (phi): 29.43%
## Checking Parameters:
## phi: 0.0050, rho: 0.1443, sig2eps: 0.1628, sig2eta: 0.6021
## beta[1]: 3.8814 beta[2]: 0.0991 beta[3]: 0.0797 beta[4]: -0.0220
## ---------------------------------------------------------------
## ---------------------------------------------------------------
## Sampled: 2000 of 2000, 100.00%.
## Batch Acceptance Rate (phi): 30.27%
## Checking Parameters:
## phi: 0.0050, rho: 0.1114, sig2eps: 0.1405, sig2eta: 0.7446
## beta[1]: 2.0783 beta[2]: 0.1443 beta[3]: 0.0865 beta[4]: 0.1397
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 . Iterations = 2000 .
## # Acceptance rate: (phi) = 30.25 %
## ##
## ##
## # Elapsed time: 1.7 Sec.
## ##
##
## # Model: GPP
## ##
## # Total time taken:: 3.95 - Sec.
M9.4 <- Bsptime(package="spTimer", model="GPP", g_size=4, n.report=2,
formula=f2, data=nysptime, validrow=vrows,
coordtype="utm", coords=4:5, scale.transform = "SQRT", mchoice=T)
##
## Output: GPP approximation models
## ---------------------------------------------------------------
## Sampled: 1000 of 2000, 50.00%.
## Batch Acceptance Rate (phi): 45.35%
## Checking Parameters:
## phi: 0.0062, rho: 0.3409, sig2eps: 0.1569, sig2eta: 1.0597
## beta[1]: 3.0971 beta[2]: 0.1351 beta[3]: 0.0258 beta[4]: 0.1065
## ---------------------------------------------------------------
## ---------------------------------------------------------------
## Sampled: 2000 of 2000, 100.00%.
## Batch Acceptance Rate (phi): 43.17%
## Checking Parameters:
## phi: 0.0071, rho: 0.3867, sig2eps: 0.1526, sig2eta: 0.8153
## beta[1]: 4.7592 beta[2]: 0.0741 beta[3]: 0.0645 beta[4]: 0.0247
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 . Iterations = 2000 .
## # Acceptance rate: (phi) = 43.15 %
## ##
## ##
## # Elapsed time: 0.89 Sec.
## ##
##
## # Model: GPP
## ##
## # Total time taken:: 3.03 - Sec.
M9.3 <- Bsptime(package="spTimer", model="GPP", g_size=3, n.report=2,
formula=f2, data=nysptime, validrow=vrows,
coordtype="utm", coords=4:5, scale.transform = "SQRT", mchoice=T)
##
## Output: GPP approximation models
## ---------------------------------------------------------------
## Sampled: 1000 of 2000, 50.00%.
## Batch Acceptance Rate (phi): 51.55%
## Checking Parameters:
## phi: 0.0034, rho: 0.7170, sig2eps: 0.1849, sig2eta: 1.6803
## beta[1]: 6.3981 beta[2]: -0.0290 beta[3]: -0.0188 beta[4]: -0.1450
## ---------------------------------------------------------------
## ---------------------------------------------------------------
## Sampled: 2000 of 2000, 100.00%.
## Batch Acceptance Rate (phi): 53.03%
## Checking Parameters:
## phi: 0.0047, rho: 0.3533, sig2eps: 0.1954, sig2eta: 0.9971
## beta[1]: 6.7283 beta[2]: -0.0019 beta[3]: 0.0442 beta[4]: -0.1023
## ---------------------------------------------------------------
## ##
## # nBurn = 1000 . Iterations = 2000 .
## # Acceptance rate: (phi) = 53 %
## ##
## ##
## # Elapsed time: 0.48 Sec.
## ##
##
## # Model: GPP
## ##
## # Total time taken:: 2.66 - Sec.
u3 <- c(3, M9.3$mchoice, unlist(M9.3$stats))
u4 <- c(4, M9.4$mchoice, unlist(M9.4$stats))
u5 <- c(5, M9.5$mchoice, unlist(M9.5$stats))
table7.11 <- rbind(u3, u4, u5)
table7.11
## gof penalty pmcc rmse mae crps cvg
## u3 3 200.13 953.99 1154.12 6.573376 5.047318 7.802914 99.79508
## u4 4 145.39 836.05 981.44 6.336695 4.820222 7.024905 99.38525
## u5 5 146.31 815.13 961.44 6.404986 4.874754 7.300427 99.38525
# xtable(table7.11, digits=2)
dput(table7.11, file=paste0(filepath, "table7.11.txt"))
coords <- nyspatial[, c("Longitude", "Latitude")]
nymap <- map_data(database="state",regions="new york")
head(nymap)
head(kcoords)
colnames(kcoords) <- c("Longitude", "Latitude")
knotplot <- ggplot() +
geom_polygon(data=nymap, aes(x=long, y=lat, group=group),
color="black", size = 0.6, fill=NA) +
geom_point(data =coords, aes(x=Longitude,y=Latitude), shape=8) +
geom_point(data =kcoords, aes(x=Longitude,y=Latitude), fill="red", col="red", size=2, shape=21) +
ggsn::scalebar(data =nymap, dist = 100, location = "topleft", transform=T, dist_unit = "km",
st.dist = .05, st.size = 5, height = .06, st.bottom=T, model="WGS84") +
ggsn::north(data=nymap, location="topright", symbol=12) +
labs(x="Longitude", y = "Latitude", size=2.5)
knotplot
ggsave(file=paste0(figurepath, "knot_sites.png"))
if (longrun) {
f2 <- y8hrmax~xmaxtemp+xwdsp+xrh
f3 <- y8hrmax~ xmaxtemp + sp(xmaxtemp)+ tp(xwdsp) + xrh
valids <- c(8,11,12,14,18,21,24,28)
vrows <- getvalidrows(sn=28, tn=62, valids=valids, allt=T)
M1.v <- Bsptime(model="lm", formula=f2, data=nysptime,
scale.transform = "SQRT", validrows=vrows, N=5000, mchoice=T)
M2.v <- Bsptime(model="separable", formula=f2, data=nysptime,
coordtype="utm", coords=4:5,
phi.s=0.005, phi.t=0.05, scale.transform = "SQRT",
validrows=vrows, N=5000, mchoice=T)
M3.v <- Bsptime(package="spTimer", model="GP", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
mchoice=T, validrows=vrows, N=5000, n.report=2)
M4.v <- Bsptime(package="stan",formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
N=5000, mchoice=T, validrows=vrows,
verbose = F)
# 8 mins
M5.v <- Bsptime(package="spTimer", model="AR", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
validrows=vrows, mchoice=T, N=5000, n.report = 2)
M6.v <- Bsptime(package="inla", model="AR", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
validrows=vrows, N=5000, mchoice=T)
#6 mins 41s
library(spTDyn)
M7.v <- Bsptime(package="sptDyn", model="GP", formula=f3, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
mchoice=T, N=5000, validrows=vrows, n.report=2)
M8.v <- Bsptime(package="spBayes", formula=f2, data=nysptime,
coordtype="utm", coords=4:5, scale.transform = "SQRT",
prior.sigma2=c(2, 25),
prior.tau2 =c(2, 25),
prior.sigma.eta =c(2, 0.001),
mchoice=T, N=5000, validrows=vrows, n.report=2)
M9.v <- Bsptime(package="spTimer", model="GPP", g_size=5,
formula=f2, data=nysptime, validrow=vrows,
coordtype="utm", coords=4:5, scale.transform = "SQRT", mchoice=T,
n.report=2, N=5000)
results <- cbind.data.frame(lm=unlist(M1.v$stats), separable=unlist(M2.v$stats),
spTimerGP=unlist(M3.v$stats), stan=unlist(M4.v$stats),inla=unlist(M6.v$stats),
spTimerAR=unlist(M5.v$stats), spTDyn=unlist(M7.v$stats),
spBayes=unlist(M8.v$stats), sptimerGPP=unlist(M9.v$stats))
mcres <- cbind.data.frame(lm=unlist(M1.v$mchoice)[9:11], separable=unlist(M2.v$mchoice)[9:11],
spTimerGP=unlist(M3.v$mchoice)[9:11], stan=unlist(M4.v$mchoice)[9:11],
inla=unlist(M6.v$mchoice)[5:7],
spTimerAR=unlist(M5.v$mchoice), spTDyn=unlist(M7.v$mchoice),
spBayes=unlist(M8.v$mchoice), sptimerGPP=unlist(M9.v$mchoice))
results
allres <- rbind.data.frame(results, mcres)
allres
table7.13 <- allres[, -8]
# xtable(table7.13, digits=2)
dput(table7.13, file=paste0(filepath, "table7.13.txt"))
} else {
table7.13 <- dget(file=paste0(filepath, "table7.13.txt"))
}
table7.13
## lm separable spTimerGP stan inla spTimerAR
## rmse 9.349713 6.491904 6.401791 6.420053 9.728621 6.461426
## mae 7.543205 5.000449 4.938585 4.851079 7.649771 4.991846
## crps 5.665111 10.561421 6.791781 3.568710 2.669567 5.969121
## cvg 98.360656 99.590164 99.590164 92.827869 65.163934 99.385246
## gof 728.907405 218.488548 181.710000 173.447236 527.418430 185.760000
## penalty 731.614443 195.369254 935.420000 266.226947 17.291486 718.470000
## pmcc 1460.521848 413.857802 1117.130000 439.674183 544.709917 904.230000
## spTDyn sptimerGPP
## rmse 6.588660 6.357974
## mae 5.114160 4.850027
## crps 5.121886 7.472485
## cvg 99.385246 99.385246
## gof 71.300000 146.690000
## penalty 467.460000 815.850000
## pmcc 538.760000 962.540000
# colnames(table7.13) <- c("M1", "M2", "M3", "M4", "M5", "M6", "M7", "M9")
# xtable(table7.13, digits=2)
# All done
end.time <- proc.time()[3]
comp.time<- (end.time-start.time)/60
# comp.time<-fancy.time(comp.time)
print(comp.time)
## elapsed
## 3.954617