Package {exposureEM}

exposureEM: combined-exposure models for binary event histories

Description

The package fits a two-component latent-cause model in which an observed binary response is the union of two independent component-specific binary outcomes. Known exposure proportions enter both component models. The implementation is domain-neutral; sport and non-sport are one motivating application rather than part of the package interface. Users can fit the model with EM, direct Marquardt-damped Newton-Raphson, or an EM-initialized hybrid through exposure_fit().

Details

The first component is called component_0 and the second is called component_1 throughout the returned objects. Users decide what those two components mean in their own study.

Author(s)

Maintainer: Wenjing Meng w2meng@ucsd.edu

Contrast between the two exposure components

Description

Computes exp(beta_component_1 - beta_component_0) with a Wald confidence interval. Its interpretation is a rate ratio for cloglog, an odds ratio for logit, and a risk ratio for log.

Usage

component_contrast(object, level = 0.95)

Arguments

Value

A one-row data frame containing the log contrast, standard error, exponentiated estimate, and confidence limits.

Fit a two-component exposure model by EM

Description

Fits the paper's combined-exposure model. For observation i, the observed event is Y_i = I(Y_0i + Y_1i >= 1), where the two latent component events are conditionally independent. The two design rows share all covariates and have separate component-specific intercepts. Exposure proportions enter as fixed offsets.

Usage

exposure_em(
  formula,
  data,
  exposure = c("a0", "a1"),
  link = c("cloglog", "logit", "log"),
  start = NULL,
  control = exposure_em_control(),
  na.action = stats::na.omit,
  contrasts = NULL
)

Arguments

Details

The outer criterion selects a preliminary stopping check. Termination also requires an absolute observed likelihood change below tol and ⁠score' information^{-1} score / length(score) <= tol⁠, with unmodified observed information having eigenvalues above 1e-10. The maximum absolute unmodified observed Newton step must also be at most sqrt(tol) to avoid treating a vanishing boundary score as a finite root. The M-step keeps E-step expectations fixed and checks Q ascent. Formula offsets are rejected: the two exposure columns supply the model's offsets. Both response values must occur, and the active component design must have full column rank. Constant cloglog exposure proportions are rejected because the separate intercepts are not identified. The rank check alone does not establish observed-model identification in other designs. For na.exclude, in-sample prediction, fitted-value, and residual methods restore excluded rows as NA; stored arrays and posterior summaries contain only analyzed rows. Model-based covariance assumes independent observations. The component contrast compares equal exposures (or the underlying rates, odds, or risks per unit exposure), not unequal observed component exposures.

Value

An object of class exposure_em_fit. It contains coefficients, vcov, the unmodified observed information, logLik, aic, the final observed score and newton_decrement, converged, convergence, iteration counts, and an optional trace. Analyzed-row probabilities, posterior expectations, residuals, exposures, and the model frame are also stored. A singular or non-positive information matrix produces an unavailable (NA) covariance and a warning; a convergence flag does not prove a global maximum or exclude weak identification.

Examples

set.seed(42)
n <- 500
a1 <- c(rep(0, 100), stats::rbeta(400, 2, 12))
a0 <- 1 - a1
age <- stats::rnorm(n)
lambda0 <- a0 * exp(-3.2 + 0.15 * age)
lambda1 <- a1 * exp(-1.4 + 0.15 * age)
y <- stats::rbinom(n, 1, -expm1(-(lambda0 + lambda1)))
dat <- data.frame(y, a0, a1, age)
fit <- exposure_em(y ~ age, dat, exposure = c("a0", "a1"))
coef(fit)
component_contrast(fit)

Control parameters for exposure-component EM fitting

Description

Control parameters for exposure-component EM fitting

Usage

exposure_em_control(
  tol = 1e-08,
  maxit = 10000L,
  mstep_tol = 1e-08,
  mstep_maxit = 1000L,
  criterion = c("parameter", "loglik", "both"),
  step_factor = 0.5,
  min_step = 2^-20,
  ridge_initial = 0,
  ridge_multiplier = 10,
  ridge_max = 1e+08,
  exposure_tol = sqrt(.Machine$double.eps),
  check_exposure_sum = TRUE,
  keep_trace = TRUE,
  trace = FALSE
)

Arguments

Value

A named list consumed by exposure_em().

