Package 'LugsailGR'

Calculate Minimum Required Effective Sample Size (M_alpha_eps_p)

Description

Computes the theoretical minimum effective sample size required to obtain a confidence region with relative volume epsilon and coverage 1 - alpha, based on Vats, Flegal and Jones (2019) and Vats and Knudson (2021).

Usage

calc_ess_bound(p = 1, alpha = 0.05, epsilon = 0.1)

Arguments

p	Integer, number of parameters.
alpha	Numeric, significance level (default 0.05 for 95% confidence).
epsilon	Numeric, relative volume tolerance (default 0.10).

Value

Numeric scalar, required minimum ESS.

Examples

calc_ess_bound(p = 1, alpha = 0.05, epsilon = 0.10)
calc_ess_bound(p = 5, alpha = 0.05, epsilon = 0.05)

---

Calculate Gelman-Rubin Termination Threshold (delta_epsilon)

Description

Computes the principled Gelman-Rubin diagnostic termination threshold delta_epsilon based on the target ESS bound M_alpha_eps_p and number of chains m.

Usage

calc_psrf_cutoff(m = 1, M_val = 1537)

Arguments

m	Integer, number of chains.
M_val	Numeric, minimum required ESS obtained from calc_ess_bound.

Value

Numeric scalar, target threshold delta_epsilon.

Examples

M_bound <- calc_ess_bound(p = 1, alpha = 0.05, epsilon = 0.10)
calc_psrf_cutoff(m = 3, M_val = M_bound)

---

Gelman-Rubin Convergence Diagnostics for Censored Data Models

Description

Evaluates lugsail Gelman-Rubin convergence diagnostics, effective sample size, and parameter estimates for MCMC chains fitted to censored data models under various censoring schemes (right, left, interval, Type-I, Type-II, progressive, and hybrid).

Usage

censor_gr_test(
  pdf_fn,
  cdf_fn,
  surv_fn = NULL,
  data,
  start_par,
  censor_type = c("right", "left", "interval", "type1", "type2", "progressive", "hybrid"),
  status = NULL,
  censor_info = NULL,
  n_iter = 2000,
  n_chains = 3,
  burn_in = 0.5,
  scale = 0.5,
  alpha = 0.05,
  epsilon = 0.1
)

Arguments

pdf_fn	Function 'function(x, par)' returning probability density values for variable 'x' given parameters 'par'.
cdf_fn	Function 'function(x, par)' returning cumulative distribution function values.
surv_fn	Optional function 'function(x, par)' returning survival function values. If 'NULL', computed as '1 - cdf_fn(x, par)'.
data	Numeric vector or matrix of observed sample data.
start_par	Initial parameter vector.
censor_type	Character string specifying censoring scheme. Options include: "right" (Right censoring), "left" (Left censoring), "interval" (Interval censoring), "type1" (Type-I censoring), "type2" (Type-II censoring), "progressive" (Progressive Type-II censoring), or "hybrid" (Hybrid censoring).
status	Numeric vector or matrix indicating event/censoring status (e.g. 1 = observed, 0 = censored).
censor_info	Optional list specifying parameters for censoring schemes (e.g. cutoff time T_censor, number of failures r_censor, or removal pattern R_pattern).
n_iter	Integer, total number of MCMC iterations per chain (default 2000).
n_chains	Integer, number of parallel chains (default 3).
burn_in	Numeric, burn-in fraction (default 0.50).
scale	Numeric or vector, proposal scale for random walk Metropolis-Hastings (default 0.5).
alpha	Numeric, significance level (default 0.05).
epsilon	Numeric, relative volume tolerance (default 0.10).

Value

An object of class "lugsail_gr" containing Gelman-Rubin diagnostics, effective sample sizes, cutoff thresholds, parameter estimates, and MCMC samples.

References

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

Examples

# Example: Right-censored exponential distribution
set.seed(123)
n_obs <- 30
true_rate <- 0.5
times <- rexp(n_obs, rate = true_rate)
censor_time <- 2.0
obs_times <- pmin(times, censor_time)
status <- as.numeric(times <= censor_time)

