sahirbhatnagar / casebase
1
#' Population Time Plot
2
#'
3
#' @description \code{plot} method for objects of class \code{popTime}
4
#'
5
#' @param x an object of class \code{popTime} or \code{popTimeExposure}.
6
#' @param ... Ignored.
7
#' @param xlab,ylab The title of the respective axis. Default: 'Follow-up time'
8
#'   for xlab and 'Population' for ylab
9
#' @param add.case.series Logical indicating if the case series should be added
10
#'   to the plot. Default: TRUE
11
#' @param add.base.series Logical indicating if the base series should be added
12
#'   to the plot. Default: FALSE
13
#' @param add.competing.event Logical indicating if the competing event should
14
#'   be added to the plot. Default: FALSE
15
#' @param casebase.theme Logical indication if the casebase theme be used. The
16
#'   casebase theme uses \code{\link{theme_minimal}}. Default: TRUE.
17
#' @param ribbon.params A list containing arguments that are passed to
18
#'   \code{\link{geom_ribbon}} which is used to plot the
19
#'   population-time area. These arguments will override the function defaults.
20
#'   For example, you can set \code{ribbon.params = list(colour = 'green')} if
21
#'   you want the area to be green.
22
#' @param case.params,base.params,competing.params A list containing arguments
23
#'   that are passed to \code{\link{geom_point}} which is used to plot
24
#'   the case series, base series, competing events. These arguments will
25
#'   override the function defaults. For example, you can set \code{case.params
26
#'   = list(size = 1.5)} if you want to increase the point size for the case
27
#'   series points. Note: do not use this argument to change the color of the
28
#'   points. Doing so will result in unexpected results for the legend. See the
29
#'   \code{color.params} and \code{fill.params} arguments, if you want to change
30
#'   the color of the points.
31
#' @param color.params A list containing arguments that are passed to
32
#'   \code{\link{scale_color_manual}} which is used to plot the legend.
33
#'   Only used if \code{legend=TRUE}. These arguments will override the function
34
#'   defaults. Use this argument if you want to change the color of the points.
35
#'   See examples for more details.
36
#' @param fill.params A list containing arguments that are passed to
37
#'   \code{\link{scale_fill_manual}} which is used to plot the legend.
38
#'   Only used if \code{legend=TRUE}. These arguments will override the function
39
#'   defaults. Use this argument if you want to change the color of the points.
40
#'   See examples for more details.
41
#' @param theme.params A list containing arguments that are passed to
42
#'   \code{\link{theme}}. For example \code{theme.params =
43
#'   list(legend.position = 'none')}.
44
#' @param facet.params A list containing arguments that are passed to
45
#'   \code{\link{facet_wrap}} which is used to create facet plots. Only
46
#'   used if plotting exposure stratified population time plots. These arguments
47
#'   will override the function defaults.
48
#' @param ratio If \code{add.base.series=TRUE}, integer, giving the ratio of the
49
#'   size of the base series to that of the case series. This argument is passed
50
#'   to the \code{\link{sampleCaseBase}} function. Default: 10.
51
#' @param censored.indicator If \code{add.base.series=TRUE}, a character string
52
#'   of length 1 indicating which value in event is the censored. This function
53
#'   will use relevel to set \code{censored.indicator} as the reference level.
54
#'   This argument is ignored if the event variable is a numeric. This argument
55
#'   is passed to the \code{\link{sampleCaseBase}} function.
56
#' @param comprisk If \code{add.base.series=TRUE}, logical indicating whether we
57
#'   have multiple event types and that we want to consider some of them as
58
#'   competing risks. This argument is passed to the
59
#'   \code{\link{sampleCaseBase}} function. Note: should be \code{TRUE} if your
60
#'   data has competing risks, even if you don't want to add competing risk
61
#'   points (\code{add.competing.event=FALSE}). Default: FALSE
62
#' @param legend Logical indicating if a legend should be added to the plot.
63
#'   Note that if you want to change the colors of the points, through the
64
#'   \code{color.params} and \code{fill.params} arguments, then set
65
#'   \code{legend=TRUE}. If you want to change the color of the points but not
66
#'   have a legend, then set \code{legend=TRUE} and \code{theme.params =
67
#'   list(legend.position = 'none'}. Default: FALSE
68
#' @param legend.position Deprecated. Specify the legend.position argument
69
#'   instead in the \code{theme.params} argument. e.g. \code{theme.params =
70
#'   list(legend.position = 'bottom')}.
71
#' @param line.width Deprecated.
72
#' @param line.colour Deprecated. specify the fill argument instead in
73
#'   \code{ribbon.params}. e.g. \code{ribbon.params = list(fill = 'red')}.
74
#' @param point.size Deprecated. specify the size argument instead in the
75
#'   \code{case.params} or \code{base.params} or \code{competing.params}
76
#'   argument. e.g. \code{case.params = list(size = 1.5)}.
77
#' @param point.colour Deprecated. Specify the values argument instead in the
78
#'   \code{color.params} and \code{fill.params} argument. See examples for
79
#'   details.
80
#' @param ncol Deprecated. Use \code{facet.params} instead.
81
#' @return The methods for \code{plot} return a population time plot, stratified
82
#'   by exposure status in the case of \code{popTimeExposure}. Note that these
83
#'   are \code{ggplot2} objects and can therefore be used in subsequent ggplot2
84
#'   type plots. See examples and vignette for details.
85
#' @details This function leverages the \code{ggplot2} package to build
86
#'   population time plots. It builds the plot by adding layers, starting with a
87
#'   layer for the area representing the population time. It then sequentially
88
#'   adds points to the plots to show the casebase sampling mechanism. This
89
#'   function gives user the flexibility to add any combination of the
90
#'   case.series, base.series and competing events. The case series and
91
#'   competing events are sampled at random vertically on the plot in order to
92
#'   visualise the incidence density using the \code{\link{popTime}} function.
93
#'   That is, imagine we draw a vertical line at a specific event time. We then
94
#'   plot the point at a randomly sampled y-coordinate along this vertical line.
95
#'   This is done to avoid having all points along the upper edge of the plot
96
#'   (because the subjects with the least amount of observation time are plotted
97
#'   at the top of the y-axis). By randomly distributing them, we can get a
98
#'   better sense of the inicidence density. The base series is sampled
99
#'   horizontally on the plot using the \code{\link{sampleCaseBase}} function.
100
#' @importFrom data.table := copy
101
#' @importFrom ggplot2 ggplot geom_point scale_fill_manual geom_ribbon aes
102
#' @importFrom ggplot2 scale_colour_manual element_blank facet_wrap
103
#' @importFrom ggplot2 theme xlab ylab
104
#' @seealso
105
#' \link{geom_point},\link{geom_ribbon},\link{theme},
106
#' \link{scale_colour_manual}, \link{scale_fill_manual},
107
#' \link{sampleCaseBase}
108
#' @examples
109
#' # change color of points
110
#' library(ggplot2)
111
#' data("bmtcrr")
112
#' popTimeData <- popTime(data = bmtcrr, time = "ftime", event = "Status")
113
#' fill_cols <- c("Case series" = "black", "Competing event" = "#009E73",
114
#'                "Base series" = "#0072B2")
115
#' color_cols <- c("Case series" = "black", "Competing event" = "black",
116
#'                 "Base series" = "black")
117
#'
118
#' plot(popTimeData,
119
#'   add.case.series = TRUE,
120
#'   add.base.series = TRUE,
121
#'   add.competing.event = FALSE,
122
#'   legend = TRUE,
123
#'   comprisk = TRUE,
124
#'   fill.params = list(
125
#'     name = element_blank(),
126
#'     breaks = c("Case series", "Competing event", "Base series"),
127
#'     values = fill_cols
128
#'   ),
129
#'   color.params = list(
130
#'     name = element_blank(),
131
#'     breaks = c("Case series", "Competing event", "Base series"),
132
#'     values = color_cols
133
#'   )
134
#' )
135
#' @export
136
#' @method plot popTime
137
#' @rdname popTime
138
plot.popTime <- function(x, ...,
139
                         xlab = "Follow-up time",
140
                         ylab = "Population",
141
                         add.case.series = TRUE,
142
                         add.base.series = FALSE,
143
                         add.competing.event = FALSE,
144
                         casebase.theme = TRUE,
145
                         ribbon.params = list(),
146
                         case.params = list(),
147
                         base.params = list(),
148
                         competing.params = list(),
149
                         color.params = list(),
150
                         fill.params = list(),
151
                         theme.params = list(),
152
                         facet.params = list(),
153
                         ratio = 1,
154
                         censored.indicator,
155
                         comprisk = FALSE,
156
                         legend = TRUE,
157
                         ncol, # Deprecated
158
                         legend.position, # Deprecated
159
                         line.width, # Deprecated
160
                         line.colour, # Deprecated
161
                         point.size, # Deprecated
162
                         point.colour) { # Deprecated
163

164
    # To prevent "no visble binding for global variable"
165 2
    ycoord <- yc <- event <- comprisk.event <- NULL
166

167
    # Okabe Ito colors
168
    # fill_cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2",
169
    #                "#999999")
170
    # colorspace::darken(fill_cols, 0.3)
171
    # color_cols <- c("#9D6C06", "#077DAA", "#026D4E", "#A39A09", "#044F7E",
172
    #                 "#696969")
173

174
    # cbbPalette
175
    # fill_cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442",
176
    #                "#0072B2", "#D55E00", "#CC79A7")
177
    # colorspace::qualitative_hcl(n = 3, palette = "Dark3") %>% put
178 2
    fill_cols <- c("#E16A86", "#50A315", "#009ADE")
179
    # colorspace::darken(fill_cols, 0.3) %>% dput
180
    # color_cols <- c("#696969", "#9D6C06", "#077DAA", "#026D4E", "#A39A09",
181
    #                 "#044F7E", "#954000", "#984B77")
182 2
    color_cols <- c("#AB3A59", "#347004", "#026A9A")
183
    # par(mfrow=c(1,2))
184
    # plot(seq_along(color_cols),col = color_cols, pch = 19, cex = 2.5)
185
    # plot(seq_along(fill_cols),col = fill_cols, pch = 19, cex = 2.5)
186

187 2
    if (!missing(line.colour)) {
188 0
        warning(paste("line.colour argument deprecated. specify the fill",
189 0
                      "argument instead in the ribbon.params.",
190 0
                      "e.g. ribbon.params = list(fill = 'red')."))
191
    }
192

193 2
    if (!missing(line.width)) {
194 0
        warning("line.width argument deprecated.")
195
    }
196

197 2
    if (!missing(point.size)) {
198 0
        warning(paste("point.size argument deprecated. specify the size",
199 0
                      "argument instead in the case.params or base.params or",
200 0
                      "competing.params argument.",
201 0
                      "e.g. case.params = list(size = 1.5)."))
202
    }
203

204 2
    if (!missing(point.colour)) {
205 0
        warning(paste("point.colour argument deprecated. specify the values",
206 0
                      "argument instead in the fill.params and color.params",
207 0
                      "arguments. see examples for details."))
208
    }
209

210 2
    if (!missing(legend.position)) {
211 0
        warning(paste("legend.position argument deprecated. specify the",
212 0
                      "legend.position argument instead in the theme.params",
213 0
                      "argument.",
214 0
                      "e.g. theme.params = list(legend.position = 'bottom')."))
215
    }
216

217 2
    if (!missing(ncol)) {
218 0
        warning(paste("ncol argument deprecated. specify the ncol argument",
219 0
                      "instead in the facet.params argument.",
220 0
                      "e.g. facet.params = list(ncol = 1)."))
221
    }
222

223 2
    if (length(fill.params) > 0 & length(color.params) == 0) {
224 0
        warning(paste("fill.params has been specified by the user but",
225 0
                      "color.params has not. Setting color.params to be equal",
226 0
                      "to fill.params."))
227 0
        color.params <- fill.params
228
    }
229

230 2
    if (length(color.params) > 0 & length(fill.params) == 0) {
231 0
        warning(paste("color.params has been specified by the user but",
232 0
                      "fill.params has not. Setting fill.params to be equal",
233 0
                      "to color.params."))
234 0
        fill.params <- color.params
235
    }
236

237 2
    exposure_variable <- attr(x, "exposure")
238

239 2
    p2 <- list(
240

241
        # Add poptime area -----------------------------------------------------
242 2
        do.call("geom_ribbon", utils::modifyList(
243 2
            list(
244 2
                data = x,
245 2
                mapping = aes(x = time, ymin = 0, ymax = ycoord),
246 2
                fill = "grey80",
247 2
                alpha = 0.5
248
            ),
249 2
            ribbon.params
250
        )),
251

252
        # Add case series ------------------------------------------------------
253 2
        if (add.case.series) {
254 2
            do.call("geom_point", utils::modifyList(
255 2
                list(
256 2
                    data = x[event == 1],
257 2
                    mapping = aes(x = time, y = yc, color = "Case series",
258 2
                                  fill = "Case series"),
259 2
                    show.legend = legend,
260 2
                    size = 1.5,
261 2
                    alpha = 0.5,
262 2
                    shape = 21
263
                ),
264 2
                case.params
265
            ))
266
        },
267

268

269
        # Add base series ------------------------------------------------------
270 2
        if (add.base.series) {
271 2
            basedata <- sampleCaseBase(
272 2
                data = x, time = "time", event = "event", ratio = ratio,
273 2
                comprisk = comprisk,
274 2
                censored.indicator = censored.indicator
275
            )
276

277 2
            do.call("geom_point", utils::modifyList(
278 2
                list(
279 2
                    data = basedata[event == 0],
280 2
                    mapping = aes(x = time, y = ycoord, colour = "Base series",
281 2
                                  fill = "Base series"),
282 2
                    show.legend = legend,
283 2
                    size = 1.5,
284 2
                    alpha = 0.5,
285 2
                    shape = 21
286
                ),
287 2
                base.params
288
            ))
289
        },
290

291

292
        # Add competing event --------------------------------------------------
293 2
        if (add.competing.event) {
294 2
            newX <- data.table::copy(x)
295 2
            newX[, `:=`(
296 2
                original.time = NULL,
297 2
                original.event = NULL,
298 2
                `event status` = NULL,
299 2
                ycoord = NULL,
300 2
                yc = NULL,
301 2
                n_available = NULL
302
            )]
303 2
            newX[event == 0, comprisk.event := 0]
304 2
            newX[event == 1, comprisk.event := 2]
305 2
            newX[event == 2, comprisk.event := 1]
306 2
            newX[, event := NULL]
307

308 2
            if (is.null(exposure_variable)) {
309 2
                compdata <- popTime(data = newX, time = "time",
310 2
                                    event = "comprisk.event")
311
            } else {
312 2
                compdata <- popTime(
313 2
                    data = newX, time = "time", event = "comprisk.event",
314 2
                    exposure = exposure_variable
315
                )
316
            }
317

318 2
            do.call("geom_point", utils::modifyList(
319 2
                list(
320 2
                    data = compdata[event == 1],
321 2
                    mapping = aes(x = time, y = yc, colour = "Competing event",
322 2
                                  fill = "Competing event"),
323 2
                    show.legend = legend,
324 2
                    size = 1.5,
325 2
                    alpha = 0.5,
326 2
                    shape = 21
327
                ),
328 2
                competing.params
329
            ))
330
        },
331

332
        # Add legend -----------------------------------------------------------
333 2
        if (legend) {
334

335
            # cols <- c("Case series" = color_cols[7],
336
            #           "Competing event" = color_cols[4],
337
            #           "Base series" = color_cols[6])
338

339 2
            cols <- c(
340 2
                "Case series" = color_cols[1],
341 2
                "Competing event" = color_cols[3],
342 2
                "Base series" = color_cols[2]
343
            )
344

345 2
            do.call("scale_colour_manual", utils::modifyList(
346 2
                list(
347 2
                    name = element_blank(),
348 2
                    breaks = c("Case series", "Competing event", "Base series"),
349 2
                    values = cols
350 2
                ), color.params
351
            ))
352
        },
353

354 2
        if (legend) {
355

356
            # cols <- c("Case series" = fill_cols[7],
357
            #           "Competing event" = fill_cols[4],
358
            #           "Base series" = fill_cols[6])
359

360 2
            cols <- c(
361 2
                "Case series" = fill_cols[1],
362 2
                "Competing event" = fill_cols[3],
363 2
                "Base series" = fill_cols[2]
364
            )
365

366 2
            do.call("scale_fill_manual", utils::modifyList(
367 2
                list(
368 2
                    name = element_blank(),
369 2
                    breaks = c("Case series", "Competing event", "Base series"),
370 2
                    values = cols
371 2
                ), fill.params
372
            ))
373
        },
374

375

376
        # Exposure stratified --------------------------------------------------
377 2
        if (!is.null(exposure_variable)) {
378 2
            do.call("facet_wrap", utils::modifyList(
379 2
                list(facets = exposure_variable, ncol = 1),
380 2
                facet.params
381
            ))
382
        },
383

384
        # Use casebase theme or not --------------------------------------------
385 2
        if (casebase.theme) {
386 2
            if (!is.null(exposure_variable)) {
387 0
                theme_cb() + panelBorder()
388
            } else {
389 2
                theme_cb()
390
            }
391
        },
392

393
        # add theme stuff  -----------------------------------------------------
394 2
        do.call("theme", utils::modifyList(
395 2
            list(legend.position = "bottom"),
396 2
            theme.params
397
        ))
398
    )
399

400 2
    p1 <- ggplot()
401

402 2
    p1 + p2 + ggplot2::xlab(xlab) + ggplot2::ylab(ylab)
403
}
404

