inst/doc/hbsaems-betalogitnorm-model.R

In hbsaems: Hierarchical Bayes Small Area Estimation Model using 'Stan'

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE
)

## -----------------------------------------------------------------------------
# library(hbsaems)
# data("data_betalogitnorm")
# data <- data_betalogitnorm
# 
# data("adjacency_matrix_car")
# M <- adjacency_matrix_car

## -----------------------------------------------------------------------------
# model_prior_pred <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    data = data,
#    sample_prior = "only",
#    prior = c(
#     set_prior("normal(-1, 0.7)", class = "Intercept"),
#     set_prior("normal(0, 0.5)", class = "b"),
#     set_prior("gamma(2, 0.1)", class = "phi")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model_prior_pred)

## -----------------------------------------------------------------------------
# model_prior_pred_phi <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    n = "n",
#    deff = "deff",
#    data = data,
#    sample_prior = "only",
#    prior = c(
#     set_prior("normal(-1, 0.7)", class = "Intercept"),
#     set_prior("normal(0, 0.5)", class = "b")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model_prior_pred_phi)

## -----------------------------------------------------------------------------
# result_hbpc <- hbpc(model_prior_pred, response_var="y")
# summary(result_hbpc)

## -----------------------------------------------------------------------------
# result_hbpc$prior_predictive_plot

## -----------------------------------------------------------------------------
# model <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    n = "n",
#    deff = "deff",
#    data = data,
#    sample_prior = "no",
#    prior = c(
#      set_prior("normal(-1, 0.7)", class = "Intercept"),
#      set_prior("normal(0, 0.5)", class = "b")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model)

## -----------------------------------------------------------------------------
# model_with_defined_re <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    n = "n",
#    deff = "deff",
#    group = "group",
#    data = data,
#    sample_prior = "no",
#    prior = c(
#      set_prior("normal(-1, 0.7)", class = "Intercept"),
#      set_prior("normal(0, 0.5)", class = "b")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model_with_defined_re)

## -----------------------------------------------------------------------------
# # Prepare Missing Data
# data_missing <- data
# data_missing$y[sample(1:30, 3)] <- NA  # 3 missing values in response

## -----------------------------------------------------------------------------
# model_deleted <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    n = "n",
#    deff = "deff",
#    group = "group",
#    data = data_missing,
#    handle_missing = "deleted",
#    sample_prior = "no",
#    prior = c(
#      set_prior("normal(-1, 0.7)", class = "Intercept"),
#      set_prior("normal(0, 0.5)", class = "b")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model_deleted)

## -----------------------------------------------------------------------------
# # Prepare missing data
# data_missing <- data
# data_missing$y[3:5] <- NA
# data_missing$x1[6:7] <- NA

## -----------------------------------------------------------------------------
# model_multiple <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    n = "n",
#    deff = "deff",
#    group = "group",
#    data = data_missing,
#    handle_missing = "multiple",
#    sample_prior = "no",
#    prior = c(
#      set_prior("normal(-1, 0.7)", class = "Intercept"),
#      set_prior("normal(0, 0.5)", class = "b")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model_multiple)

## -----------------------------------------------------------------------------
# data_missing <- data
# data_missing$x1[3:5] <- NA
# data_missing$x2[14:17] <- NA

## -----------------------------------------------------------------------------
# model_during_model <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    n = "n",
#    deff = "deff",
#    group = "group",
#    data = data_missing,
#    handle_missing = "model",
#    sample_prior = "no",
#    prior = c(
#      set_prior("normal(-1, 0.7)", class = "Intercept"),
#      set_prior("normal(0, 0.5)", class = "b")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model_during_model)

## -----------------------------------------------------------------------------
# model_spatial <- hbm_betalogitnorm(
#    response = "y",
#    predictors = c("x1", "x2", "x3"),
#    n = "n",
#    deff = "deff",
#    group = "group",
#    sre = "sre",
#    sre_type = "car",
#    car_type = "icar",
#    M = M,
#    data = data,
#    sample_prior = "no",
#    prior = c(
#      set_prior("normal(-1, 0.7)", class = "Intercept"),
#      set_prior("normal(0, 0.5)", class = "b")
#   )
# )

## -----------------------------------------------------------------------------
# summary(model_spatial)

## -----------------------------------------------------------------------------
# result_hbcc <- hbcc(model)
# summary(result_hbcc)

## -----------------------------------------------------------------------------
# model <- update_hbm(model, iter = 8000)

## -----------------------------------------------------------------------------
# summary(model)

## -----------------------------------------------------------------------------
# result_hbcc <- hbcc(model)
# summary(result_hbcc)

## -----------------------------------------------------------------------------
# result_hbcc$plots$trace

## -----------------------------------------------------------------------------
# result_hbcc$plots$dens

## -----------------------------------------------------------------------------
# result_hbcc$plots$acf

## -----------------------------------------------------------------------------
# result_hbcc$plots$nuts_energy

## -----------------------------------------------------------------------------
# result_hbcc$plots$rhat

## -----------------------------------------------------------------------------
# result_hbcc$plots$neff

## -----------------------------------------------------------------------------
# result_hbmc <- hbmc(
#       model = list(model, model_spatial),
#       comparison_metrics = c("loo", "waic", "bf"),
#       run_prior_sensitivity= TRUE,
#       sensitivity_vars = c ("b_x1")
#   )
# 
# summary(result_hbmc)

## -----------------------------------------------------------------------------
# result_hbmc$primary_model_diagnostics$pp_check_plot

## -----------------------------------------------------------------------------
# result_hbmc$primary_model_diagnostics$params_plot

## -----------------------------------------------------------------------------
# result_hbsae <- hbsae(model)
# summary(result_hbsae)