diff --git a/R-package/.Rbuildignore b/R-package/.Rbuildignore
index b37d627ba487..b1932e324589 100644
--- a/R-package/.Rbuildignore
+++ b/R-package/.Rbuildignore
@@ -4,3 +4,5 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 README.md
+^doc$
+^Meta$
diff --git a/R-package/DESCRIPTION b/R-package/DESCRIPTION
index 19b55c40181e..26a90d3254e5 100644
--- a/R-package/DESCRIPTION
+++ b/R-package/DESCRIPTION
@@ -2,7 +2,7 @@ Package: xgboost
 Type: Package
 Title: Extreme Gradient Boosting
 Version: 1.7.6.1
-Date: 2023-06-16
+Date: 2023-12-05
 Authors@R: c(
   person("Tianqi", "Chen", role = c("aut"),
          email = "tianqi.tchen@gmail.com"),
diff --git a/R-package/R/callbacks.R b/R-package/R/callbacks.R
index fa947346937d..e4690d1e73df 100644
--- a/R-package/R/callbacks.R
+++ b/R-package/R/callbacks.R
@@ -70,7 +70,7 @@ cb.print.evaluation <- function(period = 1, showsd = TRUE) {
         i == env$begin_iteration ||
         i == env$end_iteration) {
       stdev <- if (showsd) env$bst_evaluation_err else NULL
-      msg <- format.eval.string(i, env$bst_evaluation, stdev)
+      msg <- .format_eval_string(i, env$bst_evaluation, stdev)
       cat(msg, '\n')
     }
   }
