method_expTransform argument for plot_APCheatmap

R/plot_APCheatmap.R

@@ -53,6 +53,12 @@
 #' Only used if \code{y_var} is specified. Defaults to FALSE.
 #' @param plot_CI Indicator if the confidence intervals should be plotted.
 #' Only used if \code{y_var} is not specified. Defaults to TRUE.
+#' @param method_expTransform One of \code{c("simple","delta")}, stating if
+#' confidence interval limits should be transformed by
+#' a simple exp transformation or using the delta method. The delta method can
+#' be unstable in situations and lead to negative confidence interval limits.
+#' Only used when the model was estimated with a log or logit link and
+#' confidence intervals are supposed to be plotted. Defaults to \code{simple}.
 #' @param legend_limits Optional numeric vector passed as argument \code{limits}
 #' to \code{\link[ggplot2]{scale_fill_gradient2}}.
 #' 
@@ -126,6 +132,7 @@ plot_APCheatmap <- function(dat, y_var = NULL, model = NULL,
                             markLines_list = NULL,
                             markLines_displayLabels = c("age","period","cohort"),
                             y_var_logScale = FALSE, plot_CI = TRUE,
+                            method_expTransform = "simple",
                             legend_limits = NULL) {
 
   checkmate::assert_data_frame(dat)
@@ -150,6 +157,7 @@ plot_APCheatmap <- function(dat, y_var = NULL, model = NULL,
   checkmate::assert_logical(y_var_logScale, len = 1)
   checkmate::assert_logical(plot_CI, len = 1)
   checkmate::assert_numeric(legend_limits, len = 2, null.ok = TRUE)
+  checkmate::assert_choice(method_expTransform, choices = c("simple","delta"))
 
 
   # some NULL definitions to appease CRAN checks regarding use of dplyr/ggplot2
@@ -233,15 +241,33 @@ plot_APCheatmap <- function(dat, y_var = NULL, model = NULL,
     legend_title <- ifelse(used_logLink, "Mean exp effect", "Mean effect")
     y_trans      <- ifelse(used_logLink, "log", "identity")
 
+    # Create lower and upper confidence surfaces:
     if (used_logLink) {
-      plot_dat <- plot_dat %>% 
-        mutate(exp_effect = exp(effect),
-               exp_se     = sqrt((se^2) * (exp_effect^2))) %>% 
-        mutate(exp_lower  = exp_effect - qnorm(0.975) * exp_se,
-               exp_upper  = exp_effect + qnorm(0.975) * exp_se) %>% 
-        select(-effect, -se, -upper, -lower) %>% 
-        dplyr::rename(effect = exp_effect, se = exp_se,
-                      upper  = exp_upper, lower = exp_lower)
+      if (method_expTransform == "simple") {
+        plot_dat <- plot_dat %>%
+          mutate(exp_effect = exp(effect),
+                 exp_lower = exp(lower),
+                 exp_upper = exp(upper)) %>%
+          select(-effect, -se, -upper, -lower) %>% 
+          dplyr::rename(effect = exp_effect, upper  = exp_upper,
+                        lower = exp_lower)
+      }
+      else { # method_expTransform == "delta"
+        plot_dat <- plot_dat %>% 
+          mutate(exp_effect = exp(effect),
+                 exp_se     = sqrt((se^2) * (exp_effect^2))) %>% 
+          mutate(exp_lower  = exp_effect - qnorm(0.975) * exp_se,
+                 exp_upper  = exp_effect + qnorm(0.975) * exp_se) %>% 
+          select(-effect, -se, -upper, -lower) %>% 
+          dplyr::rename(effect = exp_effect, se = exp_se,
+                        upper  = exp_upper, lower = exp_lower)
+      }
+      if (any(plot_dat$lower < 0)) {
+        warning("Note: After the delta method transformation some values of the
+              lower confidence surface were negative. These
+              values were set to 0.01")
+        plot_dat$lower[plot_dat$lower < 0] <- 0.01
+      }
     }
 
   }

tests/testthat/test-plots_APCsurface.R

@@ -17,10 +17,17 @@ test_that("plot_APCheatmap", {
                                                "period" = c(1990,2010),
                                                "cohort" = c(1985,1993)),
                          apc_range = list("cohort" = 1980:2010))
+  gg4 <- plot_APCheatmap(dat = drug_deaths, y_var = "mortality_rate",
+                         markLines_list = list("age"    = c(20,70),
+                                               "period" = c(1990,2010),
+                                               "cohort" = c(1985,1993)),
+                         apc_range = list("cohort" = 1980:2010),
+                         plot_CI = FALSE)
 
   expect_s3_class(gg1, class = c("gg","ggplot"))
   expect_s3_class(gg2, class = c("gg","ggplot"))
   expect_s3_class(gg3, class = c("gg","ggplot"))
+  expect_s3_class(gg4, class = c("gg","ggplot"))
 
 
   # plot heatmap of smoothed structure
@@ -31,9 +38,12 @@ test_that("plot_APCheatmap", {
 
   gg1 <- plot_APCheatmap(dat = drug_deaths, model = model)
   gg2 <- plot_APCheatmap(dat = drug_deaths, model = model_logLink)
+  gg3 <- plot_APCheatmap(dat = drug_deaths, model = model_logLink,
+                         method_expTransform = "delta")
 
   expect_s3_class(gg1, class = c("gg","ggplot"))
   expect_s3_class(gg2, class = c("gg","ggplot"))
+  expect_s3_class(gg3, class = c("gg","ggplot"))
 })