A comprehensive R package for finding Maximum Likelihood Estimation (MLE) for any univariate distribution under various censoring and truncation schemes.
MleCensoR provides generalized functions to compute MLE under the following censoring schemes:
mle_right_truncation)mle_left_truncation)mle_random)mle_right)mle_left)mle_interval)mle_block_random)mle_bjpt2)mle_progressive_first_failure)mle_joint_type1)mle_type1)mle_type2)mle_progressive_type2)mle_joint_type2)mle_hybrid, mle_hybrid_type1,
mle_type1_hybrid)mle_hybrid_type2)mle_doubly_type2)mle_middle)mle_progressive_hybrid_type2)# Install from local directory
install.packages("path/to/MleCensoR")
# Or install using devtools
devtools::install("path/to/MleCensoR")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)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"
)The package supports the following optimization 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"
)# 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
)
)# 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))
)# 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"
)# 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
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 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)
)# 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)
)The package provides S3 methods for extracting results:
coef(fit) - Parameter estimateslogLik(fit) - Log-likelihood valuevcov(fit) - Variance-covariance matrixstdEr(fit) - Standard errorssummary(fit) - Summary table with estimates, standard
errors, z-values, and p-valuesAIC(fit) - Akaike Information CriterionnIter(fit) - Number of iterationsAdvanced 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
)GPL-3
Shikhar Tyagi
Maintainer: Shikhar Tyagi shikhar.tyagi@example.com