Brent's method for optimization

Description

Implementation for Brent's algorithm for minimizing a univariate function over an interval. The code is based on a function in the stsm package.

Usage

BrentMethod(a = 0, b, fcn, tol = .Machine$double.eps^0.25)

Arguments

Value

a list of containing the following elements:

• fx the minimum value of the input function

• x the argument that minimizes the function

• iter number of iterations to converge

• vx a vector that stores the arguments until convergence

Inverse rounding function

Description

Define the intervals associated with y = j based on the flooring function. The function returns -Inf for j = 0 (or smaller) and Inf for any j >= y_max + 1, where y_max is a known upper bound on the data y (if specified).

Usage

a_j(j, y_max = Inf)

Arguments

Value

The (lower) interval endpoint(s) associated with j.

Examples

# Standard cases:
a_j(1)
a_j(20)

# Boundary cases:
a_j(0)
a_j(20, y_max = 15)

Fit Bayesian Additive STAR Model with MCMC

Description

Run the MCMC algorithm for a STAR Bayesian additive model The transformation can be known (e.g., log or sqrt) or unknown (Box-Cox or estimated nonparametrically) for greater flexibility.

Usage

bam_star(
  y,
  X_lin,
  X_nonlin,
  splinetype = "orthogonal",
  transformation = "np",
  y_max = Inf,
  nsave = 5000,
  nburn = 5000,
  nskip = 2,
  save_y_hat = FALSE,
  verbose = TRUE
)

Arguments

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. Third, the transformation can be modeled as an unknown, monotone function using I-splines ('ispline'). The Robust Adaptive Metropolis (RAM) sampler is used for drawing the parameter of the transformation function.

Value

a list with at least 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

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

• p_waic: Effective number of parameters based on WAIC

In the case of transformation="ispline", the list also contains

• post.g: draws from the posterior distribution of the transformation g

• post.sigma.gamma: draws from the posterior distribution of sigma.gamma, the prior standard deviation of the transformation g() coefficients

Examples

# Simulate data with count-valued response y:
sim_dat = simulate_nb_friedman(n = 100, p = 5, seed=32)
y = sim_dat$y; X = sim_dat$X

# Linear and nonlinear components:
X_lin = as.matrix(X[,-(1:3)])
X_nonlin = as.matrix(X[,(1:3)])

# STAR: nonparametric transformation
fit <- bam_star(y,X_lin, X_nonlin, nburn=1000, nskip=0)

# 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')

MCMC Algorithm for BART-STAR

Description

Run the MCMC algorithm for a BART model for count-valued responses using STAR. The transformation can be known (e.g., log or sqrt) or unknown (Box-Cox or estimated nonparametrically) for greater flexibility.

Usage

bart_star(
  y,
  X,
  X_test = NULL,
  y_test = NULL,
  transformation = "np",
  y_max = Inf,
  n.trees = 200,
  sigest = NULL,
  sigdf = 3,
  sigquant = 0.9,
  k = 2,
  power = 2,
  base = 0.95,
  nsave = 5000,
  nburn = 5000,
  nskip = 2,
  save_y_hat = FALSE,
  verbose = TRUE
)

Arguments

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. Here, the model in (1) is a Bayesian additive regression tree (BART) model.

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. Third, the transformation can be modeled as an unknown, monotone function using I-splines ('ispline'). The Robust Adaptive Metropolis (RAM) sampler is used for drawing the parameter of the transformation function.

Value

a list with the following elements:

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

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

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

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

• p_waic: Effective number of parameters based on WAIC

• post.pred.test: draws from the posterior predictive distribution at the test points X_test (NULL if X_test is not given)

• post.fitted.values.test: posterior draws of the conditional mean at the test points X_test (NULL if X_test is not given)

• post.mu.test: draws of the conditional mean of z_star at the test points X_test (NULL if X_test is not given)

• post.log.pred.test: draws of the log-predictive distribution for each of the n0 test cases (NULL if X_test is not given)

• fitted.values: the posterior mean of the conditional expectation of the counts y (NULL if save_y_hat=FALSE)

• post.fitted.values: posterior draws of the conditional mean of the counts y (NULL if save_y_hat=FALSE)

In the case of transformation="ispline", the list also contains

• post.g: draws from the posterior distribution of the transformation g

• post.sigma.gamma: draws from the posterior distribution of sigma.gamma, the prior standard deviation of the transformation g() coefficients

If transformation="box-cox", then the list also contains

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

Examples

# Simulate data with count-valued response y:
sim_dat = simulate_nb_friedman(n = 100, p = 10)
y = sim_dat$y; X = sim_dat$X

# BART-STAR with log-transformation:
fit_log = bart_star(y = y, X = X, transformation = 'log',
                    save_y_hat = TRUE, nburn=1000, nskip=0)

# Fitted values
plot_fitted(y = sim_dat$Ey,
            post_y = fit_log$post.fitted.values,
            main = 'Fitted Values: BART-STAR-log')

# WAIC for BART-STAR-log:
fit_log$WAIC

# MCMC diagnostics:
plot(as.ts(fit_log$post.fitted.values[,1:10]))

# 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')

# BART-STAR with nonparametric transformation:
fit = bart_star(y = y, X = X,
                     transformation = 'np', save_y_hat = TRUE)

# Fitted values
plot_fitted(y = sim_dat$Ey,
            post_y = fit$post.fitted.values,
            main = 'Fitted Values: BART-STAR-np')

# WAIC for BART-STAR-np:
fit$WAIC

# MCMC diagnostics:
plot(as.ts(fit$post.fitted.values[,1:10]))

# 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')

MCMC sampler for BART-STAR with a monotone spline model for the transformation

Description

Run the MCMC algorithm for BART model for count-valued responses using STAR. The transformation is modeled as an unknown, monotone function using I-splines. The Robust Adaptive Metropolis (RAM) sampler is used for drawing the parameter of the transformation function.

Usage

bart_star_ispline(
  y,
  X,
  X_test = NULL,
  y_test = NULL,
  lambda_prior = 1/2,
  y_max = Inf,
  n.trees = 200,
  sigest = NULL,
  sigdf = 3,
  sigquant = 0.9,
  k = 2,
  power = 2,
  base = 0.95,
  nsave = 5000,
  nburn = 5000,
  nskip = 2,
  save_y_hat = FALSE,
  target_acc_rate = 0.3,
  adapt_rate = 0.75,
  stop_adapt_perc = 0.5,
  verbose = TRUE
)

Arguments

Value

a list with the following elements:

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

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

• post.pred.test: draws from the posterior predictive distribution at the test points X_test

• post.fitted.values.test: posterior draws of the conditional mean at the test points X_test

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

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

• post.mu.test: draws of the conditional mean of z_star at the test points

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

• post.log.pred.test: draws of the log-predictive distribution for each of the n0 test cases

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

• p_waic: Effective number of parameters based on WAIC

• post.g: draws from the posterior distribution of the transformation g

• post.sigma.gamma: draws from the posterior distribution of sigma.gamma, the prior standard deviation of the transformation g coefficients

STAR Bayesian Linear Regression

Description

Posterior inference for STAR linear model

Usage

blm_star(
  y,
  X,
  X_test = NULL,
  transformation = "np",
  y_max = Inf,
  prior = "gprior",
  use_MCMC = TRUE,
  nsave = 5000,
  nburn = 5000,
  nskip = 0,
  method_sigma = "mle",
  approx_Fz = FALSE,
  approx_Fy = FALSE,
  psi = NULL,
  compute_marg = FALSE
)

