diff --git a/DESCRIPTION b/DESCRIPTION
index f423982..866110a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -16,7 +16,7 @@ License: GPL-2
 LazyData: TRUE
 Encoding: UTF-8
 Depends:
-    ggplot2,
+    ggplot2(>= 3.4.0),
     ggpubr(>= 0.1.6)
 Imports:
     grid,
diff --git a/NAMESPACE b/NAMESPACE
index 66e0e14..3402bea 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -53,6 +53,8 @@ importFrom(grDevices,axisTicks)
 importFrom(magrittr,"%>%")
 importFrom(methods,is)
 importFrom(purrr,map)
+importFrom(rlang,"!!")
+importFrom(rlang,sym)
 importFrom(rlang,syms)
 importFrom(stats,anova)
 importFrom(stats,approx)
diff --git a/NEWS.md b/NEWS.md
index 2997f43..42721f4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,6 +7,7 @@
 ## Minor changes
 
 - R-ADDICT website is no longer live, so updating links in README (#622)
+- ggplot2 minimum version is now 3.4.0
 
 ## Bug fixes
 
diff --git a/R/ggcoxdiagnostics.R b/R/ggcoxdiagnostics.R
index 581f25a..62b64e7 100644
--- a/R/ggcoxdiagnostics.R
+++ b/R/ggcoxdiagnostics.R
@@ -99,19 +99,20 @@ ggcoxdiagnostics <- function (fit,
   else col_names <- names(stats::coef(fit))
   colnames(res) <- col_names
   res$xval <- xval
-  data2plot <- tidyr::gather_(res,
-                              key_col = "covariate", value_col = "res",
-                              gather_col = col_names)
+  data2plot <- tidyr::pivot_longer(
+    data = res, cols = dplyr::all_of(col_names),  
+    names_to = "covariate", values_to = "res"
+  )
 
   gplot <- ggplot(aes(xval, res), data = data2plot) +
            geom_point(col = point.col, shape = point.shape,
                        size = point.size, alpha = point.alpha)
 
   if (hline) gplot <- gplot + geom_hline(yintercept=hline.yintercept, col = hline.col,
-                                         size = hline.size, lty = hline.lty, alpha = hline.alpha)
+                                         linewidth = hline.size, lty = hline.lty, alpha = hline.alpha)
 
   if (sline) gplot <- gplot + geom_smooth(col = sline.col, se = sline.se, method = "loess",
-                                         size = sline.size, lty = sline.lty, alpha = sline.alpha)
+                                         linewidth = sline.size, lty = sline.lty, alpha = sline.alpha)
 
   gplot <- gplot + labs(x = xlabel, y = ylabel, title = title, subtitle = subtitle, caption = caption) + ggtheme
   # customization
