Benchmark Examples for UniCensorEM

Introduction

This vignette demonstrates the use of UniCensorEM with several benchmark distributions. The package is designed to work with any user-supplied probability density function (PDF), cumulative distribution function (CDF), and survival function.

Exponential Distribution

The exponential distribution is a simple one-parameter distribution often used for lifetime data.

Complete Data

# Define exponential distribution functions
pdf_exp <- function(x, theta) theta[1] * exp(-theta[1] * x)
cdf_exp <- function(x, theta) 1 - exp(-theta[1] * x)
survival_exp <- function(x, theta) exp(-theta[1] * x)

# Generate complete data
set.seed(123)
x_exp <- rexp(100, rate = 2)

# Fit using EM algorithm
fit_exp <- em_fit(
  data = x_exp,
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

# View results
print(fit_exp)

## EM Algorithm Fit
## ================
## Censoring scheme: complete 
## Iterations: 2 
## Convergence: Success 
## Log-likelihood: -35.1557 
## Elapsed time: 0.03 seconds
## 
## Parameter estimates:
## [1] 1.91256
## 
## Standard errors:
## [1] 0.191256

summary(fit_exp)

## EM Algorithm Fit Summary
## =========================
## Censoring scheme: complete 
## Iterations: 2 
## Convergence: Success 
## Message: Converged 
## Log-likelihood: -35.1557 
## Number of parameters: 1 
## Elapsed time: 0.03 seconds
## 
## Parameter estimates:
##   Estimate Std..Error     2.5%    97.5%
## 1  1.91256   0.191256 1.537705 2.287415
## 
## Information Criteria:
##       AIC      AICc       BIC      HQIC      CAIC       KIC       TIC    LogLik 
##  72.31145  72.35226  74.91662  73.36581  75.91662  73.31145  71.31145 -35.15572 
##      Npar      Nobs 
##   1.00000 100.00000

Right Censored Data

# Create right-censored data
set.seed(123)
x_exp_rc <- rexp(100, rate = 2)
status_exp_rc <- rep(1L, 100)
status_exp_rc[1:20] <- 0L  # Right-censor first 20 observations

# Fit with right censoring
fit_exp_rc <- em_fit(
  data = data.frame(x = x_exp_rc, status = status_exp_rc),
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "right_censor",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_exp_rc)

## EM Algorithm Fit
## ================
## Censoring scheme: right_censor 
## Iterations: 11 
## Convergence: Success 
## Log-likelihood: -45.9761 
## Elapsed time: 0.39 seconds
## 
## Parameter estimates:
## [1] 1.530048
## 
## Standard errors:
## [1] 0.171065

Weibull Distribution

The Weibull distribution is a flexible two-parameter distribution commonly used in reliability analysis.

# Define Weibull distribution functions
pdf_weibull <- function(x, theta) {
  theta[1] * theta[2] * (x^(theta[2] - 1)) * exp(-theta[1] * x^theta[2])
}
cdf_weibull <- function(x, theta) 1 - exp(-theta[1] * x^theta[2])
survival_weibull <- function(x, theta) exp(-theta[1] * x^theta[2])

# Generate Weibull data
set.seed(123)
x_weibull <- rweibull(100, shape = 2, scale = 1)

# Fit using EM algorithm
fit_weibull <- em_fit(
  data = x_weibull,
  pdf = pdf_weibull,
  cdf = cdf_weibull,
  survival = survival_weibull,
  theta0 = c(1.0, 2.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_weibull)

## EM Algorithm Fit
## ================
## Censoring scheme: complete 
## Iterations: 2 
## Convergence: Success 
## Log-likelihood: -57.5532 
## Elapsed time: 0.04 seconds
## 
## Parameter estimates:
## [1] 1.000787 2.040112
## 
## Standard errors:
## [1] 0.105432 0.156337

Gamma Distribution

The gamma distribution is a two-parameter family of continuous probability distributions.

# Define gamma distribution functions
pdf_gamma <- function(x, theta) {
  (theta[1]^theta[2] / gamma(theta[2])) * x^(theta[2] - 1) * exp(-theta[1] * x)
}
cdf_gamma <- function(x, theta) pgamma(x, shape = theta[2], rate = theta[1])
survival_gamma <- function(x, theta) pgamma(x, shape = theta[2], rate = theta[1], lower.tail = FALSE)

# Generate gamma data
set.seed(123)
x_gamma <- rgamma(100, shape = 2, rate = 1)

# Fit using EM algorithm
fit_gamma <- em_fit(
  data = x_gamma,
  pdf = pdf_gamma,
  cdf = cdf_gamma,
  survival = survival_gamma,
  theta0 = c(1.0, 2.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_gamma)

## EM Algorithm Fit
## ================
## Censoring scheme: complete 
## Iterations: 5 
## Convergence: Success 
## Log-likelihood: -139.9499 
## Elapsed time: 0.06 seconds
## 
## Parameter estimates:
## [1] 1.274009 2.193205
## 
## Standard errors:
## [1] 0.188989 0.289684

Log-Normal Distribution

The log-normal distribution is useful for modeling variables that are positively skewed.

# Define log-normal distribution functions
pdf_lnorm <- function(x, theta) {
  (1 / (x * theta[2] * sqrt(2 * pi))) * exp(-(log(x) - theta[1])^2 / (2 * theta[2]^2))
}
cdf_lnorm <- function(x, theta) plnorm(x, meanlog = theta[1], sdlog = theta[2])
survival_lnorm <- function(x, theta) plnorm(x, meanlog = theta[1], sdlog = theta[2], lower.tail = FALSE)

# Generate log-normal data
set.seed(123)
x_lnorm <- rlnorm(100, meanlog = 0, sdlog = 1)

# Fit using EM algorithm
fit_lnorm <- em_fit(
  data = x_lnorm,
  pdf = pdf_lnorm,
  cdf = cdf_lnorm,
  survival = survival_lnorm,
  theta0 = c(0.0, 1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_lnorm)

## EM Algorithm Fit
## ================
## Censoring scheme: complete 
## Iterations: 9 
## Convergence: Success 
## Log-likelihood: -141.3098 
## Elapsed time: 0.07 seconds
## 
## Parameter estimates:
## [1] 0.09040669 0.90824055
## 
## Standard errors:
## [1] 0.090824 0.064222

Normal Distribution

The normal (Gaussian) distribution is the most common continuous probability distribution.

# Define normal distribution functions
pdf_normal <- function(x, theta) {
  (1 / (theta[2] * sqrt(2 * pi))) * exp(-(x - theta[1])^2 / (2 * theta[2]^2))
}
cdf_normal <- function(x, theta) pnorm(x, mean = theta[1], sd = theta[2])
survival_normal <- function(x, theta) pnorm(x, mean = theta[1], sd = theta[2], lower.tail = FALSE)

# Generate normal data
set.seed(123)
x_normal <- rnorm(100, mean = 0, sd = 1)

# Fit using EM algorithm
fit_normal <- em_fit(
  data = x_normal,
  pdf = pdf_normal,
  cdf = cdf_normal,
  survival = survival_normal,
  theta0 = c(0.0, 1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_normal)

## EM Algorithm Fit
## ================
## Censoring scheme: complete 
## Iterations: 9 
## Convergence: Success 
## Log-likelihood: -132.2692 
## Elapsed time: 0.07 seconds
## 
## Parameter estimates:
## [1] 0.09040669 0.90824055
## 
## Standard errors:
## [1] 0.090824 0.064222

Model Selection

Compare different models using information criteria.

# Fit exponential to Weibull data
fit_exp_on_weibull <- em_fit(
  data = x_weibull,
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "complete",
  control = em.control(maxit = 500, tol = 1e-8)
)

# Compare information criteria
cat("Exponential fit on Weibull data:\n")

## Exponential fit on Weibull data:

print(model_selection(fit_exp_on_weibull))

##      AIC     AICc      BIC     HQIC     CAIC      KIC      TIC   LogLik 
## 177.6358 177.6766 180.2410 178.6902 181.2410 178.6358 176.6358 -87.8179 
##     Npar     Nobs 
##   1.0000 100.0000

cat("\nWeibull fit on Weibull data:\n")

## 
## Weibull fit on Weibull data:

print(model_selection(fit_weibull))

##       AIC      AICc       BIC      HQIC      CAIC       KIC       TIC    LogLik 
## 119.10648 119.23019 124.31682 121.21519 126.31682 121.10648 117.10648 -57.55324 
##      Npar      Nobs 
##   2.00000 100.00000

Goodness-of-Fit Statistics

Compute goodness-of-fit statistics for the fitted model.

gof <- gof_stats(fit_weibull)
print(gof)

##               KS Anderson_Darling Cramer_von_Mises           Watson 
##     0.0403382791     0.1764043503     0.0227915927     0.0206019769 
##             RMSE              MAE              MSE        R_squared 
##     0.0150968847     0.0119457630     0.0002279159     0.9972647353 
##    Adj_R_squared 
##     0.9972083381

Bootstrap Standard Errors

Compute bootstrap standard errors for parameter estimates.

# Bootstrap with fewer iterations for speed
boot_results <- bootstrap_se(
  fit_exp,
  B = 100,
  type = "parametric",
  ci_type = "percentile",
  level = 0.95
)

## Bootstrap iteration: 100 / 100

print(boot_results)

## Bootstrap Results
## =================
## Bootstrap type: parametric 
## Number of replicates: 100 
## Valid replicates: 0 
## Confidence interval type: percentile 
## Confidence level: 95 %
## 
## Bootstrap standard errors:
## [1] NA
## 
## Bootstrap bias:
## [1] NaN
## 
## Bootstrap confidence intervals:
##      2.5% 97.5%
## [1,]   NA    NA

Diagnostic Plots

Create diagnostic plots for the EM algorithm fit.

plot(fit_weibull, which = c(1, 2, 3), ask = FALSE)

Residual Analysis

Perform residual analysis to check model fit.

plot_residuals(fit_weibull, type = "cox-snell", which = c(1, 2, 3))

EM Diagnostics

Examine EM algorithm convergence diagnostics.

diag <- em_diagnostics(fit_weibull)
print(diag)

## EM Algorithm Diagnostics
## =========================
## Iterations: 2 
## Convergence status: Success 
## Convergence message: Converged 
## 
## Parameter changes:
##   Maximum change: 0.000786, 0.040112 
##   Final change: 1e-06, 0 
##   Maximum relative change: 0.000786, 0.020056 
##   Final relative change: 1e-06, 0 
## 
## Likelihood changes:
##   Maximum change: 0.037387 
##   Final change: 0 
## 
## Q-function changes:
##   Maximum change: 0 
##   Final change: 0 
## 
## Tolerance checks:
##   Absolute tolerance met: FALSE 
##   Relative tolerance met: FALSE 
##   EM tolerance met: TRUE

Interval Censoring Example

Example with interval-censored data.

# Create interval-censored data
set.seed(123)
x_interval <- rexp(100, rate = 2)
status_interval <- rep(1L, 100)
time2_interval <- rep(NA, 100)

# Make some observations interval-censored
interval_idx <- sample(1:100, 20)
status_interval[interval_idx] <- 3L
time2_interval[interval_idx] <- x_interval[interval_idx] + runif(20, 0.5, 1.5)

# Fit with interval censoring
fit_interval <- em_fit(
  data = data.frame(x = x_interval, status = status_interval, time2 = time2_interval),
  pdf = pdf_exp,
  cdf = cdf_exp,
  survival = survival_exp,
  theta0 = c(1.0),
  scheme = "interval_censor",
  control = em.control(maxit = 500, tol = 1e-8)
)

print(fit_interval)

## EM Algorithm Fit
## ================
## Censoring scheme: interval_censor 
## Iterations: 5 
## Convergence: Success 
## Log-likelihood: -22.3219 
## Elapsed time: 0.09 seconds
## 
## Parameter estimates:
## [1] 1.002217
## 
## Standard errors:
## [1] 0.233254

Conclusion

This vignette demonstrates the flexibility of UniCensorEM in handling various distributions and censoring schemes. The package can be extended to any user-defined distribution by providing the appropriate PDF, CDF, and survival functions.