Arguments

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. Here, the continuous latent data model is a linear regression.

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. The distribution-based transformations approximately preserve the mean and variance of the count data y on the latent data scale, which lends interpretability to the model parameters. Lastly, the transformation can be modeled using the Bayesian bootstrap ('bnp'), which is a Bayesian nonparametric model and incorporates the uncertainty about the transformation into posterior and predictive inference.

The Monte Carlo sampler (use_MCMC=FALSE) produces direct, discrete, and joint draws from the posterior distribution and the posterior predictive distribution of the linear regression model with a g-prior.

Value

a list with at least the following elements:

• coefficients: the posterior mean of the regression coefficients

• post.beta: posterior draws of the regression coefficients

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

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

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

• p_waic: Effective number of parameters based on WAIC

If test points are passed in, then the list will also have post.predtest, which contains draws from the posterior predictive distribution at test points.

Other elements may be present depending on the choice of prior, transformation, and sampling approach.

Note

The 'bnp' transformation (without the Fy approximation) is slower than the other transformations because of the way the TruncatedNormal sampler must be updated as the lower and upper limits change (due to the sampling of g). Thus, computational improvements are likely available.

Gibbs sampler for STAR linear regression with BNP transformation

Description

Compute MCMC samples from the posterior and predictive distributions of a STAR linear regression model with a g-prior and BNP transformation.

Usage

blm_star_bnpgibbs(
  y,
  X,
  X_test = X,
  y_max = Inf,
  psi = NULL,
  approx_Fz = FALSE,
  approx_Fy = FALSE,
  nsave = 1000,
  nburn = 1000,
  nskip = 0,
  verbose = TRUE
)

Arguments

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. Here, the continuous latent data model is a linear regression.

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'. 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. The distribution-based transformations approximately preserve the mean and variance of the count data y on the latent data scale, which lends interpretability to the model parameters. Lastly, the transformation can be modeled using the Bayesian bootstrap ('bnp'), which is a Bayesian nonparametric model and incorporates the uncertainty about the transformation into posterior and predictive inference.

Value

a list with the following elements:

• coefficients the posterior mean of the regression coefficients

• post_beta: nsave x p samples from the posterior distribution of the regression coefficients

• post_ytilde: nsave x n0 samples from the posterior predictive distribution at test points X_test

• post_g: nsave posterior samples of the transformation evaluated at the unique y values (only applies for 'bnp' transformations)

Monte Carlo sampler for STAR linear regression with a g-prior

Description

Compute direct Monte Carlo samples from the posterior and predictive distributions of a STAR linear regression model with a g-prior.

Usage

blm_star_exact(
  y,
  X,
  X_test = X,
  transformation = "np",
  y_max = Inf,
  psi = NULL,
  method_sigma = "mle",
  approx_Fz = FALSE,
  approx_Fy = FALSE,
  nsave = 5000,
  compute_marg = FALSE
)

Arguments

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. Here, the continuous latent data model is a linear regression.

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'. 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. The distribution-based transformations approximately preserve the mean and variance of the count data y on the latent data scale, which lends interpretability to the model parameters. Lastly, the transformation can be modeled using the Bayesian bootstrap ('bnp'), which is a Bayesian nonparametric model and incorporates the uncertainty about the transformation into posterior and predictive inference.

The Monte Carlo sampler produces direct, discrete, and joint draws from the posterior distribution and the posterior predictive distribution of the linear regression model with a g-prior.

Value

a list with the following elements:

• coefficients the posterior mean of the regression coefficients

• post.beta: nsave x p samples from the posterior distribution of the regression coefficients

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

• post.pred.test: nsave x n0 samples from the posterior predictive distribution at test points X_test (if given, otherwise NULL)

• sigma: The estimated latent data standard deviation

• post.g: nsave posterior samples of the transformation evaluated at the unique y values (only applies for 'bnp' transformations)

• marg.like: the marginal likelihood (if requested; otherwise NULL)

Note

The 'bnp' transformation (without the Fy approximation) is slower than the other transformations because of the way the TruncatedNormal sampler must be updated as the lower and upper limits change (due to the sampling of g). Thus, computational improvements are likely available.

Estimate the remaining time in the MCMC based on previous samples

Description

Estimate the remaining time in the MCMC based on previous samples

Usage

computeTimeRemaining(nsi, timer0, nsims, nrep = 1000)

Arguments

Value

Table of summary statistics using the function summary

Compute asymptotic confidence intervals for STAR linear regression

Description

For a linear regression model within the STAR framework, compute (asymptotic) confidence intervals for a regression coefficient of interest. Confidence intervals are computed by inverting the likelihood ratio test and profiling the log-likelihood.

Usage

## S3 method for class 'lmstar'
confint(object, parm, level = 0.95, ...)

Arguments

Value

A matrix (or vector) with columns giving lower and upper confidence limits for each parameter. These will be labelled as (1-level)/2 and 1 - (1-level)/2 in

Examples

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

#Select a transformation:
transformation = 'np'

#Estimate model
fit = lm_star(y~X, transformation=transformation)

#Confidence interval for all parameters
confint(fit)

Compute Simultaneous Credible Bands

Description

Compute (1-alpha)% credible BANDS for a function based on MCMC samples using Crainiceanu et al. (2007)

Usage

credBands(sampFuns, alpha = 0.05)

Arguments

Value

m x 2 matrix of credible bands; the first column is the lower band, the second is the upper band

Note

The input needs not be curves: the simultaneous credible "bands" may be computed for vectors. The resulting credible intervals will provide joint coverage at the (1-alpha) level across all components of the vector.

Compute the ergodic (running) mean.

Description

Compute the ergodic (running) mean.

Usage

ergMean(x)

Arguments

Value

A vector y with each element defined by y[i] = mean(x[1:i])

Examples

# Compare:
ergMean(1:10)
mean(1:10)

# Running mean for iid N(5, 1) samples:
x = rnorm(n = 10^4, mean = 5, sd = 1)
plot(ergMean(x))
abline(h=5)

Compute E(Y^2) for a STAR process

Description

Compute the conditional expectation of Y^2 for a STAR process Y under a generic link function g.

Usage

expectation2_gRcpp(g_a_j, g_a_jp1, mu, sigma, Jmax)

Arguments

Value

y2_hat m x 1 vector of conditional expectations

Note

This function uses Rcpp for computational efficiency.

Estimate the mean for a STAR process

Description

Estimate the conditional expectation for a STAR process under a generic link function g.

Usage

expectation_gRcpp(g_a_j, g_a_jp1, mu, sigma, Jmax)

Arguments

Value

y_hat m x 1 vector of conditional expectations

Note

This function uses Rcpp for computational efficiency.

Estimate the mean for a STAR process

Description

Estimate the conditional expectation for a STAR process under the identity link function.

Usage

expectation_identity(a, Jmax, Mu, sigma_t, Offset)

Arguments

Value

Zhat T x m matrix of conditional expectations

Note

This function uses Rcpp for computational efficiency.

Estimate the mean for a STAR process

Description

Estimate the conditional expectation for a STAR process under the log link function.

Usage

expectation_log(a, Jmax, Mu, sigma_t, Offset)

Arguments

Value

Zhat T x m matrix of conditional expectations

Note

This function uses Rcpp for computational efficiency.

Estimate the mean for a STAR process

Description

Estimate the conditional expectation for a STAR process under the square root link function.

Usage

expectation_sqrt(a, Jmax, Mu, sigma_t, Offset)