@@ -380,7 +380,9 @@ cb.early.stop <- function(stopping_rounds, maximize = FALSE,
     if ((maximize && score > best_score) ||
         (!maximize && score < best_score)) {
 
-      best_msg <<- format.eval.string(i, env$bst_evaluation, env$bst_evaluation_err)
+      best_msg <<- .format_eval_string(
+        i, env$bst_evaluation, env$bst_evaluation_err
+      )
       best_score <<- score
       best_iteration <<- i
       best_ntreelimit <<- best_iteration * env$num_parallel_tree
@@ -555,14 +557,18 @@ cb.cv.predict <- function(save_models = FALSE) {
 #'
 #' @examples
 #' #### Binary classification:
-#' #
+#'
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
+#'
 #' # In the iris dataset, it is hard to linearly separate Versicolor class from the rest
 #' # without considering the 2nd order interactions:
 #' x <- model.matrix(Species ~ .^2, iris)[,-1]
 #' colnames(x)
-#' dtrain <- xgb.DMatrix(scale(x), label = 1*(iris$Species == "versicolor"), nthread = 2)
+#' dtrain <- xgb.DMatrix(scale(x), label = 1*(iris$Species == "versicolor"), nthread = nthread)
 #' param <- list(booster = "gblinear", objective = "reg:logistic", eval_metric = "auc",
-#'               lambda = 0.0003, alpha = 0.0003, nthread = 2)
+#'               lambda = 0.0003, alpha = 0.0003, nthread = nthread)
 #' # For 'shotgun', which is a default linear updater, using high eta values may result in
 #' # unstable behaviour in some datasets. With this simple dataset, however, the high learning
 #' # rate does not break the convergence, but allows us to illustrate the typical pattern of
@@ -592,9 +598,9 @@ cb.cv.predict <- function(save_models = FALSE) {
 #'
 #' #### Multiclass classification:
 #' #
-#' dtrain <- xgb.DMatrix(scale(x), label = as.numeric(iris$Species) - 1, nthread = 2)
+#' dtrain <- xgb.DMatrix(scale(x), label = as.numeric(iris$Species) - 1, nthread = nthread)
 #' param <- list(booster = "gblinear", objective = "multi:softprob", num_class = 3,
-#'               lambda = 0.0003, alpha = 0.0003, nthread = 2)
+#'               lambda = 0.0003, alpha = 0.0003, nthread = nthread)
 #' # For the default linear updater 'shotgun' it sometimes is helpful
 #' # to use smaller eta to reduce instability
 #' bst <- xgb.train(param, dtrain, list(tr=dtrain), nrounds = 70, eta = 0.5,
@@ -752,7 +758,7 @@ xgb.gblinear.history <- function(model, class_index = NULL) {
 #
 
 # Format the evaluation metric string
-format.eval.string <- function(iter, eval_res, eval_err = NULL) {
+.format_eval_string <- function(iter, eval_res, eval_err = NULL) {
   if (length(eval_res) == 0)
     stop('no evaluation results')
   enames <- names(eval_res)
diff --git a/R-package/R/xgb.Booster.R b/R-package/R/xgb.Booster.R
index f31c506466f9..148863a66b61 100644
--- a/R-package/R/xgb.Booster.R
+++ b/R-package/R/xgb.Booster.R
@@ -259,11 +259,16 @@ xgb.Booster.complete <- function(object, saveraw = TRUE) {
 #'
 #' data(agaricus.train, package='xgboost')
 #' data(agaricus.test, package='xgboost')
+#'
+#' ## Keep the number of threads to 2 for examples
+#' nthread <- 2
+#' data.table::setDTthreads(nthread)
+#'
 #' train <- agaricus.train
 #' test <- agaricus.test
 #'
 #' bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-#'                eta = 0.5, nthread = 2, nrounds = 5, objective = "binary:logistic")
+#'                eta = 0.5, nthread = nthread, nrounds = 5, objective = "binary:logistic")
 #' # use all trees by default
 #' pred <- predict(bst, test$data)
 #' # use only the 1st tree
@@ -329,8 +334,14 @@ predict.xgb.Booster <- function(object, newdata, missing = NA, outputmargin = FA
                                 reshape = FALSE, training = FALSE, iterationrange = NULL, strict_shape = FALSE, ...) {
   object <- xgb.Booster.complete(object, saveraw = FALSE)
 
-  if (!inherits(newdata, "xgb.DMatrix"))
-    newdata <- xgb.DMatrix(newdata, missing = missing, nthread = NVL(object$params[["nthread"]], -1))
+  if (!inherits(newdata, "xgb.DMatrix")) {
+    config <- jsonlite::fromJSON(xgb.config(object))
+    nthread <- strtoi(config$learner$generic_param$nthread)
+    newdata <- xgb.DMatrix(
+      newdata,
+      missing = missing, nthread = NVL(nthread, -1)
+    )
+  }
   if (!is.null(object[["feature_names"]]) &&
       !is.null(colnames(newdata)) &&
       !identical(object[["feature_names"]], colnames(newdata)))
@@ -620,10 +631,15 @@ xgb.attributes <- function(object) {
 #'
 #' @examples
 #' data(agaricus.train, package='xgboost')
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
 #' train <- agaricus.train
 #'
-#' bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-#'                eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic")
+#' bst <- xgboost(
+#'   data = train$data, label = train$label, max_depth = 2,
+#'   eta = 1, nthread = nthread, nrounds = 2, objective = "binary:logistic"
+#' )
 #' config <- xgb.config(bst)
 #'
 #' @rdname xgb.config
diff --git a/R-package/R/xgb.DMatrix.R b/R-package/R/xgb.DMatrix.R
index b843d651833a..76c4029397db 100644
--- a/R-package/R/xgb.DMatrix.R
+++ b/R-package/R/xgb.DMatrix.R
@@ -18,7 +18,12 @@
 #'
 #' @examples
 #' data(agaricus.train, package='xgboost')
-#' dtrain <- with(agaricus.train, xgb.DMatrix(data, label = label, nthread = 2))
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
+#' dtrain <- with(
+#'   agaricus.train, xgb.DMatrix(data, label = label, nthread = nthread)
+#' )
 #' xgb.DMatrix.save(dtrain, 'xgb.DMatrix.data')
 #' dtrain <- xgb.DMatrix('xgb.DMatrix.data')
 #' if (file.exists('xgb.DMatrix.data')) file.remove('xgb.DMatrix.data')
diff --git a/R-package/R/xgb.load.R b/R-package/R/xgb.load.R
index d98041908eb4..6916842ec8d8 100644
--- a/R-package/R/xgb.load.R
+++ b/R-package/R/xgb.load.R
@@ -22,14 +22,23 @@
 #' @examples
 #' data(agaricus.train, package='xgboost')
 #' data(agaricus.test, package='xgboost')
+#'
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
+#'
 #' train <- agaricus.train
 #' test <- agaricus.test
-#' bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-#'                eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+#' bst <- xgboost(
+#'   data = train$data, label = train$label, max_depth = 2, eta = 1,
+#'   nthread = nthread,
+#'   nrounds = 2,
+#'   objective = "binary:logistic"
+#' )
+#'
 #' xgb.save(bst, 'xgb.model')
 #' bst <- xgb.load('xgb.model')
 #' if (file.exists('xgb.model')) file.remove('xgb.model')
-#' pred <- predict(bst, test$data)
 #' @export
 xgb.load <- function(modelfile) {
   if (is.null(modelfile))
diff --git a/R-package/R/xgb.model.dt.tree.R b/R-package/R/xgb.model.dt.tree.R
index 5411c35d2394..33e99dd0effd 100644
--- a/R-package/R/xgb.model.dt.tree.R
+++ b/R-package/R/xgb.model.dt.tree.R
@@ -46,9 +46,12 @@
 #' # Basic use:
 #'
 #' data(agaricus.train, package='xgboost')
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
 #'
 #' bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 2,
-#'                eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+#'                eta = 1, nthread = nthread, nrounds = 2,objective = "binary:logistic")
 #'
 #' (dt <- xgb.model.dt.tree(colnames(agaricus.train$data), bst))
 #'
diff --git a/R-package/R/xgb.plot.deepness.R b/R-package/R/xgb.plot.deepness.R
index 6579fb5117ab..6a182f02ff96 100644
--- a/R-package/R/xgb.plot.deepness.R
+++ b/R-package/R/xgb.plot.deepness.R
@@ -45,10 +45,13 @@
 #' @examples
 #'
 #' data(agaricus.train, package='xgboost')
+#' ## Keep the number of threads to 2 for examples
+#' nthread <- 2
+#' data.table::setDTthreads(nthread)
 #'
-#' # Change max_depth to a higher number to get a more significant result
+#' ## Change max_depth to a higher number to get a more significant result
 #' bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 6,
-#'                eta = 0.1, nthread = 2, nrounds = 50, objective = "binary:logistic",
+#'                eta = 0.1, nthread = nthread, nrounds = 50, objective = "binary:logistic",
 #'                subsample = 0.5, min_child_weight = 2)
 #'
 #' xgb.plot.deepness(bst)
diff --git a/R-package/R/xgb.plot.importance.R b/R-package/R/xgb.plot.importance.R
index ab0951463595..8bd61086270a 100644
--- a/R-package/R/xgb.plot.importance.R
+++ b/R-package/R/xgb.plot.importance.R
@@ -45,9 +45,14 @@
 #'
 #' @examples
 #' data(agaricus.train)
+#' ## Keep the number of threads to 2 for examples
+#' nthread <- 2
+#' data.table::setDTthreads(nthread)
 #'
-#' bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 3,
-#'                eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic")
+#' bst <- xgboost(
+#'   data = agaricus.train$data, label = agaricus.train$label, max_depth = 3,
+#'   eta = 1, nthread = nthread, nrounds = 2, objective = "binary:logistic"
+#' )
 #'
 #' importance_matrix <- xgb.importance(colnames(agaricus.train$data), model = bst)
 #'
diff --git a/R-package/R/xgb.plot.multi.trees.R b/R-package/R/xgb.plot.multi.trees.R
index fc54d777d2b5..cc76d4ecc19a 100644
--- a/R-package/R/xgb.plot.multi.trees.R
+++ b/R-package/R/xgb.plot.multi.trees.R
@@ -43,10 +43,15 @@
 #' @examples
 #'
 #' data(agaricus.train, package='xgboost')
-#'
-#' bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 15,
-#'                eta = 1, nthread = 2, nrounds = 30, objective = "binary:logistic",
-#'                min_child_weight = 50, verbose = 0)
+#' ## Keep the number of threads to 2 for examples
+#' nthread <- 2
+#' data.table::setDTthreads(nthread)
+#'
+#' bst <- xgboost(
+#'   data = agaricus.train$data, label = agaricus.train$label, max_depth = 15,
+#'   eta = 1, nthread = nthread, nrounds = 30, objective = "binary:logistic",
+#'   min_child_weight = 50, verbose = 0
+#' )
 #'
 #' p <- xgb.plot.multi.trees(model = bst, features_keep = 3)
 #' print(p)
diff --git a/R-package/R/xgb.plot.shap.R b/R-package/R/xgb.plot.shap.R
index 8f8e921a86b0..b7e3dca06928 100644
--- a/R-package/R/xgb.plot.shap.R
+++ b/R-package/R/xgb.plot.shap.R
@@ -74,9 +74,14 @@
 #' data(agaricus.train, package='xgboost')
 #' data(agaricus.test, package='xgboost')
 #'
-#' bst <- xgboost(agaricus.train$data, agaricus.train$label, nrounds = 50,
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
+#' nrounds <- 20
+#'
+#' bst <- xgboost(agaricus.train$data, agaricus.train$label, nrounds = nrounds,
 #'                eta = 0.1, max_depth = 3, subsample = .5,
-#'                method = "hist", objective = "binary:logistic", nthread = 2, verbose = 0)
+#'                method = "hist", objective = "binary:logistic", nthread = nthread, verbose = 0)
 #'
 #' xgb.plot.shap(agaricus.test$data, model = bst, features = "odor=none")
 #' contr <- predict(bst, agaricus.test$data, predcontrib = TRUE)
@@ -85,12 +90,11 @@
 #'
 #' # multiclass example - plots for each class separately:
 #' nclass <- 3
-#' nrounds <- 20
 #' x <- as.matrix(iris[, -5])
 #' set.seed(123)
 #' is.na(x[sample(nrow(x) * 4, 30)]) <- TRUE # introduce some missing values
 #' mbst <- xgboost(data = x, label = as.numeric(iris$Species) - 1, nrounds = nrounds,
-#'                 max_depth = 2, eta = 0.3, subsample = .5, nthread = 2,
+#'                 max_depth = 2, eta = 0.3, subsample = .5, nthread = nthread,
 #'                 objective = "multi:softprob", num_class = nclass, verbose = 0)
 #' trees0 <- seq(from=0, by=nclass, length.out=nrounds)
 #' col <- rgb(0, 0, 1, 0.5)
@@ -193,7 +197,7 @@ xgb.plot.shap <- function(data, shap_contrib = NULL, features = NULL, top_n = 1,
 #' hence allows us to see which features have a negative / positive contribution
 #' on the model prediction, and whether the contribution is different for larger
 #' or smaller values of the feature. We effectively try to replicate the
-#' \code{summary_plot} function from https://github.com/slundberg/shap.
+#' \code{summary_plot} function from https://github.com/shap/shap
 #'
 #' @inheritParams xgb.plot.shap
 #'
@@ -202,7 +206,7 @@ xgb.plot.shap <- function(data, shap_contrib = NULL, features = NULL, top_n = 1,
 #'
 #' @examples # See \code{\link{xgb.plot.shap}}.
 #' @seealso \code{\link{xgb.plot.shap}}, \code{\link{xgb.ggplot.shap.summary}},
-#'   \url{https://github.com/slundberg/shap}
+#'   \url{https://github.com/shap/shap}
 xgb.plot.shap.summary <- function(data, shap_contrib = NULL, features = NULL, top_n = 10, model = NULL,
                                   trees = NULL, target_class = NULL, approxcontrib = FALSE, subsample = NULL) {
   # Only ggplot implementation is available.
diff --git a/R-package/R/xgb.save.R b/R-package/R/xgb.save.R
index 42ecb4153a84..296dca48a201 100644
--- a/R-package/R/xgb.save.R
+++ b/R-package/R/xgb.save.R
@@ -25,14 +25,22 @@
 #' @examples
 #' data(agaricus.train, package='xgboost')
 #' data(agaricus.test, package='xgboost')
+#'
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
+#'
 #' train <- agaricus.train
 #' test <- agaricus.test
-#' bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-#'                eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+#' bst <- xgboost(
+#'   data = train$data, label = train$label, max_depth = 2, eta = 1,
+#'   nthread = nthread,
+#'   nrounds = 2,
+#'   objective = "binary:logistic"
+#' )
 #' xgb.save(bst, 'xgb.model')
 #' bst <- xgb.load('xgb.model')
 #' if (file.exists('xgb.model')) file.remove('xgb.model')
-#' pred <- predict(bst, test$data)
 #' @export
 xgb.save <- function(model, fname) {
   if (typeof(fname) != "character")
diff --git a/R-package/R/xgb.save.raw.R b/R-package/R/xgb.save.raw.R
index 48fdbca45b13..cad0fb0e01c2 100644
--- a/R-package/R/xgb.save.raw.R
+++ b/R-package/R/xgb.save.raw.R
@@ -16,13 +16,18 @@
 #' @examples
 #' data(agaricus.train, package='xgboost')
 #' data(agaricus.test, package='xgboost')
+#'
+#' ## Keep the number of threads to 2 for examples
+#' nthread <- 2
+#' data.table::setDTthreads(nthread)
+#'
 #' train <- agaricus.train
 #' test <- agaricus.test
 #' bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-#'                eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+#'                eta = 1, nthread = nthread, nrounds = 2,objective = "binary:logistic")
+#'
 #' raw <- xgb.save.raw(bst)
 #' bst <- xgb.load.raw(raw)
-#' pred <- predict(bst, test$data)
 #'
 #' @export
 xgb.save.raw <- function(model, raw_format = "deprecated") {
diff --git a/R-package/R/xgb.train.R b/R-package/R/xgb.train.R
index f23700511684..50edf6db170d 100644
--- a/R-package/R/xgb.train.R
+++ b/R-package/R/xgb.train.R
@@ -124,7 +124,8 @@
 #' than the \code{xgboost} interface.
 #'
 #' Parallelization is automatically enabled if \code{OpenMP} is present.
-#' Number of threads can also be manually specified via \code{nthread} parameter.
+#' Number of threads can also be manually specified via the \code{nthread}
+#' parameter.
 #'
 #' The evaluation metric is chosen automatically by XGBoost (according to the objective)
 #' when the \code{eval_metric} parameter is not provided.
@@ -192,17 +193,25 @@
 #' data(agaricus.train, package='xgboost')
 #' data(agaricus.test, package='xgboost')
 #'
-#' dtrain <- with(agaricus.train, xgb.DMatrix(data, label = label, nthread = 2))
-#' dtest <- with(agaricus.test, xgb.DMatrix(data, label = label, nthread = 2))
+#' ## Keep the number of threads to 1 for examples
+#' nthread <- 1
+#' data.table::setDTthreads(nthread)
+#'
+#' dtrain <- with(
+#'   agaricus.train, xgb.DMatrix(data, label = label, nthread = nthread)
+#' )
+#' dtest <- with(
+#'   agaricus.test, xgb.DMatrix(data, label = label, nthread = nthread)
+#' )
 #' watchlist <- list(train = dtrain, eval = dtest)
 #'
 #' ## A simple xgb.train example:
-#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2,
+#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread,
 #'               objective = "binary:logistic", eval_metric = "auc")
 #' bst <- xgb.train(param, dtrain, nrounds = 2, watchlist)
 #'
-#'
-#' ## An xgb.train example where custom objective and evaluation metric are used:
+#' ## An xgb.train example where custom objective and evaluation metric are
+#' ## used:
 #' logregobj <- function(preds, dtrain) {
 #'    labels <- getinfo(dtrain, "label")
 #'    preds <- 1/(1 + exp(-preds))
@@ -218,12 +227,12 @@
 #'
 #' # These functions could be used by passing them either:
 #' #  as 'objective' and 'eval_metric' parameters in the params list:
-#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2,
+#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread,
 #'               objective = logregobj, eval_metric = evalerror)
 #' bst <- xgb.train(param, dtrain, nrounds = 2, watchlist)
 #'
 #' #  or through the ... arguments:
-#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2)
+#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread)
 #' bst <- xgb.train(param, dtrain, nrounds = 2, watchlist,
 #'                  objective = logregobj, eval_metric = evalerror)
 #'
@@ -233,7 +242,7 @@
 #'
 #'
 #' ## An xgb.train example of using variable learning rates at each iteration:
-#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2,
+#' param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread,
 #'               objective = "binary:logistic", eval_metric = "auc")
 #' my_etas <- list(eta = c(0.5, 0.1))
 #' bst <- xgb.train(param, dtrain, nrounds = 2, watchlist,
@@ -245,7 +254,7 @@
 #'
 #' ## An 'xgboost' interface example:
 #' bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label,
-#'                max_depth = 2, eta = 1, nthread = 2, nrounds = 2,
+#'                max_depth = 2, eta = 1, nthread = nthread, nrounds = 2,
 #'                objective = "binary:logistic")
 #' pred <- predict(bst, agaricus.test$data)
 #'
diff --git a/R-package/man/cb.gblinear.history.Rd b/R-package/man/cb.gblinear.history.Rd
index fbbb56dbefc6..309100673584 100644
--- a/R-package/man/cb.gblinear.history.Rd
+++ b/R-package/man/cb.gblinear.history.Rd
@@ -35,14 +35,18 @@ Callback function expects the following values to be set in its calling frame:
 }
 \examples{
 #### Binary classification:
-#
+
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
+
 # In the iris dataset, it is hard to linearly separate Versicolor class from the rest
 # without considering the 2nd order interactions:
 x <- model.matrix(Species ~ .^2, iris)[,-1]
 colnames(x)
-dtrain <- xgb.DMatrix(scale(x), label = 1*(iris$Species == "versicolor"), nthread = 2)
+dtrain <- xgb.DMatrix(scale(x), label = 1*(iris$Species == "versicolor"), nthread = nthread)
 param <- list(booster = "gblinear", objective = "reg:logistic", eval_metric = "auc",
-              lambda = 0.0003, alpha = 0.0003, nthread = 2)
+              lambda = 0.0003, alpha = 0.0003, nthread = nthread)
 # For 'shotgun', which is a default linear updater, using high eta values may result in
 # unstable behaviour in some datasets. With this simple dataset, however, the high learning
 # rate does not break the convergence, but allows us to illustrate the typical pattern of
@@ -72,9 +76,9 @@ matplot(xgb.gblinear.history(bst)[[3]], type = 'l')
 
 #### Multiclass classification:
 #
-dtrain <- xgb.DMatrix(scale(x), label = as.numeric(iris$Species) - 1, nthread = 2)
+dtrain <- xgb.DMatrix(scale(x), label = as.numeric(iris$Species) - 1, nthread = nthread)
 param <- list(booster = "gblinear", objective = "multi:softprob", num_class = 3,
-              lambda = 0.0003, alpha = 0.0003, nthread = 2)
+              lambda = 0.0003, alpha = 0.0003, nthread = nthread)
 # For the default linear updater 'shotgun' it sometimes is helpful
 # to use smaller eta to reduce instability
 bst <- xgb.train(param, dtrain, list(tr=dtrain), nrounds = 70, eta = 0.5,
diff --git a/R-package/man/predict.xgb.Booster.Rd b/R-package/man/predict.xgb.Booster.Rd
index 067cbf20741e..075bcc7af193 100644
--- a/R-package/man/predict.xgb.Booster.Rd
+++ b/R-package/man/predict.xgb.Booster.Rd
@@ -128,11 +128,16 @@ of the most important features first. See below about the format of the returned
 
 data(agaricus.train, package='xgboost')
 data(agaricus.test, package='xgboost')
+
+## Keep the number of threads to 2 for examples
+nthread <- 2
+data.table::setDTthreads(nthread)
+
 train <- agaricus.train
 test <- agaricus.test
 
 bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-               eta = 0.5, nthread = 2, nrounds = 5, objective = "binary:logistic")
+               eta = 0.5, nthread = nthread, nrounds = 5, objective = "binary:logistic")
 # use all trees by default
 pred <- predict(bst, test$data)
 # use only the 1st tree
diff --git a/R-package/man/xgb.DMatrix.Rd b/R-package/man/xgb.DMatrix.Rd
index 742073faddc4..59ef0b3be976 100644
--- a/R-package/man/xgb.DMatrix.Rd
+++ b/R-package/man/xgb.DMatrix.Rd
@@ -38,7 +38,12 @@ Supported input file formats are either a LIBSVM text file or a binary file that
 }
 \examples{
 data(agaricus.train, package='xgboost')
-dtrain <- with(agaricus.train, xgb.DMatrix(data, label = label, nthread = 2))
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
+dtrain <- with(
+  agaricus.train, xgb.DMatrix(data, label = label, nthread = nthread)
+)
 xgb.DMatrix.save(dtrain, 'xgb.DMatrix.data')
 dtrain <- xgb.DMatrix('xgb.DMatrix.data')
 if (file.exists('xgb.DMatrix.data')) file.remove('xgb.DMatrix.data')
diff --git a/R-package/man/xgb.config.Rd b/R-package/man/xgb.config.Rd
index a5187c8eafe4..35545cc77a50 100644
--- a/R-package/man/xgb.config.Rd
+++ b/R-package/man/xgb.config.Rd
@@ -19,10 +19,15 @@ Accessors for model parameters as JSON string.
 }
 \examples{
 data(agaricus.train, package='xgboost')
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
 train <- agaricus.train
 
-bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-               eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic")
+bst <- xgboost(
+  data = train$data, label = train$label, max_depth = 2,
+  eta = 1, nthread = nthread, nrounds = 2, objective = "binary:logistic"
+)
 config <- xgb.config(bst)
 
 }
diff --git a/R-package/man/xgb.load.Rd b/R-package/man/xgb.load.Rd
index f644bc408125..1a406cc21d0e 100644
--- a/R-package/man/xgb.load.Rd
+++ b/R-package/man/xgb.load.Rd
@@ -27,14 +27,23 @@ not \code{xgb.load}.
 \examples{
 data(agaricus.train, package='xgboost')
 data(agaricus.test, package='xgboost')
+
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
+
 train <- agaricus.train
 test <- agaricus.test
-bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-               eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+bst <- xgboost(
+  data = train$data, label = train$label, max_depth = 2, eta = 1,
+  nthread = nthread,
+  nrounds = 2,
+  objective = "binary:logistic"
+)
+
 xgb.save(bst, 'xgb.model')
 bst <- xgb.load('xgb.model')
 if (file.exists('xgb.model')) file.remove('xgb.model')
-pred <- predict(bst, test$data)
 }
 \seealso{
 \code{\link{xgb.save}}, \code{\link{xgb.Booster.complete}}.
diff --git a/R-package/man/xgb.model.dt.tree.Rd b/R-package/man/xgb.model.dt.tree.Rd
index b89d298b6298..5a17f9d90f62 100644
--- a/R-package/man/xgb.model.dt.tree.Rd
+++ b/R-package/man/xgb.model.dt.tree.Rd
@@ -66,9 +66,12 @@ Parse a boosted tree model text dump into a \code{data.table} structure.
 # Basic use:
 
 data(agaricus.train, package='xgboost')
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
 
 bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 2,
-               eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+               eta = 1, nthread = nthread, nrounds = 2,objective = "binary:logistic")
 
 (dt <- xgb.model.dt.tree(colnames(agaricus.train$data), bst))
 
diff --git a/R-package/man/xgb.plot.deepness.Rd b/R-package/man/xgb.plot.deepness.Rd
index 39e291a811cc..9e23ac130ea0 100644
--- a/R-package/man/xgb.plot.deepness.Rd
+++ b/R-package/man/xgb.plot.deepness.Rd
@@ -61,10 +61,13 @@ This function was inspired by the blog post
 \examples{
 
 data(agaricus.train, package='xgboost')
+## Keep the number of threads to 2 for examples
+nthread <- 2
+data.table::setDTthreads(nthread)
 
-# Change max_depth to a higher number to get a more significant result
+## Change max_depth to a higher number to get a more significant result
 bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 6,
-               eta = 0.1, nthread = 2, nrounds = 50, objective = "binary:logistic",
+               eta = 0.1, nthread = nthread, nrounds = 50, objective = "binary:logistic",
                subsample = 0.5, min_child_weight = 2)
 
 xgb.plot.deepness(bst)
diff --git a/R-package/man/xgb.plot.importance.Rd b/R-package/man/xgb.plot.importance.Rd
index 691a8fdfca66..7a49061a7528 100644
--- a/R-package/man/xgb.plot.importance.Rd
+++ b/R-package/man/xgb.plot.importance.Rd
@@ -77,9 +77,14 @@ with bar colors corresponding to different clusters that have somewhat similar i
 }
 \examples{
 data(agaricus.train)
+## Keep the number of threads to 2 for examples
+nthread <- 2
+data.table::setDTthreads(nthread)
 
-bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 3,
-               eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic")
+bst <- xgboost(
+  data = agaricus.train$data, label = agaricus.train$label, max_depth = 3,
+  eta = 1, nthread = nthread, nrounds = 2, objective = "binary:logistic"
+)
 
 importance_matrix <- xgb.importance(colnames(agaricus.train$data), model = bst)
 
diff --git a/R-package/man/xgb.plot.multi.trees.Rd b/R-package/man/xgb.plot.multi.trees.Rd
index 74c4a06040f7..4fa526b90c22 100644
--- a/R-package/man/xgb.plot.multi.trees.Rd
+++ b/R-package/man/xgb.plot.multi.trees.Rd
@@ -63,10 +63,15 @@ This function is inspired by this blog post:
 \examples{
 
 data(agaricus.train, package='xgboost')
-
-bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label, max_depth = 15,
-               eta = 1, nthread = 2, nrounds = 30, objective = "binary:logistic",
-               min_child_weight = 50, verbose = 0)
+## Keep the number of threads to 2 for examples
+nthread <- 2
+data.table::setDTthreads(nthread)
+
+bst <- xgboost(
+  data = agaricus.train$data, label = agaricus.train$label, max_depth = 15,
+  eta = 1, nthread = nthread, nrounds = 30, objective = "binary:logistic",
+  min_child_weight = 50, verbose = 0
+)
 
 p <- xgb.plot.multi.trees(model = bst, features_keep = 3)
 print(p)
diff --git a/R-package/man/xgb.plot.shap.Rd b/R-package/man/xgb.plot.shap.Rd
index a55a551de286..6f2d0dfa6dcc 100644
--- a/R-package/man/xgb.plot.shap.Rd
+++ b/R-package/man/xgb.plot.shap.Rd
@@ -124,9 +124,14 @@ a meaningful thing to do.
 data(agaricus.train, package='xgboost')
 data(agaricus.test, package='xgboost')
 
-bst <- xgboost(agaricus.train$data, agaricus.train$label, nrounds = 50,
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
+nrounds <- 20
+
+bst <- xgboost(agaricus.train$data, agaricus.train$label, nrounds = nrounds,
                eta = 0.1, max_depth = 3, subsample = .5,
-               method = "hist", objective = "binary:logistic", nthread = 2, verbose = 0)
+               method = "hist", objective = "binary:logistic", nthread = nthread, verbose = 0)
 
 xgb.plot.shap(agaricus.test$data, model = bst, features = "odor=none")
 contr <- predict(bst, agaricus.test$data, predcontrib = TRUE)
@@ -135,12 +140,11 @@ xgb.ggplot.shap.summary(agaricus.test$data, contr, model = bst, top_n = 12)  # S
 
 # multiclass example - plots for each class separately:
 nclass <- 3
-nrounds <- 20
 x <- as.matrix(iris[, -5])
 set.seed(123)
 is.na(x[sample(nrow(x) * 4, 30)]) <- TRUE # introduce some missing values
 mbst <- xgboost(data = x, label = as.numeric(iris$Species) - 1, nrounds = nrounds,
-                max_depth = 2, eta = 0.3, subsample = .5, nthread = 2,
+                max_depth = 2, eta = 0.3, subsample = .5, nthread = nthread,
                 objective = "multi:softprob", num_class = nclass, verbose = 0)
 trees0 <- seq(from=0, by=nclass, length.out=nrounds)
 col <- rgb(0, 0, 1, 0.5)
diff --git a/R-package/man/xgb.plot.shap.summary.Rd b/R-package/man/xgb.plot.shap.summary.Rd
index f757fd7404a3..790237efe92d 100644
--- a/R-package/man/xgb.plot.shap.summary.Rd
+++ b/R-package/man/xgb.plot.shap.summary.Rd
@@ -67,12 +67,12 @@ Each point (observation) is coloured based on its feature value. The plot
 hence allows us to see which features have a negative / positive contribution
 on the model prediction, and whether the contribution is different for larger
 or smaller values of the feature. We effectively try to replicate the
-\code{summary_plot} function from https://github.com/slundberg/shap.
+\code{summary_plot} function from https://github.com/shap/shap
 }
 \examples{
 # See \code{\link{xgb.plot.shap}}.
 }
 \seealso{
 \code{\link{xgb.plot.shap}}, \code{\link{xgb.ggplot.shap.summary}},
-  \url{https://github.com/slundberg/shap}
+  \url{https://github.com/shap/shap}
 }
diff --git a/R-package/man/xgb.save.Rd b/R-package/man/xgb.save.Rd
index 235fc504c9ed..a7e160a12a9b 100644
--- a/R-package/man/xgb.save.Rd
+++ b/R-package/man/xgb.save.Rd
@@ -31,14 +31,22 @@ releases of XGBoost.
 \examples{
 data(agaricus.train, package='xgboost')
 data(agaricus.test, package='xgboost')
+
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
+
 train <- agaricus.train
 test <- agaricus.test
-bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-               eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+bst <- xgboost(
+  data = train$data, label = train$label, max_depth = 2, eta = 1,
+  nthread = nthread,
+  nrounds = 2,
+  objective = "binary:logistic"
+)
 xgb.save(bst, 'xgb.model')
 bst <- xgb.load('xgb.model')
 if (file.exists('xgb.model')) file.remove('xgb.model')
-pred <- predict(bst, test$data)
 }
 \seealso{
 \code{\link{xgb.load}}, \code{\link{xgb.Booster.complete}}.
diff --git a/R-package/man/xgb.save.raw.Rd b/R-package/man/xgb.save.raw.Rd
index ad188eb831ba..c7c93a734d22 100644
--- a/R-package/man/xgb.save.raw.Rd
+++ b/R-package/man/xgb.save.raw.Rd
@@ -25,12 +25,17 @@ Save xgboost model from xgboost or xgb.train
 \examples{
 data(agaricus.train, package='xgboost')
 data(agaricus.test, package='xgboost')
+
+## Keep the number of threads to 2 for examples
+nthread <- 2
+data.table::setDTthreads(nthread)
+
 train <- agaricus.train
 test <- agaricus.test
 bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-               eta = 1, nthread = 2, nrounds = 2,objective = "binary:logistic")
+               eta = 1, nthread = nthread, nrounds = 2,objective = "binary:logistic")
+
 raw <- xgb.save.raw(bst)
 bst <- xgb.load.raw(raw)
-pred <- predict(bst, test$data)
 
 }
diff --git a/R-package/man/xgb.train.Rd b/R-package/man/xgb.train.Rd
index b2c9f9edc6d4..574f0e360cdc 100644
--- a/R-package/man/xgb.train.Rd
+++ b/R-package/man/xgb.train.Rd
@@ -206,7 +206,8 @@ customized objective and evaluation metric functions, therefore it is more flexi
 than the \code{xgboost} interface.
 
 Parallelization is automatically enabled if \code{OpenMP} is present.
-Number of threads can also be manually specified via \code{nthread} parameter.
+Number of threads can also be manually specified via the \code{nthread}
+parameter.
 
 The evaluation metric is chosen automatically by XGBoost (according to the objective)
 when the \code{eval_metric} parameter is not provided.
@@ -241,17 +242,25 @@ The following callbacks are automatically created when certain parameters are se
 data(agaricus.train, package='xgboost')
 data(agaricus.test, package='xgboost')
 
-dtrain <- with(agaricus.train, xgb.DMatrix(data, label = label, nthread = 2))
-dtest <- with(agaricus.test, xgb.DMatrix(data, label = label, nthread = 2))
+## Keep the number of threads to 1 for examples
+nthread <- 1
+data.table::setDTthreads(nthread)
+
+dtrain <- with(
+  agaricus.train, xgb.DMatrix(data, label = label, nthread = nthread)
+)
+dtest <- with(
+  agaricus.test, xgb.DMatrix(data, label = label, nthread = nthread)
+)
 watchlist <- list(train = dtrain, eval = dtest)
 
 ## A simple xgb.train example:
-param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2,
+param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread,
               objective = "binary:logistic", eval_metric = "auc")
 bst <- xgb.train(param, dtrain, nrounds = 2, watchlist)
 
-
-## An xgb.train example where custom objective and evaluation metric are used:
+## An xgb.train example where custom objective and evaluation metric are
+## used:
 logregobj <- function(preds, dtrain) {
    labels <- getinfo(dtrain, "label")
    preds <- 1/(1 + exp(-preds))
@@ -267,12 +276,12 @@ evalerror <- function(preds, dtrain) {
 
 # These functions could be used by passing them either:
 #  as 'objective' and 'eval_metric' parameters in the params list:
-param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2,
+param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread,
               objective = logregobj, eval_metric = evalerror)
 bst <- xgb.train(param, dtrain, nrounds = 2, watchlist)
 
 #  or through the ... arguments:
-param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2)
+param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread)
 bst <- xgb.train(param, dtrain, nrounds = 2, watchlist,
                  objective = logregobj, eval_metric = evalerror)
 
@@ -282,7 +291,7 @@ bst <- xgb.train(param, dtrain, nrounds = 2, watchlist,
 
 
 ## An xgb.train example of using variable learning rates at each iteration:
-param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = 2,
+param <- list(max_depth = 2, eta = 1, verbose = 0, nthread = nthread,
               objective = "binary:logistic", eval_metric = "auc")
 my_etas <- list(eta = c(0.5, 0.1))
 bst <- xgb.train(param, dtrain, nrounds = 2, watchlist,
@@ -294,7 +303,7 @@ bst <- xgb.train(param, dtrain, nrounds = 25, watchlist,
 
 ## An 'xgboost' interface example:
 bst <- xgboost(data = agaricus.train$data, label = agaricus.train$label,
-               max_depth = 2, eta = 1, nthread = 2, nrounds = 2,
+               max_depth = 2, eta = 1, nthread = nthread, nrounds = 2,
                objective = "binary:logistic")
 pred <- predict(bst, agaricus.test$data)
 
diff --git a/R-package/src/xgboost_R.cc b/R-package/src/xgboost_R.cc
index fcbccb145852..a55bd07856b9 100644
--- a/R-package/src/xgboost_R.cc
+++ b/R-package/src/xgboost_R.cc
@@ -1,5 +1,5 @@
 /**
- * Copyright 2014-2022 by XGBoost Contributors
+ * Copyright 2014-2023, XGBoost Contributors
  */
 #include <dmlc/common.h>
 #include <dmlc/omp.h>
@@ -12,7 +12,6 @@
 #include <cstring>
 #include <sstream>
 #include <string>
-#include <utility>
 #include <vector>
 
 #include "../../src/c_api/c_api_error.h"
@@ -26,22 +25,24 @@
 #define R_API_BEGIN()                           \
   GetRNGstate();                                \
   try {
+
 /*!
  * \brief macro to annotate end of api
  */
-#define R_API_END()                             \
-  } catch(dmlc::Error& e) {                     \
-    PutRNGstate();                              \
-    error(e.what());                            \
-  }                                             \
+#define R_API_END()           \
+  }                           \
+  catch (dmlc::Error & e) {   \
+    PutRNGstate();            \
+    Rf_error("%s", e.what()); \
+  }                           \
   PutRNGstate();
 
 /*!
  * \brief macro to check the call.
  */
-#define CHECK_CALL(x)                           \
-  if ((x) != 0) {                               \
-    error(XGBGetLastError());                   \
+#define CHECK_CALL(x)                  \
+  if ((x) != 0) {                      \
+    Rf_error("%s", XGBGetLastError()); \
   }
 
 using dmlc::BeginPtr;
@@ -116,11 +117,25 @@ XGB_DLL SEXP XGDMatrixCreateFromMat_R(SEXP mat, SEXP missing, SEXP n_threads) {
   std::vector<float> data(nrow * ncol);
   int32_t threads = xgboost::common::OmpGetNumThreads(asInteger(n_threads));
 
-  xgboost::common::ParallelFor(nrow, threads, [&](xgboost::omp_ulong i) {
-    for (size_t j = 0; j < ncol; ++j) {
-      data[i * ncol + j] = is_int ? static_cast<float>(iin[i + nrow * j]) : din[i + nrow * j];
-    }
-  });
+  if (is_int) {
+    xgboost::common::ParallelFor(nrow, threads, [&](xgboost::omp_ulong i) {
+      for (size_t j = 0; j < ncol; ++j) {
+        auto v = iin[i + nrow * j];
+        if (v == NA_INTEGER) {
+          data[i * ncol + j] = std::numeric_limits<float>::quiet_NaN();
+        } else {
+          data[i * ncol + j] = static_cast<float>(v);
+        }
+      }
+    });
+  } else {
+    xgboost::common::ParallelFor(nrow, threads, [&](xgboost::omp_ulong i) {
+      for (size_t j = 0; j < ncol; ++j) {
+        data[i * ncol + j] = din[i + nrow * j];
+      }
+    });
+  }
+
   DMatrixHandle handle;
   CHECK_CALL(XGDMatrixCreateFromMat_omp(BeginPtr(data), nrow, ncol,
                                         asReal(missing), &handle, threads));
diff --git a/R-package/tests/helper_scripts/run-examples.R b/R-package/tests/helper_scripts/run-examples.R
new file mode 100644
index 000000000000..08dd3d2a05a1
--- /dev/null
+++ b/R-package/tests/helper_scripts/run-examples.R
@@ -0,0 +1,25 @@
+## Helper script for running individual examples.
+library(pkgload)
+library(xgboost)
+
+files <- list.files("./man")
+
+
+run_example_timeit <- function(f) {
+  path <- paste("./man/", f, sep = "")
+  print(paste("Test", f))
+  flush.console()
+  t0 <- proc.time()
+  run_example(path)
+  t1 <- proc.time()
+  list(file = f, time = t1 - t0)
+}
+
+timings <- lapply(files, run_example_timeit)
+
+for (t in timings) {
+  ratio <- t$time[1] / t$time[3]
+  if (!is.na(ratio) && !is.infinite(ratio) && ratio >= 2.5) {
+    print(paste("Offending example:", t$file, ratio))
+  }
+}
diff --git a/R-package/tests/testthat/test_basic.R b/R-package/tests/testthat/test_basic.R
index ad8c8a830955..291862e56ae6 100644
--- a/R-package/tests/testthat/test_basic.R
+++ b/R-package/tests/testthat/test_basic.R
@@ -3,24 +3,29 @@ library(Matrix)
 
 context("basic functions")
 
-data(agaricus.train, package = 'xgboost')
-data(agaricus.test, package = 'xgboost')
+data(agaricus.train, package = "xgboost")
+data(agaricus.test, package = "xgboost")
 train <- agaricus.train
 test <- agaricus.test
 set.seed(1994)
 
 # disable some tests for Win32
 windows_flag <- .Platform$OS.type == "windows" &&
-               .Machine$sizeof.pointer != 8
-solaris_flag <- (Sys.info()['sysname'] == "SunOS")
+  .Machine$sizeof.pointer != 8
+solaris_flag <- (Sys.info()["sysname"] == "SunOS")
+n_threads <- 1
+
 
 test_that("train and predict binary classification", {
   nrounds <- 2
   expect_output(
-    bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-                  eta = 1, nthread = 2, nrounds = nrounds, objective = "binary:logistic",
-                  eval_metric = "error")
-  , "train-error")
+    bst <- xgboost(
+      data = train$data, label = train$label, max_depth = 2,
+      eta = 1, nthread = n_threads, nrounds = nrounds,
+      objective = "binary:logistic", eval_metric = "error"
+    ),
+    "train-error"
+  )
   expect_equal(class(bst), "xgb.Booster")
   expect_equal(bst$niter, nrounds)
   expect_false(is.null(bst$evaluation_log))
@@ -49,26 +54,39 @@ test_that("parameter validation works", {
   d <- cbind(
     x1 = rnorm(10),
     x2 = rnorm(10),
-    x3 = rnorm(10))
+    x3 = rnorm(10)
+  )
   y <- d[, "x1"] + d[, "x2"]^2 +
     ifelse(d[, "x3"] > .5, d[, "x3"]^2, 2^d[, "x3"]) +
     rnorm(10)
-  dtrain <- xgb.DMatrix(data = d, info = list(label = y))
+  dtrain <- xgb.DMatrix(data = d, info = list(label = y), nthread = n_threads)
 
   correct <- function() {
-    params <- list(max_depth = 2, booster = "dart",
-                   rate_drop = 0.5, one_drop = TRUE,
-                   objective = "reg:squarederror")
+    params <- list(
+      max_depth = 2,
+      booster = "dart",
+      rate_drop = 0.5,
+      one_drop = TRUE,
+      nthread = n_threads,
+      objective = "reg:squarederror"
+    )
     xgb.train(params = params, data = dtrain, nrounds = nrounds)
   }
   expect_silent(correct())
   incorrect <- function() {
-    params <- list(max_depth = 2, booster = "dart",
-                   rate_drop = 0.5, one_drop = TRUE,
-                   objective = "reg:squarederror",
-                   foo = "bar", bar = "foo")
+    params <- list(
+      max_depth = 2,
+      booster = "dart",
+      rate_drop = 0.5,
+      one_drop = TRUE,
+      objective = "reg:squarederror",
+      nthread = n_threads,
+      foo = "bar",
+      bar = "foo"
+    )
     output <- capture.output(
-      xgb.train(params = params, data = dtrain, nrounds = nrounds))
+      xgb.train(params = params, data = dtrain, nrounds = nrounds)
+    )
     print(output)
   }
   expect_output(incorrect(), '\\\\"bar\\\\", \\\\"foo\\\\"')
@@ -82,15 +100,25 @@ test_that("dart prediction works", {
   d <- cbind(
     x1 = rnorm(100),
     x2 = rnorm(100),
-    x3 = rnorm(100))
+    x3 = rnorm(100)
+  )
   y <- d[, "x1"] + d[, "x2"]^2 +
     ifelse(d[, "x3"] > .5, d[, "x3"]^2, 2^d[, "x3"]) +
     rnorm(100)
 
   set.seed(1994)
-  booster_by_xgboost <- xgboost(data = d, label = y, max_depth = 2, booster = "dart",
-                                rate_drop = 0.5, one_drop = TRUE,
-                                eta = 1, nthread = 2, nrounds = nrounds, objective = "reg:squarederror")
+  booster_by_xgboost <- xgboost(
+    data = d,
+    label = y,
+    max_depth = 2,
+    booster = "dart",
+    rate_drop = 0.5,
+    one_drop = TRUE,
+    eta = 1,
+    nthread = n_threads,
+    nrounds = nrounds,
+    objective = "reg:squarederror"
+  )
   pred_by_xgboost_0 <- predict(booster_by_xgboost, newdata = d, ntreelimit = 0)
   pred_by_xgboost_1 <- predict(booster_by_xgboost, newdata = d, ntreelimit = nrounds)
   expect_true(all(matrix(pred_by_xgboost_0, byrow = TRUE) == matrix(pred_by_xgboost_1, byrow = TRUE)))
@@ -99,20 +127,20 @@ test_that("dart prediction works", {
   expect_false(all(matrix(pred_by_xgboost_0, byrow = TRUE) == matrix(pred_by_xgboost_2, byrow = TRUE)))
 
   set.seed(1994)
-  dtrain <- xgb.DMatrix(data = d, info = list(label = y))
-  booster_by_train <- xgb.train(params = list(
-                                    booster = "dart",
-                                    max_depth = 2,
-                                    eta = 1,
-                                    rate_drop = 0.5,
-                                    one_drop = TRUE,
-                                    nthread = 1,
-                                    tree_method = "exact",
-                                    objective = "reg:squarederror"
-                                ),
-                                data = dtrain,
-                                nrounds = nrounds
-                                )
+  dtrain <- xgb.DMatrix(data = d, info = list(label = y), nthread = n_threads)
+  booster_by_train <- xgb.train(
+    params = list(
+      booster = "dart",
+      max_depth = 2,
+      eta = 1,
+      rate_drop = 0.5,
+      one_drop = TRUE,
+      nthread = n_threads,
+      objective = "reg:squarederror"
+    ),
+    data = dtrain,
+    nrounds = nrounds
+  )
   pred_by_train_0 <- predict(booster_by_train, newdata = dtrain, ntreelimit = 0)
   pred_by_train_1 <- predict(booster_by_train, newdata = dtrain, ntreelimit = nrounds)
   pred_by_train_2 <- predict(booster_by_train, newdata = dtrain, training = TRUE)
@@ -126,10 +154,13 @@ test_that("train and predict softprob", {
   lb <- as.numeric(iris$Species) - 1
   set.seed(11)
   expect_output(
-    bst <- xgboost(data = as.matrix(iris[, -5]), label = lb,
-                   max_depth = 3, eta = 0.5, nthread = 2, nrounds = 5,
-                   objective = "multi:softprob", num_class = 3, eval_metric = "merror")
-  , "train-merror")
+    bst <- xgboost(
+      data = as.matrix(iris[, -5]), label = lb,
+      max_depth = 3, eta = 0.5, nthread = n_threads, nrounds = 5,
+      objective = "multi:softprob", num_class = 3, eval_metric = "merror"
+    ),
+    "train-merror"
+  )
   expect_false(is.null(bst$evaluation_log))
   expect_lt(bst$evaluation_log[, min(train_merror)], 0.025)
   expect_equal(bst$niter * 3, xgb.ntree(bst))
@@ -158,9 +189,10 @@ test_that("train and predict softprob", {
     x3 = rnorm(100)
   )
   y <- sample.int(10, 100, replace = TRUE) - 1
-  dtrain <- xgb.DMatrix(data = d, info = list(label = y))
+  dtrain <- xgb.DMatrix(data = d, info = list(label = y), nthread = n_threads)
   booster <- xgb.train(
-    params = list(tree_method = "hist"), data = dtrain, nrounds = 4, num_class = 10,
+    params = list(tree_method = "hist", nthread = n_threads),
+    data = dtrain, nrounds = 4, num_class = 10,
     objective = "multi:softprob"
   )
   predt <- predict(booster, as.matrix(d), reshape = TRUE, strict_shape = FALSE)
@@ -172,10 +204,13 @@ test_that("train and predict softmax", {
   lb <- as.numeric(iris$Species) - 1
   set.seed(11)
   expect_output(
-    bst <- xgboost(data = as.matrix(iris[, -5]), label = lb,
-                   max_depth = 3, eta = 0.5, nthread = 2, nrounds = 5,
-                   objective = "multi:softmax", num_class = 3, eval_metric = "merror")
-  , "train-merror")
+    bst <- xgboost(
+      data = as.matrix(iris[, -5]), label = lb,
+      max_depth = 3, eta = 0.5, nthread = n_threads, nrounds = 5,
+      objective = "multi:softmax", num_class = 3, eval_metric = "merror"
+    ),
+    "train-merror"
+  )
   expect_false(is.null(bst$evaluation_log))
   expect_lt(bst$evaluation_log[, min(train_merror)], 0.025)
   expect_equal(bst$niter * 3, xgb.ntree(bst))
@@ -190,16 +225,19 @@ test_that("train and predict RF", {
   set.seed(11)
   lb <- train$label
   # single iteration
-  bst <- xgboost(data = train$data, label = lb, max_depth = 5,
-                 nthread = 2, nrounds = 1, objective = "binary:logistic", eval_metric = "error",
-                 num_parallel_tree = 20, subsample = 0.6, colsample_bytree = 0.1)
+  bst <- xgboost(
+    data = train$data, label = lb, max_depth = 5,
+    nthread = n_threads,
+    nrounds = 1, objective = "binary:logistic", eval_metric = "error",
+    num_parallel_tree = 20, subsample = 0.6, colsample_bytree = 0.1
+  )
   expect_equal(bst$niter, 1)
   expect_equal(xgb.ntree(bst), 20)
 
   pred <- predict(bst, train$data)
   pred_err <- sum((pred > 0.5) != lb) / length(lb)
   expect_lt(abs(bst$evaluation_log[1, train_error] - pred_err), 10e-6)
-  #expect_lt(pred_err, 0.03)
+  # expect_lt(pred_err, 0.03)
 
   pred <- predict(bst, train$data, ntreelimit = 20)
   pred_err_20 <- sum((pred > 0.5) != lb) / length(lb)
@@ -213,11 +251,13 @@ test_that("train and predict RF with softprob", {
   lb <- as.numeric(iris$Species) - 1
   nrounds <- 15
   set.seed(11)
-  bst <- xgboost(data = as.matrix(iris[, -5]), label = lb,
-                 max_depth = 3, eta = 0.9, nthread = 2, nrounds = nrounds,
-                 objective = "multi:softprob", eval_metric = "merror",
-                 num_class = 3, verbose = 0,
-                 num_parallel_tree = 4, subsample = 0.5, colsample_bytree = 0.5)
+  bst <- xgboost(
+    data = as.matrix(iris[, -5]), label = lb,
+    max_depth = 3, eta = 0.9, nthread = n_threads, nrounds = nrounds,
+    objective = "multi:softprob", eval_metric = "merror",
+    num_class = 3, verbose = 0,
+    num_parallel_tree = 4, subsample = 0.5, colsample_bytree = 0.5
+  )
   expect_equal(bst$niter, 15)
   expect_equal(xgb.ntree(bst), 15 * 3 * 4)
   # predict for all iterations:
@@ -234,10 +274,13 @@ test_that("train and predict RF with softprob", {
 
 test_that("use of multiple eval metrics works", {
   expect_output(
-    bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-                  eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic",
-                  eval_metric = 'error', eval_metric = 'auc', eval_metric = "logloss")
-  , "train-error.*train-auc.*train-logloss")
+    bst <- xgboost(
+      data = train$data, label = train$label, max_depth = 2,
+      eta = 1, nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+      eval_metric = "error", eval_metric = "auc", eval_metric = "logloss"
+    ),
+    "train-error.*train-auc.*train-logloss"
+  )
   expect_false(is.null(bst$evaluation_log))
   expect_equal(dim(bst$evaluation_log), c(2, 4))
   expect_equal(colnames(bst$evaluation_log), c("iter", "train_error", "train_auc", "train_logloss"))
@@ -245,9 +288,11 @@ test_that("use of multiple eval metrics works", {
 
 
 test_that("training continuation works", {
-  dtrain <- xgb.DMatrix(train$data, label = train$label)
+  dtrain <- xgb.DMatrix(train$data, label = train$label, nthread = n_threads)
   watchlist <- list(train = dtrain)
-  param <- list(objective = "binary:logistic", max_depth = 2, eta = 1, nthread = 2)
+  param <- list(
+    objective = "binary:logistic", max_depth = 2, eta = 1, nthread = n_threads
+  )
 
   # for the reference, use 4 iterations at once:
   set.seed(11)
@@ -257,30 +302,33 @@ test_that("training continuation works", {
   bst1 <- xgb.train(param, dtrain, nrounds = 2, watchlist, verbose = 0)
   # continue for two more:
   bst2 <- xgb.train(param, dtrain, nrounds = 2, watchlist, verbose = 0, xgb_model = bst1)
-  if (!windows_flag && !solaris_flag)
+  if (!windows_flag && !solaris_flag) {
     expect_equal(bst$raw, bst2$raw)
+  }
   expect_false(is.null(bst2$evaluation_log))
   expect_equal(dim(bst2$evaluation_log), c(4, 2))
   expect_equal(bst2$evaluation_log, bst$evaluation_log)
   # test continuing from raw model data
   bst2 <- xgb.train(param, dtrain, nrounds = 2, watchlist, verbose = 0, xgb_model = bst1$raw)
-  if (!windows_flag && !solaris_flag)
+  if (!windows_flag && !solaris_flag) {
     expect_equal(bst$raw, bst2$raw)
+  }
   expect_equal(dim(bst2$evaluation_log), c(2, 2))
   # test continuing from a model in file
   xgb.save(bst1, "xgboost.json")
   bst2 <- xgb.train(param, dtrain, nrounds = 2, watchlist, verbose = 0, xgb_model = "xgboost.json")
-  if (!windows_flag && !solaris_flag)
+  if (!windows_flag && !solaris_flag) {
     expect_equal(bst$raw, bst2$raw)
+  }
   expect_equal(dim(bst2$evaluation_log), c(2, 2))
   file.remove("xgboost.json")
 })
 
 test_that("model serialization works", {
   out_path <- "model_serialization"
-  dtrain <- xgb.DMatrix(train$data, label = train$label)
+  dtrain <- xgb.DMatrix(train$data, label = train$label, nthread = n_threads)
   watchlist <- list(train = dtrain)
-  param <- list(objective = "binary:logistic")
+  param <- list(objective = "binary:logistic", nthread = n_threads)
   booster <- xgb.train(param, dtrain, nrounds = 4, watchlist)
   raw <- xgb.serialize(booster)
   saveRDS(raw, out_path)
@@ -295,11 +343,14 @@ test_that("model serialization works", {
 test_that("xgb.cv works", {
   set.seed(11)
   expect_output(
-    cv <- xgb.cv(data = train$data, label = train$label, max_depth = 2, nfold = 5,
-                 eta = 1., nthread = 2, nrounds = 2, objective = "binary:logistic",
-                 eval_metric = "error", verbose = TRUE)
-  , "train-error:")
-  expect_is(cv, 'xgb.cv.synchronous')
+    cv <- xgb.cv(
+      data = train$data, label = train$label, max_depth = 2, nfold = 5,
+      eta = 1., nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+      eval_metric = "error", verbose = TRUE
+    ),
+    "train-error:"
+  )
+  expect_is(cv, "xgb.cv.synchronous")
   expect_false(is.null(cv$evaluation_log))
   expect_lt(cv$evaluation_log[, min(test_error_mean)], 0.03)
   expect_lt(cv$evaluation_log[, min(test_error_std)], 0.008)
@@ -312,15 +363,19 @@ test_that("xgb.cv works", {
 })
 
 test_that("xgb.cv works with stratified folds", {
-  dtrain <- xgb.DMatrix(train$data, label = train$label)
+  dtrain <- xgb.DMatrix(train$data, label = train$label, nthread = n_threads)
   set.seed(314159)
-  cv <- xgb.cv(data = dtrain, max_depth = 2, nfold = 5,
-               eta = 1., nthread = 2, nrounds = 2, objective = "binary:logistic",
-               verbose = TRUE, stratified = FALSE)
+  cv <- xgb.cv(
+    data = dtrain, max_depth = 2, nfold = 5,
+    eta = 1., nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+    verbose = TRUE, stratified = FALSE
+  )
   set.seed(314159)
-  cv2 <- xgb.cv(data = dtrain, max_depth = 2, nfold = 5,
-                eta = 1., nthread = 2, nrounds = 2, objective = "binary:logistic",
-                verbose = TRUE, stratified = TRUE)
+  cv2 <- xgb.cv(
+    data = dtrain, max_depth = 2, nfold = 5,
+    eta = 1., nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+    verbose = TRUE, stratified = TRUE
+  )
   # Stratified folds should result in a different evaluation logs
   expect_true(all(cv$evaluation_log[, test_logloss_mean] != cv2$evaluation_log[, test_logloss_mean]))
 })
@@ -328,40 +383,57 @@ test_that("xgb.cv works with stratified folds", {
 test_that("train and predict with non-strict classes", {
   # standard dense matrix input
   train_dense <- as.matrix(train$data)
-  bst <- xgboost(data = train_dense, label = train$label, max_depth = 2,
-                 eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic", verbose = 0)
+  bst <- xgboost(
+    data = train_dense, label = train$label, max_depth = 2,
+    eta = 1, nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+    verbose = 0
+  )
   pr0 <- predict(bst, train_dense)
 
   # dense matrix-like input of non-matrix class
-  class(train_dense) <- 'shmatrix'
+  class(train_dense) <- "shmatrix"
   expect_true(is.matrix(train_dense))
   expect_error(
-    bst <- xgboost(data = train_dense, label = train$label, max_depth = 2,
-                   eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic", verbose = 0)
-    , regexp = NA)
+    bst <- xgboost(
+      data = train_dense, label = train$label, max_depth = 2,
+      eta = 1, nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+      verbose = 0
+    ),
+    regexp = NA
+  )
   expect_error(pr <- predict(bst, train_dense), regexp = NA)
   expect_equal(pr0, pr)
 
   # dense matrix-like input of non-matrix class with some inheritance
-  class(train_dense) <- c('pphmatrix', 'shmatrix')
+  class(train_dense) <- c("pphmatrix", "shmatrix")
   expect_true(is.matrix(train_dense))
   expect_error(
-    bst <- xgboost(data = train_dense, label = train$label, max_depth = 2,
-                   eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic", verbose = 0)
-    , regexp = NA)
+    bst <- xgboost(
+      data = train_dense, label = train$label, max_depth = 2,
+      eta = 1, nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+      verbose = 0
+    ),
+    regexp = NA
+  )
   expect_error(pr <- predict(bst, train_dense), regexp = NA)
   expect_equal(pr0, pr)
 
   # when someone inherits from xgb.Booster, it should still be possible to use it as xgb.Booster
-  class(bst) <- c('super.Booster', 'xgb.Booster')
+  class(bst) <- c("super.Booster", "xgb.Booster")
   expect_error(pr <- predict(bst, train_dense), regexp = NA)
   expect_equal(pr0, pr)
 })
 
 test_that("max_delta_step works", {
-  dtrain <- xgb.DMatrix(agaricus.train$data, label = agaricus.train$label)
+  dtrain <- xgb.DMatrix(
+    agaricus.train$data, label = agaricus.train$label, nthread = n_threads
+  )
   watchlist <- list(train = dtrain)
-  param <- list(objective = "binary:logistic", eval_metric = "logloss", max_depth = 2, nthread = 2, eta = 0.5)
+  param <- list(
+    objective = "binary:logistic", eval_metric = "logloss", max_depth = 2,
+    nthread = n_threads,
+    eta = 0.5
+  )
   nrounds <- 5
   # model with no restriction on max_delta_step
   bst1 <- xgb.train(param, dtrain, nrounds, watchlist, verbose = 1)
@@ -381,14 +453,16 @@ test_that("colsample_bytree works", {
   test_y <- as.numeric(rowSums(test_x) > 0)
   colnames(train_x) <- paste0("Feature_", sprintf("%03d", 1:100))
   colnames(test_x) <- paste0("Feature_", sprintf("%03d", 1:100))
-  dtrain <- xgb.DMatrix(train_x, label = train_y)
-  dtest <- xgb.DMatrix(test_x, label = test_y)
+  dtrain <- xgb.DMatrix(train_x, label = train_y, nthread = n_threads)
+  dtest <- xgb.DMatrix(test_x, label = test_y, nthread = n_threads)
   watchlist <- list(train = dtrain, eval = dtest)
   ## Use colsample_bytree = 0.01, so that roughly one out of 100 features is chosen for
   ## each tree
-  param <- list(max_depth = 2, eta = 0, nthread = 2,
-                colsample_bytree = 0.01, objective = "binary:logistic",
-                eval_metric = "auc")
+  param <- list(
+    max_depth = 2, eta = 0, nthread = n_threads,
+    colsample_bytree = 0.01, objective = "binary:logistic",
+    eval_metric = "auc"
+  )
   set.seed(2)
   bst <- xgb.train(param, dtrain, nrounds = 100, watchlist, verbose = 0)
   xgb.importance(model = bst)
@@ -398,9 +472,11 @@ test_that("colsample_bytree works", {
 })
 
 test_that("Configuration works", {
-  bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
-                 eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic",
-                 eval_metric = 'error', eval_metric = 'auc', eval_metric = "logloss")
+  bst <- xgboost(
+    data = train$data, label = train$label, max_depth = 2,
+    eta = 1, nthread = n_threads, nrounds = 2, objective = "binary:logistic",
+    eval_metric = "error", eval_metric = "auc", eval_metric = "logloss"
+  )
   config <- xgb.config(bst)
   xgb.config(bst) <- config
   reloaded_config <- xgb.config(bst)
@@ -437,22 +513,26 @@ test_that("strict_shape works", {
     y <- as.numeric(iris$Species) - 1
     X <- as.matrix(iris[, -5])
 
-    bst <- xgboost(data = X, label = y,
-                   max_depth = 2, nrounds = n_rounds,
-                   objective = "multi:softprob", num_class = 3, eval_metric = "merror")
+    bst <- xgboost(
+      data = X, label = y,
+      max_depth = 2, nrounds = n_rounds, nthread = n_threads,
+      objective = "multi:softprob", num_class = 3, eval_metric = "merror"
+    )
 
     test_strict_shape(bst, X, 3)
   }
 
 
   test_agaricus <- function() {
-    data(agaricus.train, package = 'xgboost')
+    data(agaricus.train, package = "xgboost")
     X <- agaricus.train$data
     y <- agaricus.train$label
 
-    bst <- xgboost(data = X, label = y, max_depth = 2,
-                   nrounds = n_rounds, objective = "binary:logistic",
-                   eval_metric = 'error', eval_metric = 'auc', eval_metric = "logloss")
+    bst <- xgboost(
+      data = X, label = y, max_depth = 2, nthread = n_threads,
+      nrounds = n_rounds, objective = "binary:logistic",
+      eval_metric = "error", eval_metric = "auc", eval_metric = "logloss"
+    )
 
     test_strict_shape(bst, X, 1)
   }
@@ -467,8 +547,10 @@ test_that("'predict' accepts CSR data", {
   x_csc <- as(X[1L, , drop = FALSE], "CsparseMatrix")
   x_csr <- as(x_csc, "RsparseMatrix")
   x_spv <- as(x_csc, "sparseVector")
-  bst <- xgboost(data = X, label = y, objective = "binary:logistic",
-                 nrounds = 5L, verbose = FALSE)
+  bst <- xgboost(
+    data = X, label = y, objective = "binary:logistic",
+    nrounds = 5L, verbose = FALSE, nthread = n_threads,
+  )
   p_csc <- predict(bst, x_csc)
   p_csr <- predict(bst, x_csr)
   p_spv <- predict(bst, x_spv)
diff --git a/R-package/tests/testthat/test_callbacks.R b/R-package/tests/testthat/test_callbacks.R
index 69894bd05181..f564839daa6a 100644
--- a/R-package/tests/testthat/test_callbacks.R
+++ b/R-package/tests/testthat/test_callbacks.R
@@ -11,6 +11,8 @@ data(agaricus.test, package = 'xgboost')
 train <- agaricus.train
 test <- agaricus.test
 
+n_threads <- 2
+
 # add some label noise for early stopping tests
 add.noise <- function(label, frac) {
   inoise <- sample(length(label), length(label) * frac)
@@ -20,15 +22,15 @@ add.noise <- function(label, frac) {
 set.seed(11)
 ltrain <- add.noise(train$label, 0.2)
 ltest <- add.noise(test$label, 0.2)
-dtrain <- xgb.DMatrix(train$data, label = ltrain)
-dtest <- xgb.DMatrix(test$data, label = ltest)
+dtrain <- xgb.DMatrix(train$data, label = ltrain, nthread = n_threads)
+dtest <- xgb.DMatrix(test$data, label = ltest, nthread = n_threads)
 watchlist <- list(train = dtrain, test = dtest)
 
 
 err <- function(label, pr) sum((pr > 0.5) != label) / length(label)
 
 param <- list(objective = "binary:logistic", eval_metric = "error",
-              max_depth = 2, nthread = 2)
+              max_depth = 2, nthread = n_threads)
 
 
 test_that("cb.print.evaluation works as expected", {
@@ -108,7 +110,7 @@ test_that("cb.evaluation.log works as expected", {
 
 
 param <- list(objective = "binary:logistic", eval_metric = "error",
-              max_depth = 4, nthread = 2)
+              max_depth = 4, nthread = n_threads)
 
 test_that("can store evaluation_log without printing", {
   expect_silent(
@@ -184,8 +186,10 @@ test_that("cb.save.model works as expected", {
   expect_true(file.exists('xgboost_01.json'))
   expect_true(file.exists('xgboost_02.json'))
   b1 <- xgb.load('xgboost_01.json')
+  xgb.parameters(b1) <- list(nthread = 2)
   expect_equal(xgb.ntree(b1), 1)
   b2 <- xgb.load('xgboost_02.json')
+  xgb.parameters(b2) <- list(nthread = 2)
   expect_equal(xgb.ntree(b2), 2)
 
   xgb.config(b2) <- xgb.config(bst)
@@ -269,7 +273,8 @@ test_that("early stopping works with titanic", {
     objective = "binary:logistic",
     eval_metric = "auc",
     nrounds = 100,
-    early_stopping_rounds = 3
+    early_stopping_rounds = 3,
+    nthread = n_threads
   )
 
   expect_true(TRUE)  # should not crash
@@ -310,7 +315,7 @@ test_that("prediction in xgb.cv works", {
 
 test_that("prediction in xgb.cv works for gblinear too", {
   set.seed(11)
-  p <- list(booster = 'gblinear', objective = "reg:logistic", nthread = 2)
+  p <- list(booster = 'gblinear', objective = "reg:logistic", nthread = n_threads)
   cv <- xgb.cv(p, dtrain, nfold = 5, eta = 0.5, nrounds = 2, prediction = TRUE, verbose = 0)
   expect_false(is.null(cv$evaluation_log))
   expect_false(is.null(cv$pred))
@@ -343,7 +348,7 @@ test_that("prediction in xgb.cv for softprob works", {
   set.seed(11)
   expect_warning(
     cv <- xgb.cv(data = as.matrix(iris[, -5]), label = lb, nfold = 4,
-                 eta = 0.5, nrounds = 5, max_depth = 3, nthread = 2,
+                 eta = 0.5, nrounds = 5, max_depth = 3, nthread = n_threads,
                  subsample = 0.8, gamma = 2, verbose = 0,
                  prediction = TRUE, objective = "multi:softprob", num_class = 3)
   , NA)
diff --git a/R-package/tests/testthat/test_custom_objective.R b/R-package/tests/testthat/test_custom_objective.R
index d98e7045a723..b7cea8e0edfc 100644
--- a/R-package/tests/testthat/test_custom_objective.R
+++ b/R-package/tests/testthat/test_custom_objective.R
@@ -4,10 +4,16 @@ require(xgboost)
 
 set.seed(1994)
 
+n_threads <- 2
+
 data(agaricus.train, package = 'xgboost')
 data(agaricus.test, package = 'xgboost')
-dtrain <- xgb.DMatrix(agaricus.train$data, label = agaricus.train$label)
-dtest <- xgb.DMatrix(agaricus.test$data, label = agaricus.test$label)
+dtrain <- xgb.DMatrix(
+  agaricus.train$data, label = agaricus.train$label, nthread = n_threads
+)
+dtest <- xgb.DMatrix(
+  agaricus.test$data, label = agaricus.test$label, nthread = n_threads
+)
 watchlist <- list(eval = dtest, train = dtrain)
 
 logregobj <- function(preds, dtrain) {
@@ -24,7 +30,7 @@ evalerror <- function(preds, dtrain) {
   return(list(metric = "error", value = err))
 }
 
-param <- list(max_depth = 2, eta = 1, nthread = 2,
+param <- list(max_depth = 2, eta = 1, nthread = n_threads,
               objective = logregobj, eval_metric = evalerror)
 num_round <- 2
 
@@ -69,20 +75,20 @@ test_that("custom objective using DMatrix attr works", {
 test_that("custom objective with multi-class works", {
   data <- as.matrix(iris[, -5])
   label <-  as.numeric(iris$Species) - 1
-  dtrain <- xgb.DMatrix(data = data, label = label)
-  nclasses <- 3
+  dtrain <- xgb.DMatrix(data = data, label = label, nthread = n_threads)
+  n_classes <- 3
 
   fake_softprob <- function(preds, dtrain) {
     expect_true(all(matrix(preds) == 0.5))
     grad <- rnorm(dim(as.matrix(preds))[1])
-    expect_equal(dim(data)[1] * nclasses, dim(as.matrix(preds))[1])
+    expect_equal(dim(data)[1] * n_classes, dim(as.matrix(preds))[1])
     hess <- rnorm(dim(as.matrix(preds))[1])
     return (list(grad = grad, hess = hess))
   }
   fake_merror <- function(preds, dtrain) {
-    expect_equal(dim(data)[1] * nclasses, dim(as.matrix(preds))[1])
+    expect_equal(dim(data)[1] * n_classes, dim(as.matrix(preds))[1])
   }
   param$objective <- fake_softprob
   param$eval_metric <- fake_merror
-  bst <- xgb.train(param, dtrain, 1, num_class = nclasses)
+  bst <- xgb.train(param, dtrain, 1, num_class = n_classes)
 })
diff --git a/R-package/tests/testthat/test_dmatrix.R b/R-package/tests/testthat/test_dmatrix.R
index eb83544d8d16..c8f6ecd67097 100644
--- a/R-package/tests/testthat/test_dmatrix.R
+++ b/R-package/tests/testthat/test_dmatrix.R
@@ -7,25 +7,63 @@ data(agaricus.test, package = 'xgboost')
 test_data <- agaricus.test$data[1:100, ]
 test_label <- agaricus.test$label[1:100]
 
+n_threads <- 2
+
 test_that("xgb.DMatrix: basic construction", {
   # from sparse matrix
-  dtest1 <- xgb.DMatrix(test_data, label = test_label)
+  dtest1 <- xgb.DMatrix(test_data, label = test_label, nthread = n_threads)
 
   # from dense matrix
-  dtest2 <- xgb.DMatrix(as.matrix(test_data), label = test_label)
-  expect_equal(getinfo(dtest1, 'label'), getinfo(dtest2, 'label'))
+  dtest2 <- xgb.DMatrix(as.matrix(test_data), label = test_label, nthread = n_threads)
+  expect_equal(getinfo(dtest1, "label"), getinfo(dtest2, "label"))
   expect_equal(dim(dtest1), dim(dtest2))
 
   #from dense integer matrix
   int_data <- as.matrix(test_data)
   storage.mode(int_data) <- "integer"
-  dtest3 <- xgb.DMatrix(int_data, label = test_label)
+  dtest3 <- xgb.DMatrix(int_data, label = test_label, nthread = n_threads)
   expect_equal(dim(dtest1), dim(dtest3))
 })
 
+test_that("xgb.DMatrix: NA", {
+  n_samples <- 3
+  x <- cbind(
+    x1 = sample(x = 4, size = n_samples, replace = TRUE),
+    x2 = sample(x = 4, size = n_samples, replace = TRUE)
+  )
+  x[1, "x1"] <- NA
+
+  m <- xgb.DMatrix(x, nthread = n_threads)
+  xgb.DMatrix.save(m, "int.dmatrix")
+
+  x <- matrix(as.numeric(x), nrow = n_samples, ncol = 2)
+  colnames(x) <- c("x1", "x2")
+  m <- xgb.DMatrix(x, nthread = n_threads)
+
+  xgb.DMatrix.save(m, "float.dmatrix")
+
+  iconn <- file("int.dmatrix", "rb")
+  fconn <- file("float.dmatrix", "rb")
+
+  expect_equal(file.size("int.dmatrix"), file.size("float.dmatrix"))
+
+  bytes <- file.size("int.dmatrix")
+  idmatrix <- readBin(iconn, "raw", n = bytes)
+  fdmatrix <- readBin(fconn, "raw", n = bytes)
+
+  expect_equal(length(idmatrix), length(fdmatrix))
+  expect_equal(idmatrix, fdmatrix)
+
+  close(iconn)
+  close(fconn)
+
+  file.remove("int.dmatrix")
+  file.remove("float.dmatrix")
+})
+
 test_that("xgb.DMatrix: saving, loading", {
   # save to a local file
-  dtest1 <- xgb.DMatrix(test_data, label = test_label)
+  dtest1 <- xgb.DMatrix(test_data, label = test_label, nthread = n_threads)
   tmp_file <- tempfile('xgb.DMatrix_')
   on.exit(unlink(tmp_file))
   expect_true(xgb.DMatrix.save(dtest1, tmp_file))
@@ -39,13 +77,18 @@ test_that("xgb.DMatrix: saving, loading", {
   tmp <- c("0 1:1 2:1", "1 3:1", "0 1:1")
   tmp_file <- 'tmp.libsvm'
   writeLines(tmp, tmp_file)
-  dtest4 <- xgb.DMatrix(tmp_file, silent = TRUE)
+  expect_true(file.exists(tmp_file))
+  dtest4 <- xgb.DMatrix(
+    paste(tmp_file, "?format=libsvm", sep = ""), silent = TRUE, nthread = n_threads
+  )
   expect_equal(dim(dtest4), c(3, 4))
   expect_equal(getinfo(dtest4, 'label'), c(0, 1, 0))
 
   # check that feature info is saved
   data(agaricus.train, package = 'xgboost')
-  dtrain <- xgb.DMatrix(data = agaricus.train$data, label = agaricus.train$label)
+  dtrain <- xgb.DMatrix(
+    data = agaricus.train$data, label = agaricus.train$label, nthread = n_threads
+  )
   cnames <- colnames(dtrain)
   expect_equal(length(cnames), 126)
   tmp_file <- tempfile('xgb.DMatrix_')
@@ -59,7 +102,7 @@ test_that("xgb.DMatrix: saving, loading", {
 })
 
 test_that("xgb.DMatrix: getinfo & setinfo", {
-  dtest <- xgb.DMatrix(test_data)
+  dtest <- xgb.DMatrix(test_data, nthread = n_threads)
   expect_true(setinfo(dtest, 'label', test_label))
   labels <- getinfo(dtest, 'label')
   expect_equal(test_label, getinfo(dtest, 'label'))
@@ -86,7 +129,7 @@ test_that("xgb.DMatrix: getinfo & setinfo", {
 })
 
 test_that("xgb.DMatrix: slice, dim", {
-  dtest <- xgb.DMatrix(test_data, label = test_label)
+  dtest <- xgb.DMatrix(test_data, label = test_label, nthread = n_threads)
   expect_equal(dim(dtest), dim(test_data))
   dsub1 <- slice(dtest, 1:42)
   expect_equal(nrow(dsub1), 42)
@@ -101,16 +144,20 @@ test_that("xgb.DMatrix: slice, trailing empty rows", {
   data(agaricus.train, package = 'xgboost')
   train_data <- agaricus.train$data
   train_label <- agaricus.train$label
-  dtrain <- xgb.DMatrix(data = train_data, label = train_label)
+  dtrain <- xgb.DMatrix(
+    data = train_data, label = train_label, nthread = n_threads
+  )
   slice(dtrain, 6513L)
   train_data[6513, ] <- 0
-  dtrain <- xgb.DMatrix(data = train_data, label = train_label)
+  dtrain <- xgb.DMatrix(
+    data = train_data, label = train_label, nthread = n_threads
+  )
   slice(dtrain, 6513L)
   expect_equal(nrow(dtrain), 6513)
 })
 
 test_that("xgb.DMatrix: colnames", {
-  dtest <- xgb.DMatrix(test_data, label = test_label)
+  dtest <- xgb.DMatrix(test_data, label = test_label, nthread = n_threads)
   expect_equal(colnames(dtest), colnames(test_data))
   expect_error(colnames(dtest) <- 'asdf')
   new_names <- make.names(seq_len(ncol(test_data)))
@@ -126,6 +173,6 @@ test_that("xgb.DMatrix: nrow is correct for a very sparse matrix", {
   x <- rsparsematrix(nr, 100, density = 0.0005)
   # we want it very sparse, so that last rows are empty
   expect_lt(max(x@i), nr)
-  dtest <- xgb.DMatrix(x)
+  dtest <- xgb.DMatrix(x, nthread = n_threads)
   expect_equal(dim(dtest), dim(x))
 })
diff --git a/R-package/tests/testthat/test_feature_weights.R b/R-package/tests/testthat/test_feature_weights.R
index 580f58456dbb..2757e97ea21c 100644
--- a/R-package/tests/testthat/test_feature_weights.R
+++ b/R-package/tests/testthat/test_feature_weights.R
@@ -2,6 +2,8 @@ library(xgboost)
 
 context("feature weights")
 
+n_threads <- 2
+
 test_that("training with feature weights works", {
   nrows <- 1000
   ncols <- 9
@@ -12,8 +14,12 @@ test_that("training with feature weights works", {
 
   test <- function(tm) {
     names <- paste0("f", 1:ncols)
-    xy <- xgb.DMatrix(data = x, label = y, feature_weights = weights)
-    params <- list(colsample_bynode = 0.4, tree_method = tm, nthread = 1)
+    xy <- xgb.DMatrix(
+      data = x, label = y, feature_weights = weights, nthread = n_threads
+    )
+    params <- list(
+      colsample_bynode = 0.4, tree_method = tm, nthread = n_threads
+    )
     model <- xgb.train(params = params, data = xy, nrounds = 32)
     importance <- xgb.importance(model = model, feature_names = names)
     expect_equal(dim(importance), c(ncols, 4))
diff --git a/R-package/tests/testthat/test_glm.R b/R-package/tests/testthat/test_glm.R
index 270267a0fe0c..9e0a3551f68e 100644
--- a/R-package/tests/testthat/test_glm.R
+++ b/R-package/tests/testthat/test_glm.R
@@ -1,15 +1,19 @@
 context('Test generalized linear models')
 
-require(xgboost)
+n_threads <- 2
 
 test_that("gblinear works", {
   data(agaricus.train, package = 'xgboost')
   data(agaricus.test, package = 'xgboost')
-  dtrain <- xgb.DMatrix(agaricus.train$data, label = agaricus.train$label)
-  dtest <- xgb.DMatrix(agaricus.test$data, label = agaricus.test$label)
+  dtrain <- xgb.DMatrix(
+    agaricus.train$data, label = agaricus.train$label, nthread = n_threads
+  )
+  dtest <- xgb.DMatrix(
+    agaricus.test$data, label = agaricus.test$label, nthread = n_threads
+  )
 
   param <- list(objective = "binary:logistic", eval_metric = "error", booster = "gblinear",
-                nthread = 2, eta = 0.8, alpha = 0.0001, lambda = 0.0001)
+                nthread = n_threads, eta = 0.8, alpha = 0.0001, lambda = 0.0001)
   watchlist <- list(eval = dtest, train = dtrain)
 
   n <- 5         # iterations
@@ -50,12 +54,16 @@ test_that("gblinear works", {
 test_that("gblinear early stopping works", {
   data(agaricus.train, package = 'xgboost')
   data(agaricus.test, package = 'xgboost')
-  dtrain <- xgb.DMatrix(agaricus.train$data, label = agaricus.train$label)
-  dtest <- xgb.DMatrix(agaricus.test$data, label = agaricus.test$label)
+  dtrain <- xgb.DMatrix(
+    agaricus.train$data, label = agaricus.train$label, nthread = n_threads
+  )
+  dtest <- xgb.DMatrix(
+    agaricus.test$data, label = agaricus.test$label, nthread = n_threads
+  )
 
   param <- list(
     objective = "binary:logistic", eval_metric = "error", booster = "gblinear",
-    nthread = 2, eta = 0.8, alpha = 0.0001, lambda = 0.0001,
+    nthread = n_threads, eta = 0.8, alpha = 0.0001, lambda = 0.0001,
     updater = "coord_descent"
   )
 
diff --git a/R-package/tests/testthat/test_helpers.R b/R-package/tests/testthat/test_helpers.R
index fdd0ce02b5a6..7f559753a4dc 100644
--- a/R-package/tests/testthat/test_helpers.R
+++ b/R-package/tests/testthat/test_helpers.R
@@ -161,6 +161,7 @@ test_that("SHAPs sum to predictions, with or without DART", {
     fit <- xgboost(
       params = c(
         list(
+          nthread = 2,
           booster = booster,
           objective = "reg:squarederror",
           eval_metric = "rmse"),
@@ -243,7 +244,7 @@ if (grepl('Windows', Sys.info()[['sysname']]) ||
 test_that("xgb.Booster serializing as R object works", {
   saveRDS(bst.Tree, 'xgb.model.rds')
   bst <- readRDS('xgb.model.rds')
-  dtrain <- xgb.DMatrix(sparse_matrix, label = label)
+  dtrain <- xgb.DMatrix(sparse_matrix, label = label, nthread = 2)
   expect_equal(predict(bst.Tree, dtrain), predict(bst, dtrain), tolerance = float_tolerance)
   expect_equal(xgb.dump(bst.Tree), xgb.dump(bst))
   xgb.save(bst, 'xgb.model')
@@ -345,7 +346,8 @@ test_that("xgb.importance works with and without feature names", {
   m <- xgboost::xgboost(
     data = as.matrix(data.frame(x = c(0, 1))),
     label = c(1, 2),
-    nrounds = 1
+    nrounds = 1,
+    nthread = 2
   )
   df <- xgb.model.dt.tree(model = m)
   expect_equal(df$Feature, "Leaf")
diff --git a/R-package/tests/testthat/test_interaction_constraints.R b/R-package/tests/testthat/test_interaction_constraints.R
index 7f6a8b09bd67..f3d6aa7b3844 100644
--- a/R-package/tests/testthat/test_interaction_constraints.R
+++ b/R-package/tests/testthat/test_interaction_constraints.R
@@ -2,6 +2,8 @@ require(xgboost)
 
 context("interaction constraints")
 
+n_threads <- 2
+
 set.seed(1024)
 x1 <- rnorm(1000, 1)
 x2 <- rnorm(1000, 1)
@@ -45,11 +47,18 @@ test_that("interaction constraints scientific representation", {
   d <- matrix(rexp(rows, rate = .1), nrow = rows, ncol = cols)
   y <- rnorm(rows)
 
-  dtrain <- xgb.DMatrix(data = d, info = list(label = y))
+  dtrain <- xgb.DMatrix(data = d, info = list(label = y), nthread = n_threads)
   inc <- list(c(seq.int(from = 0, to = cols, by = 1)))
 
-  with_inc <- xgb.train(data = dtrain, tree_method = 'hist',
-                        interaction_constraints = inc, nrounds = 10)
-  without_inc <- xgb.train(data = dtrain, tree_method = 'hist', nrounds = 10)
+  with_inc <- xgb.train(
+    data = dtrain,
+    tree_method = 'hist',
+    interaction_constraints = inc,
+    nrounds = 10,
+    nthread = n_threads
+  )
+  without_inc <- xgb.train(
+    data = dtrain, tree_method = 'hist', nrounds = 10, nthread = n_threads
+  )
   expect_equal(xgb.save.raw(with_inc), xgb.save.raw(without_inc))
 })
diff --git a/R-package/tests/testthat/test_interactions.R b/R-package/tests/testthat/test_interactions.R
index e90467cdcf62..ba14b1973553 100644
--- a/R-package/tests/testthat/test_interactions.R
+++ b/R-package/tests/testthat/test_interactions.R
@@ -3,6 +3,7 @@ context('Test prediction of feature interactions')
 require(xgboost)
 
 set.seed(123)
+n_threads <- 2
 
 test_that("predict feature interactions works", {
   # simulate some binary data and a linear outcome with an interaction term
@@ -21,8 +22,10 @@ test_that("predict feature interactions works", {
 
   y <- f_int(X)
 
-  dm <- xgb.DMatrix(X, label = y)
-  param <- list(eta = 0.1, max_depth = 4, base_score = mean(y), lambda = 0, nthread = 2)
+  dm <- xgb.DMatrix(X, label = y, nthread = n_threads)
+  param <- list(
+    eta = 0.1, max_depth = 4, base_score = mean(y), lambda = 0, nthread = n_threads
+  )
   b <- xgb.train(param, dm, 100)
 
   pred <- predict(b, dm, outputmargin = TRUE)
@@ -101,11 +104,13 @@ test_that("SHAP contribution values are not NAN", {
     verbose = 0,
     params = list(
       objective = "reg:squarederror",
-      eval_metric = "rmse"),
+      eval_metric = "rmse",
+      nthread = n_threads
+    ),
     data = as.matrix(subset(d, fold == 2)[, ivs]),
     label = subset(d, fold == 2)$y,
-    nthread = 1,
-    nrounds = 3)
+    nrounds = 3
+  )
 
   shaps <- as.data.frame(predict(fit,
     newdata = as.matrix(subset(d, fold == 1)[, ivs]),
@@ -118,8 +123,12 @@ test_that("SHAP contribution values are not NAN", {
 
 
 test_that("multiclass feature interactions work", {
-  dm <- xgb.DMatrix(as.matrix(iris[, -5]), label = as.numeric(iris$Species) - 1)
-  param <- list(eta = 0.1, max_depth = 4, objective = 'multi:softprob', num_class = 3)
+  dm <- xgb.DMatrix(
+    as.matrix(iris[, -5]), label = as.numeric(iris$Species) - 1, nthread = n_threads
+  )
+  param <- list(
+    eta = 0.1, max_depth = 4, objective = 'multi:softprob', num_class = 3, nthread = n_threads
+  )
   b <- xgb.train(param, dm, 40)
   pred <- t(
     array(
@@ -168,6 +177,7 @@ test_that("SHAP single sample works", {
     max_depth = 2,
     nrounds = 4,
     objective = "binary:logistic",
+    nthread = n_threads
   )
 
   predt <- predict(
diff --git a/R-package/tests/testthat/test_io.R b/R-package/tests/testthat/test_io.R
index 5b2bc4265a5f..0da40ba713b3 100644
--- a/R-package/tests/testthat/test_io.R
+++ b/R-package/tests/testthat/test_io.R
@@ -12,7 +12,8 @@ test_that("load/save raw works", {
   nrounds <- 8
   booster <- xgboost(
     data = train$data, label = train$label,
-    nrounds = nrounds, objective = "binary:logistic"
+    nrounds = nrounds, objective = "binary:logistic",
+    nthread = 2
   )
 
   json_bytes <- xgb.save.raw(booster, raw_format = "json")
diff --git a/R-package/tests/testthat/test_model_compatibility.R b/R-package/tests/testthat/test_model_compatibility.R
index 6a35ca4b31c6..5e74a5561c4f 100644
--- a/R-package/tests/testthat/test_model_compatibility.R
+++ b/R-package/tests/testthat/test_model_compatibility.R
@@ -68,7 +68,7 @@ test_that("Models from previous versions of XGBoost can be loaded", {
                 destfile = zipfile, mode = 'wb', quiet = TRUE)
   unzip(zipfile, overwrite = TRUE)
 
-  pred_data <- xgb.DMatrix(matrix(c(0, 0, 0, 0), nrow = 1, ncol = 4))
+  pred_data <- xgb.DMatrix(matrix(c(0, 0, 0, 0), nrow = 1, ncol = 4), nthread = 2)
 
   lapply(list.files(model_dir), function (x) {
     model_file <- file.path(model_dir, x)
diff --git a/R-package/tests/testthat/test_parameter_exposure.R b/R-package/tests/testthat/test_parameter_exposure.R
index 86413174b853..4e7f5333fb62 100644
--- a/R-package/tests/testthat/test_parameter_exposure.R
+++ b/R-package/tests/testthat/test_parameter_exposure.R
@@ -5,8 +5,12 @@ require(xgboost)
 data(agaricus.train, package = 'xgboost')
 data(agaricus.test, package = 'xgboost')
 
-dtrain <- xgb.DMatrix(agaricus.train$data, label = agaricus.train$label)
-dtest <- xgb.DMatrix(agaricus.test$data, label = agaricus.test$label)
+dtrain <- xgb.DMatrix(
+  agaricus.train$data, label = agaricus.train$label, nthread = 2
+)
+dtest <- xgb.DMatrix(
+  agaricus.test$data, label = agaricus.test$label, nthread = 2
+)
 
 bst <- xgboost(data = dtrain,
                max_depth = 2,
diff --git a/R-package/tests/testthat/test_poisson_regression.R b/R-package/tests/testthat/test_poisson_regression.R
index 4f3527cdb31f..214d1f218f94 100644
--- a/R-package/tests/testthat/test_poisson_regression.R
+++ b/R-package/tests/testthat/test_poisson_regression.R
@@ -5,8 +5,10 @@ set.seed(1994)
 
 test_that("Poisson regression works", {
   data(mtcars)
-  bst <- xgboost(data = as.matrix(mtcars[, -11]), label = mtcars[, 11],
-                objective = 'count:poisson', nrounds = 10, verbose = 0)
+  bst <- xgboost(
+    data = as.matrix(mtcars[, -11]), label = mtcars[, 11],
+    objective = 'count:poisson', nrounds = 10, verbose = 0, nthread = 2
+  )
   expect_equal(class(bst), "xgb.Booster")
   pred <- predict(bst, as.matrix(mtcars[, -11]))
   expect_equal(length(pred), 32)
diff --git a/R-package/tests/testthat/test_ranking.R b/R-package/tests/testthat/test_ranking.R
index 7a352bea2ec0..a9afc65b9e76 100644
--- a/R-package/tests/testthat/test_ranking.R
+++ b/R-package/tests/testthat/test_ranking.R
@@ -3,16 +3,18 @@ require(Matrix)
 
 context('Learning to rank')
 
+n_threads <- 2
+
 test_that('Test ranking with unweighted data', {
   X <- sparseMatrix(i = c(2, 3, 7, 9, 12, 15, 17, 18),
                     j = c(1, 1, 2, 2,  3,  3,  4,  4),
                     x = rep(1.0, 8), dims = c(20, 4))
   y <- c(0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0)
   group <- c(5, 5, 5, 5)
-  dtrain <- xgb.DMatrix(X, label = y, group = group)
+  dtrain <- xgb.DMatrix(X, label = y, group = group, nthread = n_threads)
 
   params <- list(eta = 1, tree_method = 'exact', objective = 'rank:pairwise', max_depth = 1,
-                 eval_metric = 'auc', eval_metric = 'aucpr')
+                 eval_metric = 'auc', eval_metric = 'aucpr', nthread = n_threads)
   bst <- xgb.train(params, dtrain, nrounds = 10, watchlist = list(train = dtrain))
   # Check if the metric is monotone increasing
   expect_true(all(diff(bst$evaluation_log$train_auc) >= 0))
@@ -26,10 +28,14 @@ test_that('Test ranking with weighted data', {
   y <- c(0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0)
   group <- c(5, 5, 5, 5)
   weight <- c(1.0, 2.0, 3.0, 4.0)
-  dtrain <- xgb.DMatrix(X, label = y, group = group, weight = weight)
+  dtrain <- xgb.DMatrix(
+    X, label = y, group = group, weight = weight, nthread = n_threads
+  )
 
-  params <- list(eta = 1, tree_method = 'exact', objective = 'rank:pairwise', max_depth = 1,
-                 eval_metric = 'auc', eval_metric = 'aucpr')
+  params <- list(
+    eta = 1, tree_method = "exact", objective = "rank:pairwise", max_depth = 1,
+    eval_metric = "auc", eval_metric = "aucpr", nthread = n_threads
+  )
   bst <- xgb.train(params, dtrain, nrounds = 10, watchlist = list(train = dtrain))
   # Check if the metric is monotone increasing
   expect_true(all(diff(bst$evaluation_log$train_auc) >= 0))
diff --git a/R-package/tests/testthat/test_unicode.R b/R-package/tests/testthat/test_unicode.R
new file mode 100644
index 000000000000..c8a225716f81
--- /dev/null
+++ b/R-package/tests/testthat/test_unicode.R
@@ -0,0 +1,22 @@
+context("Test Unicode handling")
+
+data(agaricus.train, package = 'xgboost')
+data(agaricus.test, package = 'xgboost')
+train <- agaricus.train
+test <- agaricus.test
+set.seed(1994)
+
+test_that("Can save and load models with Unicode paths", {
+  nrounds <- 2
+  bst <- xgboost(data = train$data, label = train$label, max_depth = 2,
+                 eta = 1, nthread = 2, nrounds = nrounds, objective = "binary:logistic",
+                 eval_metric = "error")
+  tmpdir <- tempdir()
+  lapply(c("모델.json", "がうる・ぐら.json", "类继承.ubj"), function(x) {
+    path <- file.path(tmpdir, x)
+    xgb.save(bst, path)
+    bst2 <- xgb.load(path)
+    xgb.parameters(bst2) <- list(nthread = 2)
+    expect_equal(predict(bst, test$data), predict(bst2, test$data))
+  })
+})
diff --git a/R-package/tests/testthat/test_update.R b/R-package/tests/testthat/test_update.R
index 541fdf68e59f..3b83846251fd 100644
--- a/R-package/tests/testthat/test_update.R
+++ b/R-package/tests/testthat/test_update.R
@@ -4,8 +4,15 @@ context("update trees in an existing model")
 
 data(agaricus.train, package = 'xgboost')
 data(agaricus.test, package = 'xgboost')
-dtrain <- xgb.DMatrix(agaricus.train$data, label = agaricus.train$label)
-dtest <- xgb.DMatrix(agaricus.test$data, label = agaricus.test$label)
+
+n_threads <- 1
+
+dtrain <- xgb.DMatrix(
+  agaricus.train$data, label = agaricus.train$label, nthread = n_threads
+)
+dtest <- xgb.DMatrix(
+  agaricus.test$data, label = agaricus.test$label, nthread = n_threads
+)
 
 # Disable flaky tests for 32-bit Windows.
 # See https://github.com/dmlc/xgboost/issues/3720
@@ -15,7 +22,10 @@ test_that("updating the model works", {
   watchlist <- list(train = dtrain, test = dtest)
 
   # no-subsampling
-  p1 <- list(objective = "binary:logistic", max_depth = 2, eta = 0.05, nthread = 2)
+  p1 <- list(
+    objective = "binary:logistic", max_depth = 2, eta = 0.05, nthread = n_threads,
+    updater = "grow_colmaker,prune"
+  )
   set.seed(11)
   bst1 <- xgb.train(p1, dtrain, nrounds = 10, watchlist, verbose = 0)
   tr1 <- xgb.model.dt.tree(model = bst1)
@@ -85,9 +95,11 @@ test_that("updating the model works", {
 })
 
 test_that("updating works for multiclass & multitree", {
-  dtr <- xgb.DMatrix(as.matrix(iris[, -5]), label = as.numeric(iris$Species) - 1)
+  dtr <- xgb.DMatrix(
+    as.matrix(iris[, -5]), label = as.numeric(iris$Species) - 1, nthread = n_threads
+  )
   watchlist <- list(train = dtr)
-  p0 <- list(max_depth = 2, eta = 0.5, nthread = 2, subsample = 0.6,
+  p0 <- list(max_depth = 2, eta = 0.5, nthread = n_threads, subsample = 0.6,
              objective = "multi:softprob", num_class = 3, num_parallel_tree = 2,
              base_score = 0)
   set.seed(121)
diff --git a/R-package/vignettes/discoverYourData.Rmd b/R-package/vignettes/discoverYourData.Rmd
index 6706be81ab14..acf5b5c1ea46 100644
--- a/R-package/vignettes/discoverYourData.Rmd
+++ b/R-package/vignettes/discoverYourData.Rmd
@@ -28,7 +28,11 @@ Package loading:
 require(xgboost)
 require(Matrix)
 require(data.table)
-if (!require('vcd')) install.packages('vcd')
+if (!require('vcd')) {
+  install.packages('vcd')
+}
+
+data.table::setDTthreads(2)
 ```
 
 > **VCD** package is used for one of its embedded dataset only.
@@ -327,10 +331,27 @@ train <- agaricus.train
 test <- agaricus.test
 
 #Random Forest - 1000 trees
-bst <- xgboost(data = train$data, label = train$label, max_depth = 4, num_parallel_tree = 1000, subsample = 0.5, colsample_bytree =0.5, nrounds = 1, objective = "binary:logistic")
+bst <- xgboost(
+    data = train$data,
+    label = train$label,
+    max_depth = 4,
+    num_parallel_tree = 1000,
+    subsample = 0.5,
+    colsample_bytree = 0.5,
+    nrounds = 1,
+    objective = "binary:logistic",
+    nthread = 2
+)
 
 #Boosting - 3 rounds
-bst <- xgboost(data = train$data, label = train$label, max_depth = 4, nrounds = 3, objective = "binary:logistic")
+bst <- xgboost(
+    data = train$data,
+    label = train$label,
+    max_depth = 4,
+    nrounds = 3,
+    objective = "binary:logistic",
+    nthread = 2
+)
 ```
 
 > Note that the parameter `round` is set to `1`.
diff --git a/R-package/vignettes/xgboost.Rnw b/R-package/vignettes/xgboost.Rnw
index c9089cd6f23b..7edf4ace3d4f 100644
--- a/R-package/vignettes/xgboost.Rnw
+++ b/R-package/vignettes/xgboost.Rnw
@@ -86,9 +86,10 @@ data(agaricus.test, package='xgboost')
 train <- agaricus.train
 test <- agaricus.test
 bst <- xgboost(data = train$data, label = train$label, max_depth = 2, eta = 1,
-               nrounds = 2, objective = "binary:logistic")
+               nrounds = 2, objective = "binary:logistic", nthread = 2)
 xgb.save(bst, 'model.save')
 bst = xgb.load('model.save')
+xgb.parameters(bst) <- list(nthread = 2)
 pred <- predict(bst, test$data)
 @
 
@@ -127,7 +128,7 @@ training from initial prediction value, weighted training instance.
 
 We can use \verb@xgb.DMatrix@ to construct an \verb@xgb.DMatrix@ object:
 <<xgb.DMatrix>>=
-dtrain <- xgb.DMatrix(train$data, label = train$label)
+dtrain <- xgb.DMatrix(train$data, label = train$label, nthread = 2)
 class(dtrain)
 head(getinfo(dtrain,'label'))
 @
@@ -161,9 +162,9 @@ evalerror <- function(preds, dtrain) {
   return(list(metric = "MSE", value = err))
 }
 
-dtest <- xgb.DMatrix(test$data, label = test$label)
+dtest <- xgb.DMatrix(test$data, label = test$label, nthread = 2)
 watchlist <- list(eval = dtest, train = dtrain)
-param <- list(max_depth = 2, eta = 1)
+param <- list(max_depth = 2, eta = 1, nthread = 2)
 
 bst <- xgb.train(param, dtrain, nrounds = 2, watchlist, logregobj, evalerror, maximize = FALSE)
 @
diff --git a/R-package/vignettes/xgboostPresentation.Rmd b/R-package/vignettes/xgboostPresentation.Rmd
index 218b12eeb103..cd389c7f18b4 100644
--- a/R-package/vignettes/xgboostPresentation.Rmd
+++ b/R-package/vignettes/xgboostPresentation.Rmd
@@ -164,7 +164,15 @@ bstSparse <- xgboost(data = train$data, label = train$label, max_depth = 2, eta
 Alternatively, you can put your dataset in a *dense* matrix, i.e. a basic **R** matrix.
 
 ```{r trainingDense, message=F, warning=F}
-bstDense <- xgboost(data = as.matrix(train$data), label = train$label, max_depth = 2, eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic")
+bstDense <- xgboost(
+    data = as.matrix(train$data),
+    label = train$label,
+    max_depth = 2,
+    eta = 1,
+    nthread = 2,
+    nrounds = 2,
+    objective = "binary:logistic"
+)
 ```
 
 ##### xgb.DMatrix
@@ -172,8 +180,15 @@ bstDense <- xgboost(data = as.matrix(train$data), label = train$label, max_depth
 **XGBoost** offers a way to group them in a `xgb.DMatrix`. You can even add other meta data in it. It will be useful for the most advanced features we will discover later.
 
 ```{r trainingDmatrix, message=F, warning=F}
-dtrain <- xgb.DMatrix(data = train$data, label = train$label)
-bstDMatrix <- xgboost(data = dtrain, max_depth = 2, eta = 1, nthread = 2, nrounds = 2, objective = "binary:logistic")
+dtrain <- xgb.DMatrix(data = train$data, label = train$label, nthread = 2)
+bstDMatrix <- xgboost(
+    data = dtrain,
+    max_depth = 2,
+    eta = 1,
+    nthread = 2,
+    nrounds = 2,
+    objective = "binary:logistic"
+)
 ```
 
 ##### Verbose option
@@ -267,8 +282,8 @@ Most of the features below have been implemented to help you to improve your mod
 For the following advanced features, we need to put data in `xgb.DMatrix` as explained above.
 
 ```{r DMatrix, message=F, warning=F}
-dtrain <- xgb.DMatrix(data = train$data, label=train$label)
-dtest <- xgb.DMatrix(data = test$data, label=test$label)
+dtrain <- xgb.DMatrix(data = train$data, label = train$label, nthread = 2)
+dtest <- xgb.DMatrix(data = test$data, label = test$label, nthread = 2)
 ```
 
 ### Measure learning progress with xgb.train
@@ -393,6 +408,7 @@ An interesting test to see how identical our saved model is to the original one
 ```{r loadModel, message=F, warning=F}
 # load binary model to R
 bst2 <- xgb.load("xgboost.model")
+xgb.parameters(bst2) <- list(nthread = 2)
 pred2 <- predict(bst2, test$data)
 
 # And now the test
@@ -417,6 +433,7 @@ print(class(rawVec))
 
 # load binary model to R
 bst3 <- xgb.load(rawVec)
+xgb.parameters(bst3) <- list(nthread = 2)
 pred3 <- predict(bst3, test$data)
 
 # pred2 should be identical to pred
diff --git a/R-package/vignettes/xgboostfromJSON.Rmd b/R-package/vignettes/xgboostfromJSON.Rmd
index 544186830d11..0946a45d5ee1 100644
--- a/R-package/vignettes/xgboostfromJSON.Rmd
+++ b/R-package/vignettes/xgboostfromJSON.Rmd
@@ -174,7 +174,7 @@ bst_preds == bst_from_json_preds
 
 None are exactly equal again.  What is going on here?  Well, since we are using the value `1` in the calculations, we have introduced a double into the calculation.  Because of this, all float values are promoted to 64-bit doubles and the 64-bit version of the exponentiation operator `exp` is also used.  On the other hand, xgboost uses the 32-bit version of the exponentiation operator in its [sigmoid function](https://github.com/dmlc/xgboost/blob/54980b8959680a0da06a3fc0ec776e47c8cbb0a1/src/common/math.h#L25-L27).
 
-How do we fix this?  We have to ensure we use the correct data types everywhere and the correct operators.  If we use only floats, the float library that we have loaded will ensure the 32-bit float exponentiation operator is applied. 
+How do we fix this?  We have to ensure we use the correct data types everywhere and the correct operators.  If we use only floats, the float library that we have loaded will ensure the 32-bit float exponentiation operator is applied.
 ```{r}
 # calculate the predictions casting doubles to floats
 bst_from_json_preds <- ifelse(fl(data$dates)<fl(node$split_condition),
diff --git a/src/common/common.h b/src/common/common.h
index b2d7211c6932..d83b58d6ce98 100644
--- a/src/common/common.h
+++ b/src/common/common.h
@@ -256,14 +256,17 @@ inline void SetDevice(std::int32_t device) {
 }
 #endif
 
-template <typename Idx, typename Container,
-          typename V = typename Container::value_type,
+template <typename Idx, typename Container, typename V = typename Container::value_type,
           typename Comp = std::less<V>>
 std::vector<Idx> ArgSort(Container const &array, Comp comp = std::less<V>{}) {
   std::vector<Idx> result(array.size());
   std::iota(result.begin(), result.end(), 0);
   auto op = [&array, comp](Idx const &l, Idx const &r) { return comp(array[l], array[r]); };
+#if defined(XGBOOST_STRICT_R_MODE) && XGBOOST_STRICT_R_MODE == 1
+  std::stable_sort(result.begin(), result.end(), op);
+#else
   XGBOOST_PARALLEL_STABLE_SORT(result.begin(), result.end(), op);
+#endif  // defined(XGBOOST_STRICT_R_MODE) && XGBOOST_STRICT_R_MODE == 1
   return result;
 }
 
diff --git a/src/gbm/gbtree_model.cc b/src/gbm/gbtree_model.cc
index 4e9cc6655eaa..6c8f0a20dee0 100644
--- a/src/gbm/gbtree_model.cc
+++ b/src/gbm/gbtree_model.cc
@@ -63,14 +63,30 @@ void GBTreeModel::Load(dmlc::Stream* fi) {
   }
 }
 
+namespace {
+std::int32_t IOThreads(Context const* ctx) {
+  CHECK(ctx);
+  std::int32_t n_threads = ctx->Threads();
+  // CRAN checks for number of threads used by examples, but we might not have the right
+  // number of threads when serializing/unserializing models as nthread is a booster
+  // parameter, which is only effective after booster initialization.
+  //
+  // The threshold ratio of CPU time to user time for R is 2.5, we set the number of
+  // threads to 2.
+#if defined(XGBOOST_STRICT_R_MODE) && XGBOOST_STRICT_R_MODE == 1
+  n_threads = std::min(2, n_threads);
+#endif
+  return n_threads;
+}
+}  // namespace
+
 void GBTreeModel::SaveModel(Json* p_out) const {
   auto& out = *p_out;
   CHECK_EQ(param.num_trees, static_cast<int>(trees.size()));
   out["gbtree_model_param"] = ToJson(param);
   std::vector<Json> trees_json(trees.size());
 
-  CHECK(ctx_);
-  common::ParallelFor(trees.size(), ctx_->Threads(), [&](auto t) {
+  common::ParallelFor(trees.size(), IOThreads(ctx_), [&](auto t) {
     auto const& tree = trees[t];
     Json tree_json{Object()};
     tree->SaveModel(&tree_json);
@@ -96,11 +112,10 @@ void GBTreeModel::LoadModel(Json const& in) {
   auto const& trees_json = get<Array const>(in["trees"]);
   trees.resize(trees_json.size());
 
-  CHECK(ctx_);
-  common::ParallelFor(trees_json.size(), ctx_->Threads(), [&](auto t) {
-    auto tree_id = get<Integer>(trees_json[t]["id"]);
-    trees.at(tree_id).reset(new RegTree());
-    trees.at(tree_id)->LoadModel(trees_json[t]);
+  common::ParallelFor(param.num_trees, IOThreads(ctx_), [&](auto t) {
+    auto tree_id = get<Integer const>(trees_json[t]["id"]);
+    trees.at(tree_id).reset(new RegTree{});
+    trees[tree_id]->LoadModel(trees_json[t]);
   });
 
   tree_info.resize(param.num_trees);
diff --git a/tests/ci_build/test_r_package.py b/tests/ci_build/test_r_package.py
index a720e7c60a99..853bf05022b1 100644
--- a/tests/ci_build/test_r_package.py
+++ b/tests/ci_build/test_r_package.py
@@ -1,91 +1,358 @@
+"""Utilities for packaging R code and running tests."""
 import argparse
 import os
+import shutil
 import subprocess
-from test_utils import DirectoryExcursion
+from io import StringIO
+from pathlib import Path
+from platform import system
 
-ROOT = os.path.normpath(
-    os.path.join(os.path.dirname(os.path.abspath(__file__)), os.path.pardir,
-                 os.path.pardir))
-r_package = os.path.join(ROOT, 'R-package')
+try:
+    import pandas as pd
+except ImportError:
+    pd = None
 
+from test_utils import R_PACKAGE, ROOT, DirectoryExcursion, cd, print_time, record_time
 
-def get_mingw_bin():
-    return os.path.join('c:/rtools40/mingw64/', 'bin')
 
+def get_mingw_bin() -> str:
+    return os.path.join("c:/rtools40/mingw64/", "bin")
 
-def test_with_autotools(args):
-    with DirectoryExcursion(r_package):
+
+@cd(ROOT)
+@record_time
+def pack_rpackage() -> Path:
+    """Compose the directory used for creating R package tar ball."""
+    dest = Path("xgboost")
+
+    def pkgroot(path: str) -> None:
+        """Change makefiles according to the package layout."""
+        with open(Path("R-package") / "src" / path, "r") as fd:
+            makefile = fd.read()
+            makefile = makefile.replace("PKGROOT=../../", "PKGROOT=.", 1)
+        with open(dest / "src" / path, "w") as fd:
+            fd.write(makefile)
+
+    output = subprocess.run(["git", "clean", "-xdf", "--dry-run"], capture_output=True)
+    if output.returncode != 0:
+        raise ValueError("Failed to check git repository status.", output)
+    would_remove = output.stdout.decode("utf-8").strip().split("\n")
+
+    if would_remove and not all(f.find("tests/ci_build") != -1 for f in would_remove):
+        raise ValueError(
+            "\n".join(would_remove) + "\nPlease cleanup the working git repository."
+        )
+
+    shutil.copytree("R-package", dest)
+    os.remove(dest / "demo" / "runall.R")
+    # core
+    shutil.copytree("src", dest / "src" / "src")
+    shutil.copytree("include", dest / "src" / "include")
+    shutil.copytree("amalgamation", dest / "src" / "amalgamation")
+    # rabit
+    rabit = Path("rabit")
+    os.mkdir(dest / "src" / rabit)
+    shutil.copytree(rabit / "src", dest / "src" / "rabit" / "src")
+    shutil.copytree(rabit / "include", dest / "src" / "rabit" / "include")
+    # dmlc-core
+    dmlc_core = Path("dmlc-core")
+    os.mkdir(dest / "src" / dmlc_core)
+    shutil.copytree(dmlc_core / "include", dest / "src" / "dmlc-core" / "include")
+    shutil.copytree(dmlc_core / "src", dest / "src" / "dmlc-core" / "src")
+    # makefile & license
+    shutil.copyfile("LICENSE", dest / "LICENSE")
+    osxmakef = dest / "src" / "Makevars.win-e"
+    if os.path.exists(osxmakef):
+        os.remove(osxmakef)
+    pkgroot("Makevars.in")
+    pkgroot("Makevars.win")
+    # misc
+    rwsp = Path("R-package") / "remove_warning_suppression_pragma.sh"
+    if system() != "Windows":
+        subprocess.check_call(rwsp)
+    rwsp = dest / "remove_warning_suppression_pragma.sh"
+    if system() != "Windows":
+        subprocess.check_call(rwsp)
+    os.remove(rwsp)
+    os.remove(dest / "CMakeLists.txt")
+    shutil.rmtree(dest / "tests" / "helper_scripts")
+    return dest
+
+
+@cd(ROOT)
+@record_time
+def build_rpackage(path: str) -> str:
+    def find_tarball() -> str:
+        found = []
+        for root, subdir, files in os.walk("."):
+            for f in files:
+                if f.endswith(".tar.gz") and f.startswith("xgboost"):
+                    found.append(os.path.join(root, f))
+        if not found:
+            raise ValueError("Failed to find output tar ball.")
+        if len(found) > 1:
+            raise ValueError("Found more than one packages:", found)
+        return found[0]
+
+    env = os.environ.copy()
+    print("Ncpus:", f"{os.cpu_count()}")
+    env.update({"MAKEFLAGS": f"-j{os.cpu_count()}"})
+    subprocess.check_call([R, "CMD", "build", path], env=env)
+
+    tarball = find_tarball()
+    return tarball
+
+
+def check_example_timing(rcheck_dir: Path, threshold: float) -> None:
+    with open(rcheck_dir / "xgboost-Ex.timings", "r") as fd:
+        timings = fd.readlines()
+        newlines = []
+        for line in timings:
+            line = line.strip()
+            newlines.append(line)
+        con_timings = "\n".join(newlines)
+        df = pd.read_csv(StringIO(con_timings), delimiter="\t")
+        ratio_n = "user/elapsed"
+        df[ratio_n] = df["user"] / df["elapsed"]
+        offending = df[df[ratio_n] > threshold]
+
+    try:
+        # requires the tabulate package
+        df.to_markdown("timings.md")
+        offending.to_markdown("offending.md")
+    except ImportError:
+        print("failed to export markdown files.")
+        pass
+
+    if offending.shape[0] == 0:
+        return
+
+    print(offending)
+    raise ValueError("There are examples using too many threads")
+
+
+@cd(ROOT)
+@record_time
+def check_rpackage(path: str) -> None:
+    env = os.environ.copy()
+    print("Ncpus:", f"{os.cpu_count()}")
+    threshold = 2.5
+    env.update(
+        {
+            "MAKEFLAGS": f"-j{os.cpu_count()}",
+            # cran specific environment variables
+            "_R_CHECK_EXAMPLE_TIMING_CPU_TO_ELAPSED_THRESHOLD_": str(threshold),
+            "_R_CHECK_TEST_TIMING_CPU_TO_ELAPSED_THRESHOLD_": str(threshold),
+            "_R_CHECK_VIGNETTE_TIMING_CPU_TO_ELAPSED_THRESHOLD_": str(threshold),
+        }
+    )
+
+    # Actually we don't run this check on windows due to dependency issue.
+    if system() == "Windows":
+        # make sure compiler from rtools is used.
         mingw_bin = get_mingw_bin()
-        CXX = os.path.join(mingw_bin, 'g++.exe')
-        CC = os.path.join(mingw_bin, 'gcc.exe')
-        cmd = ['R.exe', 'CMD', 'INSTALL', str(os.path.curdir)]
-        env = os.environ.copy()
-        env.update({'CC': CC, 'CXX': CXX})
-        subprocess.check_call(cmd, env=env)
-        subprocess.check_call([
-            'R.exe', '-q', '-e',
-            "library(testthat); setwd('tests'); source('testthat.R')"
-        ])
-        subprocess.check_call([
-            'R.exe', '-q', '-e',
-            "demo(runall, package = 'xgboost')"
-        ])
-
-
-def test_with_cmake(args):
-    os.mkdir('build')
-    with DirectoryExcursion('build'):
-        if args.compiler == 'mingw':
+        CXX = os.path.join(mingw_bin, "g++.exe")
+        CC = os.path.join(mingw_bin, "gcc.exe")
+        env.update({"CC": CC, "CXX": CXX})
+
+    status = subprocess.run(
+        [R, "CMD", "check", "--as-cran", "--timings", path], env=env
+    )
+    rcheck_dir = Path("xgboost.Rcheck")
+    with open(rcheck_dir / "00check.log", "r") as fd:
+        check_log = fd.read()
+
+    with open(rcheck_dir / "00install.out", "r") as fd:
+        install_log = fd.read()
+
+    msg = f"""
+----------------------- Install ----------------------
+{install_log}
+
+-----------------------  Check -----------------------
+{check_log}
+
+    """
+
+    if status.returncode != 0:
+        print(msg)
+        raise ValueError("Failed r package check.")
+
+    if check_log.find("WARNING") != -1:
+        print(msg)
+        raise ValueError("Has unresolved warnings.")
+    if check_log.find("Examples with CPU time") != -1:
+        print(msg)
+        raise ValueError("Suspicious NOTE.")
+    if pd is not None:
+        check_example_timing(rcheck_dir, threshold)
+
+
+@cd(R_PACKAGE)
+@record_time
+def check_rmarkdown() -> None:
+    assert system() != "Windows", "Document test doesn't support Windows."
+    env = os.environ.copy()
+    env.update({"MAKEFLAGS": f"-j{os.cpu_count()}"})
+    print("Checking R documentation.")
+    bin_dir = os.path.dirname(R)
+    rscript = os.path.join(bin_dir, "Rscript")
+    subprocess.check_call([rscript, "-e", "roxygen2::roxygenize()"], env=env)
+    output = subprocess.run(["git", "diff", "--name-only"], capture_output=True)
+    if len(output.stdout.decode("utf-8").strip()) != 0:
+        output = subprocess.run(["git", "diff"], capture_output=True)
+        raise ValueError(
+            "Please run `roxygen2::roxygenize()`. Diff:\n",
+            output.stdout.decode("utf-8"),
+        )
+
+
+@cd(R_PACKAGE)
+@record_time
+def test_with_autotools() -> None:
+    """Windows only test. No `--as-cran` check, only unittests. We don't want to manage
+    the dependencies on Windows machine.
+
+    """
+    assert system() == "Windows"
+    mingw_bin = get_mingw_bin()
+    CXX = os.path.join(mingw_bin, "g++.exe")
+    CC = os.path.join(mingw_bin, "gcc.exe")
+    cmd = [R, "CMD", "INSTALL", str(os.path.curdir)]
+    env = os.environ.copy()
+    env.update({"CC": CC, "CXX": CXX, "MAKEFLAGS": f"-j{os.cpu_count()}"})
+    subprocess.check_call(cmd, env=env)
+    subprocess.check_call(
+        ["R.exe", "-q", "-e", "library(testthat); setwd('tests'); source('testthat.R')"]
+    )
+    subprocess.check_call(["R.exe", "-q", "-e", "demo(runall, package = 'xgboost')"])
+
+
+@record_time
+def test_with_cmake(args: argparse.Namespace) -> None:
+    os.mkdir("build")
+    with DirectoryExcursion("build"):
+        if args.compiler == "mingw":
             mingw_bin = get_mingw_bin()
-            CXX = os.path.join(mingw_bin, 'g++.exe')
-            CC = os.path.join(mingw_bin, 'gcc.exe')
+            CXX = os.path.join(mingw_bin, "g++.exe")
+            CC = os.path.join(mingw_bin, "gcc.exe")
             env = os.environ.copy()
-            env.update({'CC': CC, 'CXX': CXX})
-            subprocess.check_call([
-                'cmake', os.path.pardir, '-DUSE_OPENMP=ON', '-DR_LIB=ON',
-                '-DCMAKE_CONFIGURATION_TYPES=Release', '-G', 'Unix Makefiles',
-            ],
-                                  env=env)
-            subprocess.check_call(['make', '-j', 'install'])
-        elif args.compiler == 'msvc':
-            subprocess.check_call([
-                'cmake', os.path.pardir, '-DUSE_OPENMP=ON', '-DR_LIB=ON',
-                '-DCMAKE_CONFIGURATION_TYPES=Release', '-A', 'x64'
-            ])
-            subprocess.check_call([
-                'cmake', '--build', os.path.curdir, '--target', 'install',
-                '--config', 'Release'
-            ])
+            env.update({"CC": CC, "CXX": CXX})
+            subprocess.check_call(
+                [
+                    "cmake",
+                    os.path.pardir,
+                    "-DUSE_OPENMP=ON",
+                    "-DR_LIB=ON",
+                    "-DCMAKE_CONFIGURATION_TYPES=Release",
+                    "-G",
+                    "Unix Makefiles",
+                ],
+                env=env,
+            )
+            subprocess.check_call(["make", "-j", "install"])
+        elif args.compiler == "msvc":
+            subprocess.check_call(
+                [
+                    "cmake",
+                    os.path.pardir,
+                    "-DUSE_OPENMP=ON",
+                    "-DR_LIB=ON",
+                    "-DCMAKE_CONFIGURATION_TYPES=Release",
+                    "-A",
+                    "x64",
+                ]
+            )
+            subprocess.check_call(
+                [
+                    "cmake",
+                    "--build",
+                    os.path.curdir,
+                    "--target",
+                    "install",
+                    "--config",
+                    "Release",
+                ]
+            )
         else:
-            raise ValueError('Wrong compiler')
-    with DirectoryExcursion(r_package):
-        subprocess.check_call([
-            'R.exe', '-q', '-e',
-            "library(testthat); setwd('tests'); source('testthat.R')"
-        ])
-        subprocess.check_call([
-            'R.exe', '-q', '-e',
-            "demo(runall, package = 'xgboost')"
-        ])
-
-
-def main(args):
-    if args.build_tool == 'autotools':
-        test_with_autotools(args)
+            raise ValueError("Wrong compiler")
+    with DirectoryExcursion(R_PACKAGE):
+        subprocess.check_call(
+            [
+                R,
+                "-q",
+                "-e",
+                "library(testthat); setwd('tests'); source('testthat.R')",
+            ]
+        )
+        subprocess.check_call([R, "-q", "-e", "demo(runall, package = 'xgboost')"])
+
+
+@record_time
+def main(args: argparse.Namespace) -> None:
+    if args.task == "pack":
+        pack_rpackage()
+    elif args.task == "build":
+        src_dir = pack_rpackage()
+        build_rpackage(src_dir)
+    elif args.task == "doc":
+        check_rmarkdown()
+    elif args.task == "check":
+        if args.build_tool == "autotools" and system() != "Windows":
+            src_dir = pack_rpackage()
+            tarball = build_rpackage(src_dir)
+            check_rpackage(tarball)
+        elif args.build_tool == "autotools":
+            test_with_autotools()
+        else:
+            test_with_cmake(args)
+    elif args.task == "timings":
+        check_example_timing(Path("xgboost.Rcheck"), 2.5)
     else:
-        test_with_cmake(args)
+        raise ValueError("Unexpected task.")
 
 
-if __name__ == '__main__':
-    parser = argparse.ArgumentParser()
-    parser.add_argument('--compiler',
-                        type=str,
-                        choices=['mingw', 'msvc'],
-                        help='Compiler used for compiling CXX code.')
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=(
+            "Helper script for making R package and running R tests on CI. There are"
+            " also other helper scripts in the R tests directory for installing"
+            " dependencies and running linter."
+        )
+    )
+    parser.add_argument(
+        "--task",
+        type=str,
+        choices=["pack", "build", "check", "doc", "timings"],
+        default="check",
+        required=False,
+    )
     parser.add_argument(
-        '--build-tool',
+        "--compiler",
         type=str,
-        choices=['cmake', 'autotools'],
-        help='Build tool for compiling CXX code and install R package.')
+        choices=["mingw", "msvc"],
+        help="Compiler used for compiling CXX code. Only relevant for windows build",
+        default="mingw",
+        required=False,
+    )
+    parser.add_argument(
+        "--build-tool",
+        type=str,
+        choices=["cmake", "autotools"],
+        help="Build tool for compiling CXX code and install R package.",
+        default="autotools",
+        required=False,
+    )
+    parser.add_argument(
+        "--r",
+        type=str,
+        default="R" if system() != "Windows" else "R.exe",
+        help="Path to the R executable.",
+    )
     args = parser.parse_args()
-    main(args)
+    R = args.r
+
+    try:
+        main(args)
+    finally:
+        print_time()
diff --git a/tests/ci_build/test_utils.py b/tests/ci_build/test_utils.py
index f34f8b66a62a..d4fa02629e34 100644
--- a/tests/ci_build/test_utils.py
+++ b/tests/ci_build/test_utils.py
@@ -1,14 +1,82 @@
+"""Utilities for the CI."""
 import os
-from typing import Union
+from datetime import datetime, timedelta
+from functools import wraps
+from typing import Any, Callable, Dict, TypedDict, TypeVar, Union
 
 
 class DirectoryExcursion:
-    def __init__(self, path: Union[os.PathLike, str]):
+    def __init__(self, path: Union[os.PathLike, str]) -> None:
         self.path = path
         self.curdir = os.path.normpath(os.path.abspath(os.path.curdir))
 
-    def __enter__(self):
+    def __enter__(self) -> None:
         os.chdir(self.path)
 
-    def __exit__(self, *args):
+    def __exit__(self, *args: Any) -> None:
         os.chdir(self.curdir)
+
+
+R = TypeVar("R")
+
+
+def cd(path: Union[os.PathLike, str]) -> Callable:
+    """Decorator for changing directory temporarily."""
+
+    def chdir(func: Callable[..., R]) -> Callable[..., R]:
+        @wraps(func)
+        def inner(*args: Any, **kwargs: Any) -> R:
+            with DirectoryExcursion(path):
+                return func(*args, **kwargs)
+
+        return inner
+
+    return chdir
+
+
+Record = TypedDict("Record", {"count": int, "total": timedelta})
+timer: Dict[str, Record] = {}
+
+
+def record_time(func: Callable[..., R]) -> Callable[..., R]:
+    """Decorator for recording function runtime."""
+    global timer
+
+    @wraps(func)
+    def inner(*args: Any, **kwargs: Any) -> R:
+        if func.__name__ not in timer:
+            timer[func.__name__] = {"count": 0, "total": timedelta(0)}
+        s = datetime.now()
+        try:
+            r = func(*args, **kwargs)
+        finally:
+            e = datetime.now()
+            timer[func.__name__]["count"] += 1
+            timer[func.__name__]["total"] += e - s
+        return r
+
+    return inner
+
+
+def print_time() -> None:
+    """Print all recorded items by :py:func:`record_time`."""
+    global timer
+    for k, v in timer.items():
+        print(
+            "Name:",
+            k,
+            "Called:",
+            v["count"],
+            "Elapsed:",
+            f"{v['total'].seconds} secs",
+        )
+
+
+ROOT = os.path.normpath(
+    os.path.join(
+        os.path.dirname(os.path.abspath(__file__)), os.path.pardir, os.path.pardir
+    )
+)
+R_PACKAGE = os.path.join(ROOT, "R-package")
+JVM_PACKAGES = os.path.join(ROOT, "jvm-packages")
+PY_PACKAGE = os.path.join(ROOT, "python-package")