diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index a3ac6182..bf6b41de 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -5,6 +5,7 @@ on:
     branches: [main, master]
   pull_request:
     branches: [main, master]
+  workflow_dispatch:
 
 name: R-CMD-check
 
@@ -27,6 +28,7 @@ jobs:
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
       R_KEEP_PKG_SOURCE: yes
+      SIMTRIAL_TEST_BACKWARDS_COMPATIBILITY_REF: 341f77f0a598dc6d638bd5c48746952a7db88255
 
     steps:
       - uses: actions/checkout@v3
diff --git a/.gitignore b/.gitignore
index dc3d14c7..2b1dd008 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ inst/doc
 .Rhistory
 .RData
 .Ruserdata
+tests/testthat/fixtures/backwards-compatibility
diff --git a/DESCRIPTION b/DESCRIPTION
index 0180123c..64508da4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: simtrial
 Type: Package
 Title: Clinical Trial Simulation
-Version: 0.3.0
+Version: 0.3.0.1
 Authors@R: c(
     person("Keaven", "Anderson", email = "keaven_anderson@merck.com", role = c("aut")),
     person("Yilong", "Zhang", email = "elong0527@gmail.com", role =  c("aut")),
@@ -17,6 +17,7 @@ Authors@R: c(
     person("Heng", "Zhou", role = c("ctb")),
     person("Amin", "Shirazi", role = c("ctb")),
     person("Cole", "Manschot", role = c("ctb")),
+    person("John", "Blischak", role = c("ctb")),
     person("Merck & Co., Inc., Rahway, NJ, USA and its affiliates", role = "cph")
     )
 Description: simtrial provides some basic routines for simulating a
@@ -36,6 +37,7 @@ VignetteBuilder: knitr
 Depends: R (>= 3.5.0)
 Imports:
     Rcpp,
+    data.table,
     doFuture,
     dplyr,
     foreach,
@@ -43,9 +45,11 @@ Imports:
     magrittr,
     methods,
     mvtnorm,
+    stats,
     survival,
     tibble,
-    tidyr
+    tidyr,
+    utils
 Suggests:
     Matrix,
     bshazard,
@@ -54,6 +58,7 @@ Suggests:
     gsDesign,
     knitr,
     markdown,
+    remotes,
     rmarkdown,
     stringr,
     survMisc,
diff --git a/NAMESPACE b/NAMESPACE
index 83081e8d..c2070561 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,24 +18,26 @@ export(sim_pw_surv)
 export(simfix2simpwsurv)
 export(wlr)
 importFrom(Rcpp,sourceCpp)
+importFrom(data.table,":=")
+importFrom(data.table,.N)
+importFrom(data.table,as.data.table)
+importFrom(data.table,data.table)
+importFrom(data.table,frankv)
+importFrom(data.table,last)
+importFrom(data.table,merge.data.table)
+importFrom(data.table,rbindlist)
+importFrom(data.table,setDF)
+importFrom(data.table,setDT)
+importFrom(data.table,setorderv)
 importFrom(doFuture,"%dofuture%")
-importFrom(dplyr,arrange)
-importFrom(dplyr,desc)
 importFrom(dplyr,filter)
-importFrom(dplyr,first)
 importFrom(dplyr,full_join)
 importFrom(dplyr,group_by)
-importFrom(dplyr,lag)
-importFrom(dplyr,last)
-importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
-importFrom(dplyr,n)
 importFrom(dplyr,right_join)
-importFrom(dplyr,row_number)
 importFrom(dplyr,select)
 importFrom(dplyr,starts_with)
 importFrom(dplyr,summarize)
-importFrom(dplyr,ungroup)
 importFrom(future,plan)
 importFrom(magrittr,"%>%")
 importFrom(methods,is)
@@ -43,8 +45,6 @@ importFrom(mvtnorm,GenzBretz)
 importFrom(mvtnorm,pmvnorm)
 importFrom(survival,Surv)
 importFrom(survival,is.Surv)
-importFrom(tibble,as_tibble)
 importFrom(tibble,tibble)
-importFrom(tidyr,expand)
 importFrom(tidyr,replace_na)
 useDynLib(simtrial, .registration = TRUE)
diff --git a/R/counting_process.R b/R/counting_process.R
index daf570a6..0a970a86 100644
--- a/R/counting_process.R
+++ b/R/counting_process.R
@@ -38,7 +38,7 @@
 #'   treatment group value.
 #'
 #' @return
-#' A `tibble` grouped by `stratum` and sorted within stratum by `tte`.
+#' A data frame grouped by `stratum` and sorted within stratum by `tte`.
 #' Remain rows with at least one event in the population, at least one subject
 #' is at risk in both treatment group and control group.
 #' Other variables in this represent the following within each stratum at
@@ -57,8 +57,7 @@
 #'   hypothesis)
 #' - `var_o_minus_e`: Variance of `o_minus_e` under the same assumption.
 #'
-#' @importFrom dplyr group_by arrange desc mutate
-#'   summarize first filter select lag
+#' @importFrom data.table ":=" as.data.table setDF
 #'
 #' @export
 #'
@@ -95,40 +94,38 @@ counting_process <- function(x, arm) {
     stop("counting_process: event indicator must be 0 (censoring) or 1 (event).")
   }
 
-  ans <- x %>%
-    group_by(stratum) %>%
-    arrange(desc(tte)) %>%
-    mutate(
-      one = 1,
-      n_risk_tol = cumsum(one),
-      n_risk_trt = cumsum(treatment == arm)
-    ) %>%
-    # Handling ties using Breslow's method
-    group_by(stratum, mtte = desc(tte)) %>%
-    summarize(
-      events = sum(event),
-      n_event_tol = sum((treatment == arm) * event),
-      tte = first(tte),
-      n_risk_tol = max(n_risk_tol),
-      n_risk_trt = max(n_risk_trt)
-    ) %>%
-    # Keep calculation for observed time with at least one event,
-    # at least one subject is at risk in both treatment group and control group.
-    filter(events > 0, n_risk_tol - n_risk_trt > 0, n_risk_trt > 0) %>%
-    select(-mtte) %>%
-    mutate(s = 1 - events / n_risk_tol) %>%
-    arrange(stratum, tte) %>%
-    group_by(stratum) %>%
-    mutate(
-      # Left continuous Kaplan-Meier Estimator
-      s = dplyr::lag(cumprod(s), default = 1),
-      # Observed events minus Expected events in treatment group
-      o_minus_e = n_event_tol - n_risk_trt / n_risk_tol * events,
-      # Variance of o_minus_e
-      var_o_minus_e = (n_risk_tol - n_risk_trt) *
-        n_risk_trt * events * (n_risk_tol - events) /
-        n_risk_tol^2 / (n_risk_tol - 1)
-    )
+  ans <- as.data.table(x)
+  ans <- ans[order(tte, decreasing = TRUE), ]
+  ans[, one := 1]
+  ans[, `:=`(
+    n_risk_tol = cumsum(one),
+    n_risk_trt = cumsum(treatment == arm)
+  ), by = "stratum"]
 
-  ans
+  # Handling ties using Breslow's method
+  ans[, mtte := -tte]
+  ans <- ans[, .(
+    events = sum(event),
+    n_event_tol = sum((treatment == arm) * event),
+    tte = tte[1],
+    n_risk_tol = max(n_risk_tol),
+    n_risk_trt = max(n_risk_trt)
+  ), by = c("stratum", "mtte")]
+
+  # Keep calculation for observed time with at least one event,
+  # at least one subject is at risk in both treatment group and control group.
+  ans <- ans[events > 0 & n_risk_tol - n_risk_trt > 0 & n_risk_trt > 0, ]
+  ans[, mtte := NULL]
+  ans[, s := 1 - events / n_risk_tol]
+  ans <- ans[order(stratum, tte), ]
+  # Left continuous Kaplan-Meier Estimator
+  ans[, s := c(1, cumprod(s)[-length(s)]), by = "stratum"]
+  # Observed events minus Expected events in treatment group
+  ans[, o_minus_e := n_event_tol - n_risk_trt / n_risk_tol * events]
+  # Variance of o_minus_e
+  ans[, var_o_minus_e := (n_risk_tol - n_risk_trt) *
+    n_risk_trt * events * (n_risk_tol - events) /
+    n_risk_tol^2 / (n_risk_tol - 1)]
+
+  return(setDF(ans))
 }