405

406

407
#' @title Plot Hazards and Hazard Ratios
408
#'
409
#' @description Plot method for objects returned by the \code{fitSmoothHazard}
410
#'   function. Current plot types are hazard function and hazard ratio. The
411
#'   \code{visreg} package must be installed for \code{type="hazard"}. This
412
#'   function accounts for the possible time-varying exposure effects.
413
#' @param x Fitted object of class `glm`, `gam`, `cv.glmnet` or `gbm`. This is
414
#'   the result from the [casebase::fitSmoothHazard()] function.
415
#' @param ... further arguments passed to `plot`. Only used if \code{type="hr"}.
416
#'   Any of `lwd`,`lty`,`col`,`pch`,`cex` will be applied to the hazard ratio
417
#'   line, or point (if only one time point is supplied to `newdata`).
418
#' @param type plot type. Choose one of either \code{"hazard"} for hazard
419
#'   function or \code{"hr"} for hazard ratio.  Default: \code{type = "hazard"}.
420
#' @param hazard.params Named list of arguments which will override the defaults
421
#'   passed to [visreg::visreg()], The default arguments are \code{list(fit = x,
422
#'   trans = exp, plot = TRUE, rug = FALSE, alpha = 1, partial = FALSE, overlay
423
#'   = TRUE)}. For example, if you want a 95% confidence band, specify
424
#'   \code{hazard.params = list(alpha = 0.05)}. Note that The `cond` argument
425
#'   must be provided as a named list. Each element of that list specifies the
426
#'   value for one of the terms in the model; any elements left unspecified are
427
#'   filled in with the median/most common category. Only used for
428
#'   `type="hazard"`. All other argument are used for `type="hr"`. Note that the
429
#'   `visreg` package must be installed for `type="hazard"`.
430
#' @param newdata Required for \code{type="hr"}. The \code{newdata} argument is
431
#'   the "unexposed" group, while the exposed group is defined by either: (i) a
432
#'   change (defined by the \code{increment} argument) in a variable in newdata
433
#'   defined by the \code{var} argument ; or (ii) an exposed function that takes
434
#'   a data-frame and returns the "exposed" group (e.g. \code{exposed =
435
#'   function(data) transform(data, treat=1)}). This is a generalization of the
436
#'   behavior of the rstpm2 plot function. It allows both numeric and factor
437
#'   variables to be incremented or decremented. See references for rstpm2
438
#'   package. Only used for `type="hr"`
439
#' @param var specify the variable name for the exposed/unexposed (name is given
440
#'   as a character variable). If this argument is missing, then the
441
#'   \code{exposed} argument must be specified. This is the variable which will
442
#'   be incremented by the \code{increment} argument to give the exposed
443
#'   category. If \code{var} is coded as a factor variable, then
444
#'   \code{increment=1} will return the next level of the variable in
445
#'   \code{newdata}. \code{increment=2} will return two levels above, and so on.
446
#'   If the value supplied to \code{increment} is greater than the number of
447
#'   levels, this will simply return the max level. You can also decrement the
448
#'   categorical variable by specifying a negative value, e.g.,
449
#'   \code{increment=-1} will return one level lower than the value in
450
#'   \code{newdata}. If \code{var} is a numeric, than \code{increment} will
451
#'   increment (if positive) or decrement (if negative) by the supplied value.
452
#'   Only used for `type="hr"`.
453
#' @param increment Numeric value indicating how much to increment (if positive)
454
#'   or decrement (if negative) the \code{var} variable in \code{newdata}. See
455
#'   \code{var} argument for more details. Default is 1. Only used for
456
#'   `type="hr"`.
457
#' @param exposed function that takes \code{newdata} and returns the exposed
458
#'   dataset (e.g. function(data) transform(data, treat = 1)). This argument
459
#'   takes precedence over the \code{var} argument, i.e., if both \code{var} and
460
#'   \code{exposed} are correctly specified, only the \code{exposed} argument
461
#'   will be used. Only used for `type="hr"`.
462
#' @param xvar Variable to be used on x-axis for hazard ratio plots. If NULL,
463
#'   the function defaults to using the time variable used in the call to
464
#'   \code{fitSmoothHazard}. In general, this should be any continuous variable
465
#'   which has an interaction term with another variable. Only used for
466
#'   `type="hr"`.
467
#' @param ci Logical; if TRUE confidence bands are calculated. Only available
468
#'   for `family="glm"` and `family="gam"`, and only used for `type="hr"`,
469
#'   Default: !add. Confidence intervals for hazard ratios are calculated using
470
#'   the Delta Method.
471
#' @param ci.lvl Confidence level. Must be in (0,1), Default: 0.95. Only used
472
#'   for `type="hr"`.
473
#' @param ci.col Confidence band color. Only used if argument `ci=TRUE`,
474
#'   Default: 'grey'. Only used for `type="hr"`.
475
#' @param rug Logical. Adds a rug representation (1-d plot) of the event times
476
#'   (only for `status=1`), Default: !ci. Only used for `type="hr"`.
477
#' @return a plot of the hazard function or hazard ratio. For `type="hazard"`, a
478
#'   `data.frame` (returned invisibly) of the original data used in the fitting
479
#'   along with the data used to create the plots including `predictedhazard`
480
#'   which is the predicted hazard for a given covariate pattern and time.
481
#'   `predictedloghazard` is the predicted hazard on the log scale. `lowerbound`
482
#'   and `upperbound` are the lower and upper confidence interval bounds on the
483
#'   hazard scale (i.e. used to plot the confidence bands). `standarderror` is
484
#'   the standard error of the log hazard or log hazard ratio (only if
485
#'   `family="glm"` or `family="gam"`). For `type="hr"`, `log_hazard_ratio` and
486
#'   `hazard_ratio` is returned, and if `ci=TRUE`, `standarderror` (on the log
487
#'   scale) and `lowerbound` and `upperbound` of the `hazard_ratio` are
488
#'   returned.
489
#' @details This function has only been thoroughly tested for `family="glm"`. If
490
#'   the user wants more customized plot aesthetics, we recommend saving the
491
#'   results to a `data.frame` and using  the graphical package of their choice.
492
#' @examples
493
#' if (requireNamespace("splines", quietly = TRUE)) {
494
#' data("simdat") # from casebase package
495
#' library(splines)
496
#' simdat <- transform(simdat[sample(1:nrow(simdat), size = 200),],
497
#'                     treat = factor(trt, levels = 0:1,
498
#'                     labels = c("control","treatment")))
499
#'
500
#' fit_numeric_exposure <- fitSmoothHazard(status ~ trt*bs(eventtime),
501
#'                                         data = simdat,
502
#'                                         ratio = 1,
503
#'                                         time = "eventtime")
504
#'
505
#' fit_factor_exposure <- fitSmoothHazard(status ~ treat*bs(eventtime),
506
#'                                        data = simdat,
507
#'                                        ratio = 1,
508
#'                                        time = "eventtime")
509
#'
510
#' newtime <- quantile(fit_factor_exposure[["data"]][[fit_factor_exposure[["timeVar"]]]],
511
#'                     probs = seq(0.05, 0.95, 0.01))
512
#'
513
#' par(mfrow = c(1,3))
514
#' plot(fit_numeric_exposure,
515
#'      type = "hr",
516
#'      newdata = data.frame(trt = 0, eventtime = newtime),
517
#'      exposed = function(data) transform(data, trt = 1),
518
#'      xvar = "eventtime",
519
#'      ci = TRUE)
520
#'
521
#' #by default this will increment `var` by 1 for exposed category
522
#' plot(fit_factor_exposure,
523
#'      type = "hr",
524
#'      newdata = data.frame(treat = factor("control",
525
#'               levels = c("control","treatment")), eventtime = newtime),
526
#'      var = "treat",
527
#'      increment = 1,
528
#'      xvar = "eventtime",
529
#'      ci = TRUE,
530
#'      ci.col = "lightblue",
531
#'      xlab = "Time",
532
#'      main = "Hazard Ratio for Treatment",
533
#'      ylab = "Hazard Ratio",
534
#'      lty = 5,
535
#'      lwd = 7,
536
#'      rug = TRUE)
537
#'
538
#'
539
#' # we can also decrement `var` by 1 to give hazard ratio for control/treatment
540
#' result <- plot(fit_factor_exposure,
541
#'                type = "hr",
542
#'                newdata = data.frame(treat = factor("treatment",
543
#'                                     levels = c("control","treatment")),
544
#'                                     eventtime = newtime),
545
#'                var = "treat",
546
#'                increment = -1,
547
#'                xvar = "eventtime",
548
#'                ci = TRUE)
549
#'
550
#' # see data used to create plot
551
#' head(result)
552
#' }
553
#' @seealso \code{\link[utils]{modifyList}}, [casebase::fitSmoothHazard()],
554
#'   [graphics::par()], [visreg::visreg()]
555
#' @rdname plot.singleEventCB
556
#' @references Mark Clements and Xing-Rong Liu (2019). rstpm2: Smooth Survival
557
#'   Models, Including Generalized Survival Models. R package version 1.5.1.
558
#'   https://CRAN.R-project.org/package=rstpm2
559
#'
560
#'   Breheny P and Burchett W (2017). Visualization of Regression Models Using
561
#'   visreg. The R Journal, 9: 56-71.
562
#' @export
563
#' @importFrom utils modifyList
564
plot.singleEventCB <- function(x, ...,
565
                               type = c("hazard", "hr"),
566
                               hazard.params = list(),
567
                               newdata,
568
                               exposed,
569
                               increment = 1,
570
                               var,
571
                               xvar = NULL,
572
                               ci = FALSE,
573
                               ci.lvl = 0.95,
574
                               rug = !ci,
575
                               ci.col = "grey") {
576 2
    type <- match.arg(type)
577

578 2
    if (type == "hazard") {
579 2
        if (!requireNamespace("visreg", quietly = TRUE)) {
580 0
            stop("visreg package needed for this function. please install it first.")
581
        }
582

583 2
        tt <- do.call("visreg", utils::modifyList(
584 2
            list(
585 2
                fit = x,
586 2
                trans = exp,
587 2
                plot = T,
588 2
                rug = FALSE,
589 2
                alpha = 1,
590 2
                partial = FALSE,
591 2
                overlay = TRUE,
592 2
                print.cond = TRUE
593
            ),
594 2
            hazard.params
595
        ))
596

597 2
        invisible(tt)
598

599 2
        return(tt)
600
    }
601

602

603 2
    if (type == "hr") {
604 2
        if (is.null(newdata) && type %in% c("hr")) {
605 0
            stop("Prediction using type 'hr' requires newdata to be specified.")
606
        }
607

608 2
        check_arguments_hazard(
609 2
            object = x, newdata = newdata, plot = FALSE,
610 2
            ci = ci, ci.lvl = ci.lvl
611
        )
612

613 2
        if (missing(exposed) & missing(var)) {
614 2
            stop("One of 'var' or 'exposed' arguments must be specified.")
615
        }
616

617 2
        if (missing(exposed)) {
618 2
            if (var %ni% names(newdata)) {
619 2
                stop(sprintf("%s not found in 'newdata'.", var))
620
            }
621 2
            newdata2 <- incrVar(var, increment = increment)(newdata)
622 2
        } else if (is.function(exposed)) {
623 2
            if (!missing(var)) warning("'var' argument is being ignored since 'exposure' argment has been correctly specified.")
624 2
            newdata2 <- exposed(newdata)
625 2
        } else if (!is.function(exposed) & !missing(var)) {
626 2
            if (var %ni% names(newdata)) {
627 2
                stop(sprintf("%s not found in 'newdata'.", var))
628
            }
629 2
            warning("'exposure' argument ignored since it is not correctly specified. Using 'var' argument instead.")
630 2
            newdata2 <- incrVar(var, increment = increment)(newdata)
631
        } else {
632 0
            stop("incorrect specification for 'exposed' argument. see help page for details.")
633
        }
634

635 2
        check_arguments_hazard(
636 2
            object = x, newdata = newdata2, plot = FALSE,
637 2
            ci = ci, ci.lvl = ci.lvl
638
        )
639

640 2
        plotHazardRatio(
641 2
            x = x, newdata = newdata, newdata2 = newdata2, ci = ci,
642 2
            ci.lvl = ci.lvl, ci.col = ci.col, rug = rug, xvar = xvar, ...
643
        )
644
    }
645
}
646

