diff --git a/NEWS.md b/NEWS.md
index fafbb61..3a20e96 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,13 @@
+# melt 1.11.0
+## MINOR IMPROVEMENTS
+* Updated package vignette with the publication in the Journal of Statistical Software.
+
+## DEPRECATED AND DEFUNCT
+* Removed `el_pairwise()` and associated methods.
+
+* Removed `sigTests()` for objects inheriting from `SummaryLM`.
+
+
 # melt 1.10.0
 ## NEW FEATURES
 * `el_glm()` accepts `quasipoisson` family with `"sqrt"` link function for the argument `family`.
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 5841808..dce0fcc 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -1,4 +1,4 @@
-url: https://ropensci.github.io/melt/
+destination: docs/
 reference:
 - title: Functions (Constructors)
   contents:
diff --git a/deprecated.R b/deprecated.R
deleted file mode 100644
index 3d9c2c0..0000000
--- a/deprecated.R
+++ /dev/null
@@ -1,190 +0,0 @@
-#' Pairwise comparisons for general block designs
-#'
-#' @description Tests all pairwise comparisons or comparisons with control for
-#'  general block designs. Two single step asymptotic \eqn{k}-FWER (generalized
-#'  family-wise error rate) controlling procedures are available: asymptotic
-#'  Monte Carlo (AMC) and nonparametric bootstrap (NB).
-#'
-#' \emph{This function is deprecated and will be removed in a future release.}
-#'
-#' @param formula An object of class [`formula`] (or one that can be coerced to
-#'   that class) for a symbolic description of the model to be fitted. It must
-#'   specify variables for response, treatment, and block as 'response ~
-#'   treatment | block'. Note that the use of vertical bar (|) separating
-#'   treatment and block.
-#' @param data A data frame, list or environment (or object coercible by
-#'   [as.data.frame()] to a data frame) containing the variables in `formula`.
-#' @param control An optional single character that specifies the treatment for
-#'   comparisons with control.
-#' @param k A single integer for \eqn{k} in \eqn{k}-FWER. Defaults to 1.
-#' @param alpha A single numeric for the overall significance level. Defaults to
-#'   `0.05`.
-#' @param method A single character for the procedure to be used; either `AMC`
-#'   or `NB` is supported. Defaults to `AMC`.
-#' @param B A single integer for the number of Monte Carlo samples for the AMC
-#'   (number of bootstrap replicates for the NB).
-#' @param nthreads A single integer for the number of threads for parallel
-#'   computation via OpenMP (if available). Defaults to `1`.
-#' @param maxit A single integer for the maximum number of iterations for
-#'   constrained minimization of empirical likelihood. Defaults to `10000`.
-#' @param abstol A single numeric for the the absolute convergence tolerance for
-#'   optimization. Defaults to `1e-08`.
-#' @return A list with class \code{c("pairwise", "melt")}.
-#' @references
-#'   Kim E, MacEachern SN, Peruggia M (2023).
-#'   ``Empirical likelihood for the analysis of experimental designs.''
-#'   \emph{Journal of Nonparametric Statistics}, **35**(4), 709--732.
-#'   \doi{10.1080/10485252.2023.2206919}.
-#' @examples
-#' \dontrun{
-#' # All pairwise comparisons
-#' data("clothianidin")
-#' el_pairwise(clo ~ trt | blk, data = clothianidin, B = 1000)
-#'
-#' # Comparisons with control
-#' el_pairwise(clo ~ trt | blk,
-#'   data = clothianidin, control = "Naked", method = "NB", B = 500
-#' )
-#' }
-#' @importFrom stats terms reshape
-#' @keywords internal
-#' @export
-el_pairwise <- function(formula,
-                        data,
-                        control = NULL,
-                        k = 1L,
-                        alpha = 0.05,
-                        method = c("AMC", "NB"),
-                        B,
-                        nthreads = 1L,
-                        maxit = 10000L,
-                        abstol = 1e-08) {
-  .Deprecated(msg = "`el_pairwise()` is deprecated in melt v1.5.2.")
-  alpha <- validate_alpha(alpha)
-  B <- validate_b(B)
-  max_threads <- get_max_threads()
-  nthreads <- validate_nthreads(nthreads, max_threads)
-  maxit <- validate_maxit(maxit)
-  abstol <- validate_tol(abstol)
-  method <- match.arg(method)
-  f <- attributes(terms(formula))
-  stopifnot(
-    "`formula` must be specified as `response ~ treatment | block`." =
-      (isTRUE(f$response == 1L && f$intercept == 1L)),
-    "`formula` must be specified as `response ~ treatment | block`." =
-      (isTRUE(length(f$variables) == 3L)),
-    "`formula` must be specified as `response ~ treatment | block`." =
-      (isTRUE(typeof(f$variables[[3L]]) == "language" &&
-        length(f$variables[[3L]]) == 3L)),
-    "`|` operator is missing." = (isTRUE(f$variables[[3L]][[1L]] == "|")),
-    "Transformation on variables is not allowed." =
-      (isTRUE(typeof(f$variables[[3L]][[2L]]) == "symbol" &&
-        typeof(f$variables[[3L]][[3L]]) == "symbol")),
-    "Specify distinct variables for the treatments and blocks." =
-      (isTRUE(f$variables[[3L]][[2L]] != f$variables[[3L]][[3L]]))
-  )
-  # Pseudo formula for `model.frame`
-  lhs <- f$variables[[2L]]
-  rhs <- c(f$variables[[3L]][[2L]], f$variables[[3L]][[3L]])
-  pf <- formula(paste(lhs, paste(rhs, collapse = " + "), sep = " ~ "))
-  # Extract `model.frame`
-  mf <- match.call()
-  mf <- mf[c(1L, match(c("formula", "data"), names(mf), 0L))]
-  mf$drop.unused.levels <- TRUE
-  mf[[1L]] <- quote(model.frame)
-  mf[[2L]] <- pf
-  mf <- eval(mf, parent.frame())
-  attributes(mf)$terms <- NULL
-  ## Type conversion
-  mf[[1L]] <- as.numeric(mf[[1L]])
-  mf[[2L]] <- as.factor(mf[[2L]])
-  mf[[3L]] <- as.factor(mf[[3L]])
-  if (nlevels(mf[[2L]]) >= nlevels(mf[[3L]])) {
-    stop("The number of blocks must be larger than the number of treatments.")
-  }
-  ## Construct general block design
-  # Incidence matrix
-  c <- unclass(table(mf[[3L]], mf[[2L]]))
-  # Model matrix
-  x <- reshape(mf[order(mf[[2L]]), ],
-    idvar = names(mf)[3L],
-    timevar = names(mf)[2L],
-    v.names = names(mf)[1L],
-    direction = "wide"
-  )
-  x <- x[order(x[[names(mf)[3L]]]), ]
-  # Replace `NA` with `0`
-  x[is.na(x)] <- 0
-  # Remove block variable and convert to matrix
-  x[names(mf)[3L]] <- NULL
-  x <- as.matrix(x)
-  # Name rows and columns
-  dimnames(x) <- list(levels(mf[[3L]]), levels(mf[[2L]]))
-  # General block design
-  gbd <-
-    list("model_matrix" = x, "incidence_matrix" = c, "trt" = levels(mf[[2L]]))
-  class(gbd) <- c("gbd", "melt")
-  # Check whether all pairwise comparisons or comparisons to control
-  match.arg(control, gbd$trt)
-  if (is.null(control)) {
-    ctrl <- 0
-  } else {
-    ctrl <- which(control == gbd$trt)
-  }
-  ## Pairwise comparisons
-  out <- pairwise(gbd$model_matrix, gbd$incidence_matrix,
-    control = ctrl, k, alpha, interval = TRUE,
-    method, B, nthreads, th = 50, maxit, abstol
-  )
-  out$trt <- gbd$trt
-  out$control <- control
-  out$model.matrix <- gbd$model_matrix
-  out$incidence.matrix <- gbd$incidence_matrix
-  if (!all(out$convergence)) {
-    if (method == "NB") {
-      warning(
-        "Convergence failed and switched to AMC for confidence intervals.\n"
-      )
-    } else {
-      warning("Convergence failed.\n")
-    }
-  }
-  out
-}
-
-#' @exportS3Method print pairwise
-print.pairwise <- function(x, ...) {
-  stopifnot(inherits(x, "melt"))
-  cat("\n\tEmpirical Likelihood Multiple Tests\n\n")
-  if (is.null(x$control)) {
-    cat("All pairwise comparisons\n\n")
-    rname <- vector("character", length = 0L)
-    for (i in 1L:(length(x$trt) - 1L)) {
-      for (j in (i + 1L):length(x$trt)) {
-        rname <- c(rname, paste(x$trt[i], "-", x$trt[j]))
-      }
-    }
-  } else {
-    cat("Comparisons with control\n\n")
-    diff <- setdiff(x$trt, x$control)
-    rname <- vector("character", length = length(diff))
-    for (i in seq_along(diff)) {
-      rname[i] <- paste(diff[i], "-", x$control)
-    }
-  }
-  out <- data.frame(
-    row.names = rname, Estimate = x$estimate, Chisq = x$statistic,
-    Lwr.ci = x$lower, Upr.ci = x$upper, p.adj = x$p.adj
-  )
-  printCoefmat(out,
-    digits = min(4L, getOption("digits")),
-    dig.tst = min(3L, getOption("digits")), cs.ind = c(1L, 3L, 4L),
-    tst.ind = 2L, P.values = TRUE, has.Pvalue = TRUE, eps.Pvalue = 1e-03
-  )
-  cat("\n")
-  cat(paste(c("k", "level", "method", "cutoff"),
-    c(x$k, x$level, x$method, round(x$cutoff, 4L)),
-    collapse = ", ", sep = ": "
-  ))
-  cat("\n\n")
-}
diff --git a/deprecated.cpp b/deprecated.cpp
deleted file mode 100644
index cb70e68..0000000
--- a/deprecated.cpp
+++ /dev/null
@@ -1,619 +0,0 @@
-#include "deprecated.h"
-#include <RcppEigen.h>
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-#include <cmath>
-#include <string>
-
-// [[Rcpp::export]]
-Rcpp::List pairwise(const Eigen::MatrixXd &x,
-                    const Eigen::MatrixXd &c,
-                    const int control = 0,
-                    const int k = 1,
-                    const double level = 0.05,
-                    const bool interval = true,
-                    const std::string method = "AMC",
-                    const int B = 1e+04,
-                    const int nthreads = 1,
-                    const double th = 50,
-                    const int maxit = 1e+04,
-                    const double abstol = 1e-08)
-{
-  std::vector<std::array<int, 2>> pairs = comparison_pairs(x.cols(), control);
-  const int m = pairs.size();
-  std::vector<double> estimate(m);
-  std::vector<double> statistic(m);
-  std::vector<bool> convergence(m);
-  const Eigen::VectorXd theta_hat =
-      x.array().colwise().sum() / c.array().colwise().sum();
-  for (int i = 0; i < m; ++i)
-  {
-    Rcpp::checkUserInterrupt();
-    estimate[i] = theta_hat(pairs[i][0]) - theta_hat(pairs[i][1]);
-    Eigen::MatrixXd lhs = Eigen::MatrixXd::Zero(1, x.cols());
-    lhs(pairs[i][0]) = 1;
-    lhs(pairs[i][1]) = -1;
-    minEL pairwise_result =
-        test_gbd_EL(theta_hat, x, c, lhs, Eigen::Matrix<double, 1, 1>(0),
-                    th*100, maxit, abstol);
-    statistic[i] = 2 * pairwise_result.nllr;
-    convergence[i] = pairwise_result.convergence;
-  }
-  bool anyfail = std::any_of(convergence.begin(), convergence.end(),
-                             [](bool v)
-                             { return !v; });
-  Eigen::ArrayXd bootstrap_statistics_pairwise(B);
-  if (method == "AMC")
-  {
-    bootstrap_statistics_pairwise =
-        bootstrap_statistics_pairwise_AMC(x, c, k, pairs, B, level);
-  }
-  else
-  {
-    bootstrap_statistics_pairwise = bootstrap_statistics_pairwise_NB(
-        x, c, k, pairs, B, level, nthreads, th*100, maxit, abstol);
-  }
-  std::vector<double> adj_pvalues(m);
-  for (int i = 0; i < m; ++i)
-  {
-    adj_pvalues[i] =
-        static_cast<double>(
-            (bootstrap_statistics_pairwise >= statistic[i]).count()) / B;
-  }
-  Rcpp::Function quantile("quantile");
-  double cutoff = Rcpp::as<double>(quantile(bootstrap_statistics_pairwise,
-                                            Rcpp::Named("probs") = 1 - level));
-  Rcpp::List result;
-  result["estimate"] = estimate;
-  result["statistic"] = statistic;
-  result["convergence"] = convergence;
-  result["cutoff"] = cutoff;
-  if (method == "NB" && anyfail)
-  {
-    bootstrap_statistics_pairwise =
-        bootstrap_statistics_pairwise_AMC(x, c, k, pairs, B, level);
-    cutoff = Rcpp::as<double>(quantile(bootstrap_statistics_pairwise,
-                                       Rcpp::Named("probs") = 1 - level));
-  }
-  if (interval)
-  {
-    std::vector<double> lower(m);
-    std::vector<double> upper(m);
-    for (int i = 0; i < m; ++i)
-    {
-      Rcpp::checkUserInterrupt();
-      Eigen::MatrixXd lhs = Eigen::MatrixXd::Zero(1, x.cols());
-      lhs(pairs[i][0]) = 1;
-      lhs(pairs[i][1]) = -1;
-      std::array<double, 2> ci = pair_confidence_interval_gbd(
-          theta_hat, x, c, lhs, th*100, estimate[i], cutoff);
-      lower[i] = ci[0];
-      upper[i] = ci[1];
-    }
-    result["lower"] = lower;
-    result["upper"] = upper;
-  }
-  result["p.adj"] = adj_pvalues;
-  result["k"] = k;
-  result["level"] = level;
-  result["method"] = method;
-  result["B"] = bootstrap_statistics_pairwise.size();
-  result.attr("class") = Rcpp::CharacterVector({"pairwise", "melt"});
-  return result;
-}
-
-EL_deprecated::EL_deprecated(const Eigen::Ref<const Eigen::MatrixXd> &g,
-                             const double th, const int maxit,
-                             const double abstol)
-{
-  l = (g.transpose() * g).ldlt().solve(g.colwise().sum());
-  iterations = 0;
-  convergence = false;
-  while (!convergence && iterations != maxit)
-  {
-    PseudoLog_deprecated log_tmp(Eigen::VectorXd::Ones(g.rows()) + g * l);
-    const Eigen::MatrixXd J = g.array().colwise() * log_tmp.sqrt_neg_d2plog;
-    Eigen::VectorXd step =
-        (J.transpose() * J)
-            .ldlt()
-            .solve(J.transpose() *
-                   (log_tmp.dplog / log_tmp.sqrt_neg_d2plog).matrix());
-    nllr = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g.rows()) +
-                                     g * (l + step));
-    while (nllr < log_tmp.plog_sum)
-    {
-      step /= 2;
-      nllr = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g.rows()) +
-                                       g * (l + step));
-    }
-    l += step;
-    if (nllr > th)
-    {
-      break;
-    }
-    if (nllr - log_tmp.plog_sum < abstol)
-    {
-      convergence = true;
-    }
-    else
-    {
-      ++iterations;
-    }
-  }
-}
-
-PseudoLog_deprecated::PseudoLog_deprecated(Eigen::VectorXd &&x)
-{
-  static const double n = static_cast<double>(x.size());
-  static const double a0 = 1.0 / n;
-  static const double a1 = -std::log(n) - 1.5;
-  static const double a2 = 2.0 * n;
-  static const double a3 = -0.5 * n * n;
-
-  dplog.resize(x.size());
-  sqrt_neg_d2plog.resize(x.size());
-  plog_sum = 0;
-
-  for (unsigned int i = 0; i < x.size(); ++i)
-  {
-    if (x[i] < a0)
-    {
-      dplog[i] = a2 + 2 * a3 * x[i];
-      sqrt_neg_d2plog[i] = a2 / 2;
-      plog_sum += a1 + a2 * x[i] + a3 * x[i] * x[i];
-    }
-    else
-    {
-      dplog[i] = 1.0 / x[i];
-      sqrt_neg_d2plog[i] = 1.0 / x[i];
-      plog_sum += std::log(x[i]);
-    }
-  }
-}
-
-double PseudoLog_deprecated::sum(Eigen::VectorXd &&x)
-{
-  static const double n = static_cast<double>(x.size());
-  static const double a0 = 1.0 / n;
-  static const double a1 = -std::log(n) - 1.5;
-  static const double a2 = 2.0 * n;
-  static const double a3 = -0.5 * n * n;
-  double out = 0;
-  for (unsigned int i = 0; i < x.size(); ++i)
-  {
-    out += x[i] < a0 ? a1 + a2 * x[i] + a3 * x[i] * x[i] : std::log(x[i]);
-  }
-  return out;
-}
-
-Eigen::ArrayXd PseudoLog_deprecated::dp(Eigen::VectorXd &&x)
-{
-  static const double n = static_cast<double>(x.size());
-  static const double a0 = 1.0 / n;
-  static const double a1 = 2.0 * n;
-  static const double a2 = -1.0 * n * n;
-  for (unsigned int i = 0; i < x.size(); ++i)
-  {
-    if (x[i] < a0)
-    {
-      x[i] = a1 + a2 * x[i];
-    }
-    else
-    {
-      x[i] = 1.0 / x[i];
-    }
-  }
-  return x;
-}
-
-Eigen::VectorXd linear_projection(const Eigen::Ref<const Eigen::VectorXd> &par,
-                                  const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                                  const Eigen::Ref<const Eigen::VectorXd> &rhs)
-{
-  return par - lhs.transpose() * (lhs * lhs.transpose()).inverse() *
-                   (lhs * par - rhs);
-}
-
-void linear_projection_void(Eigen::Ref<Eigen::VectorXd> par,
-                            const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                            const Eigen::Ref<const Eigen::VectorXd> &rhs)
-{
-  par -=
-      lhs.transpose() * (lhs * lhs.transpose()).inverse() * (lhs * par - rhs);
-}
-
-Eigen::MatrixXd bootstrap_sample(const Eigen::Ref<const Eigen::MatrixXd> &x,
-                                 const Eigen::Ref<const Eigen::ArrayXi> &index)
-{
-  Eigen::MatrixXd out(x.rows(), x.cols());
-  for (int i = 0; i < x.rows(); ++i)
-  {
-    out.row(i) = x.row(index(i));
-  }
-  return out;
-}
-
-Eigen::MatrixXd g_gbd(const Eigen::Ref<const Eigen::VectorXd> &par,
-                      const Eigen::Ref<const Eigen::MatrixXd> &x,
-                      const Eigen::Ref<const Eigen::MatrixXd> &c)
-{
-  return x - (c.array().rowwise() * par.array().transpose()).matrix();
-}
-
-Eigen::MatrixXd s_gbd(const Eigen::Ref<const Eigen::MatrixXd> &x,
-                        const Eigen::Ref<const Eigen::MatrixXd> &c)
-{
-  Eigen::MatrixXd g =
-      g_gbd(x.array().colwise().sum() / c.array().colwise().sum(), x, c);
-  return (g.transpose() * g) / x.rows();
-}
-
-void lambda2theta_void(const Eigen::Ref<const Eigen::VectorXd> &l,
-                       Eigen::Ref<Eigen::VectorXd> par,
-                       const Eigen::Ref<const Eigen::MatrixXd> &g,
-                       const Eigen::Ref<const Eigen::MatrixXd> &c,
-                       const double gamma)
-{
-  par += gamma *
-         ((PseudoLog_deprecated::dp(Eigen::VectorXd::Ones(g.rows()) + g * l)
-               .matrix().asDiagonal() *
-           c).array().colwise().sum().transpose() * l.array()).matrix();
-}
-
-Eigen::MatrixXd rmvn(const Eigen::MatrixXd &x, const int n)
-{
-  Eigen::MatrixXd I(n, x.cols());
-  for (int j = 0; j < x.cols(); ++j)
-  {
-    for (int i = 0; i < n; ++i)
-    {
-      I(i, j) = R::rnorm(0, 1.0);
-    }
-  }
-  const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(x);
-  return I * es.operatorSqrt();
-}
-
-minEL test_gbd_EL(const Eigen::Ref<const Eigen::VectorXd> &par0,
-                  const Eigen::Ref<const Eigen::MatrixXd> &x,
-                  const Eigen::Ref<const Eigen::MatrixXd> &c,
-                  const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                  const Eigen::Ref<const Eigen::VectorXd> &rhs,
-                  const double th,
-                  const int maxit,
-                  const double abstol)
-{
-  Eigen::VectorXd par = linear_projection(par0, lhs, rhs);
-  Eigen::MatrixXd g = g_gbd(par, x, c);
-  Eigen::VectorXd l = EL_deprecated(g, th).l;
-  double f1 =
-      PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g.rows()) + g * l);
-  double gamma = 1.0 / (c.colwise().sum().mean());
-  bool convergence = false;
-  int iterations = 0;
-  while (!convergence && iterations != maxit)
-  {
-    Eigen::VectorXd par_tmp = par;
-    lambda2theta_void(l, par_tmp, g, c, gamma);
-    linear_projection_void(par_tmp, lhs, rhs);
-    Eigen::MatrixXd g_tmp = g_gbd(par_tmp, x, c);
-    EL_deprecated eval(g_tmp, th);
-    Eigen::VectorXd l_tmp = eval.l;
-    if (!eval.convergence && iterations > 9)
-    {
-      return {par, l, f1, iterations, convergence};
-    }
-    double f0 = f1;
-    f1 = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g_tmp.rows()) +
-                                   g_tmp * l_tmp);
-    while (f0 < f1)
-    {
-      gamma /= 2;
-      par_tmp = par;
-      lambda2theta_void(l, par_tmp, g, c, gamma);
-      linear_projection_void(par_tmp, lhs, rhs);
-      g_tmp = g_gbd(par_tmp, x, c);
-      l_tmp = EL_deprecated(g_tmp, th).l;
-      if (gamma < abstol)
-      {
-        return {par, l, f0, iterations, convergence};
-      }
-      f1 = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g_tmp.rows()) +
-                                     g_tmp * l_tmp);
-    }
-    par = std::move(par_tmp);
-    l = std::move(l_tmp);
-    g = std::move(g_tmp);
-    ++iterations;
-    if (f0 - f1 < abstol)
-    {
-      convergence = true;
-    }
-  }
-  return {par, l, f1, iterations, convergence};
-}
-
-double test_nllr(const Eigen::Ref<const Eigen::VectorXd> &par0,
-                 const Eigen::Ref<const Eigen::MatrixXd> &x,
-                 const Eigen::Ref<const Eigen::MatrixXd> &c,
-                 const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                 const Eigen::Ref<const Eigen::VectorXd> &rhs,
-                 const double th,
-                 const int maxit,
-                 const double abstol)
-{
-  Eigen::VectorXd par = linear_projection(par0, lhs, rhs);
-  Eigen::MatrixXd g = g_gbd(par, x, c);
-  Eigen::VectorXd l = EL_deprecated(g, th).l;
-  double f1 =
-      PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g.rows()) + g * l);
-  double gamma = 1.0 / (c.colwise().sum().mean());
-  bool convergence = false;
-  int iterations = 0;
-  while (!convergence && iterations != maxit)
-  {
-    Eigen::VectorXd par_tmp = par;
-    lambda2theta_void(l, par_tmp, g, c, gamma);
-    linear_projection_void(par_tmp, lhs, rhs);
-    Eigen::MatrixXd g_tmp = g_gbd(par_tmp, x, c);
-    EL_deprecated eval(g_tmp, th);
-    Eigen::VectorXd l_tmp = eval.l;
-    if (!eval.convergence && iterations > 9)
-    {
-      return f1;
-    }
-    double f0 = f1;
-    f1 = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g_tmp.rows()) +
-                                   g_tmp * l_tmp);
-    while (f0 < f1)
-    {
-      gamma /= 2;
-      par_tmp = par;
-      lambda2theta_void(l, par_tmp, g, c, gamma);
-      linear_projection_void(par_tmp, lhs, rhs);
-      g_tmp = g_gbd(par_tmp, x, c);
-      l_tmp = EL_deprecated(g_tmp, th).l;
-      if (gamma < abstol)
-      {
-        return f0;
-      }
-      f1 = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g_tmp.rows()) +
-                                     g_tmp * l_tmp);
-    }
-    par = std::move(par_tmp);
-    l = std::move(l_tmp);
-    g = std::move(g_tmp);
-    ++iterations;
-    if (f0 - f1 < abstol)
-    {
-      convergence = true;
-    }
-  }
-  return f1;
-}
-
-double test_nllr(const Eigen::Ref<const Eigen::MatrixXd> &x,
-                 const Eigen::Ref<const Eigen::MatrixXd> &c,
-                 const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                 const Eigen::Ref<const Eigen::VectorXd> &rhs,
-                 const double th,
-                 const int maxit,
-                 const double abstol)
-{
-  Eigen::VectorXd par = linear_projection(
-      x.array().colwise().sum() / c.array().colwise().sum(), lhs, rhs);
-  Eigen::MatrixXd g = g_gbd(par, x, c);
-  Eigen::VectorXd l = EL_deprecated(g, th).l;
-  double f1 =
-      PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g.rows()) + g * l);
-  double gamma = 1.0 / (c.colwise().sum().mean());
-  bool convergence = false;
-  int iterations = 0;
-  while (!convergence && iterations != maxit)
-  {
-    Eigen::VectorXd par_tmp = par;
-    lambda2theta_void(l, par_tmp, g, c, gamma);
-    linear_projection_void(par_tmp, lhs, rhs);
-    Eigen::MatrixXd g_tmp = g_gbd(par_tmp, x, c);
-    EL_deprecated eval(g_tmp, th);
-    Eigen::VectorXd l_tmp = eval.l;
-    if (!eval.convergence && iterations > 9)
-    {
-      return f1;
-    }
-    double f0 = f1;
-    f1 = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g_tmp.rows()) +
-                                   g_tmp * l_tmp);
-    while (f0 < f1)
-    {
-      gamma /= 2;
-      par_tmp = par;
-      lambda2theta_void(l, par_tmp, g, c, gamma);
-      linear_projection_void(par_tmp, lhs, rhs);
-      g_tmp = g_gbd(par_tmp, x, c);
-      l_tmp = EL_deprecated(g_tmp, th).l;
-      if (gamma < abstol)
-      {
-        return f0;
-      }
-      f1 = PseudoLog_deprecated::sum(Eigen::VectorXd::Ones(g_tmp.rows()) +
-                                     g_tmp * l_tmp);
-    }
-    par = std::move(par_tmp);
-    l = std::move(l_tmp);
-    g = std::move(g_tmp);
-    ++iterations;
-    if (f0 - f1 < abstol)
-    {
-      convergence = true;
-    }
-  }
-  return f1;
-}
-
-std::vector<std::array<int, 2>> comparison_pairs(const int p,
-                                                 const int control)
-{
-  std::vector<std::array<int, 2>> pairs;
-  if (control == 0)
-  {
-    pairs.reserve(p * (p - 1) / 2);
-    for (int i = 0; i < p - 1; ++i)
-    {
-      for (int j = i + 1; j < p; ++j)
-      {
-        pairs.emplace_back(std::array<int, 2>{i, j});
-      }
-    }
-  }
-  else
-  {
-    pairs.reserve(p - 1);
-    for (int i = 0; i < p; ++i)
-    {
-      if (i == control - 1)
-      {
-        continue;
-      }
-      pairs.emplace_back(std::array<int, 2>{i, control - 1});
-    }
-  }
-  return pairs;
-}
-
-std::array<double, 2> pair_confidence_interval_gbd(
-    const Eigen::Ref<const Eigen::VectorXd> &par0,
-    const Eigen::Ref<const Eigen::MatrixXd> &x,
-    const Eigen::Ref<const Eigen::MatrixXd> &c,
-    const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-    const double th,
-    const double init,
-    const double cutoff)
-{
-  double upper_lb = init;
-  double upper_size = 1;
-  double upper_ub = init + upper_size;
-  while (2 * test_nllr(par0, x, c, lhs, Eigen::Matrix<double, 1, 1>(upper_ub),
-                       th) <=
-         cutoff)
-  {
-    upper_lb = upper_ub;
-    upper_ub += upper_size;
-  }
-  while (upper_ub - upper_lb > 1e-04)
-  {
-    if (2 * test_nllr(par0, x, c, lhs,
-                      Eigen::Matrix<double, 1, 1>((upper_lb + upper_ub) / 2),
-                      th) >
-        cutoff)
-    {
-      upper_ub = (upper_lb + upper_ub) / 2;
-    }
-    else
-    {
-      upper_lb = (upper_lb + upper_ub) / 2;
-    }
-  }
-  double lower_ub = init;
-  double lower_size = 1;
-  double lower_lb = init - lower_size;
-  while (2 * test_nllr(par0, x, c, lhs, Eigen::Matrix<double, 1, 1>(lower_lb),
-                       th) <=
-         cutoff)
-  {
-    lower_ub = lower_lb;
-    lower_lb -= lower_size;
-  }
-  while (lower_ub - lower_lb > 1e-04)
-  {
-    if (2 * test_nllr(par0, x, c, lhs,
-                      Eigen::Matrix<double, 1, 1>((lower_lb + lower_ub) / 2),
-                      th) >
-        cutoff)
-    {
-      lower_lb = (lower_lb + lower_ub) / 2;
-    }
-    else
-    {
-      lower_ub = (lower_lb + lower_ub) / 2;
-    }
-  }
-  return std::array<double, 2>{lower_ub, upper_lb};
-}
-
-Eigen::ArrayXd bootstrap_statistics_pairwise_AMC(
-    const Eigen::Ref<const Eigen::MatrixXd> &x,
-    const Eigen::Ref<const Eigen::MatrixXd> &c,
-    const int k,
-    const std::vector<std::array<int, 2>> &pairs,
-    const int B,
-    const double level)
-{
-  const Eigen::MatrixXd V_hat = s_gbd(x, c);
-  const Eigen::MatrixXd U_hat = rmvn(V_hat, B);
-  const int m = pairs.size();
-  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
-      bootstrap_statistics(B, m);
-  for (int j = 0; j < m; ++j)
-  {
-    Eigen::RowVectorXd R = Eigen::RowVectorXd::Zero(1, x.cols());
-    R(pairs[j][0]) = 1;
-    R(pairs[j][1]) = -1;
-    Eigen::MatrixXd A_hat = (R.transpose() * R) / (R * V_hat * R.transpose());
-    bootstrap_statistics.col(j) =
-        (U_hat * A_hat * U_hat.transpose()).diagonal();
-  }
-  for (int b = 0; b < B; ++b)
-  {
-    std::sort(bootstrap_statistics.row(b).data(),
-              bootstrap_statistics.row(b).data() + m);
-  }
-  return bootstrap_statistics.col(m - k);
-}
-
-Eigen::ArrayXd bootstrap_statistics_pairwise_NB(
-    const Eigen::Ref<const Eigen::MatrixXd> &x,
-    const Eigen::Ref<const Eigen::MatrixXd> &c,
-    const int k,
-    const std::vector<std::array<int, 2>> &pairs,
-    const int B,
-    const double level,
-    const int nthreads,
-    const double th,
-    const int maxit,
-    const double abstol)
-{
-  const int n = x.rows();
-  const int p = x.cols();
-  const int m = pairs.size();
-  const Eigen::MatrixXd x_centered =
-      x - (c.array().rowwise() *
-           (x.array().colwise().sum() / c.array().colwise().sum()))
-              .matrix();
-  const Eigen::ArrayXXi bootstrap_index =
-      Eigen::Map<Eigen::ArrayXXi, Eigen::Unaligned>(
-          (Rcpp::as<std::vector<int>>(Rcpp::sample(
-               Rcpp::IntegerVector(Rcpp::seq(0, n - 1)), n * B, true)))
-              .data(),
-          n, B);
-  Eigen::ArrayXd k_bootstrap_statistics = Eigen::ArrayXd::Constant(B, NA_REAL);
-  #pragma omp parallel for num_threads(nthreads) schedule(dynamic)
-  for (int b = 0; b < B; ++b)
-  {
-    std::vector<double> bootstrap_statistics(m);
-    for (int j = 0; j < m; ++j)
-    {
-      Eigen::MatrixXd lhs = Eigen::MatrixXd::Zero(1, p);
-      lhs(pairs[j][0]) = 1;
-      lhs(pairs[j][1]) = -1;
-      bootstrap_statistics[j] =
-          2 * test_nllr(bootstrap_sample(x_centered, bootstrap_index.col(b)),
-                        bootstrap_sample(c, bootstrap_index.col(b)), lhs,
-                        Eigen::Matrix<double, 1, 1>(0), th, maxit, abstol);
-    }
-    std::sort(bootstrap_statistics.begin(), bootstrap_statistics.end());
-    k_bootstrap_statistics[b] = bootstrap_statistics[m - k];
-  }
-  return k_bootstrap_statistics;
-}
diff --git a/deprecated.h b/deprecated.h
deleted file mode 100644
index 850c9ea..0000000
--- a/deprecated.h
+++ /dev/null
@@ -1,119 +0,0 @@
-#ifndef EL_deprecated_H_
-#define EL_deprecated_H_
-
-#include <RcppEigen.h>
-#include <array>
-#include <vector>
-
-struct minEL
-{
-        Eigen::VectorXd par;
-        Eigen::VectorXd l;
-        double nllr;
-        int iterations;
-        bool convergence;
-};
-
-class EL_deprecated
-{
-public:
-        Eigen::VectorXd l;
-        double nllr;
-        int iterations;
-        bool convergence;
-
-        EL_deprecated(const Eigen::Ref<const Eigen::MatrixXd> &g,
-                      const double th, const int maxit = 100,
-                      const double abstol = 1e-08);
-};
-
-class PseudoLog_deprecated
-{
-public:
-        Eigen::ArrayXd dplog;
-        Eigen::ArrayXd sqrt_neg_d2plog;
-        double plog_sum;
-
-        PseudoLog_deprecated(Eigen::VectorXd &&x);
-        static double sum(Eigen::VectorXd &&x);
-        static Eigen::ArrayXd dp(Eigen::VectorXd &&x);
-};
-
-Eigen::VectorXd linear_projection(const Eigen::Ref<const Eigen::VectorXd> &par,
-                                  const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                                  const Eigen::Ref<const Eigen::VectorXd> &rhs);
-
-void linear_projection_void(Eigen::Ref<Eigen::VectorXd> par,
-                            const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                            const Eigen::Ref<const Eigen::VectorXd> &rhs);
-
-Eigen::MatrixXd bootstrap_sample(const Eigen::Ref<const Eigen::MatrixXd> &x,
-                                 const Eigen::Ref<const Eigen::ArrayXi> &index);
-
-Eigen::MatrixXd g_gbd(const Eigen::Ref<const Eigen::VectorXd> &par,
-                      const Eigen::Ref<const Eigen::MatrixXd> &x,
-                      const Eigen::Ref<const Eigen::MatrixXd> &c);
-
-Eigen::MatrixXd s_gbd(const Eigen::Ref<const Eigen::MatrixXd> &x,
-                      const Eigen::Ref<const Eigen::MatrixXd> &c);
-
-Eigen::MatrixXd rmvn(const Eigen::MatrixXd &x, const int n);
-
-minEL test_gbd_EL(const Eigen::Ref<const Eigen::VectorXd> &par0,
-                  const Eigen::Ref<const Eigen::MatrixXd> &x,
-                  const Eigen::Ref<const Eigen::MatrixXd> &c,
-                  const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                  const Eigen::Ref<const Eigen::VectorXd> &rhs,
-                  const double th,
-                  const int maxit = 1000,
-                  const double abstol = 1e-08);
-
-double test_nllr(const Eigen::Ref<const Eigen::VectorXd> &par0,
-                 const Eigen::Ref<const Eigen::MatrixXd> &x,
-                 const Eigen::Ref<const Eigen::MatrixXd> &c,
-                 const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                 const Eigen::Ref<const Eigen::VectorXd> &rhs,
-                 const double th,
-                 const int maxit = 1000,
-                 const double abstol = 1e-08);
-
-double test_nllr(const Eigen::Ref<const Eigen::MatrixXd> &x,
-                 const Eigen::Ref<const Eigen::MatrixXd> &c,
-                 const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                 const Eigen::Ref<const Eigen::VectorXd> &rhs,
-                 const double th,
-                 const int maxit = 1000,
-                 const double abstol = 1e-08);
-
-std::vector<std::array<int, 2>> comparison_pairs(const int p,
-                                                 const int control);
-
-std::array<double, 2> pair_confidence_interval_gbd(
-                const Eigen::Ref<const Eigen::VectorXd> &par0,
-                const Eigen::Ref<const Eigen::MatrixXd> &x,
-                const Eigen::Ref<const Eigen::MatrixXd> &c,
-                const Eigen::Ref<const Eigen::MatrixXd> &lhs,
-                const double th,
-                const double init,
-                const double cutoff);
-
-Eigen::ArrayXd bootstrap_statistics_pairwise_AMC(
-                const Eigen::Ref<const Eigen::MatrixXd> &x,
-                const Eigen::Ref<const Eigen::MatrixXd> &c,
-                const int k,
-                const std::vector<std::array<int, 2>> &pairs,
-                const int B,
-                const double level);
-
-Eigen::ArrayXd bootstrap_statistics_pairwise_NB(
-    const Eigen::Ref<const Eigen::MatrixXd> &x,
-    const Eigen::Ref<const Eigen::MatrixXd> &c,
-    const int k,
-    const std::vector<std::array<int, 2>> &pairs,
-    const int B,
-    const double level,
-    const int nthreads,
-    const double th,
-    const int maxit,
-    const double abstol);
-#endif
diff --git a/deprecated.o b/deprecated.o
deleted file mode 100644
index 5c72345..0000000
Binary files a/deprecated.o and /dev/null differ