# Fit NNGP model on the meuse dataset and compare with SPDE in INLA library(sp) library(gstat) library(rstan) library(rdist) # Load the data data("meuse") coordinates(meuse) <- ~x+y proj4string(meuse) <- CRS("+init=epsg:28992") ######################################### # Fit the model in Stan ########### ######################################### coords <- coordinates(meuse) M = 15 # Number of Nearest Neighbors NN.matrix <- NNMatrix(coords = coords, n.neighbors = M, n.omp.threads = 2) str(NN.matrix) X <- as.matrix(tibble(Intercept = rep(1, nrow(meuse)), dist = meuse$dist)) meuse$log_zinc <- log(meuse$zinc) Y <- matrix(meuse$log_zinc) P <- 2 N = nrow(meuse) options(mc.cores = parallel::detectCores()) stan_data <- list(N = N, M = M, P = P, Y = Y[NN.matrix$ord], X = X[NN.matrix$ord, ], # sorted Y and X NN_ind = NN.matrix$NN_ind, NN_dist = NN.matrix$NN_dist, NN_distM = NN.matrix$NN_distM) parameters <- c("beta", "sigma", "tau", "phi", "w", "y_rep") model <- stan( file = "meuse_NNGP.stan", data = stan_data, pars = parameters, iter = 6000, chains = 1, thin = 1, seed = 123 ) rstan::traceplot(model) print(model) ####################### # Model validation #### ####################### # Has the model get rid of the spatial correlation? library(geoR) # Non spatial model m <- lm(meuse$log_zinc ~ meuse$dist) summary(m) coords_meuse <- coordinates(meuse) V2 <- variog(coords = coords_meuse, data = m$residuals) plot(V2) # Spatial model in Stan samples <- extract(model) fit <- samples$y_rep fitted <- colSums(fit) / 3000 df_fit <- tibble(obs = meuse$log_zinc, fit = fitted) df_fit$res <- df_fit$obs - df_fit$fit V3 <- variog(coords = coords_meuse, data = df_fit$res) plot(V3)