diff --git a/.Rbuildignore b/.Rbuildignore
index 9950f979..71fb8e5d 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -23,3 +23,4 @@
 ^pammtools\.Rproj$
 ^revdep$
 ^\.github$
+^CRAN-SUBMISSION$
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index cea7ca2d..a3ac6182 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -1,14 +1,10 @@
-# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag.
-# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches:
-      - main
-      - master
+    branches: [main, master]
   pull_request:
-    branches:
-      - main
-      - master
+    branches: [main, master]
 
 name: R-CMD-check
 
@@ -22,64 +18,32 @@ jobs:
       fail-fast: false
       matrix:
         config:
+          - {os: macos-latest,   r: 'release'}
           - {os: windows-latest, r: 'release'}
-          - {os: macOS-latest, r: 'release'}
-          - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
-          - {os: ubuntu-20.04,   r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest", http-user-agent: "R/4.1.0 (ubuntu-20.04) R (4.1.0 x86_64-pc-linux-gnu x86_64 linux-gnu) on GitHub Actions" }
+          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
+          - {os: ubuntu-latest,   r: 'release'}
+          - {os: ubuntu-latest,   r: 'oldrel-1'}
 
     env:
-      R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
-      RSPM: ${{ matrix.config.rspm }}
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
 
-      - uses: r-lib/actions/setup-pandoc@v1
-
-      - name: Query dependencies
-        run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
-        shell: Rscript {0}
-
-      - name: Restore R package cache
-        uses: actions/cache@v2
+      - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
-
-      - name: Install system dependencies
-        if: runner.os == 'Linux'
-        run: |
-          while read -r cmd
-          do
-            eval sudo $cmd
-          done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))')
-
-      - name: Install dependencies
-        run: |
-          remotes::install_deps(dependencies = TRUE)
-          remotes::install_cran("rcmdcheck")
-        shell: Rscript {0}
-
-      - name: Check
-        env:
-          _R_CHECK_CRAN_INCOMING_REMOTE_: false
-        run: |
-          options(crayon.enabled = TRUE)
-          rcmdcheck::rcmdcheck(args = c("--no-manual", "--as-cran"), error_on = "warning", check_dir = "check")
-        shell: Rscript {0}
+          extra-packages: any::rcmdcheck
+          needs: check
 
-      - name: Upload check results
-        if: failure()
-        uses: actions/upload-artifact@main
+      - uses: r-lib/actions/check-r-package@v2
         with:
-          name: ${{ runner.os }}-r${{ matrix.config.r }}-results
-          path: check
+          upload-snapshots: true
diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
index b79ffb49..ed7650c7 100644
--- a/.github/workflows/pkgdown.yaml
+++ b/.github/workflows/pkgdown.yaml
@@ -1,55 +1,48 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches:
-      - main
-      - master
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+  release:
+    types: [published]
+  workflow_dispatch:
 
 name: pkgdown
 
 jobs:
   pkgdown:
-    runs-on: macOS-latest
+    runs-on: ubuntu-latest
+    # Only restrict concurrency for non-PR jobs
+    concurrency:
+      group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }}
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    permissions:
+      contents: write
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-pandoc@v2
 
-      - uses: r-lib/actions/setup-pandoc@v1
-
-      - name: Query dependencies
-        run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
-        shell: Rscript {0}
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
 
-      - name: Restore R package cache
-        uses: actions/cache@v2
+      - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
+          extra-packages: any::pkgdown, local::.
+          needs: website
 
-      - name: Install dependencies
-        run: |
-          remotes::install_deps(dependencies = TRUE)
-          remotes::install_cran("pkgdown")
-          # Vignette dependencies
-          remotes::install_cran("TBFmultinomial")
-          remotes::install_cran("etm")
-          remotes::install_cran("scam")
-          remotes::install_cran("coxme")
-          remotes::install_cran("eha")
-          remotes::install_cran("timereg")
+      - name: Build site
+        run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE)
         shell: Rscript {0}
 
-      - name: Install package
-        run: R CMD INSTALL .
-
-      - name: Deploy package
-        run: |
-          git config --local user.email "actions@github.com"
-          git config --local user.name "GitHub Actions"
-          Rscript -e 'pkgdown::deploy_to_branch(new_process = FALSE)'
+      - name: Deploy to GitHub pages 🚀
+        if: github.event_name != 'pull_request'
+        uses: JamesIves/github-pages-deploy-action@v4.4.1
+        with:
+          clean: false
+          branch: gh-pages
+          folder: docs
diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index a1205cce..2c5bb502 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -1,48 +1,50 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches:
-      - main
-      - master
+    branches: [main, master]
   pull_request:
-    branches:
-      - main
-      - master
+    branches: [main, master]
 
 name: test-coverage
 
 jobs:
   test-coverage:
-    runs-on: macOS-latest
+    runs-on: ubuntu-latest
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
 
-      - uses: r-lib/actions/setup-pandoc@v1
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::covr
+          needs: coverage
 
-      - name: Query dependencies
+      - name: Test coverage
         run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
+          covr::codecov(
+            quiet = FALSE,
+            clean = FALSE,
+            install_path = file.path(Sys.getenv("RUNNER_TEMP"), "package")
+          )
         shell: Rscript {0}
 
-      - name: Restore R package cache
-        uses: actions/cache@v2
-        with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
-
-      - name: Install dependencies
+      - name: Show testthat output
+        if: always()
         run: |
-          install.packages(c("remotes"))
-          remotes::install_deps(dependencies = TRUE)
-          remotes::install_cran("covr")
-        shell: Rscript {0}
+          ## --------------------------------------------------------------------
+          find ${{ runner.temp }}/package -name 'testthat.Rout*' -exec cat '{}' \; || true
+        shell: bash
 
-      - name: Test coverage
-        run: covr::codecov(type = "all", line_exclusions = "R/zzz.R")
-        shell: Rscript {0}
+      - name: Upload test results
+        if: failure()
+        uses: actions/upload-artifact@v3
+        with:
+          name: coverage-test-failures
+          path: ${{ runner.temp }}/package
diff --git a/DESCRIPTION b/DESCRIPTION
index 9542ccc2..d4a9f62f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -35,6 +35,12 @@ Imports:
     vctrs (>= 0.3.0)
 Suggests:
     testthat
+Config/Needs/website:
+  coxme,
+  eha,
+  etm,
+  scam,
+  TBFmultinomial
 License: MIT + file LICENSE
 LazyData: true
 URL: https://adibender.github.io/pammtools/
diff --git a/NAMESPACE b/NAMESPACE
index 0ebe7c0c..82f0e784 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -61,7 +61,6 @@ export(GeomStepHazard)
 export(GeomStepribbon)
 export(GeomSurv)
 export(add_cif)
-export(add_counterfactual_transitions)
 export(add_cumu_hazard)
 export(add_hazard)
 export(add_surv_prob)
@@ -117,7 +116,6 @@ export(sample_n)
 export(select)
 export(seq_range)
 export(sim_pexp)
-export(sim_pexp_msm)
 export(slice)
 export(split_data)
 export(summarise)
@@ -160,7 +158,6 @@ importFrom(mgcv,predict.gam)
 importFrom(mvtnorm,rmvnorm)
 importFrom(pec,predictSurvProb)
 importFrom(purrr,compose)
-importFrom(purrr,cross)
 importFrom(purrr,cross_df)
 importFrom(purrr,discard)
 importFrom(purrr,flatten)
@@ -173,7 +170,6 @@ importFrom(purrr,map_int)
 importFrom(purrr,map_lgl)
 importFrom(purrr,reduce)
 importFrom(purrr,set_names)
-importFrom(purrr,transpose)
 importFrom(rlang,":=")
 importFrom(rlang,.data)
 importFrom(rlang,.env)
@@ -200,6 +196,7 @@ importFrom(stats,vcov)
 importFrom(tibble,as_tibble)
 importFrom(tidyr,complete)
 importFrom(tidyr,crossing)
+importFrom(tidyr,expand_grid)
 importFrom(tidyr,gather)
 importFrom(tidyr,nest)
 importFrom(tidyr,pivot_longer)
diff --git a/NEWS.md b/NEWS.md
index 9a3a7ef2..66e8f07c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,12 @@
+# pamtools 0.5.93
++ Maintnance (some tidyverse deprecations, link fixes, etc., smaller bugs)
+
+# pammtools 0.5.92
++ Fixed competing risks data trafo in case of more than 2 causes
+
+# pammtools 0.5.9
++ Fixes issue 154: direction argument to `geom_stepribbon`
+
 # pammtools 0.5.8
 + removed argument `methods` from `pamm`. Can be specified via `...`. Fixes #200
 + adapted `warn_about_new_time_points` when original data not stored in model object. Fixes #203
diff --git a/R/add-functions.R b/R/add-functions.R
index 368bb624..e5aa5b1e 100644
--- a/R/add-functions.R
+++ b/R/add-functions.R
@@ -63,36 +63,43 @@ add_term <- function(
     se <- unname(sqrt(rowSums( (X %*% cov.coefs) * X )))
     newdata <- newdata %>%
       mutate(
-        ci_lower = .data$fit - se_mult * se,
-        ci_upper = .data$fit + se_mult * se)
+        ci_lower = .data[["fit"]] - se_mult * se,
+        ci_upper = .data[["fit"]] + se_mult * se)
   }
 
   return(newdata)
 
 }
 
+#' Create model matrix from model object
+#' @param object a model object
 make_X <- function(object, ...) {
 
   UseMethod("make_X", object)
 
 }
 
+#' @inherit make_X
 make_X.default <- function(object, newdata, ...) {
 
   X <- model.matrix(object$formula[-2], data = newdata, ...)
 
 }
 
+#' @inherit make_X
+#' @keywords internal
 make_X.gam <- function(object, newdata, ...) {
 
   X <- predict.gam(object, newdata = newdata, type = "lpmatrix", ...)
 
 }
 
