parfm: Parametric Frailty Models
In parfm: Parametric Frailty Models

Description Usage Arguments Details Value Author(s) References See Also Examples

Fits parametric Frailty Models by maximizing the marginal likelihood.

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)

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

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) = 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) = phi ( (log(t) - mu) / sigma ) / {sigma * t * [1 - Phi( (log(t)- mu) / sigma ) )]}

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) = [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; ξ, ω, α) = \frac{2}{ω t} φ≤ft(\frac{\log(t) - ξ}{ω}\right) Φ≤ft(α\frac{\log(t)-ξ}{ω}\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) = [theta^(-1 / theta) * u^(1 / theta - 1) * exp(-u / theta)] / [Gamma(1 / theta)]

with theta > 0 and where Gamma(.) is the gamma function.

The inverse Gaussian frailty distribution (frailty="ingau") is

f(u) = 1 / sqrt(2 * pi * theta) * u^(-3 / 2) * exp[- (u - 1)^2 / (2 * theta * u)]

with theta > 0.

The positive stable frailty distribution (frailty="poosta") is

f(u) = - 1 / (pi * u) * sum_(k=1, 2, ..., infiny) [Gamma(k * (1 - nu) + 1) / k! * (-u^(nu - 1))^k * sin((1 - nu) * k * pi)]

with 0 < nu < 1.

The lognormal frailty distribution (frailty="lognormal") is

f(u) = 1 / sqrt(2 * pi * theta) * u^(-1) * exp[- (log u)^2 / (2 * theta)]

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

An object of class parfm.

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

Munda M, Rotolo F, Legrand C (2012). parfm: Parametric Frailty Models in R. Journal of Statistical Software, 51(11), 1-20. 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 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. URL http://www.jstor.org/stable/4615982

select.parfm, ci.parfm, predict.parfm

#------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="gompertz", frailty="logskewnormal")
#--------------------------

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