parfm: Parametric Frailty Models

Description

Fits parametric Frailty Models by maximizing the marginal likelihood.

Usage

parfm(formula, cluster = NULL, strata = NULL, data,
      inip = NULL, iniFpar = NULL,
      dist = c("weibull", "inweibull", "frechet", "exponential", 
               "gompertz", "loglogistic", "lognormal", "logskewnormal"),
      frailty   = c("none", "gamma", "ingau", "possta",
                    "lognormal", "loglogistic"),
      method = "nlminb", 
      maxit = 500, Fparscale = 1, showtime = FALSE, correct = 0)

Value

An object of class parfm.

Arguments

formula: A formula object, with the response on the left of a ~ operator, and the terms on the right. The response must be a survival object as returned by the Surv() function. The status indicator 'event' in the Surv object must be 0=alive, 1=dead.
cluster: The name of a cluster variable in data.
strata: The name of a strata variable in data.
data: A data.frame in which to interpret the variables named in the formula.
inip: The vector of initial values. First components are for the baseline hazard parameters according to the order given in 'details'; Other components are for the regression parameters according to the order given in 'formula'.
iniFpar: The initial value of the frailty parameter.
dist: The baseline hazard distribution. One of weibull, inweibull, frechet, exponential, gompertz, lognormal, loglogistic, or logskewnormal.
frailty: The Frailty distribution. One of: none, gamma, ingau, possta or lognormal.
method: The optimisation method from the function optimx().
maxit: Maximum number of iterations (see optimx()).
Fparscale: the scaling value for the frailty parameter in optimx(). Optimisation is performed on Fpar/Fparscale.
showtime: Show the execution time?
correct: A correction factor that does not change the marginal log-likelihood except for an additive constant given by #clusters * correct * log(10). It may be useful in order to get finite log-likelihood values in case of many events per cluster with Positive Stable frailties. Note that the value of the log-likelihood in the output is the re-adjusted value.

Author

tools:::Rd_package_author("parfm")

Details

Baseline hazards

The Weibull hazard (dist="weibull") is $$h(t; \rho, \lambda) = \rho \lambda t^{\rho - 1}$$ with $\rho, \lambda > 0$.

The inverse Weibull (or Fréchet) hazard (dist="inweibull" or dist="frechet") is $$h(t; \rho, \lambda) = \frac{\rho \lambda t^{-\rho - 1}}{\exp(\lambda t^{-\rho}) - 1}$$ with $\rho, \lambda > 0$.

The exponential hazard (dist="exponential") is $$h(t; \lambda) = \lambda$$ with $\lambda > 0$.

The Gompertz hazard (dist="gompertz") is $$h(t; \gamma, \lambda) = \lambda \exp ( \gamma t )$$ with $\gamma, \lambda > 0$.

The lognormal hazard (dist="lognormal") is $$h(t; \mu, \sigma) = \frac{\phi \left( \frac{\log ( t ) - \mu}{\sigma}\right)}{\sigma t \left[ 1 - \Phi \left( \frac{\log ( t ) - \mu}{\sigma}\right)\right]}$$ with $\mu \in {R}, \sigma > 0$, and where $\phi$ and $\Phi$ are the density and distribution functions of a standard Normal random variable.

The loglogistic hazard (dist="loglogistic") is $$h(t; \alpha, \kappa) = \frac{\exp ( \alpha ) \kappa t^{\kappa - 1}}{1 + \exp ( \alpha ) t^{\kappa}}$$ with $\alpha \in {R}, \kappa > 0$.

The log-skew-normal hazard (dist="logskewnormal") is obtained as the ratio between the density and the cumulative distribution function of a log-skew normal random variable (Azzalini, 1985), which has density $$f(t; \xi, \omega, \alpha) = \frac{2}{\omega t} \phi\left(\frac{\log(t) - \xi}{\omega}\right) \Phi\left(\alpha\frac{\log(t)-\xi}{\omega}\right)$$ with $\xi \in {R}, \omega > 0, \alpha \in {R}$, and where $\phi$ and $\Phi$ are the density and distribution functions of a standard Normal random variable. Of note, if $\alpha=0$ then the log-skew-normal boils down to the log-normal distribution, with $\mu=\xi$ and $\sigma=\omega$.

