mlefrailty_fit: Survival estimator under a gamma frailty model
In survrec: Survival Analysis for Recurrent Event Data

mlefrailty_fit

R Documentation

Survival estimator under a gamma frailty model

Description

Maximum likelihood estimation of the survival function of correlated inter-occurrence times of recurrent event data under a gamma frailty model.

Usage

mlefrailty_fit(
  x,
  tvals,
  lambda = NULL,
  alpha = NULL,
  alpha.min,
  alpha.max,
  tol = 1e-07,
  maxiter = 500,
  alpha.console = TRUE
)

Arguments

`x`	a survival recurrent event object, see `Survr()`.
`tvals`	optional vector of times at which the survival function is also evaluated.
`lambda`	optional vector of baseline hazard probabilities at the distinct event times. Defaults to the occurrence/exposure rates `n.event/colSums(AtRisk)`.
`alpha`	optional shape and scale parameter of the frailty distribution. If unknown it is estimated via the EM algorithm, starting from a seed obtained by maximising the profile likelihood (see details).
`alpha.min`	optional left bound of the interval used to find the seed of `alpha`. Default is 0.5.
`alpha.max`	optional right bound of the interval used to find the seed of `alpha`. Default is the maximum distinct event time.
`tol`	convergence tolerance of the EM algorithm. Default is 1e-7.
`maxiter`	maximum number of EM iterations. Default is 500.
`alpha.console`	if `TRUE`, prints the seed value and the final estimate of `alpha`.

Details

The product-limit estimator of Peña, Strawderman and Hollander (2001) is valid when the inter-occurrence times are an IID sample from some underlying distribution F. This assumption is clearly restrictive in biomedical applications, and one obvious generalization that allows association between inter-occurrence times is a frailty model.

A common and convenient choice of frailty distribution is a gamma with shape and scale parameters set equal to an unknown parameter \alpha. The common marginal survival function can be written as

\bar{F}(t) = \left[\frac{\alpha}{\alpha + \Lambda_0(t)} \right]^{\alpha}

The parameter \alpha controls the degree of association between inter-occurrence times within a unit. Peña, Strawderman and Hollander (2001) showed that \alpha and \Lambda_0 can be estimated by maximising the marginal likelihood function with an expectation-maximisation (EM) algorithm.

To achieve good convergence, a seed value for \alpha is estimated first by maximising the profile likelihood for \alpha over the interval (alpha.min, alpha.max) with Brent's method. If the EM algorithm does not converge from that seed, it is restarted from up to six additional seeds around it. If convergence still fails, the initial value of alpha can be used as the alpha.min argument and the fit recomputed.

Value

A list with class survfitr containing:

n: number of units or subjects observed.
m: number of events of each subject.
failed: uncensored gap times, ordered by subject.
censored: censored gap time of each subject.
time: ordered distinct event times.
n.event: number of events at each distinct time.
AtRisk: matrix of units at risk at each distinct time (rows are subjects).
status: 1 if the EM algorithm converged, -1 if it reached maxiter without converging (in that case no estimates are provided).
alpha: estimate of the gamma frailty parameter.
lambda: estimated hazard probabilities at the distinct event times.
frailties: posterior frailty estimate of each subject; values spread away from 1 indicate heterogeneity between subjects.
survfunc: survival estimate at each distinct time.
tvals: copy of the tvals argument.
MLEAttvals: survival estimate at the tvals times.

References

Peña, E.A., Strawderman, R. and Hollander, M. (2001). Nonparametric Estimation with Recurrent Event Data. Journal of the American Statistical Association 96, 1299–1315.

Examples

data(MMC)
fit <- mlefrailty_fit(Survr(MMC$id, MMC$time, MMC$event))
fit
plot(fit)

# compare with Peña-Strawderman-Hollander
fit <- psh_fit(Survr(MMC$id, MMC$time, MMC$event))
lines(fit, lty = 2)

# and with Wang-Chang
fit <- wc_fit(Survr(MMC$id, MMC$time, MMC$event))
lines(fit, lty = 3)

survrec documentation built on Sept. 22, 2026, 9:07 a.m.