647
#' @title Plot Cumulative Incidence and Survival Curves
648
#' @description Plot method for objects returned by the \code{absoluteRisk}
649
#'   function. Current plot types are cumulative incidence and survival
650
#'   functions.
651
#' @param x Fitted object of class `absRiskCB`. This is the result from the
652
#'   [casebase::absoluteRisk()] function.
653
#' @param ... further arguments passed to `matplot`. Only used if
654
#'   \code{gg=FALSE}.
655
#' @param xlab xaxis label, Default: 'time'
656
#' @param ylab yaxis label. By default, this will use the `"type"` attribute of
657
#'   the `absRiskCB` object
658
#' @param type Line type. Only used if `gg = FALSE`. This argument gets passed
659
#'   to [graphics::matplot()]. Default: 'l'
660
#' @param gg Logical for whether the `ggplot2` package should be used for
661
#'   plotting. Default: TRUE
662
#' @param id.names Optional character vector used as legend key when `gg=TRUE`.
663
#'   If missing, defaults to V1, V2, ...
664
#' @param legend.title Optional character vector of the legend title. Only used
665
#'   if `gg = FALSE`. Default is `'ID'`
666
#' @return A plot of the cumulative incidence or survival curve
667
#' @seealso \code{\link[graphics]{matplot}},
668
#'   \code{\link[casebase]{absoluteRisk}},
669
#'   \code{\link[data.table]{as.data.table}}, \code{\link[data.table]{setattr}},
670
#'   \code{\link[data.table]{melt.data.table}}
671
#' @rdname absoluteRisk
672
#' @export
673
#' @importFrom graphics matplot
674
#' @importFrom ggplot2 ggplot aes geom_line labs theme
675
#' @importFrom data.table as.data.table setnames melt
676
#' @examples
677
#' # Plot CI curves----
678
#' library(ggplot2)
679
#' data("brcancer")
680
#' mod_cb_tvc <- fitSmoothHazard(cens ~ estrec*log(time) +
681
#'                                 horTh +
682
#'                                 age +
683
#'                                 menostat +
684
#'                                 tsize +
685
#'                                 tgrade +
686
#'                                 pnodes +
687
#'                                 progrec,
688
#'                               data = brcancer,
689
#'                               time = "time", ratio = 1)
690
#' smooth_risk_brcancer <- absoluteRisk(object = mod_cb_tvc,
691
#'                                      newdata = brcancer[c(1,50),])
692
#'
693
#' class(smooth_risk_brcancer)
694
#' plot(smooth_risk_brcancer)
695
plot.absRiskCB <- function(x, ...,
696
                           xlab = "time",
697
                           ylab = ifelse(attr(x, "type") == "CI",
698
                                         "cumulative incidence",
699
                                         "survival probability"),
700
                           type = "l",
701
                           gg = TRUE,
702
                           id.names,
703
                           legend.title) {
704
    # output of absoluteRisk will always have a column named time
705
    # x equals linearRisk
706
    # ======================
707

708 2
    if (!gg) {
709 2
        graphics::matplot(x = x[, "time"],
710 2
                          y = x[, -which(colnames(x) == c("time"))],
711 2
                          type = type,
712 2
                          xlab = xlab,
713 2
                          ylab = ylab,
714
                          ...
715
        )
716

717
    } else {
718

719 2
        ID <- value <- NULL
720 2
        DT <- data.table::as.data.table(x)
721

722 2
        if (!missing(id.names)) {
723 2
            names_to_change <- grep("V", colnames(DT), value = TRUE)
724 2
            if (length(names_to_change) != length(id.names)) {
725 2
                warning("length of 'id.names' not equal to number of covariate profiles. ignoring 'id.names' argument")
726
            } else {
727 2
                data.table::setnames(DT, old = names_to_change, new = id.names)
728
            }
729
        }
730

731 2
        DT.m <- data.table::melt(DT, id.vars = "time", variable.name = "ID")
732

733 2
        ggplot(DT.m, aes(x = time, y = value, color = ID)) +
734 2
            geom_line() + theme_cb() +
735 2
            labs(color = ifelse(missing(legend.title), "ID", legend.title),
736 2
                 x = xlab, y = ylab) +
737 2
            theme(legend.position = "bottom")
738
    }
739

740
}

Read our documentation on viewing source code .

Loading