The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
In this vignette, we describe the plot method for
objects of class singleEventCB which is obtained from
running the fitSmoothHazard function. There are currently
two types of plots: hazard functions and hazard ratios. We describe each
one in detail below. Note that the plot method has only
been properly tested for family="glm".
The hazard function plots require the visreg
package.
To illustrate hazard function plots, we will use the breast cancer
dataset which contains the observations of 686 women taken from the TH.data
package. This dataset is also available from the casebase
package. In the following, we will show different hazard functions for
different combinations of continuous, binary variables as well as their
interactions.
We first fit a main effects only model with a spline on
log(time) and hormonal therapy as main effects.
mod_cb <- fitSmoothHazard(cens ~ ns(log(time), df = 3) + hormon,
data = brcancer,
time = "time")
summary(mod_cb)All arguments needed for the hazard function plots are supplied
through the hazard.params argument. This is a named list of
arguments which will override the defaults passed to
visreg::visreg(). The default arguments are
list(fit = x, trans = exp, plot = TRUE, rug = FALSE, alpha = 1, partial = FALSE, overlay = TRUE).
For example, if you want a 95% confidence band, specify
hazard.params = list(alpha = 0.05). For a complete list of
options, please see the visreg
vignettes.
We first plot the hazard as a function of time, for
hormon = 0 and hormon = 1. This is achieved by
specifying the xvar argument, as well as the
cond argument. The cond argument must be
provided as a named list. Each element of that list specifies the value
for one of the terms in the model; any elements left unspecified are
filled in with the median/most common category. Note that even though we
fit the log(time), we must specify time in the
xvar argument.
Alternatively, we can plot the hazard functions on the same plot.
This is accomplished with the by argument:
Note that if we want to extract the data used to construct the plot,
e.g. to create our own, we simply assign the call to plot
to an object (we may optionally set plot=FALSE in the
hazard.params argument as to not print any plots):
The function is flexible because you may leverage
ggplot2 just by specifying gg = TRUE, the plot
will return a ggplot object:
gg_object <- plot(mod_cb,
hazard.params = list(xvar = "time",
by = "hormon",
alpha = 0.20, # 80% CI
ylab = "Hazard",
gg = TRUE)) Now we can use it downstream for any plot while leveraging the entire
ggplot2 ecosystem of packages and functions:
Next, we fit an interaction model with a time-varying covariate, i.e. to test the hypothesis that the effect of hormonal therapy on the hazard varies with time.
mod_cb_tvc <- fitSmoothHazard(cens ~ hormon * ns(log(time), df = 3),
data = brcancer,
time = "time")
summary(mod_cb_tvc)Now we can easily plot the hazard function over time for each
hormon group:
Now we fit a model with an interaction between a continuous variable,
estrogen receptor (in fmol), and time.
mod_cb_tvc <- fitSmoothHazard(cens ~ estrec * ns(log(time), df = 3),
data = brcancer,
time = "time")
summary(mod_cb_tvc)There are now many ways to plot the time-varying effect of estrogen
receptor on the hazard function. The default is to plot the 10th, 50th
and 90th quantiles of the by variable:
# computed at the 10th, 50th and 90th quantiles of estrec
plot(mod_cb_tvc,
hazard.params = list(xvar = "time",
by = "estrec",
alpha = 1,
ylab = "Hazard")) We can also show the quartiles of estrec by specifying
the breaks argument. If breaks is a single
number, that will be the used as the number of breaks:
# computed at quartiles of estrec
plot(mod_cb_tvc,
hazard.params = list(xvar = c("time"),
by = "estrec",
alpha = 1,
breaks = 4,
ylab = "Hazard")) Alternatively, if breaks is a vector, it will be used as
the actual values to be used:
# computed where I want
plot(mod_cb_tvc,
hazard.params = list(xvar = c("time"),
by = "estrec",
alpha = 1,
breaks = c(3,2200),
ylab = "Hazard")) visreg2d(mod_cb_tvc,
xvar = "time",
yvar = "estrec",
trans = exp,
print.cond = TRUE,
zlab = "Hazard",
plot.type = "image")
visreg2d(mod_cb_tvc,
xvar = "time",
yvar = "estrec",
trans = exp,
print.cond = TRUE,
zlab = "Hazard",
plot.type = "persp")
# this can also work if 'rgl' is installed
# visreg2d(mod_cb_tvc,
# xvar = "time",
# yvar = "estrec",
# trans = exp,
# print.cond = TRUE,
# zlab = "Hazard",
# plot.type = "rgl")All the examples so far have only included two predictors in the regression equation. In this example, we fit a smooth hazard model with several predictors:
mod_cb_tvc <- fitSmoothHazard(cens ~ estrec * ns(log(time), df = 3) +
horTh +
age +
menostat +
tsize +
tgrade +
pnodes +
progrec,
data = brcancer,
time = "time")
summary(mod_cb_tvc)In the following plot, we show the time-varying effect of
estrec while controlling for all other variables. By
default, the other terms in the model are set to their median if the
term is numeric or the most common category if the term is a factor. The
values of the other variables are shown in the output:
plot(mod_cb_tvc,
hazard.params = list(xvar = "time",
by = "estrec",
alpha = 1,
breaks = 2,
ylab = "Hazard"))You can of course set the values of the other covariates as before,
i.e. by specifying the cond argument as a named list to the
hazard.params argument:
In this section we illustrate how to plot hazard ratios using the
plot method for objects of class singleEventCB
which is obtained from running the fitSmoothHazard
function. Note that these function have only been thoroughly tested with
family = "glm".
In what follows, the hazard ratio for a variable \(X\) is defined as
\[ \frac{h\left(t | X=x_1, \mathbf{Z}=\mathbf{z_1} ; \hat{\beta}\right)}{h(t | X=x_0, \mathbf{Z}=\mathbf{z_0} ; \hat{\beta})} \] where \(h(t|\cdot;\hat{\beta})\) is the hazard rate as a function of the variable \(t\) (which is usually time, but can be any other continuous variable), \(x_1\) is the value of \(X\) for the exposed group, \(x_0\) is the value of \(X\) for the unexposed group, \(\mathbf{Z}\) are other covariates in the model which are equal to \(\mathbf{z_1}\) in the exposed and \(\mathbf{z_0}\) in the unexposed group, and \(\hat{\beta}\) are the estimated regression coefficients.
As indicated by the formula above, it is most instructive to plot the hazard ratio as a function of a variable \(t\) only if there is an interaction between \(t\) and \(X\). Otherwise, the resulting plot will simply be a horizontal line across time.
We use data from the Manson trial (NEJM 2003) which is included in
the casebase package. This randomized clinical trial
investigated the effect of estrogen plus progestin (estPro)
on coronary heart disease (CHD) risk in 16,608 postmenopausal women who
were 50 to 79 years of age at base line. Participants were randomly
assigned to receive estPro or placebo. The
primary efficacy outcome of the trial was CHD (nonfatal myocardial
infarction or death due to CHD).
We fit a model with the interaction between time and treatment arm. We are therefore interested in visualizing the hazard ratio of the treatment over time.
data("eprchd")
eprchd <- transform(eprchd,
treatment = factor(treatment, levels = c("placebo","estPro")))
str(eprchd)
fit_mason <- fitSmoothHazard(status ~ treatment*time,
data = eprchd,
time = "time")
summary(fit_mason)To plot the hazard ratio, we must specify the newdata
argument with a covariate pattern for the reference group. In this
example, we treat the placebo as the reference group.
Because we have fit an interaction with time, we also provide a sequence
of times at which we would like to calculate the hazard ratio.
newtime <- quantile(fit_mason[["originalData"]][[fit_mason[["timeVar"]]]],
probs = seq(0.01, 0.99, 0.01))
# reference category
newdata <- data.frame(treatment = factor("placebo",
levels = c("placebo", "estPro")),
time = newtime)
str(newdata)
plot(fit_mason,
type = "hr",
newdata = newdata,
var = "treatment",
increment = 1,
xvar = "time",
ci = T,
rug = T)In the call to plot we specify the xvar
which is the variable plotted on the x-axis, the var
argument which specified the variable for which we want the hazard
ratio. The increment = 1 indicates that we want to
increment var by 1 level, which in this case is
estPro. Alternatively, we can specify the
exposed argument which should be a function that takes
newdata and returns the exposed dataset. The following call
is equivalent to the one above:
plot(fit_mason,
type = "hr",
newdata = newdata,
exposed = function(data) transform(data, treatment = "estPro"),
xvar = "time",
ci = T,
rug = T)Alternatively, if we want the placebo group to be the
exposed group, we can change the newdata argument to the
following:
newdata <- data.frame(treatment = factor("estPro",
levels = c("placebo", "estPro")),
time = newtime)
str(newdata)
levels(newdata$treatment)Note that the reference category in newdata is still
placebo. Therefore we must set increment = -1
in order to get the exposed dataset:
plot(fit_mason,
type = "hr",
newdata = newdata,
var = "treatment",
increment = -1,
xvar = "time",
ci = TRUE,
rug = TRUE)If the \(X\) variable has more than
two levels, than, increment works the same way,
e.g. increment = 2 will provide an exposed
group two levels above the value in newdata.
In order to save the data used to make the plot, you simply have to
assign the call to plot to a variable. This is particularly
useful if you want to really customize the plot aesthetics:
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.