Nothing
## ----set-knitr-options, cache=FALSE, echo=FALSE-------------------------------
library("knitr")
opts_chunk$set(message = FALSE, fig.width = 5.5)
## ----message=FALSE, warning=FALSE---------------------------------------------
library(glmmfields)
library(ggplot2)
library(dplyr)
## ----echo=TRUE, eval=FALSE----------------------------------------------------
# options(mc.cores = parallel::detectCores())
## ----simulate-data------------------------------------------------------------
set.seed(1)
N <- 200 # number of data points
temperature <- rnorm(N, 0, 1) # simulated temperature data
X <- cbind(1, temperature) # our design matrix
s <- sim_glmmfields(
n_draws = 1, gp_theta = 1.5, n_data_points = N,
gp_sigma = 0.2, sd_obs = 0.2, n_knots = 12, obs_error = "gamma",
covariance = "squared-exponential", X = X,
B = c(0.5, 0.2) # B represents our intercept and slope
)
d <- s$dat
d$temperature <- temperature
ggplot(s$dat, aes(lon, lat, colour = y)) +
viridis::scale_colour_viridis() +
geom_point(size = 3)
## -----------------------------------------------------------------------------
m_glm <- glm(y ~ temperature, data = d, family = Gamma(link = "log"))
m_glm
confint(m_glm)
d$m_glm_residuals <- residuals(m_glm)
ggplot(d, aes(lon, lat, colour = m_glm_residuals)) +
scale_color_gradient2() +
geom_point(size = 3)
## ----results='hide'-----------------------------------------------------------
m_spatial <- glmmfields(y ~ temperature,
data = d, family = Gamma(link = "log"),
lat = "lat", lon = "lon", nknots = 12, iter = 500, chains = 1,
prior_intercept = student_t(3, 0, 10),
prior_beta = student_t(3, 0, 3),
prior_sigma = half_t(3, 0, 3),
prior_gp_theta = half_t(3, 0, 10),
prior_gp_sigma = half_t(3, 0, 3),
seed = 123 # passed to rstan::sampling()
)
## -----------------------------------------------------------------------------
m_spatial
## -----------------------------------------------------------------------------
plot(m_spatial, type = "spatial-residual", link = TRUE) +
geom_point(size = 3)
## -----------------------------------------------------------------------------
plot(m_spatial, type = "residual-vs-fitted")
## -----------------------------------------------------------------------------
plot(m_spatial, type = "prediction", link = FALSE) +
viridis::scale_colour_viridis() +
geom_point(size = 3)
## -----------------------------------------------------------------------------
# link scale:
p <- predict(m_spatial)
head(p)
# response scale:
p <- predict(m_spatial, type = "response")
head(p)
# prediction intervals on new observations (include observation error):
p <- predict(m_spatial, type = "response", interval = "prediction")
head(p)
## -----------------------------------------------------------------------------
head(tidy(m_spatial, conf.int = TRUE, conf.method = "HPDinterval"))
## -----------------------------------------------------------------------------
pred_grid <- expand.grid(
lat = seq(min(d$lat), max(d$lat), length.out = 30),
lon = seq(min(d$lon), max(d$lon), length.out = 30)
)
pred_grid$temperature <- mean(d$temperature)
pred_grid$prediction <- predict(
m_spatial,
newdata = pred_grid,
type = "response"
)$estimate
ggplot(pred_grid, aes(lon, lat, fill = prediction)) +
geom_raster() +
viridis::scale_fill_viridis()
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.