Added vignette on censoring and updated plots

1. Added a vignette on right-censoring
2. Updated the plot functions to ggplot2

DESCRIPTION

@@ -20,7 +20,7 @@ Suggests:
     testthat (>= 3.0.0),
     knitr, 
 	  DiagrammeR, 
-	  DT 
+	  DT
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
@@ -34,7 +34,8 @@ Imports:
 	  coda (>= 0.19-4),
 	  stats,
 	  grDevices,
-	  graphics
+	  graphics,
+	  ggplot2
 NeedsCompilation: no
 Config/testthat/edition: 3
 URL: https://kuan-liu-lab.github.io/bayesmsm/, https://github.com/Kuan-Liu-Lab/bayesmsm

NAMESPACE

@@ -8,6 +8,7 @@ export(plot_ATE)
 export(plot_est_box)
 export(summary_bayesmsm)
 import(doParallel)
+import(ggplot2)
 import(parallel)
 importFrom(MCMCpack,rdirichlet)
 importFrom(R2jags,jags)

R/plot_APO.R

@@ -52,7 +52,7 @@ plot_APO <- function(input, effect_type, ...) {
   effect <- bootdata[, effect_type, drop = FALSE]
 
   # Calculate density
-  density_effect <- density(effect[[1]])
+  density_effect <- stats::density(effect[[1]])
 
   # Define titles and colors based on effect_type
   titles <- c(effect_comparator = "Comparator Level", effect_reference = "Reference Level")
@@ -62,26 +62,28 @@ plot_APO <- function(input, effect_type, ...) {
   mean_effect <- mean(effect[[1]])
 
   # Calculate CI
-  ci <- quantile(effect[[1]], probs = c(0.025, 0.975))
-  density_ci <- density(effect[[1]], from = ci[1], to = ci[2])
+  ci <- stats::quantile(effect[[1]], probs = c(0.025, 0.975))
 
-  # Plotting
-  plot(density_effect, main = paste("Average Potential Outcome (APO) of", titles[effect_type]), xlab = "Effect", ylab = "Density", col = colors[effect_type], lwd = 2, ...)
+  # Create data frame for ggplot
+  density_data <- data.frame(x = density_effect$x, y = density_effect$y)
+  ci_data <- data.frame(x = c(ci[1], ci[2]), y = c(0, 0))
 
-  # Shade the area under the curve within the 95% CI
-  polygon(c(density_ci$x, rev(density_ci$x)), c(rep(min(density_effect$y), length(density_ci$x)), rev(density_ci$y)), col = rgb(0, 0, 1, alpha = 0.3))
-
-  # Add vertical lines for the mean and 95% CI bounds
-  abline(v = mean_effect, col = "darkgrey", lty = 3)
-  abline(v = ci[1], col = "darkgreen", lty = 2)
-  abline(v = ci[2], col = "darkgreen", lty = 2)
-
-  # Legend with mean and 95% CI bounds
-  legend_text <- c(titles[effect_type],
-                   paste("Mean:", round(mean_effect, 3)),
-                   paste("95% CI: [", round(ci[1], 3), ",", round(ci[2], 3), "]"))
-
-  legend("bottom", inset=c(0,-0.5),legend = legend_text,
-         col = c(colors[effect_type], "purple", "darkgreen"),
-         lwd = 2, lty = c(1, 3, 12))
+  # ggplot2 visualization
+  ggplot2::ggplot() +
+    ggplot2::geom_line(data = density_data, ggplot2::aes(x = x, y = y), color = colors[effect_type], size = 1) +
+    ggplot2::geom_ribbon(data = density_data, ggplot2::aes(x = x, ymin = 0, ymax = y), fill = "lightblue", alpha = 0.3) +
+    ggplot2::geom_vline(xintercept = mean_effect, color = "purple", linetype = "dashed", size = 1.2) +
+    ggplot2::geom_vline(xintercept = ci, color = "darkgreen", linetype = "dotted", size = 1.2) +
+    ggplot2::labs(title = paste("Average Potential Outcome (APO) of", titles[effect_type]), x = "Effect", y = "Density") +
+    ggplot2::theme_minimal() +
+    ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
+    ggplot2::annotate("text", x = mean_effect, y = max(density_data$y) * 0.9,
+                      label = paste("Mean:", round(mean_effect, 3)),
+                      color = "purple", angle = 90, vjust = -0.5) +
+    ggplot2::annotate("text", x = ci[1], y = max(density_data$y) * 0.8,
+                      label = paste("95% CI Lower:", round(ci[1], 3)),
+                      color = "darkgreen", angle = 90, vjust = -0.5) +
+    ggplot2::annotate("text", x = ci[2], y = max(density_data$y) * 0.8,
+                      label = paste("95% CI Upper:", round(ci[2], 3)),
+                      color = "darkgreen", angle = 90, vjust = -0.5)
 }

R/plot_ATE.R

@@ -10,6 +10,7 @@
 #' @param xlim Limits for the x-axis (default is NULL).
 #' @param ylim Limits for the y-axis (default is NULL).
 #' @param ... Additional graphical parameters passed to the plot function.
+#' @import ggplot2
 #' @importFrom stats density quantile
 #' @importFrom grDevices rgb
 #' @importFrom graphics abline arrows axis legend mtext par polygon text
@@ -48,21 +49,30 @@ plot_ATE <- function(input,
     }
   }
 