diff --git a/R/ggsurvplot_combine.R b/R/ggsurvplot_combine.R
index 1ec7717..32aaca8 100644
--- a/R/ggsurvplot_combine.R
+++ b/R/ggsurvplot_combine.R
@@ -114,7 +114,7 @@ ggsurvplot_combine <- function(fit, data,
       all.levels <- c(all.levels, .levels(ss$strata))
     # convert strata into character before binding
     # avoid this warning: Unequal factor levels: coercing to character
-    grouped.d$survsummary <- map(grouped.d$survsummary,
+    grouped.d$survsummary <- purrr::map(grouped.d$survsummary,
                                  function(x){
                                    x$strata <- as.character(x$strata)
                                    x
@@ -161,9 +161,9 @@ ggsurvplot_combine <- function(fit, data,
       survtable$strata <- paste(fitname, "::", survtable$strata, sep = "") %>%
         factor(levels = strata.levels)
       survtable %>%
-        dplyr::select_( .dots = c("strata", "time", "n.risk", "pct.risk",
-                                  "n.event", "cum.n.event", "n.censor",
-                                  "cum.n.censor", "strata_size"))
+        dplyr::select( dplyr::all_of(c("strata", "time", "n.risk", "pct.risk",
+                                       "n.event", "cum.n.event", "n.censor",
+                                       "cum.n.censor", "strata_size")))
     }
 
     grouped.d <- grouped.d %>%
@@ -199,8 +199,8 @@ ggsurvplot_combine <- function(fit, data,
     #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
     if(keep.data){
-    res$data.survplot <- tibble::as.tibble(all.survsummary)
-    res$data.survtable <- tibble::as.tibble(all.survtable)
+    res$data.survplot <- tibble::as_tibble(all.survsummary)
+    res$data.survtable <- tibble::as_tibble(all.survtable)
     }
 
 
diff --git a/R/ggsurvplot_df.R b/R/ggsurvplot_df.R
index 1fb6a50..5f0eb64 100644
--- a/R/ggsurvplot_df.R
+++ b/R/ggsurvplot_df.R
@@ -176,7 +176,7 @@ ggsurvplot_df <- function(fit, fun = NULL,
   #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
   df[, .strata.var] <- factor( df[, .strata.var], levels = .levels(.strata), labels = legend.labs)
 
-  p <- ggplot2::ggplot(df, ggplot2::aes_string("time", "surv")) +
+  p <- ggplot2::ggplot(df, ggplot2::aes(x = !!sym("time"), y = !!sym("surv"))) +
     ggpubr::geom_exec(surv.geom, data = df, size = size, color = color, linetype = linetype, ...) +
     ggplot2::scale_y_continuous(breaks = y.breaks, labels = scale_labels, limits = ylim, expand = .expand) +
     ggplot2::coord_cartesian(xlim = xlim)+
diff --git a/R/ggsurvplot_facet.R b/R/ggsurvplot_facet.R
index e1f0f0d..4e00f82 100644
--- a/R/ggsurvplot_facet.R
+++ b/R/ggsurvplot_facet.R
@@ -179,10 +179,10 @@ ggsurvplot_facet <- function(fit, data, facet.by,
     pvalue <- surv_pvalue(grouped.d$fit, grouped.d$data, pval.coord = pval.coord,
                           pval.method.coord = pval.method.coord,...) %>%
       dplyr::bind_rows() %>%
-      tibble::as.tibble()
+      tibble::as_tibble()
     # Select the grouping variable columns and cbind the corresponding pvalue
     pvals.df <- grouped.d %>%
-      dplyr::select_( .dots = facet.by) %>%
+      dplyr::select(!!!syms(facet.by)) %>%
       dplyr::bind_cols(pvalue)
     pval.x <- pval.y <- pval.txt <- method.x <- method.y <- method <-  NULL
     p <- p +
diff --git a/R/surv_median.R b/R/surv_median.R
index 6ec63f6..f1caee4 100644
--- a/R/surv_median.R
+++ b/R/surv_median.R
@@ -59,7 +59,7 @@ surv_median <- function(fit, combine = FALSE){
     .table$strata <- rownames(.table)
 
     .table <- .table %>%
-      dplyr::select_(.dots = c("strata",  "median", "`0.95LCL`", "`0.95UCL`"))
+      dplyr::select(dplyr::all_of(c("strata",  "median", "0.95LCL", "0.95UCL")))
     colnames(.table) <- c("strata", "median", "lower", "upper")
     rownames(.table) <- NULL
     .table
diff --git a/R/utilities.R b/R/utilities.R
index 9049585..ebd5818 100644
--- a/R/utilities.R
+++ b/R/utilities.R
@@ -10,6 +10,7 @@
 #' @importFrom stats pchisq
 #' @importFrom survMisc ten comp
 #' @importFrom utils capture.output
+#' @importFrom rlang !! sym
 
 
 # Count the number of ggplots in a list
diff --git a/tests/testthat/setup-load_data.R b/tests/testthat/setup-load_data.R
new file mode 100644
index 0000000..175865c
--- /dev/null
+++ b/tests/testthat/setup-load_data.R
@@ -0,0 +1,2 @@
+library(survival)
+data('cancer')
diff --git a/tests/testthat/test-ggcoxdiagnostics.R b/tests/testthat/test-ggcoxdiagnostics.R
new file mode 100644
index 0000000..eb173f1
--- /dev/null
+++ b/tests/testthat/test-ggcoxdiagnostics.R
@@ -0,0 +1,16 @@
+test_that('ggcoxdiagnostics creates plot with all the observations', {
+  cph <- coxph(Surv(futime, fustat) ~ rx + age, data=ovarian)
+  p <- ggcoxdiagnostics(cph, type="deviance")
+  .build <- ggplot_build(p)
+  expect_equal(nrow(.build$data[[1]]), nrow(ovarian))
+})
+
+test_that('ggcoxdiagnostics with second type two rows for each observed event*term', {
+  cph <- coxph(Surv(futime, fustat) ~ rx + age, data=ovarian)
+  qty_terms <- length(attr(terms(cph$formula), "term.labels"))
+  qty_events <- sum(ovarian$fustat==1)
+  p <- ggcoxdiagnostics(cph, type="schoenfeld")
+  .build <- ggplot_build(p)
+  expect_equal(nrow(.build$data[[1]]), qty_terms*qty_events)
+})
+
diff --git a/tests/testthat/test-ggsurvplot_combine.R b/tests/testthat/test-ggsurvplot_combine.R
new file mode 100644
index 0000000..09e7c11
--- /dev/null
+++ b/tests/testthat/test-ggsurvplot_combine.R
@@ -0,0 +1,21 @@
+start_time <- 250
+fit1 <- survfit(Surv(time, status) ~ sex, data = lung)
+fit2 <- survfit(Surv(time, status) ~ sex, data = lung, start.time = start_time)
+
+test_that("survplot_combine plots successfully into 4 lines; second 2 fits have only (0,1) before start_time", {
+  p <- ggsurvplot_combine(list(
+      original=fit1, conditional=fit2
+    ), data = lung)
+  .build <- ggplot_build(p$plot)
+  .build_data <- .build$data[[1]]
+  expect_equal(length(unique(.build_data[['group']])), 4)
+  expect_lt(nrow(.build_data[(.build_data[['group']] >= 3) &
+                               (.build_data[['x']] < start_time), ]), 3)
+})
+
+test_that("survplot_combine includes dataframes when keep.data==TRUE", {
+  p <- ggsurvplot_combine(list(
+    original=fit1, conditional=fit2
+  ), data = lung, keep.data = TRUE)
+  expect_equal(length(names(p)), 3)
+})
diff --git a/tests/testthat/test-ggsurvplot_facet.R b/tests/testthat/test-ggsurvplot_facet.R
new file mode 100644
index 0000000..dd88862
--- /dev/null
+++ b/tests/testthat/test-ggsurvplot_facet.R
@@ -0,0 +1,21 @@
+test_that("ggsurvplot_facet creates the correct quanitty of subplots", {
+  fit <- survfit(Surv(time, status) ~ sex, data=kidney)
+  p <- ggsurvplot_facet(fit, kidney, facet.by='disease')
+  .build <- ggplot_build(p)
+  expect_equal(nrow(.build$data[[2]]),
+               length(unique(kidney[['disease']])))
+
+  fit <- survfit(Surv(time, status) ~ disease, data=kidney)
+  p <- ggsurvplot_facet(fit, kidney, facet.by='sex')
+  .build <- ggplot_build(p)
+  expect_equal(nrow(.build$data[[2]]),
+               length(unique(kidney[['sex']])))
+})
+
+test_that("ggsurvplot_facet calculates pvalue for each facet", {
+  fit <- survfit(Surv(time, status) ~ sex, data=kidney)
+  p <- ggsurvplot_facet(fit, kidney, facet.by='disease', pval = TRUE)
+  .build <- ggplot_build(p)
+  expect_equal(nrow(.build$plot$layer[[4]][['data']]),
+               length(unique(kidney[['disease']])))
+})
diff --git a/tests/testthat/test-ggsurvtable.R b/tests/testthat/test-ggsurvtable.R
index d7d6da9..d8566cd 100644
--- a/tests/testthat/test-ggsurvtable.R
+++ b/tests/testthat/test-ggsurvtable.R
@@ -1,11 +1,12 @@
 context("test-ggsurvtable")
 
-library(dplyr)
-library("survival")
-data("lung")
-fit <- survfit(Surv(time, status) ~ sex, data = lung)
+
 
 test_that("survtable y axis label colors work", {
+  library(dplyr)
+  library("survival")
+  #data("lung", package = "survival")
+  fit <- survfit(Surv(time, status) ~ sex, data = lung)
   p <- ggrisktable(fit, data = lung, color = "strata")
   .build <- ggplot_build(p)
   .build_data <- .build$data[[1]]
diff --git a/tests/testthat/test-surv_median.R b/tests/testthat/test-surv_median.R
new file mode 100644
index 0000000..560c310
--- /dev/null
+++ b/tests/testthat/test-surv_median.R
@@ -0,0 +1,6 @@
+test_that("surv_median calculates medians", {
+  df_medians <- surv_median(survfit(Surv(time, status) ~ sex,
+                                    data=lung))
+  expect_equal(df_medians[['median']], c(270,426))
+  expect_lt(df_medians[['lower']][2], 426)
+})