diff --git a/R/apes.R b/R/apes.R
index 40b7cc6..8e5f665 100644
--- a/R/apes.R
+++ b/R/apes.R
@@ -166,13 +166,13 @@ apes <- function(
 
   # Center regressor matrix (if required)
   if (control[["keep_mx"]]) {
-    MX <- object[["MX"]]
+    mx <- object[["mx"]]
   } else {
-    MX <- center_variables_r_(x, w, k_list, control[["center_tol"]], 10000L)
+    mx <- center_variables_r_(x, w, k_list, control[["center_tol"]], 10000L)
   }
 
   # Compute average partial effects, derivatives, and Jacobian
-  PX <- x - MX
+  PX <- x - mx
   Delta <- matrix(NA_real_, nt, p)
   Delta1 <- matrix(NA_real_, nt, p)
   J <- matrix(NA_real_, p, p)
@@ -193,7 +193,7 @@ apes <- function(
     } else {
       Delta[, j] <- beta[[j]] * Delta[, j]
       Delta1[, j] <- beta[[j]] * Delta1[, j]
-      J[, j] <- colSums(MX * Delta1[, j]) / nt_full
+      J[, j] <- colSums(mx * Delta1[, j]) / nt_full
       J[j, j] <- sum(mu.eta) / nt_full + J[j, j]
     }
   }
