Fit an (excess) hazard regression model using different shapes for the
baseline hazard (Weibull, piecewise constant, exponential of a Bspline
of degree 1 to 3, exponential of a restricted cubic spline), with the
possibility to include timedependent effects of variable(s) and a
random effect defined at the cluster level. The function accepts
rightcensored and counting process input styles for the followup
time. The latter allows the modelling of survival data with delayed
entries. The timedependent 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 Bspline 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 clusterspecific marginal
likelihoods. The logarithm of the full marginal likelihood, defined as
the sum of the logarithms of the clusterspecific marginal likelihoods,
is then maximised using an optimisation routine such as nlm
(default) or optim
.
1 2 3 4 5  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,
fnoptim = c("nlm", "optim"), verbose = 100, method = "NelderMead",
iterlim=10000, numHess=FALSE, print.level=1,...)

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:

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 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. 
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 to be used for estimating the
clusterspecific marginal likelihoods by adaptive GaussHermite
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 loglikelihood 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 
... 
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 followup 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. 
bounds 
vector of boundary values used to define the Bspline (or natural spline) bases. 
degree 
degree of the Bspline used to model the logarithm of the baseline hazard. 
knots 
vector of interior knots used to define the Bspline (or natural spline) bases. 
names.ph 
names of the covariables with a proportional effect. 
random 
name of the variable defining cluster membership (set to

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 variancecovariance matrix of the estimated parameters. 
mu.hat 
a 
var.mu.hat 
the covariance matrix of the clusterspecific shrinkage estimators. 
vcov.fix.mu.hat 
a 
n.par 
number of estimated parameters. 
n.gleg 
number of GaussLegendre quadrature points used to calculate the cumulative (excess) hazard (only relevant if a Bspline 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 GaussHermite quadrature points used to calculate the clusterspecific marginal likelihoods (only relevant if a multilevel 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 loglikelihood 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 nonlinear and nonproportional effects of covariates. Stat Med 2016;35:30663084 (doi: 10.1002/sim.6881)
print.mexhaz
, summary.mexhaz
, predict.mexhaz
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30  data(simdatn1)
## Fit of a mixedeffect 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 mixedeffect excess hazard model, with the baseline hazard
## described by a cubic Bspline 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 mixedeffect excess hazard model, with the baseline hazard
## described by a restricted cubic spline with two knots at the
## tertiles of the distribution of followup 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)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
Please suggest features or report bugs with the GitHub issue tracker.
All documentation is copyright its authors; we didn't write any of that.