diff --git a/R/methods.R b/R/methods.R
index 0d3fcdc43..624e518e8 100644
--- a/R/methods.R
+++ b/R/methods.R
@@ -42,16 +42,16 @@
 #' @param ... Arguments passed to [project()] if `object` is not already an
 #'   object returned by [project()].
 #'
-#' @return Let \eqn{S_{\mbox{prj}}}{S_prj} denote the number of (possibly
+#' @return Let \eqn{S_{\mathrm{prj}}}{S_prj} denote the number of (possibly
 #'   clustered) projected posterior draws (short: the number of projected draws)
 #'   and \eqn{N} the number of observations. Then, if the prediction is done for
 #'   one submodel only (i.e., `length(nterms) == 1 || !is.null(solution_terms)`
 #'   in the call to [project()]):
 #'   * [proj_linpred()] returns a `list` with elements `pred` (predictions) and
-#'   `lpd` (log predictive densities). Both elements are \eqn{S_{\mbox{prj}}
+#'   `lpd` (log predictive densities). Both elements are \eqn{S_{\mathrm{prj}}
 #'   \times N}{S_prj x N} matrices.
-#'   * [proj_predict()] returns an \eqn{S_{\mbox{prj}} \times N}{S_prj x N}
-#'   matrix of predictions where \eqn{S_{\mbox{prj}}}{S_prj} denotes
+#'   * [proj_predict()] returns an \eqn{S_{\mathrm{prj}} \times N}{S_prj x N}
+#'   matrix of predictions where \eqn{S_{\mathrm{prj}}}{S_prj} denotes
 #'   `nresample_clusters` in case of clustered projection.
 #'
 #'   If the prediction is done for more than one submodel, the output from above