f_exp <- function(x, rate) dexp(x, rate = rate)
F_exp <- function(x, rate) pexp(x, rate = rate)
S_exp <- function(x, rate) 1 - pexp(x, rate = rate)

fit <- censor_gr_test(pdf_fn = f_exp, cdf_fn = F_exp, surv_fn = S_exp,
                     data = obs_times, start_par = c(rate = 0.8),
                     censor_type = "right", status = status,
                     n_iter = 1000, n_chains = 3)
print(fit)

---

Upgraded Univariate and Multivariate Gelman-Rubin Diagnostic

Description

Computes the upgraded Gelman-Rubin convergence statistics ($\hat{R}_L$ and $\hat{R}^p_L$), Effective Sample Size ($\widehat{\mathrm{ESS}}$), and principled termination thresholds ($\delta_{\epsilon}$) proposed by Vats and Knudson (2021) using replicated lugsail batch means variance estimators.

Usage

lugsail_gr(x, alpha = 0.05, epsilon = 0.1, b = NULL, multivariate = TRUE)

Arguments

x	Input MCMC chain output. Can be a numeric vector (single univariate chain), matrix of dimension (n, p) (single multivariate chain), 3D array of dimension (n, p, m) (multiple multivariate chains), or a list of length m containing matrices or data.frames.
alpha	Numeric, significance level for the confidence region (default 0.05 for 95% confidence).
epsilon	Numeric, relative volume tolerance for target precision (default 0.10).
b	Optional batch size for lugsail batch means. If NULL, defaults to floor(sqrt(n)).
multivariate	Logical, if TRUE (default), computes multivariate Gelman-Rubin diagnostic $\hat{R}^p_L$ and multivariate ESS when p > 1.

Value

An object of class "lugsail_gr" containing:

psrf	Vector of univariate lugsail potential scale reduction factors $\hat{R}_L$ for each parameter.
mpsrf	Multivariate lugsail potential scale reduction factor $\hat{R}^p_L$.
ess	Vector of univariate effective sample sizes for each parameter.
mess	Multivariate effective sample size $\widehat{\mathrm{ESS}}_p$.
delta_eps	Principled target threshold $\delta_{\epsilon}$.
M_alpha_eps_p	Minimum required effective sample size bound $M_{\alpha, \epsilon, p}$.
converged	Logical, indicating if diagnostic has converged (mpsrf <= delta_eps).
means	Posterior sample mean vector across chains.
sd	Posterior standard deviation vector.
s2	Sample variance (scalar or matrix S).
tau_L	Replicated lugsail variance/covariance matrix estimate $\hat{T}_L$.
n	Chain length (number of iterations per chain).
m	Number of parallel chains.
p	Number of parameter dimensions.
alpha	Significance level.
epsilon	Relative volume tolerance.

References

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472. doi:10.1214/ss/1177011136.

Brooks, S. P. and Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7(4), 434–455. doi:10.1080/10618600.1998.10474787.

Examples

# Example 1: Univariate MCMC output (3 chains of 1000 iterations)
set.seed(123)
chain1 <- rnorm(1000, mean = 0, sd = 1)
chain2 <- rnorm(1000, mean = 0.05, sd = 1)
chain3 <- rnorm(1000, mean = -0.05, sd = 1)
res1 <- lugsail_gr(list(chain1, chain2, chain3))
print(res1)

# Example 2: Multivariate MCMC output (2 chains of 500 iterations, 2 parameters)
mat1 <- matrix(rnorm(1000), ncol = 2)
mat2 <- matrix(rnorm(1000), ncol = 2)
res2 <- lugsail_gr(list(mat1, mat2))
summary(res2)

---

Evaluate Gelman-Rubin Diagnostics for User-Supplied Target Density Function

Description

