Nothing
## ----setup, include=FALSE, cache=FALSE----------------------------------------
library(knitr)
# set global chunk options
# opts_chunk$set(fig.path='figure/minimal-', fig.align='center', fig.show='hold')
# options(formatR.arrow=TRUE,width=90)
knitr::opts_chunk$set(dpi=100)
## ----eval=TRUE,message=FALSE,warning=FALSE------------------------------------
library("FRK") # for carrying out FRK
library("sp") # for defining points/polygons
library("dplyr") # for easy data manipulation
library("ggplot2") # for plotting
## ----eval=TRUE,message=FALSE,warning=FALSE------------------------------------
m <- 250 # Sample size
RNGversion("3.6.0"); set.seed(1) # Fix seed
zdf <- data.frame(x = runif(m), y= runif(m)) # Generate random locs
zdf$Y <- 3 + sin(7 * zdf$x) + cos(9 * zdf$y) # Latent process
zdf$z <- rpois(m, lambda = exp(zdf$Y)) # Simulate data
coordinates(zdf) = ~x+y # Turn into sp object
## ----eval=TRUE,message=FALSE,warning=FALSE,results='hide'---------------------
S <- FRK(f = z ~ 1, # Formula to FRK
list(zdf), # All datasets are supplied in list
nres = 2, # Low-rank model to reduce run-time
response = "poisson", # data model
link = "log", # link function
nonconvex_hull = FALSE) # convex hull
pred <- predict(S) # prediction stage
## ----eval=TRUE,message=FALSE,warning=FALSE,results='hide'---------------------
plot_list <- plot(S, pred$newdata)
plot_list <- c(plot_list, plot_spatial_or_ST(zdf, "z"))
## ----echo=FALSE, warning=FALSE,fig.align='center',fig.cap="(Left) Simulated Poisson spatial data. (Centre) Prediction of the mean process. (Right) Uncertainty quantification of predictions; specifically the width of the 90\\% posterior predictive interval.\\label{fig:example1}",fig.subcap=c("",""), fig.width = 10, fig.height = 4----
ggpubr::ggarrange(plot_list$z + labs(fill = "data"),
plot_list$p_mu + labs(fill = "pred."),
plot_list$interval90_mu + labs(fill = "pred.\nuncertainty"),
nrow = 1, legend = "top")
## ----eval=TRUE,message=FALSE,warning=FALSE------------------------------------
data("Am_data")
coordinates(Am_data) = ~ Easting + Northing # convert to sp object
GZ_df <- data.frame("Easting" = 219868.09, "Northing" = 285320.8)
## ----echo=FALSE, warning=FALSE, results='hide'--------------------------------
# centre, width, and height
makeRectangle <- function(centre, w, h) {
vertices <- rbind(c(centre[, 1] - w/2, centre[, 2] - h/2),
c(centre[, 1] - w/2, centre[, 2] + h/2),
c(centre[, 1] + w/2, centre[, 2] + h/2),
c(centre[, 1] + w/2, centre[, 2] - h/2),
c(centre[, 1] - w/2, centre[, 2] - h/2))
Polygon(vertices)
}
## Following Paul and Cressie (2011),
## we predict over a series of concentric square blocks centred at Ground Zero
## (GZ), as well as a series of concentric square blocks away from GZ.
n_schemes <- 1
n_block <- 5
construct_block_scheme <- function() {
ratio <- 1.03 # width to height ratio of the blocks
w <- seq(43, 250, length.out = n_block)
h <- w / ratio
## Treat GZ as the centre, and expand relative to GZ.
blocks <- list()
for(i in 1:n_block) {
blocks[[i]] <- makeRectangle(centre = GZ_df, w = w[i], h = h[i])
blocks[[i]] <- Polygons(list(blocks[[i]]), paste0("block", i))
}
if (n_schemes == 2) {
## Now shift away from GZ
centre <- GZ_df
centre[, 1] <- GZ_df[, 1] - 153
centre[, 2] <- GZ_df[, 2] + 125
for(i in (n_block + 1):(2 * n_block)) {
blocks[[i]] <- makeRectangle(centre = centre, w = w[i - n_block], h = h[i- n_block])
blocks[[i]] <- Polygons(list(blocks[[i]]), paste0("block", i))
}
}
## (set the plotting order from largest to smallest)
pred_polygons <- SpatialPolygons(blocks, (n_schemes * n_block):1)
coordnames(pred_polygons) <- c("Easting", "Northing")
pred_polygons$Scheme <- rep(as.character(n_schemes:1), each = length(pred_polygons)/n_schemes)
return(pred_polygons)
}
blocks <- construct_block_scheme()
## ----echo=FALSE, warning=FALSE, results='hide', fig.align='center',fig.cap="Americium soil data and blocking scheme.\\label{fig:Am_data}", fig.width = 13.6, fig.height = 4.5----
nasa_palette <- c("#03006d","#02008f","#0000b6","#0001ef","#0000f6","#0428f6","#0b53f7","#0f81f3",
"#18b1f5","#1ff0f7","#27fada","#3efaa3","#5dfc7b","#85fd4e","#aefc2a","#e9fc0d","#f6da0c","#f5a009",
"#f6780a","#f34a09","#f2210a","#f50008","#d90009","#a80109","#730005")
lab1 <- xlab(as.expression(bquote("Easting /" ~ 10^5 ~ "m")))
lab2 <- ylab(as.expression(bquote("Northing /" ~ 10^5 ~ "m")))
formatter <- function(x){
x/10^5
}
x_scale <- scale_x_continuous(breaks = 10^5 *c(2.197, 2.199, 2.201), labels = formatter)
y_scale <- scale_y_continuous(breaks = 10^5 * c(2.852, 2.854, 2.856), labels = formatter)
## Basic plot to reduce code repetition
p_basic <- ggplot(data = as.data.frame(Am_data),
aes(x = Easting, y = Northing)) +
lab1 + lab2 + x_scale + y_scale + theme_bw() + coord_fixed()
## Data on the original scale
p_data <- p_basic +
geom_point(aes(colour = Am), size = 1) +
geom_point(data = GZ_df, shape = 4, size = 5) +
scale_colour_gradientn(colours = nasa_palette,
labels = scales::scientific,
breaks = c(250000, 750000))
## Data on the log scale
p_data_log_scale <- p_basic +
geom_point(aes(colour = log(Am)), size = 1) +
geom_point(data = GZ_df, shape = 4, size = 5) +
scale_colour_gradientn(colours = nasa_palette,
name = "Log-Americium",
breaks = c(9, 11, 13))
## Blocking scheme
p_Scheme_1_2 <- p_basic +
geom_point(size = 0.3) +
geom_point(data = GZ_df, shape = 4, size = 5) +
geom_polygon(data = FRK::SpatialPolygonsDataFrame_to_df(blocks),
aes(group = id, colour = Scheme), alpha = 0) +
labs(colour = "Blocking Scheme")
ggpubr::ggarrange(p_data + theme(legend.text=element_text(angle = 20)) +
theme(text = element_text(size=17)),
p_data_log_scale + theme(text = element_text(size=17)),
p_Scheme_1_2 + theme(text = element_text(size=17)),
nrow = 1, align = "hv", legend = "top")
## ----eval=TRUE,message=FALSE,warning=FALSE------------------------------------
BAUs <- auto_BAUs(manifold = plane(),
type = "grid",
data = Am_data,
nonconvex_hull = FALSE)
## Add covariates to the BAUs
d_cutoff <- 30.48
d_BAU <- distR(coordinates(BAUs), GZ_df)
BAUs$x1 <- as.numeric(d_BAU < d_cutoff)
BAUs$x2 <- d_BAU * BAUs$x1
BAUs$x3 <- as.numeric(d_BAU >= d_cutoff)
BAUs$x4 <- d_BAU * (BAUs$x3)
## ----eval=TRUE,message=FALSE,warning=FALSE,results='hide'---------------------
BAUs$fs <- 1
Am_data$std <- 1
S <- FRK(f = Am ~ -1 + x1 + x2 + x3 + x4, data = list(Am_data),
response = "gaussian",
link = "log",
BAUs = BAUs,
nres = 2,
est_error = FALSE)
## ----eval=TRUE,message=FALSE,warning=FALSE,results='hide'---------------------
pred <- predict(S, type = c("link", "mean"))
plot_list <- plot(S, pred$newdata)
## ----echo=FALSE, message=FALSE,warning=FALSE,results='hide',fig.align='center',fig.cap="Americium point predictions.\\label{fig:Am_BAU_predictions}",fig.subcap=c("",""),fig.width = 8, fig.height = 5----
plot_list <- lapply(
plot_list,
function(gg) gg + lab1 + lab2 + x_scale + y_scale)
ggpubr::ggarrange(
plot_list$p_Y + labs(fill = "Y pred.") +
scale_fill_gradientn(colours = nasa_palette, labels = scales::scientific),
plot_list$RMSPE_Y + labs(fill = "RMSPE Y"),
plot_list$p_mu + labs(fill = "mu pred.") +
scale_fill_gradientn(colours = nasa_palette, labels = scales::scientific),
plot_list$RMSPE_mu + labs(fill = "RMSPE mu"),
align = "hv", nrow = 2, ncol =2)
## ----eval=TRUE,message=FALSE,warning=FALSE,results='hide'---------------------
pred <- predict(S, newdata = blocks)
## ----echo=FALSE, warning=FALSE,fig.align='center',fig.cap="Americium block level predictions.\\label{fig:Am_block_predictions_spatial}",fig.subcap=c("",""),fig.width = 9.5, fig.height = 3.7----
block_pred_df <- FRK::SpatialPolygonsDataFrame_to_df(pred$newdata)
## Change level order to reverse the order that the blocks are plotted.
block_pred_df$block <- factor(block_pred_df$id,
levels = paste0("block", (n_schemes * n_block):1))
p_block <- p_basic +
geom_point(size = 0.3) +
geom_polygon(data = block_pred_df, aes(fill = p_mu, group = block),
alpha = 1, colour = "black") +
labs(fill = "mu pred.") +
scale_fill_gradientn(colours = nasa_palette, labels = scales::scientific)
p_block_RMSPE <- p_basic +
geom_point(size = 0.3) +
geom_polygon(data = block_pred_df, aes(fill = RMSPE_mu, group = block),
alpha = 1, colour = "black") +
scale_fill_distiller(palette = "BrBG",
labels = scales::scientific) +
labs(fill = "RMSPE mu")
ggpubr::ggarrange(p_block, p_block_RMSPE, nrow = 1, align = "hv")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.