Simultaneously estimate the regression coefficients and the baseline hazard function of proportional hazard Cox models using maximum penalised likelihood (MPL).
coxph_mpl(formula, data, subset, na.action, control, …)# S3 method for coxph_mpl
print(x,…)
a formula object, with the response on the left of a ~
operator, and
the terms on the right. The response must be a survival object as
returned by the Surv
function.
a data.frame in which to interpret the variables named in
the formula
, or in the subset
argument. If no dataset is indicated,
variables will be taken from the global environment.
expression indicating which subset of the rows of data should be used in the fit. All observations are included by default.
a missing-data filter function. This is applied to the model.frame
after any subset argument has been used. Default is options()\$na.action
.
an object inheriting from class coxph_mpl
, representing
a fitted Cox proportional hazard model.
Object of class coxph_mpl.control
specifying control options like
basis choice, smoothing parameter value and maximum number of itereations,
for example. Refer to coxph_mpl.control
to see the defaults.
Other arguments. In coxph_mpl
, these elements,
will be passed to coxph_mpl.control
. In print.coxph_mpl
,
these elements will be passed to the print
function.
an object of class coxph_mpl
representing the fit.
See coxph_mpl.object
for details.
coxph_mpl
allows to simultaneously estimate the regression
coefficients and baseline hazard function of Cox proportional hazard models,
with right censored data and independent censoring, by maximising a penalised
likelihood, in which a penalty function is used to smooth the baseline hazard
estimate.
Optimisation is achieved using a new iterative algorithm, which combines Newton's method and the multiplicative iterative algorithm proposed by Ma (2010), and respects the non-negativity constraints on the baseline hazard estimate (refer to Jun, Heritier and Lo (2014)).
The centered X matrix is used in the optimisation process to get a better shaped (penalised) log-likelihood. Baseline hazard parameter estimates and covariance matrix are then respectively corrected using a correction factor and the delta method.
When the chosen basis is not uniform, estimates of zero are possible for baseline hazard parameters and will correspond to active constraints as defined by Moore and Sadler (2008). Inference, as described by Ma, Heritier and Lo (2014), is then corrected accordingly (refer to Moore and Sadler (2008)) by adequately 'cutting' the corresponding covariance matrix.
There are currently 3 ways to perform inference on model parameters:
Let \(H\) denote the Hessian matrix of the unpenalised likelihood, \(Q\) denote the product of the first order derivative of the penalised likelihood by its transpose, and \(M_{2}\) denote the second order derivative of the penalised likelihood. Then,
'H'
refers to \(H^{-1}\), the inverse of the Hessian matrix,
'M2QM2'
, refers to the sandwich formula \(M_{2}^{-1} Q M_{2}^{-1}\),
'M2HM2'
, refers to the sandwich formula \(M_{2}^{-1} H M_{2}^{-1}\).
Simulations analysing the coverage levels of confidence intervals for the regression parameters seem to indicate that \(M_{2}^{-1} H M_{2}^{-1}\) performs better when using the uniform basis, and that \(M_{2}^{-1} Q M_{2}^{-1}\) performs when using other bases.
Ma, J. and Heritier, S. and Lo, S. (2014), On the Maximum Penalised Likelihood Approach for Proportional Hazard Models with Right Censored Survival Data. Computational Statistics and Data Analysis 74, 142-156.
Ma, J. (2010), Positively constrained multiplicative iterative algorithm for maximum penalised likelihood tomographic reconstruction. IEEE Transactions On Signal Processing 57, 181-192.
Moore, T. J. and Sadler, B. M. and Kozick R. J. (2008), Maximum-Likelihood Estimation, the Cramer-Rao Bound, and the Method of Scoring With Parameter Constraints, IEEE Transactions On Signal Processing 56, 3, 895-907.
coxph_mpl.object
, coxph_mpl.control
,
summary.coxph_mpl
and plot.coxph_mpl
.
# NOT RUN {
## data lung of the survival package. Refer to ?lung.
data(lung)
fit_mpl <- coxph_mpl(Surv(time, status == 2) ~ age + sex + ph.karno + wt.loss, data = lung)
fit_mpl
# }
Run the code above in your browser using DataLab