Frailty distributions

The gamma frailty distribution (frailty="gamma") is $$f ( u ) = \frac{\theta^{-\frac{1}{\theta}} u^{\frac{1}{\theta} - 1} \exp \left( - u / \theta \right)} {\Gamma ( 1 / \theta )}$$ with $\theta > 0$ and where $\Gamma ( \cdot )$ is the gamma function.

The inverse Gaussian frailty distribution (frailty="ingau") is $$f(u) = \frac1{\sqrt{2 \pi \theta}} u^{- \frac32} \exp \left( - \frac{(u-1)^2}{2 \theta u} \right)$$ with $\theta > 0$.

The positive stable frailty distribution (frailty="poosta") is $$f(u) = f(u) = - \frac1{\pi u} \sum_{k=1}^{\infty} \frac{\Gamma ( k (1 - \nu ) + 1 )}{k!} \left( - u^{ \nu - 1} \right)^{k} \sin ( ( 1 - \nu ) k \pi )$$ with $0 < \nu < 1$.

The lognormal frailty distribution (frailty="lognormal") is $$f(u) = \frac1{\sqrt{2 \pi \theta}} u^{-1} \exp \left( - \frac{\log(u)^2}{2 \theta} \right)$$ with $\theta > 0$. As the Laplace tranform of the lognormal frailties does not exist in closed form, the saddlepoint approximation is used (Goutis and Casella, 1999).

References

Munda M, Rotolo F, Legrand C (2012). parfm: Parametric Frailty Models in R. Journal of Statistical Software, 51(11), 1-20. DOI <tools:::Rd_expr_doi("10.18637/jss.v051.i11")>

Duchateau L, Janssen P (2008). The frailty model. Springer.

Wienke A (2010). Frailty Models in Survival Analysis. Chapman & Hall/CRC biostatistics series. Taylor and Francis.

Goutis C, Casella G (1999). Explaining the Saddlepoint Approximation. The American Statistician 53(3), 216-224. DOI <tools:::Rd_expr_doi("10.1080/00031305.1999.10474463")>

Azzalini A (1985). A class of distributions which includes the normal ones. Scandinavian Journal of Statistics, 12(2):171-178.

Examples

Run this code

#------Kidney dataset------
data(kidney) 
 # type 'help(kidney)' for a description of the data set
kidney$sex <- kidney$sex - 1

parfm(Surv(time,status) ~ sex + age, cluster = "id",
      data = kidney, dist = "exponential", frailty = "gamma")
      
      
# \donttest{
parfm(Surv(time,status) ~ sex + age, cluster = "id",
      data = kidney, dist = "exponential", frailty = "lognormal")

parfm(Surv(time,status) ~ sex + age, cluster = "id",
      data = kidney, dist = "weibull", frailty = "ingau")

parfm(Surv(time,status) ~ sex + age, cluster = "id",
      data = kidney, dist="gompertz", frailty="possta", method="CG")


parfm(Surv(time,status) ~ sex + age, cluster = "id",
      data = kidney, dist="logskewnormal", frailty="possta", method = 'BFGS')
#--------------------------

#------Asthma dataset------
data(asthma)
head(asthma)
  # type 'help(asthma)' for a description of the data set

asthma$time <- asthma$End - asthma$Begin
parfm(Surv(time, Status) ~ Drug, cluster = "Patid", data = asthma,
      dist = "weibull", frailty = "gamma")
      
parfm(Surv(time, Status) ~ Drug, cluster = "Patid", data = asthma,
      dist = "weibull", frailty = "lognormal")

parfm(Surv(Begin, End, Status) ~ Drug, cluster = "Patid", 
      data = asthma[asthma$Fevent == 0, ],
      dist = "weibull", frailty = "lognormal", method = "Nelder-Mead")
#--------------------------
# }

Run the code above in your browser using DataLab