| Title: | Model-Averaged Renewal Process |
|---|---|
| Description: | To implement a model-averaging approach with different renewal models, with a primary focus on forecasting large earthquakes. Based on six renewal models (i.e., Poisson, Gamma, Log-Logistics, Weibull, Log-Normal and BPT), model-averaged point estimates are calculated using AIC weights. Additionally, both percentile and studentized bootstrapped model-averaged confidence intervals are constructed. In comparison, point and interval estimation from the individual or "best" model (determined via model selection) can be retrieved. |
| Authors: | Jie Kang [aut], Chris Scott [aut], Vanessa Huang [aut], Veronica Tsou [aut, cre], Albert Savary [ctb] |
| Maintainer: | Veronica Tsou <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.1 |
| Built: | 2026-09-21 09:44:41 UTC |
| Source: | https://github.com/kanji709/marp |
A function to generate (double) bootstrap samples and fit BPT renewal model
bpt_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)bpt_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)
n |
Number of inter-event times generated in each bootstrap sample. |
t |
Time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
m |
Positive integer controlling repeated random-start optimizations; see [marp()]. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6; see [marp()] for model order. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards with 'length(t)' rows and six model columns. |
y |
Time point at which logit event probabilities are evaluated. |
A list of bootstrap estimates and variance/T-statistic quantities. Components beginning with 'pr_' use the logit-probability scale and those beginning with 'haz_' use the log-hazard scale. If the supplied BPT fit is unavailable, the same structure is returned with 'NA' values.
Estimated mean from bootstrapped samples
Logit event probabilities from bootstrap samples
Log-hazards from bootstrap samples
Variance of estimated mean
Variance of estimated probability
Variance of estimated hazard rates
Variance of estimated mean of bootstrapped samples (via double-bootstrapping)
Variance of estimated probability of bootstrapped samples (via double-bootstrapping)
Variance of estimated hazard rates of bootstrapped samples (via double-bootstrapping)
Pivot quantity of the estimated mean
Pivot quantity of the estimated probability
Pivot quantity of the estimated hazard rates
## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off point for probablity estimation # generate bootstrapped samples then fit renewal model res <- marp::bpt_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off point for probablity estimation # generate bootstrapped samples then fit renewal model res <- marp::bpt_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)
A function to calculate the log-likelihood of BPT model
bpt_logl(param, x)bpt_logl(param, x)
param |
Length-2 numeric vector containing log(mean) and log(alpha squared), where alpha is the BPT aperiodicity parameter. |
x |
Numeric vector of positive observations. |
returns the value of negative log-likelihood of the BPT model
set.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(292.945125794581, 0.718247184450307) # estimated parameters param <- c(log(par_hat[1]),log(par_hat[2]^2)) # input parameters for logl function # calculate log-likelihood result <- marp::bpt_logl(param, data) resultset.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(292.945125794581, 0.718247184450307) # estimated parameters param <- c(log(par_hat[1]),log(par_hat[2]^2)) # input parameters for logl function # calculate log-likelihood result <- marp::bpt_logl(param, data) result
A function to fit BPT renewal model
bpt_rp(data, t, m, y)bpt_rp(data, t, m, y)
data |
A numeric vector of positive inter-event times. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
m |
A positive integer controlling repeated random-start optimizations; the current implementation seeks 'm - 1' acceptable 'nlm()' fits and caps the number of attempted starts at '20 * m'. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
An object of class 'marp_model_fit' retaining the following eight named list components:
Estimated mean parameter (mu) of the BPT model
Estimated aperiodicity parameter (alpha) of the BPT model
Maximized log-likelihood
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event time
Logit-transformed cumulative event probability at 'y'
Log-hazard values at 't'
If the requested plausible fits cannot be obtained within the existing attempt limit, the function warns and returns the same components with unavailable estimates set to 'NA' and information criteria set to 'Inf'.
set.seed(42) data <- rgamma(30,3,0.01) # set some parameters m <- 10 # repeated random-start optimization setting t <- seq(100, 200, by=10) # time intervals y <- 304 # cut-off year for estimating probablity # fit BPT renewal model fit <- marp::bpt_rp(data, t, m, y) fit summary(fit)set.seed(42) data <- rgamma(30,3,0.01) # set some parameters m <- 10 # repeated random-start optimization setting t <- seq(100, 200, by=10) # time intervals y <- 304 # cut-off year for estimating probablity # fit BPT renewal model fit <- marp::bpt_rp(data, t, m, y) fit summary(fit)
Extract parameters from a model-specific fit
## S3 method for class 'marp_model_fit' coef(object, ...)## S3 method for class 'marp_model_fit' coef(object, ...)
object |
A 'marp_model_fit' object. |
... |
Additional arguments, currently unused. |
A named numeric vector of fitted distribution parameters.
This delegates to [marp_confint()] without changing its bootstrap calculations. The original data are deliberately not stored in 'object', so they must be supplied explicitly.
## S3 method for class 'marp_fit' confint(object, parm = NULL, level = 0.95, data, B, BB, ...)## S3 method for class 'marp_fit' confint(object, parm = NULL, level = 0.95, data, B, BB, ...)
object |
A 'marp_fit' object. |
parm |
Currently unsupported. It must be 'NULL'. |
level |
Confidence level. The default is 0.95. |
data |
The original inter-event times used to fit 'object'. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples. |
... |
Additional arguments passed to [marp_confint()]. |
The existing nested confidence-interval result returned by [marp_confint()], with class 'marp_confint'.
Density function of Log-Logistics model
dllog(x, shape = 1, scale = 1, log = FALSE)dllog(x, shape = 1, scale = 1, log = FALSE)
x |
Numeric vector of positive quantiles. |
shape |
Positive shape parameter of the log-logistic distribution. |
scale |
Positive scale parameter of the log-logistic distribution. |
log |
Logical; if 'TRUE', return log-densities. |
A numeric vector of densities, or log-densities when 'log = TRUE'.
x <- as.numeric(c(350., 450., 227., 352., 654.)) # set paramters shape <- 5 scale <- 3 log <- FALSE result_1 <- marp::dllog(x, shape, scale, log) # alternatively, set log == TRUE log <- TRUE result_2 <- marp::dllog(x, shape, scale, log)x <- as.numeric(c(350., 450., 227., 352., 654.)) # set paramters shape <- 5 scale <- 3 log <- FALSE result_1 <- marp::dllog(x, shape, scale, log) # alternatively, set log == TRUE log <- TRUE result_2 <- marp::dllog(x, shape, scale, log)
A function to generate (double) bootstrap samples and fit Gamma renewal model
gamma_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)gamma_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)
n |
Number of inter-event times generated in each bootstrap sample. |
t |
Time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
m |
Positive integer controlling repeated random-start optimizations; see [marp()]. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6; Gamma uses shape and rate. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards with 'length(t)' rows and six model columns. |
y |
Time point at which logit event probabilities are evaluated. |
A list of bootstrap estimates and variance/T-statistic quantities. Components beginning with 'pr_' use the logit-probability scale and those beginning with 'haz_' use the log-hazard scale.
Estimated mean from bootstrapped samples
Logit event probabilities from bootstrap samples
Log-hazards from bootstrap samples
Variance of estimated mean
Variance of logit event probabilities
Variance of log-hazards
Variance of estimated mean of bootstrapped samples (via double-bootstrapping)
Variance of estimated probability of bootstrapped samples (via double-bootstrapping)
Variance of estimated hazard rates of bootstrapped samples (via double-bootstrapping)
Pivot quantity of the estimated mean
Pivot quantity of the estimated probability
Pivot quantity of the estimated hazard rates
## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::gamma_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::gamma_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)
A function to calculate the log-likelihood of Gamma model
gamma_logl(param, x)gamma_logl(param, x)
param |
Length-2 numeric vector containing log(shape) and log(rate). |
x |
Numeric vector of positive observations. |
returns the value of negative log-likelihood of the Gamma model
set.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(2.7626793657057762, 0.0094307059277139432) # estimated parameters param <- log(par_hat) # input parameters for logl function # calculate log-likelihood result <- marp::gamma_logl(param, data) resultset.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(2.7626793657057762, 0.0094307059277139432) # estimated parameters param <- log(par_hat) # input parameters for logl function # calculate log-likelihood result <- marp::gamma_logl(param, data) result
A function to fit Gamma renewal model
gamma_rp(data, t, m, y)gamma_rp(data, t, m, y)
data |
A numeric vector of positive inter-event times. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
m |
A positive integer controlling repeated random-start optimizations; the current implementation seeks 'm - 1' acceptable 'nlm()' fits. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
An object of class 'marp_model_fit' retaining the following eight named list components:
Estimated shape parameter of the Gamma model
Estimated rate parameter of the Gamma model
Maximized log-likelihood
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event time
Logit-transformed cumulative event probability at 'y'
Log-hazard values at 't'
set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Gamma renewal model fit <- marp::gamma_rp(data, t, m, y) fit summary(fit)set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Gamma renewal model fit <- marp::gamma_rp(data, t, m, y) fit summary(fit)
Extract the log-likelihood from a model-specific fit
## S3 method for class 'marp_model_fit' logLik(object, ...)## S3 method for class 'marp_model_fit' logLik(object, ...)
object |
A 'marp_model_fit' object. |
... |
Additional arguments, currently unused. |
An object of class 'logLik'.
A function to generate (double) bootstrap samples and fit Log-Logistic renewal model
loglogis_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)loglogis_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)
n |
Number of inter-event times generated in each bootstrap sample. |
t |
Time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
m |
Positive integer controlling repeated random-start optimizations; see [marp()]. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6; see [marp()] for model order. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards with 'length(t)' rows and six model columns. |
y |
Time point at which logit event probabilities are evaluated. |
A list of bootstrap estimates and variance/T-statistic quantities. Components beginning with 'pr_' use the logit-probability scale and those beginning with 'haz_' use the log-hazard scale.
Estimated mean from bootstrapped samples
Logit event probabilities from bootstrap samples
Log-hazards from bootstrap samples
Variance of estimated mean
Variance of estimated probability
Variance of estimated hazard rates
Variance of estimated mean of bootstrapped samples (via double-bootstrapping)
Variance of estimated probability of bootstrapped samples (via double-bootstrapping)
Variance of estimated hazard rates of bootstrapped samples (via double-bootstrapping)
Pivot quantity of the estimated mean
Pivot quantity of the estimated probability
Pivot quantity of the estimated hazard rates
## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::loglogis_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::loglogis_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)
A function to calculate the log-likelihood of Log-Logistics model
loglogis_logl(param, x)loglogis_logl(param, x)
param |
Length-2 numeric vector containing log(scale) and log(shape). |
x |
Numeric vector of positive observations. |
returns the value of negative log-likelihood of the Log-Logistics model
set.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(2.6037079185931518, 247.59811806509711) # estimated parameters param <- c(log(par_hat[2]),log(par_hat[1])) # input parameters for logl function # calculate log-likelihood result <- marp::loglogis_logl(param, data) resultset.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(2.6037079185931518, 247.59811806509711) # estimated parameters param <- c(log(par_hat[2]),log(par_hat[1])) # input parameters for logl function # calculate log-likelihood result <- marp::loglogis_logl(param, data) result
A function to fit Log-Logistics renewal model
loglogis_rp(data, t, m, y)loglogis_rp(data, t, m, y)
data |
A numeric vector of positive inter-event times. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
m |
A positive integer controlling repeated random-start optimizations; the current implementation seeks 'm - 1' acceptable 'nlm()' fits. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
An object of class 'marp_model_fit' retaining the following eight named list components:
Estimated shape parameter of the Log-Logistics model
Estimated scale parameter of the Log-Logistics model
Maximized log-likelihood
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event time
Logit-transformed cumulative event probability at 'y'
Log-hazard values at 't'
set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Log-Logistic renewal model fit <- marp::loglogis_rp(data, t, m, y) fit summary(fit)set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Log-Logistic renewal model fit <- marp::loglogis_rp(data, t, m, y) fit summary(fit)
A function to generate (double) bootstrap samples and fit Log-Normal renewal model
lognorm_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y)lognorm_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y)
n |
Number of inter-event times generated in each bootstrap sample. |
t |
Time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6; see [marp()] for model order. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards, with 'length(t)' rows and six model columns. |
y |
Time point at which logit event probabilities are evaluated. |
A list of bootstrap estimates and variance/T-statistic quantities. Components beginning with 'pr_' use the logit-probability scale and those beginning with 'haz_' use the log-hazard scale.
Estimated mean from bootstrapped samples
Logit event probabilities from bootstrap samples
Log-hazards from bootstrap samples
Variance of estimated mean
Variance of logit event probabilities
Variance of log-hazards
Variance of estimated mean of bootstrapped samples (via double-bootstrapping)
Variance of estimated probability of bootstrapped samples (via double-bootstrapping)
Variance of estimated hazard rates of bootstrapped samples (via double-bootstrapping)
Pivot quantity of the estimated mean
Pivot quantity of the estimated probability
Pivot quantity of the estimated hazard rates
## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps # m <- 10 # number of iterations for MLE optimization par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::lognorm_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps # m <- 10 # number of iterations for MLE optimization par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::lognorm_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)
A function to fit Log-Normal renewal model
lognorm_rp(data, t, y)lognorm_rp(data, t, y)
data |
A numeric vector of positive inter-event times. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
An object of class 'marp_model_fit' retaining the following eight named list components:
Estimated mean (on the log scale) of the Log-Normal model
Estimated standard deviation (on the log scale)of the Log-Normal model
Maximized log-likelihood
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event time
Logit-transformed cumulative event probability at 'y'
Log-hazard values at 't'
set.seed(42) data <- rgamma(100,3,0.01) # set some parameters t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Log-Normal renewal model fit <- marp::lognorm_rp(data, t, y) fit summary(fit)set.seed(42) data <- rgamma(100,3,0.01) # set some parameters t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Log-Normal renewal model fit <- marp::lognorm_rp(data, t, y) fit summary(fit)
Evaluate the lower-limit equation for a model-averaged T statistic
lowerT(low, hat, sigmasq, Tstar, weights, B, alpha)lowerT(low, hat, sigmasq, Tstar, weights, B, alpha)
low |
Candidate lower confidence limit. |
hat |
Vector of model-specific estimates. |
sigmasq |
Vector of model-specific variance estimates. |
Tstar |
Matrix of bootstrap T statistics, with models in rows. |
weights |
Vector of model weights. |
B |
Number of bootstrap samples represented in 'Tstar'. |
alpha |
Significance level. |
A numeric scalar giving the lower-limit root equation value.
# set some parameters low <- 100 # lower bound hat <- rep(150, 6) # estimates obtained from each model sigmasq <- 10 # variance Tstar <- matrix(rep(100,600),6,100) # T statistics estimated from bootstrap samples weights <- rep(1/6, 6) # model weights B <- 100 # number of bootstrapped samples alpha <- 0.05 # confidence level # calculate the upper limit of T statistics res <- marp::lowerT(low, hat, sigmasq, Tstar, weights, B, alpha) res# set some parameters low <- 100 # lower bound hat <- rep(150, 6) # estimates obtained from each model sigmasq <- 10 # variance Tstar <- matrix(rep(100,600),6,100) # T statistics estimated from bootstrap samples weights <- rep(1/6, 6) # model weights B <- 100 # number of bootstrapped samples alpha <- 0.05 # confidence level # calculate the upper limit of T statistics res <- marp::lowerT(low, hat, sigmasq, Tstar, weights, B, alpha) res
A function to apply model-averaged renewal process
marp(data, t, m, y, which.model = 1)marp(data, t, m, y, which.model = 1)
data |
A numeric vector of positive inter-event times. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
m |
A positive integer controlling repeated random-start optimizations for the Gamma, log-logistic, Weibull, and BPT candidate models. The current implementation seeks 'm - 1' acceptable 'nlm()' fits for each such model. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
which.model |
Integer identifying a reference or generating model: 1 = Poisson, 2 = Gamma, 3 = log-logistic, 4 = Weibull, 5 = log-normal, and 6 = Brownian passage time (BPT). This is mainly useful when the generating model is known, such as in simulations. |
An object of class 'marp_fit': a 19-component named list of model-specific, selected-model, reference-model, and AIC-weighted estimates.
First fitted parameter for each model, in the model order listed under 'which.model'. Its meaning is model-specific.
Second fitted parameter for each model in the same order; 'NA' for the one-parameter Poisson model. Its meaning is model-specific.
Maximized log-likelihood for each model.
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event time for each model.
Logit-transformed cumulative event probability at 'y' for each model.
Log-hazard values at 't' for each model.
Model weights calculated based on AIC
Model selected based on the lowest AIC
Estimated mean obtained from the model with the lowest AIC
Estimated logit event probability from the model with the lowest AIC.
Estimated log-hazards from the model with the lowest AIC.
Estimated mean obtained from the (true or hypothetical) generating model
Estimated logit event probability from the reference model.
Estimated log-hazards from the reference model.
Estimated mean obtained from model-averaging (using AIC weights)
AIC-weighted average of the model-specific logit event probabilities.
AIC-weighted averages of the model-specific log-hazards at 't'.
set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probability which.model <- 2 # specify the generating model # model selection and averaging result <- marp::marp(data, t, m, y, which.model) result summary(result) # Bootstrap confidence intervals delegate to the existing marp_confint() # engine. The original data and bootstrap sizes are supplied explicitly. ## Not run: confint(result, data = data, B = 99, BB = 99) ## End(Not run)set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probability which.model <- 2 # specify the generating model # model selection and averaging result <- marp::marp(data, t, m, y, which.model) result summary(result) # Bootstrap confidence intervals delegate to the existing marp_confint() # engine. The original data and bootstrap sizes are supplied explicitly. ## Not run: confint(result, data = data, B = 99, BB = 99) ## End(Not run)
A function to fit model-averaged renewal process
marp_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)marp_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)
n |
Number of inter-event times generated in each bootstrap sample. |
t |
Time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
m |
Positive integer controlling repeated random-start optimizations; see [marp()]. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6; see [marp()] for model order. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards with 'length(t)' rows and six model columns. |
y |
Time point at which logit event probabilities are evaluated. |
A list combining bootstrap results from all six candidate models. Components beginning with 'pr_' use the logit-probability scale and those beginning with 'haz_' use the log-hazard scale.
Estimated mean from bootstrapped samples
A 6 by 'B' matrix of logit event probabilities
A 'length(t)' by 6 by 'B' array of log-hazards
Variance of estimated mean
Variance of model-specific logit event probabilities
Variance of model-specific log-hazards
Variance of estimated mean of bootstrapped samples (via double-bootstrapping)
Variance of estimated probability of bootstrapped samples (via double-bootstrapping)
Variance of estimated hazard rates of bootstrapped samples (via double-bootstrapping)
Pivot quantity of the estimated mean
Pivot quantity of the estimated probability
Pivot quantity of the estimated hazard rates
## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::marp_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::marp_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)
A function to apply model-averaged renewal process
marp_confint(data, m, t, B, BB, alpha, y, which.model)marp_confint(data, m, t, B, BB, alpha, y, which.model)
data |
A numeric vector of positive inter-event times. |
m |
A positive integer controlling repeated random-start optimizations; see [marp()] for the exact role of this argument. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples. |
alpha |
Significance level; the confidence level is '1 - alpha'. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
which.model |
Integer identifying a reference or generating model: 1 = Poisson, 2 = Gamma, 3 = log-logistic, 4 = Weibull, 5 = log-normal, and 6 = BPT. |
The returned object has three top-level components: 'out', containing the fitted/model-averaged estimates and bootstrap-weighted estimates; 'percent_CI', containing percentile intervals; and 'student_CI', containing studentized intervals. The object has class 'marp_confint'; its three existing list components remain directly accessible with '$'. Probability quantities are stored on the logit scale and hazard quantities on the log-hazard scale throughout.
The 19 fitted/model-averaged components returned by [marp()], followed by 'mu_bstrp', 'pr_bstrp', and 'haz_bstrp'.
Bootstrap model weights and percentile summaries for reference-model and selected-model means, logit probabilities, and log-hazards.
Studentized lower and upper limits for reference-model, selected-model, and model-averaged quantities.
First fitted parameter for each model; its meaning is model-specific.
Second fitted parameter for each model; its meaning is model-specific.
Maximized log-likelihood
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event times
Logit-transformed cumulative event probabilities at 'y'
Log-hazard values at 't'
Model weights calculated based on AIC
Model selected based on the lowest AIC
Estimated mean obtained from the model with the lowest AIC
Estimated logit event probability from the model with the lowest AIC
Estimated log-hazards from the model with the lowest AIC
Estimated mean obtained from the (true or hypothetical) generating model
Estimated logit event probability from the reference model
Estimated log-hazards from the reference model
Estimated mean obtained from model-averaging (using AIC weights)
AIC-weighted average of logit event probabilities
AIC-weighted averages of log-hazards
Estimated mean obtained from model-averaging (using bootstrapped weights)
Bootstrap-weighted average of logit event probabilities
Bootstrap-weighted averages of log-hazards
Model weights calculated by bootstrapping, that is, the frequency of each model being selected as the best model is divided by the total number of bootstraps
Median of the percentile bootstrap confidence interval of the estimated mean based on the generating model
Lower limit of the percentile bootstrap confidence interval of the estimated mean based on the generating model
Upper limit of the percentile bootstrap confidence interval of the estimated mean based on the generating model
Median of the percentile bootstrap confidence interval of the estimated mean based on the best model
Lower limit of the percentile bootstrap confidence interval of the estimated mean based on the best model
Upper limit of the percentile bootstrap confidence interval of the estimated mean based on the best model
Median reference-model logit event probability
Lower limit for the reference-model logit event probability
Upper limit for the reference-model logit event probability
Median best-model logit event probability
Lower limit for the best-model logit event probability
Upper limit for the best-model logit event probability
Median reference-model log-hazards
Lower limits for reference-model log-hazards
Upper limits for reference-model log-hazards
Median best-model log-hazards
Lower limits for best-model log-hazards
Upper limits for best-model log-hazards
Lower limit of the studentized bootstrap confidence interval of the estimated mean based on the generating model
Upper limit of the studentized bootstrap confidence interval of the estimated mean based on the generating model
Lower limit of the studentized bootstrap confidence interval of the estimated mean based on the best model
Upper limit of the studentized bootstrap confidence interval of the estimated mean based on the best model
Lower studentized limit for the reference-model logit probability
Upper studentized limit for the reference-model logit probability
Lower studentized limit for the best-model logit probability
Upper studentized limit for the best-model logit probability
Lower studentized limits for reference-model log-hazards
Upper studentized limits for reference-model log-hazards
Lower studentized limits for best-model log-hazards
Upper studentized limits for best-model log-hazards
Lower limit of model-averaged studentized bootstrap confidence interval of the estimated mean
Upper limit of model-averaged studentized bootstrap confidence interval of the estimated mean
Lower model-averaged limit on the logit-probability scale
Upper model-averaged limit on the logit-probability scale
Lower model-averaged limits on the log-hazard scale
Upper model-averaged limits on the log-hazard scale
## Not run: # generate random data set.seed(42) data <- rgamma(30, 3, 0.01) # set some parameters m <- 10 # repeated random-start optimization setting t <- seq(100,200,by=10) # time intervals alpha <- 0.05 # confidence level y <- 304 # cut-off year for estimating probability B <- 100 # number of bootstraps BB <- 100 # number of double bootstraps which.model <- 2 # specify the generating model # Construct confidence intervals. Production choices of B and BB may be # computationally expensive because BB double bootstraps are generated for # each of B bootstrap samples. res <- marp::marp_confint(data,m,t,B,BB,alpha,y,which.model) res ## End(Not run)## Not run: # generate random data set.seed(42) data <- rgamma(30, 3, 0.01) # set some parameters m <- 10 # repeated random-start optimization setting t <- seq(100,200,by=10) # time intervals alpha <- 0.05 # confidence level y <- 304 # cut-off year for estimating probability B <- 100 # number of bootstraps BB <- 100 # number of double bootstraps which.model <- 2 # specify the generating model # Construct confidence intervals. Production choices of B and BB may be # computationally expensive because BB double bootstraps are generated for # each of B bootstrap samples. res <- marp::marp_confint(data,m,t,B,BB,alpha,y,which.model) res ## End(Not run)
Number of observations in a model-averaged fit
## S3 method for class 'marp_fit' nobs(object, ...)## S3 method for class 'marp_fit' nobs(object, ...)
object |
A 'marp_fit' object. |
... |
Additional arguments, currently unused. |
The number of observations used for fitting.
Number of observations in a model-specific fit
## S3 method for class 'marp_model_fit' nobs(object, ...)## S3 method for class 'marp_model_fit' nobs(object, ...)
object |
A 'marp_model_fit' object. |
... |
Additional arguments, currently unused. |
The number of observations used for fitting.
A function to calculate percentile bootstrap confidence interval
percent_confint(data, B, t, m, y, which.model = 1)percent_confint(data, B, t, m, y, which.model = 1)
data |
A numeric vector of positive inter-event times. |
B |
Number of nonparametric bootstrap samples. |
t |
Time points at which log-hazards are evaluated. |
m |
A positive integer controlling repeated random-start optimizations; see [marp()]. |
y |
Time point at which logit event probabilities are evaluated. |
which.model |
Integer identifying the reference model using the mapping 1 = Poisson, 2 = Gamma, 3 = log-logistic, 4 = Weibull, 5 = log-normal, and 6 = BPT. |
A list of 95 percent percentile-bootstrap summaries. Components beginning with 'pr_' are on the logit-probability scale and components beginning with 'haz_' are on the log-hazard scale.
Model weights calculated by bootstrapping, that is, the frequency of each model being selected as the best model is divided by the total number of bootstraps
Median of the percentile bootstrap confidence interval of the estimated mean based on the generating model
Lower limit of the percentile bootstrap confidence interval of the estimated mean based on the generating model
Upper limit of the percentile bootstrap confidence interval of the estimated mean based on the generating model
Median of the percentile bootstrap confidence interval of the estimated mean based on the best model
Lower limit of the percentile bootstrap confidence interval of the estimated mean based on the best model
Upper limit of the percentile bootstrap confidence interval of the estimated mean based on the best model
Median logit event probability based on the reference model
Lower limit for the reference-model logit event probability
Upper limit for the reference-model logit event probability
Median logit event probability based on the best model
Lower limit for the best-model logit event probability
Upper limit for the best-model logit event probability
Median log-hazards based on the reference model
Lower limits for reference-model log-hazards
Upper limits for reference-model log-hazards
Median log-hazards based on the best model
Lower limits for best-model log-hazards
Upper limits for best-model log-hazards
## Not run: # generate random data set.seed(42) data <- rgamma(30, 3, 0.01) # set some parameters m <- 10 # repeated random-start optimization setting t <- seq(100,200,by=10) # time intervals y <- 304 # time point for estimating event probability B <- 100 # number of bootstraps BB <- 100 # number of double bootstraps which.model <- 2 # specify the generating model # construct percentile bootstrap confidence intervals marp::percent_confint(data, B, t, m, y, which.model) ## End(Not run)## Not run: # generate random data set.seed(42) data <- rgamma(30, 3, 0.01) # set some parameters m <- 10 # repeated random-start optimization setting t <- seq(100,200,by=10) # time intervals y <- 304 # time point for estimating event probability B <- 100 # number of bootstraps BB <- 100 # number of double bootstraps which.model <- 2 # specify the generating model # construct percentile bootstrap confidence intervals marp::percent_confint(data, B, t, m, y, which.model) ## End(Not run)
Probability function of Log-Logistics model
pllog(q, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE)pllog(q, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE)
q |
Numeric vector of positive quantiles. |
shape |
Positive shape parameter of the log-logistic distribution. |
scale |
Positive scale parameter of the log-logistic distribution. |
lower.tail |
Logical; if 'TRUE', return lower-tail probabilities, otherwise return upper-tail probabilities. |
log.p |
Logical; if 'TRUE', return probabilities on the log scale. |
A numeric vector of probabilities, or log-probabilities when 'log.p = TRUE'.
q <- c(1, 2, 3, 4) # set paramters shape <- 5 scale <- 3 result_1 <- marp::pllog(q, shape, scale) # alternatively, return log-probabilities result_2 <- marp::pllog(q, shape, scale, log.p = TRUE)q <- c(1, 2, 3, 4) # set paramters shape <- 5 scale <- 3 result_1 <- marp::pllog(q, shape, scale) # alternatively, return log-probabilities result_2 <- marp::pllog(q, shape, scale, log.p = TRUE)
A function to generate (double) bootstrap samples and fit Poisson renewal model
poisson_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y)poisson_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y)
n |
Number of inter-event times generated in each bootstrap sample. |
t |
Time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6; see [marp()] for model order. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards, with 'length(t)' rows and six columns in the model order documented by [marp()]. |
y |
Time point at which logit event probabilities are evaluated. |
A list of bootstrap estimates and variance/T-statistic quantities. Components beginning with 'pr_' use the logit-probability scale and those beginning with 'haz_' use the log-hazard scale.
Estimated mean from bootstrapped samples
Logit event probabilities from bootstrap samples
Log-hazards from bootstrap samples
Variance of estimated mean
Variance of logit event probabilities
Variance of log-hazards
Variance of estimated mean of bootstrapped samples (via double-bootstrapping)
Variance of estimated probability of bootstrapped samples (via double-bootstrapping)
Variance of estimated hazard rates of bootstrapped samples (via double-bootstrapping)
Pivot quantity of the estimated mean
Pivot quantity of the estimated probability
Pivot quantity of the estimated hazard rates
## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps # m <- 10 # number of iterations for MLE optimization par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::poisson_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps # m <- 10 # number of iterations for MLE optimization par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01 ) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::poisson_bstrp(n, t, B, BB, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)
A function to fit Poisson renewal model
poisson_rp(data, t, y)poisson_rp(data, t, y)
data |
A numeric vector of positive inter-event times. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
An object of class 'marp_model_fit' retaining the following eight named list components:
Estimated exponential rate parameter of the Poisson renewal model
N/A, only keep it as a place holder for output formatting purpose
Maximized log-likelihood
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event time
Logit-transformed cumulative event probability at 'y'
Log-hazard values at 't'
set.seed(42) data <- rgamma(100,3,0.01) # set some parameters t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Poisson renewal model fit <- marp::poisson_rp(data, t, y) fit summary(fit)set.seed(42) data <- rgamma(100,3,0.01) # set some parameters t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Poisson renewal model fit <- marp::poisson_rp(data, t, y) fit summary(fit)
Print a model-averaged confidence-interval result
## S3 method for class 'marp_confint' print(x, ...)## S3 method for class 'marp_confint' print(x, ...)
x |
An object returned by [marp_confint()] or 'confint.marp_fit()'. |
... |
Additional arguments, currently unused. |
The input object, invisibly.
Print a model-averaged renewal-process fit
## S3 method for class 'marp_fit' print(x, ...)## S3 method for class 'marp_fit' print(x, ...)
x |
An object returned by [marp()]. |
... |
Additional arguments, currently unused. |
The input object, invisibly.
Print a model-specific renewal-process fit
## S3 method for class 'marp_model_fit' print(x, ...)## S3 method for class 'marp_model_fit' print(x, ...)
x |
An object returned by one of the model-specific fitting functions. |
... |
Additional arguments, currently unused. |
The input object, invisibly.
Print a summary of a model-averaged renewal-process fit
## S3 method for class 'summary_marp_fit' print(x, ...)## S3 method for class 'summary_marp_fit' print(x, ...)
x |
A 'summary_marp_fit' object. |
... |
Additional arguments, currently unused. |
The input object, invisibly.
Print a summary of a model-specific renewal-process fit
## S3 method for class 'summary_marp_model_fit' print(x, ...)## S3 method for class 'summary_marp_model_fit' print(x, ...)
x |
A 'summary_marp_model_fit' object. |
... |
Additional arguments, currently unused. |
The input object, invisibly.
A function to calculate Studentized bootstrap confidence interval
student_confint( n, B, t, m, BB, par_hat, mu_hat, pr_hat, haz_hat, weights, alpha, y, best.model, which.model = 1 )student_confint( n, B, t, m, BB, par_hat, mu_hat, pr_hat, haz_hat, weights, alpha, y, best.model, which.model = 1 )
n |
Number of inter-event times generated in each bootstrap sample. |
B |
Number of bootstrap samples. |
t |
Time points at which log-hazards are evaluated. |
m |
Positive integer controlling repeated random-start optimizations; see [marp()]. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards with 'length(t)' rows and six model columns. |
weights |
Length-6 vector of model weights. |
alpha |
Significance level; the confidence level is '1 - alpha'. |
y |
Time point at which logit event probabilities are evaluated. |
best.model |
Integer identifying the lowest-AIC model. |
which.model |
Integer identifying the reference model using the mapping 1 = Poisson, 2 = Gamma, 3 = log-logistic, 4 = Weibull, 5 = log-normal, and 6 = BPT. |
A list of studentized-bootstrap confidence limits. Components beginning with 'pr_' are on the logit-probability scale and components beginning with 'haz_' are on the log-hazard scale.
Lower limit of the studentized bootstrap confidence interval of the estimated mean based on the generating model
Upper limit of the studentized bootstrap confidence interval of the estimated mean based on the generating model
Lower limit of the studentized bootstrap confidence interval of the estimated mean based on the best model
Upper limit of the studentized bootstrap confidence interval of the estimated mean based on the best model
Lower limit for the reference-model logit event probability
Upper limit for the reference-model logit event probability
Lower limit for the best-model logit event probability
Upper limit for the best-model logit event probability
Lower limits for reference-model log-hazards
Upper limits for reference-model log-hazards
Lower limits for best-model log-hazards
Upper limits for best-model log-hazards
Lower limit of model-averaged studentized bootstrap confidence interval of the estimated mean
Upper limit of model-averaged studentized bootstrap confidence interval of the estimated mean
Lower model-averaged limit on the logit-probability scale
Upper model-averaged limit on the logit-probability scale
Lower model-averaged limits on the log-hazard scale
Upper model-averaged limits on the log-hazard scale
## Not run: # generate random data set.seed(42) data <- rgamma(30, 3, 0.01) # set some parameters n <- 30 # sample size m <- 10 # repeated random-start optimization setting t <- seq(100,200,by=10) # time intervals y <- 304 # cut-off year for estimating probablity B <- 100 # number of bootstraps BB <- 100 # number of double bootstraps par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) weights <- c(0.00000, 0.21000, 0.02000, 0.55000, 0.00000, 0.22000) # model weights alpha <- 0.05 # confidence level y <- 304 # cut-off year for estimating probablity best.model <- 2 which.model <- 2 # specify the generating model#' # construct Studentized bootstrap confidence interval marp::student_confint( n,B,t,m,BB,par_hat,mu_hat,pr_hat,haz_hat,weights,alpha,y,best.model,which.model ) ## End(Not run)## Not run: # generate random data set.seed(42) data <- rgamma(30, 3, 0.01) # set some parameters n <- 30 # sample size m <- 10 # repeated random-start optimization setting t <- seq(100,200,by=10) # time intervals y <- 304 # cut-off year for estimating probablity B <- 100 # number of bootstraps BB <- 100 # number of double bootstraps par_hat <- c( 3.41361e-03, 2.76268e+00, 2.60370e+00, 3.30802e+02, 5.48822e+00, 2.92945e+02, NA, 9.43071e-03, 2.47598e+02, 1.80102e+00, 6.50845e-01, 7.18247e-01) mu_hat <- c(292.94512, 292.94513, 319.72017, 294.16945, 298.87286, 292.94512) pr_hat <- c(0.60039, 0.42155, 0.53434, 0.30780, 0.56416, 0.61795) haz_hat <- matrix(c( -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -5.67999, -6.09420, -5.99679, -5.91174, -5.83682, -5.77031, -5.71085, -5.65738, -5.60904, -5.56512, -5.52504, -5.48833, -6.09902, -5.97017, -5.85769, -5.75939, -5.67350, -5.59856, -5.53336, -5.47683, -5.42805, -5.38621, -5.35060, -6.17146, -6.09512, -6.02542, -5.96131, -5.90194, -5.84668, -5.79498, -5.74642, -5.70064, -5.65733, -5.61624, -5.92355, -5.80239, -5.70475, -5.62524, -5.55994, -5.50595, -5.46106, -5.42359, -5.39222, -5.36591, -5.34383, -5.79111, -5.67660, -5.58924, -5.52166, -5.46879, -5.42707, -5.39394, -5.36751, -5.34637, -5.32946, -5.31596 ),length(t),6) weights <- c(0.00000, 0.21000, 0.02000, 0.55000, 0.00000, 0.22000) # model weights alpha <- 0.05 # confidence level y <- 304 # cut-off year for estimating probablity best.model <- 2 which.model <- 2 # specify the generating model#' # construct Studentized bootstrap confidence interval marp::student_confint( n,B,t,m,BB,par_hat,mu_hat,pr_hat,haz_hat,weights,alpha,y,best.model,which.model ) ## End(Not run)
Summarise a model-averaged renewal-process fit
## S3 method for class 'marp_fit' summary(object, ...)## S3 method for class 'marp_fit' summary(object, ...)
object |
An object returned by [marp()]. |
... |
Additional arguments, currently unused. |
A structured object of class 'summary_marp_fit'.
Summarise a model-specific renewal-process fit
## S3 method for class 'marp_model_fit' summary(object, ...)## S3 method for class 'marp_model_fit' summary(object, ...)
object |
An object returned by one of the model-specific fitting functions. |
... |
Additional arguments, currently unused. |
A structured object of class 'summary_marp_model_fit'.
Evaluate the upper-limit equation for a model-averaged T statistic
upperT(up, hat, sigmasq, Tstar, weights, B, alpha)upperT(up, hat, sigmasq, Tstar, weights, B, alpha)
up |
Candidate upper confidence limit. |
hat |
Vector of model-specific estimates. |
sigmasq |
Vector of model-specific variance estimates. |
Tstar |
Matrix of bootstrap T statistics, with models in rows. |
weights |
Vector of model weights. |
B |
Number of bootstrap samples represented in 'Tstar'. |
alpha |
Significance level. |
A numeric scalar giving the upper-limit root equation value.
# set some parameters up <- 100 # upper bound hat <- rep(150, 6) # estimates obtained from each model sigmasq <- 10 # variance Tstar <- matrix(rep(100,600),6,100) # T statistics estimated from bootstrap samples weights <- rep(1/6, 6) # model weights B <- 100 # number of bootstrapped samples alpha <- 0.05 # confidence level # calculate the upper limit of T statistics res <- marp::upperT(up, hat, sigmasq, Tstar, weights, B, alpha) res# set some parameters up <- 100 # upper bound hat <- rep(150, 6) # estimates obtained from each model sigmasq <- 10 # variance Tstar <- matrix(rep(100,600),6,100) # T statistics estimated from bootstrap samples weights <- rep(1/6, 6) # model weights B <- 100 # number of bootstrapped samples alpha <- 0.05 # confidence level # calculate the upper limit of T statistics res <- marp::upperT(up, hat, sigmasq, Tstar, weights, B, alpha) res
A function to generate (double) bootstrap samples and fit Weibull renewal model
weibull_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)weibull_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y)
n |
Number of inter-event times generated in each bootstrap sample. |
t |
Time points at which log-hazards are evaluated. |
B |
Number of bootstrap samples. |
BB |
Number of double-bootstrap samples per bootstrap sample. |
m |
Positive integer controlling repeated random-start optimizations; see [marp()]. |
par_hat |
Length-12 vector containing 'par1' for models 1–6 followed by 'par2' for models 1–6; see [marp()] for model order. |
mu_hat |
Length-6 vector of model-specific mean estimates. |
pr_hat |
Length-6 vector of model-specific logit event probabilities. |
haz_hat |
Matrix of model-specific log-hazards with 'length(t)' rows and six model columns. |
y |
Time point at which logit event probabilities are evaluated. |
A list of bootstrap estimates and variance/T-statistic quantities. Components beginning with 'pr_' use the logit-probability scale and those beginning with 'haz_' use the log-hazard scale.
Estimated mean from bootstrapped samples
Logit event probabilities from bootstrap samples
Log-hazards from bootstrap samples
Variance of estimated mean
Variance of estimated probability
Variance of estimated hazard rates
Variance of estimated mean of bootstrapped samples (via double-bootstrapping)
Variance of estimated probability of bootstrapped samples (via double-bootstrapping)
Variance of estimated hazard rates of bootstrapped samples (via double-bootstrapping)
Pivot quantity of the estimated mean
Pivot quantity of the estimated probability
Pivot quantity of the estimated hazard rates
## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.4136086430979953e-03, 2.7626793657057762e+00, 2.6037039674870583e+00, 3.3080162440951688e+02, 5.4882183788378658e+00, 2.9294512422957860e+02, NA, 9.4307059277139432e-03, 2.4759796859031687e+02, 1.8010183507666513e+00, 6.5084541680686814e-01, 7.1824719073918109e-01 ) mu_hat <- c( 292.94512187913182, 292.94512912200048, 319.72017228620746, 294.16945213908519, 298.87285747700128, 292.94512422957860 ) pr_hat <- c( 0.60038574701819891, 0.42154974433034809, 0.53433568234281148, 0.30779792692414687, 0.56416103510057725, 0.61794524610544410 ) haz_hat <- matrix(c( -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -6.0942031084732298, -5.9967873794574516, -5.9117418563554684, -5.8368230853439300, -5.7703089176306639, -5.7108525626839901, -5.6573839062669986, -5.6090408956082456, -5.5651206740587922, -5.5250440506799734, -5.4883291920475745, -6.0990192429336094, -5.9701664705134210, -5.8576899644670348, -5.7593884711134971, -5.6734972529860741, -5.5985621349393231, -5.5333565788683616, -5.4768259914915305, -5.4280496904694857, -5.3862145095364315, -5.3505961502861927, -6.1714638710963881, -6.0951186680582552, -6.0254209583640863, -5.9613052806725335, -5.9019434350392981, -5.8466788789061646, -5.7949823391436279, -5.7464209045603756, -5.7006359661738628, -5.6573271297614109, -5.6162402596857071, -5.9235521978533958, -5.8023896004395645, -5.7047473880293342, -5.6252373537796752, -5.5599409055534252, -5.5059486025117375, -5.4610610586440487, -5.4235891601883868, -5.3922173604047572, -5.3659081375131672, -5.3438339586221275, -5.7911126719889303, -5.6765973314326752, -5.5892417143301261, -5.5216608261560411, -5.4687921205249133, -5.4270729562323066, -5.3939387902533049, -5.3675067327627373, -5.3463701567645607, -5.3294619641245422, -5.3159614865560094 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::weibull_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)## Not run: # set some parameters n <- 30 # sample size t <- seq(100, 200, by = 10) # time intervals B <- 100 # number of bootstraps BB <- 100 # number of double-bootstraps m <- 10 # repeated random-start optimization setting par_hat <- c( 3.4136086430979953e-03, 2.7626793657057762e+00, 2.6037039674870583e+00, 3.3080162440951688e+02, 5.4882183788378658e+00, 2.9294512422957860e+02, NA, 9.4307059277139432e-03, 2.4759796859031687e+02, 1.8010183507666513e+00, 6.5084541680686814e-01, 7.1824719073918109e-01 ) mu_hat <- c( 292.94512187913182, 292.94512912200048, 319.72017228620746, 294.16945213908519, 298.87285747700128, 292.94512422957860 ) pr_hat <- c( 0.60038574701819891, 0.42154974433034809, 0.53433568234281148, 0.30779792692414687, 0.56416103510057725, 0.61794524610544410 ) haz_hat <- matrix(c( -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -5.6799852941338829, -6.0942031084732298, -5.9967873794574516, -5.9117418563554684, -5.8368230853439300, -5.7703089176306639, -5.7108525626839901, -5.6573839062669986, -5.6090408956082456, -5.5651206740587922, -5.5250440506799734, -5.4883291920475745, -6.0990192429336094, -5.9701664705134210, -5.8576899644670348, -5.7593884711134971, -5.6734972529860741, -5.5985621349393231, -5.5333565788683616, -5.4768259914915305, -5.4280496904694857, -5.3862145095364315, -5.3505961502861927, -6.1714638710963881, -6.0951186680582552, -6.0254209583640863, -5.9613052806725335, -5.9019434350392981, -5.8466788789061646, -5.7949823391436279, -5.7464209045603756, -5.7006359661738628, -5.6573271297614109, -5.6162402596857071, -5.9235521978533958, -5.8023896004395645, -5.7047473880293342, -5.6252373537796752, -5.5599409055534252, -5.5059486025117375, -5.4610610586440487, -5.4235891601883868, -5.3922173604047572, -5.3659081375131672, -5.3438339586221275, -5.7911126719889303, -5.6765973314326752, -5.5892417143301261, -5.5216608261560411, -5.4687921205249133, -5.4270729562323066, -5.3939387902533049, -5.3675067327627373, -5.3463701567645607, -5.3294619641245422, -5.3159614865560094 ),length(t),6) y <- 304 # cut-off year for estimating probablity # generate bootstrapped samples then fit renewal model res <- marp::weibull_bstrp(n, t, B, BB, m, par_hat, mu_hat, pr_hat, haz_hat, y) ## End(Not run)
A function to calculate the log-likelihood of Weibull model
weibull_logl(param, x)weibull_logl(param, x)
param |
Length-2 numeric vector containing log(scale) and log(shape). |
x |
Numeric vector of positive observations. |
returns the value of negative log-likelihood of the Weibull model
set.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(330.801103808081, 1.80101338777944) # estimated parameters param <- log(par_hat) # input parameters for logl function # calculate log-likelihood result <- marp::weibull_logl(param, data) resultset.seed(42) data <- rgamma(30,3,0.01) # set some parameters par_hat <- c(330.801103808081, 1.80101338777944) # estimated parameters param <- log(par_hat) # input parameters for logl function # calculate log-likelihood result <- marp::weibull_logl(param, data) result
A function to fit Weibull renewal model
weibull_rp(data, t, m, y)weibull_rp(data, t, m, y)
data |
A numeric vector of positive inter-event times. |
t |
A numeric vector of time points at which log-hazards are evaluated. |
m |
A positive integer controlling repeated random-start optimizations; the current implementation seeks 'm - 1' acceptable 'nlm()' fits. |
y |
A time point at which the logit-transformed cumulative event probability is evaluated. |
An object of class 'marp_model_fit' retaining the following eight named list components:
Estimated scale parameter of the Weibull model
Estimated shape parameter of the Weibull model
Maximized log-likelihood
Akaike information criterion (AIC)
Bayesian information criterion (BIC)
Estimated mean inter-event time
Logit-transformed cumulative event probability at 'y'
Log-hazard values at 't'
set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Weibull renewal model fit <- marp::weibull_rp(data, t, m, y) fit summary(fit)set.seed(42) data <- rgamma(100,3,0.01) # set some parameters m = 10 # repeated random-start optimization setting t = seq(100, 200, by=10) # time intervals y = 304 # cut-off year for estimating probablity # fit Weibull renewal model fit <- marp::weibull_rp(data, t, m, y) fit summary(fit)