glmmfields
### https://cran.r-project.org/web/packages/glmmfields/vignettes/spatial-glms.html
library(glmmfields)
library(ggplot2)
library(dplyr)
library(viridis)
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
head(d)
## time pt y lon lat station_id temperature
## 1 1 1 2.5622221 6.588776 8.142518 1 -0.6264538
## 2 1 2 1.1557119 1.850700 9.287772 2 0.1836433
## 3 1 3 0.8233918 9.543781 1.474810 3 -0.8356286
## 4 1 4 2.2892430 8.978485 7.498217 4 1.5952808
## 5 1 5 1.7609267 9.436971 9.756573 5 0.3295078
## 6 1 6 1.2091675 7.236908 9.747925 6 -0.8204684
ggplot(s$dat, aes(lon, lat, colour = y)) +
viridis::scale_colour_viridis() +
geom_point(size = 3)+theme_bw()
m_glm <- glm(y ~ temperature, data = d, family = Gamma(link = "log"))
m_glm
##
## Call: glm(formula = y ~ temperature, family = Gamma(link = "log"),
## data = d)
##
## Coefficients:
## (Intercept) temperature
## 0.5369 0.1967
##
## Degrees of Freedom: 199 Total (i.e. Null); 198 Residual
## Null Deviance: 23.27
## Residual Deviance: 16.72 AIC: 280.9
## 2.5 % 97.5 %
## (Intercept) 0.4964271 0.5780115
## temperature 0.1522553 0.2413277
d$m_glm_residuals <- residuals(m_glm)
ggplot(d, aes(lon, lat, colour = m_glm_residuals)) +
scale_color_gradient2() +
geom_point(size = 3)+theme_bw()
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()
)
##
## SAMPLING FOR MODEL 'glmmfields' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.002 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 20 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 500 [ 0%] (Warmup)
## Chain 1: Iteration: 50 / 500 [ 10%] (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%] (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%] (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%] (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%] (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%] (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%] (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%] (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%] (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%] (Sampling)
## Chain 1: Iteration: 500 / 500 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 3.673 seconds (Warm-up)
## Chain 1: 2.545 seconds (Sampling)
## Chain 1: 6.218 seconds (Total)
## Chain 1:
## Inference for Stan model: glmmfields.
## 1 chains, each with iter=500; warmup=250; thin=1;
## post-warmup draws per chain=250, total post-warmup draws=250.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## gp_sigma 0.26 0.00 0.07 0.16 0.21 0.25 0.30 0.41 219 1.00
## gp_theta 1.65 0.01 0.21 1.28 1.51 1.63 1.78 2.12 221 1.00
## B[1] 0.46 0.01 0.07 0.33 0.42 0.46 0.50 0.60 33 1.01
## B[2] 0.19 0.00 0.02 0.16 0.18 0.19 0.20 0.22 164 1.02
## CV[1] 0.21 0.00 0.01 0.19 0.20 0.21 0.22 0.23 369 1.00
## lp__ -62.94 0.24 2.91 -68.98 -64.78 -62.46 -60.65 -58.23 145 1.00
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 03 21:30:01 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
plot(m_spatial, type = "spatial-residual", link = TRUE) +
geom_point(size = 3)+theme_bw()
plot(m_spatial, type = "residual-vs-fitted")+theme_bw()
plot(m_spatial, type = "prediction", link = FALSE) +
viridis::scale_colour_viridis() +
geom_point(size = 3)+theme_bw()
p <- predict(m_spatial)
head(p)
## # A tibble: 6 x 3
## estimate conf_low conf_high
## <dbl> <dbl> <dbl>
## 1 0.724 0.642 0.830
## 2 0.559 0.463 0.674
## 3 0.188 0.0797 0.326
## 4 0.922 0.795 1.05
## 5 0.619 0.498 0.730
## 6 0.496 0.394 0.600
p <- predict(m_spatial, type = "response")
head(p)
## # A tibble: 6 x 3
## estimate conf_low conf_high
## <dbl> <dbl> <dbl>
## 1 2.06 1.90 2.29
## 2 1.75 1.59 1.96
## 3 1.21 1.08 1.38
## 4 2.52 2.22 2.85
## 5 1.86 1.65 2.07
## 6 1.64 1.48 1.82
p <- predict(m_spatial, type = "response", interval = "prediction")
head(p)
## # A tibble: 6 x 3
## estimate conf_low conf_high
## <dbl> <dbl> <dbl>
## 1 2.06 1.38 3.09
## 2 1.75 1.05 2.66
## 3 1.21 0.747 1.79
## 4 2.52 1.52 3.88
## 5 1.86 1.21 2.58
## 6 1.64 0.997 2.49
head(tidy(m_spatial, conf.int = TRUE, conf.method = "HPDinterval"))
## # A tibble: 6 x 5
## term estimate std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 gp_sigma 0.252 0.0700 0.160 0.407
## 2 gp_theta 1.63 0.212 1.18 2.00
## 3 B[1] 0.456 0.0665 0.330 0.605
## 4 B[2] 0.192 0.0166 0.158 0.222
## 5 CV[1] 0.206 0.0119 0.190 0.233
## 6 spatialEffectsKnots[1,1] 0.358 0.0786 0.193 0.498
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()+theme_bw()
