diff --git a/.Rbuildignore b/.Rbuildignore
index 823e4f0..8247b0a 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -11,3 +11,4 @@ dev
 ^CONTRIBUTING.md$
 ^DOCKERFILE$
 ^README.Rmd$
+^tests/script
diff --git a/DESCRIPTION b/DESCRIPTION
index bbdf720..8fc6d1c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -3,13 +3,13 @@ Title: Perform productivity analyses in a web browser
 Version: 0.3.0.9000
 Authors@R: c(
     person("Ove Haugland", "Jakobsen", role = c("aut", "cre"), email = "ohj@mac.com"),
+    person("Aleksander", "Eilertsen", role = "aut", email = "ale@riksrevisjonen.no"),
     person("Jan Roar", "Beckstrøm", role = "ctb", email = "jrb@riksrevisjonen.no"),
     person("Lars Skaage", "Engebretsen", role = "ctb", email = "lse@riksrevisjonen.no"),
     person("Jonas", "Månsson", role = "ctb", email = "jonas.mansson@riksrevisjonen.no"),
-    person("Aleksander", "Eilertsen", role = "ctb", email = "ale@riksrevisjonen.no"),
     person("Joachim", "Sandnes", role = "ctb", email = "jsn@riksrevisjonen.no"),
     person(family = "Riksrevisjonen", role = "cph"))
-Description: What the package does (one paragraph).
+Description: Application to conduct productivity and efficiency (DEA) analysis. 
 License: GPL-3 | file LICENSE
 URL: https://riksrevisjonen.github.io/pioneeR/, https://github.com/Riksrevisjonen/pioneeR
 BugReports: https://github.com/Riksrevisjonen/pioneeR/issues
@@ -17,7 +17,7 @@ Depends:
     R (>= 4.1.0),
     shiny (>= 1.7.0)
 Imports:
-    Benchmarking (>= 0.30),
+    lpSolveAPI (>= 5.5.2),
     haven (>= 2.5.0),
     markdown,
     productivity (>= 1.1.0),
@@ -25,22 +25,22 @@ Imports:
     writexl,
     bslib (>= 0.5.1),
     ggplot2 (>= 3.3.0),
-    reactable (>= 0.4.0)
+    reactable (>= 0.4.0),
+    cli (>= 3.6.0)
 Suggests:
-    cli,
+    Benchmarking (>= 0.30),
     data.table,
-    lpSolve,
-    lpSolveAPI,
     rmarkdown,
     knitr,
     scales,
     rlang,
     testthat (>= 3.0.0),
-    tibble
+    tibble,
+    withr
 Language: en-GB
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 VignetteBuilder: knitr
 Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index e702b5c..a0a84f4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,6 +2,8 @@
 
 S3method(summary_tbl_dea,Farrell)
 S3method(summary_tbl_dea,numeric)
+S3method(summary_tbl_dea,pioneer_dea)
+export(compute_dea)
 export(compute_scale_efficiency)
 export(create_matrix)
 export(runPioneeR)
diff --git a/NEWS.md b/NEWS.md
index 33f40ca..d26c1d2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,11 +1,18 @@
 # pioneeR 0.3.0.9000
 