@@ -689,18 +689,18 @@ print.vsel <- function(x, ...) {
 #'   the lower or upper bound (depending on argument `type`) of the
 #'   normal-approximation confidence interval (with nominal coverage `1 -
 #'   alpha`; see argument `alpha` of [summary.vsel()]) for \eqn{U_k -
-#'   U_{\mbox{base}}}{U_k - U_base} (with \eqn{U_k} denoting the \eqn{k}-th
-#'   submodel's true utility and \eqn{U_{\mbox{base}}}{U_base} denoting the
+#'   U_{\mathrm{base}}}{U_k - U_base} (with \eqn{U_k} denoting the \eqn{k}-th
+#'   submodel's true utility and \eqn{U_{\mathrm{base}}}{U_base} denoting the
 #'   baseline model's true utility) falls above (or is equal to)
-#'   \deqn{\texttt{pct} \cdot (u_0 - u_{\mbox{base}})}{pct * (u_0 - u_base)}
+#'   \deqn{\texttt{pct} \cdot (u_0 - u_{\mathrm{base}})}{pct * (u_0 - u_base)}
 #'   where \eqn{u_0} denotes the null model's estimated utility and
-#'   \eqn{u_{\mbox{base}}}{u_base} the baseline model's estimated utility. The
+#'   \eqn{u_{\mathrm{base}}}{u_base} the baseline model's estimated utility. The
 #'   baseline is either the reference model or the best submodel found (see
 #'   argument `baseline` of [summary.vsel()]).
 #'
 #'   For example, `alpha = 0.32`, `pct = 0`, and `type = "upper"` means that we
 #'   select the smallest model size for which the upper bound of the 68%
-#'   confidence interval for \eqn{U_k - U_{\mbox{base}}}{U_k - U_base} exceeds
+#'   confidence interval for \eqn{U_k - U_{\mathrm{base}}}{U_k - U_base} exceeds
 #'   (or is equal to) zero, that is, for which the submodel's utility estimate
 #'   is at most one standard error smaller than the baseline model's utility
 #'   estimate.
@@ -1056,8 +1056,8 @@ get_subparams.gamm4 <- function(x, ...) {
 #'   uses `"rstanarm"` if the reference model fit is of an unknown class).
 #' @param ... Currently ignored.
 #'
-#' @return An \eqn{S_{\mbox{prj}} \times Q}{S_prj x Q} matrix of projected
-#'   draws, with \eqn{S_{\mbox{prj}}}{S_prj} denoting the number of projected
+#' @return An \eqn{S_{\mathrm{prj}} \times Q}{S_prj x Q} matrix of projected
+#'   draws, with \eqn{S_{\mathrm{prj}}}{S_prj} denoting the number of projected
 #'   draws and \eqn{Q} the number of parameters.
 #'
 #' @examples
diff --git a/R/misc.R b/R/misc.R
index 8d528ed25..57f00cefc 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -181,24 +181,24 @@ bootstrap <- function(x, fun = mean, B = 2000,
 #   subsampled (without replacement).
 #
 # @return Let \eqn{y} denote the response (vector), \eqn{N} the number of
-#   observations, and \eqn{S_{\mbox{prj}}}{S_prj} the number of projected draws
-#   (= either `nclusters` or `ndraws`, depending on which one is used). Then the
-#   return value is a list with elements:
+#   observations, and \eqn{S_{\mathrm{prj}}}{S_prj} the number of projected
+#   draws (= either `nclusters` or `ndraws`, depending on which one is used).
+#   Then the return value is a list with elements:
 #
-#   * `mu`: An \eqn{N \times S_{\mbox{prj}}}{N x S_prj} matrix of expected
+#   * `mu`: An \eqn{N \times S_{\mathrm{prj}}}{N x S_prj} matrix of expected
 #   values for \eqn{y} for each draw/cluster.
-#   * `var`: An \eqn{N \times S_{\mbox{prj}}}{N x S_prj} matrix of predictive
+#   * `var`: An \eqn{N \times S_{\mathrm{prj}}}{N x S_prj} matrix of predictive
 #   variances for \eqn{y} for each draw/cluster which are needed for projecting
 #   the dispersion parameter (the predictive variances are NA for those families
 #   that do not have a dispersion parameter).
-#   * `dis`: A vector of length \eqn{S_{\mbox{prj}}}{S_prj} containing the
+#   * `dis`: A vector of length \eqn{S_{\mathrm{prj}}}{S_prj} containing the
 #   reference model's dispersion parameter value for each draw/cluster (NA for
 #   those families that do not have a dispersion parameter).
-#   * `weights`: A vector of length \eqn{S_{\mbox{prj}}}{S_prj} containing the
+#   * `weights`: A vector of length \eqn{S_{\mathrm{prj}}}{S_prj} containing the
 #   weights for the draws/clusters.
 #   * `cl`: Cluster assignment for each posterior draw, that is, a vector that
 #   has length equal to the number of posterior draws and each value is an
-#   integer between 1 and \eqn{S_{\mbox{prj}}}{S_prj}.
+#   integer between 1 and \eqn{S_{\mathrm{prj}}}{S_prj}.
 .get_refdist <- function(refmodel, ndraws = NULL, nclusters = NULL,
                          thinning = TRUE) {
   # Number of draws in the reference model:
diff --git a/R/refmodel.R b/R/refmodel.R
index f7a85978d..66f228106 100644
--- a/R/refmodel.R
+++ b/R/refmodel.R
@@ -101,8 +101,8 @@
 #' Arguments `ref_predfun`, `proj_predfun`, and `div_minimizer` may be `NULL`
 #' for using an internal default. Otherwise, let \eqn{N} denote the number of
 #' observations (in case of CV, these may be reduced to each fold),
-#' \eqn{S_{\mbox{ref}}}{S_ref} the number of posterior draws for the reference
-#' model's parameters, and \eqn{S_{\mbox{prj}}}{S_prj} the number of (possibly
+#' \eqn{S_{\mathrm{ref}}}{S_ref} the number of posterior draws for the reference
+#' model's parameters, and \eqn{S_{\mathrm{prj}}}{S_prj} the number of (possibly
 #' clustered) parameter draws for projection (short: the number of projected
 #' draws). Then the functions supplied to these arguments need to have the
 #' following prototypes:
@@ -114,7 +114,7 @@
 #'     typically stored in `fit`) or data for new observations (at least in the
 #'     form of a `data.frame`).
 #' * `proj_predfun`: `proj_predfun(fits, newdata)` where:
-#'     + `fits` accepts a `list` of length \eqn{S_{\mbox{prj}}}{S_prj}
+#'     + `fits` accepts a `list` of length \eqn{S_{\mathrm{prj}}}{S_prj}
 #'     containing this number of submodel fits. This `list` is the same as that
 #'     returned by [project()] in its output element `submodl` (which in turn is
 #'     the same as the return value of `div_minimizer`, except if [project()]
@@ -125,15 +125,15 @@
 #' * `div_minimizer` does not need to have a specific prototype, but it needs to
 #' be able to be called with the following arguments:
 #'     + `formula` accepts either a standard [`formula`] with a single response
-#'     (if \eqn{S_{\mbox{prj}} = 1}{S_prj = 1}) or a [`formula`] with
-#'     \eqn{S_{\mbox{prj}} > 1}{S_prj > 1} response variables [cbind()]-ed on
+#'     (if \eqn{S_{\mathrm{prj}} = 1}{S_prj = 1}) or a [`formula`] with
+#'     \eqn{S_{\mathrm{prj}} > 1}{S_prj > 1} response variables [cbind()]-ed on
 #'     the left-hand side in which case the projection has to be performed for
 #'     each of the response variables separately.
 #'     + `data` accepts a `data.frame` to be used for the projection.
 #'     + `family` accepts a [`family`] object.
 #'     + `weights` accepts either observation weights (at least in the form of a
 #'     numeric vector) or `NULL` (for using a vector of ones as weights).
-#'     + `projpred_var` accepts an \eqn{N \times S_{\mbox{prj}}}{N x S_prj}
+#'     + `projpred_var` accepts an \eqn{N \times S_{\mathrm{prj}}}{N x S_prj}
 #'     matrix of predictive variances (necessary for \pkg{projpred}'s internal
 #'     GLM fitter).
 #'     + `projpred_regul` accepts a single numeric value as supplied to argument
@@ -141,10 +141,10 @@
 #'     + `...` accepts further arguments specified by the user.
 #'
 #' The return value of these functions needs to be:
-#' * `ref_predfun`: an \eqn{N \times S_{\mbox{ref}}}{N x S_ref} matrix.
-#' * `proj_predfun`: an \eqn{N \times S_{\mbox{prj}}}{N x S_prj} matrix.
-#' * `div_minimizer`: a `list` of length \eqn{S_{\mbox{prj}}}{S_prj} containing
-#' this number of submodel fits.
+#' * `ref_predfun`: an \eqn{N \times S_{\mathrm{ref}}}{N x S_ref} matrix.
+#' * `proj_predfun`: an \eqn{N \times S_{\mathrm{prj}}}{N x S_prj} matrix.
+#' * `div_minimizer`: a `list` of length \eqn{S_{\mathrm{prj}}}{S_prj}
+#' containing this number of submodel fits.
 #'
 #' # Argument `extract_model_data`
 #'
diff --git a/man/as.matrix.projection.Rd b/man/as.matrix.projection.Rd
index d835af0e1..567d4e63c 100644
--- a/man/as.matrix.projection.Rd
+++ b/man/as.matrix.projection.Rd
@@ -18,8 +18,8 @@ uses \code{"rstanarm"} if the reference model fit is of an unknown class).}
 \item{...}{Currently ignored.}
 }
 \value{
-An \eqn{S_{\mbox{prj}} \times Q}{S_prj x Q} matrix of projected
-draws, with \eqn{S_{\mbox{prj}}}{S_prj} denoting the number of projected
+An \eqn{S_{\mathrm{prj}} \times Q}{S_prj x Q} matrix of projected
+draws, with \eqn{S_{\mathrm{prj}}}{S_prj} denoting the number of projected
 draws and \eqn{Q} the number of parameters.
 }
 \description{
diff --git a/man/pred-projection.Rd b/man/pred-projection.Rd
index f55151ffd..9f91f4c8d 100644
--- a/man/pred-projection.Rd
+++ b/man/pred-projection.Rd
@@ -82,17 +82,17 @@ the set of clustered posterior draws after projection (with this set being
 determined by argument \code{nclusters} of \code{\link[=project]{project()}}).}
 }
 \value{
-Let \eqn{S_{\mbox{prj}}}{S_prj} denote the number of (possibly
+Let \eqn{S_{\mathrm{prj}}}{S_prj} denote the number of (possibly
 clustered) projected posterior draws (short: the number of projected draws)
 and \eqn{N} the number of observations. Then, if the prediction is done for
 one submodel only (i.e., \code{length(nterms) == 1 || !is.null(solution_terms)}
 in the call to \code{\link[=project]{project()}}):
 \itemize{
 \item \code{\link[=proj_linpred]{proj_linpred()}} returns a \code{list} with elements \code{pred} (predictions) and
-\code{lpd} (log predictive densities). Both elements are \eqn{S_{\mbox{prj}}
+\code{lpd} (log predictive densities). Both elements are \eqn{S_{\mathrm{prj}}
   \times N}{S_prj x N} matrices.
-\item \code{\link[=proj_predict]{proj_predict()}} returns an \eqn{S_{\mbox{prj}} \times N}{S_prj x N}
-matrix of predictions where \eqn{S_{\mbox{prj}}}{S_prj} denotes
+\item \code{\link[=proj_predict]{proj_predict()}} returns an \eqn{S_{\mathrm{prj}} \times N}{S_prj x N}
+matrix of predictions where \eqn{S_{\mathrm{prj}}}{S_prj} denotes
 \code{nresample_clusters} in case of clustered projection.
 }
 
diff --git a/man/refmodel-init-get.Rd b/man/refmodel-init-get.Rd
index 223fc1b8a..8bcafa28b 100644
--- a/man/refmodel-init-get.Rd
+++ b/man/refmodel-init-get.Rd
@@ -155,8 +155,8 @@ analogously for higher-order joint effects, e.g., of three predictors).
 Arguments \code{ref_predfun}, \code{proj_predfun}, and \code{div_minimizer} may be \code{NULL}
 for using an internal default. Otherwise, let \eqn{N} denote the number of
 observations (in case of CV, these may be reduced to each fold),
-\eqn{S_{\mbox{ref}}}{S_ref} the number of posterior draws for the reference
-model's parameters, and \eqn{S_{\mbox{prj}}}{S_prj} the number of (possibly
+\eqn{S_{\mathrm{ref}}}{S_ref} the number of posterior draws for the reference
+model's parameters, and \eqn{S_{\mathrm{prj}}}{S_prj} the number of (possibly
 clustered) parameter draws for projection (short: the number of projected
 draws). Then the functions supplied to these arguments need to have the
 following prototypes:
@@ -172,7 +172,7 @@ form of a \code{data.frame}).
 }
 \item \code{proj_predfun}: \code{proj_predfun(fits, newdata)} where:
 \itemize{
-\item \code{fits} accepts a \code{list} of length \eqn{S_{\mbox{prj}}}{S_prj}
+\item \code{fits} accepts a \code{list} of length \eqn{S_{\mathrm{prj}}}{S_prj}
 containing this number of submodel fits. This \code{list} is the same as that
 returned by \code{\link[=project]{project()}} in its output element \code{submodl} (which in turn is
 the same as the return value of \code{div_minimizer}, except if \code{\link[=project]{project()}}
@@ -185,15 +185,15 @@ as with \code{refit_prj = FALSE}).
 be able to be called with the following arguments:
 \itemize{
 \item \code{formula} accepts either a standard \code{\link{formula}} with a single response
-(if \eqn{S_{\mbox{prj}} = 1}{S_prj = 1}) or a \code{\link{formula}} with
-\eqn{S_{\mbox{prj}} > 1}{S_prj > 1} response variables \code{\link[=cbind]{cbind()}}-ed on
+(if \eqn{S_{\mathrm{prj}} = 1}{S_prj = 1}) or a \code{\link{formula}} with
+\eqn{S_{\mathrm{prj}} > 1}{S_prj > 1} response variables \code{\link[=cbind]{cbind()}}-ed on
 the left-hand side in which case the projection has to be performed for
 each of the response variables separately.
 \item \code{data} accepts a \code{data.frame} to be used for the projection.
 \item \code{family} accepts a \code{\link{family}} object.
 \item \code{weights} accepts either observation weights (at least in the form of a
 numeric vector) or \code{NULL} (for using a vector of ones as weights).
-\item \code{projpred_var} accepts an \eqn{N \times S_{\mbox{prj}}}{N x S_prj}
+\item \code{projpred_var} accepts an \eqn{N \times S_{\mathrm{prj}}}{N x S_prj}
 matrix of predictive variances (necessary for \pkg{projpred}'s internal
 GLM fitter).
 \item \code{projpred_regul} accepts a single numeric value as supplied to argument
@@ -204,10 +204,10 @@ GLM fitter).
 
 The return value of these functions needs to be:
 \itemize{
-\item \code{ref_predfun}: an \eqn{N \times S_{\mbox{ref}}}{N x S_ref} matrix.
-\item \code{proj_predfun}: an \eqn{N \times S_{\mbox{prj}}}{N x S_prj} matrix.
-\item \code{div_minimizer}: a \code{list} of length \eqn{S_{\mbox{prj}}}{S_prj} containing
-this number of submodel fits.
+\item \code{ref_predfun}: an \eqn{N \times S_{\mathrm{ref}}}{N x S_ref} matrix.
+\item \code{proj_predfun}: an \eqn{N \times S_{\mathrm{prj}}}{N x S_prj} matrix.
+\item \code{div_minimizer}: a \code{list} of length \eqn{S_{\mathrm{prj}}}{S_prj}
+containing this number of submodel fits.
 }
 }
 
