if (FALSE) {
library(sf)
library(vmsae)
# this function is time consuming for the first run
install_environment()
load_environment()
acs_data <- read_sf(system.file("example", "mo_county.shp", package = "vmsae"))
y <- readRDS(system.file("example", "y.Rds", package = "vmsae"))
y_sigma <- readRDS(system.file("example", "y_sigma.Rds", package = "vmsae"))
X <- readRDS(system.file("example", "X.Rds", package = "vmsae"))
W <- readRDS(system.file("example", "W.Rds", package = "vmsae"))
num_samples <- 1000 # set to larger values in practice, e.g. 10000.
model <- vgmsfh_numpyro(y, y_sigma, X, W,
GEOID = acs_data$GEOID,
model_name = "mo_county", save_dir = NULL,
num_samples = num_samples, num_warmup = num_samples)
y_hat_np <- model@yhat_samples
y_hat_mean_np <- apply(y_hat_np, c(2, 3), mean)
y_hat_lower_np <- apply(y_hat_np, c(2, 3), quantile, 0.025)
y_hat_upper_np <- apply(y_hat_np, c(2, 3), quantile, 0.975)
plot(model, shp = acs_data, type = "compare", var_idx = 2)
}
Run the code above in your browser using DataLab