| Title: | Mixed Effect Excess Hazard Models |
|---|---|
| Description: | Fit flexible (excess) hazard regression models with the possibility of including non-proportional effects of covariables and of adding a random effect at the cluster level (corresponding to a shared frailty). A detailed description of the package functionalities is provided in Charvat and Belot (2021) <doi: 10.18637/jss.v098.i14>. |
| Authors: | Hadrien Charvat, Aurelien Belot, Juste Goungounga |
| Maintainer: | Hadrien Charvat <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 2.6 |
| Built: | 2024-09-27 06:17:13 UTC |
| Source: | CRAN |
Fit an (excess) hazard regression model using different shapes for the
baseline hazard (Weibull, piecewise constant, exponential of a B-spline
of degree 1 to 3, exponential of a restricted cubic spline), with the
possibility to include time-dependent effects of variable(s) and a
random effect defined at the cluster level. Follow-up time can be
entered in the right-censored or counting process input style. The
latter allows the modelling of survival data with delayed entries. The
time-dependent effect of a covariable is modelled by adding interaction
terms between the covariable and a function of time of the same class as
the one used for the baseline hazard (in particular, with the same knots
for piecewise constant hazards; and with the same degree and the same
knots for B-spline or restricted cubic spline functions). The random
effect is assumed to be normally distributed with mean 0 and standard
deviation sigma. The optimisation process uses adaptive Gaussian
quadrature to calculate the cluster-specific marginal likelihoods. The
logarithm of the full marginal likelihood, defined as the sum of the
logarithms of the cluster-specific marginal likelihoods, is then
maximised using an optimisation routine such as nlm (default) or
optim. Functions to compute and plot the predicted (excess)
hazard and (net) survival are provided. In the case of a random
intercept model, three types of predictions can be obtained : marginal
(population-averaged), cluster-specific posterior (cluster-averaged), or
conditional (for a given quantile of the ditribution of the random effect).
Hadrien Charvat, Aurelien Belot
Charvat H, Remontet L, Bossard N, Roche L, Dejardin O, Rachet B, Launoy G, Belot A; CENSUR Working Survival Group. A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates. Stat Med 2016;35(18):3066-3084 (doi: 10.1002/sim.6881)
Charvat H, Belot A. An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. J Stat Softw 2021;98(14):1-36 (doi: 10.18637/jss.v098.i14)
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") ## Examples of syntax for various models (not run) ## Fit of a fixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution and with effects of age (agecr), ## deprivation index (depindex) and sex (IsexH) using the optim ## procedure and the BFGS method (see help of optim). # Mod_weib <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="weibull", # expected="popmrate", verbose=1000, fnoptim="optim", # method="BFGS") ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2mix <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", # degree=3, knots=c(1,5), expected="popmrate", random="clust", # verbose=1000) ## Fit of a fixed-effect overall hazard model, with the baseline hazard ## described by a piecewise constant function with the following vector ## of knots (defining the endpoints of the intervals on which the hazard ## is constant): (1,3,5,8), and with effects of age (agecr), deprivation ## index (depindex) and sex (IsexH) # Mod_pw <- mexhaz(formula=Surv(time=timesurv, event=vstat)~ # agecr+depindex+IsexH, data= simdatn1, base="pw.cst", knots=c(1,3,5,8), # verbose=1000) ## Fit of a fixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2 <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", # degree=3, knots=c(1,5), expected="popmrate", verbose=1000) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2mix <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", # degree=3, knots=c(1,5), expected="popmrate", random="clust", # verbose=1000) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year, with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) ## and with a time-dependent effect for age (agecr) and sex (IsexH). # Mod_bs3_2mixnph <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH + nph(agecr+IsexH), data=simdatn1, # base="exp.bs", degree=3, knots=c(1,5), expected="popmrate", # random="clust", verbose=1000)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") ## Examples of syntax for various models (not run) ## Fit of a fixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution and with effects of age (agecr), ## deprivation index (depindex) and sex (IsexH) using the optim ## procedure and the BFGS method (see help of optim). # Mod_weib <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="weibull", # expected="popmrate", verbose=1000, fnoptim="optim", # method="BFGS") ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2mix <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", # degree=3, knots=c(1,5), expected="popmrate", random="clust", # verbose=1000) ## Fit of a fixed-effect overall hazard model, with the baseline hazard ## described by a piecewise constant function with the following vector ## of knots (defining the endpoints of the intervals on which the hazard ## is constant): (1,3,5,8), and with effects of age (agecr), deprivation ## index (depindex) and sex (IsexH) # Mod_pw <- mexhaz(formula=Surv(time=timesurv, event=vstat)~ # agecr+depindex+IsexH, data= simdatn1, base="pw.cst", knots=c(1,3,5,8), # verbose=1000) ## Fit of a fixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2 <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", # degree=3, knots=c(1,5), expected="popmrate", verbose=1000) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2mix <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", # degree=3, knots=c(1,5), expected="popmrate", random="clust", # verbose=1000) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year, with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) ## and with a time-dependent effect for age (agecr) and sex (IsexH). # Mod_bs3_2mixnph <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH + nph(agecr+IsexH), data=simdatn1, # base="exp.bs", degree=3, knots=c(1,5), expected="popmrate", # random="clust", verbose=1000)
Function for computing direct adjusted survival estimates
from a model fitted with the mexhaz. It can be used to obtain
direct adjusted survival estimates for one or two populations. In the
latter case, survival difference estimates are also
computed. Corresponding variance estimates are based on the Delta Method
(based on the assumption of multivariate normality of the model
parameter estimates). When the model includes a random effect, three types of
predictions can be made: (i) marginal predictions (obtained by
integration over the random effect distribution), (ii) cluster-specific
posterior predictions for an existing cluster, or (iii) conditional
predictions for a given quantile of the random effect distribution (by
default, for the median value, that is, 0).
adjsurv(object, time.pts, data, data.0 = NULL, weights = NULL, marginal = TRUE, quant.rdm = 0.5, cluster = NULL, quant.rdm.0 = 0.5, cluster.0 = NULL, level = 0.95, dataset = NULL)adjsurv(object, time.pts, data, data.0 = NULL, weights = NULL, marginal = TRUE, quant.rdm = 0.5, cluster = NULL, quant.rdm.0 = 0.5, cluster.0 = NULL, level = 0.95, dataset = NULL)
object |
an object of class |
time.pts |
a vector of numerical values representing the time points at which predictions are requested. Time values greater than the maximum follow-up time on which the model estimation was based are discarded. |
data |
a |
data.0 |
an optional |
weights |
optional argument specifying the weights to be
associated with each row of |
marginal |
logical value controlling the type of predictions
returned by the function when the model includes a random
intercept. When |
quant.rdm |
numerical value (between 0 and 1) specifying the
quantile of the random effect distribution that should be used when
requesting conditional predictions. The default value is set to 0.5
(corresponding to the median, that is a value of the random effect
of 0). This argument is ignored if the model is a fixed effect
model, when the |
cluster |
a single value corresponding to the name of the cluster for
which posterior predictions should be calculated. These predictions
are obtained by integrating over the cluster-specific posterior
distribution of the random effect and thus require the original
dataset. The dataset can either be provided as part of the
|
quant.rdm.0 |
random effect distribution quantile value to be used with |
cluster.0 |
cluster value to be used with |
level |
a number in (0,1) specifying the level of confidence for
computing the confidence intervals of the hazard and the
survival. By default, |
dataset |
original dataset used to fit the |
An object of class resMexhaz that can be used by the function
plot.resMexhaz to produce graphics
of the direct adjusted survival curve. It contains the following
elements:
results |
a |
type |
type of results returned by the function. The value is
used by |
multiobs |
value used by
|
ci.method |
method used to compute confidence limits. Currently set
to |
level |
level of confidence used to compute confidence limits. |
Hadrien Charvat, Aurelien Belot
Charvat H, Remontet L, Bossard N, Roche L, Dejardin O, Rachet B, Launoy G, Belot A; CENSUR Working Survival Group. A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates. Stat Med 2016;35:3066-3084 (doi: 10.1002/sim.6881)
Skrondal A, Rabe-Hesketh S. Prediction in multilevel generalized linear models. J R Stat Soc A Stat Soc 2009;172(3):659-687 (doi: 10.1111/j.1467-985X.2009.00587.x).
plot.resMexhaz, lines.resMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Direct adjusted survival for the simdatn1 population DAS_Modbs1_2 <- adjsurv(Mod_bs1_2, time.pts=seq(1,10), data=simdatn1)data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Direct adjusted survival for the simdatn1 population DAS_Modbs1_2 <- adjsurv(Mod_bs1_2, time.pts=seq(1,10), data=simdatn1)
This is a generic function.
fixef(x, ...)fixef(x, ...)
x |
a fitted object from which fixed effects can be extracted. |
... |
may be required by some methods using this generic function. |
see the documentation for the specific class.
## See the specific class documentation.## See the specific class documentation.
Display a vector containing the fixed effects of a mexhaz model.
## S3 method for class 'mexhaz' fixef(x, ...)## S3 method for class 'mexhaz' fixef(x, ...)
x |
an object of class |
... |
not used. |
A vector containing the fixed effect estimates.
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") fixef(Mod_weib_mix)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") fixef(Mod_weib_mix)
Function for adding to an already existing graphical window
the predicted (excess) hazard or (net) survival based on a predMexhaz object.
## S3 method for class 'predMexhaz' lines(x, which = c("surv", "hazard"), conf.int = TRUE, lty.pe = "solid", lty.ci = "dashed", ...)## S3 method for class 'predMexhaz' lines(x, which = c("surv", "hazard"), conf.int = TRUE, lty.pe = "solid", lty.ci = "dashed", ...)
x |
an object of class |
which |
type of curve to be plotted. Selection can be made
between |
conf.int |
logical values allowing the user to decide whether to plot the confidence limits of the survival (or hazard). |
lty.pe |
type of line used for drawing the hazard/survival estimate. |
lty.ci |
type of line used for drawing the confidence limits. |
... |
additional parameters that are directly passed to the
|
predict.mexhaz, plot.predMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs2_2A <- predict(Mod_bs2_2, time.pts=seq(0,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) plot(Pred_Modbs2_2A, which="hazard", col="red") lines(Pred_Modbs2_2A, which="hazard", conf.int=FALSE)data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs2_2A <- predict(Mod_bs2_2, time.pts=seq(0,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) plot(Pred_Modbs2_2A, which="hazard", col="red") lines(Pred_Modbs2_2A, which="hazard", conf.int=FALSE)
Function for adding to an already existing graphical window
various quantities calculated from a mexhaz model.
## S3 method for class 'resMexhaz' lines(x, conf.int = TRUE, lty.pe = "solid", lty.ci = "blank", col.ci = "blue", alpha.col.ci = 0.25, ...)## S3 method for class 'resMexhaz' lines(x, conf.int = TRUE, lty.pe = "solid", lty.ci = "blank", col.ci = "blue", alpha.col.ci = 0.25, ...)
x |
an object of class |
conf.int |
logical value allowing the user to decide whether to plot the confidence limits of the survival (or hazard). |
lty.pe |
type of line used for drawing the hazard/survival estimate. |
lty.ci |
type of line used for drawing the confidence limits. |
col.ci |
color used to fill in the polygon defined by the confidence limits. |
alpha.col.ci |
parameter used internally by the |
... |
additional parameters that are directly passed to the
|
adjsurv, riskfunc,
plot.resMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Relative risk (ratio of cumulative incidence curves) for men versus ## women for two different values of age RR_Modbs2a <- riskfunc(Mod_bs2, time.pts=seq(0,10, by=0.1), data=data.frame(agecr=0.1, IsexH=1, depindex=0), data.0=data.frame(agecr=0.1, IsexH=0, depindex=0), conf.int="delta", type="rr") RR_Modbs2b <- riskfunc(Mod_bs2, time.pts=seq(0,10, by=0.1), data=data.frame(agecr=-0.1, IsexH=1, depindex=0), data.0=data.frame(agecr=-0.1, IsexH=0, depindex=0), conf.int="delta", type="rr") plot(RR_Modbs2a) lines(RR_Modbs2a, col.ci="green") abline(h=1, lty="dashed")data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Relative risk (ratio of cumulative incidence curves) for men versus ## women for two different values of age RR_Modbs2a <- riskfunc(Mod_bs2, time.pts=seq(0,10, by=0.1), data=data.frame(agecr=0.1, IsexH=1, depindex=0), data.0=data.frame(agecr=0.1, IsexH=0, depindex=0), conf.int="delta", type="rr") RR_Modbs2b <- riskfunc(Mod_bs2, time.pts=seq(0,10, by=0.1), data=data.frame(agecr=-0.1, IsexH=1, depindex=0), data.0=data.frame(agecr=-0.1, IsexH=0, depindex=0), conf.int="delta", type="rr") plot(RR_Modbs2a) lines(RR_Modbs2a, col.ci="green") abline(h=1, lty="dashed")
Fit an (excess) hazard regression model using different shapes for the
baseline hazard (Weibull, piecewise constant, exponential of a B-spline
of degree 1 to 3, exponential of a restricted cubic spline), with the
possibility to include time-dependent effects of variable(s) and a
random effect defined at the cluster level. The function accepts
right-censored and counting process input styles for the follow-up
time. The latter allows the modelling of survival data with delayed
entries. The time-dependent effect of a covariable is modelled by adding
interaction terms between the covariable and a function of time of the
same class as the one used for the baseline hazard (in particular, with
the same knots for piecewise constant hazards; and with the same degree
and the same knots for B-spline or restricted cubic spline
functions). The random effect is assumed to be normally distributed with
mean 0 and standard deviation sigma. The optimisation process uses
adaptive Gaussian quadrature to calculate the cluster-specific marginal
likelihoods. The logarithm of the full marginal likelihood, defined as
the sum of the logarithms of the cluster-specific marginal likelihoods,
is then maximised using an optimisation routine such as nlm
(default) or optim.
mexhaz(formula, data, expected = NULL, base = c("weibull", "exp.bs", "exp.ns", "pw.cst"), degree = 3, knots = NULL, bound = NULL, n.gleg = 20, init = NULL, random = NULL, n.aghq = 10, recurrent = FALSE, fnoptim = c("nlm", "optim"), verbose = 0, method = "Nelder-Mead", iterlim = 10000, numHess = FALSE, print.level = 1, exactGradHess = TRUE, gradtol = ifelse(exactGradHess, 1e-8, 1e-6), testInit = TRUE, keep.data = FALSE, ...)mexhaz(formula, data, expected = NULL, base = c("weibull", "exp.bs", "exp.ns", "pw.cst"), degree = 3, knots = NULL, bound = NULL, n.gleg = 20, init = NULL, random = NULL, n.aghq = 10, recurrent = FALSE, fnoptim = c("nlm", "optim"), verbose = 0, method = "Nelder-Mead", iterlim = 10000, numHess = FALSE, print.level = 1, exactGradHess = TRUE, gradtol = ifelse(exactGradHess, 1e-8, 1e-6), testInit = TRUE, keep.data = FALSE, ...)
formula |
a formula object, with the response on the left of the In case |
data |
a |
expected |
name of the variable (must be given in quotes) representing the
population (i.e., expected) hazard. By default, |
base |
functional form that should be used to model the baseline
hazard. Selection can be made between the following options:
For the Weibull hazard model, the cumulative hazard is given by the following relationship: H(t,x,z) = lambda*exp(x'b)*t^(rho*exp(z'g)) where lambda and rho are the parameters of the Weibull baseline hazard,
x represent variables with proportional effect (with corresponding
regression coefficients 'b') and z represent variables with
time-dependent effects (with corresponding regression coefficients
'g'). The For the spline-based hazards, it should be noted that the B-spline
and restricted cubic spline bases created internally in the
|
degree |
if |
knots |
if |
bound |
a vector of two numerical values corresponding to the boundary knots
of the spline functions. If |
n.gleg |
if |
init |
vector of initial values. By default for the baseline hazard: if if - if - if - if - if the parameters describing the effects of the covariables are all set to 0; the parameter representing the standard deviation of the random
effect is set to 0.1. In case of failed convergence, if
|
random |
name of the variable to be entered as a random effect (must be given
between quotes), representing the cluster membership. By default,
|
n.aghq |
number of quadrature points used for estimating the
cluster-specific marginal likelihoods by adaptive Gauss-Hermite
quadrature. By default, |
recurrent |
logical value indicating that the dataset
corresponds to recurrent events expressed in calendar time, in which
case the |
fnoptim |
name of the R optimisation procedure used to maximise the
likelihood. Selection can be made between |
verbose |
integer parameter representing the frequency at which the current state
of the optimisation process is displayed. Internally, an 'evaluation' is
defined as an estimation of the log-likelihood for a given vector of
parameters. This means that the number of evaluations is increased each
time the optimisation procedure updates the value of any of the
parameters to be estimated. If |
method |
if |
iterlim |
if |
numHess |
logical value allowing the user to choose between the Hessian returned
by the optimization algorithm (default) or the Hessian estimated by
the |
print.level |
this argument is only used if |
exactGradHess |
logical value allowing the user to decide
whether maximisation of the likelihood should be based on the
analytic gradient and Hessian computed internally (default,
corresponding to |
gradtol |
this argument is only used if |
testInit |
this argument is used only when
|
keep.data |
logical argument determining whether the dataset
should be kept in the object returned by the function: this can be
useful in certain contexts (e.g., to calculate cluster-specific
posterior predictions from a random intercept model) but might create
unnecessarily voluminous objects. The default
value is set to |
... |
represents additional parameters directly passed to |
An object of class mexhaz containing the following elements:
dataset |
name of the dataset used to fit the model. |
call |
function call on which the model is based. |
formula |
formula part of the call. |
expected |
name of the variable corresponding to the expected
hazard (takes the value |
xlevels |
information concerning the levels of the categorical
variables used in the model (used by |
n.obs.tot |
total number of observations in the dataset. |
n.obs |
number of observations used to fit the model (after exclusion of missing values). |
n.events |
number of events (after exclusion of missing values). |
n.clust |
number of clusters. |
n.time.0 |
number of observations for which the observed follow-up time was equal to 0 (only for right censored type data). |
base |
function used to model the baseline hazard. |
max.time |
maximal observed time in the dataset. |
boundary.knots |
vector of boundary values used to define the B-spline (or natural spline) bases. |
degree |
degree of the B-spline used to model the logarithm of the baseline hazard. |
knots |
vector of interior knots used to define the B-spline (or natural spline) bases. |
names.ph |
names of the covariables with a proportional effect. |
random |
name of the variable defining cluster membership (set to
|
recurrent |
logical value indicating whether the dataset corresponds to recurrent events. |
init |
a vector containing the initial values of the parameters. |
coefficients |
a vector containing the parameter estimates. |
std.errors |
a vector containing the standard errors of the parameter estimates. |
vcov |
the variance-covariance matrix of the estimated parameters. |
gradient |
the gradient of the log-likelihood function evaluated at the estimated parameters. |
hessian |
the Hessian of the log-likelihood function evaluated at the estimated parameters. |
mu.hat |
a |
var.mu.hat |
the covariance matrix of the cluster-specific shrinkage estimators. |
vcov.fix.mu.hat |
a |
data |
original dataset used to fit the model (if
|
n.par |
number of estimated parameters. |
n.gleg |
number of Gauss-Legendre quadrature points used to calculate the cumulative (excess) hazard (only relevant if a B-spline of degree 2 or 3 or a cubic restricted spline was used to model the logarithm of the baseline hazard). |
n.aghq |
number of adaptive Gauss-Hermite quadrature points used to calculate the cluster-specific marginal likelihoods (only relevant if a multi-level model is fitted). |
fnoptim |
name of the R optimisation procedure used to maximise the likelihood. |
method |
optimisation method used by |
code |
code (integer) indicating the status of the optimisation
process (this code has a different meaning for |
loglik |
value of the log-likelihood at the end of the optimisation procedure. |
iter |
number of iterations used in the optimisation process. |
eval |
number of evaluations used in the optimisation process. |
time.elapsed |
total time required to reach convergence. |
Hadrien Charvat, Aurelien Belot
Charvat H, Remontet L, Bossard N, Roche L, Dejardin O, Rachet B, Launoy G, Belot A; CENSUR Working Survival Group. A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates. Stat Med 2016;35(18):3066-3084 (doi: 10.1002/sim.6881)
Charvat H, Belot A. An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. J Stat Softw 2021;98(14):1-36 (doi: 10.18637/jss.v098.i14)
fixef.mexhaz, predict.mexhaz, print.mexhaz,
ranef.mexhaz, summary.mexhaz, update.mexhaz, vcov.mexhaz
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") ## More complex examples (not run) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2mix_nph <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH+nph(agecr), data=simdatn1, # base="exp.bs", degree=3, knots=c(1,5), expected="popmrate", # random="clust", verbose=1000) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a restricted cubic spline with two knots at the ## tertiles of the distribution of follow-up times for individuals who ## presented the event and with effects of age (agecr) and sex (IsexH) # Mod_ns3_2mix_nph <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+nph(agecr), data=simdatn1, base="exp.ns", degree=3, # knots=quantile(simdatn1[simdatn1$vstat==1,]$timesurv, probs=c(1:2/3)), # expected="popmrate", random="clust", verbose=1000)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") ## More complex examples (not run) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) # Mod_bs3_2mix_nph <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+depindex+IsexH+nph(agecr), data=simdatn1, # base="exp.bs", degree=3, knots=c(1,5), expected="popmrate", # random="clust", verbose=1000) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a restricted cubic spline with two knots at the ## tertiles of the distribution of follow-up times for individuals who ## presented the event and with effects of age (agecr) and sex (IsexH) # Mod_ns3_2mix_nph <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+nph(agecr), data=simdatn1, base="exp.ns", degree=3, # knots=quantile(simdatn1[simdatn1$vstat==1,]$timesurv, probs=c(1:2/3)), # expected="popmrate", random="clust", verbose=1000)
Function for plotting the predicted (excess) hazard or
(net) survival based on a predMexhaz object.
## S3 method for class 'predMexhaz' plot(x, which = c("surv", "hazard"), conf.int = TRUE, lty.pe = "solid", lty.ci = "dashed", ...)## S3 method for class 'predMexhaz' plot(x, which = c("surv", "hazard"), conf.int = TRUE, lty.pe = "solid", lty.ci = "dashed", ...)
x |
an object of class |
which |
type of curve to be plotted. Selection can be made
between |
conf.int |
logical value allowing the user to decide whether to plot the confidence limits of the survival (or hazard). |
lty.pe |
type of line used for drawing the hazard/survival estimate. |
lty.ci |
type of line used for drawing the confidence limits. |
... |
additional parameters that are directly passed to the
|
predict.mexhaz, points.predMexhaz, lines.predMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs2_2A <- predict(Mod_bs2_2, time.pts=seq(0,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) plot(Pred_Modbs2_2A, which="hazard")data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs2_2A <- predict(Mod_bs2_2, time.pts=seq(0,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) plot(Pred_Modbs2_2A, which="hazard")
Function for plotting various quantities calculated from a
mexhaz model.
## S3 method for class 'resMexhaz' plot(x, conf.int = TRUE, lty.pe = "solid", lty.ci = "blank", col.ci = "blue", alpha.col.ci = 0.25, ylim = NULL, ...)## S3 method for class 'resMexhaz' plot(x, conf.int = TRUE, lty.pe = "solid", lty.ci = "blank", col.ci = "blue", alpha.col.ci = 0.25, ylim = NULL, ...)
x |
an object of class |
conf.int |
logical value allowing the user to decide whether to plot the confidence limits of the survival (or hazard). |
lty.pe |
type of line used for drawing the hazard/survival estimate. |
lty.ci |
type of line used for drawing the confidence limits. |
col.ci |
color used to fill in the polygon defined by the confidence limits. |
alpha.col.ci |
parameter used internally by the |
ylim |
when set to |
... |
additional parameters that are directly passed to the
|
adjsurv, riskfunc,
lines.resMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Relative risk (ratio of cumulative incidence curves) for men versus ## women RR_Modbs2 <- riskfunc(Mod_bs2, time.pts=seq(0,10, by=0.1), data=data.frame(agecr=0.1, IsexH=1, depindex=0), data.0=data.frame(agecr=0.1, IsexH=0, depindex=0), conf.int="delta", type="rr") plot(RR_Modbs2) abline(h=1, lty="dashed")data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Relative risk (ratio of cumulative incidence curves) for men versus ## women RR_Modbs2 <- riskfunc(Mod_bs2, time.pts=seq(0,10, by=0.1), data=data.frame(agecr=0.1, IsexH=1, depindex=0), data.0=data.frame(agecr=0.1, IsexH=0, depindex=0), conf.int="delta", type="rr") plot(RR_Modbs2) abline(h=1, lty="dashed")
Function for adding to an already existing graphical window
the predicted (excess) hazard or (net) survival based on a predMexhaz object.
## S3 method for class 'predMexhaz' points(x, which = c("surv", "hazard"), conf.int = TRUE, lty.pe = "solid", lty.ci = "dashed", ...)## S3 method for class 'predMexhaz' points(x, which = c("surv", "hazard"), conf.int = TRUE, lty.pe = "solid", lty.ci = "dashed", ...)
x |
an object of class |
which |
type of curve to be plotted. Selection can be made
between |
conf.int |
logical values allowing the user to decide whether to plot the confidence limits of the survival (or hazard). |
lty.pe |
type of line used for drawing the hazard/survival estimate (used when 'type="l"'). |
lty.ci |
type of line used for drawing the confidence limits (used when 'type="l"'). |
... |
additional parameters that are directly passed to the
|
predict.mexhaz, plot.predMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs2_2A <- predict(Mod_bs2_2, time.pts=seq(0,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) plot(Pred_Modbs2_2A, which="hazard", col="red") points(Pred_Modbs2_2A, which="hazard", type="l", conf.int=FALSE)data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs2_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=2, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs2_2A <- predict(Mod_bs2_2, time.pts=seq(0,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) plot(Pred_Modbs2_2A, which="hazard", col="red") points(Pred_Modbs2_2A, which="hazard", type="l", conf.int=FALSE)
Function for predicting the (excess) hazard and the
corresponding (net) survival from a model fitted with the mexhaz
function for a particular vector of covariates. If the survival model
was fitted with an expected hazard (excess hazard model), the estimates
obtained are excess hazard and net survival estimates. Corresponding
variance estimates are based on the Delta Method or Monte Carlo
simulation (based on the assumption of multivariate normality of the
model parameter estimates). This function allows the computation of the
hazard and the survival at one time point for several vectors of
covariates or for one vector of covariates at several time
points. When the model includes a random effect, three types of
predictions can be made: (i) marginal predictions (obtained by
integration over the random effect distribution), (ii) cluster-specific
posterior predictions for an existing cluster, or (iii) conditional
predictions for a given quantile of the random effect distribution (by
default, for the median value, that is, 0).
## S3 method for class 'mexhaz' predict(object, time.pts, data.val = data.frame(.NotUsed = NA), marginal = FALSE, quant.rdm = 0.5, cluster = NULL, conf.int = c("delta", "simul", "none"), level = 0.95, delta.type.h = c("log", "plain"), delta.type.s = c("log-log", "log", "plain"), nb.sim = 10000, keep.sim = FALSE, include.gradient = FALSE, dataset = NULL, ...)## S3 method for class 'mexhaz' predict(object, time.pts, data.val = data.frame(.NotUsed = NA), marginal = FALSE, quant.rdm = 0.5, cluster = NULL, conf.int = c("delta", "simul", "none"), level = 0.95, delta.type.h = c("log", "plain"), delta.type.s = c("log-log", "log", "plain"), nb.sim = 10000, keep.sim = FALSE, include.gradient = FALSE, dataset = NULL, ...)
object |
an object of class |
time.pts |
a vector of numerical values representing the time points at which predictions are requested. Time values greater than the maximum follow-up time on which the model estimation was based are discarded. |
data.val |
a |
marginal |
logical value controlling the type of predictions
returned by the function when the model includes a random
intercept. When |
quant.rdm |
numerical value (between 0 and 1) specifying the
quantile of the random effect distribution that should be used when
requesting conditional predictions. The default value is set to 0.5
(corresponding to the median, that is a value of the random effect
of 0). This argument is ignored if the model is a fixed effect
model, when the |
cluster |
a single value corresponding to the name of the cluster for
which posterior predictions should be calculated. These predictions
are obtained by integrating over the cluster-specific posterior
distribution of the random effect and thus require the original
dataset. The dataset can either be provided as part of the
|
conf.int |
method to be used to compute confidence
limits. Selection can be made between the following options:
|
level |
a number in (0,1) specifying the level of confidence for
computing the confidence intervals of the hazard and the survival. The
default value is set to |
delta.type.h |
type of confidence limits for the hazard when
using the Delta Method. With the default value ( |
delta.type.s |
type of confidence limits for the survival when
using the Delta Method. With the default value ( |
nb.sim |
integer value representing the number of simulations
used to estimate the confidence limits for the (excess) hazard
and the (net) survival. This argument is used only if
|
keep.sim |
logical value determining if the simulated hazard and
survival values should be returned (only used when
|
include.gradient |
logical value allowing the function to return
the components of the gradient of the logarithm of the hazard and of
the logarithm of the cumulative hazard for each prediction. This
argument is used only if |
dataset |
original dataset used to fit the |
... |
for potential additional parameters. |
An object of class predMexhaz that can be used by
the function plot.predMexhaz to produce graphics of the (excess) hazard and
the (net) survival. It contains the following elements:
call |
the |
results |
a |
variances |
a |
grad.loghaz |
a |
grad.logcum |
a |
vcov |
a matrix corresponding to the covariance matrix used to compute the confidence intervals. |
type |
this item can take the value
|
type.me |
the type of predictions produced in case of a model
including a random intercept. Can take the values |
ci.method |
the method used to compute confidence limits. |
level |
level of confidence used to compute confidence limits. |
delta.type |
type of confidence limits for the hazard and the survival when using the Delta Method. |
nb.sim |
number of simulations used to estimate the confidence
limits when |
sim.haz |
matrix containing the simulated hazards (each column
representing a simulated vector of values); only returned when
|
sim.surv |
matrix containing the simulated survival probabilities (each column
representing a simulated vector of values); only returned when
|
Hadrien Charvat, Aurelien Belot
Charvat H, Remontet L, Bossard N, Roche L, Dejardin O, Rachet B, Launoy G, Belot A; CENSUR Working Survival Group. A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates. Stat Med 2016;35:3066-3084 (doi: 10.1002/sim.6881).
Skrondal A, Rabe-Hesketh S. Prediction in multilevel generalized linear models. J R Stat Soc A Stat Soc 2009;172(3):659-687 (doi: 10.1111/j.1467-985X.2009.00587.x).
print.predMexhaz, plot.predMexhaz,
points.predMexhaz, lines.predMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and ## with effects of age (agecr), deprivation index (depindex) and sex ## (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs1_2A <- predict(Mod_bs1_2, time.pts=seq(0.1,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) ## Prediction for several vectors of covariates at one time point Pred_Modbs1_2B <- predict(Mod_bs1_2, time.pts=10, data.val=data.frame(agecr=c(-0.2,-0.1,0), depindex=c(0.5,0.5,0.5), IsexH=c(1,1,1))) ## Prediction for all individuals of the study population ## at one time point Pred_Modbs1_2C <- predict(Mod_bs1_2, time.pts=10, data.val=simdatn1) # Example of cluster-specific posterior prediction (not run) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year # Mod_bs3_2mix <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+IsexH, data=simdatn1, base="exp.bs", degree=3, # knots=c(1,5), expected="popmrate", random="clust", verbose=1000) ## Posterior predictions at several time points for an individual ## in cluster 15 with a specific vector of covariates # Pred_Modbs3_2A <- predict(Mod_bs3_2mix, # time.pts=seq(0.1,10,by=0.1), data.val=data.frame(agecr=0.2, IsexH=1), # cluster=15, dataset=simdatn1)data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and ## with effects of age (agecr), deprivation index (depindex) and sex ## (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs1_2A <- predict(Mod_bs1_2, time.pts=seq(0.1,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1)) ## Prediction for several vectors of covariates at one time point Pred_Modbs1_2B <- predict(Mod_bs1_2, time.pts=10, data.val=data.frame(agecr=c(-0.2,-0.1,0), depindex=c(0.5,0.5,0.5), IsexH=c(1,1,1))) ## Prediction for all individuals of the study population ## at one time point Pred_Modbs1_2C <- predict(Mod_bs1_2, time.pts=10, data.val=simdatn1) # Example of cluster-specific posterior prediction (not run) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a cubic B-spline with two knots at 1 and 5 year # Mod_bs3_2mix <- mexhaz(formula=Surv(time=timesurv, # event=vstat)~agecr+IsexH, data=simdatn1, base="exp.bs", degree=3, # knots=c(1,5), expected="popmrate", random="clust", verbose=1000) ## Posterior predictions at several time points for an individual ## in cluster 15 with a specific vector of covariates # Pred_Modbs3_2A <- predict(Mod_bs3_2mix, # time.pts=seq(0.1,10,by=0.1), data.val=data.frame(agecr=0.2, IsexH=1), # cluster=15, dataset=simdatn1)
Display the model call as well as the values of the estimated model parameters.
## S3 method for class 'mexhaz' print(x, ...)## S3 method for class 'mexhaz' print(x, ...)
x |
an object of class |
... |
represents additional parameters directly passed to |
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") print(Mod_weib_mix)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") print(Mod_weib_mix)
Display the first lines of the data.frame containing the
predictions provided by the predMexhaz function.
## S3 method for class 'predMexhaz' print(x, ...)## S3 method for class 'predMexhaz' print(x, ...)
x |
an object of class |
... |
represents additional parameters directly passed to |
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs1_2A <- predict(Mod_bs1_2, time.pts=seq(0.1,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1), conf.int="delta") print(Pred_Modbs1_2A)data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Prediction at several time points for one vector of covariates Pred_Modbs1_2A <- predict(Mod_bs1_2, time.pts=seq(0.1,10,by=0.1), data.val=data.frame(agecr=0,depindex=0.5,IsexH=1), conf.int="delta") print(Pred_Modbs1_2A)
Display the model call, the values of the estimated model parameters, as well as the corresponding hazard ratios (only for proportional effects).
## S3 method for class 'summary.mexhaz' print(x, ...)## S3 method for class 'summary.mexhaz' print(x, ...)
x |
an object of class |
... |
represents additional parameters directly passed to |
mexhaz, print.mexhaz, summary.mexhaz
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", verbose=0) summary(Mod_weib)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", verbose=0) summary(Mod_weib)
This is a generic function.
ranef(x, ...)ranef(x, ...)
x |
a fitted object from which random effects can be extracted. |
... |
might be used by some methods using this generic function. |
see the documentation for the specific class.
## See the specific class documentation.## See the specific class documentation.
Display a data frame containing the cluster-specific random effects with their standard errors.
## S3 method for class 'mexhaz' ranef(x, ...)## S3 method for class 'mexhaz' ranef(x, ...)
x |
an object of class |
... |
not used. |
A data frame with three columns containing the cluster names, the random effect estimates, and their standard errors.
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") ranef(Mod_weib_mix)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") ranef(Mod_weib_mix)
Function for computing hazard ratio and risk ratio (ratio
of cumulative probabilities of failure) estimates from a model fitted
with the mexhaz function. Corresponding confidence intervals
are based on the Delta Method or Monte Carlo simulation (based on the
assumption of multivariate normality of the model parameter
estimates). This function allows the computation of estimates at one
time point for several vectors of covariates or for one vector of
covariates at several time points. When the model includes a random effect, three types of
predictions can be made: (i) marginal predictions (obtained by
integration over the random effect distribution), (ii) cluster-specific
posterior predictions for an existing cluster, or (iii) conditional
predictions for a given quantile of the random effect distribution (by
default, for the median value, that is, 0).
riskfunc(object, time.pts, data, data.0, marginal = TRUE, quant.rdm = 0.5, cluster = NULL, quant.rdm.0 = 0.5, cluster.0 = NULL, type = c("hr", "rr"), conf.int = c("delta", "simul"), level = 0.95, nb.sim = 10000, seed = NULL, dataset = NULL)riskfunc(object, time.pts, data, data.0, marginal = TRUE, quant.rdm = 0.5, cluster = NULL, quant.rdm.0 = 0.5, cluster.0 = NULL, type = c("hr", "rr"), conf.int = c("delta", "simul"), level = 0.95, nb.sim = 10000, seed = NULL, dataset = NULL)
object |
an object of class |
time.pts |
a vector of numerical values representing the time points at which predictions are requested. Time values greater than the maximum follow-up time on which the model estimation was based are discarded. |
data |
a |
data.0 |
a |
marginal |
logical value controlling the type of predictions
returned by the function when the model includes a random
intercept. When |
quant.rdm |
numerical value (between 0 and 1) specifying the
quantile of the random effect distribution that should be used when
requesting conditional predictions. The default value is set to 0.5
(corresponding to the median, that is a value of the random effect
of 0). This argument is ignored if the model is a fixed effect
model, when the |
cluster |
a single value corresponding to the name of the cluster for
which posterior predictions should be calculated. These predictions
are obtained by integrating over the cluster-specific posterior
distribution of the random effect and thus require the original
dataset. The dataset can either be provided as part of the
|
quant.rdm.0 |
random effect distribution quantile value to be used with |
cluster.0 |
cluster value to be used with |
type |
argument specifying the type of predictions to be
calculated. Selection can be made between |
conf.int |
method to be used to compute confidence
limits. Selection can be made between the following options:
|
level |
a number in (0,1) specifying the level of confidence for
computing the confidence intervals of the hazard and the
survival. By default, this argument is set to |
nb.sim |
integer value representing the number of simulations
used to estimate the confidence limits for the (excess) hazard
and the (net) survival. This argument is used only if
|
seed |
argument allowing the user to set a random seed for
simulations (only relevant when |
dataset |
original dataset used to fit the |
An object of class resMexhaz that can be used by the function
plot.resMexhaz to produce graphics
of the hazard ratio or risk ratio curve. It contains the following
elements:
results |
a |
type |
type of results returned by the function. The value is
used by |
multiobs |
value used by
|
ci.method |
method used to compute confidence limits. |
level |
level of confidence used to compute confidence limits. |
Hadrien Charvat, Aurelien Belot
Charvat H, Remontet L, Bossard N, Roche L, Dejardin O, Rachet B, Launoy G, Belot A; CENSUR Working Survival Group. A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates. Stat Med 2016;35:3066-3084 (doi: 10.1002/sim.6881)
Skrondal A, Rabe-Hesketh S. Prediction in multilevel generalized linear models. J R Stat Soc A Stat Soc 2009;172(3):659-687 (doi: 10.1111/j.1467-985X.2009.00587.x).
plot.resMexhaz, lines.resMexhaz
data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Risk ratio along time for agecr=0.2 compared to agecr=0.1 RR_Modbs1_2 <- riskfunc(Mod_bs1_2, time.pts=seq(0,10,le=101), data=data.frame(agecr=0.2,depindex=0,IsexH=1), data.0=data.frame(agecr=0.1,depindex=0,IsexH=1),type="rr", conf.int="delta") plot(RR_Modbs1_2)data(simdatn1) ## Fit of a fixed-effect hazard model, with the baseline hazard ## described by a linear B-spline with two knots at 1 and 5 year and with ## effects of age (agecr), deprivation index (depindex) and sex (IsexH) Mod_bs1_2 <- mexhaz(formula=Surv(time=timesurv, event=vstat)~agecr+depindex+IsexH, data=simdatn1, base="exp.bs", degree=1, knots=c(1,5), verbose=0) ## Risk ratio along time for agecr=0.2 compared to agecr=0.1 RR_Modbs1_2 <- riskfunc(Mod_bs1_2, time.pts=seq(0,10,le=101), data=data.frame(agecr=0.2,depindex=0,IsexH=1), data.0=data.frame(agecr=0.1,depindex=0,IsexH=1),type="rr", conf.int="delta") plot(RR_Modbs1_2)
The simdatn1 dataset has 4000 rows and 8 columns.
This dataset contains the following columns:
Age at diagnosis (continuous).
Centred and rescaled age variable (age-70)/100.
Deprivation index (continuous).
Sex (0 = Female, 1 = Male).
ID number of the cluster.
Vital status (0 = Alive, 1 = Dead).
Follow-up time (years).
Population (expected) mortality rate at the time of censoring.
The simdatn2 dataset has 4000 rows and 8 columns. The
IsexH variable is simulated with a non-proportional effect.
This dataset contains the following columns:
Age at diagnosis (continuous).
Centred and rescaled age variable (age-70)/100.
Deprivation index (continuous).
Sex (0 = Female, 1 = Male).
ID number of the cluster.
Vital status (0 = Alive, 1 = Dead).
Follow-up time (years).
Population (expected) mortality rate at the time of censoring.
Produces a summary of a mexhaz object.
## S3 method for class 'mexhaz' summary(object, ...)## S3 method for class 'mexhaz' summary(object, ...)
object |
an object of class |
... |
represents additional parameters directly passed to |
mexhaz, print.mexhaz, print.summary.mexhaz
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", verbose=0) summary(Mod_weib)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", verbose=0) summary(Mod_weib)
Function allowing the user to update an existing mexhaz model. All the arguments of the model can be updated. If the argument 'init' is not provided, the function uses the estimated values of the existing model as starting values for the corresponding parameters of the new model.
## S3 method for class 'mexhaz' update(object, formula, data, expected = NULL, base = NULL, degree = 3, knots = NULL, bound = NULL, n.gleg = 20, init = NULL, random = NULL, n.aghq = 10, fnoptim = c("nlm", "optim"), verbose = 0, method = "Nelder-Mead", iterlim = 10000, numHess = FALSE, print.level = 1, exactGradHess = TRUE, gradtol = ifelse(exactGradHess, 1e-8, 1e-6), envir = parent.frame(), ...)## S3 method for class 'mexhaz' update(object, formula, data, expected = NULL, base = NULL, degree = 3, knots = NULL, bound = NULL, n.gleg = 20, init = NULL, random = NULL, n.aghq = 10, fnoptim = c("nlm", "optim"), verbose = 0, method = "Nelder-Mead", iterlim = 10000, numHess = FALSE, print.level = 1, exactGradHess = TRUE, gradtol = ifelse(exactGradHess, 1e-8, 1e-6), envir = parent.frame(), ...)
object |
an object of class |
formula |
a formula object, with the response on the left of the In case See |
data |
a |
expected |
name of the variable (must be given in quotes) representing the
population (i.e., expected) hazard. By default, |
base |
functional form that should be used to model the baseline
hazard. Selection can be made between the following options:
|
degree |
if |
knots |
if |
bound |
If |
n.gleg |
if |
init |
vector of initial values. By default for the baseline hazard: if if - if - if - if - if the parameters describing the effects of the covariates are all set to 0; the parameter representing the standard deviation of the random
effect is set to 0.1. (if |
random |
name of the variable to be entered as a random effect (must be given
between quotes), representing the cluster membership. By default,
|
n.aghq |
number of quadrature points used for estimating the
cluster-specific marginal likelihoods by adaptive Gauss-Hermite
quadrature. By default, |
fnoptim |
name of the R optimisation procedure used to maximise the
likelihood. Selection can be made between |
verbose |
integer parameter representing the frequency at which the current state
of the optimisation process is displayed. Internally, an 'evaluation' is
defined as an estimation of the log-likelihood for a given vector of
parameters. This means that the number of evaluations is increased each
time the optimisation procedure updates the value of any of the
parameters to be estimated. If |
method |
if |
iterlim |
if |
numHess |
logical value allowing the user to choose between the Hessian returned
by the optimization algorithm (default) or the Hessian estimated by
the |
print.level |
this argument is only used if |
exactGradHess |
logical value allowing the user to decide
whether maximisation of the likelihood should be based on the
analytic gradient and Hessian computed internally (default,
corresponding to |
gradtol |
this argument is only used if |
envir |
environment in which the objects' names given as arguments to the updated model are to be found. |
... |
represents additional parameters directly passed to |
An object of class mexhaz. See mexhaz for more details.
Hadrien Charvat, Aurelien Belot
Charvat H, Remontet L, Bossard N, Roche L, Dejardin O, Rachet B, Launoy G, Belot A; CENSUR Working Survival Group. A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates. Stat Med 2016;35:3066-3084 (doi: 10.1002/sim.6881)
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", verbose=0) ## Add an effect of gender Mod_weib_2 <- update(Mod_weib, formula=~.+IsexH)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", verbose=0) ## Add an effect of gender Mod_weib_2 <- update(Mod_weib, formula=~.+IsexH)
Display a the covariance matrix of the fixed effects of a mexhaz model.
## S3 method for class 'mexhaz' vcov(object, ...)## S3 method for class 'mexhaz' vcov(object, ...)
object |
an object of class |
... |
not used. |
The estimated covariance matrix of the fixed effects.
data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") vcov(Mod_weib_mix)data(simdatn1) ## Fit of a mixed-effect excess hazard model, with the baseline hazard ## described by a Weibull distribution (without covariables) Mod_weib_mix <- mexhaz(formula=Surv(time=timesurv, event=vstat)~1, data=simdatn1, base="weibull", expected="popmrate", verbose=0, random="clust") vcov(Mod_weib_mix)