update.mexhaz: Update of a mexhaz model

update.mexhazR Documentation

Update of a mexhaz model

Description

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.

Usage

## 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(), ...)

Arguments

object

an object of class mexhaz, corresponding to a survival model fitted with the mexhaz function.

formula

a formula object, with the response on the left of the ~ operator, and the linear predictor on the right. The response must be of the form Surv(time, event) for right censored data or Surv(time, time2, event) for counting process style data. The linear predictor accepts a special instruction nph() for specifying variables for which a time-dependent effect should be modelled (if several variables are modelled with time-dependent effects, separate these variables inside the nph() instruction with a + sign).

In case time takes the value 0 for some observations when data are entered in the right censored style, it is assumed that these observations refer to events (or censoring) that occurred on the first day of follow-up. Consequently, a value of 1/730.5 (half a day) is substituted in order to make computations possible. However, it should be stressed that this is just a convention and that it does not make much sense if the time scale is not expressed in years. We therefore advise the analyst to deal with 0 time values during the dataset preparation stage.

See update.formula for more details.

data

a data.frame containing the variables referred to in the formula, as well as in the expected and random arguments if these arguments are used.

expected

name of the variable (must be given in quotes) representing the population (i.e., expected) hazard. By default, expected=NULL, which means that the function estimates the overall hazard (and not the excess hazard).

base

functional form that should be used to model the baseline hazard. Selection can be made between the following options: "weibull" for a Weibull hazard, "exp.bs" for a hazard described by the exponential of a B-spline (only B-splines of degree 1, 2 or 3 are accepted), "exp.ns" for a hazard described by the exponential of a restricted cubic B-spline (also called 'natural spline'), "pw.cst" for a piecewise constant hazard. By default, base=NULL.

degree

if base="exp.bs", degree represents the degree of the B-spline used. Only integer values between 1 and 3 are accepted, and 3 is the default.

knots

if base="exp.bs" or "exp.ns", knots is the vector of interior knots of the spline. If base="pw.cst", knots is the vector defining the endpoints of the time intervals on which the hazard is assumed to be constant. By default, knots=NULL (that is, it produces a B-spline with no interior knots if base="exp.bs", a linear B-spline with no interior knots if base="exp.ns", or a constant hazard over the whole follow-up period if base="pw.cst").

bound

If base="exp.bs" or base="exp.ns", computation of the B-spline basis requires that boundary knots be given. The bound argument allows the user to specify these boundary knots. If base="exp.bs", the interval defined by the boundary knots must at least include the interval c(0,max(time)) (otherwise, there could be problems with ill-conditioned bases). If base="exp.ns", the boundary knots correspond to the knots beyond which the spline is contrained to be linear (in that case, the boundary knots can be values contained in c(0,max(time))). By default, the boundary knots are set to c(0,max(time)).

n.gleg

if base="exp.bs" and degree is equal to 2 or 3, or if base="exp.ns", the cumulative hazard is computed via Gauss-Legendre quadrature and n.gleg is the number of quadrature nodes to be used to compute the cumulative hazard. By default, n.gleg=20.

init

vector of initial values. By default init=NULL and the initial values are internally set to the following values:

for the baseline hazard:

if exactGradHess=TRUE (except for the excess hazard random effect model for which this argument is ignored), the intercept is set to 0.5*log(Number of events/Person-years of observation) and all other parameters set to 0. In case of failed convergence, several trials are run with an adaptation of the value of the intercept.

if exactGradHess=FALSE, the following values are used:

- if base="weibull", the scale and shape parameters are set to 0.1;

- if base="exp.bs", the parameters of the B-spline are all set to -1;

- if base="exp.ns", the parameters of the restricted cubic B-spline are all set to -1;

- if base="pw.cst", the logarithm of the piecewise-constant hazards are set to -1;

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 exactGradHess=TRUE, several trials are run with an adaptation of the value in case of failed convergence).

random

name of the variable to be entered as a random effect (must be given between quotes), representing the cluster membership. By default, random=NULL which means that the function fits a fixed effects model.

n.aghq

number of quadrature points used for estimating the cluster-specific marginal likelihoods by adaptive Gauss-Hermite quadrature. By default, n.aghq=10.

fnoptim

name of the R optimisation procedure used to maximise the likelihood. Selection can be made between "nlm" (by default) and "optim". Note: if exactGradHess=TRUE, this argument will be ignored (fnoptim will be set automatically to "nlm").

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 verbose=n (with n an integer), the function will display the current values of the parameters, the log-likelihood and the time elapsed every n evaluations. If verbose=0, nothing is displayed.

method

if fnoptim="optim", method represents the optimisation method to be used by optim. By default, method="Nelder-Mead". This parameter is not used if fnoptim="nlm".

iterlim

if fnoptim="nlm", iterlim represents the maximum number of iterations before the nlm optimisation procedure is terminated. By default, iterlim is set to 10000. This parameter is not used if fnoptim="optim" (in this case, the maximum number of iterations must be given as part of a list of control parameters via the control argument: see the help page of optim for further details).

numHess

logical value allowing the user to choose between the Hessian returned by the optimization algorithm (default) or the Hessian estimated by the hessian function from the numDeriv package. The latter might be more accurate but its estimation is more time-consuming. We suggest to use the default Hessian estimation procedure during model building and estimate the numDeriv-based Hessian only on the final model. Note: if exactGradHess=TRUE, this argument is ignored.

print.level

this argument is only used if fnoptim="nlm". It determines the level of printing during the optimisation process. The default value (for the mexhaz function) is set to '1' which means that details on the initial and final step of the optimisation procedure are printed (see the help page of nlm for further details).

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 exactGradHess=TRUE). In that case, optimisation is performed with the nlm function. Note: even if set to TRUE, this argument is ignored when the user wants to fit an excess hazard model including a random effect because in that case, there is no simple way to obtain the analytic gradient and Hessian.

gradtol

this argument is only used if fnoptim="nlm". It corresponds to the tolerance at which the scaled gradient is considered close enough to zero to terminate the algorithm. The default value depends on the value of the argument exactGradHess.

envir

environment in which the objects' names given as arguments to the updated model are to be found.

...

represents additional parameters directly passed to nlm or optim to control the optimisation process.

Value

An object of class mexhaz. See mexhaz for more details.

Author(s)

Hadrien Charvat, Aurelien Belot

References

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)

See Also

mexhaz

Examples


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)


mexhaz documentation built on Oct. 31, 2022, 5:08 p.m.