Fit a two-component exposure model using EM, Marquardt, or both

Description

Fits the same two-component combined-exposure likelihood through one of three computational paths. method = "em" uses Bernoulli-indicator EM. method = "marquardt" directly maximizes the observed-data likelihood with Marquardt-stabilized Newton-Raphson. method = "hybrid" first performs a limited number of EM outer iterations and then passes those coefficients to the direct optimizer.

Usage

exposure_fit(
  formula,
  data,
  exposure = c("a0", "a1"),
  link = c("cloglog", "logit", "log"),
  method = c("em", "marquardt", "hybrid"),
  start = NULL,
  hybrid_em_iterations = 10L,
  em_control = exposure_em_control(),
  marquardt_control = exposure_marquardt_control(),
  na.action = stats::na.omit,
  contrasts = NULL
)

Arguments

Details

All three methods support complementary log-log, logit, and log component links. Hybrid fitting uses the selected link in both phases.

Direct termination requires the sum of squared coefficient changes, absolute likelihood change, and the unmodified Newton decrement per coefficient all to be at most tol. At an exactly zero score the proposed change is zero. Diagonal inflation uses 0.99 * abs(diag(H)) + 0.01 * abs(sum(diag(H))), with a machine-precision floor for flat diagonals. Backtracking enforces Armijo sufficient increase up to a rounding allowance of 1e-12 * (1 + abs(logLik)). Final covariance never uses inflation. Log-link steps must keep every component probability below one; invalid proposals are rejected during backtracking. For logit, proposals whose component probabilities round to one are also rejected numerically. These are interior optimizers: reaching a probability boundary does not certify a regular finite maximum or valid Wald inference. Hybrid fitting switches after at most hybrid_em_iterations, or earlier if EM terminates. Each phase respects its own keep_trace setting; the combined trace contains both phases when both are retained.

Value

An object inheriting from exposure_fit and exposure_em_fit, with the fields described in exposure_em(). Direct and hybrid fits additionally retain em_trace and marquardt_trace when requested. iterations is the sum of completed EM and accepted direct updates; em_iterations and marquardt_iterations give the separate counts.

Examples

set.seed(12)
n <- 500
a1 <- ifelse(stats::runif(n) < 0.2, 0, stats::rbeta(n, 2, 15))
a0 <- 1 - a1
x <- stats::rnorm(n)
p <- -expm1(-(a0 * exp(-3.2 + 0.2 * x) +
               a1 * exp(-1.5 + 0.2 * x)))
dat <- data.frame(y = stats::rbinom(n, 1, p), a0, a1, x)

fit_em <- exposure_fit(y ~ x, dat, method = "em")
fit_direct <- exposure_fit(y ~ x, dat, method = "marquardt")
fit_hybrid <- exposure_fit(
  y ~ x, dat, method = "hybrid", hybrid_em_iterations = 5
)
component_contrast(fit_hybrid)
fit_logit <- exposure_fit(y ~ x, dat, link = "logit", method = "marquardt")
fit_log <- exposure_fit(y ~ x, dat, link = "log", method = "hybrid")

Control parameters for direct Marquardt fitting

Description

These settings control direct maximization of the observed-data likelihood. They are used by method = "marquardt" and by the direct phase of method = "hybrid" in exposure_fit().

Usage

exposure_marquardt_control(
  tol = 1e-08,
  maxit = 200L,
  step_factor = 0.5,
  min_step = 1e-04,
  armijo = 1e-04,
  inflation_initial = 0.01,
  inflation_multiplier = 10,
  inflation_attempts = 20L,
  pd_tol = 1e-10,
  keep_trace = TRUE,
  trace = FALSE
)

Arguments

Value

A named list consumed by exposure_fit().

Posterior latent-component probabilities

Description

Returns the converged E-step quantities. component_0 and component_1 are the posterior expectations of the two latent indicators; because both events can occur, they are not mutually exclusive and can sum to more than one.

Usage

posterior_components(object, joint = TRUE)

Arguments

Value

A data frame with one row per analyzed observation.

Predict from a fitted exposure model

Description

Predict from a fitted exposure model

Usage

## S3 method for class 'exposure_em_fit'
predict(object, newdata = NULL, type = c("response", "component", "link"), ...)

Arguments

Value

A numeric vector for type = "response", otherwise a two-column matrix.