Title: | Weighted Estimation in Cox Regression |
---|---|
Description: | Implements weighted estimation in Cox regression as proposed by Schemper, Wakounig and Heinze (Statistics in Medicine, 2009, <doi:10.1002/sim.3623>) and as described in Dunkler, Ploner, Schemper and Heinze (Journal of Statistical Software, 2018, <doi:10.18637/jss.v084.i02>). Weighted Cox regression provides unbiased average hazard ratio estimates also in case of non-proportional hazards. Approximated generalized concordance probability an effect size measure for clear-cut decisions can be obtained. The package provides options to estimate time-dependent effects conveniently by including interactions of covariates with arbitrary functions of time, with or without making use of the weighting option. |
Authors: | Daniela Dunkler [aut, cre], Georg Heinze [aut], Meinhard Ploner [aut] |
Maintainer: | Daniela Dunkler <[email protected]> |
License: | GPL-3 |
Version: | 4.0.3 |
Built: | 2024-11-23 06:32:14 UTC |
Source: | CRAN |
This package implements weighted estimation in Cox regression as proposed by Schemper, Wakounig and Heinze (Statistics in Medicine, 2009, doi:10.1002/sim.3623). Weighted Cox regression provides unbiased average hazard ratio estimates also in case of non-proportional hazards. The package provides options to estimate time-dependent effects conveniently by including interactions of covariates with arbitrary functions of time, with or without making use of the weighting option. For more details we refer to Dunkler, Ploner, Schemper and Heinze (Journal of Statistical Software, 2018, doi:10.18637/jss.v084.i02).
Package: | coxphw
|
Type: | Package |
Version: | 4.0.2 |
Date: | 2020-06-16 |
License: | GPL-2 |
Main functions included in the coxphw package are
coxphw |
weighted estimation of Cox regression: either (recommended) estimation of |
average hazard ratios (Schemper et al., 2009), estimation of average regression | |
effects (Xu and O'Quigley, 2000), or proportional hazards regression. | |
plot |
plots the weights used in a weighted Cox regression analysis against time. |
concord |
obtains generalized concordance probabilities with confidence intervalls. |
predict |
obtains the effect estimates (of e.g. a nonlinear or a time-dependent effect) |
at specified values of a continuous covariable. With plot.coxphw.predict
|
|
these relative or log relative hazard versus values of the continuous covariable | |
can be plotted. | |
wald |
obtain Wald chi-squared test statistics and p-values for one or more regression |
coefficients given their variance-covariance matrix. | |
Data sets included in the coxphw package are
biofeedback |
biofeedback treatment data |
gastric |
gastric cancer data |
The SAS
macro WCM
with similar functionality can be obtained at
https://cemsiis.meduniwien.ac.at/en/kb/science-research/software/statistical-software/wcmcoxphw/.
Important version changes:
Up to Version 2.13 coxphw used a slightly different syntax (arguments: AHR
, AHR.norobust
,
ARE
, PH
, normalize
, censcorr
, prentice
, breslow
,
taroneware
).
From Version 3.0.0 on the old syntax is disabled.
From Version 4.0.0 fractional polynomials are disabled and plotshape
is replaced with predict
and plot.coxphw.predict
.
Georg Heinze, Meinhard Ploner, Daniela Dunkler
Maintainer: [email protected]
Dunkler D, Ploner M, Schemper M, Heinze G. (2018) Weighted Cox Regression Using the R Package coxphw. JSS 84, 1–26, doi:10.18637/jss.v084.i02.
Dunkler D, Schemper M, Heinze G. (2010) Gene Selection in Microarray Survival Studies Under Possibly Non-Proportional Hazards. Bioinformatics 26:784-90.
Lin D and Wei L (1989). The Robust Inference for the Cox Proportional Hazards Model. J AM STAT ASSOC 84, 1074-1078.
Lin D (1991). Goodness-of-Fit Analysis for the Cox Regression Model Based on a Class of Parameter Estimators. J AM STAT ASSOC 86, 725-728.
Royston P and Altman D (1994). Regression Using Fractional Polynomials of Continuous Covariates: Parsimonious Parametric Modelling. J R STAT SOC C-APPL 43, 429-467.
Royston P and Sauerbrei W (2008). Multivariable Model-Building. A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables. Wiley, Chichester, UK.
Sasieni P (1993). Maximum Weighted Partial Likelihood Estimators for the Cox Model. J AM STAT ASSOC 88, 144-152.
Schemper M (1992). Cox Analysis of Survival Data with Non-Proportional Hazard Functions. J R STAT SOC D 41, 455-465.
Schemper M, Wakounig S and Heinze G (2009). The Estimation of Average Hazard Ratios by Weighted Cox Regression. Stat Med 28, 2473-2489, doi:10.1002/sim.3623.
Xu R and O'Quigley J (2000). Estimating Average Regression Effect Under Non-Proportional Hazards. Biostatistics 1, 423-439.
coxphw
, concord
, plot.coxphw
, predict.coxphw
, plot.coxphw.predict
, wald
## for examples see coxphw
## for examples see coxphw
In this study the effect of biofeedback treatment on time until treatment success
was evaluated in patients suffering from aspiration after head and neck
surgery. The outcome of interest was the time from start of treatment until the patient achieved full
swallowing rehabilitation (thdur
). Patients were randomized into two groups (bfb
):
one group of patients received videoendoscopic biofeedback treatment; the other group received the
conservative treatment including thermal stimulation with ice and exercises for the lips, tongue,
laryngeal closure and elevation. Treatment was started as soon as the healing process after surgery
was finished (thealing
).
data(biofeedback)
data(biofeedback)
A data frame with 33 observations on the following 6 variables:
id
the patient id.
success
of treatment within the first 100 days; either 0 = no success or 1 = success.
thdur
the duration of therapy in days.
bfb
indicates the treatment group; either 0 = conservative or 1 = biofeedback.
theal
time from surgery to start of therapy in days.
log2heal
log2-transformed time from surgery to start of therapy.
Data were supplied by Dr. D.-M. Denk, who gave permission to freely distribute the data and use them for non-commercial purposes.
Denk, D.-M. & Kaider, A. (1997). Videoendoscopic Biofeedback: A Simple Method to Improve the Efficacy of Swallowing Rehabilitation of Patients After Head and Neck Surgery. ORL J OTO-RHINO-LARY 59, 100-105.
data("biofeedback") plot(survfit(Surv(thdur, success) ~ bfb, data = biofeedback), lty = 1:2, las = 1, xlab = "time (days)", ylab = "propability of success") coxphw(Surv(thdur, success) ~ bfb, data = biofeedback, template = "AHR")
data("biofeedback") plot(survfit(Surv(thdur, success) ~ bfb, data = biofeedback), lty = 1:2, las = 1, xlab = "time (days)", ylab = "propability of success") coxphw(Surv(thdur, success) ~ bfb, data = biofeedback, template = "AHR")
coxphw
This class of objects is returned by the coxphw
function.
Objects of this class have methods for the functions summary
, print
, coef
,
vcov
, plot
and confint
.
## S3 method for class 'coxphw' coef(object, ...)
## S3 method for class 'coxphw' coef(object, ...)
object |
object of class |
... |
further arguments. |
Daniela Dunkler
coxphw
or
coxph
Compute generalized concordance probabilities with accompanying
confidence intervals for objects of class coxphw
or coxph
.
concord(fit, digits = 4)
concord(fit, digits = 4)
fit |
an object of class |
digits |
integer indicating the number of decimal places to be used. Default is 4. |
The generalized concordance probability is defined as
with
and
as survival times of randomly chosen subjects with covariate
values
and
, respectively. Assuming that
and
are
1 and 0, respectively, this definition includes a two-group comparison.
If proportional hazards can be assumed, the generalized concordance probability can also
be derived from Cox proportional hazards regression (coxphw
with template = "PH"
or coxph
) or weighted Cox regression as suggested by Xu and O'Quigley (2000)
(coxphw
with template = "ARE"
).
If in a fit to coxphw
the betafix
argument was used, then for the
fixed parameters only the point estimates are given.
A matrix with estimates of the generalized concordance probability with accompanying confidence intervalls for each explanatory variable in the model.
Daniela Dunkler
Dunkler D, Schemper M, Heinze G. (2010) Gene Selection in Microarray Survival Studies Under Possibly Non-Proportional Hazards. Bioinformatics 26:784-90.
Xu R and O'Quigley J (2000). Estimating Average Regression Effect Under Non-Proportional Hazards. Biostatistics 1, 423-439.
data("gastric") fit <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR") concord(fit)
data("gastric") fit <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR") concord(fit)
Computes confidence intervals for one or more parameters in a model fitted by
coxphw
. Objects of this class have methods for the functions summary
,
print
, coef
, vcov
, plot
, and confint
.
## S3 method for class 'coxphw' confint(object, parm, level = 0.95, ...)
## S3 method for class 'coxphw' confint(object, parm, level = 0.95, ...)
object |
a fitted model object of class |
parm |
a specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names. If missing, all parameters are considered. |
level |
the confidence level required. |
... |
additional argument(s) for methods. |
A matrix (or vector) with columns giving lower and upper confidence limits for each parameter. These will be labeled as (1-level)/2 and 1 - (1-level)/2 in % (by default 2.5% and 97.5%).
Daniela Dunkler
Weighted Cox regression as proposed by Schemper et al. (2009) doi:10.1002/sim.3623 provides unbiased estimates of average hazard ratios also in case of non-proportional hazards. Time-dependent effects can be conveniently estimated by including interactions of covariates with arbitrary functions of time, with or without making use of the weighting option.
coxphw(formula, data, template = c("AHR", "ARE", "PH"), subset, na.action, robust = TRUE, jack = FALSE, betafix = NULL, alpha = 0.05, trunc.weights = 1, control, caseweights, x = TRUE, y = TRUE, verbose = FALSE, sorted = FALSE, id = NULL, clusterid = NULL, ...)
coxphw(formula, data, template = c("AHR", "ARE", "PH"), subset, na.action, robust = TRUE, jack = FALSE, betafix = NULL, alpha = 0.05, trunc.weights = 1, control, caseweights, x = TRUE, y = TRUE, verbose = FALSE, sorted = FALSE, id = NULL, clusterid = NULL, ...)
formula |
a formula object with the response on the left of the operator and the model terms
on the right. The response must be a survival object as returned by
|
data |
a data frame in which to interpret the variables named in |
template |
choose among three pre-defined templates: |
subset |
expression indicating which subset of the rows of data should be used in the fit. All observations are included by default. |
na.action |
missing-data filtering. Defaults to |
robust |
if set to TRUE, the robust covariance estimate (Lin-Wei) is used; otherwise the Lin-Sasieni covariance estimate is applied. Default is TRUE. |
jack |
if set to TRUE, the variance is based on a complete jackknife. Each individual (as
identified by |
betafix |
can be used to restrict the estimation of one or more regression coefficients to pre-defined
values. A vector with one element for each model term as given in |
alpha |
the significance level (1- |
trunc.weights |
specifies a quantile at which the (combined normalized) weights are to be truncated.
It can be used to increase the precision of the estimates, particularly if
|
control |
Object of class |
caseweights |
vector of case weights, equivalent to |
x |
requests copying explanatory variables into the output object. Default is TRUE. |
y |
requests copying survival information into the output object. Default is TRUE. |
verbose |
requests echoing of intermediate results. Default is FALSE. |
sorted |
if set to TRUE, the data set will not be sorted prior to passing it to FORTRAN. This may speed up computations. Default is FALSE. |
id |
a vector of subject identification integer numbers starting from 1 used only if the data are
in the counting process format. These IDs are used to compute the robust covariance matrix.
If |
clusterid |
a vector of cluster identification integer numbers starting from 1. These IDs are used
to compute the robust covariance matrix. If |
... |
additional arguments. |
If Cox's proportional hazards regression is used in the presence of non-proportional hazards, i.e., with underlying time-dependent hazard ratios of prognostic factors, the average relative risk for such a factor is under- or overestimated and testing power for the corresponding regression parameter is reduced. In such a situation weighted estimation provides a parsimonious alternative to more elaborate modelling of time-dependent effects. Weighted estimation in Cox regression extends the tests by Breslow and Prentice to a multi-covariate situation as does the Cox model to Mantel's logrank test. Weighted Cox regression can also be seen as a robust alternative to the standard Cox estimator, reducing the influence of outlying survival times on parameter estimates.
Three pre-defined templates can be requested:
1) "AHR"
, i.e., estimation of average hazard ratios
(Schemper et al., 2009) using Prentice weights with censoring correction and robust variance estimation;
2) "ARE"
, i.e., estimation of average regression effects (Xu and O'Quigley, 2000) using censoring
correction and robust variance estimation; or
3) "PH"
, i.e., Cox proportional hazards regression
using robust variance estimation.
Breslow's tie-handling method is used by the program, other methods to handle ties are currently not available.
A fit of coxphw
with template = "PH"
will yield identical estimates as a fit of
coxph
using Breslow's tie handling method and robust variance estimation
(using cluster
).
If robust = FALSE
, the program estimates the covariance matrix using the Lin (1991)
and Sasieni (1993) sandwich estimate with
and
denoting the sum
of contributions to the second derivative of the log likelihood, weighted by
and
,
respectively. This estimate is independent from the scaling of the weights and reduces to the inverse
of the information matrix in case of no weighting. However, it is theoretically valid only in case of
proportional hazards. Therefore, since application of weighted Cox regression usually implies a violated
proportional hazards assumption, the robust Lin-Wei covariance estimate is used by default (
robust =
TRUE
).
If some regression coefficients are held constant using betafix
, no standard errors
are given for these coefficients as they are not estimated in the model. The global Wald test
only relates to those variables for which regression coefficients were estimated.
An offset
term can be included in the formula
of coxphw
.
In this way a variable can be specified which is included in the model but its parameter estimate is
fixed at 1.
A list with the following components:
coefficients |
the parameter estimates. |
var |
the estimated covariance matrix. |
df |
the degrees of freedom. |
ci.lower |
the lower confidence limits of exp(beta). |
ci.upper |
the upper confidence limits of exp(beta). |
prob |
the p-values. |
linear.predictors |
the linear predictors. |
n |
the number of observations. |
dfbeta.resid |
matrix of DFBETA residuals. |
iter |
the number of iterations needed to converge. |
method.ties |
the ties handling method. |
PTcoefs |
matrix with scale and shift used for pretransformation of |
cov.j |
the covariance matrix computed by the jackknife method (only computed if |
cov.lw |
the covariance matrix computed by the Lin-Wei method (robust covariance) |
cov.ls |
the covariance matrix computed by the Lin-Sasieni method. |
cov.method |
the method used to compute the (displayed) covariance matrix and the standard errors.
This method is either "jack" if |
w.matrix |
a matrix with four columns according to the number of uncensored failure
times. The first column contains the failure times, the remaining columns (labeled
|
caseweights |
if |
Wald |
Wald-test statistics. |
means |
the means of the covariates. |
offset.values |
offset values. |
dataline |
the first dataline of the input data set (required for |
x |
if |
y |
the response. |
alpha |
the significance level = 1 - confidence level. |
template |
the requested template. |
formula |
the model formula. |
betafix |
the |
call |
the function call. |
The SAS
macro WCM
with similar functionality is offered for download at
https://cemsiis.meduniwien.ac.at/en/kb/science-research/software/statistical-software/wcmcoxphw/ .
Up to Version 2.13 coxphw used a slightly different syntax (arguments: AHR
,
AHR.norobust
, ARE
, PH
, normalize
, censcorr
, prentice
,
breslow
, taroneware
). From Version 3.0.0 on the old syntax is disabled.
From Version 4.0.0 estimation of fractional polynomials is disabled.
Georg Heinze, Meinhard Ploner, Daniela Dunkler
Dunkler D, Ploner M, Schemper M, Heinze G. (2018) Weighted Cox Regression Using the R Package coxphw. JSS 84, 1–26, doi:10.18637/jss.v084.i02.
Lin D and Wei L (1989). The Robust Inference for the Cox Proportional Hazards Model. J AM STAT ASSOC 84, 1074-1078.
Lin D (1991). Goodness-of-Fit Analysis for the Cox Regression Model Based on a Class of Parameter Estimators. J AM STAT ASSOC 86, 725-728.
Royston P and Altman D (1994). Regression Using Fractional Polynomials of Continuous Covariates: Parsimonious Parametric Modelling. J R STAT SOC C-APPL 43, 429-467.
Royston P and Sauerbrei W (2008). Multivariable Model-Building. A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables. Wiley, Chichester, UK.
Sasieni P (1993). Maximum Weighted Partial Likelihood Estimators for the Cox Model. J AM STAT ASSOC 88, 144-152.
Schemper M (1992). Cox Analysis of Survival Data with Non-Proportional Hazard Functions. J R STAT SOC D 41, 455-465.
Schemper M, Wakounig S and Heinze G (2009). The Estimation of Average Hazard Ratios by Weighted Cox Regression. STAT MED 28, 2473-2489. doi:10.1002/sim.3623
Xu R and O'Quigley J (2000). Estimating Average Regression Effect Under Non-Proportional Hazards. Biostatistics 1, 423-439.
concord
, plot.coxphw
, predict.coxphw
, plot.coxphw.predict
, coxph
data("gastric") # weighted estimation of average hazard ratio fit1 <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR") summary(fit1) fit1$cov.lw # robust covariance fit1$cov.ls # Lin-Sasieni covariance # unweighted estimation, include interaction with years # ('radiation' must be included in formula!) gastric$years <- gastric$time / 365.25 fit2 <- coxphw(Surv(years, status) ~ radiation + years : radiation, data = gastric, template = "PH") summary(fit2) # unweighted estimation with a function of time data("gastric") gastric$yrs <- gastric$time / 365.25 fun <- function(t) { (t > 1) * 1 } fit3 <- coxphw(Surv(yrs, status) ~ radiation + fun(yrs):radiation, data = gastric, template = "PH") # for more examples see vignette or predict.coxphw
data("gastric") # weighted estimation of average hazard ratio fit1 <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR") summary(fit1) fit1$cov.lw # robust covariance fit1$cov.ls # Lin-Sasieni covariance # unweighted estimation, include interaction with years # ('radiation' must be included in formula!) gastric$years <- gastric$time / 365.25 fit2 <- coxphw(Surv(years, status) ~ radiation + years : radiation, data = gastric, template = "PH") summary(fit2) # unweighted estimation with a function of time data("gastric") gastric$yrs <- gastric$time / 365.25 fun <- function(t) { (t > 1) * 1 } fit3 <- coxphw(Surv(yrs, status) ~ radiation + fun(yrs):radiation, data = gastric, template = "PH") # for more examples see vignette or predict.coxphw
This is used to set various numeric parameters controlling a Cox model fit using coxphw
.
Typically it would only be used in a call to coxphw
.
coxphw.control(iter.max = 200, maxhs = 5, xconv = 1e-4, gconv = 1e-4, maxstep = 1, round.times.to = 0.00001, add.constant = 0, pc = TRUE, pc.time = TRUE, normalize = TRUE)
coxphw.control(iter.max = 200, maxhs = 5, xconv = 1e-4, gconv = 1e-4, maxstep = 1, round.times.to = 0.00001, add.constant = 0, pc = TRUE, pc.time = TRUE, normalize = TRUE)
iter.max |
maximum number of iterations to attempt for convergence. Default is 200. |
maxhs |
maximum number of step-halvenings per iterations. Default is 5. The increments of the
parameter vector in one Newton-Rhaphson iteration step are halved, unless the new
pseudo-likelihood is greater than the old one, maximally doing |
xconv |
specifies the maximum allowed change in standardized parameter estimates to declare convergence. Default is 0.0001. |
gconv |
specifies the maximum allowed change in score function to declare convergence. Default is 0.0001. |
maxstep |
specifies the maximum change of (standardized) parameter values allowed in one iteration. Default is 1. |
round.times.to |
rounds survival times to the nearest number that is a multiple of
|
add.constant |
this number will be added to all times. It can be used if some survival times are exactly 0. Default is 0. |
pc |
if set to TRUE, it will orthogonalize the model matrix to speed up convergence. Default is TRUE. |
pc.time |
if set to TRUE, it will orthogonalize time-dependent covariates (i.e., interactions of covariates with functions of time) to speed up convergence. Default is TRUE. |
normalize |
if set to TRUE, weights are normalized such that their sum is equal to the number of events. This may speed up or enable convergence, but has no consequences on the estimated regression coefficients. |
A list containing the values of each of the above constants
Daniela Dunkler
Provides fractional polynomials as accessible function.
fp.power(z, a, b = NULL)
fp.power(z, a, b = NULL)
z |
a scalar or vector of positive numerical values. |
a |
first power. |
b |
optional second power. |
The function returns fp(a
) of z
(and optionally fp(b)
of z
).
A matrix with one or two columns (if a second power b
was specified), and number of
rows equal to the length of z
. The columns are sorted by descending power.
Georg Heinze
Royston P and Altman D (1994). Regression Using Fractional Polynomials of Continuous Covariates: Parsimonious Parametric Modelling. J R STAT SOC C-APPL 43, 429-467.
Royston P and Sauerbrei W (2008). Multivariable Model-Building. A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables. Wiley, Chichester, UK.
fp.power(z = c(1, 4, 6), a = 1) fp.power(z = c(1, 4, 6), a = 0.5) fp.power(z = c(1, 4, 6), a = 0.5, b = 0.5) fp.power(z = c(1, 4, 6), a = 0, b = 2)
fp.power(z = c(1, 4, 6), a = 1) fp.power(z = c(1, 4, 6), a = 0.5) fp.power(z = c(1, 4, 6), a = 0.5, b = 0.5) fp.power(z = c(1, 4, 6), a = 0, b = 2)
A data set of survival times of patients with locally advanced, nonresectable gastric carcinoma. The patients were either treated with chemotherapy plus radiation or chemotherapy alone.
data(gastric)
data(gastric)
A data frame with 90 observations on the following 4 variables:
id
unique patient id.
radiation
treatment of either 0 = chemotherapy alone or 1 = chemotherapy plus radiation.
time
survival time in days.
status
0 = censored or 1 = death.
Stablein DM, Carter J, Novak JW. (1981) Analysis of Survival Data with Nonproportional Hazard Functions. Controlled Clinical Trials 2:149-159. OR
https://www.mayo.edu/research/documents/gastrichtml/DOC-10027680
Gastrointestinal Tumor Study Group. (1982) A Comparison of Combination Chemotherapy and Combined Modality Therapy for Locally Advanced Gastric Carcinoma. Cancer 49:1771-7.
data("gastric") plot(survfit(Surv(time, status) ~ radiation, data = gastric), lty = 1:2, las = 1, xscale = 365.25, xlab = "time (years)", ylab = "survival distribution function") coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR")
data("gastric") plot(survfit(Surv(time, status) ~ radiation, data = gastric), lty = 1:2, las = 1, xscale = 365.25, xlab = "time (years)", ylab = "survival distribution function") coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR")
This function plots the weights used in a weighted Cox regression analysis against time.
## S3 method for class 'coxphw' plot(x, rank = FALSE, log = FALSE, legendxy = NULL,...)
## S3 method for class 'coxphw' plot(x, rank = FALSE, log = FALSE, legendxy = NULL,...)
x |
a |
rank |
if set to TRUE, plots the weights against ranked time. Default is FALSE. |
log |
if set to TRUE, shows logarithm of weights. Default is FALSE. |
legendxy |
an optional vector of length 2 of the x and y co-ordinates to be used to position the legend. |
... |
additional arguments for plotting |
The function plots the survival weights, i.e., the left-continuous survivor function estimates, the censoring weights, i.e., estimates of the follow-up distribution obtained by Kaplan-Meier estimation with reversed meaning of the status indicator and the combined normalized weights, i.e. the product of the survival and the censoring weights, rescaled to a mean of 1.
No output value.
In coxphw version 4.0.0 the new plot
function replaces the old plotw
function.
Georg Heinze, Daniela Dunkler
data("gastric") # weighted estimation of average hazard ratio fit1 <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR") plot(fit1) # estimation of average regression effect by inverse probability of censoring weights; # truncate weights at 95th percentile fit2 <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "ARE", trunc.weights = 0.95) plot(fit2)
data("gastric") # weighted estimation of average hazard ratio fit1 <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "AHR") plot(fit1) # estimation of average regression effect by inverse probability of censoring weights; # truncate weights at 95th percentile fit2 <- coxphw(Surv(time, status) ~ radiation, data = gastric, template = "ARE", trunc.weights = 0.95) plot(fit2)
This function visualizes a nonlinear or a time-dependent effect of a predict.coxphw
object.
## S3 method for class 'coxphw.predict' plot(x, addci = TRUE, xlab = NULL, ylab = NULL, ...)
## S3 method for class 'coxphw.predict' plot(x, addci = TRUE, xlab = NULL, ylab = NULL, ...)
x |
an output object of |
addci |
confidence intervalls are obtained. Default is TRUE. |
xlab |
label for x-axis of plot, uses variable specified in |
ylab |
label for y-axis of plot, uses as appropriate either "relative hazard" or "log relative hazard" as default. |
... |
further parameters, to be used for plots (e.g., scaling of axes). |
This function can be used to depict the estimated nonlinear or time-dependent
effect of an object of class predict.coxphw
. It supports simple nonlinear
effects as well as interaction effects of continuous variables with binary
covariates (see examples section in predict.coxphw
. ).
No output value.
In coxphw version 4.0.0 the old plotshape
function is replaced with
predict.coxphw
and plot.coxphw.predict
.
Georg Heinze, Meinhard Ploner, Daniela Dunkler
Royston P and Altman D (1994). Regression Using Fractional Polynomials of Continuous Covariates: Parsimonious Parametric Modelling. Applied Statistics J R STAT SOC C-APPL 43, 429-467.
Royston P and Sauerbrei W (2008). Multivariable Model-building. A pragmatic approach to regression analysis based on fractional polynomials for modelling continuous variables. Wiley, Chichester, UK.
set.seed(30091) n <- 300 x <- 1:n true.func <- function(x) 3 * (x / 100)^{2} - log(x / 100) - 3 * x / 100 x <- round(rnorm(n = x) * 10 + 40, digits = 0) time <- rexp(n = n, rate = 1) / exp(true.func(x)) event <- rep(x = 1, times = n) futime <- runif(n = n, min = 0, max = 309000) event <- (time < futime) + 0 time[event == 0] <- futime[event == 0] my.data <- data.frame(x, time, event) fitahr <- coxphw(Surv(time, event) ~ x, data = my.data, template = "AHR") # estimated function plotx <- quantile(x, probs = 0.05):quantile(x, probs = 0.95) plot(predict(fitahr, type = "shape", newx = plotx, refx = median(x), x = "x", verbose = FALSE)) # true function lines(x = plotx, true.func(plotx) - true.func(median(plotx)), lty = 2) legend("topright", legend=c("AHR estimates", "true"), bty = "n", lty = 1:2, inset = 0.05) # for more examples see predict.coxphw
set.seed(30091) n <- 300 x <- 1:n true.func <- function(x) 3 * (x / 100)^{2} - log(x / 100) - 3 * x / 100 x <- round(rnorm(n = x) * 10 + 40, digits = 0) time <- rexp(n = n, rate = 1) / exp(true.func(x)) event <- rep(x = 1, times = n) futime <- runif(n = n, min = 0, max = 309000) event <- (time < futime) + 0 time[event == 0] <- futime[event == 0] my.data <- data.frame(x, time, event) fitahr <- coxphw(Surv(time, event) ~ x, data = my.data, template = "AHR") # estimated function plotx <- quantile(x, probs = 0.05):quantile(x, probs = 0.95) plot(predict(fitahr, type = "shape", newx = plotx, refx = median(x), x = "x", verbose = FALSE)) # true function lines(x = plotx, true.func(plotx) - true.func(median(plotx)), lty = 2) legend("topright", legend=c("AHR estimates", "true"), bty = "n", lty = 1:2, inset = 0.05) # for more examples see predict.coxphw
This function obtains the effect estimates (e.g. of a nonlinear or a
time-dependent effect) at specified values of a continuous
covariable for a model fitted by coxphw
. It prints the
relative or log relative hazard. Additionally, the linear predictor lp
or the risk score exp(lp) can be obtained.
## S3 method for class 'coxphw' predict(object, type = c("shape", "slice.time", "slice.z", "slice.x", "lp", "risk"), x = NULL, newx = NA, refx = NA, z = NULL, at = NULL, exp = FALSE, se.fit = FALSE, pval = FALSE, digits = 4, verbose = FALSE, ...)
## S3 method for class 'coxphw' predict(object, type = c("shape", "slice.time", "slice.z", "slice.x", "lp", "risk"), x = NULL, newx = NA, refx = NA, z = NULL, at = NULL, exp = FALSE, se.fit = FALSE, pval = FALSE, digits = 4, verbose = FALSE, ...)
object |
an output object of |
type |
the type of predicted value. Choices are:
|
x |
name of the continuous or time variable (use "") for |
newx |
the data values for |
refx |
the reference value for variable |
z |
variable which is in interaction with |
at |
if |
exp |
if set to TRUE (default), the log relative hazard is given, otherwise the relative hazard
is requested for |
se.fit |
if set to TRUE, pointwise standard errors are produced for the predictions for
|
pval |
if set to TRUE add Wald-test p-values to effect estimates at values of
|
digits |
number of printed digits. Default is 4. |
verbose |
if set to TRUE (default), results are printed. |
... |
further parameters. |
This function can be used to depict the estimated nonlinear or time-dependent
effect of an object of class coxphw
. It supports simple nonlinear
effects as well as interaction effects of a continuous variable with a binary
covariate or with time (see examples section).
If the effect estimates of a simple nonlinear effect of x
without
interaction is requested with type = "shape"
, then x
(the usually
continuous covariate), refx
(the reference value of x
) and newx
(for these values of x
the effect estimates are obtained) must be given.
If the effect estimates of an interaction of z
with x
are requested
with type = "slice.x"
, then x
(the usually continuous variable),
z
(the categorical variable) and newx
(for these values of x
the effect estimates are obtained) must be given.
If the effect estimates of an interaction of z
with x
for one level of z
are requested with type = "slice.z"
), then x
(the usually continuous variable),
z
(the categorical variable), at
(at which level of z
),
refx
(the reference value of x
), and newx
(for these values of x
the effect estimates are obtained) must be given.
If the effect estimates of an interaction of z
with time
are requested
with type = "slice.time"
, then x
(the time
), z
(the
categorical variable) and newx
(for these values of x
the effect
estimates are obtained) must be given.
Note that if the model formula contains time-by-covariate interactions, then the linear predictor and the risk score are obtained for the failure or censoring time of each subject.
If type = "shape"
, "slice.time"
, "slice.x"
, or "slice.z"
a list with the following components:
estimates |
a matrix with estimates of (log) relative hazard and corresponding confidence limits. |
se |
pointwise standard errors, only if |
p |
p-value, only if |
alpha |
the significance level = 1 - confidence level. |
exp |
an indicator if log relative hazard or relative hazard was obtained. |
x |
name of |
If type = "lp"
or "risk"
, a vector.
In coxphw version 4.0.0 the old plotshape
function is replaced with
predict.coxphw
and plot.coxphw.predict
.
Georg Heinze, Meinhard Ploner, Daniela Dunkler
Royston P and Altman D (1994). Regression Using Fractional Polynomials of Continuous Covariates: Parsimonious Parametric Modelling. Applied Statistics 43, 429-467.
Royston P and Sauerbrei W (2008). Multivariable Model-building. A pragmatic approach to regression analysis based on fractional polynomials for modelling continuous variables. Wiley, Chichester, UK.
### Example for type = "slice.time" data("gastric") gastric$yrs <- gastric$time / 365.25 # check proportional hazards fitcox <- coxph(Surv(yrs, status) ~ radiation + cluster(id), data = gastric, x = TRUE, method = "breslow") fitcox.ph <- cox.zph(fit = fitcox, transform = "identity") ## compare and visualize linear and log-linear time-dependent effects of radiation fit1 <- coxphw(Surv(yrs, status) ~ yrs * radiation, data = gastric, template = "PH") summary(fit1) predict(fit1, type = "slice.time", x = "yrs", z = "radiation", newx = c(0.5, 1, 2), verbose = TRUE, exp = TRUE, pval = TRUE) fit2 <- coxphw(Surv(yrs, status) ~ log(yrs) * radiation, data = gastric, template = "PH") summary(fit2) predict(fit2, type = "slice.time", x = "yrs", z = "radiation", newx = c(0.5, 1, 2), verbose = TRUE, exp = TRUE, pval = TRUE) plotx <- seq(from = quantile(gastric$yrs, probs = 0.05), to = quantile(gastric$yrs, probs = 0.95), length = 100) y1 <- predict(fit1, type = "slice.time", x = "yrs", z = "radiation", newx = plotx) y2 <- predict(fit2, type = "slice.time", x = "yrs", z = "radiation", newx = plotx) plot(x = fitcox.ph, se = FALSE, xlim = c(0, 3), las = 1, lty = 3) abline(a = 0, b = 0, lty = 3) lines(x = plotx, y = y1$estimates[, "coef"], col = "red", lty = 1, lwd = 2) lines(x = plotx, y = y2$estimates[, "coef"], col = "blue", lty = 2, lwd = 2) legend(x = 1.7, y = 1.6, title = "time-dependent effect", title.col = "black", legend = c("LOWESS", "linear", "log-linear"), col = c("black", "red", "blue"), lty = c(3, 1:2), bty = "n", lwd = 2, text.col = c("black", "red", "blue")) ### Example for type = "shape" set.seed(512364) n <- 200 x <- 1:n true.func <- function(x) 2.5 * log(x) - 2 x <- round(runif(x) * 60 + 10, digits = 0) time <- round(100000 * rexp(n= n, rate = 1) / exp(true.func(x)), digits = 1) event <- rep(x = 1, times = n) my.data <- data.frame(x,time,event) fit <- coxphw(Surv(time, event) ~ log(x) + x, data = my.data, template = "AHR") predict(fit, type = "shape", newx = c(30, 50), refx = 40, x = "x", verbose = TRUE) plotx <- seq(from = quantile(x, probs = 0.05), to = quantile(x, probs = 0.95), length = 100) plot(predict(fit, type = "shape", newx = plotx, refx = 40, x = "x")) ### Example for type = "slice.x" and "slice.z" set.seed(75315) n <- 200 trt <- rbinom(n = n, size = 1, prob = 0.5) x <- 1:n true.func <- function(x) 2.5 * log(x) - 2 x <- round(runif(n = x) * 60 + 10, digits = 0) time <- 100 * rexp(n = n, rate = 1) / exp(true.func(x) / 4 * trt - (true.func(x) / 4)^2 * (trt==0)) event <- rep(x = 1, times = n) my.data <- data.frame(x, trt, time, event) fun<-function(x) x^(-2) fit <- coxphw(Surv(time, event) ~ x * trt + fun(x) * trt , data = my.data, template = "AHR", verbose = FALSE) # plots the interaction of trt with x (the effect of trt dependent on the values of x) plotx <- quantile(x, probs = 0.05):quantile(x, probs = 0.95) plot(predict(fit, type = "slice.x", x = "x", z = "trt", newx = plotx, verbose = FALSE), main = "interaction of trt with x") # plot the effect of x in subjects with trt = 0 y0 <- predict(fit, type = "slice.z", x = "x", z = "trt", at = 0, newx = plotx, refx = median(x), verbose = FALSE) plot(y0, main = "effect of x in subjects with trt = 0") # plot the effect of x in subjects with trt = 1 y1 <- predict(fit, type = "slice.z", x = "x", z = "trt", at = 1, newx = plotx, refx = median(x), verbose = FALSE) plot(y1, main = "effect of x in subjects with trt = 1") # Example for type = "slice.time" set.seed(23917) time <- 100 * rexp(n = n, rate = 1) / exp((true.func(x) / 10)^2 / 2000 * trt + trt) event <- rep(x = 1, times = n) my.data <- data.frame(x, trt, time, event) plot.x <- seq(from = 1, to = 100, by = 1) fun <- function(t) { PT(t)^-2 * log(PT(t)) } fun2 <- function(t) { PT(t)^-2 } fitahr <- coxphw(Surv(time, event) ~ fun(time) * trt + fun2(time) * trt + x, data = my.data, template = "AHR") yahr <- predict(fitahr, type = "slice.time", x = "time", z = "trt", newx = plot.x) fitph <- coxphw(Surv(time, event) ~ fun(time) * trt + fun2(time) * trt + x, data = my.data, template = "PH") yph <- predict(fitph, type = "slice.time", x = "time", z = "trt", newx = plot.x) plot(yahr, addci = FALSE) lines(yph$estimates$time, yph$estimates$coef, lty = 2) legend("bottomright", legend = c("AHR", "PH"), bty = "n", lty = 1:2, inset = 0.05)
### Example for type = "slice.time" data("gastric") gastric$yrs <- gastric$time / 365.25 # check proportional hazards fitcox <- coxph(Surv(yrs, status) ~ radiation + cluster(id), data = gastric, x = TRUE, method = "breslow") fitcox.ph <- cox.zph(fit = fitcox, transform = "identity") ## compare and visualize linear and log-linear time-dependent effects of radiation fit1 <- coxphw(Surv(yrs, status) ~ yrs * radiation, data = gastric, template = "PH") summary(fit1) predict(fit1, type = "slice.time", x = "yrs", z = "radiation", newx = c(0.5, 1, 2), verbose = TRUE, exp = TRUE, pval = TRUE) fit2 <- coxphw(Surv(yrs, status) ~ log(yrs) * radiation, data = gastric, template = "PH") summary(fit2) predict(fit2, type = "slice.time", x = "yrs", z = "radiation", newx = c(0.5, 1, 2), verbose = TRUE, exp = TRUE, pval = TRUE) plotx <- seq(from = quantile(gastric$yrs, probs = 0.05), to = quantile(gastric$yrs, probs = 0.95), length = 100) y1 <- predict(fit1, type = "slice.time", x = "yrs", z = "radiation", newx = plotx) y2 <- predict(fit2, type = "slice.time", x = "yrs", z = "radiation", newx = plotx) plot(x = fitcox.ph, se = FALSE, xlim = c(0, 3), las = 1, lty = 3) abline(a = 0, b = 0, lty = 3) lines(x = plotx, y = y1$estimates[, "coef"], col = "red", lty = 1, lwd = 2) lines(x = plotx, y = y2$estimates[, "coef"], col = "blue", lty = 2, lwd = 2) legend(x = 1.7, y = 1.6, title = "time-dependent effect", title.col = "black", legend = c("LOWESS", "linear", "log-linear"), col = c("black", "red", "blue"), lty = c(3, 1:2), bty = "n", lwd = 2, text.col = c("black", "red", "blue")) ### Example for type = "shape" set.seed(512364) n <- 200 x <- 1:n true.func <- function(x) 2.5 * log(x) - 2 x <- round(runif(x) * 60 + 10, digits = 0) time <- round(100000 * rexp(n= n, rate = 1) / exp(true.func(x)), digits = 1) event <- rep(x = 1, times = n) my.data <- data.frame(x,time,event) fit <- coxphw(Surv(time, event) ~ log(x) + x, data = my.data, template = "AHR") predict(fit, type = "shape", newx = c(30, 50), refx = 40, x = "x", verbose = TRUE) plotx <- seq(from = quantile(x, probs = 0.05), to = quantile(x, probs = 0.95), length = 100) plot(predict(fit, type = "shape", newx = plotx, refx = 40, x = "x")) ### Example for type = "slice.x" and "slice.z" set.seed(75315) n <- 200 trt <- rbinom(n = n, size = 1, prob = 0.5) x <- 1:n true.func <- function(x) 2.5 * log(x) - 2 x <- round(runif(n = x) * 60 + 10, digits = 0) time <- 100 * rexp(n = n, rate = 1) / exp(true.func(x) / 4 * trt - (true.func(x) / 4)^2 * (trt==0)) event <- rep(x = 1, times = n) my.data <- data.frame(x, trt, time, event) fun<-function(x) x^(-2) fit <- coxphw(Surv(time, event) ~ x * trt + fun(x) * trt , data = my.data, template = "AHR", verbose = FALSE) # plots the interaction of trt with x (the effect of trt dependent on the values of x) plotx <- quantile(x, probs = 0.05):quantile(x, probs = 0.95) plot(predict(fit, type = "slice.x", x = "x", z = "trt", newx = plotx, verbose = FALSE), main = "interaction of trt with x") # plot the effect of x in subjects with trt = 0 y0 <- predict(fit, type = "slice.z", x = "x", z = "trt", at = 0, newx = plotx, refx = median(x), verbose = FALSE) plot(y0, main = "effect of x in subjects with trt = 0") # plot the effect of x in subjects with trt = 1 y1 <- predict(fit, type = "slice.z", x = "x", z = "trt", at = 1, newx = plotx, refx = median(x), verbose = FALSE) plot(y1, main = "effect of x in subjects with trt = 1") # Example for type = "slice.time" set.seed(23917) time <- 100 * rexp(n = n, rate = 1) / exp((true.func(x) / 10)^2 / 2000 * trt + trt) event <- rep(x = 1, times = n) my.data <- data.frame(x, trt, time, event) plot.x <- seq(from = 1, to = 100, by = 1) fun <- function(t) { PT(t)^-2 * log(PT(t)) } fun2 <- function(t) { PT(t)^-2 } fitahr <- coxphw(Surv(time, event) ~ fun(time) * trt + fun2(time) * trt + x, data = my.data, template = "AHR") yahr <- predict(fitahr, type = "slice.time", x = "time", z = "trt", newx = plot.x) fitph <- coxphw(Surv(time, event) ~ fun(time) * trt + fun2(time) * trt + x, data = my.data, template = "PH") yph <- predict(fitph, type = "slice.time", x = "time", z = "trt", newx = plot.x) plot(yahr, addci = FALSE) lines(yph$estimates$time, yph$estimates$coef, lty = 2) legend("bottomright", legend = c("AHR", "PH"), bty = "n", lty = 1:2, inset = 0.05)
coxphw
This class of objects is returned by the coxphw
function.
Objects of this class have methods for the functions summary
, print
, coef
,
vcov
, plot
, and confint
.
## S3 method for class 'coxphw' print(x, ...)
## S3 method for class 'coxphw' print(x, ...)
x |
object of class |
... |
further arguments. |
If some regression coefficients were held fixed by betafix
, then no standard errors are
given for these coefficients as they are not estimated in the model. The global Wald test
only relates to those variables for which regression coefficients were estimated.
Georg Heinze, Daniela Dunkler
predict.coxphw
This class of objects is returned by the predict.coxphw
function.
Objects of this class have methods for the functions print
and plot
.
## S3 method for class 'coxphw.predict' print(x, ...)
## S3 method for class 'coxphw.predict' print(x, ...)
x |
object of class |
... |
further arguments. |
Daniela Dunkler
coxphw
, predict.coxphw
plot.coxphw.predict
Provides automatic pretransformation of variables (to well-scaled and nonzero values).
PT(z)
PT(z)
z |
a vector of numerical values. |
The function transforms a variable by shifting to positive values, and dividing by scaling factor (a power of 10) such that the standard deviation is approximately equal to 1.
(z
+ shift) / scale
Georg Heinze
PT(z = c(-6, -1, 4, 6))
PT(z = c(-6, -1, 4, 6))
coxphw
This class of objects is returned by the coxphw
function.
Objects of this class have methods for the functions summary
, print
, coef
,
vcov
, plot
, and confint
.
## S3 method for class 'coxphw' summary(object, print = TRUE, ...)
## S3 method for class 'coxphw' summary(object, print = TRUE, ...)
object |
object of class |
print |
print summary. Default is TRUE. |
... |
further arguments. |
If some regression coefficients are held fixed by betafix
, no standard errors and related quantities
are given for these coefficients as they are not estimated in the model. The global Wald test
only relates to those variables for which regression coefficients were estimated.
Georg Heinze, Daniela Dunkler
coxphw
This class of objects is returned by the coxphw
function.
Objects of this class have methods for the functions summary
, print
, coef
,
vcov
, plot
, and confint
.
## S3 method for class 'coxphw' vcov(object, ...)
## S3 method for class 'coxphw' vcov(object, ...)
object |
object of class |
... |
further arguments. |
Daniela Dunkler
Obtain Wald chi-squared test statistics and p-values for one or more regression coefficients given their variance-covariance matrix.
wald(coeff, cov, index = NULL, h0 = NULL)
wald(coeff, cov, index = NULL, h0 = NULL)
coeff |
a vector with regression coefficients. |
cov |
a variance-covariance matrix. |
index |
an integer specifying which parameters should be jointly tested. Default is to test
all parameters given in |
h0 |
a vector with the same length as |
The test is based on the assumption that the coefficients asymptotically follow
a (multivariate) normal distribution with mean coeff
and a variance-covariance
matrix cov
.
A vector with the following components:
chi2 |
the Wald-test statistic. |
df |
degrees of freedom. |
p |
p-value. |
Daniela Dunkler