Arguments

Value

Zhat T x m matrix of conditional expectations

Note

This function uses Rcpp for computational efficiency.

Box-Cox transformation

Description

Evaluate the Box-Cox transformation, which is a scaled power transformation to preserve continuity in the index lambda at zero. Negative values are permitted.

Usage

g_bc(t, lambda)

Arguments

Value

The evaluation(s) of the Box-Cox function at the given input(s) t.

Note

Special cases include the identity transformation (lambda = 1), the square-root transformation (lambda = 1/2), and the log transformation (lambda = 0).

Examples

# Log-transformation:
g_bc(1:5, lambda = 0); log(1:5)

# Square-root transformation: note the shift and scaling
g_bc(1:5, lambda = 1/2); sqrt(1:5)

Bayesian bootstrap-based transformation

Description

Compute one posterior draw from the smoothed transformation implied by (separate) Bayesian bootstrap models for the CDFs of y and X.

Usage

g_bnp(
  y,
  xtSigmax = rep(0, length(y)),
  zgrid = NULL,
  sigma_epsilon = 1,
  approx_Fz = FALSE
)

Arguments

Value

A smooth monotone function which can be used for evaluations of the transformation at each posterior draw.

Examples

# Sample some data:
y = rpois(n = 200, lambda = 5)
# Compute 200 draws of g on a grid:
t = seq(0, max(y), length.out = 100) # grid
g_post = t(sapply(1:500, function(s) g_bnp(y, approx_Fz = TRUE)(t)))
# Plot together:
plot(t, t, ylim = range(g_post), type='n', ylab = 'g(t)',  main = 'Bayesian bootstrap posterior: g')
apply(g_post, 1, function(g) lines(t, g, col='gray'))
# And the posterior mean of g:
lines(t, colMeans(g_post), lwd=3)

Cumulative distribution function (CDF)-based transformation

Description

Compute a CDF-based transformation using the observed count data. The CDF can be estimated nonparametrically or parametrically based on the Poisson or Negative Binomial distributions. In the parametric case, the parameters are determined based on the moments of y. Note that this is a fixed quantity and does not come with uncertainty quantification.

Usage

g_cdf(y, distribution = "np")

Arguments

Value

A smooth monotone function which can be used for evaluations of the transformation.

Examples

# Sample some data:
y = rpois(n = 500, lambda = 5)

# Empirical CDF version:
g_np = g_cdf(y, distribution = 'np')

# Poisson version:
g_pois = g_cdf(y, distribution = 'pois')

# Negative binomial version:
g_negbin = g_cdf(y, distribution = 'neg-bin')

# Plot together:
t = 1:max(y) # grid
plot(t, g_np(t), type='l')
lines(t, g_pois(t), lty = 2)
lines(t, g_negbin(t), lty = 3)

Approximate inverse transformation

Description

Compute the inverse function of a transformation g based on a grid search.

Usage

g_inv_approx(g, t_grid)

Arguments

Value

A function which can be used for evaluations of the (approximate) inverse transformation function.

Examples

# Sample some data:
y = rpois(n = 500, lambda = 5)

# Empirical CDF transformation:
g_np = g_cdf(y, distribution = 'np')

# Grid for approximation:
t_grid = seq(1, max(y), length.out = 100)

# Approximate inverse:
g_inv = g_inv_approx(g = g_np, t_grid = t_grid)

# Check the approximation:
plot(t_grid, g_inv(g_np(t_grid)), type='p')
lines(t_grid, t_grid)

Inverse Box-Cox transformation

Description

Evaluate the inverse Box-Cox transformation. Negative values are permitted.

Usage

g_inv_bc(s, lambda)

Arguments

Value

The evaluation(s) of the inverse Box-Cox function at the given input(s) s.

Note

Special cases include the identity transformation (lambda = 1), the square-root transformation (lambda = 1/2), and the log transformation (lambda = 0).

#' @examples # (Inverse) log-transformation: g_inv_bc(1:5, lambda = 0); exp(1:5)

# (Inverse) square-root transformation: note the shift and scaling g_inv_bc(1:5, lambda = 1/2); (1:5)^2

Fitting STAR Gradient Boosting Machines via EM algorithm

Description

Compute the MLEs and log-likelihood for the Gradient Boosting Machines (GBM) STAR model. The STAR model requires a *transformation* and an *estimation function* for the conditional mean given observed data. The transformation can be known (e.g., log or sqrt) or unknown (Box-Cox or estimated nonparametrically) for greater flexibility. The estimator in this case is a GBM. Standard function calls including fitted and residuals apply.

Usage

gbm_star(
  y,
  X,
  X.test = NULL,
  transformation = "np",
  y_max = Inf,
  sd_init = 10,
  tol = 10^-10,
  max_iters = 1000,
  n.trees = 100,
  interaction.depth = 1,
  shrinkage = 0.1,
  bag.fraction = 1
)

Arguments

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. The Gaussian model in this case is a GBM.

Value

a list with the following elements:

• fitted.values: the fitted values at the MLEs (training)

• fitted.values.test: the fitted values at the MLEs (testing)

• g.hat a function containing the (known or estimated) transformation

• sigma.hat the MLE of the standard deviation

• mu.hat the MLE of the conditional mean (on the transformed scale)

• z.hat the estimated latent data (on the transformed scale) at the MLEs

• residuals the Dunn-Smyth residuals (randomized)

• residuals_rep the Dunn-Smyth residuals (randomized) for 10 replicates

• logLik the log-likelihood at the MLEs

• logLik0 the log-likelihood at the MLEs for the *unrounded* initialization

• lambda the Box-Cox nonlinear parameter

• gbmObj: the object returned by gbm() at the MLEs

• and other parameters that (1) track the parameters across EM iterations and (2) record the model specifications

Note

Infinite latent data values may occur when the transformed Gaussian model is highly inadequate. In that case, the function returns the *indices* of the data points with infinite latent values, which are significant outliers under the model. Deletion of these indices and re-running the model is one option, but care must be taken to ensure that (i) it is appropriate to treat these observations as outliers and (ii) the model is adequate for the remaining data points.

References

Kowal, D. R., & Wu, B. (2021). Semiparametric count data regression for self‐reported mental health. Biometrics. doi:10.1111/biom.13617

Examples

# Simulate data with count-valued response y:
sim_dat = simulate_nb_friedman(n = 100, p = 10)
y = sim_dat$y; X = sim_dat$X

# EM algorithm for STAR (using the log-link)
fit_em = gbm_star(y = y, X = X,
                 transformation = 'log')

# Evaluate convergence:
plot(fit_em$logLik_all, type='l', main = 'GBM-STAR-log', xlab = 'Iteration', ylab = 'log-lik')

# Fitted values:
y_hat = fitted(fit_em)
plot(y_hat, y);

# Residuals:
plot(residuals(fit_em))
qqnorm(residuals(fit_em)); qqline(residuals(fit_em))

# Log-likelihood at MLEs:
fit_em$logLik

Generalized EM estimation for STAR

Description

Compute MLEs and log-likelihood for a generalized STAR model. The STAR model requires a *transformation* and an *estimation function* for the conditional mean given observed data. The transformation can be known (e.g., log or sqrt) or unknown (Box-Cox or estimated nonparametrically) for greater flexibility. The estimator can be any least squares estimator, including nonlinear models. Standard function calls including coefficients(), fitted(), and residuals() apply.

Usage