-- Fixed bug where checkbox for timeseries would be available before any data is uploaded (#74)
+## Changes & improvements
+
 - Updated documentation for `run_pioneer()`
+- Added custom functions for DEA calculations in order to remove `{Benchmarking}` as a dependency (#80)
+
+## Bug fixes
+
+- Fixed bug where checkbox for timeseries would be available before any data is uploaded (#74)
+- Fixed a few bugs where the UI did not correctly reflect the user supplied input for rounding (#84)
 
 # pioneeR 0.3.0
 
-## Breaking changes:
+## Breaking changes
 
 - pioneeR now depends on bslib 0.5.1 or higher. Version 0.6.0 or higher is recommended, and support for version 0.5.1 will be deprecated in the next release
 
diff --git a/R/compute_dea.R b/R/compute_dea.R
new file mode 100644
index 0000000..26814ef
--- /dev/null
+++ b/R/compute_dea.R
@@ -0,0 +1,91 @@
+#' Compute DEA
+#'
+#' DEA analysis.
+#'
+#' @param data Dataset to analyse.
+#' @param id A string with the DMU id or name variable.
+#' @param input A character vector with input variables.
+#' @param output A character vector with output variables.
+#' @param rts Returns to scale.
+#' \tabular{rlll}{
+#'    \tab crs \tab\tab Constant returns to scale, convexity and free disposability. \cr
+#'    \tab vrs \tab\tab Variable returns to scale, convexity and free disposability. \cr
+#'    \tab drs \tab\tab Decreasing returns to scale, convexity, down-scaling and free disposability. \cr
+#'    \tab irs \tab\tab Increasing returns to scale, (up-scaling, but not down-scaling), convexity and free disposability.
+#' }
+#' @param orientation Model orientation.
+#' @param super If `TRUE` super efficiency scores are calculated.
+#' @param slack If `TRUE` slack values are calculated.
+#' @param peers If `TRUE` peers are added to the response.
+#' @return list
+#' @export
+compute_dea <- function(
+    data,
+    id,
+    input,
+    output,
+    rts = c('crs', 'vrs', 'drs', 'irs'),
+    orientation = c('in', 'out'),
+    super = FALSE,
+    slack = FALSE,
+    peers = FALSE) {
+
+  rts <- match.arg(rts)
+  orientation <- match.arg(orientation)
+
+  # Prepare and validate input
+  x <- as.matrix(data[input])
+  y <- as.matrix(data[output])
+  ids <- data[[id]]
+  check_data(x, y)
+
+  # Calculate DEA metrics
+  res <- compute_dea_(
+    x, y, ids,
+    rts = rts,
+    orientation = orientation,
+    super = super,
+    slack = slack,
+    peers = peers
+  )
+
+  res
+
+}
+
+#' compute_dea_ (internal helper)
+#' @inheritParams compute_dea
+#' @inheritParams compute_efficiency
+#' @return list
+#' @noRd
+compute_dea_ <- function(x, y, ids, rts, orientation, super, slack, peers) {
+
+  # Set initial values
+  super_res <- NULL
+  slack_res <- NULL
+  peers_res <- NULL
+
+  eff_res <- compute_efficiency(x, y, rts = rts, orientation = orientation)
+  if (super) super_res <- compute_super_efficiency(x, y, rts = rts, orientation = orientation)
+  if (slack) slack_res <- compute_slack(x, y, eff_res)
+  if (peers) {
+    peers_res <- get_peers(eff_res, ids)
+    peers_res <- peers_res[2:ncol(peers_res)]
+  }
+
+  res <- list(
+    efficiency = eff_res$values,
+    super_efficiency = super_res$values,
+    slack = slack_res$values
+  )
+  res <- do.call('cbind', res)
+  if (peers) res <- cbind(res, peers_res)
+
+  out <- structure(
+    list(
+      results = data.frame(dmu = ids, res),
+      info = eff_res$info
+    ), class = 'pioneer_dea'
+  )
+
+}
diff --git a/R/dea-methods.R b/R/dea-methods.R
new file mode 100644
index 0000000..90fd5c3
--- /dev/null
+++ b/R/dea-methods.R
@@ -0,0 +1,182 @@
+#' Compute efficiency
+#'
+#' Compute DEA efficiency scores.
+#'
+#' @inheritParams compute_dea
+#' @param x A matrix with input values.
+#' @param y A matrix with output values.
+#' @param xref A matrix with input values for the reference technology. Defaults to `NULL`.
+#' @param yref A matrix with output values for the reference technology. Defaults to `NULL`.
+#' @param values_only Only return the calculated values.
+#' @return list
+#' @noRd
+compute_efficiency <- function(
+    x, y,
+    xref = NULL,
+    yref = NULL,
+    rts,
+    orientation,
+    values_only = FALSE) {
+  # Get dims
+  dims <- get_dims(x, y, xref, yref, rts = rts, slack = FALSE)
+  # Create linear program model
+  lp <- create_lp(
+    x = x, y = y,
+    xref = xref, yref = yref,
+    rts = rts,
+    orientation = orientation,
+    m = dims$n_inputs,
+    n = dims$n_outputs,
+    n_units = dims$n_units,
+    n_constraints = dims$n_constraints,
+    n_vars = dims$n_vars,
+    slack = FALSE)
+  # Solve model
+  res <- solve_lp(
+    lp, x, y,
+    orientation = orientation,
+    m = dims$n_inputs,
+    n = dims$n_outputs)
+  # Adjust values
+  res <- adjust_values(res$values, res$lambda, res$eps)
+  # Return
+  res <- create_dea_output(
+    res,
+    rts = rts,
+    orientation = orientation,
+    dims = dims,
+    values_only = values_only)
+  res
+}
+
+#' Compute super efficiency
+#'
+#' Compute DEA super efficiency scores.
+#'
+#' @inheritParams compute_efficiency
+#' @noRd
+compute_super_efficiency <- function(
+    x, y,
+    rts,
+    orientation) {
+  dims <- get_dims(x, y, rts = rts, slack = FALSE)
+  lambda <- matrix(0, dims$n_units, dims$n_units)
+  values <- rep(NA_real_, dims$n_units)
+  for (i in seq_len(dims$n_units)) {
+    # Create linear program model
+    lp <- create_lp(
+      x[-i,,drop=FALSE],
+      y[-i,,drop=FALSE],
+      rts = rts,
+      orientation = orientation,
+      m = dims$n_inputs,
+      n = dims$n_outputs,
+      n_units = dims$n_units - 1,
+      n_constraints = dims$n_constraints,
+      n_vars = dims$n_vars - 1,
+      slack = FALSE)
+    # Solve model
+    tmp <- solve_lp_single(
+      x[i,,drop=FALSE],
+      y[i,,drop=FALSE],
+      lp = lp,
+      orientation = orientation,
+      m = dims$n_inputs,
+      n = dims$n_outputs)[[1]]
+    values[i] <- tmp$value
+    lambda[i,-i] <- tmp$solution
+  }
+  eps <- lpSolveAPI::lp.control(lp)$epsilon['epsint']
+  res <- adjust_values(values, lambda, eps)
+  # Return
+  res <- create_dea_output(
+    res,
+    rts = rts,
+    orientation = orientation,
+    dims = dims,
+    values_only = FALSE)
+  res
+}
+
+#' Compute slack
+#'
+#' Compute slack for a DEA efficiency analysis.
+#'
+#' @inheritParams compute_efficiency
+#' @param model Output of `compute_efficiency()`.
+#' @return list
+#' @noRd
+compute_slack <- function(x, y, model) {
+  # Get dims
+  dims <- get_dims(x, y, rts = model$info$rts, slack = TRUE)
+  # Scale inputs and outputs if needed
+  scale <- scale_vars(
+    x, y,
+    m = dims$n_inputs,
+    n = dims$n_outputs,
+    n_units = dims$n_units)
+  # Create linear program model
+  lp <- create_lp(
+    scale$x, scale$y,
+    rts = model$info$rts,
+    orientation = model$info$orientation,
+    m = dims$n_inputs,
+    n = dims$n_outputs,
+    n_units = dims$n_units,
+    n_constraints = dims$n_constraints,
+    n_vars = dims$n_vars,
+    n_lambda = dims$n_lambda,
+    slack = TRUE)
+  # Solve model
+  res <- solve_lp_slack(
+    lp,
+    scale$x, scale$y,
+    values = model$unadj_values, # slack is calculated based on unadjusted eff.scores
+    orientation = model$info$orientation,
+    m = dims$n_inputs,
+    n = dims$n_outputs,
+    n_units = dims$n_units,
+    n_vars = dims$n_vars)
+  # Adjust values
+  res <- adjust_slack(
+    lp,
+    values = res$values,
+    sx = res$sx,
+    sy = res$sy,
+    lambda = res$lambda,
+    scaling = scale,
+    m = dims$n_inputs,
+    n = dims$n_outputs,
+    n_units = dims$n_units)
+  # Return
+  res <- create_dea_output(
+    res,
+    rts = model$info$rts,
+    orientation = model$info$orientation,
+    dims = dims,
+    values_only = FALSE)
+  res
+}
+
+#' Get peers
+#'
+#' Get peers for each DMU.
+#'
+#' @param model Output of `compute_efficiency()`.
+#' @param ids A vector with DMU names or ids.
+#' @param threshold Minimum weight for extracted peers. Defaults to 0.
+#' @return data.frame
+#' @noRd
+get_peers <- function(model, ids, threshold = 0) {
+  lambda <- model$lambda
+  pt_ <- apply(lambda, 1, function(x) {which(x > threshold)})
+  if (dim(lambda)[1] == 1) pt_ <- list(c(pt_))  # Only 1 DMU
+  bench <- t(mapply(function(x) x[1:max(sapply(pt_, length))], pt_))
+  maxp <- max(sapply(pt_, length))
+  if (maxp == 0 | is.na(maxp)) maxp <- 1
+  dim(bench) <- c(dim(lambda)[1], maxp)
+  bench <- matrix(ids[bench], nrow=nrow(bench))
+  colnames(bench) <- paste0('peer', 1:ncol(bench))
+  bench <- data.frame(cbind(dmu = ids, bench))
+  bench
+}
diff --git a/R/dea-utils.R b/R/dea-utils.R
index 9dd30a9..36142fa 100644
--- a/R/dea-utils.R
+++ b/R/dea-utils.R
@@ -42,8 +42,7 @@ compute_scale_efficiency <- function(
     x,
     y,
     orientation = c('in', 'out'),
-    digits = NULL)
-{
+    digits = NULL) {
 
   if (!is.matrix(x)) {
     cli::cli_abort('inputs must be a matrix of input values')
@@ -71,26 +70,26 @@ compute_scale_efficiency <- function(
   orientation <- match.arg(orientation)
 
   # Run DEA models
-  crs_mod <- Benchmarking::dea(x, y, RTS = 'crs', ORIENTATION = orientation)
-  vrs_mod <- Benchmarking::dea(x, y, RTS = 'vrs', ORIENTATION = orientation)
-  nirs_mod <- Benchmarking::dea(x, y, RTS = 'drs', ORIENTATION = orientation)
+  crs_mod <- compute_efficiency(x, y, rts = 'crs', orientation = orientation, values_only = TRUE)$values
+  vrs_mod <- compute_efficiency(x, y, rts = 'vrs', orientation = orientation, values_only = TRUE)$values
+  nirs_mod <- compute_efficiency(x, y, rts = 'drs', orientation = orientation, values_only = TRUE)$values
 
   # If efficiency scores for a unit differs in the CRS and VRS models and the ratio
   # of the NIRS and VRS models is equal to 1, the unit should decrease its size. When
   # the CRS and VRS models differ, and the ratio of the NIRS and VRS models is *not*
   # equal to 1, the unit should increase its size.
-  equal_crs_vrs <- round(crs_mod$eff, 6L) == round(vrs_mod$eff, 6L)
+  equal_crs_vrs <- round(crs_mod, 6L) == round(vrs_mod, 6L)
   optimal_scale <- ifelse(
     !equal_crs_vrs,
-    ifelse(vrs_mod$eff / nirs_mod$eff == 1, 'Decrease', 'Increase'),
+    ifelse(vrs_mod / nirs_mod == 1, 'Decrease', 'Increase'),
     '-'
   )
 
   out_mod <- data.frame(
-    'CRS' = crs_mod$eff,
-    'VRS' = vrs_mod$eff,
-    'Scale.eff.' = crs_mod$eff / vrs_mod$eff,
-    'VRS.NIRS.ratio' = vrs_mod$eff / nirs_mod$eff
+    'CRS' = crs_mod,
+    'VRS' = vrs_mod,
+    'Scale.eff.' = crs_mod / vrs_mod,
+    'VRS.NIRS.ratio' = vrs_mod / nirs_mod
   )
 
   if (!is.null(digits) || is.numeric(digits)) {
@@ -107,3 +106,83 @@ compute_scale_efficiency <- function(
   out_mod
 
 }
+
+#' check_data
+#' @noRd
+check_data <- function(x, y, xref = NULL, yref = NULL) {
+  check_numeric(x, y)
+  check_nunits(x, y)
+  if (!is.null(xref) || !is.null(yref)) {
+    if (ncol(xref) != ncol(x)) {
+      msg <- c('Inconsistent number of input variables:' ,
+               'x' = "You've supplied are {ncol(x)} column(s) for `x` and {ncol(xref)} column(s) for `xref`.")
+      cli::cli_abort(msg)
+    }
+    if (ncol(yref) != ncol(y)) {
+      msg <- c('Inconsistent number of output variables:' ,
+               'x' = "You've supplied are {ncol(y)} column(s) for `y` and {ncol(yref)} column(s) for `yref`.")
+      cli::cli_abort(msg)
+    }
+    check_numeric(xref, yref)
+    check_nunits(xref, yref)
+  }
+}
+
+#' check_data
+#' @noRd
+check_numeric <- function(x, y) {
+  chk <- apply(cbind(x, y), 2, is.numeric)
+  if (any(!chk)) {
+    cli::cli_abort('All variables must be numeric.')
+  }
+}
+
+#' check_data
+#' @noRd
+check_nunits <- function(x, y, ref = FALSE) {
+  nx <- nrow(x)
+  ny <- nrow(y)
+  if (nx != ny) {
+    if (ref) {
+      msg <- "You've supplied {nrows(x)} rows for `xref` and {nrow(y)} column(s) for `yref`."
+    } else {
+      msg <- "You've supplied {nrows(x)} rows for `x` and {nrow(y)} column(s) for `y`."
+    }
+    cli::cli_abort(c('Inconsistent number of units: ', 'x' = msg))
+  }
+}
+
+#' create_dea_output
+#' @noRd
+create_dea_output <- function(res, rts, orientation, dims, values_only) {
+  if (values_only) {
+    res <- list(values = res$values)
+  } else {
+    res$info <- create_model_info(rts, orientation, dims)
+  }
+  res
+}
+
+#' create_model_info
+#' @noRd
+create_model_info <- function(rts, orientation, dims) {
+  list(
+    rts = rts,
+    orientation = orientation,
+    dims = dims
+  )
+}
+
+#' round_numeric
+#' @noRd
+round_numeric <- function(df, digits = 4L) {
+  df[] <- lapply(df, function(x) if(is.numeric(x)) round(x, digits) else x)
+  df
+}
+
+#' pioneer_dea print method
+#' @noRd
+print.pioneer_dea <- function(x, ...) {
+  cat("DEA result:\n")
+  utils::str(x)
+}
diff --git a/R/lp-utils.R b/R/lp-utils.R
new file mode 100644
index 0000000..366fdd3
--- /dev/null
+++ b/R/lp-utils.R
@@ -0,0 +1,330 @@
+#' Create linear program
+#' @inheritParams compute_efficiency
+#' @param m Number of inputs.
+#' @param n Number of outputs.
+#' @param n_units Number of units (DMUs).
+#' @param n_constraints Number of constraints.
+#' @param n_vars Number of decision variables.
+#' @param n_lambda Number of lambda variables.
+#' @param slack If `TRUE` creates a model for slack analysis.
+#' @noRd
+create_lp <- function(x, y,
+                      xref = NULL,
+                      yref = NULL,
+                      rts,
+                      orientation,
+                      m, n,
+                      n_units,
+                      n_constraints,
+                      n_vars,
+                      n_lambda = NULL,
+                      slack = FALSE) {
+  if (is.null(xref) && is.null(yref)) {
+    lp <- create_lp_(
+      x, y,
+      rts = rts,
+      orientation = orientation,
+      m = m,
+      n = n,
+      n_units = n_units,
+      n_constraints = n_constraints,
+      n_vars = n_vars,
+      n_lambda = n_lambda,
+      slack = slack)
+  } else {
+    lp <- create_lp_(
+      xref, yref,
+      rts = rts,
+      orientation = orientation,
+      m = m,
+      n = n,
+      n_units = n_units,
+      n_constraints = n_constraints,
+      n_vars = n_vars,
+      n_lambda = n_lambda,
+      slack = slack)
+  }
+  lp
+}
+
+#' Create linear program
+#' @inheritParams create_lp
+#' @noRd
+create_lp_ <- function(x, y, rts, orientation, m, n, n_units, n_constraints,
+                       n_vars, n_lambda, slack) {
+  lp <- lpSolveAPI::make.lp(n_constraints, n_vars)
+  if (slack) {
+    lpSolveAPI::set.objfn(lp, rep(1, m + n), 1:(m + n))
+    lpSolveAPI::set.constr.type(lp, c(rep('=', m + n), rep('<=', n_lambda)))
+    lpSolveAPI::lp.control(lp, sense='max')
+  } else {
+    lpSolveAPI::set.objfn(lp, 1, 1)
+    lpSolveAPI::set.constr.type(lp, rep('>=', n_constraints))
+    if (orientation == 'in') lpSolveAPI::lp.control(lp, sense = 'min')
+    if (orientation == 'out') lpSolveAPI::lp.control(lp, sense = 'max')
+  }
+  lpSolveAPI::lp.control(lp, scaling = c('range', 'equilibrate', 'integers'))
+  lp <- set_lp_constraints(
+    lp, x, y, rts, m = m, n = n, n_units = n_units,
+    n_vars = n_vars, slack = slack)
+  lp
+}
+
+#' set_lp_constraints
+#' @param lp An lpSolve linear program model object.
+#' @inheritParams create_lp
+#' @noRd
+set_lp_constraints <- function(lp, x, y, rts, m, n, n_units, n_vars, slack) {
+  if (slack) {
+    l1 <- m+n+1
+    # Set restrictions on matrix
+    for (h in seq_len(m)) lpSolveAPI::set.row(lp, h, x[,h], l1:n_vars) # (m+n+1):(m+n+n_units)
+    for (h in seq_len(n)) lpSolveAPI::set.row(lp, m+h, -y[,h], l1:n_vars) # (m+n+1):(m+n+n_units)
+    for (h in seq_len(m+n)) lpSolveAPI::set.mat(lp, h, h, 1)
+  } else {
+    # Set restrictions on input/output
+    for (h in seq_len(m)) lpSolveAPI::set.row(lp, h, c(0, -x[,h])) # inputs
+    for (h in seq_len(n)) lpSolveAPI::set.row(lp, m+h, c(0, y[,h])) # outputs
+  }
+  lp <- set_lambda_constraints(
+    lp, rts = rts, m = m, n = n,
+    n_units = n_units, slack = slack)
+  lp
+}
+
+#' set_lambda_constraints
+#' @noRd
+set_lambda_constraints <- function(lp, rts, m, n, n_units, slack){
+  if (slack) q <- rep(0, m+n) else q <- 0 # Set number of starting zeros
+  if (rts != 'crs') {
+    l1 <- m+n+1; l2 <- m+n+2
+    if (rts == 'vrs') {
+      lpSolveAPI::set.row(lp, l1, c(q, rep(-1, n_units)))
+      lpSolveAPI::set.row(lp, l2, c(q, rep(1, n_units)))
+      lpSolveAPI::set.rhs(lp, c(-1, 1), l1:l2)
+    } else if (rts == 'irs') {
+      lpSolveAPI::set.row(lp, l1, c(q, rep(1, n_units)))
+      lpSolveAPI::set.rhs(lp, 1, l1)
+      if (slack) lpSolveAPI::set.constr.type(lp, '>=', l1)
+    } else if (rts == 'drs') {
+      lpSolveAPI::set.row(lp, l1, c(q, rep(-1, n_units)))
+      lpSolveAPI::set.rhs(lp, -1, l1)
+      if (slack) lpSolveAPI::set.constr.type(lp, '>=', l1)
+    }
+  }
+  lp
+}
+
+#' solve_lp
+#' @inheritParams create_lp
+#' @noRd
+solve_lp <- function(lp, x, y, orientation, m, n) {
+  # Solve LP for each row
+  dl <- mapply(\(x, y) {
+    solve_lp_single(x, y, lp, orientation, m = m, n = n)
+  }, x = asplit(x, 1), y = asplit(y, 1))
+  values <- vapply(dl, \(x) x$value, FUN.VALUE = numeric(1))
+  lambda <- do.call('rbind', lapply(dl, \(x) x$solution)) |> as.matrix()
+  eps <- lpSolveAPI::lp.control(lp)$epsilon['epsint']
+  list(
+    values = values,
+    lambda = lambda,
+    eps = eps
+  )
+}
+
+#' solve_lp_single
+#' @noRd
+solve_lp_single <- function(x, y, lp, orientation, m, n){
+  if (orientation == 'in') {
+    lpSolveAPI::set.column(lp, 1, c(1, x), 0:m)
+    lpSolveAPI::set.rhs(lp, y, (m + 1):(m + n))
+  } else {
+    lpSolveAPI::set.column(lp, 1, c(1, -y), c(0, (m + 1):(m + n)))
+    lpSolveAPI::set.rhs(lp, -x, 1:m)
+  }
+  lpSolveAPI::set.basis(lp, default = TRUE)
+  status <- lpSolveAPI::solve.lpExtPtr(lp)
+  value <- get_lp_objective(lp, status, orientation)
+  solution <- get_lp_solution(lp, status)[-1]
+  list(list(value = value, solution = solution))
+}
+
+#' solve_lp_slack
+#' @inheritParams create_lp
+#' @noRd
+solve_lp_slack <- function(lp, x, y, values, orientation, m, n, n_units, n_vars) {
+  if (orientation == 'in') {
+    e_values <- values
+    f_values <- rep(1, n_units)
+  } else if (orientation == 'out') {
+    f_values <- values
+    e_values <- rep(1, n_units)
+  }
+  rhs <- cbind(e_values * x, -f_values * y)
+  mn <- m+n
+  mn1 <- mn+1
+  values <- rep(NA_real_, n_units)
+  sx <- matrix(NA_real_, n_units, m)
+  sy <- matrix(NA_real_, n_units, n)
+  lambda <- matrix(NA_real_, n_units, n_units)
+  for (i in seq_len(n_units)) {
+    lpSolveAPI::set.rhs(lp, rhs[i,], 1:mn)
+    lpSolveAPI::set.basis(lp, default=TRUE)
+    status <- solve(lp)
+    if (status != 0) {
+      if (status == 2 || status == 3) {
+        values[i] <- 0
+        solution <- rep(0, n_vars) #m+n+Kr
+        solution[mn+i] <- 1
+      } else {
+        values[i] <- NA_real_
+        solution <- NA_real_
+      }
+    } else {
+      values[i] <- lpSolveAPI::get.objective(lp)
+      solution <- lpSolveAPI::get.variables(lp)
+    }
+    sx[i,] <- solution[1:m]
+    sy[i,] <- solution[(m+1):mn]
+    lambda[i,] <- solution[mn1:n_vars]
+  }
+  colnames(sx) <- paste('sx',1:m,sep='')
+  colnames(sy) <- paste('sy',1:n,sep='')
+  list(
+    lp = lp,
+    values = values,
+    sx = sx,
+    sy = sy,
+    lambda = lambda
+  )
+}
+
+#' scale_vars (for slack analysis)
+#' @inheritParams create_lp
+#' @noRd
+scale_vars <- function(x, y, m, n, n_units) {
+  mmm <- colMeans(x)
+  nnn <- colMeans(y)
+  if (min(mmm) < 1e-4 || max(mmm) > 1e4 ||
+      min(nnn) < 1e-4 || max(nnn) > 1e4) {
+    x <- x / matrix(mmm, nrow=n_units, ncol=m, byrow=TRUE)
+    y <- y / matrix(nnn, nrow=n_units, ncol=n, byrow=TRUE)
+    scale <- TRUE
+  } else {
+    scale = FALSE
+  }
+  list(x=x, y=y, mmm = mmm, nnn = nnn, scale = scale)
+}
+
+#' adjust_values
+#' @noRd
+adjust_values <- function(unadj_values, lambda, eps) {
+  # efficiency scores close to 1 are converted to 1
+  values <- adjust_efficiency(unadj_values, eps)
+  # lambda values close to 0 and 1 are converted to 0 and 1 respectively
+  lambda <- adjust_lambda(lambda, eps)
+  list(
+    values = values,
+    unadj_values = unadj_values,
+    lambda = lambda)
+}
+
+#' adjust_slack
+#' @noRd
+adjust_slack <- function(lp, values, sx, sy, lambda, scaling, m, n, n_units){
+  lpcontr <- lpSolveAPI::lp.control(lp)
+  eps <- lpcontr$epsilon['epsint']
+  lambda <- adjust_lambda(lambda, eps)
+  sx[abs(sx) < eps] <- 0
+  sy[abs(sy) < eps] <- 0
+  if (scaling$scale) {
+    sx <- sx * matrix(scaling$mmm, nrow=n_units, ncol=m, byrow=TRUE)
+    sy <- sy * matrix(scaling$nnn, nrow=n_units, ncol=n, byrow=TRUE)
+  }
+  sum <- rowSums(sx) + rowSums(sy)
+  is_slack <- sum > eps
+  list(values = sum,
+       unadj_values = values,
+       lambda = lambda,
+       data = data.frame(sum = sum, is_slack, sx, sy)
+  )
+}
+
+#' adjust_efficiency
+#' @noRd
+adjust_efficiency <- function(x, eps){
+  x[abs(x-1) < eps] <- 1
+  x
+}
+
+#' adjust_lambda
+#' @noRd
+adjust_lambda <- function(x, eps) {
+  x[abs(x-1) < eps] <- 1   # 'lambda' close to 1 should be 1
+  x[abs(x) < eps] <- 0     # 'lambda' close to 0 should be 0
+  x
+}
+
+#' Get model dimensions
+#' @inheritParams create_lp
+#' @noRd
+get_dims <- function(x, y, xref = NULL, yref = NULL, rts, slack = FALSE) {
+  n_inputs <- ncol(x)
+  n_outputs <- ncol(y)
+  n_units <- if (!is.null(xref)) nrow(xref) else nrow(x)
+  if (rts == 'crs') {
+    n_lambda <- 0
+  } else if (rts %in% c('irs', 'drs')) {
+    n_lambda <- 1
+  } else {
+    n_lambda <- 2
+  }
+  if (slack) {
+    n_vars <- n_inputs + n_outputs + n_units
+  } else {
+    n_vars <- 1 + n_units
+  }
+  list(
+    n_inputs = n_inputs,
+    n_outputs = n_outputs,
+    n_units = n_units,
+    n_constraints = n_inputs + n_outputs + n_lambda,
+    n_vars = n_vars,
+    n_lambda = n_lambda
+  )
+}
+
+#' get_invalid_objective
+#' @noRd
+get_invalid_objective <- function(status, orientation) {
+  if (status != 0) {
+    if (status == 2 || status == 3) {
+      # if the unit is not in the technology set the value to Inf
+      if (orientation == 'in') return(Inf) else return(-Inf)
+    } else {
+      return(NA_real_)
+    }
+  }
+}
+
+#' get_lp_objective
+#' @noRd
+get_lp_objective <- function(lp, status, orientation) {
+  if (status != 0) {
+    get_invalid_objective(status, orientation)
+  } else {
+    lpSolveAPI::get.objective(lp)
+  }
+}
+
+#' get_lp_solution
+#' @noRd
+get_lp_solution <- function(lp, status) {
+  if (status != 0) {
+    n <- lpSolveAPI::dim.lpExtPtr(lp)[2]
+    rep(NA_real_, n)
+  } else {
+    lpSolveAPI::get.variables(lp)
+  }
+}
diff --git a/R/summary-dea.R b/R/summary-dea.R
index 5e8b8bf..210d507 100644
--- a/R/summary-dea.R
+++ b/R/summary-dea.R
@@ -25,6 +25,13 @@ summary_tbl_dea.Farrell <- function(x) {
   summary_tbl_dea(eff)
 }
 
+#' @method summary_tbl_dea pioneer_dea
+#' @export
+summary_tbl_dea.pioneer_dea <- function(x) {
+  eff <- x$results$efficiency
+  summary_tbl_dea(eff)
+}
+
 #' @method summary_tbl_dea numeric
 #' @export
 summary_tbl_dea.numeric <- function(x) {
diff --git a/inst/app/R/fct_dea-boostrap.R b/inst/app/R/fct_dea-bootstrap.R
similarity index 89%
rename from inst/app/R/fct_dea-boostrap.R
rename to inst/app/R/fct_dea-bootstrap.R
index 264274f..81b110a 100644
--- a/inst/app/R/fct_dea-boostrap.R
+++ b/inst/app/R/fct_dea-bootstrap.R
@@ -38,11 +38,11 @@ perform_boot <- function(x, y, rts, orientation, i, h, theta, boot) {
   beta <- bootstrap_sample(theta, h = h)
   if (orientation == 'in') {
     x_ref <- (theta / beta) * x
-    boot[, i] <- Benchmarking::dea(x, y, RTS = rts, ORIENTATION = orientation, XREF = x_ref, FAST = TRUE)
+    boot[, i] <- pioneeR:::compute_efficiency(x, y, rts = rts, orientation = orientation, xref = x_ref, yref = y, values_only = TRUE)$values
   } else if (orientation == 'out') {
     beta <- 1 / beta
     y_ref <- (theta / beta) * y
-    boot[, i] <- Benchmarking::dea(x, y, RTS = rts, ORIENTATION = orientation, YREF = y_ref, FAST = TRUE)
+    boot[, i] <- pioneeR:::compute_efficiency(x, y, rts = rts, orientation = orientation, xref = x, yref = y_ref, values_only = TRUE)$values
   }
   return(invisible(boot))
 }
diff --git a/inst/app/dea_analysis.Rmd b/inst/app/dea_analysis.Rmd
index 80ec932..24ce93f 100644
--- a/inst/app/dea_analysis.Rmd
+++ b/inst/app/dea_analysis.Rmd
@@ -64,12 +64,12 @@ knitr::kable(df)
 The technology for the model is `r mod$rts` and the orientation is `r mod$orient`.
 
 ```{r, echo=FALSE}
-tbl <- summary_tbl_dea(params$model$eff)
+tbl <- summary_tbl_dea(params$model$values)
 knitr::kable(tbl)
 ```
 
 ```{r, echo=FALSE}
-x <- summary(params$model$eff)
+x <- summary(params$model$values)
 x <- data.frame(as.list(x))
 colnames(x) <- gsub('^X', '', colnames(x))
 knitr::kable(x)
@@ -81,7 +81,7 @@ The efficiency scores have been estimated using the `dea` function from the `Ben
 
 ```{r, echo=FALSE}
 d <- data.frame(
-  `Efficiency score` = round(params$model$eff, params$settings$digits),
+  `Efficiency score` = round(params$model$values, params$settings$digits),
   check.names = FALSE
 )
 knitr::kable(d)
@@ -173,26 +173,35 @@ idvar <- params$params$id
 inputs <- c(params$inputs)
 outputs <- c(params$outputs)
 
-x <- create_matrix(df, inputs, idvar, normalize = params$params$normalize)
-y <- create_matrix(df, outputs, idvar, normalize = params$params$normalize)
-
 # Estimate the DEA model
-model <- Benchmarking::dea(
-  x, y, RTS = rts, ORIENTATION = orientation
+model <- compute_dea(
+  df,
+  idvar,
+  inputs,
+  outputs,
+  rts = rts,
+  orientation = orientation
 )
 
 # Print out the summary of the DEA model
 summary_tbl_dea(model)
 
-# Summary statistics for the efficiency scores
-summary(model$eff)
-
-# Estimate slack for all DMUs using the DEA model
-slack <- Benchmarking::slack(x, y, model)
+# Estimate the model with slack and super efficiency
+model <- compute_dea(
+  df,
+  idvar,
+  inputs,
+  outputs,
+  rts = rts,
+  orientation = orientation,
+  slack = TRUE,
+  super = TRUE
+)
 
-# Estimate super efficiency
-super_eff <- Benchmarking::sdea(x, y, RTS = rts, ORIENTATION = orientation)
+# Create input and output matrises for scale efficiency
+x <- create_matrix(df, inputs, idvar, normalize = params$params$normalize)
+y <- create_matrix(df, outputs, idvar, normalize = params$params$normalize)
 
 # Get table for scale efficiency
-compute_scale_efficiency(inputs, outputs, orientation, digits = 4L)
+compute_scale_efficiency(x, y, orientation, digits = 4L)
 ```
diff --git a/inst/app/global.R b/inst/app/global.R
index 6486c81..6ab1034 100644
--- a/inst/app/global.R
+++ b/inst/app/global.R
@@ -1,4 +1,5 @@
 # Load packages
+library(pioneeR)
 library(shiny)
 library(bslib)
 source('R/utils.R')
diff --git a/inst/app/server.R b/inst/app/server.R
index ec139ee..62782eb 100644
--- a/inst/app/server.R
+++ b/inst/app/server.R
@@ -1,9 +1,6 @@
 # Load required packages
 require(readxl)
 require(ucminf)
-require(lpSolveAPI)
-require(lpSolve)
-require(Benchmarking)
 require(productivity)
 require(ggplot2)
 require(scales)
@@ -254,9 +251,9 @@ shinyServer(function(input, output, session) {
 
     req(data(), dea.in(), dea.out())
 
-    d <- Benchmarking::dea(
-      dea.in(), dea.out(), RTS = model_params$rts,
-      ORIENTATION = model_params$orientation)
+    d <- pioneeR:::compute_efficiency(
+      dea.in(), dea.out(), rts = model_params$rts,
+      orientation = model_params$orientation)
 
     d
 
@@ -264,7 +261,7 @@ shinyServer(function(input, output, session) {
 
   dea.slack <- reactive({
     x <- tryCatch({
-      Benchmarking::slack(dea.in(), dea.out(), dea.prod())
+      pioneeR:::compute_slack(dea.in(), dea.out(), dea.prod())
     }, warning = function(e) {
       NULL
     }, error = function(e) {
@@ -277,9 +274,10 @@ shinyServer(function(input, output, session) {
 
     req(data(), dea.in(), dea.out())
 
-    d <- Benchmarking::sdea(
-      dea.in(), dea.out(), RTS = model_params$rts,
-      ORIENTATION = model_params$orientation)
+    d <- pioneeR:::compute_super_efficiency(
+      dea.in(), dea.out(), rts = model_params$rts,
+      orientation = model_params$orientation
+    )
 
     d
 
@@ -296,7 +294,7 @@ shinyServer(function(input, output, session) {
     if (is.matrix(x)) x <- rowSums(x)
     if (is.matrix(y)) y <- rowSums(y)
 
-    data.frame(dmu = rownames(dea.in()), x = unname(x), y = unname(y), eff = unname(prod$eff))
+    data.frame(dmu = rownames(dea.in()), x = unname(x), y = unname(y), eff = unname(prod$values))
   })
 
   dea_plot <- reactive({
@@ -380,7 +378,7 @@ shinyServer(function(input, output, session) {
 
   dea_scaleeff <- reactive({
 
-    out <- compute_scale_efficiency(dea.in(), dea.out(), model_params$orientation, 4L)
+    out <- compute_scale_efficiency(dea.in(), dea.out(), model_params$orientation, input$dea_round)
     out
 
   })
@@ -393,7 +391,7 @@ shinyServer(function(input, output, session) {
       dmu = selection()[[params()$id]],
       inputs = sapply(seq_len(nrow(dea.in())), function(i) sum(dea.in()[i,])),
       outputs = sapply(seq_len(nrow(dea.out())), function(i) sum(dea.out()[i,])),
-      eff = dea.prod()$eff,
+      eff = dea.prod()$values,
       stringsAsFactors = FALSE
     )
 
@@ -474,16 +472,15 @@ shinyServer(function(input, output, session) {
 
     req(dea.prod())
 
-    eff_tbl <- summary_tbl_dea(dea.prod())
-    eff <- dea.prod()$eff
-
+    eff <- dea.prod()$values
+    eff_tbl <- summary_tbl_dea(eff)
     if (model_params$orientation == 'in')
       sum_eff <- sum(dea.in() * eff) / sum(dea.in())
     else if (model_params$orientation == 'out')
       sum_eff <- sum(dea.out() * eff) / sum(dea.out())
 
     rts <- switch(
-      dea.prod()$RTS,
+      input$dea_rts,
       crs = 'Constant',
       vrs = 'Variable',
       drs = 'Non-increasing',
@@ -491,7 +488,7 @@ shinyServer(function(input, output, session) {
     )
 
     orient <- switch(
-      dea.prod()$ORIENTATION,
+      input$dea_orientation,
       'in' = 'Input',
       'out' = 'Output'
     )
@@ -582,7 +579,7 @@ shinyServer(function(input, output, session) {
 
   dea.tbl <- reactive({
 
-    deff <- matrix(dea.prod()$eff, ncol = 1, dimnames = list(NULL, 'Efficiency'))
+    deff <- matrix(dea.prod()$values, ncol = 1, dimnames = list(NULL, 'Efficiency'))
 
     # Initialize to NULL
     ins <- outs <- sl <- seff <- NULL
@@ -598,13 +595,13 @@ shinyServer(function(input, output, session) {
       outs <- apply(dea.out(), 1, sum)
 
     if (input$out.slack)
-      sl <- matrix(dea.slack()$sum, ncol = 1, dimnames = list(NULL, 'Slack'))
+      sl <- matrix(dea.slack()$values, ncol = 1, dimnames = list(NULL, 'Slack'))
 
     if (input$out.sdea)
-      seff <- matrix(sdea.prod()$eff, ncol = 1, dimnames = list(NULL, 'sDEA'))
+      seff <- matrix(sdea.prod()$values, ncol = 1, dimnames = list(NULL, 'sDEA'))
 
     df <- data.frame(
-      DMU = names(dea.prod()$eff),
+      DMU = selection()[, input$dea_id],
       round(cbind(ins, outs, deff, sl, seff), input$dea_round),
       stringsAsFactors = FALSE, row.names = NULL
     )
@@ -627,8 +624,7 @@ shinyServer(function(input, output, session) {
   })
 
   output$dea.slack <- renderReactable({
-    sl <- dea.slack()
-    df <- data.frame(round(cbind(sl$sx, sl$sy, sl$sum), 3))
+    df <- round(dea.slack()$data, input$dea_round)
     colnames(df)[ncol(df)] <- 'Total'
 
     opts <- list2(!!!reactable_opts, data = df)
@@ -636,11 +632,7 @@ shinyServer(function(input, output, session) {
   })
 
   output$peers.table <- renderReactable({
-    pr <- dea.prod()
-    pe <- Benchmarking::peers(pr, NAMES = TRUE)
-    df <- data.frame(
-      cbind(selection()[, input$dea_id], pe),
-      stringsAsFactors = FALSE)
+    df <- pioneeR:::get_peers(dea.prod(), ids = selection()[, input$dea_id], threshold = 0)
     colnames(df)[1] <- 'DMU'
 
     opts <- list2(!!!reactable_opts, data = df, columns = list(
@@ -692,13 +684,14 @@ shinyServer(function(input, output, session) {
 
   observeEvent(input$save_model, {
     if (length(models()) >= 10) return()
+    # EDIT HERE!
     mod <- dea.prod()
     mod_save <- list(
       id = rand_id(),
       data = data.frame(
-        idx = seq_len(length(mod$eff)),
-        dmu = names(mod$eff),
-        eff = round(unname(mod$eff), input$dea_round)
+        idx = seq_len(length(mod$values)),
+        dmu = selection()[, input$dea_id], #names(mod$eff),
+        eff = round(unname(mod$values), input$dea_round)
       ),
       params = list(
         rts = mod$RTS,
@@ -860,7 +853,7 @@ shinyServer(function(input, output, session) {
     y <- isolate(dea.out())
 
     withProgress(message = 'Running', value = 0, {
-      theta <- isolate(as.vector(dea.prod()$eff))
+      theta <- isolate(as.vector(dea.prod()$values))
       # theta >= 1 if 'out', <= if 'in'
       h <- if (is.numeric(bw_rule)) bw_rule else bw_rule(theta, rule = bw_rule)
       boot <- matrix(NA, nrow = nrow(x), ncol = b)
@@ -906,16 +899,16 @@ shinyServer(function(input, output, session) {
       return()
 
     # Add DMU names and round inputs
-    df <- cbind(data.frame(DMU = names(dea.prod()$eff)), round(res$tbl, input$boot_round))
+    df <- cbind(data.frame(DMU = names(dea.prod()$values)), round(res$tbl, input$boot_round))
 
     opts <- list2(!!!reactable_opts, data = df, columns = list(
-        eff = colDef(show = input$boot_show_eff, name = 'Efficiency'),
-        bias = colDef(show = input$boot_show_bias, name = 'Bias'),
-        eff_bc = colDef(name = 'Bias corr. score'),
-        lower = colDef(name = 'Lower bound'),
-        upper = colDef(name = 'Upper bound'),
-        range = colDef(name = 'CI range')
-      )
+      eff = colDef(show = input$boot_show_eff, name = 'Efficiency'),
+      bias = colDef(show = input$boot_show_bias, name = 'Bias'),
+      eff_bc = colDef(name = 'Bias corr. score'),
+      lower = colDef(name = 'Lower bound'),
+      upper = colDef(name = 'Upper bound'),
+      range = colDef(name = 'CI range')
+    )
     )
     tbl <- do.call(reactable, opts)
 
@@ -945,7 +938,7 @@ shinyServer(function(input, output, session) {
     },
     content = function(file) {
       res <- dea_boostrap()
-      df <- cbind(data.frame(DMU = names(dea.prod()$eff)), round(res$tbl, input$boot_round))
+      df <- cbind(data.frame(DMU = names(dea.prod()$values)), round(res$tbl, input$boot_round))
       # Export based on chosen file format
       if (input$boot_fileformat == 'dta') {
         colnames(df) <- gsub('\\s', '_', colnames(df))
diff --git a/man/compute_dea.Rd b/man/compute_dea.Rd
new file mode 100644
index 0000000..0f87a52
--- /dev/null
+++ b/man/compute_dea.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/compute_dea.R
+\name{compute_dea}
+\alias{compute_dea}
+\title{Compute DEA}
+\usage{
+compute_dea(
+  data,
+  id,
+  input,
+  output,
+  rts = c("crs", "vrs", "drs", "irs"),
+  orientation = c("in", "out"),
+  super = FALSE,
+  slack = FALSE,
+  peers = FALSE
+)
+}
+\arguments{
+\item{data}{Dataset to analyse.}
+
+\item{id}{A string with the DMU id or name variable.}
+
+\item{input}{A character vector with input variables.}
+
+\item{output}{A character vector with output variables.}
+
+\item{rts}{Returns to scale.
+\tabular{rlll}{
+\tab crs \tab\tab Constant returns to scale, convexity and free disposability. \cr
+\tab vrs \tab\tab Variable returns to scale, convexity and free disposability. \cr
+\tab drs \tab\tab Decreasing returns to scale, convexity, down-scaling and free disposability. \cr
+\tab irs \tab\tab Increasing returns to scale, (up-scaling, but not down-scaling), convexity and free disposability.
+}}
+
+\item{orientation}{Model orientation.}
+
+\item{super}{If \code{TRUE} super efficiency scores are calculated.}
+
+\item{slack}{If \code{TRUE} slack values are calculated.}
+
+\item{peers}{If \code{TRUE} peers are added to the response.}
+}
+\value{
+list
+}
+\description{
+DEA analysis.
+}
diff --git a/tests/script/create-testdata-benchmarking.R b/tests/script/create-testdata-benchmarking.R
new file mode 100644
index 0000000..a533394
--- /dev/null
+++ b/tests/script/create-testdata-benchmarking.R
@@ -0,0 +1,89 @@
+# ---- Load datasets ---- #
+
+# Simple example (taken from https://pystoned.readthedocs.io/en/latest/examples/DEA/dea_ref.html)
+x <- as.matrix(c(100, 200, 300, 500, 100, 200, 600, 400, 550, 600))
+y <- as.matrix(c(75, 100, 300, 400, 25, 50, 400, 260, 180, 240))
+xref <- as.matrix(c(100, 300, 500, 100, 600))
+yref <- as.matrix(c(75, 300, 400, 25, 400))
+
+# Frontier 4.1
+frontier41 <- readRDS('tests/testdata/frontier41.RDS')
+f41_x <- as.matrix(frontier41[c('capital', 'labour')])
+f41_y <- as.matrix(frontier41[c('output')])
+
+# Norwegian District Courts
+norCourts2018 <- readRDS('tests/testdata/norCourts2018.RDS')
+nc_x <- as.matrix(norCourts2018[c('judges', 'case_workers', 'costs')])
+nc_y <- as.matrix(norCourts2018[c(
+  'civil_cases', 'criminal_case_single',
+  'criminal_case_full_bench',
+  'other_civil_cases')])
+
+# Electric Plants
+electricPlants <- readRDS('tests/testdata/electricPlants.RDS')
+ecp_x <- as.matrix(electricPlants[c('labor', 'fuel', 'capital')])
+ecp_y <- as.matrix(electricPlants[c('output')])
+
+# Hospitals
+hospitals <- readRDS('tests/testdata/hospitals.RDS')
+hospitals_x <- as.matrix(hospitals[c('labor', 'capital')])
+hospitals_y <- as.matrix(hospitals[c('inpatients', 'outpatients')])
+
+
+# ---- Helper functions ---- #
+
+#' Create DEA test data (using Benchmarking)
+#' @param x A matrix with input values.
+#' @param y A matrix with output values.
+create_testdata <- function(x, y, xref = NULL, yref = NULL, ids) {
+  orientations <- c('in', 'out')
+  rts_ <- c('crs', 'vrs', 'drs', 'irs')
+  dl <- list()
+  for (a in seq_along(orientations)) {
+    for (b in seq_along(rts_)) {
+      tmp <- create_testdata_single(
+        x, y, xref, yref, ids,
+        rts = rts_[b],
+        orientation = orientations[a])
+      tmp <- list(tmp)
+      names(tmp) <- paste0(orientations[a], '_', rts_[b])
+      dl <- append(dl, tmp)
+    }
+  }
+  dl
+}
+
+#' Create DEA test data (single)
+#' @param x A matrix with input values.
+#' @param y A matrix with output values.
+#' @param rts Returns to scale
+#' @param orientation Model orientation
+create_testdata_single <- function(x, y, xref, yref, ids, rts, orientation) {
+  eff <- Benchmarking::dea(x, y, XREF = xref, YREF = yref, RTS = rts, ORIENTATION = orientation)
+  super_eff <- Benchmarking::sdea(x, y, RTS = rts, ORIENTATION = orientation)
+  slack <- Benchmarking::slack(x, y, e = eff)
+  colnames(eff$lambda) <- ids
+  peers <- Benchmarking::peers(eff, NAMES = TRUE)
+
+  dimnames(eff$lambda) <- NULL
+  dimnames(super_eff$lambda) <- NULL
+  dimnames(slack$lambda) <- NULL
+
+  dl <- list(efficiency = list(eff = eff$eff, lambda = eff$lambda, objval = eff$objval, peers = peers),
+             super_efficiency = list(eff = super_eff$eff, lambda = super_eff$lambda, objval = super_eff$objval),
+             slack = list(eff = slack$eff, lambda = slack$lambda, objval = slack$objval, sum = slack$sum)
+  )
+  dl
+}
+
+
+# ---- Create testdata ---- #
+
+results <- list(
+  'simple' = create_testdata(x, y, xref, yref, ids = rownames(x)),
+  'frontier41' = create_testdata(f41_x, f41_y, ids = frontier41$firm),
+  'norCourts2018' = create_testdata(nc_x, nc_y, ids = norCourts2018$district_court),
+  'electricPlants' = create_testdata(ecp_x, ecp_y, ids = electricPlants$plant),
+  'hospitals' = create_testdata(hospitals_x, hospitals_y, ids = hospitals$firm_id)
+)
+saveRDS(results, 'tests/testdata/benchmarking-results.RDS')
diff --git a/tests/testdata/benchmarking-results.RDS b/tests/testdata/benchmarking-results.RDS
new file mode 100644
index 0000000..0f9f2d3
Binary files /dev/null and b/tests/testdata/benchmarking-results.RDS differ
diff --git a/tests/testdata/electricPlants.RDS b/tests/testdata/electricPlants.RDS
new file mode 100644
index 0000000..9e7ea41
Binary files /dev/null and b/tests/testdata/electricPlants.RDS differ
diff --git a/tests/testdata/frontier41.RDS b/tests/testdata/frontier41.RDS
new file mode 100644
index 0000000..8a1f33c
Binary files /dev/null and b/tests/testdata/frontier41.RDS differ
diff --git a/tests/testdata/hospitals.RDS b/tests/testdata/hospitals.RDS
new file mode 100644
index 0000000..9f0a270
Binary files /dev/null and b/tests/testdata/hospitals.RDS differ
diff --git a/tests/testdata/norCourts2018.RDS b/tests/testdata/norCourts2018.RDS
new file mode 100644
index 0000000..cb0fa6c
Binary files /dev/null and b/tests/testdata/norCourts2018.RDS differ
diff --git a/tests/testthat/test-compute_efficiency.R b/tests/testthat/test-compute_efficiency.R
new file mode 100644
index 0000000..5dc49b8
--- /dev/null
+++ b/tests/testthat/test-compute_efficiency.R
@@ -0,0 +1,343 @@
+# Benchmarking results
+benchmarking_results <- readRDS('../testdata/benchmarking-results.RDS')
+
+# Frontier 4.1
+frontier41 <- readRDS('../testdata/frontier41.RDS')
+f41_x <- as.matrix(frontier41[c('capital', 'labour')])
+f41_y <- as.matrix(frontier41[c('output')])
+
+# Norwegian District Courts
+norCourts2018 <- readRDS('../testdata/norCourts2018.RDS')
+nc_x <- as.matrix(norCourts2018[c('judges', 'case_workers', 'costs')])
+nc_y <- as.matrix(norCourts2018[c(
+  'civil_cases', 'criminal_case_single',
+  'criminal_case_full_bench',
+  'other_civil_cases')])
+
+# Electric Plants
+electricPlants <- readRDS('../testdata/electricPlants.RDS')
+ecp_x <- as.matrix(electricPlants[c('labor', 'fuel', 'capital')])
+ecp_y <- as.matrix(electricPlants[c('output')])
+
+# Hospitals
+hospitals <- readRDS('../testdata/hospitals.RDS')
+hp_x <- as.matrix(hospitals[c('labor', 'capital')])
+hp_y <- as.matrix(hospitals[c('inpatients', 'outpatients')])
+
+test_that('compute_efficiency() returns the correct structure', {
+
+  res <- compute_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'out')
+  # class
+  expect_identical(class(res), 'list')
+  # object names
+  expect_identical(names(res), c('values', 'unadj_values', 'lambda', 'info'))
+  expect_identical(names(res$info), c('rts', 'orientation', 'dims'))
+  expect_identical(names(res$info$dims), c('n_inputs', 'n_outputs', 'n_units', 'n_constraints', 'n_vars', 'n_lambda'))
+  expect_identical(res$info$rts, 'vrs')
+  expect_identical(res$info$orientation, 'out')
+  # dimensions (dim object)
+  expect_equal(res$info$dims$n_units, nrow(f41_x))
+  expect_equal(res$info$dims$n_inputs, ncol(f41_x))
+  expect_equal(res$info$dims$n_outputs, ncol(f41_y))
+  # dimensions (values)
+  expect_equal(length(res$values), nrow(f41_x))
+  expect_equal(length(res$unadj_values), nrow(f41_x))
+  expect_equal(nrow(res$lambda), nrow(f41_x))
+  expect_equal(ncol(res$lambda), nrow(f41_x))
+
+  # values_only = TRUE
+  res <- compute_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'out', values_only = TRUE)
+  expect_identical(class(res), 'list')
+  expect_identical(names(res), c('values'))
+
+})
+
+test_that('compute_efficiency() works for CRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_crs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_crs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_crs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_crs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_crs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_crs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_crs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_crs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
+test_that('compute_efficiency() works for VRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_vrs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_vrs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_vrs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_vrs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_vrs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_vrs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_vrs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_vrs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+
+})
+
+test_that('compute_efficiency() works for IRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_irs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_irs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_irs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_irs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_irs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_irs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_irs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_irs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+
+})
+
+test_that('compute_efficiency() works for DRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_drs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_drs$efficiency
+  res <- compute_efficiency(f41_x, f41_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_drs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_drs$efficiency
+  res <- compute_efficiency(nc_x, nc_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_drs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_drs$efficiency
+  res <- compute_efficiency(hp_x, hp_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_drs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_drs$efficiency
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
+test_that('compute_efficiency() works when a reference techology set is used', {
+
+  # Simple example (taken from https://pystoned.readthedocs.io/en/latest/examples/DEA/dea_ref.html)
+  x <- as.matrix(c(100, 200, 300, 500, 100, 200, 600, 400, 550, 600))
+  y <- as.matrix(c(75, 100, 300, 400, 25, 50, 400, 260, 180, 240))
+  xref <- as.matrix(c(100, 300, 500, 100, 600))
+  yref <- as.matrix(c(75, 300, 400, 25, 400))
+
+  bench_res <- benchmarking_results$simple$out_vrs$efficiency
+  res <- compute_efficiency(x, y, xref, yref, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
diff --git a/tests/testthat/test-compute_scale_efficiency.R b/tests/testthat/test-compute_scale_efficiency.R
new file mode 100644
index 0000000..6e3f147
--- /dev/null
+++ b/tests/testthat/test-compute_scale_efficiency.R
@@ -0,0 +1,27 @@
+df <- data.frame(
+  id = letters[1:20],
+  a = 1:20,
+  b = 21:40,
+  c = 41:60,
+  d = 61:80
+)
+
+test_that('compute_scale_efficiency returns expected results', {
+  withr::local_options(lifecycle_verbosity = 'quiet')
+  x <- create_matrix(df, id = 'id', columns = c('a', 'b', 'c'))
+  y <- create_matrix(df, id = 'id', columns = 'd')
+  res <- compute_scale_efficiency(x, y, orientation = 'in')
+  expect_true(is.data.frame(res))
+  res <- compute_scale_efficiency(x, y, orientation = 'in', digits = 4L)
+  expect_true(is.data.frame(res))
+})
+
+test_that('compute_scale_efficiency returns errors as expected', {
+  withr::local_options(lifecycle_verbosity = 'quiet')
+  x <- create_matrix(df, id = 'id', columns = c('a', 'b', 'c'))
+  y <- create_matrix(df, id = 'id', columns = 'd')
+  expect_error(compute_scale_efficiency(x, 1:3))
+  expect_error(compute_scale_efficiency(1:3, x))
+  expect_error(compute_scale_efficiency(matrix(1:4, ncol = 2), matrix(1:6, ncol = 2)))
+  expect_warning(compute_scale_efficiency(x, y, digits = 'a'))
+})
diff --git a/tests/testthat/test-compute_slack.R b/tests/testthat/test-compute_slack.R
new file mode 100644
index 0000000..b08989d
--- /dev/null
+++ b/tests/testthat/test-compute_slack.R
@@ -0,0 +1,357 @@
+# Benchmarking results
+benchmarking_results <- readRDS('../testdata/benchmarking-results.RDS')
+
+# Frontier 4.1
+frontier41 <- readRDS('../testdata/frontier41.RDS')
+f41_x <- as.matrix(frontier41[c('capital', 'labour')])
+f41_y <- as.matrix(frontier41[c('output')])
+
+# Norwegian District Courts
+norCourts2018 <- readRDS('../testdata/norCourts2018.RDS')
+nc_x <- as.matrix(norCourts2018[c('judges', 'case_workers', 'costs')])
+nc_y <- as.matrix(norCourts2018[c(
+  'civil_cases', 'criminal_case_single',
+  'criminal_case_full_bench',
+  'other_civil_cases')])
+
+# Electric Plants
+electricPlants <- readRDS('../testdata/electricPlants.RDS')
+ecp_x <- as.matrix(electricPlants[c('labor', 'fuel', 'capital')])
+ecp_y <- as.matrix(electricPlants[c('output')])
+
+# Hospitals
+hospitals <- readRDS('../testdata/hospitals.RDS')
+hp_x <- as.matrix(hospitals[c('labor', 'capital')])
+hp_y <- as.matrix(hospitals[c('inpatients', 'outpatients')])
+
+
+test_that('compute_slack() returns the correct structure', {
+  res <- compute_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'out')
+  res <- compute_slack(f41_x, f41_y, model = res)
+  # class
+  expect_identical(class(res), 'list')
+  expect_identical(class(res$data), 'data.frame')
+  # object names
+  expect_identical(names(res), c('values', 'unadj_values', 'lambda', 'data', 'info'))
+  expect_identical(names(res$data), c('sum', 'is_slack', 'sx1', 'sx2', 'sy1'))
+  expect_identical(names(res$info), c('rts', 'orientation', 'dims'))
+  expect_identical(names(res$info$dims), c('n_inputs', 'n_outputs', 'n_units', 'n_constraints', 'n_vars', 'n_lambda'))
+  expect_identical(res$info$rts, 'vrs')
+  expect_identical(res$info$orientation, 'out')
+  # dimensions (dim object)
+  expect_equal(res$info$dims$n_units, nrow(f41_x))
+  expect_equal(res$info$dims$n_inputs, ncol(f41_x))
+  expect_equal(res$info$dims$n_outputs, ncol(f41_y))
+  # dimensions (values)
+  expect_equal(length(res$values), nrow(f41_x))
+  expect_equal(length(res$unadj_values), nrow(f41_x))
+  expect_equal(nrow(res$lambda), nrow(f41_x))
+  expect_equal(ncol(res$lambda), nrow(f41_x))
+
+})
+
+test_that('compute_slack() works for CRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_crs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'crs', orientation = 'in')
+  res <- compute_slack(f41_x, f41_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_crs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'crs', orientation = 'out')
+  res <- compute_slack(f41_x, f41_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_crs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'crs', orientation = 'in')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_crs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'crs', orientation = 'out')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_crs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'crs', orientation = 'in')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_crs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'crs', orientation = 'out')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_crs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'crs', orientation = 'in')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_crs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'crs', orientation = 'out')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
+test_that('compute_slack() works for VRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_vrs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'in')
+  res <- compute_slack(f41_x, f41_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_vrs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'out')
+  res <- compute_slack(f41_x, f41_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_vrs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'vrs', orientation = 'in')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_vrs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'vrs', orientation = 'out')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_vrs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'vrs', orientation = 'in')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_vrs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'vrs', orientation = 'out')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_vrs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'vrs', orientation = 'in')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_vrs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'vrs', orientation = 'out')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
+test_that('compute_slack() works for IRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_irs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'irs', orientation = 'in')
+  res <- compute_slack(f41_x, f41_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_irs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'irs', orientation = 'out')
+  res <- compute_slack(f41_x, f41_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_irs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'irs', orientation = 'in')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_irs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'irs', orientation = 'out')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_irs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'irs', orientation = 'in')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_irs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'irs', orientation = 'out')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_irs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'irs', orientation = 'in')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_irs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'irs', orientation = 'out')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+
+})
+
+test_that('compute_slack() works for DRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_drs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'drs', orientation = 'in')
+  res <- compute_slack(f41_x, f41_y, model = res )
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_drs$slack
+  res <- compute_efficiency(f41_x, f41_y, rts = 'drs', orientation = 'out')
+  res <- compute_slack(f41_x, f41_y, model = res )
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_drs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'drs', orientation = 'in')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_drs$slack
+  res <- compute_efficiency(nc_x, nc_y, rts = 'drs', orientation = 'out')
+  res <- compute_slack(nc_x, nc_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_drs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'drs', orientation = 'in')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_drs$slack
+  res <- compute_efficiency(hp_x, hp_y, rts = 'drs', orientation = 'out')
+  res <- compute_slack(hp_x, hp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_drs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'drs', orientation = 'in')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_drs$slack
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'drs', orientation = 'out')
+  res <- compute_slack(ecp_x, ecp_y, model = res)
+  expect_equal(res$values, bench_res$sum) # slack sum
+  expect_equal(res$unadj_values, bench_res$objval) # unadjusted slack sum
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
diff --git a/tests/testthat/test-compute_super_efficiency.R b/tests/testthat/test-compute_super_efficiency.R
new file mode 100644
index 0000000..7c68edd
--- /dev/null
+++ b/tests/testthat/test-compute_super_efficiency.R
@@ -0,0 +1,291 @@
+# Benchmarking results
+benchmarking_results <- readRDS('../testdata/benchmarking-results.RDS')
+
+# Frontier 4.1
+frontier41 <- readRDS('../testdata/frontier41.RDS')
+f41_x <- as.matrix(frontier41[c('capital', 'labour')])
+f41_y <- as.matrix(frontier41[c('output')])
+
+# Norwegian District Courts
+norCourts2018 <- readRDS('../testdata/norCourts2018.RDS')
+nc_x <- as.matrix(norCourts2018[c('judges', 'case_workers', 'costs')])
+nc_y <- as.matrix(norCourts2018[c(
+  'civil_cases', 'criminal_case_single',
+  'criminal_case_full_bench',
+  'other_civil_cases')])
+
+# Electric Plants
+electricPlants <- readRDS('../testdata/electricPlants.RDS')
+ecp_x <- as.matrix(electricPlants[c('labor', 'fuel', 'capital')])
+ecp_y <- as.matrix(electricPlants[c('output')])
+
+# Hospitals
+hospitals <- readRDS('../testdata/hospitals.RDS')
+hp_x <- as.matrix(hospitals[c('labor', 'capital')])
+hp_y <- as.matrix(hospitals[c('inpatients', 'outpatients')])
+
+
+test_that('compute_super_efficiency() returns the correct structure', {
+
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'out')
+  # class
+  expect_identical(class(res), 'list')
+  # object names
+  expect_identical(names(res), c('values', 'unadj_values', 'lambda', 'info'))
+  expect_identical(names(res$info), c('rts', 'orientation', 'dims'))
+  expect_identical(names(res$info$dims), c('n_inputs', 'n_outputs', 'n_units', 'n_constraints', 'n_vars', 'n_lambda'))
+  expect_identical(res$info$rts, 'vrs')
+  expect_identical(res$info$orientation, 'out')
+  # dimensions (dim object)
+  expect_equal(res$info$dims$n_units, nrow(f41_x))
+  expect_equal(res$info$dims$n_inputs, ncol(f41_x))
+  expect_equal(res$info$dims$n_outputs, ncol(f41_y))
+  # dimensions (values)
+  expect_equal(length(res$values), nrow(f41_x))
+  expect_equal(length(res$unadj_values), nrow(f41_x))
+  expect_equal(nrow(res$lambda), nrow(f41_x))
+  expect_equal(ncol(res$lambda), nrow(f41_x))
+
+})
+
+test_that('compute_super_efficiency() works for CRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_crs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_crs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_crs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_crs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_crs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_crs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_crs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'crs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_crs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'crs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
+test_that('compute_super_efficiency() works for VRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_vrs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_vrs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_vrs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_vrs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_vrs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_vrs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_vrs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'vrs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_vrs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'vrs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+
+})
+
+test_that('compute_super_efficiency() works for IRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_irs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_irs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_irs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_irs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_irs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_irs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_irs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'irs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_irs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'irs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
+test_that('compute_super_efficiency() works for DRS', {
+
+  # --- Frontier 4.1 data --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$frontier41$in_drs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$frontier41$out_drs$super_efficiency
+  res <- compute_super_efficiency(f41_x, f41_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Norwegian District Courts --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$norCourts2018$in_drs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$norCourts2018$out_drs$super_efficiency
+  res <- compute_super_efficiency(nc_x, nc_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Hospitals --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$hospitals$in_drs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$hospitals$out_drs$super_efficiency
+  res <- compute_super_efficiency(hp_x, hp_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # --- Electric Plants --- #
+
+  # orientation in
+  bench_res <- benchmarking_results$electricPlants$in_drs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'drs', orientation = 'in')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+  # orientation out
+  bench_res <- benchmarking_results$electricPlants$out_drs$super_efficiency
+  res <- compute_super_efficiency(ecp_x, ecp_y, rts = 'drs', orientation = 'out')
+  expect_equal(res$values, bench_res$eff) # efficiency
+  expect_equal(res$lambda, bench_res$lambda) # lambda
+
+})
+
diff --git a/tests/testthat/test-dea-utils.R b/tests/testthat/test-dea-utils.R
index b0230da..15f08df 100644
--- a/tests/testthat/test-dea-utils.R
+++ b/tests/testthat/test-dea-utils.R
@@ -6,38 +6,50 @@ df <- data.frame(
   d = 61:80
 )
 
-test_that('create_matrix returns expected results', {
+test_that('create_matrix() returns expected results', {
+  withr::local_options(lifecycle_verbosity = 'quiet')
   x <- create_matrix(df, id = 'id', columns = c('a', 'b', 'c'))
   expect_true(is.matrix(x))
   expect_equal(c('a', 'b', 'c'), colnames(x))
   expect_equal(letters[1:20], rownames(x))
 })
 
-test_that('normalization works in create_matrix', {
+test_that('normalization works in create_matrix()', {
+  withr::local_options(lifecycle_verbosity = 'quiet')
   x <- create_matrix(df, id = 'id', columns = c('a', 'b', 'c'), normalize = TRUE)
   expect_equal(mean(x), 1)
 })
 
-test_that('create_matrix returns errors as expected', {
+test_that('create_matrix() returns errors as expected', {
+  withr::local_options(lifecycle_verbosity = 'quiet')
   expect_error(create_matrix(1:10, 1, 1))
   expect_error(create_matrix(df, id = 'id', columns = 'x'))
   expect_error(create_matrix(df, id = 'x', columns = 'd'))
 })
 
-test_that('compute_scale_efficiency returns expected results', {
-  x <- create_matrix(df, id = 'id', columns = c('a', 'b', 'c'))
-  y <- create_matrix(df, id = 'id', columns = 'd')
-  res <- compute_scale_efficiency(x, y, orientation = 'in')
-  expect_true(is.data.frame(res))
-  res <- compute_scale_efficiency(x, y, orientation = 'in', digits = 4L)
-  expect_true(is.data.frame(res))
+
+test_that('check_numeric() works correctly', {
+  # numeric data
+  numeric_x <- data.frame(a = c(1, 2, 3), b = c(4, 5, 6))
+  numeric_y <- data.frame(c = c(7, 8, 9), d = c(10, 11, 12))
+  expect_silent(check_numeric(numeric_x, numeric_y))
+
+  # non-numeric data
+  non_numeric_x <- data.frame(a = c(1, 2, 3), b = c('a', 'b', 'c'))
+  numeric_y <- data.frame(a = c(1, 2, 3), b = c(4, 5, 6))
+  expect_error(check_numeric(non_numeric_x, numeric_y), 'All variables must be numeric.')
+  expect_error(check_numeric(numeric_y, non_numeric_x), 'All variables must be numeric.')
 })
 
-test_that('compute_scale_efficiency returns errors as expected', {
-  x <- create_matrix(df, id = 'id', columns = c('a', 'b', 'c'))
-  y <- create_matrix(df, id = 'id', columns = 'd')
-  expect_error(compute_scale_efficiency(x, 1:3))
-  expect_error(compute_scale_efficiency(1:3, x))
-  expect_error(compute_scale_efficiency(matrix(1:4, ncol = 2), matrix(1:6, ncol = 2)))
-  expect_warning(compute_scale_efficiency(x, y, digits = 'a'))
+test_that('check_nunits() works correctly', {
+  # correct
+  x <- data.frame(a = c(1, 2, 3), b = c(4, 5, 6))
+  y <- data.frame(c = c(7, 8, 9), d = c(10, 11, 12))
+  expect_silent(check_nunits(x, y))
+
+  # incorrect
+  x <- data.frame(a = c(1, 2, 3), b = c(4, 5, 6))
+  y <- data.frame(c = c(7, 8, 9, 10), d = c(11, 12, 13, 14))
+  expect_error(check_nunits(x, y))
 })
+
diff --git a/tests/testthat/test-get_peers.r b/tests/testthat/test-get_peers.r
new file mode 100644
index 0000000..2f3dc5c
--- /dev/null
+++ b/tests/testthat/test-get_peers.r
@@ -0,0 +1,68 @@
+# Benchmarking results
+benchmarking_results <- readRDS('../testdata/benchmarking-results.RDS')
+
+# Frontier 4.1
+frontier41 <- readRDS('../testdata/frontier41.RDS')
+f41_x <- as.matrix(frontier41[c('capital', 'labour')])
+f41_y <- as.matrix(frontier41[c('output')])
+
+# Norwegian District Courts
+norCourts2018 <- readRDS('../testdata/norCourts2018.RDS')
+nc_x <- as.matrix(norCourts2018[c('judges', 'case_workers', 'costs')])
+nc_y <- as.matrix(norCourts2018[c(
+  'civil_cases', 'criminal_case_single',
+  'criminal_case_full_bench',
+  'other_civil_cases')])
+
+# Electric Plants
+electricPlants <- readRDS('../testdata/electricPlants.RDS')
+ecp_x <- as.matrix(electricPlants[c('labor', 'fuel', 'capital')])
+ecp_y <- as.matrix(electricPlants[c('output')])
+
+# Hospitals
+hospitals <- readRDS('../testdata/hospitals.RDS')
+hp_x <- as.matrix(hospitals[c('labor', 'capital')])
+hp_y <- as.matrix(hospitals[c('inpatients', 'outpatients')])
+
+test_that('get_peers() gives the correct results', {
+
+  # --- Frontier 4.1 data --- #
+
+  bench_res <- benchmarking_results$frontier41$out_crs$efficiency$peers
+  res <- compute_efficiency(f41_x, f41_y, rts = 'crs', orientation = 'out')
+  df <- get_peers(res, ids = as.character(frontier41$firm))
+  expect_identical(df$peer1, bench_res[,1])
+  expect_identical(df$peer2, bench_res[,2])
+
+  # --- Norwegian District Courts --- #
+
+  bench_res <- benchmarking_results$norCourts2018$out_vrs$efficiency$peers
+  res <- compute_efficiency(nc_x, nc_y, rts = 'vrs', orientation = 'out')
+  df <- get_peers(res, ids = norCourts2018$district_court)
+  expect_identical(df$peer1, bench_res[,1])
+  expect_identical(df$peer2, bench_res[,2])
+  expect_identical(df$peer3, bench_res[,3])
+  expect_identical(df$peer4, bench_res[,4])
+  expect_identical(df$peer5, bench_res[,5])
+  expect_identical(df$peer6, bench_res[,6])
+  expect_identical(df$peer7, bench_res[,7])
+
+  # --- Hospitals --- #
+
+  bench_res <- benchmarking_results$hospitals$out_vrs$efficiency$peers
+  res <- compute_efficiency(hp_x, hp_y, rts = 'vrs', orientation = 'out')
+  df <- get_peers(res, ids = as.character(hospitals$firm_id))
+  expect_identical(df$peer1, bench_res[,1])
+  expect_identical(df$peer2, bench_res[,2])
+  expect_identical(df$peer3, bench_res[,3])
+  expect_identical(df$peer4, bench_res[,4])
+
+  # --- Electric Plants --- #
+
+  bench_res <- benchmarking_results$electricPlants$out_drs$efficiency$peers
+  res <- compute_efficiency(ecp_x, ecp_y, rts = 'drs', orientation = 'out')
+  df <- get_peers(res, ids = electricPlants$plant)
+  expect_identical(df$peer1, bench_res[,1])
+  expect_identical(df$peer2, bench_res[,2])
+
+})
diff --git a/tests/testthat/tests-lp-utils.R b/tests/testthat/tests-lp-utils.R
new file mode 100644
index 0000000..a19ddaa
--- /dev/null
+++ b/tests/testthat/tests-lp-utils.R
@@ -0,0 +1,81 @@
+test_that('adjust_efficiency() works correctly', {
+  eps <- 0.02
+  x <- c(0.99, 1.01, 0.5)
+  expected <- c(1, 1, 0.5)
+  result <- adjust_efficiency(x, eps)
+  expect_equal(result, expected)
+
+  eps <- 1e-06
+  x <- c(0.9999999, 1.000001, 0.5)
+  expected <- c(1, 1, 0.5)
+  result <- adjust_efficiency(x, eps)
+  expect_equal(result, expected)
+})
+
+test_that('adjust_lambda() works correctly', {
+  eps <- 0.02
+  x <- c(0.99, 1.01, 0.5, 0.009)
+  expected <- c(1, 1, 0.5, 0)
+  result <- adjust_lambda(x, eps)
+  expect_equal(result, expected)
+
+  eps <- 1e-06
+  x <- c(0.9999999, 1.000001, 0.5, 0.0000009)
+  expected <- c(1, 1, 0.5, 0)
+  result <- adjust_lambda(x, eps)
+  expect_equal(result, expected)
+})
+
+test_that('get_dims() works correctly', {
+  x <- matrix(1:6, nrow=3, ncol=2)
+  y <- matrix(7:9, nrow=3, ncol=1)
+  xref <- matrix(1:4, nrow=2, ncol=2)
+  yref <- matrix(5:6, nrow=2, ncol=1)
+
+  # crs
+  result <- get_dims(x, y, rts = 'crs')
+  expect_equal(result$n_inputs, 2L)
+  expect_equal(result$n_outputs, 1L)
+  expect_equal(result$n_units, 3L)
+  expect_equal(result$n_constraints, 3L)
+  expect_equal(result$n_lambda, 0L)
+  expect_equal(result$n_vars, 4L)
+
+  # irs/drs
+  result <- get_dims(x, y, rts = 'irs')
+  expect_equal(result$n_inputs, 2L)
+  expect_equal(result$n_outputs, 1L)
+  expect_equal(result$n_units, 3L)
+  expect_equal(result$n_constraints, 4L)
+  expect_equal(result$n_lambda, 1L)
+  expect_equal(result$n_vars, 4L)
+
+  # vrs
+  result <- get_dims(x, y, rts = 'vrs')
+  expect_equal(result$n_inputs, 2L)
+  expect_equal(result$n_outputs, 1L)
+  expect_equal(result$n_units, 3L)
+  expect_equal(result$n_constraints, 5L)
+  expect_equal(result$n_lambda, 2L)
+  expect_equal(result$n_vars, 4L)
+
+  # xref/yref (only two units here)
+  result <- get_dims(x, y, xref, yref, rts = 'crs')
+  expect_equal(result$n_units, 2L)
+  expect_equal(result$n_vars, 3L)
+
+  # slack (LP matrix is transposed in this calculation, hence n_vars = n_inputs + n_outputs + n_units)
+  result <- get_dims(x, y, rts = 'crs', slack = TRUE)
+  expect_equal(result$n_vars, 6L)
+
+})
+
+test_that('get_invalid_objective() works correctly', {
+  expect_equal(get_invalid_objective(0, 'in'), NULL)
+  expect_equal(get_invalid_objective(2, 'in'), Inf)
+  expect_equal(get_invalid_objective(3, 'in'), Inf)
+  expect_equal(get_invalid_objective(2, 'out'), -Inf)
+  expect_equal(get_invalid_objective(3, 'out'), -Inf)
+  expect_equal(get_invalid_objective(1, 'in'), NA_real_)
+  expect_equal(get_invalid_objective(1, 'out'), NA_real_)
+})