Combined-Exposure Models by EM and Marquardt Optimization

Fits general two-component combined-exposure models for binary event histories when the event setting is not observed. The observed binary event is represented as the union of two latent component-specific binary events. Known exposure proportions enter as offsets. Parameters can be estimated by expectation-maximization, direct Marquardt-damped Newton-Raphson maximization of the observed likelihood, or a hybrid that uses several expectation-maximization iterations before direct optimization. Uncertainty is estimated with Louis' formula for the expectation-maximization estimator and the inverse observed Hessian for direct and hybrid fits. Complementary log-log, logit, and log component links are available for all three estimation methods.

exposureEM

exposureEM fits two-component exposure models for binary outcomes when the setting in which an event occurred is not observed. The motivating method is described in A Combined Exposure Model for Binary Injury Histories: An EM Approach to Sport and Non-Sport Risk in Children, but the package itself is not specific to sports or health research.

Users can choose pure EM, direct Marquardt-damped Newton-Raphson, or a hybrid that runs several EM iterations before switching to direct optimization. The package does not implement expectation-solution estimation or application-specific data processing.

Model

For observation i, the observed binary event is generated by two latent component events:

Y_i = I(Y_0i + Y_1i >= 1)

The components have known exposures a0 and a1, separate intercepts, and shared covariate effects. In the primary censored-Poisson/complementary-log-log specification,

lambda_ki = a_ki * exp(x_ki' beta)
P(Y_ki = 1) = 1 - exp(-lambda_ki).

Logit and log component links are also available as sensitivity models.

Installation

From the directory containing the package source archive:

install.packages("exposureEM_0.3.0.tar.gz", repos = NULL, type = "source")

During development, the package can be loaded from its source directory:

pkgload::load_all("exposureEM")

Choose an estimation method

library(exposureEM)

fit_em <- exposure_fit(
  event ~ age + group + baseline_score,
  data = analysis_data,
  exposure = c("a0", "a1"),
  link = "cloglog",
  method = "em"
)

fit_direct <- exposure_fit(
  event ~ age + group + baseline_score,
  data = analysis_data,
  exposure = c("a0", "a1"),
  link = "cloglog",
  method = "marquardt"
)

fit_hybrid <- exposure_fit(
  event ~ age + group + baseline_score,
  data = analysis_data,
  exposure = c("a0", "a1"),
  link = "cloglog",
  method = "hybrid",
  hybrid_em_iterations = 10
)

summary(fit_hybrid)
component_contrast(fit_hybrid)
posterior_components(fit_hybrid)

The first exposure column belongs to component 0 and the second to component 1. By default, both must lie in [0, 1] and sum to one. Zero exposure to either component is valid.

exposure_em() remains available as a backward-compatible shortcut for a pure EM fit. All three methods support cloglog, logit, and log component links. Set link = "logit" or link = "log" in any call above; a hybrid uses the same link in both phases.

For a hybrid fit, hybrid_em_iterations is the maximum number of complete EM outer iterations before switching; EM convergence can trigger an earlier switch. A value such as 5 or 10 is a practical starting point, but the best choice depends on starting values and the dataset.

Main outputs

• coef(fit) returns the component-specific intercepts and shared covariate effects.
• vcov(fit) returns the Louis-formula covariance for EM or the inverse unmodified observed Hessian for direct and hybrid fits.
• component_contrast(fit) returns exp(beta_component_1 - beta_component_0) and its Wald confidence interval.
• posterior_components(fit) returns final E-step marginal and joint probabilities.
• predict(fit, newdata, type = "response") returns the probability of any observed event.
• fit$trace records the observed log-likelihood and numerical diagnostics. A hybrid trace labels its EM and Marquardt phases separately.

Numerical behavior

The EM M-step uses full Newton-Raphson iterations. Step-halving prevents invalid probabilities and requires each accepted inner step to improve the expected complete-data log-likelihood.

The direct method instead updates the observed-data likelihood. When the negative observed Hessian is not positive definite, Marquardt diagonal inflation supplies a stable ascent direction. Backtracking shortens the step until it sufficiently increases the observed likelihood. The inflation is removed before calculating the final covariance matrix. With the log link, backtracking rejects any proposal with a component probability at or above one. Boundary solutions and failed backtracking are reported as nonconvergence; the optimizer does not clip probabilities or silently change the link.

The default outer tolerance is 1e-8. EM checks its selected change criterion, observed likelihood change, and the score scaled by unmodified observed information. Direct fitting checks squared parameter change, likelihood change, and the same score measure. Inspect fit$score, fit$newton_decrement, and fit$converged; small steps alone do not establish identification or a global maximum. Invalid observed information produces NA covariance with a warning. Model-based standard errors assume independent observations. No elapsed-time advantage of the hybrid is claimed.

Run the standalone smoke test

After installing the package, run:

source(system.file("examples", "test_exposureEM.R", package = "exposureEM"))

The script simulates a reproducible dataset, demonstrates EM, direct Marquardt, and hybrid fitting, checks all three EM component links, exercises predictions and posterior probabilities, and stops if any numerical or API check fails.

Check against the saved ABCD results

On the authorized computer that contains the original ABCD files, first install the current package and then run this command from the package source folder:

Rscript tests/abcd-reference.R

This quick integration test rebuilds the analysis data with the same canonical runner used for the earlier fits. It checks the adjusted and unadjusted primary models through EM, direct Marquardt, and hybrid fitting against EM_BL_to_Y3_results_all_20260226.rds, including sample counts, component coefficients, the component rate ratio and confidence interval, and AIC. The console prints the saved and newly fitted estimates side by side, and then prints the largest numerical differences and the final pass/fail result. It does not create a separate report file. The test never copies participant-level ABCD data into the package or prints participant records.

To check all saved primary cloglog models or all usable saved EM models, use one of these longer commands:

EXPOSUREEM_ABCD_MODE=primary Rscript tests/abcd-reference.R
EXPOSUREEM_ABCD_MODE=full Rscript tests/abcd-reference.R

The two input paths can be changed without editing the script by setting EXPOSUREEM_ABCD_RUNNER and EXPOSUREEM_ABCD_RESULTS.