knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This vignette demonstrates the use of the Wcompo
package
in the proportional means regression of weighted composite endpoint
of recurrent hospitalizations and death (Mao and Lin, 2016, Biostatistics).
Let $D$ denote the survival time and write $N_D(t)=I(D\leq t)$. Let $N_1(t),\ldots, N_{K-1}(t)$ denote the counting processes for $K-1$ different types of possibly recurrent nonfatal event. We are interested in a weighted composite event process of the form \begin{equation} \mathcal R(t)=\sum_{k=1}^{K-1} w_kN_k(t)+w_DN_D(t), \end{equation} where $w_1,\ldots, w_K$ and $w_D$ are prespecified weights. To reflect the greater importance of death, we typically set $w_D>w_k$ $(k=1,\ldots, K-1)$.
Let $\boldsymbol Z$ denote a vector of baseline covariates. We model the conditional mean of $\mathcal R(t)$ given $\boldsymbol Z$ by \begin{equation}\tag{1} E{\mathcal R(t)\mid\boldsymbol Z}=\exp(\boldsymbol\beta^{\rm T}\boldsymbol Z)\mu_0(t), \end{equation} where $\boldsymbol\beta$ is a vector of regression coefficients and $\mu_0(t)$ is a nonparametric baseline mean function. We call model (1) the proportional means (PM) model because it implies that the conditional mean ratio between any two covariate groups is proportional over time, i.e., [\frac{E{\mathcal R(t)\mid\boldsymbol Z_i}}{E{\mathcal R(t)\mid\boldsymbol Z_j}} =\exp{\beta^{\rm T}(\boldsymbol Z_i-\boldsymbol Z_j)}.] This also means that the components of $\boldsymbol\beta$ can be interpreted as the log-mean ratios associated with unit increases in the corresponding covariates. With censored data, $\boldsymbol\beta$ can be estimated using the inverse probability censoring weighting (IPCW) technique.
The basic function to fit the PM model is CompoML()
.
To use the function, the input data must be
in the "long" format. Specifically, we need an id
variable containing
the unique patient identifiers, a time
variable containing the event times,
a status
variable labeling the event type (1
for death, 2, ..., K
for nonfatal event types $1,\ldots, K-1$, and 0
for censoring),
and a covariate matrix Z
. To fit an unweighted PM model (i.e., $w_D=w_1=\cdots=w_{K-1}=1$),
run the function in the default mode
obj <- CompoML(id,time,status,Z)
To weight the components differently, use an optional argument w
to specify the $K$-vector of weights $(w_D,w_1,\ldots, w_{K-1})^{\rm T}$.
Among the output of the CompoML()
function are beta
, the estimated $\boldsymbol\beta$,
and var
, its estimated covariance matrix. For a summary of analysis results, simply print
the returned object. To predict the model-based mean function for a new covariate z
, use
plot(obj, z)
Heart Failure: A Controlled Trial Investigating Outcomes of Exercise Training (HF-ACTION) was a randomized controlled clinical trial to evaluate the efficacy and safety of exercise training among patients with heart failure (O’Connor and others, 2009). A total of 2331 medically stable outpatients with heart failure and reduced ejection fraction were recruited between April 2003 and February 2007 at 82 centers in the USA, Canada, and France. Patients were randomly assigned to usual care alone or usual care plus aerobic exercise training that consists of 36 supervised sessions followed by home-based training. The usual care group consisted of 1172 patients (follow-up data not available for 1 patient), and the exercise training group consisted of 1159 patients. There were a large number of hospital admissions (nonfatal event) and a considerable number of deaths in each treatment arm.
To illustrate the PM model, we analyze a mock dataset hfmock
consisting of 963 patients,
461 in exercise training and the remaining 502 in usual care.
This dataset is contained in the Wcompo
package.
## load the package and data library(Wcompo) data(hfmock) head(hfmock)
The variables id
, time
(in units of months), and status
(2: hospitalization, 1: death)
are already in the required forms for CompoML()
. The binary variables Training
and HF.etiology
are indicators for exercise training vs usual care and for ischemic vs non-ischemic patients, respectively.
We fit a PM model to the composite of recurrent hospitalizations and death weighted by $1:2$
against the treatment and etiology indicators:
## fit a weighted PM (w_D=2, w_1=1) obj <- CompoML(hfmock$id,hfmock$time,hfmock$status,hfmock[,c("Training","HF.etiology")], w=c(2,1)) ## print out the result obj
The above output shows that adding exercise training reduces the average number of composite events, with death counted twice as much as hospitalization, by 1-0.592=40.8\% compared to usual care lone. This effect is highly significant ($p$-value $1.357e-08$). To display the differences visually, we plot the model-based mean functions by treatment for each etiology.
oldpar <- par(mfrow = par("mfrow")) par(mfrow=c(1,2)) ## plot the estimated mean function for ## non-ischemic patients by treatment plot(obj,c(1,0),ylim=c(0,1.5),xlim=c(0,50), main="Non-ischemic", xlab="Time (months)",cex.main=1.2,lwd=2) plot(obj,c(0,0),add=TRUE,cex.main=1.2,lwd=2,lty=2) legend("topleft",lty=1:2,lwd=2,c("Exercise training","Usual care")) ## plot the estimated mean function for ## ischemic patients by treatment plot(obj,c(1,1),ylim=c(0,1.5),xlim=c(0,50), main="Ischemic", xlab="Time (months)",cex.main=1.2,lwd=2) plot(obj,c(0,1),add=TRUE,cex.main=1.2,lwd=2,lty=2) legend("topleft",lty=1:2,lwd=2,c("Exercise training","Usual care")) par(oldpar)
We can see that the mean function for exercise training is substantially lower than that for usual care in both ischemic and non-ischemic patients.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.