parfm {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)

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.

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).

Value

An object of class parfm.

Author(s)

Federico Rotolo [aut, cre], Marco Munda [aut], Andrea Callegaro [ctb]

References

Munda M, Rotolo F, Legrand C (2012). parfm: Parametric Frailty Models in R. Journal of Statistical Software, 51(11), 1-20. DOI <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 <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.

See Also

select.parfm, ci.parfm, predict.parfm

Examples

#------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")
      
      

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")
#--------------------------