hrr: Estimate a hierarchical related regression model

In chrishanretty/hrr: Hierarchical Related Regression for Categorical Outcomes

Estimate a hierarchical related regression model

Description

hrr returns a list containing a Stan fit, a tidy data frame of coefficients, and pre-calculated tables showing the popularity of different options amongst different groups defined by the categorical predictors in the model.

Usage

hrr(
  f,
  data,
  ps,
  aux,
  res,
  areavar,
  weightvar,
  testing = FALSE,
  adjust = FALSE,
  overdispersed = FALSE,
  threading = FALSE,
  probs = c(0.025, 0.5, 0.975),
  mrp_only = FALSE,
  stan_control = list(adapt_delta = 0.99, max_treedepth = 12),
  ...
)

Arguments

f	a two-part formula, with categorical predictors found in the post-stratification frame followed by a vertical bar and (continuous or categorical) predictors found in the auxiliary area-level data.
data	a data frame containing the outcome variable and the categorical predictors specified in the first part of formula f. The data frame must also contain the area identifier named in areavar.
ps	a data frame containing the categorical predictors specified in the first part of formula f. The data frame must also contain the area identifier named in areavar, and the post-stratification weights named in weightvar.
aux	a data frame containing the continuous or categorical predictors specified in the second part of formula f. The data frame must also contain the area identified named in areavar.
res	a data frame containing the aggregate-level outcomes for each area. The data frame must include the area identifier named in areavar, and variables names equivalent to the levels of the response variable found in data.
areavar	a character string giving the name of an area identifier found in the different input data frames (data, ps, aux, res)
weightvar	a character string giving the name of the variable in data frame ps which contains the post-stratification weights. This is not the name of any variable in data containing survey weights: the package does not support such weights.
testing	a logical value. If testing is equal to TRUE, then the return value is a list containing Stan code and a list of data to pass to Stan. If testing is equal to FALSE, the model is estimated. Setting testing to TRUE can be useful to prototype Stan code which is then manually tweaked.
adjust	a logical value. If adjust is equal to TRUE, the model for individual level outcomes includes a parameter not found in the aggregate model. This parameter can capture survey effects. If adjust is equal to FALSE (the default) the standard hierarchical related regression model is estimated.
overdispersed	a logical value. If overdispersed is equal to TRUE, binary or categorical counts are modelled using a Beta-Binomial or Dirichlet-Multinomial distribution rather than a binomial or multinomial distribution.
threading	an integer. If threading is equal to zero (the default), hrr estimates the model using RStan, with each chain operating on at most one core. If threading is a non-zero number, hrr estimates the model using cmdstanr, with within-chain multi-threading and threading threads per chain.
probs	a vector of probabilities. This parameter affects the summary statistics returned.
mrp_only	a logical value. If mrp_only is equal to TRUE, then only individual responses are modelled. If mrp_only is equal to FALSE (the default), aggregate outcomes are also modelled.
stan_control	a list of arguments passed to RStan or cmdstanr's control argument.
...	other arguments passed on RStan or (if threading > 0) to cmdstanr.

Value

A list with entries area_smry, grp_smry, and fit.

Examples

data("toydat")
data("toypsf")

f <- cat_y ~ k_1 + k_2 + k_3 | x_1 + x_2

aux <- unique(toydat[,c("area", "x_1", "x_2")])
res <- unique(toydat[,c("area", "red", "green", "blue")])

## Generate Stan code and data for later use
mod <- hrr(f, data = toydat, ps = toypsf, aux = aux,
    res = res, areavar = "area", weightvar = "count",
    testing = TRUE, adjust = FALSE, overdispersed = TRUE)

## Computationally intensive bit
## Not run: 
mod <- hrr(f, data = toydat, ps = toypsf, aux = aux,
    res = res, areavar = "area", weightvar = "count",
    testing = FALSE, adjust = FALSE, overdispersed = TRUE,
    iter = 320, chains = 4, cores = 4)

## End(Not run)

chrishanretty/hrr documentation built on Oct. 23, 2024, 9:23 p.m.