@@ -209,18 +209,22 @@ apes <- function(
 
   # Compute analytical bias correction of average partial effects
   if (bias_corr) {
-    b <- apes_bias_correction_(eta, family, x, beta, binary, nt, p, PPsi, z,
-      w, k_list, panel_structure, L, k, MPsi, v)
+    b <- apes_bias_correction_(
+      eta, family, x, beta, binary, nt, p, PPsi, z,
+      w, k_list, panel_structure, L, k, MPsi, v
+    )
     delta <- delta - b
   }
   rm(eta, w, z, MPsi)
 
   # Compute covariance matrix
-  Gamma <- gamma_(MX, object[["hessian"]], J, PPsi, v, nt_full)
+  Gamma <- gamma_(mx, object[["hessian"]], J, PPsi, v, nt_full)
   V <- crossprod(Gamma)
 
-  V <- apes_adjust_covariance_(V, Delta, Gamma, k_list, adj, sampling_fe,
-    weak_exo, panel_structure)
+  V <- apes_adjust_covariance_(
+    V, Delta, Gamma, k_list, adj, sampling_fe,
+    weak_exo, panel_structure
+  )
 
   # Add names
   names(delta) <- nms.sp
@@ -268,8 +272,9 @@ apes_set_adj_ <- function(n_pop, nt_full) {
   return(adj)
 }
 
-apes_adjust_covariance_ <- function(V, Delta, Gamma, k_list, adj, sampling_fe,
-  weak_exo, panel_structure) {
+apes_adjust_covariance_ <- function(
+    V, Delta, Gamma, k_list, adj, sampling_fe,
+    weak_exo, panel_structure) {
   if (adj > 0.0) {
     # Simplify covariance if sampling assumptions are imposed
     if (sampling_fe == "independence") {
@@ -298,8 +303,9 @@ apes_adjust_covariance_ <- function(V, Delta, Gamma, k_list, adj, sampling_fe,
   return(V)
 }
 
-apes_bias_correction_ <- function(eta, family, x, beta, binary, nt, p, PPsi,
-  z, w, k_list, panel_structure, L, k, MPsi, v) {
+apes_bias_correction_ <- function(
+    eta, family, x, beta, binary, nt, p, PPsi,
+    z, w, k_list, panel_structure, L, k, MPsi, v) {
   # Compute second-order partial derivatives
   Delta2 <- matrix(NA_real_, nt, p)
   Delta2[, !binary] <- partial_mu_eta_(eta, family, 3L)
diff --git a/R/feglm.R b/R/feglm.R
index aba25cf..3d4a2e6 100644
--- a/R/feglm.R
+++ b/R/feglm.R
@@ -178,10 +178,10 @@ feglm <- function(
   x <- NULL
   eta <- NULL
 
-  # Add names to beta, hessian, and MX (if provided) ----
+  # Add names to beta, hessian, and mx (if provided) ----
   names(fit[["coefficients"]]) <- nms_sp
   if (control[["keep_mx"]]) {
-    colnames(fit[["MX"]]) <- nms_sp
+    colnames(fit[["mx"]]) <- nms_sp
   }
   dimnames(fit[["hessian"]]) <- list(nms_sp, nms_sp)
 
diff --git a/R/feglm_helpers.R b/R/feglm_helpers.R
index b1e4a63..66204bc 100644
--- a/R/feglm_helpers.R
+++ b/R/feglm_helpers.R
@@ -468,7 +468,7 @@ get_score_matrix_ <- function(object) {
 
   # Center regressor matrix (if required)
   if (control[["keep_mx"]]) {
-    MX <- object[["MX"]]
+    mx <- object[["mx"]]
   } else {
     # Extract additional required quantities from result list
     formula <- object[["formula"]]
@@ -483,24 +483,24 @@ get_score_matrix_ <- function(object) {
     attr(x, "dimnames") <- NULL
 
     # Center variables
-    MX <- center_variables_r_(x, w, k.list, control[["center_tol"]], 10000L)
-    colnames(MX) <- nms_sp
+    mx <- center_variables_r_(x, w, k.list, control[["center_tol"]], 10000L)
+    colnames(mx) <- nms_sp
   }
 
   # Return score matrix
-  MX * (nu * w)
+  mx * (nu * w)
 }
 
 #' @title Gamma computation
 #' @description Computes the gamma matrix for the APES function
-#' @param MX Regressor matrix
+#' @param mx Regressor matrix
 #' @param H Hessian matrix
 #' @param J Jacobian matrix
 #' @param PPsi Psi matrix
 #' @param v Vector of weights
 #' @param nt Number of observations
 #' @noRd
-gamma_ <- function(MX, H, J, PPsi, v, nt) {
+gamma_ <- function(mx, H, J, PPsi, v, nt) {
   inv_nt <- 1.0 / nt
-  (MX %*% solve(H * inv_nt, J) - PPsi) * v * inv_nt
+  (mx %*% solve(H * inv_nt, J) - PPsi) * v * inv_nt
 }
diff --git a/R/fenegbin.R b/R/fenegbin.R
index 57379a4..2fc6389 100644
--- a/R/fenegbin.R
+++ b/R/fenegbin.R
@@ -201,15 +201,17 @@ fenegbin <- function(
   # Information if convergence failed ----
   if (!conv && trace) cat("Algorithm did not converge.\n")
 
-  # Add names to beta, hessian, and MX (if provided) ----
+  # Add names to beta, hessian, and mx (if provided) ----
   names(fit[["coefficients"]]) <- nms_sp
   if (control[["keep_mx"]]) {
-    colnames(fit[["MX"]]) <- nms_sp
+    colnames(fit[["mx"]]) <- nms_sp
   }
   dimnames(fit[["hessian"]]) <- list(nms_sp, nms_sp)
 
-  fenegbin_result_list_(fit, theta, iter, conv, nobs, lvls_k, nms_fe,
-    formula, data, family, control)
+  fenegbin_result_list_(
+    fit, theta, iter, conv, nobs, lvls_k, nms_fe,
+    formula, data, family, control
+  )
 }
 
 # Convergence Check ----
@@ -222,8 +224,9 @@ fenegbin_check_convergence_ <- function(dev, dev_old, theta, theta_old, tol) {
 
 # Generate result list ----
 
-fenegbin_result_list_ <- function(fit, theta, iter, conv, nobs, lvls_k,
-  nms_fe, formula, data, family, control) {
+fenegbin_result_list_ <- function(
+    fit, theta, iter, conv, nobs, lvls_k,
+    nms_fe, formula, data, family, control) {
   reslist <- c(
     fit, list(
       theta      = theta,
diff --git a/R/generics_vcov.R b/R/generics_vcov.R
index 5eac8a6..6b306c2 100644
--- a/R/generics_vcov.R
+++ b/R/generics_vcov.R
@@ -69,7 +69,7 @@ vcov.feglm <- function(
   # it is totally fine not to have a cluster variable
   cl_vars <- vcov_feglm_vars_(object)
   k <- length(cl_vars)
-  
+
   if (isTRUE(k >= 1L) && type != "clustered") {
     type <- "clustered"
   }
@@ -77,7 +77,7 @@ vcov.feglm <- function(
   # Compute requested type of covariance matrix
   H <- object[["hessian"]]
   p <- ncol(H)
-  
+
   if (type == "hessian") {
     # If the hessian is invertible, compute its inverse
     V <- vcov_feglm_hessian_covariance_(H, p)
@@ -87,8 +87,10 @@ vcov.feglm <- function(
       # Check if the OP is invertible and compute its inverse
       V <- vcov_feglm_outer_covariance_(G, p)
     } else {
-      V <- vcov_feglm_clustered_sandwich_covariance_(object, type, H, G,
-        cl_vars, sp_vars, k, p)
+      V <- vcov_feglm_clustered_sandwich_covariance_(
+        object, type, H, G,
+        cl_vars, sp_vars, k, p
+      )
     }
   }
 
@@ -117,8 +119,9 @@ vcov_feglm_outer_covariance_ <- function(G, p) {
   V
 }
 
-vcov_feglm_clustered_sandwich_covariance_ <- function(object, type, H, G,
-  cl_vars, sp_vars, k, p) {
+vcov_feglm_clustered_sandwich_covariance_ <- function(
+    object, type, H, G,
+    cl_vars, sp_vars, k, p) {
   # Check if the hessian is invertible and compute its inverse
   V <- try(solve(H), silent = TRUE)
   if (inherits(V, "try-error")) {
@@ -198,7 +201,7 @@ vcov_feglm_clustered_covariance_ <- function(G, cl_vars, sp_vars, p) {
         )
       )
     }
-    
+
     # Update outer product
     if (i %% 2L) {
       B <- B + B_r
diff --git a/tests/testthat/test-kendall.R b/tests/testthat/test-kendall.R
deleted file mode 100644
index 1cf1077..0000000
--- a/tests/testthat/test-kendall.R
+++ /dev/null
@@ -1,79 +0,0 @@
-test_that("kendall", {
-  x <- 1:2
-  expect_equal(kendall_cor(x, x), cor(x, x, method = "kendall"))
-
-  x <- 1:3
-  expect_equal(kendall_cor(x, x), 1)
-
-  x <- rep(1, 3)
-  expect_warning(kendall_cor(x, x), "zero variance")
-
-  x <- c(1, 0, 2)
-  y <- c(5, 3, 4)
-  expect_equal(kendall_cor(x, y), cor(x, y, method = "kendall"))
-
-  k1 <- kendall_cor_test(x, y, alternative = "two.sided")
-  k2 <- cor.test(x, y, method = "kendall", alternative = "two.sided")
-  expect_equal(k1$statistic, unname(k2$estimate))
-  expect_equal(k1$p_value, k2$p.value)
-
-  x <- 1:3
-  y <- 3:1
-  expect_equal(kendall_cor(x, y), cor(x, y, method = "kendall"))
-
-  x <- c(1, NA, 2)
-  y <- 3:1
-  expect_equal(kendall_cor(x, y), cor(x, y, method = "kendall", use = "pairwise.complete.obs"))
-
-  set.seed(123)
-  x <- rnorm(100)
-  y <- rpois(100, 2)
-  expect_equal(kendall_cor(x, y), cor(x, y, method = "kendall"))
-
-  k1 <- kendall_cor_test(x, y, alternative = "two.sided")
-  k2 <- cor.test(x, y, method = "kendall", alternative = "two.sided")
-  expect_equal(k1$statistic, unname(k2$estimate))
-  expect_equal(k1$p_value, k2$p.value)
-
-  k1 <- kendall_cor_test(x, y, alternative = "greater")
-  k2 <- cor.test(x, y, method = "kendall", alternative = "greater")
-  expect_equal(k1$statistic, unname(k2$estimate))
-  expect_equal(k1$p_value, k2$p.value)
-
-  k1 <- kendall_cor_test(x, y, alternative = "less")
-  k2 <- cor.test(x, y, method = "kendall", alternative = "less")
-  expect_equal(k1$statistic, unname(k2$estimate))
-  expect_equal(k1$p_value, k2$p.value)
-
-  x <- rnorm(1e3)
-  y <- rpois(1e3, 2)
-  expect_equal(kendall_cor(x, y), cor(x, y, method = "kendall"))
-
-  t_kendall <- c()
-  for (i in 1:100) {
-    t1 <- Sys.time()
-    kendall_cor(x, y)
-    t2 <- Sys.time()
-    t_kendall[i] <- t2 - t1
-  }
-  t_kendall <- median(t_kendall)
-
-  t_cor <- c()
-  for (i in 1:100) {
-    t1 <- Sys.time()
-    cor(x, y, method = "kendall")
-    t2 <- Sys.time()
-    t_cor[i] <- t2 - t1
-  }
-  t_cor <- median(t_cor)
-
-  expect_lt(t_kendall, t_cor)
-
-  x <- 1:3
-  y <- NA
-  expect_error(kendall_cor(x, y))
-
-  x <- 1:3
-  y <- c(1, NA, NA)
-  expect_error(kendall_cor(x, y))
-})