MleCensoR: Maximum Likelihood Estimation under Censoring Schemes

A comprehensive R package for finding Maximum Likelihood Estimation (MLE) for any univariate distribution under various censoring and truncation schemes.

Features

MleCensoR provides generalized functions to compute MLE under the following censoring schemes:

• Right Truncation (mle_right_truncation)
• Left Truncation (mle_left_truncation)
• Random Censoring (mle_random)
• Right Censoring (mle_right)
• Left Censoring (mle_left)
• Interval Censoring (mle_interval)
• Block Random Censoring (mle_block_random)
• Balanced Joint Progressive Type-II (BJPT-II) Censoring (mle_bjpt2)
• Progressive First Failure Censoring (mle_progressive_first_failure)
• Joint Type-I Censoring (mle_joint_type1)
• Type-I Censoring (mle_type1)
• Type-II Censoring (mle_type2)
• Progressive Type-II Censoring (mle_progressive_type2)
• Joint Type-II Censoring (mle_joint_type2)
• Hybrid Censoring (Type-I Hybrid) (mle_hybrid, mle_hybrid_type1, mle_type1_hybrid)
• Hybrid Type-II Censoring (mle_hybrid_type2)
• Doubly Type-II Censoring (mle_doubly_type2)
• Middle Censoring (mle_middle)
• Type-II Progressively Hybrid Censoring (mle_progressive_hybrid_type2)

Installation

# Install from local directory
install.packages("path/to/MleCensoR")

# Or install using devtools
devtools::install("path/to/MleCensoR")

Quick Start

Basic Usage Example

library(MleCensoR)

# Define distribution functions (e.g., Weibull distribution)
pdf_weibull <- function(x, theta) {
  shape <- theta[1]
  scale <- theta[2]
  dweibull(x, shape = shape, scale = scale)
}

cdf_weibull <- function(x, theta) {
  shape <- theta[1]
  scale <- theta[2]
  pweibull(x, shape = shape, scale = scale)
}

surv_weibull <- function(x, theta) {
  1 - cdf_weibull(x, theta)
}

# Example: Right Censoring
set.seed(123)
x <- rweibull(50, shape = 2, scale = 1)
status <- rep(1, 50)  # All observed

# MLE estimation
fit <- mle_right(
  x = x,
  status = status,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 1.5, scale = 0.8),
  method = "NR"
)

# View results
summary(fit)
coef(fit)
vcov(fit)
stdEr(fit)
logLik(fit)
AIC(fit)

Using Custom Log-Likelihood

You can also provide a custom log-likelihood function directly:

# Custom log-likelihood for right censoring
custom_loglik <- function(theta, x, status) {
  shape <- theta[1]
  scale <- theta[2]
  pdf_val <- dweibull(x, shape = shape, scale = scale)
  surv_val <- pweibull(x, shape = shape, scale = scale, lower.tail = FALSE)
  ll <- status * log(pdf_val) + (1 - status) * log(surv_val)
  return(ll)
}

fit <- mle_right(
  x = x,
  status = status,
  logLik = custom_loglik,
  start = c(shape = 1.5, scale = 0.8),
  method = "NR"
)

Optimization Methods

The package supports the following optimization methods:

• NR - Newton-Raphson (default for unconstrained problems)
• BFGS - Broyden-Fletcher-Goldfarb-Shanno
• BFGSR - BFGS algorithm implemented in R
• BHHH - Berndt-Hall-Hall-Hausman (recommended for observation-wise likelihood)
• SANN - Simulated Annealing
• CG - Conjugate Gradients
• NM - Nelder-Mead (default for constrained problems without gradient)

Example with Different Methods

# Using BHHH method
fit_bhhh <- mle_right(
  x = x,
  status = status,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 1.5, scale = 0.8),
  method = "BHHH"
)

# Using BFGS method
fit_bfgs <- mle_right(
  x = x,
  status = status,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 1.5, scale = 0.8),
  method = "BFGS"
)

Advanced Features

Parameter Constraints

# Constrained optimization
fit_constrained <- mle_right(
  x = x,
  status = status,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 1.5, scale = 0.8),
  method = "BFGS",
  constraints = list(
    ineqA = matrix(c(1, 0, 0, 1), nrow = 2, byrow = TRUE),
    ineqB = c(0, 0)  # theta > 0
  )
)

