The marp package fits six parametric renewal-process
models to positive inter-event times: Poisson (exponential waiting
times), Gamma, log-logistic, Weibull, log-normal, and Brownian passage
time (BPT). It compares the candidate models using AIC and calculates
AIC-weighted model-averaged estimates.
The main fitted quantities are the mean inter-event time, the
cumulative event probability at a user-supplied time y
stored on the logit scale, and the hazard evaluated at
user-supplied times t stored on the log
scale.
This vignette demonstrates the standard fitted-object workflow. The final confidence-interval step is shown but not evaluated because the existing studentized procedure uses nested bootstrap sampling and is intentionally computationally expensive.
marp() expects a numeric vector of positive inter-event
times. The following small reproducible example is simulated from a
Gamma distribution.
Install the current version of marp from GitHub
with:
if (requireNamespace("marp", quietly = TRUE)) {
library(marp)
} else {
devtools::load_all("..")
}
set.seed(42)
dat <- rgamma(50, shape = 3, rate = 0.01)
summary(dat)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 48.94 164.41 270.19 292.30 367.20 670.99The remaining inputs control the displayed estimates. t
gives the times for log-hazard evaluation, y gives the time
for the logit event probability, and m controls repeated
random-start optimizations in candidate models that use
nlm().
The reference/generating-model codes are 1 Poisson, 2 Gamma, 3
log-logistic, 4 Weibull, 5 log-normal, and 6 BPT. A reference model is
most naturally known in a simulation; here code 2 matches the
distribution used to generate dat.
The primary interface fits all six candidates, selects the lowest-AIC model, and calculates AIC-weighted estimates.
set.seed(42)
fit <- marp(dat, t, m, y, which.model = 2)
fit
#> Model-averaged renewal-process fit
#> Observations: 50
#>
#> Model comparison:
#> Model AIC AIC weight Status
#> Poisson renewal (exponential waiting times) 669.7785 2.931522e-07 ok
#> Gamma 641.3777 4.307620e-01 ok
#> Log-logistic 645.1522 6.525531e-02 ok
#> Weibull 641.9885 3.173946e-01 ok
#> Log-normal 644.0323 1.142342e-01 ok
#> Brownian passage time 644.9457 7.235364e-02 ok
#>
#> Best model: Gamma
#> Best-model mean: 292.3017
#> Model-averaged mean: 294.4973
#> Model-averaged logit event probability: 0.3966815
#> Model-averaged log-hazard evaluations: 6 storedPrinting the fitted object gives a concise model-comparison table, AIC weights, the selected model, and the principal model-averaged estimates. A structured summary provides the model-specific parameters and fit criteria.
fit_summary <- summary(fit)
fit_summary
#> Model-averaged renewal-process fit summary
#> Observations: 50
#>
#> Model comparison:
#> model par1 par2
#> Poisson renewal (exponential waiting times) 3.421129e-03 NA
#> Gamma 3.449353e+00 0.01180066
#> Log-logistic 3.005431e+00 258.95021871
#> Weibull 3.309584e+02 2.02295201
#> Log-normal 5.525885e+00 0.58013126
#> Brownian passage time 2.923001e+02 0.63165008
#> logLik AIC BIC AIC_weight status
#> -333.8892 669.7785 671.6905 2.931522e-07 ok
#> -318.6889 641.3777 645.2018 4.307620e-01 ok
#> -320.5761 645.1522 648.9762 6.525531e-02 ok
#> -318.9943 641.9885 645.8126 3.173946e-01 ok
#> -320.0162 644.0323 647.8564 1.142342e-01 ok
#> -320.4728 644.9457 648.7697 7.235364e-02 ok
#>
#> Selected model: Gamma
#> Selected-model mean: 292.3017
#> Selected-model logit event probability: 0.4054542
#>
#> Model-averaged estimates:
#> Mean: 294.4973
#> Logit event probability at y = 304 : 0.3966815
#> Log-hazard evaluations: 6 stored
#>
#> Reference model: GammaThe six rows are always in the documented candidate-model order. Parameters are model-specific, so they should be interpreted using the documentation for the corresponding model-specific fitter. For example, the Gamma parameters are shape and rate, whereas the Weibull parameters are scale and shape.
The S3 print and summary methods are the usual entry points. Existing named list components remain available for programmatic use and backward compatibility.
fit$weights_AIC
#> [1] 2.931522e-07 4.307620e-01 6.525531e-02 3.173946e-01 1.142342e-01
#> [6] 7.235364e-02
fit$model_best
#> [1] 2
fit$mu_best
#> [1] 292.3017
fit$mu_aic
#> [1] 294.4973
fit$pr_aic
#> [1] 0.3966815
fit$haz_aic
#> [1] -6.279153 -6.031304 -5.842676 -5.693959 -5.573597 -5.474192mu_best is the mean estimate from the lowest-AIC model,
while mu_aic is the AIC-weighted mean estimate.
pr_aic is an AIC-weighted average on the logit
event-probability scale. haz_aic contains
AIC-weighted values on the log-hazard scale, with one
value for every element of t.
The supplied reference model is also retained. These quantities are useful for simulation comparisons and do not imply that the true model is known in an observational application.
The standard confint() method delegates to the package’s
existing marp_confint() implementation. Because the fitted
object deliberately does not store the original observations, the
original data and bootstrap sizes must be supplied explicitly.
The public function can equivalently be called directly with
significance level alpha = 1 - level:
ci_direct <- marp_confint(
data = dat,
m = m,
t = t,
B = 99,
BB = 99,
alpha = 0.05,
y = y,
which.model = 2
)Both forms use the same numerical engine. The studentized interval
calculation generates BB double-bootstrap samples for each
of B bootstrap samples and fits multiple renewal models
repeatedly. Production runs can therefore take substantially longer than
the point-estimation workflow. Set a random seed for reproducibility and
choose B and BB appropriate to the analysis;
the values shown above illustrate the API rather than prescribing
settings for every use.
A typical analysis uses marp() followed by
print() and summary() to compare candidate
models and inspect selected and model-averaged estimates. Named
components support programmatic extraction, while confint()
provides a conventional interface to the existing bootstrap
confidence-interval engine.