diff --git a/R/cut_data_by_date.R b/R/cut_data_by_date.R
index cafdd0a4..75473759 100644
--- a/R/cut_data_by_date.R
+++ b/R/cut_data_by_date.R
@@ -24,7 +24,7 @@
 #'
 #' @return A dataset ready for survival analysis.
 #'
-#' @importFrom dplyr filter mutate select
+#' @importFrom data.table ":=" as.data.table setDF
 #'
 #' @export
 #'
@@ -33,13 +33,11 @@
 #' # cut at calendar time 5 after start of randomization
 #' sim_pw_surv(n = 20) %>% cut_data_by_date(5)
 cut_data_by_date <- function(x, cut_date) {
-  ans <- x %>%
-    filter(enroll_time <= cut_date) %>%
-    mutate(
-      tte = pmin(cte, cut_date) - enroll_time,
-      event = fail * (cte <= cut_date)
-    ) %>%
-    select(tte, event, stratum, treatment)
+  ans <- as.data.table(x)
+  ans <- ans[enroll_time <= cut_date, ]
+  ans[, tte := pmin(cte, cut_date) - enroll_time]
+  ans[, event := fail * (cte <= cut_date)]
+  ans <- ans[, c("tte", "event", "stratum", "treatment")]
 
-  ans
+  return(setDF(ans))
 }
diff --git a/R/cut_data_by_event.R b/R/cut_data_by_event.R
index fbdaf1d7..4578576d 100644
--- a/R/cut_data_by_event.R
+++ b/R/cut_data_by_event.R
@@ -24,9 +24,11 @@
 #' @param x A time-to-event dataset, for example, generated by [sim_pw_surv()].
 #' @param event Event count at which data cutoff is to be made.
 #'
-#' @return A tibble ready for survival analysis, including columns
+#' @return A data frame ready for survival analysis, including columns
 #'   time to event (`tte`), `event`, the `stratum`, and the `treatment`.
 #'
+#' @importFrom data.table setDF
+#'
 #' @export
 #'
 #' @examples
@@ -36,5 +38,5 @@
 cut_data_by_event <- function(x, event) {
   cut_date <- get_cut_date_by_event(x, event)
   ans <- x %>% cut_data_by_date(cut_date = cut_date)
-  ans
+  return(setDF(ans))
 }
diff --git a/R/get_cut_date_by_event.R b/R/get_cut_date_by_event.R
index 4bfed1f0..20817910 100644
--- a/R/get_cut_date_by_event.R
+++ b/R/get_cut_date_by_event.R
@@ -25,7 +25,7 @@
 #'   at which the targeted event count is reached, or if the final event count
 #'   is never reached, the final `cte` at which an event occurs.
 #'
-#' @importFrom dplyr ungroup select filter arrange mutate row_number last
+#' @importFrom data.table ":=" as.data.table frankv last
 #'
 #' @export
 #'
@@ -62,14 +62,12 @@
 #' y <- cut_data_by_date(x, cut_date = d)
 #' table(y$stratum, y$event)
 get_cut_date_by_event <- function(x, event) {
-  y <- x %>%
-    ungroup() %>%
-    select(cte, fail) %>%
-    filter(fail == 1) %>%
-    select(cte) %>%
-    arrange(cte) %>%
-    mutate(eventCount = row_number()) %>%
-    subset(eventCount <= event)
+  y <- as.data.table(x)
+  y <- y[fail == 1, ]
+  y <- y[, .(cte)]
+  y <- y[order(cte), ]
+  y[, eventCount := frankv(y, "cte", ties.method = "first")]
+  y <- y[eventCount <= event, ]
 
-  last(y$cte)
+  return(last(y$cte))
 }
diff --git a/R/global.R b/R/global.R
index 644abd80..218d9aa3 100644
--- a/R/global.R
+++ b/R/global.R
@@ -21,43 +21,46 @@
 
 utils::globalVariables(
   c(
-    "atrisk",
-    "Count",
+    ".",
+    "Ex1delayedEffect",
+    "N",
     "cte",
-    "dropoutRate",
-    "dropoutTime",
-    "enrollTime",
+    "dropout_rate",
+    "dropout_time",
+    "duration",
+    "enroll_time",
+    "event",
     "eventCount",
     "events",
-    "Ex1delayedEffect",
     "fail",
-    "fail_rate",
     "fail_time",
-    "hr",
-    "duration",
-    "enrollTime",
-    "event",
     "finish",
+    "hr",
+    "i",
     "lambda",
+    "max_weight",
     "mtte",
-    "N",
-    "OminusE",
+    "n_event_tol",
+    "n_risk_tol",
+    "n_risk_trt",
+    "nbrOfWorkers",
+    "o_minus_e",
     "one",
     "origin",
     "period",
     "rate",
     "s",
-    "S",
     "status",
     "stratum",
     "time",
     "treatment",
     "tte",
-    "txevents",
-    "Var",
-    "w",
-    "wOminusE",
-    "wVar",
-    "txatrisk"
+    "var_o_minus_e"
   )
 )
+
+# Workaround to remove `R CMD check` NOTE "All declared Imports should be used."
+# https://r-pkgs.org/dependencies-in-practice.html#how-to-not-use-a-package-in-imports
+ignore_unused_imports <- function() {
+  utils::globalVariables
+}
diff --git a/R/rpw_enroll.R b/R/rpw_enroll.R
index ccc549e9..c8c3fed1 100644
--- a/R/rpw_enroll.R
+++ b/R/rpw_enroll.R
@@ -35,9 +35,7 @@
 #'
 #' @return A vector of random enrollment times.
 #'
-#' @importFrom tibble tibble
-#' @importFrom dplyr mutate row_number group_by
-#' @importFrom tidyr expand
+#' @importFrom data.table ":=" .N as.data.table last
 #'
 #' @export
 #'
@@ -89,15 +87,12 @@ rpw_enroll <- function(
   }
 
   # Build `y` summarizes the start/end time, period order, etc.
