# parfm: Parametric Frailty Models In parfm: Parametric Frailty Models

## Description

Fits parametric Frailty Models by maximizing the marginal likelihood.

## Usage

 1 2 3 4 5 6 7 8 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) = 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).

## 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 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
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 #------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") #--------------------------