diff --git a/R-package/R/xgb.Booster.R b/R-package/R/xgb.Booster.R
index 77d75fa9c2a5..77b33f16db44 100644
--- a/R-package/R/xgb.Booster.R
+++ b/R-package/R/xgb.Booster.R
@@ -249,7 +249,7 @@ xgb.get.handle <- function(object) {
 #' summary(rowSums(pred_contr) - qlogis(pred))
 #' # for the 1st record, let's inspect its features that had non-zero contribution to prediction:
 #' contr1 <- pred_contr[1,]
-#' contr1 <- contr1[-length(contr1)]    # drop BIAS
+#' contr1 <- contr1[-length(contr1)]    # drop intercept
 #' contr1 <- contr1[contr1 != 0]        # drop non-contributing features
 #' contr1 <- contr1[order(abs(contr1))] # order by contribution magnitude
 #' old_mar <- par("mar")
@@ -473,7 +473,7 @@ predict.xgb.Booster <- function(object, newdata, missing = NA, outputmargin = FA
 
   .Call(XGSetArrayDimInplace_R, arr, rev(shape))
 
-  cnames <- if (!is.null(colnames(newdata))) c(colnames(newdata), "BIAS") else NULL
+  cnames <- if (!is.null(colnames(newdata))) c(colnames(newdata), "(Intercept)") else NULL
   n_groups <- shape[2]
 
   ## Needed regardless of whether strict shape is being used.
diff --git a/R-package/man/predict.xgb.Booster.Rd b/R-package/man/predict.xgb.Booster.Rd
index 88a2f203efcd..9c2e434d0625 100644
--- a/R-package/man/predict.xgb.Booster.Rd
+++ b/R-package/man/predict.xgb.Booster.Rd
@@ -211,7 +211,7 @@ str(pred_contr)
 summary(rowSums(pred_contr) - qlogis(pred))
 # for the 1st record, let's inspect its features that had non-zero contribution to prediction:
 contr1 <- pred_contr[1,]
-contr1 <- contr1[-length(contr1)]    # drop BIAS
+contr1 <- contr1[-length(contr1)]    # drop intercept
 contr1 <- contr1[contr1 != 0]        # drop non-contributing features
 contr1 <- contr1[order(abs(contr1))] # order by contribution magnitude
 old_mar <- par("mar")
diff --git a/R-package/tests/testthat/test_helpers.R b/R-package/tests/testthat/test_helpers.R
index 38b5ca0667bf..c619bc50b5da 100644
--- a/R-package/tests/testthat/test_helpers.R
+++ b/R-package/tests/testthat/test_helpers.R
@@ -101,7 +101,7 @@ test_that("predict feature contributions works", {
   # gbtree binary classifier
   expect_error(pred_contr <- predict(bst.Tree, sparse_matrix, predcontrib = TRUE), regexp = NA)
   expect_equal(dim(pred_contr), c(nrow(sparse_matrix), ncol(sparse_matrix) + 1))
-  expect_equal(colnames(pred_contr), c(colnames(sparse_matrix), "BIAS"))
+  expect_equal(colnames(pred_contr), c(colnames(sparse_matrix), "(Intercept)"))
   pred <- predict(bst.Tree, sparse_matrix, outputmargin = TRUE)
   expect_lt(max(abs(rowSums(pred_contr) - pred)), 1e-5)
   # must work with data that has no column names
@@ -114,14 +114,14 @@ test_that("predict feature contributions works", {
   # gbtree binary classifier (approximate method)
   expect_error(pred_contr <- predict(bst.Tree, sparse_matrix, predcontrib = TRUE, approxcontrib = TRUE), regexp = NA)
   expect_equal(dim(pred_contr), c(nrow(sparse_matrix), ncol(sparse_matrix) + 1))
-  expect_equal(colnames(pred_contr), c(colnames(sparse_matrix), "BIAS"))
+  expect_equal(colnames(pred_contr), c(colnames(sparse_matrix), "(Intercept)"))
   pred <- predict(bst.Tree, sparse_matrix, outputmargin = TRUE)
   expect_lt(max(abs(rowSums(pred_contr) - pred)), 1e-5)
 
   # gblinear binary classifier
   expect_error(pred_contr <- predict(bst.GLM, sparse_matrix, predcontrib = TRUE), regexp = NA)
   expect_equal(dim(pred_contr), c(nrow(sparse_matrix), ncol(sparse_matrix) + 1))
-  expect_equal(colnames(pred_contr), c(colnames(sparse_matrix), "BIAS"))
+  expect_equal(colnames(pred_contr), c(colnames(sparse_matrix), "(Intercept)"))
   pred <- predict(bst.GLM, sparse_matrix, outputmargin = TRUE)
   expect_lt(max(abs(rowSums(pred_contr) - pred)), 1e-5)
   # manual calculation of linear terms
@@ -137,7 +137,7 @@ test_that("predict feature contributions works", {
   expect_is(pred_contr, "list")
   expect_length(pred_contr, 3)
   for (g in seq_along(pred_contr)) {
-    expect_equal(colnames(pred_contr[[g]]), c(colnames(iris[, -5]), "BIAS"))
+    expect_equal(colnames(pred_contr[[g]]), c(colnames(iris[, -5]), "(Intercept)"))
     expect_lt(max(abs(rowSums(pred_contr[[g]]) - pred[, g])), 1e-5)
   }
 
@@ -151,7 +151,7 @@ test_that("predict feature contributions works", {
     byrow = TRUE
   )
   for (g in seq_along(pred_contr)) {
-    expect_equal(colnames(pred_contr[[g]]), c(colnames(iris[, -5]), "BIAS"))
+    expect_equal(colnames(pred_contr[[g]]), c(colnames(iris[, -5]), "(Intercept)"))
     expect_lt(max(abs(rowSums(pred_contr[[g]]) - pred[, g])), float_tolerance)
     # manual calculation of linear terms
     coefs <- c(coefs_all[-1, g], coefs_all[1, g]) # intercept needs to be the last
diff --git a/R-package/tests/testthat/test_interactions.R b/R-package/tests/testthat/test_interactions.R
index 645efc12a14c..60cf9d80039a 100644
--- a/R-package/tests/testthat/test_interactions.R
+++ b/R-package/tests/testthat/test_interactions.R
@@ -48,7 +48,7 @@ test_that("predict feature interactions works", {
   intr <- predict(b, dm, predinteraction = TRUE)
   expect_equal(dim(intr), c(N, P + 1, P + 1))
   # check assigned colnames
-  cn <- c(letters[1:P], "BIAS")
+  cn <- c(letters[1:P], "(Intercept)")
   expect_equal(dimnames(intr), list(NULL, cn, cn))
 
   # check the symmetry
@@ -60,7 +60,7 @@ test_that("predict feature interactions works", {
   # diagonal terms for features 3,4,5 must be close to zero
   expect_lt(Reduce(max, sapply(3:P, function(i) max(abs(intr[, i, i])))), 0.05)
 
-  # BIAS must have no interactions
+  # Intercept must have no interactions
   expect_lt(max(abs(intr[, 1:P, P + 1])), 0.00001)
 
   # interactions other than 2 x 3 must be close to zero