-  y <- enroll_rate %>%
-    mutate(
-      period = row_number(),
-      finish = cumsum(duration),
-      lambda = duration * rate,
-      origin = dplyr::lag(finish, default = 0)
-    ) %>%
-    group_by(period) %>%
-    mutate(N = stats::rpois(n = 1, lambda = lambda))
+  y <- as.data.table(enroll_rate)
+  y[, period := seq_len(.N)]
+  y[, finish := cumsum(duration)]
+  y[, lambda := duration * rate]
+  y[, origin := c(0, finish[-length(finish)])]
+  y[, N := stats::rpois(n = .N, lambda = lambda)]
 
   # Deal with extreme cases where none randomized in fixed intervals
   if (sum(y$N) == 0) {
@@ -106,18 +101,18 @@ rpw_enroll <- function(
       return(ans)
     }
 
-    if (dplyr::last(enroll_rate$rate) <= 0) {
+    if (last(enroll_rate$rate) <= 0) {
       # Stop with error message if enrollment has not finished but enrollment rate for last period is less or equal with 0
       stop("rpw_enroll: please specify > 0 enrollment rate for the last period; otherwise enrollment cannot finish.")
     } else {
       # Otherwise, return inter-arrival exponential times
-      ans <- cumsum(stats::rexp(n = n, rate = dplyr::last(enroll_rate$rate))) + dplyr::last(y$finish)
+      ans <- cumsum(stats::rexp(n = n, rate = last(enroll_rate$rate))) + last(y$finish)
       return(ans)
     }
   }
 
   # Generate sorted uniform observations for Poisson count for each interval
-  z <- expand(y, enroll_time = sort(stats::runif(n = N, min = origin, max = finish)))
+  z <- y[, .(enroll_time = sort(stats::runif(n = N, min = origin, max = finish))), by = "period"]
 
   # If n not specified, return generated times
   if (is.null(n)) {
@@ -136,14 +131,14 @@ rpw_enroll <- function(
   n_add <- n - nrow(z)
   # Stop with error message if enrollment has not finished but
   # enrollment rate for last period is less or equal with 0
-  if (dplyr::last(enroll_rate$rate) <= 0) {
+  if (last(enroll_rate$rate) <= 0) {
     stop("rpw_enroll: please specify > 0 enrollment rate for the last period; otherwise enrollment cannot finish.")
   }
   # Otherwise, return inter-arrival exponential times
   else {
     ans <- c(
       z$enroll_time,
-      cumsum(stats::rexp(n_add, rate = dplyr::last(enroll_rate$rate))) + dplyr::last(y$finish)
+      cumsum(stats::rexp(n_add, rate = last(enroll_rate$rate))) + last(y$finish)
     )
     return(ans)
   }
diff --git a/R/sim_fixed_n.R b/R/sim_fixed_n.R
index 02400227..bbefff57 100644
--- a/R/sim_fixed_n.R
+++ b/R/sim_fixed_n.R
@@ -53,7 +53,7 @@
 #'   (2 and 3).
 #'
 #' @return
-#' A `tibble` including columns:
+#' A data frame including columns:
 #' - `event`: Event count.
 #' - `ln_hr`: Log-hazard ratio.
 #' - `z`: Normal test statistic; < 0 favors experimental.
@@ -67,6 +67,7 @@
 #' then columns `rho` and `gamma` are also included.
 #'
 #' @importFrom tibble tibble
+#' @importFrom data.table ":=" rbindlist setDF
 #' @importFrom doFuture "%dofuture%"
 #' @importFrom future plan
 #' @importFrom methods is
@@ -239,28 +240,6 @@ sim_fixed_n <- function(
 
   n_stratum <- nrow(stratum)
 
-  # Build a function to calculate z and log-hr
-  doAnalysis <- function(d, rho_gamma, n_stratum) {
-    if (nrow(rho_gamma) == 1) {
-      z <- tibble(z = (d %>% counting_process(arm = "experimental") %>% wlr(rho_gamma = rho_gamma))$z)
-    } else {
-      z <- d %>%
-        counting_process(arm = "experimental") %>%
-        wlr(rho_gamma = rho_gamma, return_corr = TRUE)
-    }
-
-    ans <- tibble(
-      event = sum(d$event),
-      ln_hr = ifelse(n_stratum > 1,
-        survival::coxph(survival::Surv(tte, event) ~ (treatment == "experimental") + survival::strata(stratum), data = d)$coefficients,
-        survival::coxph(survival::Surv(tte, event) ~ (treatment == "experimental"), data = d)$coefficients
-      ) %>% as.numeric()
-    )
-
-    ans <- cbind(ans, z)
-    return(ans)
-  }
-
   # Compute minimum planned follow-up time
   minFollow <- max(0, total_duration - sum(enroll_rate$duration))
 
@@ -299,7 +278,7 @@ sim_fixed_n <- function(
     )
 
     # Study date that targeted event rate achieved
-    tedate <- sim %>% get_cut_date_by_event(target_event)
+    tedate <- get_cut_date_by_event(sim, target_event)
 
     # Study data that targeted minimum follow-up achieved
     tmfdate <- max(sim$enroll_time) + minFollow
@@ -342,100 +321,91 @@ sim_fixed_n <- function(
 
     # Total duration cutoff
     if (tests[1]) {
-      d <- sim %>% cut_data_by_date(total_duration)
-      r1 <- d %>% doAnalysis(rho_gamma, n_stratum)
+      d <- cut_data_by_date(sim, total_duration)
+      r1 <- doAnalysis(d, rho_gamma, n_stratum)
     }
 
     # Targeted events cutoff
     if (tests[2]) {
-      d <- sim %>% cut_data_by_date(tedate)
-      r2 <- d %>% doAnalysis(rho_gamma, n_stratum)
+      d <- cut_data_by_date(sim, tedate)
+      r2 <- doAnalysis(d, rho_gamma, n_stratum)
     }
 
     # Minimum follow-up cutoff
     if (tests[3]) {
-      d <- sim %>% cut_data_by_date(tmfdate)
-      r3 <- d %>% doAnalysis(rho_gamma, n_stratum)
+      d <- cut_data_by_date(sim, tmfdate)
+      r3 <- doAnalysis(d, rho_gamma, n_stratum)
     }
 
-    addit <- NULL
+    addit <- list()
     # Planned duration cutoff
     if (1 %in% timing_type) {
-      addit <- rbind(
-        addit,
-        r1 %>% mutate(
-          cut = "Planned duration",
-          duration = total_duration
-        )
-      )
+      rtemp <- cbind(r1, cut = "Planned duration", duration = total_duration)
+      addit <- c(addit, list(rtemp))
     }
 
     # Targeted events cutoff
     if (2 %in% timing_type) {
-      addit <- rbind(
-        addit,
-        r2 %>% mutate(
-          cut = "Targeted events",
-          duration = tedate
-        )
-      )
+      rtemp <- cbind(r2, cut = "Targeted events", duration = tedate)
+      addit <- c(addit, list(rtemp))
     }
 
     # Minimum follow-up duration target
     if (3 %in% timing_type) {
-      addit <- rbind(
-        addit,
-        r3 %>% mutate(
-          cut = "Minimum follow-up",
-          duration = tmfdate
-        )
-      )
+      rtemp <- cbind(r3, cut = "Minimum follow-up", duration = tmfdate)
+      addit <- c(addit, list(rtemp))
     }
 
     # Max of planned duration, targeted events
     if (4 %in% timing_type) {
       if (tedate > total_duration) {
-        addit <- rbind(
-          addit,
-          r2 %>% mutate(
-            cut = "Max(planned duration, event cut)",
-            duration = tedate
-          )
-        )
+        rtemp <- cbind(r2, cut = "Max(planned duration, event cut)", duration = tedate)
+        addit <- c(addit, list(rtemp))
       } else {
-        addit <- rbind(
-          addit,
-          r1 %>% mutate(
-            cut = "Max(planned duration, event cut)",
-            duration = total_duration
-          )
-        )
+        rtemp <- cbind(r1, cut = "Max(planned duration, event cut)", duration = total_duration)
+        addit <- c(addit, list(rtemp))
       }
     }
 
     # Max of minimum follow-up, targeted events
     if (5 %in% timing_type) {
       if (tedate > tmfdate) {
-        addit <- rbind(
-          addit,
-          r2 %>% mutate(
-            cut = "Max(min follow-up, event cut)",
-            duration = tedate
-          )
-        )
+        rtemp <- cbind(r2, cut = "Max(min follow-up, event cut)", duration = tedate)
+        addit <- c(addit, list(rtemp))
       } else {
-        addit <- rbind(
-          addit,
-          r3 %>% mutate(
-            cut = "Max(min follow-up, event cut)",
-            duration = tmfdate
-          )
-        )
+        rtemp <- cbind(r3, cut = "Max(min follow-up, event cut)", duration = tmfdate)
+        addit <- c(addit, list(rtemp))
       }
     }
 
-    addit %>% mutate(sim = i)
+    results_sim <- rbindlist(addit)
+    results_sim[, sim := i]
+    setDF(results_sim)
+    return(results_sim)
   }
 
   return(results)
 }
+
+# Build a function to calculate z and log-hr
+doAnalysis <- function(d, rho_gamma, n_stratum) {
+  if (nrow(rho_gamma) == 1) {
+    tmp <- counting_process(d, arm = "experimental")
+    tmp <- wlr(tmp, rho_gamma = rho_gamma)
+    z <- data.frame(z = tmp$z)
+  } else {
+    tmp <- counting_process(d, arm = "experimental")
+    z <- wlr(tmp, rho_gamma = rho_gamma, return_corr = TRUE)
+  }
+
+  event <- sum(d$event)
+  if (n_stratum > 1) {
+    ln_hr <- survival::coxph(survival::Surv(tte, event) ~ (treatment == "experimental") + survival::strata(stratum), data = d)$coefficients
+  } else {
+    ln_hr <- survival::coxph(survival::Surv(tte, event) ~ (treatment == "experimental"), data = d)$coefficients
+  }
+  ln_hr <- as.numeric(ln_hr)
+
+  ans <- cbind(event, ln_hr, z)
+  return(ans)
+}
diff --git a/R/sim_pw_surv.R b/R/sim_pw_surv.R
index 573a7667..5482ad0c 100644
--- a/R/sim_pw_surv.R
+++ b/R/sim_pw_surv.R
@@ -44,7 +44,7 @@
 #' @param dropout_rate Dropout rates; see details and examples;
 #'   note that treatments need to be the same as input in block.
 #'
-#' @return A `tibble` with the following variables for each observation:
+#' @return A data frame with the following variables for each observation:
 #' - `stratum`.
 #' - `enroll_time`: Enrollment time for the observation.
 #' - `Treatment`: Treatment group; this will be one of the values
@@ -56,6 +56,7 @@
 #' - `fail`: Indicator that `cte` was set using failure time;
 #'   i.e., 1 is a failure, 0 is a dropout.
 #'
+#' @importFrom data.table ":=" .N data.table setDF setDT setorderv
 #' @importFrom dplyr group_by mutate
 #' @importFrom tibble tibble
 #'
@@ -142,24 +143,29 @@ sim_pw_surv <- function(
       duration = rep(100, 2),
       rate = rep(.001, 2)
     )) {
-  # Start tibble by generating stratum and enrollment times
-  x <- tibble(stratum = sample(
+  # Start table by generating stratum and enrollment times
+  x <- data.table(stratum = sample(
     x = stratum$stratum,
     size = n,
     replace = TRUE,
     prob = stratum$p
-  )) %>%
-    mutate(enroll_time = rpw_enroll(n, enroll_rate)) %>%
-    group_by(stratum) %>%
-    # Assign treatment
-    mutate(treatment = randomize_by_fixed_block(n = n(), block = block)) %>%
-    # Generate time to failure and time to dropout
-    group_by(stratum, treatment)
+  ))
+  x[, enroll_time := rpw_enroll(n, enroll_rate)]
+  # The awkward back and forth ordering is to maintain 1:1 parity with
+  # dplyr::group_by() for backwards compatibility. group_by() sorts by the
+  # grouping variable and then returns the rows to their original positions.
+  # This is mainly for testing for backwards compatibility. Since the
+  # treatments are assigned randomly by group, it would still be statistically
+  # valid without this ordering
+  setorderv(x, "stratum")
+  x[, treatment := randomize_by_fixed_block(n = .N, block = block), by = "stratum"]
+  setorderv(x, "enroll_time")
 
+  # Generate time to failure and time to dropout
   unique_stratum <- unique(x$stratum)
   unique_treatment <- unique(x$treatment)
-  x$fail_time <- 0
-  x$dropout_time <- 0
+  x[, fail_time := 0]
+  x[, dropout_time := 0]
 
   for (sr in unique_stratum) {
     for (tr in unique_treatment) {
@@ -176,11 +182,9 @@ sim_pw_surv <- function(
   }
 
   # Set calendar time-to-event and failure indicator
-  ans <- x %>%
-    mutate(
-      cte = pmin(dropout_time, fail_time) + enroll_time,
-      fail = (fail_time <= dropout_time) * 1
-    )
+  ans <- setDT(x)
+  ans[, cte := pmin(dropout_time, fail_time) + enroll_time]
+  ans[, fail := (fail_time <= dropout_time) * 1]
 
-  ans
+  return(setDF(ans))
 }
diff --git a/R/simfix2simPWSurv.R b/R/simfix2simPWSurv.R
index d8da7c9f..f6fe0300 100644
--- a/R/simfix2simPWSurv.R
+++ b/R/simfix2simPWSurv.R
@@ -33,11 +33,11 @@
 #'   hazard ratio for experimental vs. control,
 #'   and dropout rates by stratum and time period.
 #'
-#' @return A list of two `tibble` components formatted for
+#' @return A list of two data frame components formatted for
 #'   [sim_pw_surv()]: `fail_rate` and `dropout_rate`.
 #'
 #' @importFrom tibble tibble
-#' @importFrom dplyr group_by mutate n ungroup select
+#' @importFrom data.table ":=" .N as.data.table setDF
 #'
 #' @export
 #'
@@ -92,40 +92,36 @@ simfix2simpwsurv <- function(
       dropout_rate = rep(.001, 2)
     )) {
   # Put failure rates into sim_pw_surv format
-  fr <- rbind(
-    fail_rate %>%
-      group_by(stratum) %>%
-      mutate(
-        treatment = "control",
-        rate = fail_rate, period = 1:n()
-      ) %>%
-      ungroup(),
-    fail_rate %>%
-      group_by(stratum) %>%
-      mutate(
-        treatment = "experimental",
-        rate = fail_rate * hr,
-        period = 1:n()
-      ) %>%
-      ungroup()
-  ) %>%
-    select("stratum", "period", "treatment", "duration", "rate")
+  fr_control <- as.data.table(fail_rate)
+  fr_control[, `:=`(
+    treatment = "control",
+    rate = fail_rate,
+    period = seq_len(.N)
+  ), by = "stratum"]
+
+  fr_experimental <- as.data.table(fail_rate)
+  fr_experimental[, `:=`(
+    treatment = "experimental",
+    rate = fail_rate * hr,
+    period = seq_len(.N)
+  ), by = "stratum"]
+
+  fr <- rbind(fr_control, fr_experimental)
+  fr <- fr[, c("stratum", "period", "treatment", "duration", "rate")]
 
   # Put dropout rates into sim_pw_surv format
-  dr <- fail_rate %>%
-    group_by(stratum) %>%
-    mutate(
-      treatment = "control",
-      rate = dropout_rate,
-      period = 1:n()
-    ) %>%
-    select("stratum", "period", "treatment", "duration", "rate") %>%
-    ungroup()
+  dr_control <- as.data.table(fail_rate)
+  dr_control[, `:=`(
+    treatment = "control",
+    rate = dropout_rate,
+    period = seq_len(.N)
+  ), by = "stratum"]
+
+  dr_experimental <- dr_control
+  dr_experimental$treatment <- "experimental"
 
-  dr <- rbind(
-    dr,
-    dr %>% mutate(treatment = "experimental")
-  )
+  dr <- rbind(dr_control, dr_experimental)
+  dr <- dr[, c("stratum", "period", "treatment", "duration", "rate")]
 
-  list(fail_rate = fr, dropout_rate = dr)
+  list(fail_rate = setDF(fr), dropout_rate = setDF(dr))
 }
diff --git a/R/wlr.R b/R/wlr.R
index a37eea3c..7202a6ac 100644
--- a/R/wlr.R
+++ b/R/wlr.R
@@ -31,7 +31,7 @@
 #'   estimated correlation for weighted sum of observed minus expected;
 #'   see details; Default: `FALSE`.
 #'
-#' @return A `tibble` with `rho_gamma` as input and the FH test statistic
+#' @return A data frame with `rho_gamma` as input and the FH test statistic
 #'   for the data in `x`. (`z`, a directional square root of the usual
 #'   weighted logrank test); if variance calculations are specified
 #'   (for example, to be used for covariances in a combination test),
@@ -78,8 +78,7 @@
 #'   The stratified Fleming-Harrington weighted logrank test is then computed as:
 #'   \deqn{z = \sum_i X_i/\sqrt{\sum_i V_i}.}
 #'
-#' @importFrom tibble tibble as_tibble
-#' @importFrom dplyr left_join select
+#' @importFrom data.table data.table merge.data.table setDF
 #'
 #' @export
 #'
@@ -178,43 +177,6 @@ wlr <- function(
     stop("wlr: can't report the correlation for a single WLR test.")
   }
 
-  # Build an internal function to compute the Z statistics
-  # under a sequence of rho and gamma of WLR.
-  wlr_z_stat <- function(x, rho_gamma, return_variance) {
-    ans <- rho_gamma
-
-    xx <- x %>%
-      ungroup() %>%
-      select(s, o_minus_e, var_o_minus_e)
-
-    ans$z <- rep(0, nrow(rho_gamma))
-
-    if (return_variance) {
-      ans$var <- rep(0, nrow(rho_gamma))
-    }
-
-    for (i in 1:nrow(rho_gamma)) {
-      y <- xx %>%
-        mutate(
-          weight = s^rho_gamma$rho[i] * (1 - s)^rho_gamma$gamma[i],
-          weighted_o_minus_e = weight * o_minus_e,
-          weighted_var = weight^2 * var_o_minus_e
-        ) %>%
-        summarize(
-          weighted_var = sum(weighted_var),
-          weighted_o_minus_e = sum(weighted_o_minus_e)
-        )
-
-      ans$z[i] <- y$weighted_o_minus_e / sqrt(y$weighted_var)
-
-      if (return_variance) {
-        ans$var[i] <- y$weighted_var
-      }
-    }
-
-    ans
-  }
-
   if (n_weight == 1) {
     ans <- wlr_z_stat(x, rho_gamma = rho_gamma, return_variance = return_variance)
   } else {
@@ -228,16 +190,19 @@ wlr <- function(
       matrix(rho_gamma$gamma, nrow = n_weight, ncol = n_weight, byrow = TRUE)
     ) / 2
 
-    # Convert back to tibble
-    rg_new <- tibble(rho = as.numeric(ave_rho), gamma = as.numeric(ave_gamma))
+    # Convert to data.table
+    rg_new <- data.table(rho = as.numeric(ave_rho), gamma = as.numeric(ave_gamma))
     # Get unique values of rho, gamma
-    rg_unique <- rg_new %>% unique()
+    rg_unique <- unique(rg_new)
 
     # Compute FH statistic for unique values
     # and merge back to full set of pairs
-    rg_fh <- rg_new %>% left_join(
-      wlr_z_stat(x, rho_gamma = rg_unique, return_variance = TRUE),
-      by = c("rho" = "rho", "gamma" = "gamma")
+    rg_fh <- merge.data.table(
+      x = rg_new,
+      y = wlr_z_stat(x, rho_gamma = rg_unique, return_variance = TRUE),
+      by = c("rho", "gamma"),
+      all.x = TRUE,
+      sort = FALSE
     )
 
     # Get z statistics for input rho, gamma combinations
@@ -249,16 +214,48 @@ wlr <- function(
     if (return_corr) {
       corr_mat <- stats::cov2cor(cov_mat)
       colnames(corr_mat) <- paste("v", 1:ncol(corr_mat), sep = "")
-      ans <- cbind(rho_gamma, z, as_tibble(corr_mat))
+      ans <- cbind(rho_gamma, z, as.data.frame(corr_mat))
     } else if (return_variance) {
       corr_mat <- cov_mat
       colnames(corr_mat) <- paste("v", 1:ncol(corr_mat), sep = "")
-      ans <- cbind(rho_gamma, z, as_tibble(corr_mat))
+      ans <- cbind(rho_gamma, z, as.data.frame(corr_mat))
     } else if (return_corr + return_corr == 0) {
       corr_mat <- NULL
       ans <- cbind(rho_gamma, z)
     }
   }
 
-  ans
+  return(setDF(ans))
+}
+
+# Build an internal function to compute the Z statistics
+# under a sequence of rho and gamma of WLR.
+wlr_z_stat <- function(x, rho_gamma, return_variance) {
+  ans <- rho_gamma
+
+  xx <- x[, c("s", "o_minus_e", "var_o_minus_e")]
+
+  ans$z <- rep(0, nrow(rho_gamma))
+
+  if (return_variance) {
+    ans$var <- rep(0, nrow(rho_gamma))
+  }
+
+  for (i in 1:nrow(rho_gamma)) {
+    weight <- xx$s^rho_gamma$rho[i] * (1 - xx$s)^rho_gamma$gamma[i]
+    weighted_o_minus_e <- weight * xx$o_minus_e
+    weighted_var <- weight^2 * xx$var_o_minus_e
+
+    weighted_var_total <- sum(weighted_var)
+    weighted_o_minus_e_total <- sum(weighted_o_minus_e)
+
+
+    ans$z[i] <- weighted_o_minus_e_total / sqrt(weighted_var_total)
+
+    if (return_variance) {
+      ans$var[i] <- weighted_var_total
+    }
+  }
+
+  return(ans)
 }
diff --git a/man/counting_process.Rd b/man/counting_process.Rd
index 9ecc9454..1ec52bfb 100644
--- a/man/counting_process.Rd
+++ b/man/counting_process.Rd
@@ -20,7 +20,7 @@ counting_process(x, arm)
 treatment group value.}
 }
 \value{
-A \code{tibble} grouped by \code{stratum} and sorted within stratum by \code{tte}.
+A data frame grouped by \code{stratum} and sorted within stratum by \code{tte}.
 Remain rows with at least one event in the population, at least one subject
 is at risk in both treatment group and control group.
 Other variables in this represent the following within each stratum at
diff --git a/man/cut_data_by_event.Rd b/man/cut_data_by_event.Rd
index b0554810..23b4a25e 100644
--- a/man/cut_data_by_event.Rd
+++ b/man/cut_data_by_event.Rd
@@ -12,7 +12,7 @@ cut_data_by_event(x, event)
 \item{event}{Event count at which data cutoff is to be made.}
 }
 \value{
-A tibble ready for survival analysis, including columns
+A data frame ready for survival analysis, including columns
 time to event (\code{tte}), \code{event}, the \code{stratum}, and the \code{treatment}.
 }
 \description{
diff --git a/man/sim_fixed_n.Rd b/man/sim_fixed_n.Rd
index eac14e8b..87961cc7 100644
--- a/man/sim_fixed_n.Rd
+++ b/man/sim_fixed_n.Rd
@@ -48,7 +48,7 @@ see details. Default is to include all available cutoff methods.}
 to specify one Fleming-Harrington weighted logrank test per row.}
 }
 \value{
-A \code{tibble} including columns:
+A data frame including columns:
 \itemize{
 \item \code{event}: Event count.
 \item \code{ln_hr}: Log-hazard ratio.
@@ -128,7 +128,7 @@ mean(p < .025)
 # Example 3
 # Use two cores
 set.seed(2023)
-plan("multisession", workers=2)
+plan("multisession", workers = 2)
 sim_fixed_n(n_sim = 10)
 plan("sequential")
 }
diff --git a/man/sim_pw_surv.Rd b/man/sim_pw_surv.Rd
index e835c8a0..1143010b 100644
--- a/man/sim_pw_surv.Rd
+++ b/man/sim_pw_surv.Rd
@@ -34,7 +34,7 @@ note that treatments need to be the same as input in block.}
 note that treatments need to be the same as input in block.}
 }
 \value{
-A \code{tibble} with the following variables for each observation:
+A data frame with the following variables for each observation:
 \itemize{
 \item \code{stratum}.
 \item \code{enroll_time}: Enrollment time for the observation.
diff --git a/man/simfix2simpwsurv.Rd b/man/simfix2simpwsurv.Rd
index d17a306d..698ae2d6 100644
--- a/man/simfix2simpwsurv.Rd
+++ b/man/simfix2simpwsurv.Rd
@@ -16,7 +16,7 @@ hazard ratio for experimental vs. control,
 and dropout rates by stratum and time period.}
 }
 \value{
-A list of two \code{tibble} components formatted for
+A list of two data frame components formatted for
 \code{\link[=sim_pw_surv]{sim_pw_surv()}}: \code{fail_rate} and \code{dropout_rate}.
 }
 \description{
diff --git a/man/simtrial-package.Rd b/man/simtrial-package.Rd
index 789f7ce4..c8bd6e1a 100644
--- a/man/simtrial-package.Rd
+++ b/man/simtrial-package.Rd
@@ -39,6 +39,7 @@ Other contributors:
   \item Heng Zhou [contributor]
   \item Amin Shirazi [contributor]
   \item Cole Manschot [contributor]
+  \item John Blischak [contributor]
   \item Merck & Co., Inc., Rahway, NJ, USA and its affiliates [copyright holder]
 }
 
