sahirbhatnagar / casebase
1
#' @rdname plot.singleEventCB
2
#' @param newdata2 `data.frame` for exposed group. calculated and passed
3
#'   internally to `plotHazardRatio` function
4
#' @importFrom graphics polygon lines arrows points
5
#' @importFrom stats terms delete.response vcov qnorm
6
#' @importFrom utils modifyList
7
plotHazardRatio <- function(x, newdata, newdata2, ci, ci.lvl, ci.col,
8
                            rug, xvar, ...) {
9 2
    tt <- stats::terms(x)
10 2
    Terms <- stats::delete.response(tt)
11 2
    beta2 <- coef(x)
12 2
    gradient <- hrJacobian(
13 2
        object = x, newdata = newdata,
14 2
        newdata2 = newdata2, term = Terms
15
    )
16

17 2
    log_hazard_ratio <- gradient %*% beta2
18

19 2
    if (is.null(xvar)) {
20 0
        xvar <- x[["timeVar"]]
21
    } else {
22 2
        if (length(xvar) > 1) {
23 0
            warning(paste("more than one xvar supplied. Only plotting hazard",
24 0
                          "ratio for first element."))
25
            }
26 2
        xvar <- xvar[1]
27
    }
28

29 2
    if (xvar %ni% colnames(newdata)) {
30 0
        stop(sprintf(paste("%s column (which you supplied to 'xvar' argument)",
31 0
                           "not found in newdata"), xvar))
32
    } else {
33 2
        xvar_values <- newdata[[xvar]]
34
    }
35

36
    # sorting indices for ploting
37 2
    i.backw <- order(xvar_values, decreasing = TRUE)
38 2
    i.forw <- order(xvar_values)
39

40 2
    other_plot_args <- list(...)
41 2
    line_args <- grep("^lwd$|^lty$|^col$|^pch$|^cex$", names(other_plot_args))
42

43 2
    if (length(line_args) == 0) {
44 2
        line_args <- list(NULL)
45
    } else {
46 2
        line_args <- other_plot_args[line_args]
47
    }
48

49
    # plot CI as polygon shade - if 'se = TRUE' (default)
50 2
    if (ci) {
51 2
        v2 <- stats::vcov(x)
52 2
        SE_log_hazard_ratio <- sqrt(diag(gradient %*% tcrossprod(v2, gradient)))
53

54 2
        hazard_ratio_lower <- exp(stats::qnorm(p = (1 - ci.lvl) / 2,
55 2
                                               mean = log_hazard_ratio,
56 2
                                               sd = SE_log_hazard_ratio))
57 2
        hazard_ratio_upper <- exp(stats::qnorm(p = 1 - (1 - ci.lvl) / 2,
58 2
                                               mean = log_hazard_ratio,
59 2
                                               sd = SE_log_hazard_ratio))
60 2
        x.poly <- c(xvar_values[i.forw], xvar_values[i.backw])
61 2
        y.poly <- c(hazard_ratio_lower[i.forw], hazard_ratio_upper[i.backw])
62

63 2
        do.call("plot", utils::modifyList(
64 2
            list(
65 2
                x = range(x.poly),
66 2
                y = range(y.poly) * c(0.99, 1.01),
67 2
                type = "n",
68 2
                ylab = "hazard ratio",
69 2
                xlab = xvar
70
            ),
71 2
            other_plot_args
72
        ))
73

74 2
        if (length(unique(x.poly)) == 1) {
75 2
            graphics::arrows(x0 = x.poly[1], y0 = y.poly[1], y1 = y.poly[2],
76 2
                             angle = 90, length = 0.5, code = 3)
77 2
            do.call("points", utils::modifyList(
78 2
                list(
79 2
                    x = xvar_values[i.forw],
80 2
                    y = exp(log_hazard_ratio[i.forw]),
81 2
                    pch = 19,
82 2
                    col = "black",
83 2
                    cex = 2
84
                ),
85 2
                line_args
86
            ))
87
        } else {
88 2
            graphics::polygon(x.poly, y.poly, col = ci.col, border = NA)
89 2
            do.call("lines", utils::modifyList(
90 2
                list(
91 2
                    x = xvar_values[i.forw],
92 2
                    y = exp(log_hazard_ratio[i.forw]),
93 2
                    lwd = 2,
94 2
                    lty = 1,
95 2
                    col = "black"
96
                ),
97 2
                line_args
98
            ))
99
        }
100

101 2
        results <- transform(newdata,
102 2
                             log_hazard_ratio = log_hazard_ratio,
103 2
                             standarderror = SE_log_hazard_ratio,
104 2
                             hazard_ratio = exp(log_hazard_ratio),
105 2
                             lowerbound = hazard_ratio_lower,
106 2
                             upperbound = hazard_ratio_upper
107
        )
108

109
    } else {
110

111 2
        if (length(xvar_values) == 1) {
112 2
            do.call("plot", utils::modifyList(
113 2
                list(
114 2
                    x = xvar_values, y = exp(log_hazard_ratio), lwd = 2,
115 2
                    lty = 1, pch = 19, cex = 2, ylab = "hazard ratio",
116 2
                    xlab = xvar
117
                ),
118 2
                other_plot_args
119
            ))
120

121
        } else {
122

123 2
            do.call("plot", utils::modifyList(
124 2
                list(
125 2
                    x = xvar_values, y = exp(log_hazard_ratio), lwd = 2,
126 2
                    lty = 1, type = "l", ylab = "hazard ratio", xlab = xvar
127
                ),
128 2
                other_plot_args
129
            ))
130
        }
131

132 2
        results <- transform(newdata,
133 2
                             log_hazard_ratio = log_hazard_ratio,
134 2
                             hazard_ratio = exp(log_hazard_ratio)
135
        )
136
    }
137

138 2
    if (rug) {
139 2
        events <- x[["originalData"]][[x[["eventVar"]]]]
140 2
        rug(x[["originalData"]][which(events == 1), , drop = FALSE][[xvar]],
141 2
            quiet = TRUE
142 2
        ) # Silence warnings about clipped values
143
    }
144

145 2
    invisible(results)
146
}

Read our documentation on viewing source code .

Loading