-  ate_density <- density(ate_values)
-  ci <- quantile(ate_values, probs = c(0.025, 0.975))
-  density_ci <- density(ate_values, from = ci[1], to = ci[2])
+  # Calculate density and CI
+  ate_density <- stats::density(ate_values)
+  ci <- stats::quantile(ate_values, probs = c(0.025, 0.975))
 
-  plot(ate_density, col = col_density, main = main, xlab = xlab, ylab = ylab, xlim = xlim, ylim = ylim, ...)
-  polygon(c(density_ci$x, rev(density_ci$x)), c(rep(min(ate_density$y), length(density_ci$x)), rev(density_ci$y)), col = rgb(0, 0, 1, alpha = 0.3))
-  abline(v = mean(ate_values), col = "purple", lwd = 2, lty = 3)
-  abline(v = ci[1], col = "darkgreen", lty = 2)
-  abline(v = ci[2], col = "darkgreen", lty = 2)
+  # Create data frame for ggplot
+  density_data <- data.frame(x = ate_density$x, y = ate_density$y)
+  ci_data <- data.frame(x = c(ci[1], ci[2]), y = c(0, 0))
 
-  legend_text <- c(paste(ATE, "Density",sep = " "),
-                   paste("Mean:", round(mean(ate_values), 3)),
-                   paste("95% CI: [", round(ci[1], 3), ",", round(ci[2], 3), "]"))
-
-  legend("bottom", inset=c(0,-0.5), legend = legend_text,
-         col = c(col_density, "purple", "darkgreen"),
-         lwd = 2, lty = c(1, 3, 2))
+  # ggplot2 visualization
+  ggplot2::ggplot() +
+    ggplot2::geom_line(data = density_data, ggplot2::aes(x = x, y = y), color = col_density, size = 1) +
+    ggplot2::geom_ribbon(data = density_data, ggplot2::aes(x = x, ymin = 0, ymax = y), fill = fill_density, alpha = 0.3) +
+    ggplot2::geom_vline(xintercept = mean(ate_values), color = "purple", linetype = "dashed", size = 1.2) +
+    ggplot2::geom_vline(xintercept = ci, color = "darkgreen", linetype = "dotted", size = 1.2) +
+    ggplot2::labs(title = main, x = xlab, y = ylab) +
+    ggplot2::theme_minimal() +
+    ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
+    ggplot2::annotate("text", x = mean(ate_values), y = max(density_data$y) * 0.9,
+                      label = paste("Mean:", round(mean(ate_values), 3)),
+                      color = "purple", angle = 90, vjust = -0.5) +
+    ggplot2::annotate("text", x = ci[1], y = max(density_data$y) * 0.8,
+                      label = paste("95% CI Lower:", round(ci[1], 3)),
+                      color = "darkgreen", angle = 90, vjust = -0.5) +
+    ggplot2::annotate("text", x = ci[2], y = max(density_data$y) * 0.8,
+                      label = paste("95% CI Upper:", round(ci[2], 3)),
+                      color = "darkgreen", angle = 90, vjust = -0.5)
 }