genEM_star(
  y,
  estimator,
  transformation = "np",
  y_max = Inf,
  sd_init = 10,
  tol = 10^-10,
  max_iters = 1000
)

Arguments

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.

The expectation-maximization (EM) algorithm is used to produce maximum likelihood estimators (MLEs) for the parameters defined in the estimator function, such as linear regression coefficients, which define the Gaussian model for the continuous latent data. Fitted values (point predictions), residuals, and log-likelihood values are also available. Inference for the estimators proceeds via classical maximum likelihood. Initialization of the EM algorithm can be randomized to monitor convergence. However, the log-likelihood is concave for all transformations (except 'box-cox'), so global convergence is guaranteed.

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 estimated within the EM algorithm ('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 MLEs of the coefficients

• fitted.values the fitted values at the MLEs

• g.hat a function containing the (known or estimated) transformation

• sigma.hat the MLE of the standard deviation

• mu.hat the MLE of the conditional mean (on the transformed scale)

• z.hat the estimated latent data (on the transformed scale) at the MLEs

• residuals the Dunn-Smyth residuals (randomized)

• residuals_rep the Dunn-Smyth residuals (randomized) for 10 replicates

• logLik the log-likelihood at the MLEs

• logLik0 the log-likelihood at the MLEs for the *unrounded* initialization

• lambda the Box-Cox nonlinear parameter

• and other parameters that (1) track the parameters across EM iterations and (2) record the model specifications

Note

Infinite latent data values may occur when the transformed Gaussian model is highly inadequate. In that case, the function returns the *indices* of the data points with infinite latent values, which are significant outliers under the model. Deletion of these indices and re-running the model is one option, but care must be taken to ensure that (i) it is appropriate to treat these observations as outliers and (ii) the model is adequate for the remaining data points.

References

Kowal, D. R., & Wu, B. (2021). Semiparametric count data regression for self‐reported mental health. Biometrics. doi:10.1111/biom.13617

Examples

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

# Select a transformation:
transformation = 'np'

# Example using GAM as underlying estimator (for illustration purposes only)
if(require("mgcv")){
  fit_em = genEM_star(y = y,
                      estimator = function(y) gam(y ~ s(X1)+s(X2),
                      data=data.frame(y,X)),
                      transformation = transformation)
}

# Fitted coefficients:
coef(fit_em)

# Fitted values:
y_hat = fitted(fit_em)
plot(y_hat, y);

# Log-likelihood at MLEs:
fit_em$logLik

Generalized 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

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

Arguments

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 at least the following elements:

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

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

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

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

• p_waic: Effective number of parameters based on WAIC

• post.lambda: draws from the posterior distribution of lambda (NULL unless transformation='box-cox')

• fitted.values: the posterior mean of the conditional expectation of the counts y (NULL if save_y_hat=FALSE)

• post.fitted.values: posterior draws of the conditional mean of the counts y (NULL if save_y_hat=FALSE)

If the coefficients list from init_params and sample_params contains a named element beta, e.g. for linear regression, then the function output contains

• coefficients: the posterior mean of the beta coefficients

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

• post.othercoefs: draws from the posterior distribution of any other sampled coefficients, e.g. variance terms

If no beta exists in the parameter coefficients, then the output list just contains

• coefficients: the posterior mean of all coefficients

• post.beta: draws from the posterior distribution of all coefficients

Additionally, if init_params and sample_params have output mu_test, then the sampler will output post.predtest, which contains draws from the posterior predictive distribution at test points.

Examples

# 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 = genMCMC_star(y = y,
                         sample_params = function(y, params) sample_lm_gprior(y, X, params),
                         init_params = function(y) init_lm_gprior(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.beta[,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')

MCMC sampler for STAR with a monotone spline model for the transformation

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 is modeled as an unknown, monotone function using I-splines. The Robust Adaptive Metropolis (RAM) sampler is used for drawing the parameter of the transformation function.

Usage

genMCMC_star_ispline(
  y,
  sample_params,
  init_params,
  lambda_prior = 1/2,
  y_max = Inf,
  nsave = 5000,
  nburn = 5000,
  nskip = 0,
  save_y_hat = FALSE,
  target_acc_rate = 0.3,
  adapt_rate = 0.75,
  stop_adapt_perc = 0.5,
  verbose = TRUE
)

Arguments

Details

If the coefficients list from init_params and sample_params contains a named element beta, e.g. for linear regression, then the function output contains

• coefficients: the posterior mean of the beta coefficients

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

• post.othercoefs: draws from the posterior distribution of any other sampled coefficients, e.g. variance terms

If no beta exists in the parameter coefficients, then the output list just contains

• coefficients: the posterior mean of all coefficients

• post.beta: draws from the posterior distribution of all coefficients

Additionally, if init_params and sample_params have output mu_test, then the sampler will output post.predtest, which contains draws from the posterior predictive distribution at test points.

Value

A list with at least the following elements:

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

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

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

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

• p_waic: Effective number of parameters based on WAIC

• post.g: draws from the posterior distribution of the transformation g

• post.sigma.gamma: draws from the posterior distribution of sigma.gamma, the prior standard deviation of the transformation g() coefficients

• fitted.values: the posterior mean of the conditional expectation of the counts y (NULL if save_y_hat=FALSE)

• post.fitted.values: posterior draws of the conditional mean of the counts y (NULL if save_y_hat=FALSE)

along with other elements depending on the nature of the initialization and sampling functions. See details for more info.

Summarize of effective sample size

Description

Compute the summary statistics for the effective sample size (ESS) across posterior samples for possibly many variables

Usage

getEffSize(postX)

Arguments

Value

Table of summary statistics using the function summary().

Examples

# ESS for iid simulations:
rand_iid = rnorm(n = 10^4)
getEffSize(rand_iid)

# ESS for several AR(1) simulations with coefficients 0.1, 0.2,...,0.9:
rand_ar1 = sapply(seq(0.1, 0.9, by = 0.1), function(x) arima.sim(n = 10^4, list(ar = x)))
getEffSize(rand_ar1)

Initialize the parameters for an additive model

Description

Initialize the parameters for an additive model, which may contain both linear and nonlinear predictors. The nonlinear terms are modeled using orthogonalized splines.

Usage

init_bam_orthog(y, X_lin, X_nonlin, B_all = NULL)

Arguments

Value

a named list params containing

1. mu: vector of conditional means (fitted values)

2. sigma: the conditional standard deviation

3. coefficients: a named list of parameters that determine mu

Note

The parameters in coefficients are:

• beta_lin: the p x 1 linear coefficients, including the linear terms from X_nonlin

• f_j: the n x pNL matrix of fitted values for each nonlinear function

• theta_j: the pNL-dimensional of nonlinear basis coefficients

• sigma_beta: p x 1 vector of linear regression coefficient standard deviations

• sigma_theta_j: pNL x 1 vector of nonlinear coefficient standard deviations

Initialize the parameters for an additive model

Description

Initialize the parameters for an additive model, which may contain both linear and nonlinear predictors. The nonlinear terms are modeled using low-rank thin plate splines.

Usage

init_bam_thin(y, X_lin, X_nonlin, B_all = NULL)

Arguments

Value

a named list params containing

1. mu: vector of conditional means (fitted values)

2. sigma: the conditional standard deviation

3. coefficients: a named list of parameters that determine mu

Note

The parameters in coefficients are:

• beta_lin: the p x 1 linear coefficients, including the linear terms from X_nonlin

• f_j: the n x pNL matrix of fitted values for each nonlinear function

• theta_j: the pNL-dimensional of nonlinear basis coefficients

• sigma_beta: p x 1 vector of linear regression coefficient standard deviations

• sigma_theta_j: pNL x 1 vector of nonlinear coefficient standard deviations

Initialize linear regression parameters assuming a g-prior

Description

Initialize the parameters for a linear regression model assuming a g-prior for the coefficients.

Usage

init_lm_gprior(y, X, X_test = NULL)

Arguments

Value

a named list params containing at least

1. mu: vector of conditional means (fitted values)

2. sigma: the conditional standard deviation

3. coefficients: a named list of parameters that determine mu

Additionally, if X_test is not NULL, then the list includes an element mu_test, the vector of conditional means at the test points

Note

The parameters in coefficients are:

• beta: the p x 1 vector of regression coefficients components of beta

Examples

# Simulate data for illustration:
sim_dat = simulate_nb_lm(n = 100, p = 5)
y = sim_dat$y; X = sim_dat$X

# Initialize:
params = init_lm_gprior(y = y, X = X)
names(params)
names(params$coefficients)

Initialize linear regression parameters assuming a horseshoe prior

Description

Initialize the parameters for a linear regression model assuming a horseshoe prior for the (non-intercept) coefficients. The number of predictors p may exceed the number of observations n.

Usage

init_lm_hs(y, X, X_test = NULL)

Arguments

Value

a named list params containing at least

1. mu: vector of conditional means (fitted values)

2. sigma: the conditional standard deviation

3. coefficients: a named list of parameters that determine mu

Additionally, if X_test is not NULL, then the list includes an element mu_test, the vector of conditional means at the test points

Note

The parameters in coefficients are:

• beta: the p x 1 vector of regression coefficients

• sigma_beta: the p x 1 vector of regression coefficient standard deviations (local scale parameters)

• xi_sigma_beta: the p x 1 vector of parameter-expansion variables for sigma_beta

• lambda_beta: the global scale parameter

• xi_lambda_beta: the parameter-expansion variable for lambda_beta components of beta

Initialize linear regression parameters assuming a ridge prior

Description

Initialize the parameters for a linear regression model assuming a ridge prior for the (non-intercept) coefficients. The number of predictors p may exceed the number of observations n.

Usage

init_lm_ridge(y, X, X_test = NULL)

Arguments

Value

a named list params containing at least

1. mu: vector of conditional means (fitted values)

2. sigma: the conditional standard deviation

3. coefficients: a named list of parameters that determine mu

Additionally, if X_test is not NULL, then the list includes an element mu_test, the vector of conditional means at the test points

Note

The parameters in coefficients are:

• beta: the p x 1 vector of regression coefficients

• sigma_beta: the prior standard deviation for the (non-intercept) components of beta

Initialize the parameters for a simple mean-only model

Description

Initialize the parameters for the model y ~ N(mu0, sigma^2) with a flat prior on mu0.

Usage

init_params_mean(y)

Arguments

Value

a named list params containing

1. mu: vector of conditional means (fitted values)

2. sigma: the conditional standard deviation

3. coefficients: a named list of parameters that determine mu

Note

The only parameter in coefficients is mu0. Although redundant here, this parametrization is useful in other functions.

Estimate confidence intervals/bands for a STAR process

Description

Compute confidence intervals/bands for the expected value of the count-valued STAR process y based on intervals/bands for the Gaussian process mu.

Usage

interval_gRcpp(g_a_j, g_a_jp1, L_mu, U_mu, sigma, Jmax)

Arguments

Value

LU_y m x 2 vector of intervals for y.

Note

This function uses Rcpp for computational efficiency.

Compute the inverse log-odds

Description

Compute the inverse log-odds

Usage

invlogit(x)

Arguments

Value

A scalar or vector of values in (0,1)

Fitting frequentist STAR linear model via EM algorithm

Description

Compute the MLEs and log-likelihood for the STAR linear model. The regression coefficients are estimated using least squares within an EM algorithm.

Usage

lm_star(
  formula,
  data = NULL,
  transformation = "np",
  y_max = Inf,
  sd_init = 10,
  tol = 10^-10,
  max_iters = 1000
)

Arguments

Details

Standard function calls including coefficients, fitted, and residuals apply. Fitted values are the expectation at the MLEs, and as such are not necessarily count-valued.

Value

an object of class "lmstar", which is a list with the following elements:

• coefficients the MLEs of the coefficients

• fitted.values the fitted values at the MLEs

• g.hat a function containing the (known or estimated) transformation

• ginv.hat a function containing the inverse of the transformation

• sigma.hat the MLE of the standard deviation

• mu.hat the MLE of the conditional mean (on the transformed scale)

• z.hat the estimated latent data (on the transformed scale) at the MLEs

• residuals the Dunn-Smyth residuals (randomized)

• residuals_rep the Dunn-Smyth residuals (randomized) for 10 replicates

• logLik the log-likelihood at the MLEs

• logLik0 the log-likelihood at the MLEs for the *unrounded* initialization

• lambda the Box-Cox nonlinear parameter

• and other parameters that (1) track the parameters across EM iterations and (2) record the model specifications

Note

Infinite latent data values may occur when the transformed Gaussian model is highly inadequate. In that case, the function returns the *indices* of the data points with infinite latent values, which are significant outliers under the model. Deletion of these indices and re-running the model is one option, but care must be taken to ensure that (i) it is appropriate to treat these observations as outliers and (ii) the model is adequate for the remaining data points.

References

Kowal, D. R., & Wu, B. (2021). Semiparametric count data regression for self‐reported mental health. Biometrics. doi:10.1111/biom.13617

Examples

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

# Fit model
fit_em = lm_star(y~X)

# Fitted coefficients:
coef(fit_em)
# Fitted values:
y_hat = fitted(fit_em)
plot(y_hat, y);

# Residuals:
plot(residuals(fit_em))
qqnorm(residuals(fit_em)); qqline(residuals(fit_em))

Compute the pointwise log-likelihood for STAR

Description

Compute the pointwise log-likelihood for a STAR model. The code here assumes that the transformed real-valued process (z_star) has conditionally independent components with means mu and standard deviations sigma.

Usage

logLikePointRcpp(g_a_j, g_a_jp1, mu, sigma)

Arguments

Value

loglike m x 1 log-likelihood value

Note

This function uses Rcpp for computational efficiency.

Compute the log-likelihood for STAR

Description

Compute the log-likelihood for a STAR model. The code here assumes that the transformed real-valued process (z_star) has conditionally independent components with means mu and standard deviations sigma.

Usage

logLikeRcpp(g_a_j, g_a_jp1, mu, sigma)

Arguments

Value

loglike scalar log-likelihood value

Note

This function uses Rcpp for computational efficiency.

Compute the log-odds

Description

Compute the log-odds

Usage

logit(x)

Arguments

Value

A scalar or vector of log-odds

Plot the estimated regression coefficients and credible intervals

Description

Plot the estimated regression coefficients and credible intervals for the linear effects in up to two models.

Usage

plot_coef(
  post_coefficients_1,
  post_coefficients_2 = NULL,
  alpha = 0.05,
  labels = NULL
)

Arguments

Value

A plot of regression coefficients and credible intervals for 1-2 models

Plot the fitted values and the data

Description

Plot the fitted values, plus pointwise credible intervals, against the data. For simulations, one may use the true values in place of the data.

Usage

plot_fitted(y, post_y, y_hat = NULL, alpha = 0.05, ...)

Arguments

Value

A plot with the fitted values and the credible intervals against the data

Plot the empirical and model-based probability mass functions

Description

Plot the empirical probability mass function, i.e., the proportion of data values y that equal j for each j=0,1,..., together with the model-based estimate of the probability mass function based on the posterior predictive distribution.

Usage

plot_pmf(y, post.pred, error.bars = FALSE, alpha = 0.05)

Arguments

Value

A plot of the empirical PMF of y along with a PMF estimate from the model posterior predictive distribution

pmax() in Rcpp

Description

Compute the pointwise max for two vectors of equal length

Usage

pmaxRcpp(v1, v2)

Arguments

Value

vm m x 1 vector of pointwise maxima

Note

This function uses Rcpp for computational efficiency.

pmin() in Rcpp

Description

Compute the pointwise min for two vectors of equal length

Usage

pminRcpp(v1, v2)

Arguments

Value

vm m x 1 vector of pointwise minima

Note

This function uses Rcpp for computational efficiency.

Predict method for response in STAR linear model

Description

Outputs predicted values based on an lmstar fit and optionally prediction intervals based on the the (plug-in) predictive distribution for the STAR linear model

Usage

## S3 method for class 'lmstar'
predict(object, newdata = NULL, interval = FALSE, level = 0.95, N = 1000, ...)

Arguments

Details

If interval=TRUE, then predict.lmstar uses a Monte Carlo approach to estimating the (plug-in) predictive distribution for the STAR linear model. The algorithm iteratively samples (i) the latent data given the observed data, (ii) the latent predictive data given the latent data from (i), and (iii) (inverse) transforms and rounds the latent predictive data to obtain a draw from the integer-valued predictive distribution.

The appropriate quantiles of these Monte Carlo draws are computed and reported as the prediction interval.

Value

Either a a vector of predictions (if interval=FALSE) or a matrix of predictions and bounds with column names fit, lwr, and upr

Note

The “plug-in" predictive distribution is a crude approximation. Better approaches are available using the Bayesian models, e.g. blm_star, which provide samples from the posterior predictive distribution.

For highly skewed responses, prediction intervals especially at lower levels may not include the predicted value itself, since the mean is often much larger than the median.

Examples

# Simulate data with count-valued response y:
x = seq(0, 1, length.out = 100)
y = rpois(n = length(x), lambda = exp(1.5 + 5*(x -.5)^2))

# Estimate model--assume a quadratic effect (better for illustration purposes)
fit = lm_star(y~x+I(x^2), transformation = 'sqrt')

#Compute the predictive draws for the test points (same as observed points here)
#Also compute intervals using plug-in predictive distribution
y_pred = predict(fit, interval=TRUE)

# Plot the results
plot(x, y, ylim = range(y, y_pred), main = 'STAR: Predictions and 95% PI')
lines(x,y_pred[,"fit"], col='black', type='s', lwd=4)
lines(x, y_pred[,"lwr"], col='darkgray', type='s', lwd=4)
lines(x, y_pred[,"upr"], col='darkgray', type='s', lwd=4)

Compute coefficient p-values for STAR linear regression using likelihood ratio test

Description

For a linear regression model within the STAR framework, compute p-values for regression coefficients using a likelihood ratio test. It also computes a p-value for excluding all predictors, akin to a (partial) F test.

Usage

pvals(object)

Arguments

Value

a list of p+1 p-values, one for each predictor as well as the joint p-value excluding all predictors

Examples

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

# Select a transformation:
transformation = 'np'

#Estimate model
fit = lm_star(y~X, transformation = transformation)

#Compute p-values
pvals(fit)

Fit Random Forest STAR with EM algorithm

Description

Compute the MLEs and log-likelihood for the Random Forest STAR model. The STAR model requires a *transformation* and an *estimation function* for the conditional mean given observed data. The transformation can be known (e.g., log or sqrt) or unknown (Box-Cox or estimated nonparametrically) for greater flexibility. The estimator in this case is a random forest. Standard function calls including fitted and residuals apply.

Usage

randomForest_star(
  y,
  X,
  X.test = NULL,
  transformation = "np",
  y_max = Inf,
  sd_init = 10,
  tol = 10^-10,
  max_iters = 1000,
  ntree = 500,
  mtry = max(floor(ncol(X)/3), 1),
  nodesize = 5
)

Arguments

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.

The expectation-maximization (EM) algorithm is used to produce maximum likelihood estimators (MLEs) for the parameters defined in the The fitted values are computed using out-of-bag samples. As a result, the log-likelihood is based on out-of-bag prediction, and it is similarly straightforward to compute out-of-bag squared and absolute errors.

Value

a list with the following elements:

• fitted.values: the fitted values at the MLEs based on out-of-bag samples (training)

• fitted.values.test: the fitted values at the MLEs (testing)

• g.hat a function containing the (known or estimated) transformation

• sigma.hat the MLE of the standard deviation

• mu.hat the MLE of the conditional mean (on the transformed scale)

• z.hat the estimated latent data (on the transformed scale) at the MLEs

• residuals the Dunn-Smyth residuals (randomized)

• residuals_rep the Dunn-Smyth residuals (randomized) for 10 replicates

• logLik the log-likelihood at the MLEs

• logLik0 the log-likelihood at the MLEs for the *unrounded* initialization

• lambda the Box-Cox nonlinear parameter

• rfObj: the object returned by randomForest() at the MLEs

• and other parameters that (1) track the parameters across EM iterations and (2) record the model specifications

Note

Since the random forest produces random predictions, the EM algorithm will never converge exactly.

Infinite latent data values may occur when the transformed Gaussian model is highly inadequate. In that case, the function returns the *indices* of the data points with infinite latent values, which are significant outliers under the model. Deletion of these indices and re-running the model is one option, but care must be taken to ensure that (i) it is appropriate to treat these observations as outliers and (ii) the model is adequate for the remaining data points.

References

Kowal, D. R., & Wu, B. (2021). Semiparametric count data regression for self‐reported mental health. Biometrics. doi:10.1111/biom.13617

Examples

# Simulate data with count-valued response y:
sim_dat = simulate_nb_friedman(n = 100, p = 10)
y = sim_dat$y; X = sim_dat$X

# EM algorithm for STAR (using the log-link)
fit_em = randomForest_star(y = y, X = X,
                 transformation = 'log',
                 max_iters = 100)

# Fitted values (out-of-bag)
y_hat = fitted(fit_em)
plot(y_hat, y);

# Residuals:
plot(residuals(fit_em))
qqnorm(residuals(fit_em)); qqline(residuals(fit_em))

# Log-likelihood at MLEs (out-of-bag):
fit_em$logLik

Data on the efficacy of a pest management system at reducing the number of roaches in urban apartments.

Description

Data on the efficacy of a pest management system at reducing the number of roaches in urban apartments.

Usage

roaches

Format

## 'roaches' A data frame with 262 obs. of 6 variables:

y

Number of roaches caught

roach1

Pretreatment number of roaches

treatment

Treatment indicator

senior

Indicator for only elderly residents in building

exposure2

Number of days for which the roach traps were used

Source

Gelman and Hill (2007); package 'rstanarm'

Rounding function

Description

Define the rounding operator associated with the floor function. The function also returns zero whenever the input is negative and caps the value at y_max, where y_max is a known upper bound on the data y (if specified).

Usage

round_floor(z, y_max = Inf)

Arguments

Value

The count-valued output(s) from the rounding function.

Examples

# Floor function:
round_floor(1.5)
round_floor(0.5)

# Special treatmeant of negative numbers:
round_floor(-1)

Sample from a truncated normal distribution

Description

Sample from a truncated normal distribution. Samples are drawn componentwise, so each component of the vector is allowed its own mean, standard deviation, and upper and lower limits. The components are assumed to be independent.

Usage

rtruncnormRcpp(y_lower, y_upper, mu, sigma, u_rand)

Arguments

Value

z_star m x 1 draw from the truncated normal distribution

Note

This function uses Rcpp for computational efficiency.

Sample a Gaussian vector using the fast sampler of BHATTACHARYA et al.

Description

Sample from N(mu, Sigma) where Sigma = solve(crossprod(Phi) + solve(D)) and mu = Sigma*crossprod(Phi, alpha):

Usage

sampleFastGaussian(Phi, Ddiag, alpha)

Arguments

Value

Draw from N(mu, Sigma), which is p x 1, and is computed in O(n^2*p)

Note

Assumes D is diagonal, but extensions are available

Sample the parameters for an additive model

Description

Sample the parameters for an additive model, which may contain both linear and nonlinear predictors. The nonlinear terms are modeled using orthogonalized splines. The sampler draws the linear terms jointly and then samples each vector of nonlinear coefficients using Bayesian backfitting (i.e., conditional on all other nonlinear and linear terms).

Usage

sample_bam_orthog(
  y,
  X_lin,
  X_nonlin,
  params,
  A = 10^4,
  B_all = NULL,
  diagBtB_all = NULL,
  XtX = NULL
)

Arguments

Value

The updated named list params with draws from the full conditional distributions of sigma and coefficients (and updated mu).

Note

The parameters in coefficients are:

• beta_lin: the p x 1 linear coefficients, including the linear terms from X_nonlin

• f_j: the n x pNL matrix of fitted values for each nonlinear function

• theta_j: the pNL-dimensional of nonlinear basis coefficients

• sigma_beta: p x 1 vector of linear regression coefficient standard deviations

• sigma_theta_j: pNL x 1 vector of nonlinear coefficient standard deviations

Sample the parameters for an additive model

Description

Sample the parameters for an additive model, which may contain both linear and nonlinear predictors. The nonlinear terms are modeled using low-rank thin plate splines. The sampler draws the linear terms jointly and then samples each vector of nonlinear coefficients using Bayesian backfitting (i.e., conditional on all other nonlinear and linear terms).

Usage

sample_bam_thin(
  y,
  X_lin,
  X_nonlin,
  params,
  A = 10^4,
  B_all = NULL,
  BtB_all = NULL,
  XtX = NULL
)

Arguments

Value

The updated named list params with draws from the full conditional distributions of sigma and coefficients (and updated mu).

Note

The parameters in coefficients are:

• beta_lin: the p x 1 linear coefficients, including the linear terms from X_nonlin

• f_j: the n x pNL matrix of fitted values for each nonlinear function

• theta_j: the pNL-dimensional of nonlinear basis coefficients

• sigma_beta: p x 1 vector of linear regression coefficient standard deviations

• sigma_theta_j: pNL x 1 vector of nonlinear coefficient standard deviations

Sample the linear regression parameters assuming a g-prior

Description

Sample the parameters for a linear regression model assuming a g-prior for the coefficients.

Usage

sample_lm_gprior(y, X, params, psi = NULL, XtX = NULL, X_test = NULL)

Arguments

Value

The updated named list params with draws from the full conditional distributions of sigma and coefficients (along with updated mu and mu_test if applicable).

Note

The parameters in coefficients are:

• beta: the p x 1 vector of regression coefficients components of beta

Examples

# Simulate data for illustration:
sim_dat = simulate_nb_lm(n = 100, p = 5)
y = sim_dat$y; X = sim_dat$X
# Initialize:
params = init_lm_gprior(y = y, X = X)
# Sample:
params = sample_lm_gprior(y = y, X = X, params = params)
names(params)
names(params$coefficients)

Sample linear regression parameters assuming horseshoe prior

Description

Sample the parameters for a linear regression model assuming a horseshoe prior for the (non-intercept) coefficients. The number of predictors p may exceed the number of observations n.

Usage

sample_lm_hs(y, X, params, XtX = NULL, X_test = NULL)

Arguments

Value

The updated named list params with draws from the full conditional distributions of sigma and coefficients (along with updated mu and mu_test if applicable).

Note

The parameters in coefficients are:

• beta the p x 1 vector of regression coefficients

• sigma_beta p x 1 vector of regression coefficient standard deviations (local scale parameters)

• xi_sigma_beta p x 1 vector of parameter-expansion variables for sigma_beta

• lambda_beta the global scale parameter

• xi_lambda_beta parameter-expansion variable for lambda_beta components of beta

Sample linear regression parameters assuming a ridge prior

Description

Sample the parameters for a linear regression model assuming a ridge prior for the (non-intercept) coefficients. The number of predictors p may exceed the number of observations n.

Usage

sample_lm_ridge(y, X, params, A = 10^4, XtX = NULL, X_test = NULL)

Arguments

Value

The updated named list params with draws from the full conditional distributions of sigma and coefficients (along with updated mu and mu_test if applicable).

Note

The parameters in coefficients are:

• beta: the p x 1 vector of regression coefficients

• sigma_beta: the prior standard deviation for the (non-intercept) components of beta

Sample the parameters for a simple mean-only model

Description

Sample the parameters for the model y ~ N(mu0, sigma^2) with a flat prior on mu0 and sigma ~ Unif(0, A).

Usage

sample_params_mean(y, params)

Arguments

Value

The updated named list params with draws from the full conditional distributions of sigma and coefficients (and updated mu).

Note

The only parameter in coefficients is mu0. Although redundant here, this parametrization is useful in other functions.

Compute Simultaneous Band Scores (SimBaS)

Description

Compute simultaneous band scores (SimBaS) from Meyer et al. (2015, Biometrics). SimBaS uses MC(MC) simulations of a function of interest to compute the minimum alpha such that the joint credible bands at the alpha level do not include zero. This quantity is computed for each grid point (or observation point) in the domain of the function.

Usage

simBaS(sampFuns)

Arguments

Value

m x 1 vector of simBaS

Note

The input needs not be curves: the simBaS may be computed for vectors to achieve a multiplicity adjustment.

The minimum of the returned value, PsimBaS_t, over the domain t is the Global Bayesian P-Value (GBPV) for testing whether the function is zero everywhere.

Simulate count data from Friedman's nonlinear regression

Description

Simulate data from a negative-binomial distribution with nonlinear mean function.

Usage

simulate_nb_friedman(
  n = 100,
  p = 10,
  r_nb = 1,
  b_int = log(1.5),
  b_sig = log(5),
  sigma_true = sqrt(2 * log(1)),
  seed = NULL
)

Arguments

Details

The log-expected counts are modeled using the Friedman (1991) nonlinear function with interactions, possibly with additional Gaussian noise (on the log-scale). We assume that half of the predictors are associated with the response, i.e., true signals. For sufficiently large dispersion parameter r_nb, the distribution will approximate a Poisson distribution. Here, the predictor variables are simulated from independent uniform distributions.

Value

A named list with the simulated count response y, the simulated design matrix X, and the true expected counts Ey.

Note

Specifying sigma_true = sqrt(2*log(1 + a)) implies that the expected counts are inflated by 100*a% (relative to exp(X*beta)), in addition to providing additional overdispersion.

Examples

# Simulate and plot the count data:
sim_dat = simulate_nb_friedman(n = 100, p = 10);
plot(sim_dat$y)

Simulate count data from a linear regression

Description

Simulate data from a negative-binomial distribution with linear mean function.

Usage

simulate_nb_lm(
  n = 100,
  p = 10,
  r_nb = 1,
  b_int = log(1.5),
  b_sig = log(2),
  sigma_true = sqrt(2 * log(1)),
  ar1 = 0,
  intercept = FALSE,
  seed = NULL
)

Arguments

Details

The log-expected counts are modeled as a linear function of covariates, possibly with additional Gaussian noise (on the log-scale). We assume that half of the predictors are associated with the response, i.e., true signals. For sufficiently large dispersion parameter r_nb, the distribution will approximate a Poisson distribution. Here, the predictor variables are simulated from independent standard normal distributions.

Value

A named list with the simulated count response y, the simulated design matrix X (possibly including the intercept), the true expected counts Ey, and the true regression coefficients beta_true.

Note

Specifying sigma_true = sqrt(2*log(1 + a)) implies that the expected counts are inflated by 100*a% (relative to exp(X*beta)), in addition to providing additional overdispersion.

Examples

# Simulate and plot the count data:
sim_dat = simulate_nb_lm(n = 100, p = 10);
plot(sim_dat$y)

Initialize and reparametrize a spline basis matrix

Description

Following Wand and Ormerod (2008), compute a low-rank thin plate spline basis which is diagonalized such that the prior variance for the nonlinear component is a scalar times a diagonal matrix. Knot locations are determined by quantiles and the penalty is the integrated squared second derivative.

Usage

splineBasis(tau, sumToZero = TRUE, rescale01 = TRUE)

Arguments

Value

B_nl: the nonlinear component of the spline basis matrix

Note

To form the full spline basis matrix, compute cbind(1, tau, B_nl). The sum-to-zero constraint implicitly assumes that the linear term is centered and scaled, i.e., scale(tau).

Estimation for Bayesian STAR spline regression

Description

Compute samples from the predictive distributions of a STAR spline regression model using either a Gibbs sampling approach or exact Monte Carlo sampling (default is Gibbs sampling which scales better for large n)

Usage

spline_star(
  y,
  tau = NULL,
  transformation = "np",
  y_max = Inf,
  psi = NULL,
  approx_Fz = FALSE,
  approx_Fy = FALSE,
  nsave = 1000,
  use_MCMC = TRUE,
  nburn = 1000,
  nskip = 0,
  verbose = TRUE,
  method_sigma = "mle"
)

Arguments

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. Here, the continuous latent data model is a spline regression.

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'. 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. The distribution-based transformations approximately preserve the mean and variance of the count data y on the latent data scale, which lends interpretability to the model parameters. Lastly, the transformation can be modeled using the Bayesian bootstrap ('bnp'), which is a Bayesian nonparametric model and incorporates the uncertainty about the transformation into posterior and predictive inference.

Value

a list with the following elements:

• post_ytilde: nsave x n samples from the posterior predictive distribution at the observation points tau

• marg_like: the marginal likelihood (only if use_MCMC=FALSE; otherwise NULL)

Note

For the 'bnp' transformation (without the Fy approximation), there are numerical stability issues when psi is modeled as unknown. In this case, it is better to fix psi at some positive number.

Examples

# Simulate some data:
n = 100
tau = seq(0,1, length.out = n)
y = round_floor(exp(1 + rnorm(n)/4 + poly(tau, 4)%*%rnorm(n=4, sd = 4:1)))

# Sample from the predictive distribution of a STAR spline model:
fit = spline_star(y = y, tau = tau)

# Compute 90% prediction intervals:
pi_y = t(apply(fit$post_ytilde, 2, quantile, c(0.05, .95)))

# Plot the results: intervals, median, and smoothed mean
plot(tau, y, ylim = range(pi_y, y))
polygon(c(tau, rev(tau)),c(pi_y[,2], rev(pi_y[,1])),col='gray', border=NA)
lines(tau, apply(fit$post_ytilde, 2, median), lwd=5, col ='black')
lines(tau, smooth.spline(tau, apply(fit$post_ytilde, 2, mean))$y, lwd=5, col='blue')
lines(tau, y, type='p')

Monte Carlo predictive sampler for spline regression

Description

Compute direct Monte Carlo samples from the posterior predictive distribution of a STAR spline regression model.

Usage

spline_star_exact(
  y,
  tau = NULL,
  transformation = "np",
  y_max = Inf,
  psi = 1000,
  method_sigma = "mle",
  approx_Fz = FALSE,
  approx_Fy = FALSE,
  nsave = 1000,
  compute_marg = TRUE
)

Arguments

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. Here, the continuous latent data model is a spline regression.

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'. 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. The distribution-based transformations approximately preserve the mean and variance of the count data y on the latent data scale, which lends interpretability to the model parameters. Lastly, the transformation can be modeled using the Bayesian bootstrap ('bnp'), which is a Bayesian nonparametric model and incorporates the uncertainty about the transformation into posterior and predictive inference.

The Monte Carlo sampler produces direct, discrete, and joint draws from the posterior predictive distribution of the spline regression model at the observed tau points.

Value

a list with the following elements:

• post_ytilde: nsave x n samples from the posterior predictive distribution at the observation points tau

• marg_like: the marginal likelihood (if requested; otherwise NULL)

Compute the first and second moment of a truncated normal

Description

Given lower and upper endpoints and the mean and standard deviation of a (non-truncated) normal distribution, compute the first and second moment of the truncated normal distribution. All inputs may be scalars or vectors.

Usage

truncnorm_mom(a, b, mu, sig)

Arguments

Value

a list containing the first moment m1 and the second moment m2

Univariate Slice Sampler from Neal (2008)

Description

Compute a draw from a univariate distribution using the code provided by Radford M. Neal. The documentation below is also reproduced from Neal (2008).

Usage

uni.slice(x0, g, w = 1, m = Inf, lower = -Inf, upper = +Inf, gx0 = NULL)

Arguments

Value

The point sampled, with its log density attached as an attribute.

Note

The log density function may return -Inf for points outside the support of the distribution. If a lower and/or upper bound is specified for the support, the log density function will not be called outside such limits.

Update parameters for warpDLM model with trend DLM

Description

This function serves to update the warpDLM variance parameters when the underlying DLM is a structural model (i.e. local level or local linear trend). It assumes a Unif(0,A=10^4) prior on all standard deviations.

Usage

update_struct(fit, z_star, theta)

Arguments

Value

A KFAS model object (of class SSModel) updated with the newly sampled variance parameters

Posterior Inference for warpDLM model with latent structural DLM

Description

This function outputs posterior quantities and forecasts from a univariate warpDLM model. Currently two latent DLM specifications are supported: local level and the local linear trend.

Usage

warpDLM(
  y,
  type = c("level", "trend"),
  transformation = c("np", "identity", "log", "sqrt", "pois", "neg-bin"),
  y_max = Inf,
  R0 = 10,
  nsave = 5000,
  nburn = 5000,
  nskip = 1,
  n.ahead = 1
)

Arguments

Value

A list with the following elements:

• V_post: posterior draws of the observation variance

• W_post: posterior draws of the state update variance(s)

• fc_post: draws from the forecast distribution (of length n.ahead)

• post_pred: draws from the posterior predictive distribution of y

• g_func: transformation function

• g_inv_func: inverse transformation function

• KFAS_mod: the final KFAS model representing the latent DLM