Fitting dynamic frailty models with the EM algorithm
dynfrail(formula, data, distribution = dynfrail_dist(),
control = dynfrail_control(), ...)A formula that contains on the left hand side an object of the type Surv
and on the right hand side a +cluster(id) statement.
A data frame in which the formula argument can be evaluated
An object as created by dynfrail_dist
An object as created by dynfrail_control
Other arguments, currently used to warn about deprecated argument names
A dynfrail object that contains the following fields:
A named vector of the estimated regression coefficients
The breslow estimate of the baseline hazard at each event time point, in chronological order
Fisher's information matrix corresponding to the coefficients and hazard, assuming \(\theta, \lambda\) constant
The point estimate of the logarithm of the frailty parameter \(\theta\). See details.
The point estimate of the logarithm of the autocorrelation parameter \(\lambda\). See details.
A data.frame containing the variables: id (cluster id), interval (for piecewise constant frailty, the label of the interval
on which the frailty is constant), Y (a Surv object which determines a starting and a stopping time for each row), frail (the empirical
Bayes estimates of the piecewise constant frailty corresponding to that specific cluster and that specific time period)
The time points of the events in the data set, this is the same length as hazard
A vector of length two with the log-likelihood of the starting Cox model and the maximized log-likelihood
The original formula argument
The original distribution argument
The original control argument
This function fits dynamic frailty models where the intensity of the process is described by $$\lambda(t) = Z(t) \exp(\beta^\top x) \lambda_0(t).$$ As in regular frailty models, the random effect is shared by observations from a cluster, or by recurrent event episodes within an individual. This implementation generally follows the lines of Putter & van Houwelingen (2015). The maximum likelihood estimates are obtained with an exact E step.
\(Z(t)\) has two parameters: \(\theta\) plays the role of the spread of the frailty distribution. For the frailty distributons with finite
variance (all except the positive stable) this is the inverse of the variance, so that 0 corresponds to infinite variance and infinity to
variance 0. The second parameter \(\lambda\) determines how much variation in time is in \(Z(t)\), so that
$$cor(Z(t_1), Z(t_2)) = exp(-\lambda (t_2 - t_1)).$$ Note that this heavily depends on the time scale, so
the starting value in the distribution should reflect that.
By default, the program must calculate \(Z(t)\) for each cluster and for each event time point in the data.
This is computationally challenging. An option is to use the nints argument in the control argument.
This considers \(Z(t)\) to be piecewise constant over nints + 1 intervals. These intervals are determined
automatically so that there are roughly an equal number of observations for each interval. Using nints = 0
is equivalent to fitting a shared frailty model with the frailtyEM package.
It is recommended that the user starts with nints = 0 and then slowly increase the number of intervals.
Other options for performance may be set within the control argument. Also, this could be tried out first
on a subset of the data.
For computational reasons, the standard errors of \(\theta\) and \(\lambda\) are not calculated, and the standard errors of the regression coefficients are obtained under the assumption that the frailty distribution is fixed.
Putter, H., & Van Houwelingen, H. C. (2015). Dynamic frailty models based on compound birth<U+2013>death processes. Biostatistics, 16(3), 550-564.
# NOT RUN {
# 5 piecewise constant intervals
m2 <- dynfrail(Surv(time, status) ~ rx + sex + cluster(litter),
data = rats,
distribution = dynfrail_dist(n_ints = 4))
# }
# NOT RUN {
#' # essentially a gamma frailty fit
m1 <- dynfrail(Surv(time, status) ~ rx + sex + cluster(litter),
data = rats,
distribution = dynfrail_dist(n_ints = 0))
# completely semiparametric gamma frailty
m2 <- dynfrail(Surv(time, status) ~ rx + sex + cluster(litter),
data = rats)
# }
Run the code above in your browser using DataLab