diff --git a/man/wlr.Rd b/man/wlr.Rd
index d5166109..1ca05a91 100644
--- a/man/wlr.Rd
+++ b/man/wlr.Rd
@@ -28,7 +28,7 @@ estimated correlation for weighted sum of observed minus expected;
 see details; Default: \code{FALSE}.}
 }
 \value{
-A \code{tibble} with \code{rho_gamma} as input and the FH test statistic
+A data frame with \code{rho_gamma} as input and the FH test statistic
 for the data in \code{x}. (\code{z}, a directional square root of the usual
 weighted logrank test); if variance calculations are specified
 (for example, to be used for covariances in a combination test),
diff --git a/tests/testthat/generate-backwards-compatible-test-data.R b/tests/testthat/generate-backwards-compatible-test-data.R
new file mode 100644
index 00000000..8499e217
--- /dev/null
+++ b/tests/testthat/generate-backwards-compatible-test-data.R
@@ -0,0 +1,272 @@
+# Generates values from a previous version of simtrial for testing backwards
+# compatibility
+
+# Generate backwards compatible test data
+#
+# Used by tests/testthat/test-backwards-compatibility.R
+#
+# reference - the Git commit ID or tag name to pass to the arg `ref` of install_github()
+#
+# outdir - directory to save the output RDS files
+generate_test_data <- function(
+    reference = Sys.getenv(
+      x = "SIMTRIAL_TEST_BACKWARDS_COMPATIBILITY_REF",
+      unset = "341f77f0a598dc6d638bd5c48746952a7db88255"
+    ),
+    outdir = "tests/testthat/fixtures/backwards-compatibility") {
+  # Setup ----------------------------------------------------------------------
+
+  if (nchar(reference) < 1) {
+    stop("Invalid Git reference to pass to remotes::install_github()")
+  }
+  message("Installing simtrial from reference ", reference)
+
+  tmplib <- tempfile()
+  dir.create(tmplib)
+  old_lib_paths <- .libPaths()
+  .libPaths(c(tmplib, old_lib_paths))
+  on.exit(.libPaths(old_lib_paths), add = TRUE)
+  on.exit(unlink(tmplib, recursive = TRUE), add = TRUE)
+
+  remotes::install_github(
+    repo = "Merck/simtrial",
+    ref = reference,
+    dependencies = FALSE,
+    upgrade = FALSE
+  )
+  library("simtrial")
+  packageVersion("simtrial")
+  on.exit(detach("package:simtrial"), add = TRUE)
+
+  message("Saving output RDS files to ", outdir)
+  dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
+
+  # sim_fixed_n() --------------------------------------------------------------
+
+  set.seed(12345)
+  ex1 <- sim_fixed_n(n_sim = 2)
+  saveRDS(ex1, file.path(outdir, "sim_fixed_n_ex1.rds"))
+  set.seed(12345)
+  ex2 <- sim_fixed_n(n_sim = 2, rho_gamma = data.frame(rho = 0, gamma = c(0, 1)))
+  saveRDS(ex2, file.path(outdir, "sim_fixed_n_ex2.rds"))
+  set.seed(12345)
+  ex3 <- sim_fixed_n(
+    n_sim = 2,
+    timing_type = c(2, 5),
+    rho_gamma = data.frame(rho = 0, gamma = c(0, 1))
+  )
+  saveRDS(ex3, file.path(outdir, "sim_fixed_n_ex3.rds"))
+
+  # cut_data_by_date() ---------------------------------------------------------
+
+  set.seed(12345)
+  ex1 <- cut_data_by_date(
+    x = sim_pw_surv(n = 20),
+    cut_date = 5
+  )
+  saveRDS(ex1, file.path(outdir, "cut_data_by_date_ex1.rds"))
+
+  # get_cut_date_by_event() ----------------------------------------------------
+
+  set.seed(12345)
+  x <- sim_pw_surv(
+    n = 200,
+    stratum = data.frame(
+      stratum = c("Positive", "Negative"),
+      p = c(.5, .5)
+    ),
+    fail_rate = data.frame(
+      stratum = rep(c("Positive", "Negative"), 2),
+      period = rep(1, 4),
+      treatment = c(rep("control", 2), rep("experimental", 2)),
+      duration = rep(1, 4),
+      rate = log(2) / c(6, 9, 9, 12)
+    ),
+    dropout_rate = data.frame(
+      stratum = rep(c("Positive", "Negative"), 2),
+      period = rep(1, 4),
+      treatment = c(rep("control", 2), rep("experimental", 2)),
+      duration = rep(1, 4),
+      rate = rep(.001, 4)
+    )
+  )
+  ex1 <- get_cut_date_by_event(subset(x, stratum == "Positive"), event = 50)
+  saveRDS(ex1, file.path(outdir, "get_cut_date_by_event_ex1.rds"))
+
+  # counting_process() ---------------------------------------------------------
+
+  # Example 1
+  x <- data.frame(
+    stratum = c(rep(1, 10), rep(2, 6)),
+    treatment = rep(c(1, 1, 0, 0), 4),
+    tte = 1:16,
+    event = rep(c(0, 1), 8)
+  )
+  ex1 <- counting_process(x, arm = 1)
+  saveRDS(ex1, file.path(outdir, "counting_process_ex1.rds"))
+
+  # Example 2
+  set.seed(12345)
+  x <- sim_pw_surv(n = 400)
+  y <- cut_data_by_event(x, 150)
+  ex2 <- counting_process(y, arm = "experimental")
+  saveRDS(ex2, file.path(outdir, "counting_process_ex2.rds"))
+
+  # Example 3
+  # Counting Process Format with ties
+  x <- data.frame(
+    stratum = c(rep(1, 10), rep(2, 6)),
+    treatment = rep(c(1, 1, 0, 0), 4),
+    tte = c(rep(1:4, each = 4)),
+    event = rep(c(0, 1), 8)
+  )
+  arm <- 1
+  ex3 <- counting_process(x, arm)
+  saveRDS(ex3, file.path(outdir, "counting_process_ex3.rds"))
+
+  # wlr() ----------------------------------------------------------------------
+
+  # Example 1
+  # Use default enrollment and event rates at cut at 100 events
+  set.seed(12345)
+  x <- sim_pw_surv(n = 200)
+  x <- cut_data_by_event(x, 100)
+  x <- counting_process(x, arm = "experimental")
+
+  # Compute the corvariance between FH(0, 0), FH(0, 1) and FH(1, 0)
+  ex1 <- wlr(x, rho_gamma = data.frame(rho = c(0, 0, 1), gamma = c(0, 1, 0)))
+  saveRDS(ex1, file.path(outdir, "wlr_ex1.rds"))
+  ex1_var <- wlr(x, rho_gamma = data.frame(rho = c(0, 0, 1), gamma = c(0, 1, 0)), return_variance = TRUE)
+  saveRDS(ex1_var, file.path(outdir, "wlr_ex1_var.rds"))
+  ex1_cor <- wlr(x, rho_gamma = data.frame(rho = c(0, 0, 1), gamma = c(0, 1, 0)), return_corr = TRUE)
+  saveRDS(ex1_cor, file.path(outdir, "wlr_ex1_cor.rds"))
+
+  # Example 2
+  # Use default enrollment and event rates at cut of 100 events
+  set.seed(12345)
+  x <- sim_pw_surv(n = 200)
+  x <- cut_data_by_event(x, 100)
+  x <- counting_process(x, arm = "experimental")
+  ex2 <- wlr(x, rho_gamma = data.frame(rho = c(0, 0), gamma = c(0, 1)), return_corr = TRUE)
+  saveRDS(ex2, file.path(outdir, "wlr_ex2.rds"))
+
+  # rpw_enroll() ---------------------------------------------------------------
+
+  # Example 1
+  # Piecewise uniform (piecewise exponential inter-arrival times) for 10k patients
+  # enrollment Enrollment rates of 5 for time 0-100, 15 for 100-300, and 30
+  # thereafter
+  set.seed(12345)
+  ex1 <- rpw_enroll(
+    n = 1e5,
+    enroll_rate = data.frame(
+      rate = c(5, 15, 30),
+      duration = c(100, 200, 100)
+    )
+  )
+  saveRDS(ex1, file.path(outdir, "rpw_enroll_ex1.rds"))
+
+  # Example 2
+  # Exponential enrollment
+  set.seed(12345)
+  ex2 <- rpw_enroll(
+    n = 1e5,
+    enroll_rate = data.frame(rate = .03, duration = 1)
+  )
+  saveRDS(ex2, file.path(outdir, "rpw_enroll_ex2.rds"))
+
+  # simfix2simpwsurv() ---------------------------------------------------------
+
+  # Example 1
+  # Convert standard input
+  ex1 <- simfix2simpwsurv()
+  saveRDS(ex1, file.path(outdir, "simfix2simpwsurv_ex1.rds"))
+
+  # Example 2
+  # Stratified example
+  fail_rate <- data.frame(
+    stratum = c(rep("Low", 3), rep("High", 3)),
+    duration = rep(c(4, 10, 100), 2),
+    fail_rate = c(
+      .04, .1, .06,
+      .08, .16, .12
+    ),
+    hr = c(
+      1.5, .5, 2 / 3,
+      2, 10 / 16, 10 / 12
+    ),
+    dropout_rate = .01
+  )
+  ex2 <- simfix2simpwsurv(fail_rate)
+  saveRDS(ex2, file.path(outdir, "simfix2simpwsurv_ex2.rds"))
+
+  # sim_pw_surv() --------------------------------------------------------------
+
+  # Example 1
+  set.seed(12345)
+  ex1 <- sim_pw_surv(n = 20)
+  saveRDS(ex1, file.path(outdir, "sim_pw_surv_ex1.rds"))
+
+  # Example 2
+  # 3:1 randomization
+  set.seed(12345)
+  ex2 <- sim_pw_surv(
+    n = 20,
+    block = c(rep("experimental", 3), "control")
+  )
+  saveRDS(ex2, file.path(outdir, "sim_pw_surv_ex2.rds"))
+
+  # Example 3
+  # Simulate 2 stratum; will use defaults for blocking and enrollRates
+  set.seed(12345)
+  ex3 <- sim_pw_surv(
+    n = 20,
+    # 2 stratum,30% and 70% prevalence
+    stratum = data.frame(stratum = c("Low", "High"), p = c(.3, .7)),
+    fail_rate = data.frame(
+      stratum = c(rep("Low", 4), rep("High", 4)),
+      period = rep(1:2, 4),
+      treatment = rep(c(
+        rep("control", 2),
+        rep("experimental", 2)
+      ), 2),
+      duration = rep(c(3, 1), 4),
+      rate = c(.03, .05, .03, .03, .05, .08, .07, .04)
+    ),
+    dropout_rate = data.frame(
+      stratum = c(rep("Low", 2), rep("High", 2)),
+      period = rep(1, 4),
+      treatment = rep(c("control", "experimental"), 2),
+      duration = rep(1, 4),
+      rate = rep(.001, 4)
+    )
+  )
+  saveRDS(ex3, file.path(outdir, "sim_pw_surv_ex3.rds"))
+
+  # Example 4
+  # If you want a more rectangular entry for a tibble
+  fail_rate <- list(
+    data.frame(stratum = "Low", period = 1, treatment = "control", duration = 3, rate = .03),
+    data.frame(stratum = "Low", period = 1, treatment = "experimental", duration = 3, rate = .03),
+    data.frame(stratum = "Low", period = 2, treatment = "experimental", duration = 3, rate = .02),
+    data.frame(stratum = "High", period = 1, treatment = "control", duration = 3, rate = .05),
+    data.frame(stratum = "High", period = 1, treatment = "experimental", duration = 3, rate = .06),
+    data.frame(stratum = "High", period = 2, treatment = "experimental", duration = 3, rate = .03)
+  )
+  fail_rate <- do.call(rbind, fail_rate)
+  dropout_rate <- list(
+    data.frame(stratum = "Low", period = 1, treatment = "control", duration = 3, rate = .001),
+    data.frame(stratum = "Low", period = 1, treatment = "experimental", duration = 3, rate = .001),
+    data.frame(stratum = "High", period = 1, treatment = "control", duration = 3, rate = .001),
+    data.frame(stratum = "High", period = 1, treatment = "experimental", duration = 3, rate = .001)
+  )
+  dropout_rate <- do.call(rbind, dropout_rate)
+  set.seed(12345)
+  ex4 <- sim_pw_surv(
+    n = 12,
+    stratum = data.frame(stratum = c("Low", "High"), p = c(.3, .7)),
+    fail_rate = fail_rate,
+    dropout_rate = dropout_rate
+  )
+  saveRDS(ex4, file.path(outdir, "sim_pw_surv_ex4.rds"))
+}
diff --git a/tests/testthat/test-backwards-compatibility.R b/tests/testthat/test-backwards-compatibility.R
new file mode 100644
index 00000000..4446dbee
--- /dev/null
+++ b/tests/testthat/test-backwards-compatibility.R
@@ -0,0 +1,281 @@
+skip_if_offline()
+skip_on_cran()
+
+# see generate-backwards-compatible-test-data.R for how test data is generated
+# with a previous version of simtrial
+reference <- Sys.getenv("SIMTRIAL_TEST_BACKWARDS_COMPATIBILITY_REF")
+if (nchar(reference) > 0) {
+  source("generate-backwards-compatible-test-data.R")
+  generate_test_data(reference = reference, outdir = "fixtures/backwards-compatibility")
+  library("simtrial")
+} else {
+  skip(message = "Not testing backwards compatibility")
+}
+
+test_that("cut_data_by_date()", {
+  set.seed(12345)
+  observed <- cut_data_by_date(
+    x = sim_pw_surv(n = 20),
+    cut_date = 5
+  )
+  expected <- readRDS("fixtures/backwards-compatibility/cut_data_by_date_ex1.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+})
+
+test_that("get_cut_date_by_event()", {
+  set.seed(12345)
+  x <- sim_pw_surv(
+    n = 200,
+    stratum = data.frame(
+      stratum = c("Positive", "Negative"),
+      p = c(.5, .5)
+    ),
+    fail_rate = data.frame(
+      stratum = rep(c("Positive", "Negative"), 2),
+      period = rep(1, 4),
+      treatment = c(rep("control", 2), rep("experimental", 2)),
+      duration = rep(1, 4),
+      rate = log(2) / c(6, 9, 9, 12)
+    ),
+    dropout_rate = data.frame(
+      stratum = rep(c("Positive", "Negative"), 2),
+      period = rep(1, 4),
+      treatment = c(rep("control", 2), rep("experimental", 2)),
+      duration = rep(1, 4),
+      rate = rep(.001, 4)
+    )
+  )
+  observed <- get_cut_date_by_event(subset(x, stratum == "Positive"), event = 50)
+  expected <- readRDS("fixtures/backwards-compatibility/get_cut_date_by_event_ex1.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+})
+
+test_that("counting_process()", {
+  # Example 1
+  x <- data.frame(
+    stratum = c(rep(1, 10), rep(2, 6)),
+    treatment = rep(c(1, 1, 0, 0), 4),
+    tte = 1:16,
+    event = rep(c(0, 1), 8)
+  )
+  observed <- counting_process(x, arm = 1)
+  expected <- readRDS("fixtures/backwards-compatibility/counting_process_ex1.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+
+  # Example 2
+  set.seed(12345)
+  x <- sim_pw_surv(n = 400)
+  y <- cut_data_by_event(x, 150)
+  observed <- counting_process(y, arm = "experimental")
+  expected <- readRDS("fixtures/backwards-compatibility/counting_process_ex2.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+
+  # Example 3
+  # Counting Process Format with ties
+  x <- data.frame(
+    stratum = c(rep(1, 10), rep(2, 6)),
+    treatment = rep(c(1, 1, 0, 0), 4),
+    tte = c(rep(1:4, each = 4)),
+    event = rep(c(0, 1), 8)
+  )
+  arm <- 1
+  observed <- counting_process(x, arm)
+  expected <- readRDS("fixtures/backwards-compatibility/counting_process_ex3.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+})
+
+test_that("wlr()", {
+  # Example 1
+  # Use default enrollment and event rates at cut at 100 events
+  set.seed(12345)
+  x <- sim_pw_surv(n = 200)
+  x <- cut_data_by_event(x, 100)
+  x <- counting_process(x, arm = "experimental")
+
+  # Compute the corvariance between FH(0, 0), FH(0, 1) and FH(1, 0)
+  observed <- wlr(x, rho_gamma = data.frame(rho = c(0, 0, 1), gamma = c(0, 1, 0)))
+  expected <- readRDS("fixtures/backwards-compatibility/wlr_ex1.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+  observed <- wlr(x, rho_gamma = data.frame(rho = c(0, 0, 1), gamma = c(0, 1, 0)), return_variance = TRUE)
+  expected <- readRDS("fixtures/backwards-compatibility/wlr_ex1_var.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+  observed <- wlr(x, rho_gamma = data.frame(rho = c(0, 0, 1), gamma = c(0, 1, 0)), return_corr = TRUE)
+  expected <- readRDS("fixtures/backwards-compatibility/wlr_ex1_cor.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+
+  # Example 2
+  # Use default enrollment and event rates at cut of 100 events
+  set.seed(12345)
+  x <- sim_pw_surv(n = 200)
+  x <- cut_data_by_event(x, 100)
+  x <- counting_process(x, arm = "experimental")
+  observed <- wlr(x, rho_gamma = data.frame(rho = c(0, 0), gamma = c(0, 1)), return_corr = TRUE)
+  expected <- readRDS("fixtures/backwards-compatibility/wlr_ex2.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+})
+
+test_that("rpw_enroll()", {
+  set.seed(12345)
+  observed <- rpw_enroll(
+    n = 1e5,
+    enroll_rate = data.frame(
+      rate = c(5, 15, 30),
+      duration = c(100, 200, 100)
+    )
+  )
+  expected <- readRDS("fixtures/backwards-compatibility/rpw_enroll_ex1.rds")
+  expect_equal(observed, expected)
+
+  # Example 2
+  # Exponential enrollment
+  set.seed(12345)
+  observed <- rpw_enroll(
+    n = 1e5,
+    enroll_rate = data.frame(rate = .03, duration = 1)
+  )
+  expected <- readRDS("fixtures/backwards-compatibility/rpw_enroll_ex2.rds")
+  expect_equal(observed, expected)
+})
+
+test_that("simfix2simpwsurv()", {
+  # Example 1
+  # Convert standard input
+  observed <- simfix2simpwsurv()
+  expected <- readRDS("fixtures/backwards-compatibility/simfix2simpwsurv_ex1.rds")
+  expect_equivalent(
+    as.data.frame(observed$fail_rate),
+    as.data.frame(expected$fail_rate)
+  )
+  expect_equivalent(
+    as.data.frame(observed$dropout_rate),
+    as.data.frame(expected$dropout_rate)
+  )
+
+  # Example 2
+  # Stratified example
+  fail_rate <- data.frame(
+    stratum = c(rep("Low", 3), rep("High", 3)),
+    duration = rep(c(4, 10, 100), 2),
+    fail_rate = c(
+      .04, .1, .06,
+      .08, .16, .12
+    ),
+    hr = c(
+      1.5, .5, 2 / 3,
+      2, 10 / 16, 10 / 12
+    ),
+    dropout_rate = .01
+  )
+  observed <- simfix2simpwsurv(fail_rate)
+  expected <- readRDS("fixtures/backwards-compatibility/simfix2simpwsurv_ex2.rds")
+  expect_equivalent(
+    as.data.frame(observed$fail_rate),
+    as.data.frame(expected$fail_rate)
+  )
+  expect_equivalent(
+    as.data.frame(observed$dropout_rate),
+    as.data.frame(expected$dropout_rate)
+  )
+})
+
+test_that("sim_pw_surv()", {
+  # Example 1
+  set.seed(12345)
+  observed <- sim_pw_surv(n = 20)
+  expected <- readRDS("fixtures/backwards-compatibility/sim_pw_surv_ex1.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+
+  # Example 2
+  # 3:1 randomization
+  set.seed(12345)
+  observed <- sim_pw_surv(
+    n = 20,
+    block = c(rep("experimental", 3), "control")
+  )
+  expected <- readRDS("fixtures/backwards-compatibility/sim_pw_surv_ex2.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+
+  # Example 3
+  # Simulate 2 stratum; will use defaults for blocking and enrollRates
+  set.seed(12345)
+  observed <- sim_pw_surv(
+    n = 20,
+    # 2 stratum,30% and 70% prevalence
+    stratum = data.frame(stratum = c("Low", "High"), p = c(.3, .7)),
+    fail_rate = data.frame(
+      stratum = c(rep("Low", 4), rep("High", 4)),
+      period = rep(1:2, 4),
+      treatment = rep(c(
+        rep("control", 2),
+        rep("experimental", 2)
+      ), 2),
+      duration = rep(c(3, 1), 4),
+      rate = c(.03, .05, .03, .03, .05, .08, .07, .04)
+    ),
+    dropout_rate = data.frame(
+      stratum = c(rep("Low", 2), rep("High", 2)),
+      period = rep(1, 4),
+      treatment = rep(c("control", "experimental"), 2),
+      duration = rep(1, 4),
+      rate = rep(.001, 4)
+    )
+  )
+  expected <- readRDS("fixtures/backwards-compatibility/sim_pw_surv_ex3.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+
+  # Example 4
+  # If you want a more rectangular entry for a tibble
+  fail_rate <- list(
+    data.frame(stratum = "Low", period = 1, treatment = "control", duration = 3, rate = .03),
+    data.frame(stratum = "Low", period = 1, treatment = "experimental", duration = 3, rate = .03),
+    data.frame(stratum = "Low", period = 2, treatment = "experimental", duration = 3, rate = .02),
+    data.frame(stratum = "High", period = 1, treatment = "control", duration = 3, rate = .05),
+    data.frame(stratum = "High", period = 1, treatment = "experimental", duration = 3, rate = .06),
+    data.frame(stratum = "High", period = 2, treatment = "experimental", duration = 3, rate = .03)
+  )
+  fail_rate <- do.call(rbind, fail_rate)
+  dropout_rate <- list(
+    data.frame(stratum = "Low", period = 1, treatment = "control", duration = 3, rate = .001),
+    data.frame(stratum = "Low", period = 1, treatment = "experimental", duration = 3, rate = .001),
+    data.frame(stratum = "High", period = 1, treatment = "control", duration = 3, rate = .001),
+    data.frame(stratum = "High", period = 1, treatment = "experimental", duration = 3, rate = .001)
+  )
+  dropout_rate <- do.call(rbind, dropout_rate)
+  set.seed(12345)
+  observed <- sim_pw_surv(
+    n = 12,
+    stratum = data.frame(stratum = c("Low", "High"), p = c(.3, .7)),
+    fail_rate = fail_rate,
+    dropout_rate = dropout_rate
+  )
+  expected <- readRDS("fixtures/backwards-compatibility/sim_pw_surv_ex4.rds")
+  expect_equivalent(as.data.frame(observed), as.data.frame(expected))
+})
+
+test_that("sim_fixed_n()", {
+  # Example 1
+  # Show output structure
+  set.seed(12345)
+  observed <- sim_fixed_n(n = 2)
+  expected <- readRDS("fixtures/backwards-compatibility/sim_fixed_n_ex1.rds")
+  expect_equal(observed, expected)
+
+  # Example 2
+  # Example with 2 tests: logrank and FH(0,1)
+  set.seed(12345)
+  observed <- sim_fixed_n(n = 2, rho_gamma = data.frame(rho = 0, gamma = c(0, 1)))
+  expected <- readRDS("fixtures/backwards-compatibility/sim_fixed_n_ex2.rds")
+  expect_equal(observed, expected)
+
+  # Example 3
+  # Power by test
+  # Only use cuts for events, events + min follow-up
+  set.seed(12345)
+  observed <- sim_fixed_n(
+    n_sim = 2,
+    timing_type = c(2, 5),
+    rho_gamma = data.frame(rho = 0, gamma = c(0, 1))
+  )
+  expected <- readRDS("fixtures/backwards-compatibility/sim_fixed_n_ex3.rds")
+  expect_equal(observed, expected)
+})
diff --git a/tests/testthat/test-independent_test_counting_process.R b/tests/testthat/test-independent_test_counting_process.R
index 3aa042b5..811e5461 100644
--- a/tests/testthat/test-independent_test_counting_process.R
+++ b/tests/testthat/test-independent_test_counting_process.R
@@ -42,7 +42,7 @@ surv_to_count <- function(time, status, trt, strats) {
 
   # Log Rank Expectation Difference and Variance
   res <- merge(km, km_by_trt, all = TRUE) %>%
-    arrange(trt, strats, time) %>%
+    dplyr::arrange(trt, strats, time) %>%
     mutate(
       OminusE = tn.event - tn.risk / n.risk * n.event,
       Var = (n.risk - tn.risk) * tn.risk * n.event * (n.risk - n.event) / n.risk^2 / (n.risk - 1)
@@ -61,7 +61,7 @@ testthat::test_that("Counting Process Format without ties", {
   res_counting_process <- simtrial::counting_process(x, arm)
   res_test <- surv_to_count(time = x$tte, status = x$event, trt = x$treatment, strats = x$stratum)
 
-  res_test <- as_tibble(subset(res_test, trt == 1)) %>%
+  res_test <- tibble::as_tibble(subset(res_test, trt == 1)) %>%
     subset(n.event > 0 & n.risk - tn.risk > 0 & tn.risk > 0)
 
   testthat::expect_equal(res_counting_process$o_minus_e, res_test$OminusE)
@@ -80,7 +80,7 @@ testthat::test_that("Counting Process Format with ties", {
   res_counting_process <- counting_process(x, arm)
   res_test <- surv_to_count(time = x$tte, status = x$event, trt = x$treatment, strats = x$stratum)
 
-  res_test <- as_tibble(subset(res_test, trt == 1)) %>%
+  res_test <- tibble::as_tibble(subset(res_test, trt == 1)) %>%
     subset(n.event > 0 & n.risk - tn.risk > 0 & tn.risk > 0)
 
   testthat::expect_equal(res_counting_process$o_minus_e, res_test$OminusE)
diff --git a/tests/testthat/test-independent_test_getCutDateForCount.R b/tests/testthat/test-independent_test_getCutDateForCount.R
index 4effa1e8..17e01ca6 100644
--- a/tests/testthat/test-independent_test_getCutDateForCount.R
+++ b/tests/testthat/test-independent_test_getCutDateForCount.R
@@ -4,9 +4,9 @@ test_that("get_cut_date_by_event returns the correct cut date", {
   ycutdate <- get_cut_date_by_event(y, event)
 
   x <- y %>%
-    ungroup() %>%
+    dplyr::ungroup() %>%
     filter(fail == 1) %>%
-    arrange(cte) %>%
+    dplyr::arrange(cte) %>%
     dplyr::slice(event)
   expect_equal(x$cte, ycutdate)
 })