Parameter Ranges

# Specify parameter ranges
fit_range <- mle_right(
  x = x,
  status = status,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 1.5, scale = 0.8),
  param_range = list(shape = c(0.1, 10), scale = c(0.1, 10))
)

Custom Gradient and Hessian

# Provide analytic gradient
custom_grad <- function(theta, x, status) {
  shape <- theta[1]
  scale <- theta[2]
  # ... compute gradient ...
  grad <- c(d_shape, d_scale)
  return(grad)
}

fit_grad <- mle_right(
  x = x,
  status = status,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 1.5, scale = 0.8),
  grad = custom_grad,
  method = "NR"
)

Censoring Scheme Examples

Interval Censoring

# Interval censoring data
l <- c(0.5, 1.0, 1.5, 2.0)  # Lower bounds
r <- c(1.0, 1.5, 2.0, Inf)  # Upper bounds

fit_interval <- mle_interval(
  l = l,
  r = r,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 2, scale = 1)
)

Progressive Type-II Censoring

# Progressive Type-II censoring
x <- c(0.5, 0.8, 1.2, 1.5, 2.0)  # Observed failure times
r_removals <- c(1, 1, 2, 1, 0)   # Removals at each failure

fit_prog <- mle_progressive_type2(
  x = x,
  r_removals = r_removals,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 2, scale = 1)
)

Joint Type-I Censoring

# Joint Type-I censoring for two populations
x <- c(0.5, 0.8, 1.2, 1.5, 2.0)  # Failure times
z <- c(1, 1, 0, 0, 1)             # Group membership (1 = pop A, 0 = pop B)
tc <- 2.0                          # Censoring time
n_A <- 3                           # Initial units in pop A
n_B <- 2                           # Initial units in pop B

fit_joint1 <- mle_joint_type1(
  x = x,
  z = z,
  tc = tc,
  n_A = n_A,
  n_B = n_B,
  pdf_A = pdf_weibull,
  cdf_A = cdf_weibull,
  surv_A = surv_weibull,
  pdf_B = pdf_weibull,
  cdf_B = cdf_weibull,
  surv_B = surv_weibull,
  start = c(shape = 2, scale = 1)
)

Type-II Progressively Hybrid Censoring

# Suppose n = 20, m = 10, tc = 2.0
# 7 failures observed before tc: Case II
x <- c(0.12, 0.35, 0.62, 0.89, 1.15, 1.48, 1.82)
r_removals <- c(1, 0, 1, 0, 1, 0, 0, 1, 0, 0) # full plan for m = 10

fit_prog_hybrid <- mle_progressive_hybrid_type2(
  x = x,
  r_removals = r_removals,
  tc = 2.0,
  n = 20,
  m = 10,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 2, scale = 1)
)

Output Methods

The package provides S3 methods for extracting results:

• coef(fit) - Parameter estimates
• logLik(fit) - Log-likelihood value
• vcov(fit) - Variance-covariance matrix
• stdEr(fit) - Standard errors
• summary(fit) - Summary table with estimates, standard errors, z-values, and p-values
• AIC(fit) - Akaike Information Criterion
• nIter(fit) - Number of iterations

Control Parameters

Advanced control parameters can be passed via ...:

• tol - Stopping tolerance (default: 1e-8)
• reltol - Relative convergence tolerance (default: sqrt(.Machine$double.eps))
• gradtol - Gradient norm tolerance (default: 1e-6)
• steptol - Step size tolerance (default: 1e-10)
• iterlim - Maximum iterations (default: 100 for NR, 200 for BFGSR)
• printLevel - Printing level (0 = none, 1 = initial/final, 2 = all iterations)

fit <- mle_right(
  x = x,
  status = status,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  surv = surv_weibull,
  start = c(shape = 1.5, scale = 0.8),
  method = "NR",
  iterlim = 500,
  tol = 1e-10,
  printLevel = 2
)

License

GPL-3

Author

Shikhar Tyagi

Contact

Maintainer: Shikhar Tyagi shikhar.tyagi@example.com