diff --git a/man/suggest_size.Rd b/man/suggest_size.Rd
index 0aa1a18a0..3549a58dc 100644
--- a/man/suggest_size.Rd
+++ b/man/suggest_size.Rd
@@ -52,18 +52,18 @@ The suggested model size is the smallest model size \eqn{k \in \{0,
   1, ..., \texttt{nterms\_max\}}}{k = 0, 1, ..., nterms_max} for which either
 the lower or upper bound (depending on argument \code{type}) of the
 normal-approximation confidence interval (with nominal coverage \code{1 - alpha}; see argument \code{alpha} of \code{\link[=summary.vsel]{summary.vsel()}}) for \eqn{U_k -
-  U_{\mbox{base}}}{U_k - U_base} (with \eqn{U_k} denoting the \eqn{k}-th
-submodel's true utility and \eqn{U_{\mbox{base}}}{U_base} denoting the
+  U_{\mathrm{base}}}{U_k - U_base} (with \eqn{U_k} denoting the \eqn{k}-th
+submodel's true utility and \eqn{U_{\mathrm{base}}}{U_base} denoting the
 baseline model's true utility) falls above (or is equal to)
-\deqn{\texttt{pct} \cdot (u_0 - u_{\mbox{base}})}{pct * (u_0 - u_base)}
+\deqn{\texttt{pct} \cdot (u_0 - u_{\mathrm{base}})}{pct * (u_0 - u_base)}
 where \eqn{u_0} denotes the null model's estimated utility and
-\eqn{u_{\mbox{base}}}{u_base} the baseline model's estimated utility. The
+\eqn{u_{\mathrm{base}}}{u_base} the baseline model's estimated utility. The
 baseline is either the reference model or the best submodel found (see
 argument \code{baseline} of \code{\link[=summary.vsel]{summary.vsel()}}).
 
 For example, \code{alpha = 0.32}, \code{pct = 0}, and \code{type = "upper"} means that we
 select the smallest model size for which the upper bound of the 68\%
-confidence interval for \eqn{U_k - U_{\mbox{base}}}{U_k - U_base} exceeds
+confidence interval for \eqn{U_k - U_{\mathrm{base}}}{U_k - U_base} exceeds
 (or is equal to) zero, that is, for which the submodel's utility estimate
 is at most one standard error smaller than the baseline model's utility
 estimate.