Runs Markov chain Monte Carlo (MCMC) sampling for a user-supplied probability density function (PDF) or log-likelihood function of variable and parameters with data, and returns upgraded lugsail Gelman-Rubin convergence test statistics, effective sample sizes, parameter estimates, and termination decisions.

Usage

mcmc_gr_test(
  target_pdf,
  data = NULL,
  start_par,
  n_iter = 2000,
  n_chains = 3,
  burn_in = 0.5,
  scale = 0.5,
  alpha = 0.05,
  epsilon = 0.1,
  ...
)

Arguments

target_pdf	Function taking parameter vector par as its first argument and optional data. It should return the probability density value or log-density value for par.
data	Optional data vector, matrix, data.frame, or list passed to target_pdf.
start_par	Initial parameter vector. If named, parameter names will be preserved.
n_iter	Integer, total number of MCMC iterations per chain (default 2000).
n_chains	Integer, number of parallel MCMC chains to generate (default 3).
burn_in	Numeric, fraction of initial burn-in samples to discard (default 0.50).
scale	Numeric or vector, proposal standard deviation for random walk Metropolis-Hastings sampler (default 0.5).
alpha	Numeric, significance level for convergence bound (default 0.05).
epsilon	Numeric, relative volume tolerance (default 0.10).
...	Additional parameters passed to target_pdf.

Value

An object of class "lugsail_gr" containing Gelman-Rubin test statistics, effective sample sizes, cutoff thresholds, parameter estimates, standard errors, and MCMC chain samples.

References

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

Examples

# Example 1: Normal distribution target PDF with data
set.seed(42)
true_mean <- 3
my_data <- rnorm(50, mean = true_mean, sd = 1)

# Log-posterior density function for parameter mu
log_post <- function(mu, data) {
  sum(dnorm(data, mean = mu, sd = 1, log = TRUE))
}

# Run diagnostic test
fit <- mcmc_gr_test(target_pdf = log_post, data = my_data, start_par = c(mu = 0), 
                    n_iter = 1000, n_chains = 3)
print(fit)

---

Plot Method for lugsail_gr Objects

Description

Produces diagnostic plots for MCMC chains, including trace plots, running Gelman-Rubin statistic $\hat{R}_L$ plots, and density estimates.

Usage

## S3 method for class 'lugsail_gr'
plot(x, type = c("all", "trace", "running_gr", "density"), ...)

Arguments

x	An object of class "lugsail_gr".
type	Character string specifying plot type: "all" (default), "trace", "running_gr", or "density".
...	Additional graphical parameters.

Value

No return value, called for side effects (generating plots).

Examples

set.seed(123)
chain1 <- rnorm(500)
chain2 <- rnorm(500)
fit <- lugsail_gr(list(chain1, chain2))
plot(fit, type = "running_gr")

---

Print Method for lugsail_gr Objects

Description

Print Method for lugsail_gr Objects

Usage

## S3 method for class 'lugsail_gr'
print(x, ...)

Arguments

x	An object of class "lugsail_gr".
...	Further arguments passed to print.

Value

Invisibly returns x.

---

Summary Method for lugsail_gr Objects

Description

Summary Method for lugsail_gr Objects

Usage

## S3 method for class 'lugsail_gr'
summary(object, ...)

Arguments

object	An object of class "lugsail_gr".
...	Further arguments passed to summary.

Value

A summary table object.

---

Sample Titanic Survival Dataset

Description

Subset of the Titanic passenger dataset used in Vats and Knudson (2021) for demonstrating multivariate Gelman-Rubin convergence diagnostics in Bayesian logistic regression.

Usage

titanic_sub

Format

A data frame with 100 rows and 6 variables:

Survived

Binary indicator of survival (1 = survived, 0 = died)

Pclass

Passenger class (1 = 1st, 2 = 2nd, 3 = 3rd)

Sex

Sex of passenger (1 = female, 0 = male)

Age

Age in years

SibSp

Number of siblings/spouses aboard

Fare

Passenger fare

Source

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

---