GenfrailtyPenal: Fit a Shared or a Joint Frailty Generalized Survival Model
In frailtypack: Shared, Joint (Generalized) Frailty Models; Surrogate Endpoints
GenfrailtyPenal R Documentation
Fit a Shared or a Joint Frailty Generalized Survival Model
Description
I. SHARED FRAILTY GENERALIZED SURVIVAL MODELS
Fit a gamma Shared Frailty Generalized Survival Model using
a parametric estimation, or a semi-parametric penalized likelihood estimation.
Right-censored data and strata (up to 6 levels) are allowed.
It allows to obtain a parametric or flexible semi-parametric smooth
hazard and survival functions.
Each frailty term ui is assumed
to act multiplicatively on the hazard function, and to be drawn from a
Gamma distribution with unit mean and variance \theta.
Conditional on the frailty term, the hazard function for the
jth subject in the ith group
is then expressed by
where xij
is a collection of baseline covariates,
\xi is a vector of parameters, and
\lambdaij
(t | xij ; \xi)
is the hazard function for an average value of the frailty.
The associated conditional survival function writes
where
Sij
(t | xij ; \xi)
designates the survival function for an average value of the frailty.
Following Liu et al. (2017, 2018), the latter function is expressed in terms of
a link function g(.) and a linear predictor
\etaij
(t, xij; \xi)
such that
g[Sij
(t | xij ; \xi)]
=
\etaij
(t, xij; \xi),
i.e.
Sij
(t | xij ; \xi)
=
h[\etaij
(t, xij; \xi)]
with h() = g-1().
The conditional survival function is finally modeled by
The table below summarizes the most commonly used (inverse) link functions and
their associated conditional survival, hazard and cumulative hazard functions.
PHM stands for "Proportional Hazards Model",
POM for "Proportional Odds Model,
PROM for "Probit Model" and AHM for "Additive Hazards Model".
I.(a) Fully parametric case
In the fully parametric case, linear predictors considered are of the form
where \rho > 0 is a shape parameter,
\gamma > 0 a scale parameter,
\beta a vector of regression coefficients,
and \xi = (\rho ,\gamma, \beta).
With the appropriate link function, such linear parametric predictors
make it possible to recover
a Weibull baseline survival function for PHMs and AHMs,
a log-logistic baseline survival function for POMs,
and a log-normal one for PROMs.
I. (b) Flexible semi-parametric case
For PHM and AHM, a more flexible splines-based approach is proposed for
modeling the baseline hazard function and time-varying regression coefficients.
In this case, conditional on the frailty term ui,
the hazard function for the jth subject
in the ith group is still expressed by
\lambdaij
(t | xij, ui ;
\xi)
= ui
\lambdaij
(t | xij ; \xi),
but we have this time
The smoothness of baseline hazard function \lambda0()
is ensured by penalizing the log-likelihood by a term which has
large values for rough functions.
Moreover, for parametric and flexible semi-parametric AHMs, the
log-likelihood is constrained to ensure the strict positivity of the hazards,
since the latter is not naturally guaranteed by the model.
II. JOINT FRAILTY GENERALIZED SURVIVAL MODELS
Fit a gamma Joint Frailty Generalized Survival Model for recurrent and
terminal events using a parametric estimation,
or a semi-parametric penalized likelihood estimation.
Right-censored data and strata (up to 6 levels) for the recurrent event part
are allowed.
Joint frailty models allow studying, jointly, survival processes
of recurrent and terminal events, by considering the terminal event as an
informative censoring.
This model includes an common patient-specific frailty term
ui for the two survival functions which will
take into account the unmeasured heterogeneity in the data,
associated with unobserved covariates.
The frailty term acts differently for the two survival functions
(ui for the recurrent survival function and
uiα for the death one).
The covariates could be different for the recurrent and terminal event parts.
II.(a) Fully parametric case
For the jth recurrence (j=1,...,ni)
and the ith patient (i=1,...,N),
the gamma Joint Frailty Generalized Survival Model
for recurrent event survival function
SRij(.) and death survival function
SDi(.) is
- \etaRij (resp. \etaDi)
is the linear predictor for the recurrent (resp. terminal) event process.
The form of these linear predictors is the same as the one presented in I.(a).
- hR(.) (resp. hD(.))
is the inverse link function associated with
recurrent events (resp. terminal event).
- xRij and xDi
are two vectors of baseline covariates associated with
recurrent and terminal events.
- \xiR and \xiD
are the parameter vectors for recurrent and terminal events.
- \alpha is a parameter allowing more flexibility in the association
between recurrent and terminal events processes.
- The random frailties ui are still assumed iid and
drown from a \Gamma(1/\theta,1/\theta).
II.(b) Flexible semi-parametric case
If one chooses to fit a PHM or an AHM for recurrent and/or terminal events,
a splines-based approach for modeling baseline hazard functions
and time-varying regression coefficients is still available.
In this approach, the submodel for recurrent events is expressed as
\lambdaRij
(t | xRij, ui ;
\xiR)
= ui
\lambdaRij
(t | xRij ;
\xiR), where
The submodel for terminal event is expressed as
\lambdaDi
(t | xDi, ui ;
\xiD)
= uiα
\lambdaDi
(t | xDi ;
\xiD), where
Baseline hazard functions
\lambdaR0(.) and \lambdaD0(.)
are estimated using cubic M-splines (of order 4)
with positive coefficients, and the time-varying coefficients
\betaR(.) and \betaD(.)
are estimated using B-splines of order q.
The smoothness of baseline hazard functions
is ensured by penalizing the log-likelihood by two terms
which has large values for rough functions.
Moreover,
if one chooses an AHM for recurrent and/or terminal event submodel,
the log-likelihood is constrained to ensure
the strict positivity of the hazards,
since the latter is not naturally guaranteed by the model.
Usage
GenfrailtyPenal(formula, formula.terminalEvent, data, recurrentAG = FALSE,
family, hazard = "Splines", n.knots, kappa, betaknots = 1, betaorder = 3,
RandDist = "Gamma", init.B, init.Theta, init.Alpha, Alpha, maxit = 300,
nb.gh, nb.gl, LIMparam = 1e-3, LIMlogl = 1e-3, LIMderiv = 1e-3, print.times = TRUE,
cross.validation, jointGeneral, nb.int, initialize, init.Ksi, Ksi, init.Eta)
Arguments
formula
A formula object, with the response on the left of a
\sim operator, and the terms on the right. The response must be a
survival object as returned by the 'Surv' function
like in survival package. Interactions are possible using ' * ' or ' : '.
formula.terminalEvent
Only for joint frailty models: a formula object,
only requires terms on the right to indicate which variables are used for
the terminal event. Interactions are possible using ' * ' or ' : '.
data
A 'data.frame' with the variables used in 'formula'.
recurrentAG
Logical value. Is Andersen-Gill model fitted? If so
indicates that recurrent event times with the counting process approach of
Andersen and Gill is used. This formulation can be used for dealing with
time-dependent covariates. The default is FALSE.
family
Type of Generalized Survival Model to fit.
"PH" for a proportional hazards model,
"AH" for an additive hazards model,
"PO" for a proportional odds model and
"probit" for a probit model.
A vector of length 2 is expected for joint models
(e.g., family=c("PH","PH")).
hazard
Type of hazard functions:
"Splines" for semi-parametric hazard functions using equidistant
intervals, or "parametric" for parametric distribution functions.
In case of family="PH" or family="AH",
the "parametric" option corresponds to a Weibull distribution.
In case of family="PO" and family="probit",
the "parametric" option corresponds to a log-logistic
and a log-normal distribution, respectively.
So far, the "Splines" option is only available for PH and AH submodels.
Default is "Splines".
n.knots
Integer giving the number of knots to use. Value required in
the penalized likelihood estimation. It corresponds to the n.knots+2
splines functions for the approximation of the hazard or the survival
functions. We estimate I- or M-splines of order 4. When the user set a
number of knots equals to k (i.e. n.knots=k),
then the number of interior knots is k-2 and the number of splines is
(k-2)+order. Number of knots must be between 4 and 20. (See Note)
kappa
Positive smoothing parameter in the penalized likelihood
estimation. The coefficient kappa tunes the intensity of the penalization
(the integral of the squared second derivative of hazard function).
In a stratified shared model, this argument must be a vector
with kappas for both strata. In a stratified joint model, this argument
must be a vector with kappas for both strata for recurrent events plus one
kappa for terminal event.
We advise the user to identify several possible tuning
parameters, note their defaults and look at the sensitivity of the results
to varying them. Value required. (See Note).
betaknots
Number of inner knots used for the
B-splines time-varying coefficient estimation. Default is 1.
See 'timedep' function for more details.
betaorder
Order of the B-splines used for the
time-varying coefficient estimation.
Default is cubic B-splines (order=3).
See 'timedep' function for more details.
Not implemented for Proportional Odds and Probit submodels.
RandDist
Type of random effect distribution:
"Gamma" for a gamma distribution,
and "LogN" for a log-normal distribution (not implemented yet).
Default is "Gamma".
init.B
A vector of initial values for regression coefficients. This
vector should be of the same size as the whole vector of covariates with the
first elements for the covariates related to the recurrent events and then
to the terminal event (interactions in the end of each component). Default
is 0.1 for each (for Generalized Survival and Shared Frailty Models)
or 0.5 (for Generalized Joint Frailty Models).
init.Theta
Initial value for frailty variance.
init.Alpha
Only for Generalized Joint Frailty Models:
initial value for parameter alpha.
Alpha
Only for Generalized Joint Frailty Models:
input "None" so as to fit a joint model without the parameter alpha.
maxit
Maximum number of iterations for the Marquardt algorithm.
Default is 300
nb.gh
Number of nodes for the Gaussian-Hermite quadrature.
It can be chosen among 5, 7, 9, 12, 15, 20 and 32.
The default is 20 if hazard="Splines", 32 otherwise.
nb.gl
Number of nodes for the Gaussian-Laguerre quadrature.
It can be chosen between 20 and 32.
The default is 20 if hazard="Splines", 32 otherwise.
LIMparam
Convergence threshold of the Marquardt algorithm for the
parameters (see Details), 10-3 by default.
LIMlogl
Convergence threshold of the Marquardt algorithm for the
log-likelihood (see Details), 10-3 by default.
LIMderiv
Convergence threshold of the Marquardt algorithm for the
gradient (see Details), 10-3 by default.
print.times
A logical parameter to print iteration process. Default
is TRUE.
cross.validation
Not implemented yet for the generalized settings.
jointGeneral
Not implemented yet for the generalized settings.
nb.int
Not implemented yet for the generalized settings.
initialize
Not implemented yet for the generalized settings.
init.Ksi
Not implemented yet for the generalized settings.
Ksi
Not implemented yet for the generalized settings.
init.Eta
Not implemented yet for the generalized settings.
Details
TYPICAL USES
For a Generalized Survival Model:
GenfrailtyPenal(
formula=Surv(time,event)~var1+var2,
data, family, \dots)
For a Shared Frailty Generalized Survival Model:
GenfrailtyPenal(
formula=Surv(time,event)~cluster(group)+var1+var2,
data, family, \dots)
For a Joint Frailty Generalized Survival Model:
GenfrailtyPenal(
formula=Surv(time,event)~cluster(group)+var1+var2+var3+terminal(death),
formula.terminalEvent= ~var1+var4,
data, family, \dots)
OPTIMIZATION ALGORITHM
The estimated parameters are obtained using the robust Marquardt algorithm
(Marquardt, 1963) which is a combination between a Newton-Raphson algorithm
and a steepest descent algorithm. The iterations are stopped when
the difference between two consecutive log-likelihoods is small
(<10-3),
the estimated coefficients are stable
(consecutive values <10-3,
and the gradient small enough (<10-3).
When the frailty variance is small, numerical problems may arise.
To solve this problem, an alternative formula of the penalized log-likelihood
is used (see Rondeau, 2003 for further details).
For Proportional Hazards and Additive Hazards submodels,
cubic M-splines of order 4 can be used to estimate the hazard function.
In this case, I-splines (integrated M-splines) are used to compute the
cumulative hazard function.
The inverse of the Hessian matrix is the variance estimator.
To deal with the positivity constraint of the variance component and the
spline coefficients, a squared transformation is used and the standard errors
are computed by the \Delta-method (Knight & Xekalaki, 2000).
The integrations in the full log likelihood are evaluated using
Gaussian quadrature. Laguerre polynomials with 20 points are used to treat
the integrations on [0,\infty[.
INITIAL VALUES
In case of a shared frailty model,
the splines and the regression coefficients are initialized to 0.1.
The program fits, firstly, an adjusted Cox model to give new initial values
for the splines and the regression coefficients.
The variance of the frailty term \theta is initialized to 0.1.
Then, a shared frailty model is fitted.
In case of a joint frailty model,
the splines and the regression coefficients are initialized to 0.5.
The program fits firstly, an adjusted Cox model to have new initial values
for the splines and the regression coefficients.
The variance of the frailty term \theta and the association parameter
\alpha are initialized to 1.
Then, a joint frailty model is fitted.
Value
The following components are included in a 'frailtyPenal' object for each
model.
b
Sequence of the corresponding estimation of the coefficients for
the hazard functions (parametric or semiparametric), the random effects
variances and the regression coefficients.
call
The code used for the model.
formula
The formula part of the code used for the model.
n
The number of observations used in the fit.
groups
The maximum number of groups used in the fit.
n.events
The number of events observed in the fit.
n.eventsbygrp
A vector of length the number of groups
giving the number of observed events in each group.
loglik
The marginal log-likelihood in the parametric case.
loglikPenal
The marginal penalized log-likelihood
in the semiparametric case.
coef
The regression coefficients.
varH
The variance matrix of
the regression coefficients before positivity constraint transformation.
Then, the delta method is needed to obtain the estimated variance parameters.
That is why some variances don't match with
the printed values at the end of the model.
varHtotal
The variance matrix of
all the parameters before positivity constraint transformation.
Then, the delta method is needed to obtain the estimated variance parameters.
That is why some variances don't match with
the printed values at the end of the model.
varHIH
The robust estimation of
the variance matrix of the regression coefficients
varHIHtotal
The robust estimation of
the variance matrix of all parameters.
x
Matrix of times where the hazard functions are estimated.
xSu
Matrix of times where the survival functions are estimated.
lam
Array (dim=3) of baseline hazard estimates
and confidence bands.
surv
Array (dim=3) of baseline survival estimates
and confidence bands.
type
Character string specifying the type of censoring,
see the Surv function for more details.
n.strat
Number of strata.
n.iter
Number of iterations needed to converge.
median
The value of the median survival and its confidence bands.
If there are two strata or more, the first value corresponds to the value
for the first strata, etc.
LCV
The approximated likelihood cross-validation criterion in the
semiparametric case.
With H (resp. Hpen) the hessian matrix
of log-likelihood (resp. penalized log-likelihood),
EDF = Hpen-1 H
the effective degrees of freedom,
L(\xi,\theta) the log-likelihood and
n the number of observations,
LCV = 1/n x (trace(EDF) - L(\xi,\theta)).
AIC
The Akaike information Criterion for the parametric case.
With p the number of parameters,
n the number of observations and L(\xi,\theta) the log-likelihood,
AIC = 1/n x (p - L(\xi,\theta)).
npar
Number of parameters.
nvar
Number of explanatory variables.
typeof
Indicator of the type of hazard functions computed :
0 for "Splines", 2 for "parametric".
istop
Convergence indicator:
1 if convergence is reached,
2 if convergence is not reached,
3 if the hessian matrix is not positive definite,
4 if a numerical problem has occurred in the likelihood calculation
shape.param
Shape parameter for the parametric hazard function
(a Weibull distribution is used for proportional and additive hazards models,
a log-logistic distribution is used for proportional odds models,
a log-normal distribution is used for probit models).
scale.param
Scale parameter for the parametric hazard function.
Names.data
Name of the dataset.
Frailty
Logical value. Was model with frailties fitted ?
linear.pred
Linear predictor:
\beta'X in the generalized survival models or
\beta'X + log(ui)
in the shared frailty generalized survival models.
BetaTpsMat
Matrix of time varying-effects and confidence bands
(the first column used for abscissa of times).
nvartimedep
Number of covariates with time-varying effects.
Names.vardep
Name of the covariates with time-varying effects.
EPS
Convergence criteria concerning
the parameters, the likelihood and the gradient.
family
Type of Generalized Survival Model fitted
(0 for PH, 1 for PO, 2 for probit, 3 for AH).
global_chisq.test
A binary variable equals to 0 when no multivariate
Wald is given, 1 otherwise.
beta_p.value
p-values of the Wald test for the estimated
regression coefficients.
cross.Val
Logical value. Is cross validation procedure used for
estimating the smoothing parameters in the penalized likelihood estimation?
DoF
Degrees of freedom associated with the smoothing parameter
kappa.
kappa
A vector with the smoothing parameters in the penalized
likelihood estimation corresponding to each baseline function as components.
n.knots
Number of knots for estimating the baseline functions in the
penalized likelihood estimation.
n.knots.temp
Initial value for the number of knots.
global_chisq
A vector with the values of each multivariate Wald test.
dof_chisq
A vector with the degree of freedom for each multivariate
Wald test.
p.global_chisq
A vector with the p-values for each global multivariate
Wald test.
names.factor
Names of the "as.factor" variables.
Xlevels
Vector of the values that factor might have taken.
The following components are specific to shared models.
equidistant
Indicator for the intervals used in the spline estimation
of baseline hazard functions :
1 for equidistant intervals ; 0 for intervals using percentile
(note: equidistant = 2 in case of parametric estimation).
Names.cluster
Cluster names.
theta
Variance of the gamma frailty parameter, i.e.
Var(ui).
varTheta
Variance of parameter theta.
theta_p.value
p-value of the Wald test for
the estimated variance of the gamma frailty.
The following components are specific to joint models.
formula
The formula part of the code
used for the recurrent events.
formula.terminalEvent
The formula part of the code
used for the terminal model.
n.deaths
Number of observed deaths.
n.censored
Number of censored individuals.
theta
Variance of the gamma frailty parameter, i.e.
Var(ui).
indic_alpha
Indicator if a joint frailty model with
\alpha parameter was fitted.
alpha
The coefficient \alpha associated
with the frailty parameter in the terminal hazard function.
nvar
A vector with the number of covariates
of each type of hazard function as components.
nvarnotdep
A vector with the number of constant effect covariates
of each type of hazard function as components.
nvarRec
Number of recurrent explanatory variables.
nvarEnd
Number of death explanatory variables.
noVar1
Indicator of recurrent explanatory variables.
noVar2
Indicator of death explanatory variables.
Names.vardep
Name of the covariates with time-varying effects
for the recurrent events.
Names.vardepdc
Name of the covariates with time-varying effects
for the terminal event.
xR
Matrix of times where both survival and hazard function
are estimated for the recurrent event.
xD
Matrix of times for the terminal event.
lamR
Array (dim=3) of hazard estimates and confidence bands
for recurrent event.
lamD
The same value as lamR for the terminal event.
survR
Array (dim=3) of baseline survival estimates and
confidence bands for recurrent event.
survD
The same value as survR for the terminal event.
nb.gh
Number of nodes for the Gaussian-Hermite quadrature.
nb.gl
Number of nodes for the Gaussian-Laguerre quadrature.
medianR
The value of the median survival for the recurrent events
and its confidence bands.
medianD
The value of the median survival for the terminal event
and its confidence bands.
names.factor
Names of the "as.factor" variables
for the recurrent events.
names.factordc
Names of the "as.factor" variables
for the terminal event.
Xlevels
Vector of the values that factor might have taken
for the recurrent events.
Xlevels2
Vector of the values that factor might have taken
for the terminal event.
linear.pred
Linear predictor for the recurrent part:
\beta'X + log(ui).
lineardeath.pred
Linear predictor for the terminal part:
\beta'X + \alpha x log(ui).
Xlevels
Vector of the values that factor might have taken
for the recurrent part.
Xlevels2
vector of the values that factor might have taken
for the death part.
BetaTpsMat
Matrix of time varying-effects and confidence bands for
recurrent event (the first column used for abscissa of times of recurrence).
BetaTpsMatDc
Matrix of time varying-effects and confidence bands for
terminal event (the first column used for abscissa of times of death).
alpha_p.value
p-value of the Wald test for the estimated \alpha.
Note
In the flexible semiparametric case, smoothing parameters kappa and
number of knots n.knots are the arguments that the user have to change
if the fitted model does not converge.
n.knots takes integer values between 4 and 20.
But with n.knots=20, the model would take a long time to converge.
So, usually, begin first with n.knots=7, and increase it step by step
until it converges.
kappa only takes positive values. So, choose a value for kappa (for
instance 10000), and if it does not converge, multiply or divide this value
by 10 or 5 until it converges.
References
J. Chauvet and V. Rondeau (2021). A flexible class of generalized
joint frailty models for the analysis of survival endpoints. In revision.
Liu XR, Pawitan Y, Clements M. (2018)
Parametric and penalized generalized survival models.
Statistical Methods in Medical Research 27(5), 1531-1546.
Liu XR, Pawitan Y, Clements MS. (2017)
Generalized survival models for correlated time-to-event data.
Statistics in Medicine 36(29), 4743-4762.
A. Krol, A. Mauguen, Y. Mazroui, A. Laurent, S. Michiels and V. Rondeau
(2017). Tutorial in Joint Modeling and Prediction: A Statistical Software
for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event.
Journal of Statistical Software 81(3), 1-52.
V. Rondeau, Y. Mazroui and J. R. Gonzalez (2012). Frailtypack: An R package
for the analysis of correlated survival data with frailty models using
penalized likelihood estimation or parametric estimation. Journal of
Statistical Software 47, 1-28.
V. Rondeau, J.P. Pignon, S. Michiels (2011). A joint model for the
dependance between clustered times to tumour progression and deaths: A
meta-analysis of chemotherapy in head and neck cancer. Statistical
methods in medical research 897, 1-19.
V. Rondeau, S. Mathoulin-Pellissier, H. Jacqmin-Gadda, V. Brouste, P.
Soubeyran (2007). Joint frailty models for recurring events and death using
maximum penalized likelihood estimation:application on cancer events.
Biostatistics 8,4, 708-721.
V. Rondeau, D. Commenges, and P. Joly (2003). Maximum penalized likelihood
estimation in a gamma-frailty model. Lifetime Data Analysis 9,
139-153.
C.A. McGilchrist, and C.W. Aisbett (1991). Regression with frailty in
survival analysis. Biometrics 47, 461-466.
D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear
parameters. SIAM Journal of Applied Mathematics, 431-441.
See Also
Surv,
terminal,
timedep
Examples
#############################################################################
# ----- GENERALIZED SURVIVAL MODELS (without frailties) ----- #
#############################################################################
adult.retino = retinopathy[retinopathy$type == "adult", ]
adult.retino[adult.retino$futime >= 50, "status"] = 0
adult.retino[adult.retino$futime >= 50, "futime"] = 50
### --- Parametric PH, AH, PO and probit models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="PH")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="AH")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="PO")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="probit")
### --- Semi-parametric PH and AH models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino,
family="PH", hazard="Splines", n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino,
family="AH", hazard="Splines", n.knots=8, kappa=10^10, betaknots=1, betaorder=2)
#############################################################################
# ----- SHARED FRAILTY GENERALIZED SURVIVAL MODELS ----- #
#############################################################################
adult.retino = retinopathy[retinopathy$type == "adult", ]
adult.retino[adult.retino$futime >= 50, "status"] = 0
adult.retino[adult.retino$futime >= 50, "futime"] = 50
### --- Parametric PH, AH, PO and probit models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="PH")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="AH")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="PO")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="probit")
### --- Semi-parametric PH and AH models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt),
data=adult.retino, family="PH", hazard="Splines",
n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt),
data=adult.retino, family="AH", hazard="Splines",
n.knots=8, kappa=10^10, betaknots=1, betaorder=2)
#############################################################################
# ----- JOINT FRAILTY GENERALIZED SURVIVAL MODELS ----- #
#############################################################################
data("readmission")
readmission[, 3:5] = readmission[, 3:5]/365.25
### --- Parametric dual-PH, AH, PO and probit models --- ###
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("PH","PH"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("AH","AH"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("PO","PO"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("probit","probit"))
### --- Semi-parametric dual-PH and AH models --- ###
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE,
hazard="Splines", family=c("PH","PH"),
n.knots=5, kappa=c(100,100), betaknots=1, betaorder=3)
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE,
hazard="Splines", family=c("AH","AH"),
n.knots=5, kappa=c(600,600), betaknots=1, betaorder=3)
frailtypack documentation built on Nov. 25, 2023, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| GenfrailtyPenal | R Documentation |
Fit a Shared or a Joint Frailty Generalized Survival Model
Description
I. SHARED FRAILTY GENERALIZED SURVIVAL MODELS
Fit a gamma Shared Frailty Generalized Survival Model using a parametric estimation, or a semi-parametric penalized likelihood estimation. Right-censored data and strata (up to 6 levels) are allowed. It allows to obtain a parametric or flexible semi-parametric smooth hazard and survival functions.
Each frailty term ui is assumed
to act multiplicatively on the hazard function, and to be drawn from a
Gamma distribution with unit mean and variance \theta.
Conditional on the frailty term, the hazard function for the
jth subject in the ith group
is then expressed by
where xij
is a collection of baseline covariates,
\xi is a vector of parameters, and
\lambdaij
(t | xij ; \xi)
is the hazard function for an average value of the frailty.
The associated conditional survival function writes
where
Sij
(t | xij ; \xi)
designates the survival function for an average value of the frailty.
Following Liu et al. (2017, 2018), the latter function is expressed in terms of
a link function g(.) and a linear predictor
\etaij
(t, xij; \xi)
such that
g[Sij
(t | xij ; \xi)]
=
\etaij
(t, xij; \xi),
i.e.
Sij
(t | xij ; \xi)
=
h[\etaij
(t, xij; \xi)]
with h() = g-1().
The conditional survival function is finally modeled by
The table below summarizes the most commonly used (inverse) link functions and their associated conditional survival, hazard and cumulative hazard functions. PHM stands for "Proportional Hazards Model", POM for "Proportional Odds Model, PROM for "Probit Model" and AHM for "Additive Hazards Model".
I.(a) Fully parametric case
In the fully parametric case, linear predictors considered are of the form
where \rho > 0 is a shape parameter,
\gamma > 0 a scale parameter,
\beta a vector of regression coefficients,
and \xi = (\rho ,\gamma, \beta).
With the appropriate link function, such linear parametric predictors
make it possible to recover
a Weibull baseline survival function for PHMs and AHMs,
a log-logistic baseline survival function for POMs,
and a log-normal one for PROMs.
I. (b) Flexible semi-parametric case
For PHM and AHM, a more flexible splines-based approach is proposed for
modeling the baseline hazard function and time-varying regression coefficients.
In this case, conditional on the frailty term ui,
the hazard function for the jth subject
in the ith group is still expressed by
\lambdaij
(t | xij, ui ;
\xi)
= ui
\lambdaij
(t | xij ; \xi),
but we have this time
The smoothness of baseline hazard function \lambda0()
is ensured by penalizing the log-likelihood by a term which has
large values for rough functions.
Moreover, for parametric and flexible semi-parametric AHMs, the log-likelihood is constrained to ensure the strict positivity of the hazards, since the latter is not naturally guaranteed by the model.
II. JOINT FRAILTY GENERALIZED SURVIVAL MODELS
Fit a gamma Joint Frailty Generalized Survival Model for recurrent and terminal events using a parametric estimation, or a semi-parametric penalized likelihood estimation. Right-censored data and strata (up to 6 levels) for the recurrent event part are allowed. Joint frailty models allow studying, jointly, survival processes of recurrent and terminal events, by considering the terminal event as an informative censoring.
This model includes an common patient-specific frailty term
ui for the two survival functions which will
take into account the unmeasured heterogeneity in the data,
associated with unobserved covariates.
The frailty term acts differently for the two survival functions
(ui for the recurrent survival function and
uiα for the death one).
The covariates could be different for the recurrent and terminal event parts.
II.(a) Fully parametric case
For the jth recurrence (j=1,...,ni)
and the ith patient (i=1,...,N),
the gamma Joint Frailty Generalized Survival Model
for recurrent event survival function
SRij(.) and death survival function
SDi(.) is
- \etaRij (resp. \etaDi)
is the linear predictor for the recurrent (resp. terminal) event process.
The form of these linear predictors is the same as the one presented in I.(a).
- hR(.) (resp. hD(.))
is the inverse link function associated with
recurrent events (resp. terminal event).
- xRij and xDi
are two vectors of baseline covariates associated with
recurrent and terminal events.
- \xiR and \xiD
are the parameter vectors for recurrent and terminal events.
- \alpha is a parameter allowing more flexibility in the association
between recurrent and terminal events processes.
- The random frailties ui are still assumed iid and
drown from a \Gamma(1/\theta,1/\theta).
II.(b) Flexible semi-parametric case
If one chooses to fit a PHM or an AHM for recurrent and/or terminal events,
a splines-based approach for modeling baseline hazard functions
and time-varying regression coefficients is still available.
In this approach, the submodel for recurrent events is expressed as
\lambdaRij
(t | xRij, ui ;
\xiR)
= ui
\lambdaRij
(t | xRij ;
\xiR), where
The submodel for terminal event is expressed as
\lambdaDi
(t | xDi, ui ;
\xiD)
= uiα
\lambdaDi
(t | xDi ;
\xiD), where
Baseline hazard functions
\lambdaR0(.) and \lambdaD0(.)
are estimated using cubic M-splines (of order 4)
with positive coefficients, and the time-varying coefficients
\betaR(.) and \betaD(.)
are estimated using B-splines of order q.
The smoothness of baseline hazard functions is ensured by penalizing the log-likelihood by two terms which has large values for rough functions.
Moreover, if one chooses an AHM for recurrent and/or terminal event submodel, the log-likelihood is constrained to ensure the strict positivity of the hazards, since the latter is not naturally guaranteed by the model.
Usage
GenfrailtyPenal(formula, formula.terminalEvent, data, recurrentAG = FALSE,
family, hazard = "Splines", n.knots, kappa, betaknots = 1, betaorder = 3,
RandDist = "Gamma", init.B, init.Theta, init.Alpha, Alpha, maxit = 300,
nb.gh, nb.gl, LIMparam = 1e-3, LIMlogl = 1e-3, LIMderiv = 1e-3, print.times = TRUE,
cross.validation, jointGeneral, nb.int, initialize, init.Ksi, Ksi, init.Eta)
Arguments
formula |
A formula object, with the response on the left of a
|
formula.terminalEvent |
Only for joint frailty models: a formula object, only requires terms on the right to indicate which variables are used for the terminal event. Interactions are possible using ' * ' or ' : '. |
data |
A 'data.frame' with the variables used in ' |
recurrentAG |
Logical value. Is Andersen-Gill model fitted? If so indicates that recurrent event times with the counting process approach of Andersen and Gill is used. This formulation can be used for dealing with time-dependent covariates. The default is FALSE. |
family |
Type of Generalized Survival Model to fit.
|
hazard |
Type of hazard functions:
|
n.knots |
Integer giving the number of knots to use. Value required in
the penalized likelihood estimation. It corresponds to the |
kappa |
Positive smoothing parameter in the penalized likelihood estimation. The coefficient kappa tunes the intensity of the penalization (the integral of the squared second derivative of hazard function). In a stratified shared model, this argument must be a vector with kappas for both strata. In a stratified joint model, this argument must be a vector with kappas for both strata for recurrent events plus one kappa for terminal event. We advise the user to identify several possible tuning parameters, note their defaults and look at the sensitivity of the results to varying them. Value required. (See Note). |
betaknots |
Number of inner knots used for the
B-splines time-varying coefficient estimation. Default is 1.
See ' |
betaorder |
Order of the B-splines used for the
time-varying coefficient estimation.
Default is cubic B-splines ( |
RandDist |
Type of random effect distribution:
|
init.B |
A vector of initial values for regression coefficients. This vector should be of the same size as the whole vector of covariates with the first elements for the covariates related to the recurrent events and then to the terminal event (interactions in the end of each component). Default is 0.1 for each (for Generalized Survival and Shared Frailty Models) or 0.5 (for Generalized Joint Frailty Models). |
init.Theta |
Initial value for frailty variance. |
init.Alpha |
Only for Generalized Joint Frailty Models: initial value for parameter alpha. |
Alpha |
Only for Generalized Joint Frailty Models: input "None" so as to fit a joint model without the parameter alpha. |
maxit |
Maximum number of iterations for the Marquardt algorithm. Default is 300 |
nb.gh |
Number of nodes for the Gaussian-Hermite quadrature.
It can be chosen among 5, 7, 9, 12, 15, 20 and 32.
The default is 20 if |
nb.gl |
Number of nodes for the Gaussian-Laguerre quadrature.
It can be chosen between 20 and 32.
The default is 20 if |
LIMparam |
Convergence threshold of the Marquardt algorithm for the
parameters (see Details), |
LIMlogl |
Convergence threshold of the Marquardt algorithm for the
log-likelihood (see Details), |
LIMderiv |
Convergence threshold of the Marquardt algorithm for the
gradient (see Details), |
print.times |
A logical parameter to print iteration process. Default is TRUE. |
cross.validation |
Not implemented yet for the generalized settings. |
jointGeneral |
Not implemented yet for the generalized settings. |
nb.int |
Not implemented yet for the generalized settings. |
initialize |
Not implemented yet for the generalized settings. |
init.Ksi |
Not implemented yet for the generalized settings. |
Ksi |
Not implemented yet for the generalized settings. |
init.Eta |
Not implemented yet for the generalized settings. |
Details
TYPICAL USES
For a Generalized Survival Model:
GenfrailtyPenal( formula=Surv(time,event)~var1+var2, data, family, \dots)
For a Shared Frailty Generalized Survival Model:
GenfrailtyPenal( formula=Surv(time,event)~cluster(group)+var1+var2, data, family, \dots)
For a Joint Frailty Generalized Survival Model:
GenfrailtyPenal( formula=Surv(time,event)~cluster(group)+var1+var2+var3+terminal(death), formula.terminalEvent= ~var1+var4, data, family, \dots)
OPTIMIZATION ALGORITHM
The estimated parameters are obtained using the robust Marquardt algorithm
(Marquardt, 1963) which is a combination between a Newton-Raphson algorithm
and a steepest descent algorithm. The iterations are stopped when
the difference between two consecutive log-likelihoods is small
(<10-3),
the estimated coefficients are stable
(consecutive values <10-3,
and the gradient small enough (<10-3).
When the frailty variance is small, numerical problems may arise.
To solve this problem, an alternative formula of the penalized log-likelihood
is used (see Rondeau, 2003 for further details).
For Proportional Hazards and Additive Hazards submodels,
cubic M-splines of order 4 can be used to estimate the hazard function.
In this case, I-splines (integrated M-splines) are used to compute the
cumulative hazard function.
The inverse of the Hessian matrix is the variance estimator.
To deal with the positivity constraint of the variance component and the
spline coefficients, a squared transformation is used and the standard errors
are computed by the \Delta-method (Knight & Xekalaki, 2000).
The integrations in the full log likelihood are evaluated using
Gaussian quadrature. Laguerre polynomials with 20 points are used to treat
the integrations on [0,\infty[.
INITIAL VALUES
In case of a shared frailty model,
the splines and the regression coefficients are initialized to 0.1.
The program fits, firstly, an adjusted Cox model to give new initial values
for the splines and the regression coefficients.
The variance of the frailty term \theta is initialized to 0.1.
Then, a shared frailty model is fitted.
In case of a joint frailty model,
the splines and the regression coefficients are initialized to 0.5.
The program fits firstly, an adjusted Cox model to have new initial values
for the splines and the regression coefficients.
The variance of the frailty term \theta and the association parameter
\alpha are initialized to 1.
Then, a joint frailty model is fitted.
Value
The following components are included in a 'frailtyPenal' object for each model.
b |
Sequence of the corresponding estimation of the coefficients for the hazard functions (parametric or semiparametric), the random effects variances and the regression coefficients. |
call |
The code used for the model. |
formula |
The formula part of the code used for the model. |
n |
The number of observations used in the fit. |
groups |
The maximum number of groups used in the fit. |
n.events |
The number of events observed in the fit. |
n.eventsbygrp |
A vector of length the number of groups giving the number of observed events in each group. |
loglik |
The marginal log-likelihood in the parametric case. |
loglikPenal |
The marginal penalized log-likelihood in the semiparametric case. |
coef |
The regression coefficients. |
varH |
The variance matrix of the regression coefficients before positivity constraint transformation. Then, the delta method is needed to obtain the estimated variance parameters. That is why some variances don't match with the printed values at the end of the model. |
varHtotal |
The variance matrix of all the parameters before positivity constraint transformation. Then, the delta method is needed to obtain the estimated variance parameters. That is why some variances don't match with the printed values at the end of the model. |
varHIH |
The robust estimation of the variance matrix of the regression coefficients |
varHIHtotal |
The robust estimation of the variance matrix of all parameters. |
x |
Matrix of times where the hazard functions are estimated. |
xSu |
Matrix of times where the survival functions are estimated. |
lam |
Array (dim=3) of baseline hazard estimates and confidence bands. |
surv |
Array (dim=3) of baseline survival estimates and confidence bands. |
type |
Character string specifying the type of censoring,
see the |
n.strat |
Number of strata. |
n.iter |
Number of iterations needed to converge. |
median |
The value of the median survival and its confidence bands. If there are two strata or more, the first value corresponds to the value for the first strata, etc. |
LCV |
The approximated likelihood cross-validation criterion in the
semiparametric case.
With H (resp. Hpen) the hessian matrix
of log-likelihood (resp. penalized log-likelihood),
EDF = Hpen-1 H
the effective degrees of freedom,
L(
|
AIC |
The Akaike information Criterion for the parametric case.
With p the number of parameters,
n the number of observations and L(
|
npar |
Number of parameters. |
nvar |
Number of explanatory variables. |
typeof |
Indicator of the type of hazard functions computed : 0 for "Splines", 2 for "parametric". |
istop |
Convergence indicator: 1 if convergence is reached, 2 if convergence is not reached, 3 if the hessian matrix is not positive definite, 4 if a numerical problem has occurred in the likelihood calculation |
shape.param |
Shape parameter for the parametric hazard function (a Weibull distribution is used for proportional and additive hazards models, a log-logistic distribution is used for proportional odds models, a log-normal distribution is used for probit models). |
scale.param |
Scale parameter for the parametric hazard function. |
Names.data |
Name of the dataset. |
Frailty |
Logical value. Was model with frailties fitted ? |
linear.pred |
Linear predictor:
|
BetaTpsMat |
Matrix of time varying-effects and confidence bands (the first column used for abscissa of times). |
nvartimedep |
Number of covariates with time-varying effects. |
Names.vardep |
Name of the covariates with time-varying effects. |
EPS |
Convergence criteria concerning the parameters, the likelihood and the gradient. |
family |
Type of Generalized Survival Model fitted (0 for PH, 1 for PO, 2 for probit, 3 for AH). |
global_chisq.test |
A binary variable equals to 0 when no multivariate Wald is given, 1 otherwise. |
beta_p.value |
p-values of the Wald test for the estimated regression coefficients. |
cross.Val |
Logical value. Is cross validation procedure used for estimating the smoothing parameters in the penalized likelihood estimation? |
DoF |
Degrees of freedom associated with the smoothing parameter
|
kappa |
A vector with the smoothing parameters in the penalized likelihood estimation corresponding to each baseline function as components. |
n.knots |
Number of knots for estimating the baseline functions in the penalized likelihood estimation. |
n.knots.temp |
Initial value for the number of knots. |
global_chisq |
A vector with the values of each multivariate Wald test. |
dof_chisq |
A vector with the degree of freedom for each multivariate Wald test. |
p.global_chisq |
A vector with the p-values for each global multivariate Wald test. |
names.factor |
Names of the "as.factor" variables. |
Xlevels |
Vector of the values that factor might have taken. |
The following components are specific to shared models.
equidistant |
Indicator for the intervals used in the spline estimation
of baseline hazard functions :
1 for equidistant intervals ; 0 for intervals using percentile
(note: |
Names.cluster |
Cluster names. |
theta |
Variance of the gamma frailty parameter, i.e.
Var( |
varTheta |
Variance of parameter |
theta_p.value |
p-value of the Wald test for the estimated variance of the gamma frailty. |
The following components are specific to joint models.
formula |
The formula part of the code used for the recurrent events. |
formula.terminalEvent |
The formula part of the code used for the terminal model. |
n.deaths |
Number of observed deaths. |
n.censored |
Number of censored individuals. |
theta |
Variance of the gamma frailty parameter, i.e.
Var( |
indic_alpha |
Indicator if a joint frailty model with
|
alpha |
The coefficient |
nvar |
A vector with the number of covariates of each type of hazard function as components. |
nvarnotdep |
A vector with the number of constant effect covariates of each type of hazard function as components. |
nvarRec |
Number of recurrent explanatory variables. |
nvarEnd |
Number of death explanatory variables. |
noVar1 |
Indicator of recurrent explanatory variables. |
noVar2 |
Indicator of death explanatory variables. |
Names.vardep |
Name of the covariates with time-varying effects for the recurrent events. |
Names.vardepdc |
Name of the covariates with time-varying effects for the terminal event. |
xR |
Matrix of times where both survival and hazard function are estimated for the recurrent event. |
xD |
Matrix of times for the terminal event. |
lamR |
Array (dim=3) of hazard estimates and confidence bands for recurrent event. |
lamD |
The same value as |
survR |
Array (dim=3) of baseline survival estimates and confidence bands for recurrent event. |
survD |
The same value as |
nb.gh |
Number of nodes for the Gaussian-Hermite quadrature. |
nb.gl |
Number of nodes for the Gaussian-Laguerre quadrature. |
medianR |
The value of the median survival for the recurrent events and its confidence bands. |
medianD |
The value of the median survival for the terminal event and its confidence bands. |
names.factor |
Names of the "as.factor" variables for the recurrent events. |
names.factordc |
Names of the "as.factor" variables for the terminal event. |
Xlevels |
Vector of the values that factor might have taken for the recurrent events. |
Xlevels2 |
Vector of the values that factor might have taken for the terminal event. |
linear.pred |
Linear predictor for the recurrent part:
|
lineardeath.pred |
Linear predictor for the terminal part:
|
Xlevels |
Vector of the values that factor might have taken for the recurrent part. |
Xlevels2 |
vector of the values that factor might have taken for the death part. |
BetaTpsMat |
Matrix of time varying-effects and confidence bands for recurrent event (the first column used for abscissa of times of recurrence). |
BetaTpsMatDc |
Matrix of time varying-effects and confidence bands for terminal event (the first column used for abscissa of times of death). |
alpha_p.value |
p-value of the Wald test for the estimated |
Note
In the flexible semiparametric case, smoothing parameters kappa and
number of knots n.knots are the arguments that the user have to change
if the fitted model does not converge.
n.knots takes integer values between 4 and 20.
But with n.knots=20, the model would take a long time to converge.
So, usually, begin first with n.knots=7, and increase it step by step
until it converges.
kappa only takes positive values. So, choose a value for kappa (for
instance 10000), and if it does not converge, multiply or divide this value
by 10 or 5 until it converges.
References
J. Chauvet and V. Rondeau (2021). A flexible class of generalized joint frailty models for the analysis of survival endpoints. In revision.
Liu XR, Pawitan Y, Clements M. (2018) Parametric and penalized generalized survival models. Statistical Methods in Medical Research 27(5), 1531-1546.
Liu XR, Pawitan Y, Clements MS. (2017) Generalized survival models for correlated time-to-event data. Statistics in Medicine 36(29), 4743-4762.
A. Krol, A. Mauguen, Y. Mazroui, A. Laurent, S. Michiels and V. Rondeau (2017). Tutorial in Joint Modeling and Prediction: A Statistical Software for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event. Journal of Statistical Software 81(3), 1-52.
V. Rondeau, Y. Mazroui and J. R. Gonzalez (2012). Frailtypack: An R package for the analysis of correlated survival data with frailty models using penalized likelihood estimation or parametric estimation. Journal of Statistical Software 47, 1-28.
V. Rondeau, J.P. Pignon, S. Michiels (2011). A joint model for the dependance between clustered times to tumour progression and deaths: A meta-analysis of chemotherapy in head and neck cancer. Statistical methods in medical research 897, 1-19.
V. Rondeau, S. Mathoulin-Pellissier, H. Jacqmin-Gadda, V. Brouste, P. Soubeyran (2007). Joint frailty models for recurring events and death using maximum penalized likelihood estimation:application on cancer events. Biostatistics 8,4, 708-721.
V. Rondeau, D. Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gamma-frailty model. Lifetime Data Analysis 9, 139-153.
C.A. McGilchrist, and C.W. Aisbett (1991). Regression with frailty in survival analysis. Biometrics 47, 461-466.
D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431-441.
See Also
Surv,
terminal,
timedep
Examples
#############################################################################
# ----- GENERALIZED SURVIVAL MODELS (without frailties) ----- #
#############################################################################
adult.retino = retinopathy[retinopathy$type == "adult", ]
adult.retino[adult.retino$futime >= 50, "status"] = 0
adult.retino[adult.retino$futime >= 50, "futime"] = 50
### --- Parametric PH, AH, PO and probit models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="PH")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="AH")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="PO")
GenfrailtyPenal(formula=Surv(futime,status)~trt, data=adult.retino,
hazard="parametric", family="probit")
### --- Semi-parametric PH and AH models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino,
family="PH", hazard="Splines", n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
GenfrailtyPenal(formula=Surv(futime,status)~timedep(trt), data=adult.retino,
family="AH", hazard="Splines", n.knots=8, kappa=10^10, betaknots=1, betaorder=2)
#############################################################################
# ----- SHARED FRAILTY GENERALIZED SURVIVAL MODELS ----- #
#############################################################################
adult.retino = retinopathy[retinopathy$type == "adult", ]
adult.retino[adult.retino$futime >= 50, "status"] = 0
adult.retino[adult.retino$futime >= 50, "futime"] = 50
### --- Parametric PH, AH, PO and probit models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="PH")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="AH")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="PO")
GenfrailtyPenal(formula=Surv(futime,status)~trt+cluster(id), data=adult.retino,
hazard="parametric", family="probit")
### --- Semi-parametric PH and AH models --- ###
GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt),
data=adult.retino, family="PH", hazard="Splines",
n.knots=8, kappa=10^6, betaknots=1, betaorder=2)
GenfrailtyPenal(formula=Surv(futime,status)~cluster(id)+timedep(trt),
data=adult.retino, family="AH", hazard="Splines",
n.knots=8, kappa=10^10, betaknots=1, betaorder=2)
#############################################################################
# ----- JOINT FRAILTY GENERALIZED SURVIVAL MODELS ----- #
#############################################################################
data("readmission")
readmission[, 3:5] = readmission[, 3:5]/365.25
### --- Parametric dual-PH, AH, PO and probit models --- ###
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("PH","PH"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("AH","AH"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("PO","PO"))
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+chemo,
formula.terminalEvent=~sex+dukes+chemo, data=readmission, recurrentAG=TRUE,
hazard="parametric", family=c("probit","probit"))
### --- Semi-parametric dual-PH and AH models --- ###
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE,
hazard="Splines", family=c("PH","PH"),
n.knots=5, kappa=c(100,100), betaknots=1, betaorder=3)
GenfrailtyPenal(
formula=Surv(t.start,t.stop,event)~cluster(id)+terminal(death)+sex+dukes+timedep(chemo),
formula.terminalEvent=~sex+dukes+timedep(chemo), data=readmission, recurrentAG=TRUE,
hazard="Splines", family=c("AH","AH"),
n.knots=5, kappa=c(600,600), betaknots=1, betaorder=3)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.