+#' @inherit make_X
+#' @keywords internal
 make_X.scam <- function(object, newdata, ...) {
-  
+
   X <- predict.scam(object, newdata = newdata, type = "lpmatrix", ...)
-  
+
 }
 
 prep_X <- function(object, newdata, reference = NULL, ...) {
@@ -258,7 +265,7 @@ get_hazard.default <- function(
       add_ci(object, X, type = type, ci_type = ci_type, se_mult = se_mult, ...)
   }
   if (type == "response") {
-    newdata <- newdata %>% mutate(hazard = exp(.data$hazard))
+    newdata <- newdata %>% mutate(hazard = exp(.data[["hazard"]]))
   }
 
   newdata %>% arrange(.data[[time_var]], .by_group = TRUE)
@@ -329,7 +336,7 @@ get_cumu_hazard <- function(
 
   interval_length <- sym(interval_length)
 
-  mutate_args  <- list(cumu_hazard = quo(cumsum(.data$hazard *
+  mutate_args  <- list(cumu_hazard = quo(cumsum(.data[["hazard"]] *
     (!!interval_length))))
   haz_vars_in_data <- map(c("hazard", "se", "ci_lower", "ci_upper"),
     ~ grep(.x, colnames(newdata), value = TRUE, fixed = TRUE)) %>% flatten_chr()
@@ -343,8 +350,8 @@ get_cumu_hazard <- function(
       if (ci_type == "default") {
         mutate_args <- mutate_args %>%
           append(list(
-            cumu_lower = quo(cumsum(.data$ci_lower * (!!interval_length))),
-            cumu_upper = quo(cumsum(.data$ci_upper * (!!interval_length)))))
+            cumu_lower = quo(cumsum(.data[["ci_lower"]] * (!!interval_length))),
+            cumu_upper = quo(cumsum(.data[["ci_upper"]] * (!!interval_length)))))
       } else {
         # ci delta rule
         newdata <- split(newdata, group_indices(newdata)) %>%
@@ -443,7 +450,7 @@ get_surv_prob <- function(
 
   interval_length <- sym(interval_length)
 
-  mutate_args  <- list(surv_prob = quo(exp(-cumsum(.data$hazard *
+  mutate_args  <- list(surv_prob = quo(exp(-cumsum(.data[["hazard"]] *
     (!!interval_length)))))
   haz_vars_in_data <- map(c("hazard", "se", "ci_lower", "ci_upper"),
     ~grep(.x, colnames(newdata), value = TRUE, fixed = TRUE)) %>% flatten_chr()
@@ -457,8 +464,8 @@ get_surv_prob <- function(
       if (ci_type == "default") {
         mutate_args <- mutate_args %>%
           append(list(
-            surv_upper = quo(exp(-cumsum(.data$ci_lower * (!!interval_length)))),
-            surv_lower = quo(exp(-cumsum(.data$ci_upper * (!!interval_length))))))
+            surv_upper = quo(exp(-cumsum(.data[["ci_lower"]] * (!!interval_length)))),
+            surv_lower = quo(exp(-cumsum(.data[["ci_upper"]] * (!!interval_length))))))
       } else {
         # ci delta rule
         newdata <- split(newdata, group_indices(newdata)) %>%
@@ -512,16 +519,16 @@ add_ci <- function(
   if (type == "link") {
     newdata <- newdata %>%
       mutate(
-        ci_lower = .data$hazard - se_mult * .data$se,
-        ci_upper = .data$hazard + se_mult * .data$se)
+        ci_lower = .data[["hazard"]] - se_mult * .data[["se"]],
+        ci_upper = .data[["hazard"]] + se_mult * .data[["se"]])
   }
 
   if (type != "link") {
     if (ci_type == "default") {
       newdata <- newdata %>%
         mutate(
-          ci_lower = exp(.data$hazard - se_mult * .data$se),
-          ci_upper = exp(.data$hazard + se_mult * .data$se))
+          ci_lower = exp(.data[["hazard"]] - se_mult * .data[["se"]]),
+          ci_upper = exp(.data[["hazard"]] + se_mult * .data[["se"]]))
     } else {
       if (ci_type == "delta") {
         newdata <- split(newdata, group_indices(newdata)) %>%
@@ -545,8 +552,8 @@ add_delta_ci <- function(newdata, object, se_mult = 2, ...) {
   newdata %>%
     mutate(
       se       = sqrt(rowSums( (Jacobi %*% V) * Jacobi )),
-      ci_lower = exp(.data$hazard) - .data$se * se_mult,
-      ci_upper = exp(.data$hazard) + .data$se * se_mult)
+      ci_lower = exp(.data[["hazard"]]) - .data[["se"]] * se_mult,
+      ci_upper = exp(.data[["hazard"]]) + .data[["se"]] * se_mult)
 
 }
 
@@ -560,8 +567,8 @@ add_delta_ci_cumu <- function(newdata, object, se_mult = 2, ...) {
   newdata %>%
     mutate(
       se       = sqrt(rowSums( (LHS %*% V) * LHS )),
-      cumu_lower = cumsum(.data$intlen * .data$hazard) - .data$se * se_mult,
-      cumu_upper = cumsum(.data$intlen * .data$hazard) + .data$se * se_mult)
+      cumu_lower = cumsum(.data[["intlen"]] * .data[["hazard"]]) - .data[["se"]] * se_mult,
+      cumu_upper = cumsum(.data[["intlen"]] * .data[["hazard"]]) + .data[["se"]] * se_mult)
 
 }
 
@@ -576,8 +583,8 @@ add_delta_ci_surv <- function(newdata, object, se_mult = 2, ...) {
   newdata %>%
     mutate(
       se       = sqrt(rowSums( (LHS %*% V) * LHS)),
-      surv_lower = exp(-cumsum(.data$hazard * .data$intlen)) - .data$se * se_mult,
-      surv_upper = exp(-cumsum(.data$hazard * .data$intlen)) + .data$se * se_mult)
+      surv_lower = exp(-cumsum(.data[["hazard"]] * .data[["intlen"]])) - .data[["se"]] * se_mult,
+      surv_upper = exp(-cumsum(.data[["hazard"]] * .data[["intlen"]])) + .data[["se"]] * se_mult)
 
 }
 
@@ -768,7 +775,7 @@ get_trans_prob <- function(
     tend = "tend",
     cumu_hazard = "cumu_hazard",
     ...) {
-  
+
   # interval_length
   assert_character(interval_length)
   assert_subset(interval_length, colnames(newdata))
@@ -782,49 +789,49 @@ get_trans_prob <- function(
   assert_character(cumu_hazard)
   assert_subset(cumu_hazard, colnames(newdata))
   assert_data_frame(newdata, all.missing = FALSE)
-  
+
   # assert_class(object, classes = "glm")
-  
+
   interval_length <- sym(interval_length)
   transition <- sym(transition)
-  
+
   # include from and to, to obtain transition probability in multidim array
-  newdata <- newdata %>% 
+  newdata <- newdata %>%
     mutate(from = as.numeric(gsub("->.*", "", !!transition))
            , to = as.numeric(gsub(".*->", "", !!transition)))
-  
+
   # get unique transitions to build transition matrix
   unique_transition <- data.frame(unique(newdata %>% select(!!transition, from, to)))
   # get unique time points
-  unique_tend <- data.frame(unique(newdata %>% 
-                                     ungroup(!!transition) %>% 
+  unique_tend <- data.frame(unique(newdata %>%
+                                     ungroup(!!transition) %>%
                                      select(!!tend)))
-  
+
   # transition matrix
   m <- sapply(unique_transition[,c(2,3)], max) + 1 #transition starts at 0, integer of matrix at 1
-  M <- array(0, dim=c(max(m), max(m), nrow(unique_transition))) 
-  
-  
+  M <- array(0, dim=c(max(m), max(m), nrow(unique_transition)))
+
+
   # create transition matrices to be used at every time point,
   # multiply matrices with "scalar" alpha_ij_k which is the delta cumu hazard at time t_k for transition i->j
-  
+
   for (iter in 1:nrow(unique_transition)){
     M[unique_transition$from[iter] + 1, unique_transition$to[iter] + 1,iter] <- 1
     M[unique_transition$from[iter] + 1, unique_transition$from[iter] + 1,iter] <- -1
   }
-  
+
   # add cumu hazards to dataset
-  newdata <- newdata %>% 
+  newdata <- newdata %>%
     # group_by(!!transition) %>%
     mutate(delta_cumu_hazard = cumu_hazard - ifelse(is.na(lag(cumu_hazard)), 0, lag(cumu_hazard)))
-  
+
   # create dA array, to calculate transition probabilities
   alpha <- array(rep(0, nrow(unique_tend)*nrow(unique_transition)), dim=c(nrow(unique_tend), nrow(unique_transition)))
   I <- array(rep(diag(max(m)), nrow(unique_tend))
              , dim=c( max(m), max(m), nrow(unique_tend)))
   A <- array(0, dim=c(max(m), max(m), nrow(unique_tend)))
   cum_A <- array(0, dim=c(max(m), max(m), nrow(unique_tend)))
-  
+
   # calculate differences in hazards
   alpha <- sapply(1:nrow(unique_transition), function(iter){
     val <- newdata %>% ungroup() %>% filter(transition == unique_transition[iter,1]) %>% arrange(tend)
@@ -835,10 +842,10 @@ get_trans_prob <- function(
       A[,,t] <- A[,,t] + M[,,trans] * alpha[t, trans]
     }
   }
-  
+
   # prepare transition probabilities
   A <- I + A
-  
+
   for (iter in 1:nrow(unique_tend)) {
     if (iter == 1) {
       cum_A[,,iter] = A[,,iter]
@@ -846,21 +853,21 @@ get_trans_prob <- function(
       cum_A[,,iter] = round(cum_A[,,iter-1] %*% A[,,iter],10) #use matrix multiplikation
     }
   }
-  
+
   # transform array so that transition probability can be joined via tend and transition
   tmp <- cbind(unique_tend
                , sapply(1:nrow(unique_transition), function(row) cum_A[unique_transition$from[row] + 1, unique_transition$to[row] + 1, ]))
   colnames(tmp) <- c("tend", as.character(unique_transition$transition))
   trans_prob_df <- tmp %>%
     pivot_longer(cols = c(as.character(unique_transition$transition)), names_to = "transition", values_to = "trans_prob")
-  
+
   # join probabilities and return matrix
   newdata <- newdata %>%
     left_join(trans_prob_df, by=c("tend", "transition")) %>%
     select(-delta_cumu_hazard, -from, -to)
-  
+
   return(newdata)
-  
+
 }
 
 #' Add transition probabilities
@@ -868,15 +875,15 @@ get_trans_prob <- function(
 add_trans_prob <- function(
     newdata
     , object
-    , overwrite       = FALSE 
+    , overwrite       = FALSE
     , alpha           = 0.05
     , n_sim           = 100L
     , time_var        = NULL
     , interval_length = "intlen",
     ...
 ) {
-  
-  
+
+
   if (!overwrite) {
     if ("trans_prob" %in% names(newdata)) {
       stop("Data set already contains 'trans_prob' column.
@@ -891,69 +898,69 @@ add_trans_prob <- function(
       names(newdata))
     newdata <- newdata %>% select(-one_of(rm.vars))
   }
-  
+
   # extract to simulate ci
   coefs         <- coef(object)
   V             <- object$Vp
   sim_coef_mat  <- mvtnorm::rmvnorm(n_sim, mean = coefs, sigma = V)
-  
-  
+
+
   newdata <- newdata %>% add_cumu_hazard(object)
-  
+
   old_groups <- dplyr::groups(newdata)
   res_data <- newdata %>% ungroup(transition)
-  out_data <- group_split(res_data) |> 
-    map(res_data, .f = ~ group_by(.x, transition))|> 
+  out_data <- group_split(res_data) |>
+    map(res_data, .f = ~ group_by(.x, transition))|>
     map(res_data, .f = ~ get_trans_prob(.x)) |>
     map(res_data, .f = ~ group_by(.x, !!!old_groups)) |>
     bind_rows()
-  
+
   return(out_data)
-  
+
 }
 
 #' helper function for add_trans_ci
 #' @keywords internal
 get_sim_cumu <- function(newdata, ...) {
-  
+
   newdata$cumu_hazard <- cumsum(newdata$intlen * newdata$hazard)
-  
+
   newdata
-  
+
 }
 
 #' Add transition probabilities confidence intervals
 #' @export
 add_trans_ci <- function(newdata, object, n_sim=100L, alpha=0.05, ...) {
-  
+
   X             <- predict.gam(object, newdata = newdata, type = "lpmatrix")
   coefs         <- coef(object)
   V             <- object$Vp
-  
+
   # define groups: 1. all grouping variables -> cumu hazards, 2. all but transition -> trans_prob
   groups_array <- group_indices(newdata)
   groups_trans <- newdata %>% ungroup(transition) %>% group_indices()
-  
+
   sim_coef_mat <- mvtnorm::rmvnorm(n_sim, mean = coefs, sigma = V)
   sim_fit_mat <- apply(sim_coef_mat, 1, function(z)
     exp(X %*% z))
-  
+
   # create list with replicated newdata
   nlst <- as.list(replicate(n_sim, newdata[,c("tend", "transition", "intlen")], simplify=F))
-  
+
   # add cumu-hazard in each element and calculate trans_prob with perturbed hazards
   nlst <- lapply(1:n_sim, function(i) {
     nlst[[i]] <- cbind(nlst[[i]], hazard = sim_fit_mat[, i]) # add hazard
     # split by group and calculate cumu hazard
     nlst[[i]] <- split(nlst[[i]], groups_array) %>% map_dfr(get_sim_cumu)
     nlst[[i]] <- split(nlst[[i]], groups_trans) %>% map_dfr(get_trans_prob)
-    
+
     nlst[[i]]
   })
-  
+
   sim_trans_probs <- do.call(cbind, lapply(nlst, function(df) df$trans_prob))
   newdata$trans_lower <- apply(sim_trans_probs, 1, quantile, probs = alpha / 2)
   newdata$trans_upper <- apply(sim_trans_probs, 1, quantile, probs = 1 - alpha / 2)
-  
+
   newdata
 }
diff --git a/R/as-ped.R b/R/as-ped.R
index fac4bb1d..dc9709ba 100644
--- a/R/as-ped.R
+++ b/R/as-ped.R
@@ -203,7 +203,7 @@ is.ped <- function(x) inherits(x, "ped")
 
 #' @rdname as_ped
 #' @param newdata A new data set (\code{data.frame}) that contains the same
-#' variables that were used to create the PED object (code{data}).
+#' variables that were used to create the PED object (\code{data}).
 #' @export
 as_ped.ped <- function(data, newdata, ...) {
 
@@ -251,7 +251,9 @@ as_ped_cr <- function(
   ...) {
 
   lhs_vars <- get_lhs_vars(formula)
+  n_lhs <- length(lhs_vars)
   event_types <- get_event_types(data, formula, censor_code)
+  n_events <- sum(event_types != censor_code)
 
   cut <- map2(
     event_types,
@@ -269,7 +271,7 @@ as_ped_cr <- function(
     cut,
     function(.event, .cut) {
       ped_i <- data %>%
-        mutate(!!lhs_vars[2] := 1L * (.data[[lhs_vars[2]]] == .env[[".event"]])) %>%
+        mutate(!!lhs_vars[n_lhs] := 1L * (.data[[lhs_vars[n_lhs]]] == .env[[".event"]])) %>%
         as_ped(
           formula      = formula,
           cut          = .cut,
diff --git a/R/convenience-plots.R b/R/convenience-plots.R
index b97624dd..223996be 100644
--- a/R/convenience-plots.R
+++ b/R/convenience-plots.R
@@ -24,10 +24,10 @@ gg_smooth.default <- function(x, fit, ...) {
 
   sobj <- get_terms(data = x, fit = fit, ...)
 
-  ggsmooth <- ggplot(sobj, aes_string(x = "x", y = "eff", group = "term")) +
+  ggsmooth <- ggplot(sobj, aes(x = .data[["x"]], y = .data[["eff"]], group = .data[["term"]])) +
     geom_hline(yintercept = 0, lty = 3) +
     geom_line() +
-    geom_ribbon(aes_string(ymin = "ci_lower", ymax = "ci_upper"), alpha = 0.2) +
+    geom_ribbon(aes(ymin = .data[["ci_lower"]], ymax = .data[["ci_upper"]]), alpha = 0.2) +
     facet_wrap(~term, scales = "free_x") +
     ylab(expression(f[p](x[p]))) + xlab(expression(x[p]))
 
@@ -60,8 +60,8 @@ gg_tensor <- function(x, ci = FALSE, ...) {
         levels = c("ci_lower", "fit", "ci_upper")))
   }
 
-  gg2d <- ggplot(df2d, aes_string(x = "x", y = "y", z = "fit")) +
-    geom_raster(aes_string(fill = "fit")) +
+  gg2d <- ggplot(df2d, aes(x = .data[["x"]], y = .data[["y"]], z = .data[["fit"]])) +
+    geom_raster(aes(fill = .data[["fit"]])) +
     scale_fill_gradient2(
       name = expression(f(list(x, y))),
       low  = "steelblue", high = "firebrick2") +
@@ -96,9 +96,9 @@ gg_tensor <- function(x, ci = FALSE, ...) {
 gg_re <- function(x, ...) {
 
   re <- tidy_re(x, ...)
-  ggplot(re, aes_string(sample = "fit")) +
-    geom_abline(aes_string(intercept = "qqintercept", slope = "qqslope")) +
-    geom_qq(distribution = "qnorm") +
+  ggplot(re, aes(sample = .data[["fit"]])) +
+    geom_abline(aes(intercept = .data[["qqintercept"]], slope = .data[["qqslope"]])) +
+    geom_qq(distribution = stats::qnorm) +
     facet_wrap(~main) +
     theme_set(theme_bw())
 
@@ -122,8 +122,8 @@ gg_fixed <- function(x, intercept=FALSE, ...) {
 
   fixed_df <- tidy_fixed(x, intercept = intercept, ...)
 
-  ggplot(fixed_df, aes_string(x = "variable", y = "coef", ymin = "ci_lower",
-      ymax = "ci_upper")) +
+  ggplot(fixed_df, aes(x = .data[["variable"]], y = .data[["coef"]], ymin = .data[["ci_lower"]],
+      ymax = .data[["ci_upper"]])) +
     geom_hline(yintercept = 0, lty = 3) +
     geom_pointrange() +
     coord_flip() +
diff --git a/R/cumulative-coefficient.R b/R/cumulative-coefficient.R
index 1f737e57..64a733ab 100644
--- a/R/cumulative-coefficient.R
+++ b/R/cumulative-coefficient.R
@@ -44,11 +44,11 @@ get_cumu_coef.aalen <- function(model, data = NULL, terms, ci = TRUE, ...) {
   cumu_coef <- model[["cum"]] %>%
     as_tibble() %>%
     select(one_of(terms)) %>%
-    gather("variable", "cumu_hazard", -.data$time)
+    gather("variable", "cumu_hazard", -.data[["time"]])
   cumu_var <- model[["var.cum"]] %>%
     as_tibble() %>%
     select(terms) %>%
-    gather("variable", "cumu_var", -.data$time)
+    gather("variable", "cumu_var", -.data[["time"]])
 
     suppressMessages(
       left_join(cumu_coef, cumu_var) %>%
diff --git a/R/formula-specials.R b/R/formula-specials.R
index 3b04b71a..2d765e61 100644
--- a/R/formula-specials.R
+++ b/R/formula-specials.R
@@ -310,35 +310,39 @@ make_mat_names <- function(x, ...) {
 
 #' @keywords internal
 make_mat_names.default <- function(
-  col_vars,
+  x,
   latency_var = NULL,
   tz_var      = NULL,
   suffix      = NULL,
-  nfunc       = 1) {
+  nfunc       = 1,
+  ...) {
 
   if (!is.null(suffix)) {
-    return(paste(col_vars, suffix, sep = "_"))
+    return(paste(x, suffix, sep = "_"))
   } else {
     if (!is.null(tz_var) & nfunc > 1)  {
-      tz_ind <- col_vars == tz_var
-      col_vars[!tz_ind] <- paste(col_vars[!tz_ind], tz_var,  sep = "_")
+      tz_ind <- x == tz_var
+      x[!tz_ind] <- paste(x[!tz_ind], tz_var,  sep = "_")
     }
     if (!is.null(latency_var)) {
-      latency_ind <- col_vars == latency_var
-      col_vars[latency_ind] <- paste(col_vars[latency_ind], "latency",
+      latency_ind <- x == latency_var
+      x[latency_ind] <- paste(x[latency_ind], "latency",
         sep = "_")
     }
   }
 
-  return(col_vars)
+  return(x)
 
 }
 
 #' @keywords internal
-make_mat_names.list <- function(func_list, time_var) {
-  hist_names <- map(func_list, ~ make_mat_names(c(.x[["col_vars"]], "LL"),
+make_mat_names.list <- function(
+  x,
+  time_var,
+  ...) {
+  hist_names <- map(x, ~ make_mat_names(c(.x[["col_vars"]], "LL"),
     .x[["latency_var"]], .x[["tz_var"]], .x[["suffix"]],
-    nfunc = length(func_list)))
+    nfunc = length(x)))
 
   time_mat_ind <- map(hist_names, ~grepl(time_var, .))
   for (i in seq_along(time_mat_ind)) {
diff --git a/R/get-cut-points.R b/R/get-cut-points.R
index 168f5ab1..1f026be4 100644
--- a/R/get-cut-points.R
+++ b/R/get-cut-points.R
@@ -40,7 +40,8 @@ get_cut.default <- function(
 
 }
 
-
+#' @rdname get_cut
+#' @inherit get_cut
 get_cut.list <- function (
   data,
   formula,
diff --git a/R/ggplot-extensions.R b/R/ggplot-extensions.R
index 54fda7e6..5e6b9d2b 100644
--- a/R/ggplot-extensions.R
+++ b/R/ggplot-extensions.R
@@ -4,12 +4,15 @@
 #'
 #' \code{geom_stepribbon} is an extension of the \code{geom_ribbon}, and
 #' is optimized for Kaplan-Meier plots with pointwise confidence intervals
-#' or a confidence band.
+#' or a confidence band. The default \code{direction}-argument \code{"hv"} is
+#' appropriate for right-continuous step functions like the hazard rates etc
+#' returned by \code{pammtools}.
 #'
 #' @seealso
 #'   \code{\link[ggplot2]{geom_ribbon}} \code{geom_stepribbon}
 #'   inherits from \code{geom_ribbon}.
-#' @inheritParams ggplot2:::geom_ribbon
+#' @inheritParams ggplot2::geom_ribbon
+#' @inheritParams ggplot2::geom_step
 #' @examples
 #' library(ggplot2)
 #' huron <- data.frame(year = 1875:1972, level = as.vector(LakeHuron))
@@ -26,6 +29,7 @@ geom_stepribbon <- function(
   data        = NULL,
   stat        = "identity",
   position    = "identity",
+  direction   = "hv",
   na.rm       = FALSE,
   show.legend = NA,
   inherit.aes = TRUE, ...) {
@@ -38,7 +42,7 @@ geom_stepribbon <- function(
     position    = position,
     show.legend = show.legend,
     inherit.aes = inherit.aes,
-    params      = list(na.rm = na.rm, ... )
+    params      = list(na.rm = na.rm, direction = direction, ... )
   )
 
 }
@@ -48,21 +52,59 @@ geom_stepribbon <- function(
 #' @format NULL
 #' @usage NULL
 #' @export
-
 GeomStepribbon <- ggproto(
   "GeomStepribbon", GeomRibbon,
 
   extra_params = c("na.rm"),
 
-  draw_group = function(data, panel_scales, coord, na.rm = FALSE) {
+  draw_group = function(data, panel_scales, coord, na.rm = FALSE, direction = "hv") {
 
     if (na.rm) data <- data[complete.cases(data[c("x", "ymin", "ymax")]), ]
     data   <- rbind(data, data)
     data   <- data[order(data$x), ]
-    data$x <- c(data$x[2:nrow(data)], NA)
-    data   <- data[complete.cases(data["x"]), ]
-    GeomRibbon$draw_group(data, panel_scales, coord, na.rm = FALSE)
-
+    data   <- ggplot2_stairstep(data[complete.cases(data["x"]), ],
+                               direction = direction)
+    GeomRibbon$draw_group(data, panel_scales, coord, na.rm = na.rm)
   }
 
 )
+
+
+# code adapted from
+# https://github.com/tidyverse/ggplot2/blob/9741da5050f81b7b5c012c56d02f45fc93d68f89/R/geom-path.r#L320
+ggplot2_stairstep <- function(data, direction =  c("hv", "vh", "mid")) {
+  direction <- match.arg(direction)
+  data <- as.data.frame(data)[order(data$x), ]
+  n <- nrow(data)
+  if (n <= 1) {
+    return(data[0, , drop = FALSE])
+  }
+  if (direction == "vh") {
+    xs <- rep(1:n, each = 2)[-2 * n]
+    ys <- c(1, rep(2:n, each = 2))
+  }
+  if (direction == "hv") {
+    xs <- c(1, rep(2:n, each = 2))
+    ys <- rep(1:n, each = 2)[-2 * n]
+  }
+  if (direction == "mid") {
+    xs <- rep(1:(n - 1), each = 2)
+    ys <- rep(1:n, each = 2)
+  }
+
+  ymin <- c(data$ymin[ys])
+  ymax <- c(data$ymax[ys])
+  if (direction == "mid") {
+    gaps <- data$x[-1] - data$x[-n]
+    mid_x <- data$x[-n] + gaps/2
+    x <- c(data$x[1], mid_x[xs], data$x[n])
+    data_attr <- data[c(1, xs, n),
+                      setdiff(names(data), c("x", "ymin", "ymax"))]
+  } else {
+    x <- data$x[xs]
+    ymin <- data$ymin[ys]
+    ymax <- data$ymax[ys]
+    data_attr <- data[xs, setdiff(names(data), c("x", "ymin", "ymax"))]
+  }
+  cbind(data.frame(x = x, ymin = ymin, ymax = ymax), data_attr)
+}
diff --git a/R/lag-lead-utils.R b/R/lag-lead-utils.R
index 8e72bced..97e81834 100644
--- a/R/lag-lead-utils.R
+++ b/R/lag-lead-utils.R
@@ -103,8 +103,8 @@ gg_laglead.LL_df <- function(
       tz = factor(.data$tz, levels = sort(unique(.data$tz),
         decreasing = FALSE)),
       interval = factor(.data$interval, levels = levels(.data$interval)) )
-  gg_ll <- ggplot(x, aes_string(x = "interval", y = "tz")) +
-    geom_tile(aes_string(fill = "LL"), colour = grid_col) +
+  gg_ll <- ggplot(x, aes(x = .data[["interval"]], y = .data[["tz"]])) +
+    geom_tile(aes(fill = .data[["LL"]]), colour = grid_col) +
     scale_fill_gradient(low = low_col, high = high_col) +
     scale_x_discrete(expand = c(0, 0)) +
     scale_y_discrete(expand = c(0, 0)) +
diff --git a/R/make-newdata.R b/R/make-newdata.R
index 237a942e..c1af0212 100644
--- a/R/make-newdata.R
+++ b/R/make-newdata.R
@@ -83,7 +83,7 @@ sample_info.fped <- function(x) {
 #'
 #' @importFrom dplyr slice bind_cols
 #' @importFrom vctrs vec_c
-#' @importFrom purrr map map_lgl map2 transpose cross
+#' @importFrom purrr map map_lgl
 #' @importFrom checkmate test_data_frame
 #' @param ... Data frames that should be combined to one data frame.
 #' Elements of first df vary fastest, elements of last df vary slowest.
@@ -102,10 +102,10 @@ combine_df <- function(...) {
   }
   ind_seq   <- map(dots, ~ seq_len(nrow(.x)))
   not_empty <- map_lgl(ind_seq, ~ length(.x) > 0)
-  ind_list  <- ind_seq[not_empty] %>% cross() %>% transpose() %>% map(function(x) vec_c(!!!x))
 
-  map2(dots[not_empty], ind_list, function(.x, .y) slice(.x, .y)) %>%
-    bind_cols()
+  ord <- lapply(dots[not_empty], function(z) colnames(z)) |> unlist()
+  out <- do.call(expand_grid, rev(dots[not_empty]))
+  out <- out[, ord]
 
 }
 
@@ -121,6 +121,7 @@ combine_df <- function(...) {
 #'
 #' @rdname newdata
 #' @aliases make_newdata
+#' @importFrom tidyr expand_grid
 #' @inheritParams sample_info
 #' @param ... Covariate specifications (expressions) that will be evaluated
 #' by looking for variables in \code{x}. Must be of the form \code{z = f(z)}
@@ -131,7 +132,8 @@ combine_df <- function(...) {
 #' the time argument, but see "Details" an "Examples" below.
 #' @import dplyr
 #' @importFrom checkmate assert_data_frame assert_character
-#' @importFrom purrr map cross_df
+#' @importFrom purrr map
+#' @importFrom tidyr expand_grid
 #' @details Depending on the type of variables in \code{x}, mean or modus values
 #' will be used for variables not specified in ellipsis
 #' (see also \code{\link[pammtools]{sample_info}}). If \code{x} is an object
@@ -159,9 +161,6 @@ combine_df <- function(...) {
 #' # mean/modus values of unspecified variables are calculated over whole data
 #' tumor %>% make_newdata(sex=unique(sex))
 #' tumor %>% group_by(sex) %>% make_newdata()
-#' # You can also pass a part of the data sets as data frame to make_newdata
-#' purrr::cross_df(list(days = c(0, 500, 1000), sex = c("male", "female"))) %>%
-#'   make_newdata(x=tumor)
 #'
 #' # Examples for PED data
 #' ped <- tumor %>% slice(1:3) %>% as_ped(Surv(days, status)~., cut = c(0, 500, 1000))
@@ -193,11 +192,13 @@ make_newdata.default <- function(x, ...) {
   orig_names <- names(x)
 
   expressions    <- quos(...)
-  expr_evaluated <- map(expressions, lazyeval::f_eval, data = x)
+  expr_evaluated <- map(expressions, lazyeval::f_eval, data = x) |>
+    map(c)
 
   # construct data parts depending on input type
   lgl_atomic <- map_lgl(expr_evaluated, is_atomic)
-  part1 <- expr_evaluated[lgl_atomic] %>% cross_df()
+  # part1 <- expr_evaluated[lgl_atomic] |> cross_df()
+  part1 <- do.call(tidyr::expand_grid, rev(expr_evaluated[lgl_atomic]))
   part2 <- do.call(combine_df, expr_evaluated[!lgl_atomic])
 
   ndf  <- combine_df(part1, part2)
diff --git a/R/model-evaluation.R b/R/model-evaluation.R
index b356e146..9be594c0 100644
--- a/R/model-evaluation.R
+++ b/R/model-evaluation.R
@@ -16,7 +16,7 @@ as.data.frame.crps <- function(x, row.names = NULL, optional = FALSE, ...) {
   m$method <- attr(x, "dimnames")[[1]]
 
   m <- m %>%
-    pivot_longer(cols = -.data$method, values_to = "IBS") %>%
-    dplyr::rename(time = .data$name)
+    pivot_longer(cols = -one_of("method"), values_to = "IBS") %>%
+    dplyr::rename(time = "name")
 
 }
diff --git a/R/multi-state-helpers.R b/R/multi-state-helpers.R
deleted file mode 100644
index 8a4fbcc2..00000000
--- a/R/multi-state-helpers.R
+++ /dev/null
@@ -1,130 +0,0 @@
-#' @examples
-#' t_mat <- matrix(data = NA, nrow = 4, ncol = 4)
-#' t_mat[1,2] <- "-1 - sin(x2)"
-#' t_mat[1,3] <- "-1.3 + 0.1*x2*x3"
-#' t_mat[2,3] <- "-1.5 + 0.8*x1 + 0.3*x2 - 0.1*x3"
-#' @keywords internal
-from_to_pairs <- function(t_mat) {
-
-
-
-}
-
-#' @examples
-#' ftp <- from_to_pairs2(t_mat)
-from_to_pairs2 <- function(t_mat) {
-
-  res <- apply(t_mat, 1, function(x) which(x) - 1)
-  names(res) <- seq_len(nrow(t_mat)) - 1
-  res <- res[vapply(res, length, 0) != 0]
-
-  res
-
-}
-
-from_to_pairs.data.frame <- function(data, from_col = "from", to_col = "to") {
-
-  map(
-    .x = sort(unique(data[[from_col]])),
-    .f = ~{
-      data %>%
-        filter(.data[[from_col]] == .x) %>%
-        pull(to_col) %>%
-        unique()
-    })
-
-}
-
-
-#' Add counterfactual observations for possible transitions
-#'
-#' If data only contains one row per transition that took place, this function
-#' adds additional rows for each transition that was possible at that time
-#' (for each subject in the data).
-#' @param data Data set that only contains rows for transitions that took place.
-#' @param from_to_pairs A list with one element for each possible initial state.
-#' The values of each list element indicate possible transitions from that state.
-#' Will be calculated from the data if unspecified.
-#' @param from_col Name of the column that stores initial state.
-#' @param to_col Name of the column that stores end state.
-#' @param transition_col Name of the column that contains the transition identifier (factor variable).
-#' @examples
-#'
-#' res <- add_counterfactual_transitions(data, list(c(2,3), 3))
-#' res <- add_counterfactual_transitions(data)
-#' res %>% select(id, from, to, status, tstart, tstop, transition, trans)
-#' @export
-add_counterfactual_transitions <- function(
-  data,
-  from_to_pairs  = list(),
-  from_col       = "from",
-  to_col         = "to",
-  transition_col = "transition") {
-
-
-  if(length(from_to_pairs) == 0) {
-    from_to_pairs <- from_to_pairs.data.frame(data, from_col, to_col) %>%
-      discard(~length(.x) == 0)
-  }
-
-  l_from <- split(data, data[[from_col]])
-
-  orig_status <- data %>%
-    select(id, from, to, tstart, status) %>%
-    rename(
-      orig_to = to,
-      orig_status = status)
-
-  res <- map2_dfr(
-    from_to_pairs,
-    l_from, ~{
-      n_to <- length(.x)
-      .y %>%
-        group_by(id) %>%
-        mutate(initial_id = seq_len(n())) %>%
-        ungroup() %>%
-        slice(rep(row_number(), n_to)) %>%
-        arrange(id, from, to) %>%
-        group_by(id, initial_id) -> temp
-        temp %>%
-        mutate(
-          to = .x,
-          transition = paste0(from, "->", to)
-        ) %>%
-        ungroup() %>%
-        mutate(
-          initial_id = NULL,
-          transition = as.factor(transition))
-    }
-  )
-
-  res %>%
-    ungroup() %>%
-    left_join(orig_status) %>%
-    mutate(status = status * (to == orig_to)) %>%
-    select(-orig_status, -orig_to) %>%
-    arrange(id, tstart, tstop, from, to, status)
-
-
-}
-
-# ped_msm <- res %>% as_ped(
-#   formula = Surv(tstart, tstop, status)~.,
-#   transition = "transition",
-#   timescale = "calendar")
-# data %>%
-#   select(id, tstart, tstop, status, transition)
-# res %>%
-#   select(id, tstart, tstop, status, from, to, transition)
-# ped_msm %>%
-#   group_by(id, transition) %>%
-#   slice(1,n()) %>%
-#   select(id, tstart, tend, ped_status, transition) %>%
-#   filter(id %in% c(1, 4, 5))
-
-# split_df  %>%
-#   mutate(row_number = row_number()) %>%
-#   group_by(id, transition) %>%
-#   slice(1,n()) %>%
-#   select(row_number, id, tstart, tend, ped_status, transition) %>%
-#   filter(id %in% c(1, 4, 5))
diff --git a/R/nest-utils.R b/R/nest-utils.R
index b92cae63..f66603ce 100644
--- a/R/nest-utils.R
+++ b/R/nest-utils.R
@@ -22,8 +22,6 @@ nest_tdc <- function(data, formula, ...) {
 }
 
 #' @rdname nest_tdc
-#' @param vars A character vector of TDCs that will be nested.
-#' @param id A character giving the name of the ID column.
 #' @export
 nest_tdc.default <- function(data, formula, ...) {
 
diff --git a/R/pammtools.R b/R/pammtools.R
index 7c335d7d..7a7d00af 100644
--- a/R/pammtools.R
+++ b/R/pammtools.R
@@ -38,7 +38,8 @@
 #' }
 #'
 #' @name pammtools
-#' @docType package
+#' @keywords internal
+"_PACKAGE"
 #' @references
 #' Bender, Andreas, Andreas Groll, and Fabian Scheipl. 2018.
 #' “A Generalized Additive Model Approach to Time-to-Event Analysis”
diff --git a/R/predict.R b/R/predict.R
index 772d3f7f..0b18aae1 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -8,7 +8,8 @@
 predictSurvProb.pamm <- function(
   object,
   newdata,
-  times) {
+  times,
+  ...) {
 
   if (!is.ped(newdata)) {
 
diff --git a/R/sim-pexp.R b/R/sim-pexp.R
index d0f0c09a..eab54cfd 100644
--- a/R/sim-pexp.R
+++ b/R/sim-pexp.R
@@ -155,9 +155,9 @@ sim_pexp <- function(formula, data, cut) {
       status = 1L * (.data$time <= max(cut)),
       time   = pmin(.data$time, max(cut)))
 
-    suppressMessages(
+  suppressMessages(
     sim_df <- sim_df %>%
-      left_join(select(data, -.data$time, -.data$status))
+      left_join(select(data, -all_of(c("time", "status"))))
   )
 
   attr(sim_df, "id_var")     <- "id"
@@ -166,7 +166,7 @@ sim_pexp <- function(formula, data, cut) {
   attr(sim_df, "tz_var")     <- tz_vars
   attr(sim_df, "cens_value") <- 0
   attr(sim_df, "breaks")     <- cut
-  attr(sim_df, "tz")         <- imap(tz_vars, ~select(sim_df, .x) %>%
+  attr(sim_df, "tz")         <- imap(tz_vars, ~select(sim_df, all_of(.x)) %>%
     pull(.x) %>% unique()) %>% flatten()
   if (exists("ll_funs")) attr(sim_df, "ll_funs") <- ll_funs
   if (exists("cumu_funs")) attr(sim_df, "cumu_funs") <- cumu_funs
@@ -314,7 +314,7 @@ sim_pexp_cr <- function(formula, data, cut) {
     cut     = cut,
     id      = "id") %>%
     mutate(
-      t = t + tstart
+      t = t + .data$tstart
     )
 
   # calculate cause specific hazards
@@ -348,149 +348,3 @@ sim_pexp_cr <- function(formula, data, cut) {
     select(-one_of(c("tstart", "tend", "interval", "offset", "ped_status", "rate")))
 
 }
-
-
-#' Simulate data for multi-state models
-#'
-#' @rdname sim_pexp
-#' @inherit sim_pexp
-#' @param tmat Transition matrix containing the log-hazard definition for each transition.
-#' Transitions that can not occurr should contain a \code{NA}.
-#' @param proportion_censoring The proportion of censored observations
-#' @param keep_transitions_at_risk Logical. If \code{TRUE} (default), all transitions
-#' for which a subject was at risk will be included in the data set.
-#' @export
-#' @keywords internal
-#' @examples
-#' t_mat <- matrix(data = NA, nrow = 3, ncol = 3)
-#' t_mat[1,2] <- "log(0.7) + x1"
-#' t_mat[1,3] <- "log(0.5) - 0.25 * x2"
-#' t_mat[2,3] <- "log(0.9)"
-#'
-#' n = 100
-#' data <- cbind.data.frame(
-#'  id = 1:n,
-#'  x1 = runif(n, -3, 3),
-#'  x2 = runif(n, 0, 6),
-#'  from = 1,
-#'  t = 0)
-#'
-#' cut =  seq(0, 3, by = 0.01)
-#'
-#' msm_df <- sim_pexp_msm(
-#'  t_mat = t_mat,
-#'  data = data,
-#'  cut = cut,
-#'  proportoin_censoring = 0.3,
-#'  keep_transitions_at_risk = TRUE)
-#' head(msm_df)
-#'
-sim_pexp_msm <- function(
-  t_mat,
-  data,
-  cut,
-  proportion_censoring = 0.3,
-  keep_transitions_at_risk = FALSE ) {
-
-  # Aus der Transition Matrix können die Transition-States abgelesen werden
-  trans_states <- which(apply(t_mat,1, function(x) !all(is.na(x))))
-
-  # cens_time <- c(
-  #   rep(max(cut), times = floor( (1 - proportoin_censoring)*nrow(data))),
-  #   runif((nrow(data) - floor((1 - proportoin_censoring)*nrow(data))), 0, max(cut))
-  # )
-  # data$cens_time <- sample(cens_time)
-
-  # Initialize object where results are stored
-  results <- NULL
-  sim_df <- NULL
-  trans_not_made <- NULL
-
-
-  while (nrow(data) > 0) {
-
-    # Für jeden Transition state
-    for (i in seq_along(trans_states)) {
-      # Welche Beobachtungen befinden sich in diesem Status
-      if (nrow(data[data$from == i,]) == 0) next
-      # Welche Hazards führen von diesem Status weg
-      form <- as.Formula(paste("~ ", paste(t_mat[i,which(!is.na(t_mat[i,]))], collapse = "|")))
-      # Führe für jede der Beobachtungen in dem state ein CR Experiment durch
-      # Speichere die Ergebnisse in einem Dataframe
-      results <- rbind(results,
-        sim_pexp_cr(form, data[data$from == i,], cut = cut) %>%
-          mutate(
-            time = time + t,
-            to =  which(!is.na(t_mat[i,]))[type]
-          ) %>%
-          select(-one_of(paste0("hazard",seq_len(length(attr(form,"rhs"))))))
-      )
-
-
-      if (keep_transitions_at_risk) {
-        # Alle Transitions, die nicht gemacht wurden, aber theoretisch möglich gewesen wären, sollen als zensiert im Dataframe auftauchen (Vgl. mstate)
-        trans_not_made <- rbind(trans_not_made,
-          results %>%
-            slice(rep(row_number(),length(attr(form, "rhs")) )) %>%
-            mutate(
-              pos_types = seq_len(length(attr(form, "rhs"))),
-              status = if_else(type != pos_types, 0, 1),
-              to =  which(!is.na(t_mat[i,]))[pos_types]
-            ) %>%
-            filter(
-              type != pos_types
-            ) %>%
-            mutate(
-              pos_types = NULL
-            )
-
-        )
-
-      }
-
-    }
-
-    sim_df <- rbind(sim_df, results, trans_not_made)
-
-    data <- data %>%
-      left_join(
-        results[,c("id", "to", "time","status")], by = "id") %>%
-      filter(
-        status != 0 & to %in% trans_states) %>%
-      mutate(
-        from   = to,
-        to     = NULL,
-        t      = t + time,
-        time   = NULL,
-        status = NULL
-      )
-    # results und trans_not_made wieder leeren
-    results <- NULL
-    trans_not_made <- NULL
-
-    # Ergebnisframe sortieren, ids untereinander und nach Transitionzeiten sortieren
-    sim_df <- sim_df %>%
-      group_by(id) %>%
-      arrange(time, .by_group = TRUE)
-    # Solange, wie sich Individuen in data befinden gibt es id's in transient States, und das Vorgehen wird wiederholt
-  }
-
-  # sim_df aufräumen
-  sim_df <- sim_df %>%
-    mutate(
-      status = status * (time <= max(cut)),
-      time   = pmin(time, max(cut)),
-      # to     = ifelse(status, to, from),
-      tstart = t,
-      tstop  = time,
-      gap    = tstop - tstart,
-      time   = NULL,
-      t      = NULL,
-      type   = NULL,
-      transition = as.factor(paste0(from, "->", to))
-    ) %>%
-    relocate(to, .after = from)
-
-  return(sim_df)
-
-}
diff --git a/R/tidiers.R b/R/tidiers.R
index db113897..7fef4c2e 100644
--- a/R/tidiers.R
+++ b/R/tidiers.R
@@ -46,7 +46,7 @@ tidy_fixed.gam <- function(x, intercept=FALSE, ...) {
 
   ftab <- summary(x)[["p.table"]][, 1:2]
   if (!intercept) {
-    ftab <- ftab[!grepl("Intercept", rownames(ftab)), ]
+    ftab <- ftab[!grepl("Intercept", rownames(ftab)), , drop = FALSE]
   }
   calc_ci(ftab)
 
diff --git a/R/tidyverse-methods.R b/R/tidyverse-methods.R
index ea664c57..2cbf674e 100644
--- a/R/tidyverse-methods.R
+++ b/R/tidyverse-methods.R
@@ -66,9 +66,7 @@ fped_attr <- function(fped) {
 #' @param .data an  object of class \code{ped}, see \code{\link{as_ped}}.
 #' @param tbl an  object of class \code{ped}, see \code{\link{as_ped}}.
 #' @param x an  object of class \code{ped}, see \code{\link{as_ped}}.
-#' @param funs see \code{\link[dplyr]{summarize_all}}
 #' @param ... see \code{dplyr} documentation
-#' @param .dots see \code{dplyr} documentation
 #' @description See \code{dplyr} documentation of the respective functions for
 #'   description and examples.
 #' @return a modified \code{ped} object (except for \code{do})
@@ -218,7 +216,6 @@ select.ped <- function(.data, ...) {
 
 }
 
-#' @param keep_attributes conserve attributes? defaults to \code{TRUE}
 #' @export
 #' @export mutate
 #' @rdname dplyr_verbs
diff --git a/R/viz-elra.R b/R/viz-elra.R
index f3494a61..209a82a6 100644
--- a/R/viz-elra.R
+++ b/R/viz-elra.R
@@ -32,16 +32,16 @@ gg_partial <- function(data, model, term, ..., reference = NULL, ci = TRUE) {
   # ndf      <- ndf %>% mutate_at(vars[-1], ~as.factor(.x))
   n_vars   <- length(vars)
 
-  gg_base <- ggplot(ndf, aes_string(x = vars[1])) + xlab(vars[1])
+  gg_base <- ggplot(ndf, aes(x = .data[[vars[1]]])) + xlab(vars[1])
   if (n_vars == 1) {
     gg_out <- gg_base +
-      geom_ribbon(aes_string(ymin = "ci_lower", ymax = "ci_upper"),
+      geom_ribbon(aes(ymin = .data[["ci_lower"]], ymax = .data[["ci_upper"]]),
         alpha = 0.3) +
-      geom_line(aes_string(y = "fit"))
+      geom_line(aes(y = .data[["fit"]]))
   } else {
     # if (n_vars == 2) {
-      gg_out <- gg_base + aes_string(y = vars[2], z = "fit") +
-        geom_tile(aes_string(fill = "fit")) +
+      gg_out <- gg_base + aes(y = .data[[vars[2]]], z = .data[["fit"]]) +
+        geom_tile(aes(fill = .data[["fit"]])) +
         geom_contour(col = "grey30") +
         scale_y_continuous(expand = c(0, 0)) +
         scale_x_continuous(expand = c(0, 0)) +
@@ -94,8 +94,8 @@ gg_partial_ll <- function(
         levels = c("ci_lower", "fit", "ci_upper")))
   }
 
-  gg_base <- ggplot(ll_df, aes_string(x = "intmid", y =  tz_var)) +
-    geom_tile(aes_string(fill = "fit"), colour = "grey30") +
+  gg_base <- ggplot(ll_df, aes(x = .data[["intmid"]], y =  tz_var)) +
+    geom_tile(aes(fill = .data[["fit"]]), colour = "grey30") +
     scale_fill_gradient2(high = "firebrick2", low = "steelblue",
       na.value = "grey30") +
     scale_x_continuous(expand = c(0, 0)) +
@@ -142,19 +142,19 @@ gg_slice <- function(data, model, term, ..., reference = NULL, ci = TRUE) {
   ndf      <- ndf %>% mutate_at(vars[-1], ~as.factor(.x))
   n_vars   <- length(vars)
 
-  gg_out <- ggplot(ndf, aes_string(x = vars[1], y = "fit"))
+  gg_out <- ggplot(ndf, aes(x = .data[[vars[1]]], y = .data[["fit"]]))
   if (ci) {
     gg_out <- gg_out +
-      geom_ribbon(aes_string(ymin = "ci_lower", ymax = "ci_upper"), alpha = 0.3)
+      geom_ribbon(aes(ymin = .data[["ci_lower"]], ymax = .data[["ci_upper"]]), alpha = 0.3)
   }
   gg_out <- gg_out + geom_line()
   if (n_vars > 1) {
     if(ci) {
-      gg_out <- gg_out + aes_string(group = vars[2], fill = vars[2]) +
-      geom_line(aes_string(col = vars[2]))
+      gg_out <- gg_out + aes(group = .data[[vars[2]]], fill = .data[[vars[2]]]) +
+      geom_line(aes(col = .data[[vars[2]]]))
     } else {
-      gg_out <- gg_out + aes_string(group = vars[2]) +
-        geom_line(aes_string(col = vars[2]))
+      gg_out <- gg_out + aes(group = .data[[vars[2]]]) +
+        geom_line(aes(col = .data[[vars[2]]]))
     }
     if (n_vars > 2) {
       form   <- as.formula(paste0("~", vars[-1:-2], collapse = "+"))
@@ -176,10 +176,10 @@ gg_cumu_eff <- function(data, model, term, z1, z2=NULL, se_mult = 2, ci = TRUE)
 
   cumu_eff_df <- get_cumu_eff(data, model, term, z1, z2, se_mult)
 
-  gg_out <- ggplot(cumu_eff_df, aes_string(x = "tend", y = "cumu_eff"))
+  gg_out <- ggplot(cumu_eff_df, aes(x = .data[["tend"]], y = .data[["cumu_eff"]]))
   if (ci) {
     gg_out <- gg_out +
-      geom_ribbon(aes_string(ymin = "cumu_eff_lower", ymax = "cumu_eff_upper"),
+      geom_ribbon(aes(ymin = .data[["cumu_eff_lower"]], ymax = .data[["cumu_eff_upper"]]),
         alpha = 0.3)
   }
 
diff --git a/R/warnings.R b/R/warnings.R
index acd06f91..899705c7 100644
--- a/R/warnings.R
+++ b/R/warnings.R
@@ -7,7 +7,7 @@ warn_about_new_time_points <- function(object, newdata, ...) {
 
 }
 
-
+#' @inherit warn_about_new_time_points
 warn_about_new_time_points.glm <- function(object, newdata, time_var, ...) {
 
   is_pam <- (inherits(object, "gam" ) | inherits( object, "scam"))
@@ -41,7 +41,7 @@ warn_about_new_time_points.pamm <- function(object, newdata, ...) {
     int_original <- int_info(object)
     if ("interval" %in% colnames(newdata)) {
       int_new <- unique(newdata[["interval"]])
-      if(!all(int_new %in% int_original)) {
+      if(!all(int_new %in% int_original[["interval"]])) {
         warning(
           paste0(
             "Time points/intervals in new data not equivalent to time points/intervals during model fit.",
diff --git a/README.Rmd b/README.Rmd
index 230df33f..cf2fa7fb 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -9,19 +9,18 @@ opts_chunk$set(warning = FALSE)
 
 <!-- badges: start -->
 
-[![R-CMD-check](https://github.com/adibender/pammtools/workflows/R-CMD-check/badge.svg)](https://github.com/adibender/pammtools/actions)
+[![R-CMD-check](https://github.com/adibender/pammtools/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/adibender/pammtools/actions/workflows/R-CMD-check.yaml)
 [![cran
 checks](https://cranchecks.info/badges/worst/pammtools)](https://cran.r-project.org/web/checks/check_results_pammtools.html)
 [![lifecycle](https://lifecycle.r-lib.org/articles/figures/lifecycle-stable.svg)](https://lifecycle.r-lib.org/articles/stages.html)
 [![Build
 Status](https://travis-ci.org/adibender/pammtools.svg?branch=master)](https://travis-ci.org/adibender/pammtools)
-[![Build
 Status](https://ci.appveyor.com/api/projects/status/github/adibender/pammtools?branch=master&svg=true)](https://ci.appveyor.com/project/adibender/pammtools/branch/master)
 [![codecov.io](https://codecov.io/github/adibender/pammtools/coverage.svg?branch=master)](https://codecov.io/github/adibender/pammtools/branch/master)
 [![CRAN_Status_Badge](https://www.r-pkg.org/badges/version-ago/pammtools)](https://cran.r-project.org/package=pammtools)
 [![CRAN_Download_Badge](https://cranlogs.r-pkg.org/badges/pammtools)](https://cran.r-project.org/package=pammtools)
 [![MIT
-license](http://img.shields.io/badge/license-MIT-brightgreen.svg)]( https://opensource.org/licenses/MIT)
+license](https://img.shields.io/badge/license-MIT-brightgreen.svg)]( https://opensource.org/license/mit)
 <!-- badges: end -->
 
 # **`pammtools`**: Piece-Wise Exponential Additive Mixed Modeling Tools
diff --git a/README.md b/README.md
index e6d14eea..608e2d6f 100644
--- a/README.md
+++ b/README.md
@@ -3,17 +3,17 @@
 
 [![R-CMD-check](https://github.com/adibender/pammtools/workflows/R-CMD-check/badge.svg)](https://github.com/adibender/pammtools/actions)
 [![cran
-checks](https://cranchecks.info/badges/worst/pammtools)](https://cran.r-project.org/web/checks/check_results_pammtools.html)
+checks](https://badges.cranchecks.info/worst/pammtools.svg)](https://cran.r-project.org/web/checks/check_results_pammtools.html)
 [![lifecycle](https://lifecycle.r-lib.org/articles/figures/lifecycle-stable.svg)](https://lifecycle.r-lib.org/articles/stages.html)
 [![Build
-Status](https://travis-ci.org/adibender/pammtools.svg?branch=master)](https://travis-ci.org/adibender/pammtools)
+Status](https://travis-ci.org/adibender/pammtools.svg?branch=master)](https://app.travis-ci.com/adibender/pammtools)
 [![Build
 Status](https://ci.appveyor.com/api/projects/status/github/adibender/pammtools?branch=master&svg=true)](https://ci.appveyor.com/project/adibender/pammtools/branch/master)
-[![codecov.io](https://codecov.io/github/adibender/pammtools/coverage.svg?branch=master)](https://codecov.io/github/adibender/pammtools/branch/master)
+[![codecov.io](https://codecov.io/github/adibender/pammtools/coverage.svg?branch=master)](https://app.codecov.io/github/adibender/pammtools/branch/master)
 [![CRAN_Status_Badge](https://www.r-pkg.org/badges/version-ago/pammtools)](https://cran.r-project.org/package=pammtools)
 [![CRAN_Download_Badge](https://cranlogs.r-pkg.org/badges/pammtools)](https://cran.r-project.org/package=pammtools)
 [![MIT
-license](http://img.shields.io/badge/license-MIT-brightgreen.svg)]( https://opensource.org/licenses/MIT)
+license](https://img.shields.io/badge/license-MIT-brightgreen.svg)]( https://opensource.org/license/mit)
 <!-- badges: end -->
 
 # **`pammtools`**: Piece-Wise Exponential Additive Mixed Modeling Tools
diff --git a/cran-comments.md b/cran-comments.md
index 1605fd09..3f68afd9 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,5 +1,5 @@
 ## Test environments
-* Local Ubuntu Linux R 4.1.2
+* Local Ubuntu Linux R 4.2.1
 * R release with Travis CI
 * Appveyor (patched)
 * win-builder
@@ -12,4 +12,4 @@
 No problems.
 
 ## Comments
-Some minor updates
+- Maintenance updates
diff --git a/inst/CITATION b/inst/CITATION
index 490cb6e1..cdc267c5 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -1,8 +1,9 @@
 citHeader("Citation information for 1) the Software 2) the tutorial on PAMMs 3)
-paper on cumulative effects/exposure-lag-response associations:")
+paper on cumulative effects/exposure-lag-response associations, 4) Recurrent events:")
 
 
-citEntry(entry = "Article",
+bibentry(
+  bibtype = "Article",
   author = "Andreas Bender and Fabian Scheipl",
   title = "pammtools: Piece-wise exponential Additive Mixed Modeling tools",
   year = "2018",
@@ -11,21 +12,30 @@ citEntry(entry = "Article",
   exponential additive mixed modeling tools. arXiv:1806.01042 [stat]",
   url = "https://arxiv.org/abs/1806.01042")
 
-citEntry(
-  entry = "Article",
+bibentry(
+  bibtype = "Article",
   author = "Andreas Bender and Andreas Groll and Fabian Scheipl",
   title = "A generalized additive model approach to time-to-event analysis",
   textVersion = "Andreas Bender, Andreas Groll and Fabian Scheipl. 2018. A Generalized Additive Model Approach to Time-to-Event Analysis. Statistical Modelling.",
   journal = "Statistical Modelling",
   year = "2018",
-  url = "https://journals.sagepub.com/doi/10.1177/1471082X17748083"
+  url = "https://doi.org/10.1177/1471082X17748083"
 )
 
-citEntry(
-  entry = "Article",
+bibentry(
+  bibtype = "Article",
   author = "Bender, Andreas and Scheipl, Fabian and Hartl, Wolfgang and Day, Andrew G and K\u00fcchenhoff, Helmut",
   title = "Penalized estimation of complex, non-linear exposure-lag-response associations",
   textVersion = "Andreas Bender, Fabian Scheipl, Wolfgang Hartl, Andrew G Day, Helmut K\u00fcchenhoff; Penalized estimation of complex, non-linear exposure-lag-response associations, Biostatistics, , kxy003",
   journal = "Biostatistics",
   year = "2019"
 )
+
+bibentry(
+  bibtype = "Article",
+  author = "Ramjith Jordache, and Bender, Andreas and Roes, Kit C. B. and Jonker, Marianne A.",
+  title = "Recurrent events analysis with piece-wise exponential additive mixed models",
+  textVersion = "Ramjith J, Bender A, Roes KCB, Jonker MA. Recurrent events analysis with piece-wise exponential additive mixed models. 2022. Statistical Modelling., 2022",
+  journal = "Statistical Modelling",
+  year = "2022"
+)
diff --git a/man/add_counterfactual_transitions.Rd b/man/add_counterfactual_transitions.Rd
deleted file mode 100644
index 3f1bd012..00000000
--- a/man/add_counterfactual_transitions.Rd
+++ /dev/null
@@ -1,38 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/multi-state-helpers.R
-\name{add_counterfactual_transitions}
-\alias{add_counterfactual_transitions}
-\title{Add counterfactual observations for possible transitions}
-\usage{
-add_counterfactual_transitions(
-  data,
-  from_to_pairs = list(),
-  from_col = "from",
-  to_col = "to",
-  transition_col = "transition"
-)
-}
-\arguments{
-\item{data}{Data set that only contains rows for transitions that took place.}
-
-\item{from_to_pairs}{A list with one element for each possible initial state.
-The values of each list element indicate possible transitions from that state.
-Will be calculated from the data if unspecified.}
-
-\item{from_col}{Name of the column that stores initial state.}
-
-\item{to_col}{Name of the column that stores end state.}
-
-\item{transition_col}{Name of the column that contains the transition identifier (factor variable).}
-}
-\description{
-If data only contains one row per transition that took place, this function
-adds additional rows for each transition that was possible at that time
-(for each subject in the data).
-}
-\examples{
-
-res <- add_counterfactual_transitions(data, list(c(2,3), 3))
-res <- add_counterfactual_transitions(data)
-res \%>\% select(id, from, to, status, tstart, tstop, transition, trans)
-}
diff --git a/man/as_ped.Rd b/man/as_ped.Rd
index c542bbbc..ab15a9cf 100644
--- a/man/as_ped.Rd
+++ b/man/as_ped.Rd
@@ -88,7 +88,7 @@ Often this will be \code{0}, which is the default.}
 \item{x}{any R object.}
 
 \item{newdata}{A new data set (\code{data.frame}) that contains the same
-variables that were used to create the PED object (code{data}).}
+variables that were used to create the PED object (\code{data}).}
 }
 \value{
 A data frame class \code{ped} in piece-wise exponential data format.
diff --git a/man/dplyr_verbs.Rd b/man/dplyr_verbs.Rd
index b68fd11d..ca8e9d2b 100644
--- a/man/dplyr_verbs.Rd
+++ b/man/dplyr_verbs.Rd
@@ -134,12 +134,6 @@ it is a potentially expensive operation so you must opt into it.}
 \item{suffix}{If there are non-joined duplicate variables in \code{x} and
 \code{y}, these suffixes will be added to the output to disambiguate them.
 Should be a character vector of length 2.}
-
-\item{funs}{see \code{\link[dplyr]{summarize_all}}}
-
-\item{.dots}{see \code{dplyr} documentation}
-
-\item{keep_attributes}{conserve attributes? defaults to \code{TRUE}}
 }
 \value{
 a modified \code{ped} object (except for \code{do})
diff --git a/man/geom_stepribbon.Rd b/man/geom_stepribbon.Rd
index 83f814dd..4207c3f3 100644
--- a/man/geom_stepribbon.Rd
+++ b/man/geom_stepribbon.Rd
@@ -11,16 +11,112 @@ geom_stepribbon(
   data = NULL,
   stat = "identity",
   position = "identity",
+  direction = "hv",
   na.rm = FALSE,
   show.legend = NA,
   inherit.aes = TRUE,
   ...
 )
 }
+\arguments{
+\item{mapping}{Set of aesthetic mappings created by \code{\link[ggplot2:aes]{aes()}}. If specified and
+\code{inherit.aes = TRUE} (the default), it is combined with the default mapping
+at the top level of the plot. You must supply \code{mapping} if there is no plot
+mapping.}
+
+\item{data}{The data to be displayed in this layer. There are three
+options:
+
+If \code{NULL}, the default, the data is inherited from the plot
+data as specified in the call to \code{\link[ggplot2:ggplot]{ggplot()}}.
+
+A \code{data.frame}, or other object, will override the plot
+data. All objects will be fortified to produce a data frame. See
+\code{\link[ggplot2:fortify]{fortify()}} for which variables will be created.
+
+A \code{function} will be called with a single argument,
+the plot data. The return value must be a \code{data.frame}, and
+will be used as the layer data. A \code{function} can be created
+from a \code{formula} (e.g. \code{~ head(.x, 10)}).}
+
+\item{stat}{The statistical transformation to use on the data for this layer.
+When using a \verb{geom_*()} function to construct a layer, the \code{stat}
+argument can be used the override the default coupling between geoms and
+stats. The \code{stat} argument accepts the following:
+\itemize{
+\item A \code{Stat} ggproto subclass, for example \code{StatCount}.
+\item A string naming the stat. To give the stat as a string, strip the
+function name of the \code{stat_} prefix. For example, to use \code{stat_count()},
+give the stat as \code{"count"}.
+\item For more information and other ways to specify the stat, see the
+\link[ggplot2:layer_stats]{layer stat} documentation.
+}}
+
+\item{position}{A position adjustment to use on the data for this layer. This
+can be used in various ways, including to prevent overplotting and
+improving the display. The \code{position} argument accepts the following:
+\itemize{
+\item The result of calling a position function, such as \code{position_jitter()}.
+This method allows for passing extra arguments to the position.
+\item A string naming the position adjustment. To give the position as a
+string, strip the function name of the \code{position_} prefix. For example,
+to use \code{position_jitter()}, give the position as \code{"jitter"}.
+\item For more information and other ways to specify the position, see the
+\link[ggplot2:layer_positions]{layer position} documentation.
+}}
+
+\item{direction}{direction of stairs: 'vh' for vertical then horizontal,
+'hv' for horizontal then vertical, or 'mid' for step half-way between
+adjacent x-values.}
+
+\item{na.rm}{If \code{FALSE}, the default, missing values are removed with
+a warning. If \code{TRUE}, missing values are silently removed.}
+
+\item{show.legend}{logical. Should this layer be included in the legends?
+\code{NA}, the default, includes if any aesthetics are mapped.
+\code{FALSE} never includes, and \code{TRUE} always includes.
+It can also be a named logical vector to finely select the aesthetics to
+display.}
+
+\item{inherit.aes}{If \code{FALSE}, overrides the default aesthetics,
+rather than combining with them. This is most useful for helper functions
+that define both data and aesthetics and shouldn't inherit behaviour from
+the default plot specification, e.g. \code{\link[ggplot2:borders]{borders()}}.}
+
+\item{...}{Other arguments passed on to \code{\link[ggplot2:layer]{layer()}}'s \code{params} argument. These
+arguments broadly fall into one of 4 categories below. Notably, further
+arguments to the \code{position} argument, or aesthetics that are required
+can \emph{not} be passed through \code{...}. Unknown arguments that are not part
+of the 4 categories below are ignored.
+\itemize{
+\item Static aesthetics that are not mapped to a scale, but are at a fixed
+value and apply to the layer as a whole. For example, \code{colour = "red"}
+or \code{linewidth = 3}. The geom's documentation has an \strong{Aesthetics}
+section that lists the available options. The 'required' aesthetics
+cannot be passed on to the \code{params}. Please note that while passing
+unmapped aesthetics as vectors is technically possible, the order and
+required length is not guaranteed to be parallel to the input data.
+\item When constructing a layer using
+a \verb{stat_*()} function, the \code{...} argument can be used to pass on
+parameters to the \code{geom} part of the layer. An example of this is
+\code{stat_density(geom = "area", outline.type = "both")}. The geom's
+documentation lists which parameters it can accept.
+\item Inversely, when constructing a layer using a
+\verb{geom_*()} function, the \code{...} argument can be used to pass on parameters
+to the \code{stat} part of the layer. An example of this is
+\code{geom_area(stat = "density", adjust = 0.5)}. The stat's documentation
+lists which parameters it can accept.
+\item The \code{key_glyph} argument of \code{\link[ggplot2:layer]{layer()}} may also be passed on through
+\code{...}. This can be one of the functions described as
+\link[ggplot2:draw_key]{key glyphs}, to change the display of the layer in the legend.
+}}
+}
 \description{
 \code{geom_stepribbon} is an extension of the \code{geom_ribbon}, and
 is optimized for Kaplan-Meier plots with pointwise confidence intervals
-or a confidence band.
+or a confidence band. The default \code{direction}-argument \code{"hv"} is
+appropriate for right-continuous step functions like the hazard rates etc
+returned by \code{pammtools}.
 }
 \examples{
 library(ggplot2)
diff --git a/man/get_cut.Rd b/man/get_cut.Rd
index 98df0edd..a3a2b5b2 100644
--- a/man/get_cut.Rd
+++ b/man/get_cut.Rd
@@ -3,11 +3,22 @@
 \name{get_cut}
 \alias{get_cut}
 \alias{get_cut.default}
+\alias{get_cut.list}
 \title{Obtain interval break points}
 \usage{
 get_cut(data, formula, cut = NULL, ...)
 
 \method{get_cut}{default}(data, formula, cut = NULL, max_time = NULL, event = 1L, ...)
+
+\method{get_cut}{list}(
+  data,
+  formula,
+  cut = NULL,
+  max_time = NULL,
+  event = 1L,
+  timescale = "gap",
+  ...
+)
 }
 \description{
 Default method words for data frames.
diff --git a/man/make_X.Rd b/man/make_X.Rd
new file mode 100644
index 00000000..250fb456
--- /dev/null
+++ b/man/make_X.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/add-functions.R
+\name{make_X}
+\alias{make_X}
+\title{Create model matrix from model object}
+\usage{
+make_X(object, ...)
+}
+\arguments{
+\item{object}{a model object}
+}
+\description{
+Create model matrix from model object
+}
diff --git a/man/make_X.default.Rd b/man/make_X.default.Rd
new file mode 100644
index 00000000..3a4c1a28
--- /dev/null
+++ b/man/make_X.default.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/add-functions.R
+\name{make_X.default}
+\alias{make_X.default}
+\title{Create model matrix from model object}
+\usage{
+\method{make_X}{default}(object, newdata, ...)
+}
+\arguments{
+\item{object}{a model object}
+}
+\description{
+Create model matrix from model object
+}
diff --git a/man/make_X.gam.Rd b/man/make_X.gam.Rd
new file mode 100644
index 00000000..621e75cc
--- /dev/null
+++ b/man/make_X.gam.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/add-functions.R
+\name{make_X.gam}
+\alias{make_X.gam}
+\title{Create model matrix from model object}
+\usage{
+\method{make_X}{gam}(object, newdata, ...)
+}
+\arguments{
+\item{object}{a model object}
+}
+\description{
+Create model matrix from model object
+}
+\keyword{internal}
diff --git a/man/make_X.scam.Rd b/man/make_X.scam.Rd
new file mode 100644
index 00000000..aa9354ba
--- /dev/null
+++ b/man/make_X.scam.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/add-functions.R
+\name{make_X.scam}
+\alias{make_X.scam}
+\title{Create model matrix from model object}
+\usage{
+\method{make_X}{scam}(object, newdata, ...)
+}
+\arguments{
+\item{object}{a model object}
+}
+\description{
+Create model matrix from model object
+}
+\keyword{internal}
diff --git a/man/nest_tdc.Rd b/man/nest_tdc.Rd
index e943d707..b3cebe41 100644
--- a/man/nest_tdc.Rd
+++ b/man/nest_tdc.Rd
@@ -22,10 +22,6 @@ Only TDCs present in \code{formula} will be returned.}
 part indicates the structure of the TDC structure.}
 
 \item{...}{Further arguments passed to methods.}
-
-\item{vars}{A character vector of TDCs that will be nested.}
-
-\item{id}{A character giving the name of the ID column.}
 }
 \description{
 Provides methods to nest data with time-dependent covariates (TDCs).
diff --git a/man/newdata.Rd b/man/newdata.Rd
index acfc8bc7..46f33f7c 100644
--- a/man/newdata.Rd
+++ b/man/newdata.Rd
@@ -63,9 +63,6 @@ tumor \%>\% make_newdata(days=seq_range(days, 3), status=unique(status), age=c(5
 # mean/modus values of unspecified variables are calculated over whole data
 tumor \%>\% make_newdata(sex=unique(sex))
 tumor \%>\% group_by(sex) \%>\% make_newdata()
-# You can also pass a part of the data sets as data frame to make_newdata
-purrr::cross_df(list(days = c(0, 500, 1000), sex = c("male", "female"))) \%>\%
-  make_newdata(x=tumor)
 
 # Examples for PED data
 ped <- tumor \%>\% slice(1:3) \%>\% as_ped(Surv(days, status)~., cut = c(0, 500, 1000))
diff --git a/man/pammtools.Rd b/man/pammtools.Rd
index 82527775..ff333a9e 100644
--- a/man/pammtools.Rd
+++ b/man/pammtools.Rd
@@ -44,19 +44,6 @@ Estimation for left-truncated data.
 Competing risks analysis.
 }
 }
-\references{
-Bender, Andreas, Andreas Groll, and Fabian Scheipl. 2018.
-“A Generalized Additive Model Approach to Time-to-Event Analysis”
-Statistical Modelling, February. https://doi.org/10.1177/1471082X17748083.
-
-Bender, Andreas, Fabian Scheipl, Wolfgang Hartl, Andrew G. Day, and Helmut Küchenhoff. 2019.
-“Penalized Estimation of Complex, Non-Linear Exposure-Lag-Response Associations.”
-Biostatistics 20 (2): 315–31. https://doi.org/10.1093/biostatistics/kxy003.
-
-Bender, Andreas, and Fabian Scheipl. 2018.
-“pammtools: Piece-Wise Exponential Additive Mixed Modeling Tools.”
-ArXiv:1806.01042 \link{Stat}, June. https://arxiv.org/abs/1806.01042.
-}
 \seealso{
 Useful links:
 \itemize{
@@ -80,3 +67,4 @@ Other contributors:
 }
 
 }
+\keyword{internal}
diff --git a/man/predictSurvProb.pamm.Rd b/man/predictSurvProb.pamm.Rd
index eb6715f2..a1b8158c 100644
--- a/man/predictSurvProb.pamm.Rd
+++ b/man/predictSurvProb.pamm.Rd
@@ -4,7 +4,7 @@
 \alias{predictSurvProb.pamm}
 \title{S3 method for pamm objects for compatibility with package pec}
 \usage{
-\method{predictSurvProb}{pamm}(object, newdata, times)
+\method{predictSurvProb}{pamm}(object, newdata, times, ...)
 }
 \arguments{
 \item{object}{A fitted model from which to extract predicted survival
@@ -16,6 +16,8 @@ which to compute predicted survival probabilities.}
 \item{times}{A vector of times in the range of the response variable, e.g.
 times when the response is a survival object, at which to return the
 survival probabilities.}
+
+\item{...}{Additional arguments that are passed on to the current method.}
 }
 \description{
 S3 method for pamm objects for compatibility with package pec
diff --git a/man/sim_pexp.Rd b/man/sim_pexp.Rd
index f81fe87c..1edba3c0 100644
--- a/man/sim_pexp.Rd
+++ b/man/sim_pexp.Rd
@@ -2,18 +2,9 @@
 % Please edit documentation in R/sim-pexp.R
 \name{sim_pexp}
 \alias{sim_pexp}
-\alias{sim_pexp_msm}
 \title{Simulate survival times from the piece-wise exponential distribution}
 \usage{
 sim_pexp(formula, data, cut)
-
-sim_pexp_msm(
-  t_mat,
-  data,
-  cut,
-  proportion_censoring = 0.3,
-  keep_transitions_at_risk = FALSE
-)
 }
 \arguments{
 \item{formula}{An extended formula that specifies the linear predictor.
@@ -26,19 +17,9 @@ below.}
 \item{data}{A data set with variables specified in \code{formula}.}
 
 \item{cut}{A sequence of time-points starting with 0.}
-
-\item{proportion_censoring}{The proportion of censored observations}
-
-\item{keep_transitions_at_risk}{Logical. If \code{TRUE} (default), all transitions
-for which a subject was at risk will be included in the data set.}
-
-\item{tmat}{Transition matrix containing the log-hazard definition for each transition.
-Transitions that can not occurr should contain a \code{NA}.}
 }
 \description{
 Simulate survival times from the piece-wise exponential distribution
-
-Simulate data for multi-state models
 }
 \examples{
 library(survival)
@@ -115,28 +96,4 @@ sim_df <- sim_pexp(
   data = df,
   cut = 0:10)
 }
-t_mat <- matrix(data = NA, nrow = 3, ncol = 3)
-t_mat[1,2] <- "log(0.7) + x1"
-t_mat[1,3] <- "log(0.5) - 0.25 * x2"
-t_mat[2,3] <- "log(0.9)"
-
-n = 100
-data <- cbind.data.frame(
- id = 1:n,
- x1 = runif(n, -3, 3),
- x2 = runif(n, 0, 6),
- from = 1,
- t = 0)
-
-cut =  seq(0, 3, by = 0.01)
-
-msm_df <- sim_pexp_msm(
- t_mat = t_mat,
- data = data,
- cut = cut,
- proportoin_censoring = 0.3,
- keep_transitions_at_risk = TRUE)
-head(msm_df)
-
 }
-\keyword{internal}
diff --git a/man/warn_about_new_time_points.glm.Rd b/man/warn_about_new_time_points.glm.Rd
new file mode 100644
index 00000000..a4e420a1
--- /dev/null
+++ b/man/warn_about_new_time_points.glm.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/warnings.R
+\name{warn_about_new_time_points.glm}
+\alias{warn_about_new_time_points.glm}
+\title{Warn if new t_j are used}
+\usage{
+\method{warn_about_new_time_points}{glm}(object, newdata, time_var, ...)
+}
+\description{
+Warn if new t_j are used
+}
diff --git a/tests/testthat/test-as-ped.R b/tests/testthat/test-as-ped.R
index 1a9e9569..bab22009 100644
--- a/tests/testthat/test-as-ped.R
+++ b/tests/testthat/test-as-ped.R
@@ -137,7 +137,7 @@ test_that("Trafo works for recurrent events data", {
 
 
 test_that("Trafo works for multi-state data without recurrent events", {
-  
+
   test_df <- data.frame(
     id     = c(1, 2,2),
     tstart = c(0, 0,1.2),
@@ -147,11 +147,11 @@ test_that("Trafo works for multi-state data without recurrent events", {
     to = c(3, 2,3),
     transition   = c("1->3", "1->2","2->3"),
     age    = c(24, 36,36))
-  
+
   test_df <- test_df %>% filter(status == 1) %>% add_counterfactual_transitions()
-  
+
   print(test_df)
-  
+
   # CALENDAR timescale
   cal_df <- as_ped_multistate(
     data       = test_df,
@@ -159,9 +159,9 @@ test_that("Trafo works for multi-state data without recurrent events", {
     transition = "transition",
     id         = "id",
     timescale  = "calendar")
-  
+
   print(cal_df)
-  
+
   # according to code: order by transition -> id -> tstart
   expect_data_frame(cal_df, nrows = 7L, ncols = 10L)
   expect_identical(cal_df$transition,
@@ -187,6 +187,6 @@ test_that("Trafo works for multi-state data without recurrent events", {
     cal_df$to,
     c(2,2,2, 3,3,3, 3)
   )
-  
+
   }
 )
diff --git a/vignettes/baseline.Rmd b/vignettes/baseline.Rmd
index 564ef93c..68073c01 100644
--- a/vignettes/baseline.Rmd
+++ b/vignettes/baseline.Rmd
@@ -3,7 +3,7 @@ title: "Baseline Hazard"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
- %\VignetteIndexEntry{baseline}
+ %\VignetteIndexEntry{Baseline Hazard}
  %\VignetteEngine{knitr::rmarkdown}
 ---
 
diff --git a/vignettes/basics.Rmd b/vignettes/basics.Rmd
index fff1c470..f5a6ee6f 100644
--- a/vignettes/basics.Rmd
+++ b/vignettes/basics.Rmd
@@ -56,7 +56,7 @@ library(mgcv)
 library(pammtools)
 library(dplyr)
 
-data("veteran", package = "survival") # load veteran data
+veteran <- survival::veteran
 # remove ties to illustrate equivalence with Cox approach
 vetu <- filter(veteran, !duplicated(time))
 ped_vetu <- vetu %>%
diff --git a/vignettes/competing-risks.Rmd b/vignettes/competing-risks.Rmd
index e47c6c5d..1dc7a36e 100644
--- a/vignettes/competing-risks.Rmd
+++ b/vignettes/competing-risks.Rmd
@@ -4,7 +4,7 @@ author: "Andreas Bender"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{competing-risks}
+  %\VignetteIndexEntry{Competing Risks}
   %\VignetteEngine{knitr::rmarkdown}
 ---
 
@@ -156,79 +156,3 @@ ggplot(ndf, aes(x = tend, y = cif)) +
     alpha = .3) +
   facet_wrap(~cause, labeller = label_both)
 ```
-
-
-# `mgus2` Data
-
-```{r}
-library(survival)
-etime <- with(mgus2, ifelse(pstat == 0, futime, ptime))
-event <- with(mgus2, ifelse(pstat == 0, 2 * death, 1))
-event <- factor(event, 0:2, labels=c("censor", "pcm", "death"))
-mfit2 <- survfit(Surv(etime, event) ~ sex, data=mgus2)
-plot(mfit2, col=c(1,2,1,2), lty=c(2,2,1,1), mark.time=FALSE, lwd=2, xscale=12,
-  xlab="Years post diagnosis", ylab="Probability in State")
-legend(240, .6, c("death:female", "death:male", "pcm:female", "pcm:male"),
-  col=c(1,2,1,2), lty=c(1,1,2,2), lwd=2, bty='n')
-
-dat <- mgus2
-dat$event <- as.numeric(event) - 1
-dat$time <- etime
-dat$sex <- as.factor(sex)
-ped <- dat %>%
-  as_ped(Surv(time, event)~. - ptime - pstat - futime - death) %>%
-  mutate(cause = as.factor(cause))
-pam <- pamm(ped_status ~ cause + s(tend, by = cause) + sex, data = ped, family = poisson(), offset = offset, discrete = TRUE, method = "fREML", engine = "bam")
-summary(pam)
-ndf <- ped %>%
-  make_newdata(tend = unique(tend), sex = unique(sex), cause = unique(cause)) %>%
-  group_by(cause, sex) %>%
-  add_cif(pam)
-df_cox <- as.data.frame(mfit2$pstate)
-colnames(df_cox) <- c("0", "1", "2")
-df_cox$sex <- as.factor(rep(c("F", "M"), each = nrow(df_cox)/2))
-df_cox$tend <- mfit2$time
-df_cox <- df_cox %>%
-  select(-0) %>%
-  tidyr::pivot_longer(cols = c("1", "2"), names_to="cause", values_to = "cif") %>%
-  mutate(cause = as.factor(cause))
-ggplot(ndf, aes(x = tend, y = cif, col = sex)) +
-  geom_line() +
-  # geom_ribbon(
-  #   aes(ymin = cif_lower, ymax = cif_upper, fill = sex),
-  #   alpha = .3) +
-  facet_wrap(~cause, labeller = label_both) +
-  geom_step(data = df_cox, lty = 2, ) +
-  scale_y_continuous(breaks = seq(0, .8, by = .2))
-
-# multi-state formulation
-ptemp <- with(mgus2, ifelse(ptime==futime & pstat==1, ptime-.1, ptime))
-data3 <- tmerge(mgus2, mgus2, id=id, death=event(futime, death),pcm = event(ptemp, pstat))
-data3 <- tmerge(data3, data3, id, enum=cumtdc(tstart))
-temp <- with(data3, ifelse(death==1, 2, pcm))
-data3$event <- temp
-data3$event2 <- factor(temp, 0:2, labels=c("censor", "pcm", "death"))
-mfit3 <- survfit(Surv(tstart, tstop, event2) ~ sex, data=data3, id=id)
-
-datt <- data3 %>%
-  mutate(from = ifelse(enum == 1, 0, 1)) %>%
-  mutate(
-    to = event,
-    status = 1L * (event != 0),
-    to = ifelse(to == 0, 2, to),
-    transition = paste0(from, "->", to)
-  )
-table(datt$transition)
-
-dat_ped <- add_counterfactual_transitions(datt)
-ped <- as_ped(dat_ped, Surv(tstart, tstop, status)~., transition = "transition",
-  timescale = "calendar")
-pam <- pamm(ped_status ~ transition + s(tend, by = transition) + sex,
-  data = ped, engine = "bam", method = "fREML", discrete = TRUE)
-summary(pam)
-
-ndf <- ped %>%
-  make_newdata(tend = unique(tend), transition = unique(transition), sex = unique(sex)) %>%
-  group_by(transition, sex) %>%
-  add_cif(pam)
-```
diff --git a/vignettes/convenience.Rmd b/vignettes/convenience.Rmd
index a003a671..60a57449 100644
--- a/vignettes/convenience.Rmd
+++ b/vignettes/convenience.Rmd
@@ -3,7 +3,7 @@ title: "Workflow and convenience functions"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{convenience}
+  %\VignetteIndexEntry{Workflow and convenience functions}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 link-citations: yes
diff --git a/vignettes/cumulative-effects.Rmd b/vignettes/cumulative-effects.Rmd
index 565408f8..b3a113de 100644
--- a/vignettes/cumulative-effects.Rmd
+++ b/vignettes/cumulative-effects.Rmd
@@ -3,7 +3,7 @@ title: "Cumulative Effects"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{cumulative-effects}
+  %\VignetteIndexEntry{Cumulative Effects}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 link-citations: yes
diff --git a/vignettes/data-transformation.Rmd b/vignettes/data-transformation.Rmd
index 8aa9e0c7..3e5c8c37 100644
--- a/vignettes/data-transformation.Rmd
+++ b/vignettes/data-transformation.Rmd
@@ -3,7 +3,7 @@ title: "Data Transformation"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{data-transformation}
+  %\VignetteIndexEntry{Data Transformation}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 ---
diff --git a/vignettes/frailty.Rmd b/vignettes/frailty.Rmd
index 34ef07e1..3ed944ef 100644
--- a/vignettes/frailty.Rmd
+++ b/vignettes/frailty.Rmd
@@ -3,7 +3,7 @@ title: "Frailty and random effects"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{frailty}
+  %\VignetteIndexEntry{Frailty and random effects}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 ---
diff --git a/vignettes/left-truncation.Rmd b/vignettes/left-truncation.Rmd
index 0cba2fa8..78daac1a 100644
--- a/vignettes/left-truncation.Rmd
+++ b/vignettes/left-truncation.Rmd
@@ -3,7 +3,7 @@ title: "Left Truncation"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{left-truncation}
+  %\VignetteIndexEntry{Left Truncation}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 link-citations: yes
diff --git a/vignettes/model-evaluation.Rmd b/vignettes/model-evaluation.Rmd
index cc0bd11a..49dfccd5 100644
--- a/vignettes/model-evaluation.Rmd
+++ b/vignettes/model-evaluation.Rmd
@@ -3,7 +3,7 @@ title: "Model evaluation"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{model-evaluation}
+  %\VignetteIndexEntry{Model evaluation}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 link-citations: yes
diff --git a/vignettes/recurrent-events.Rmd b/vignettes/recurrent-events.Rmd
index 53bbc5e1..90822018 100644
--- a/vignettes/recurrent-events.Rmd
+++ b/vignettes/recurrent-events.Rmd
@@ -5,7 +5,7 @@ date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 number_sections: true
 vignette: >
-  %\VignetteIndexEntry{recurrent-events}
+  %\VignetteIndexEntry{Recurrent Events}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 link-citations: yes
diff --git a/vignettes/smj-tutorial.Rmd b/vignettes/smj-tutorial.Rmd
index 5f8c3883..e2572018 100644
--- a/vignettes/smj-tutorial.Rmd
+++ b/vignettes/smj-tutorial.Rmd
@@ -3,7 +3,7 @@ title: "Paper: A generalized additive model approach to time-to-event analysis"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{smj-tutorial}
+  %\VignetteIndexEntry{Paper: A generalized additive model approach to time-to-event analysis}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---
diff --git a/vignettes/splines.Rmd b/vignettes/splines.Rmd
index 1c23fe44..3374be7c 100644
--- a/vignettes/splines.Rmd
+++ b/vignettes/splines.Rmd
@@ -3,7 +3,7 @@ title: "Non-linear effects (penalized splines)"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{splines.html}
+  %\VignetteIndexEntry{Non-linear effects (penalized splines)}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 ---
@@ -35,7 +35,7 @@ We first illustrate the estimation on non-linear smooth effects using
 the `veteran` data from the survival package:
 
 ```{r}
-data("veteran", package="survival") # load veteran package
+veteran <- survival::veteran # load veteran package
 veteran <- veteran %>%
   mutate(
     trt   = 1+(trt == 2),
diff --git a/vignettes/strata.Rmd b/vignettes/strata.Rmd
index 28018362..659d786d 100644
--- a/vignettes/strata.Rmd
+++ b/vignettes/strata.Rmd
@@ -3,7 +3,7 @@ title: "Stratified baseline"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{strata}
+  %\VignetteIndexEntry{Stratified baseline}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 ---
@@ -43,7 +43,7 @@ from the `survival` package,
 where the cell type of the lung cancer had four levels:
 
 ```{r}
-data("veteran") # load veteran data
+veteran <- survival::veteran # load veteran data
 veteran <- veteran %>%
   mutate(
     trt   = 1*(trt == 2),
diff --git a/vignettes/tdcovar.Rmd b/vignettes/tdcovar.Rmd
index ba306a62..7b3d35c2 100644
--- a/vignettes/tdcovar.Rmd
+++ b/vignettes/tdcovar.Rmd
@@ -3,7 +3,7 @@ title: "Time-dependent covariates"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{tdcovar}
+  %\VignetteIndexEntry{Time-dependent covariates}
   %\VignetteEngine{knitr::rmarkdown}
 bibliography: Remote.bib
 ---
diff --git a/vignettes/tveffects.Rmd b/vignettes/tveffects.Rmd
index b08d2cf5..6921fb9b 100644
--- a/vignettes/tveffects.Rmd
+++ b/vignettes/tveffects.Rmd
@@ -3,7 +3,7 @@ title: "Time-varying effects"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{tveffects}
+  %\VignetteIndexEntry{Time-varying effects}
   %\VignetteEngine{knitr::rmarkdown}
 ---
 
@@ -71,7 +71,7 @@ Karnofsky performance scores `karno` (the higher the better), `age` and whether
 
 ```{r}
 # for some reason the prior variable is coded 0/10 instead of 0/1
-data("veteran", package = "survival")
+veteran <- survival::veteran
 veteran <- veteran %>%
   mutate(
     trt   = 1L * (trt == 2),