star_MCMC: MCMC Algorithm for STAR

In drkowal/rSTAR: MCMC and EM algorithms for Simultaneous Transformation and Rounding (STAR) Models

MCMC Algorithm for STAR

Description

Run the MCMC algorithm for STAR given

1. a function to initialize model parameters; and

2. a function to sample (i.e., update) model parameters.

The transformation can be known (e.g., log or sqrt) or unknown (Box-Cox or estimated nonparametrically) for greater flexibility.

Usage

star_MCMC(
  y,
  sample_params,
  init_params,
  transformation = "np",
  y_max = Inf,
  nsave = 5000,
  nburn = 5000,
  nskip = 2,
  save_y_hat = FALSE,
  verbose = TRUE
)

Arguments

y	n x 1 vector of observed counts
sample_params	a function that inputs data y and a named list params containing mu: the n x 1 vector of conditional means (fitted values) sigma: the conditional standard deviation coefficients: a named list of parameters that determine mu and outputs an updated list params of samples from the full conditional posterior distribution of coefficients and sigma (and updates mu)
init_params	an initializing function that inputs data y and initializes the named list params of mu, sigma, and coefficients
transformation	transformation to use for the latent data; must be one of "identity" (identity transformation) "log" (log transformation) "sqrt" (square root transformation) "np" (nonparametric transformation estimated from empirical CDF) "pois" (transformation for moment-matched marginal Poisson CDF) "neg-bin" (transformation for moment-matched marginal Negative Binomial CDF) "box-cox" (box-cox transformation with learned parameter)
y_max	a fixed and known upper bound for all observations; default is Inf
nsave	number of MCMC iterations to save
nburn	number of MCMC iterations to discard
nskip	number of MCMC iterations to skip between saving iterations, i.e., save every (nskip + 1)th draw
save_y_hat	logical; if TRUE, compute and save the posterior draws of the expected counts, E(y), which may be slow to compute
verbose	logical; if TRUE, print time remaining

Details

STAR defines a count-valued probability model by (1) specifying a Gaussian model for continuous *latent* data and (2) connecting the latent data to the observed data via a *transformation and rounding* operation.

Posterior and predictive inference is obtained via a Gibbs sampler that combines (i) a latent data augmentation step (like in probit regression) and (ii) an existing sampler for a continuous data model.

There are several options for the transformation. First, the transformation can belong to the *Box-Cox* family, which includes the known transformations 'identity', 'log', and 'sqrt', as well as a version in which the Box-Cox parameter is inferred within the MCMC sampler ('box-cox'). Second, the transformation can be estimated (before model fitting) using the empirical distribution of the data y. Options in this case include the empirical cumulative distribution function (CDF), which is fully nonparametric ('np'), or the parametric alternatives based on Poisson ('pois') or Negative-Binomial ('neg-bin') distributions. For the parametric distributions, the parameters of the distribution are estimated using moments (means and variances) of y.

Value

a list with the following elements:

• coefficients: the posterior mean of the coefficients

• fitted.values: the posterior mean of the conditional expectation of the counts y

• post.coefficients: posterior draws of the coefficients

• post.fitted.values: posterior draws of the conditional mean of the counts y

• post.pred: draws from the posterior predictive distribution of y

• post.lambda: draws from the posterior distribution of lambda

• post.sigma: draws from the posterior distribution of sigma

• post.log.like.point: draws of the log-likelihood for each of the n observations

• logLik: the log-likelihood evaluated at the posterior means

• WAIC: Widely-Applicable/Watanabe-Akaike Information Criterion

• p_waic: Effective number of parameters based on WAIC

Examples

## Not run: 
# Simulate data with count-valued response y:
sim_dat = simulate_nb_lm(n = 100, p = 5)
y = sim_dat$y; X = sim_dat$X

# STAR: log-transformation:
fit_log = star_MCMC(y = y,
                         sample_params = function(y, params) sample_params_lm(y, X, params),
                         init_params = function(y) init_params_lm(y, X),
                         transformation = 'log')
# Posterior mean of each coefficient:
coef(fit_log)

# WAIC for STAR-log:
fit_log$WAIC

# MCMC diagnostics:
plot(as.ts(fit_log$post.coefficients[,1:3]))

# Posterior predictive check:
hist(apply(fit_log$post.pred, 1,
           function(x) mean(x==0)), main = 'Proportion of Zeros', xlab='');
abline(v = mean(y==0), lwd=4, col ='blue')

# STAR: nonparametric transformation
fit = star_MCMC(y = y,
                sample_params = function(y, params) sample_params_lm(y, X, params),
                init_params = function(y) init_params_lm(y, X),
                transformation = 'np')

# Posterior mean of each coefficient:
coef(fit)

# WAIC:
fit$WAIC

# MCMC diagnostics:
plot(as.ts(fit$post.coefficients[,1:3]))

# Posterior predictive check:
hist(apply(fit$post.pred, 1,
           function(x) mean(x==0)), main = 'Proportion of Zeros', xlab='');
abline(v = mean(y==0), lwd=4, col ='blue')


## End(Not run)

drkowal/rSTAR documentation built on July 5, 2023, 2:18 p.m.