R/plot_est_box.R

@@ -47,43 +47,27 @@ plot_est_box <- function(input, ...) {
     stop("bootdata must contain 'effect_comparator', 'effect_reference', and 'RD' columns.")
   }
 
-  # Calculate means and standard deviations
+  # Calculate means and confidence intervals
   means <- sapply(bootdata[required_columns], mean)
   lowerbd <- sapply(bootdata[required_columns], function(x) quantile(x, probs = 0.025))
   upperbd <- sapply(bootdata[required_columns], function(x) quantile(x, probs = 0.975))
 
-  # Define the position for each point
-  position <- 1:length(means)
+  # Create data frame for ggplot
+  plot_data <- data.frame(
+    Treatment = factor(names(means), levels = names(means)),
+    Mean = means,
+    LowerCI = lowerbd,
+    UpperCI = upperbd
+  )
 
-  # Define some offsets for text placement
-  text_offset <- (max(upperbd) - min(lowerbd)) * 0.05
-
-  old_par <- par(no.readonly = TRUE) # Save current graphical parameters
-  on.exit(par(old_par)) # Restore graphical parameters on exit
-  par(mar = c(5, 4, 4, 3) + 0.1) # Adjust margins
-
-  # Plotting
-  plot(position, means, ylim = range(lowerbd - text_offset, upperbd + 3*text_offset),
-       xlim = range(0.5, length(means)+0.5),
-       pch = 19, xaxt = "n", # round down vs round up;
-       xlab = "Treatment Level", ylab = "Effect", main = "Treatment Effect Estimates", ...)
-  axis(1, at = position, labels = if (is_binomial) c("Comparator Level", "Reference Level", "RD", "RR", "OR") else c("Comparator Level", "Reference Level", "RD"))
-
-  # Error bars
-  arrows(position, lowerbd, position, upperbd, angle = 90, code = 3, length = 0.1, ...)
-
-  # Adding text for means and CIs
-  for (i in seq_along(means)) {
-    text(position[i], upperbd[i] + text_offset, labels = paste("Mean:", round(means[i], 2)), cex = 0.8, pos = 3)
-    text(position[i], upperbd[i] + 2 * text_offset, labels = paste("95% CI: [", round(lowerbd[i], 2), ", ", round(upperbd[i], 2), "]"), cex = 0.8, pos = 3)
-  }
-
-  # Check if the input is a model and extract treatment sequences if they exist
-  has_treatment_info <- "reference" %in% names(input) && "comparator" %in% names(input)
-
-  # Conditional treatment sequence information below x-axis labels
-  if (has_treatment_info) {
-    mtext(paste("(", paste(input$reference, collapse = ", "), ")", sep = ""), side = 1, at = position[2], line = 2)
-    mtext(paste("(", paste(input$comparator, collapse = ", "), ")", sep = ""), side = 1, at = position[1], line = 2)
-  }
+  # ggplot2 visualization
+  ggplot2::ggplot(plot_data, ggplot2::aes(x = Treatment, y = Mean)) +
+    ggplot2::geom_point(size = 3) +
+    ggplot2::geom_errorbar(ggplot2::aes(ymin = LowerCI, ymax = UpperCI), width = 0.2, color = "blue") +
+    ggplot2::labs(title = "Treatment Effect Estimates", x = "Treatment Level", y = "Effect") +
+    ggplot2::theme_minimal() +
+    ggplot2::geom_text(ggplot2::aes(y = UpperCI + 0.05, label = paste0("Mean: ", round(Mean, 2))), vjust = -0.5) +
+    ggplot2::geom_text(ggplot2::aes(y = UpperCI + 0.25, label = paste0("95% CI: [", round(LowerCI, 2), ", ", round(UpperCI, 2), "]")), vjust = -0.5) +
+    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) +
+    ggplot2::expand_limits(y = max(plot_data$UpperCI) + 0.15)
 }