diff --git a/R/fit.R b/R/fit.R
index 620d99af..51b05b8f 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -1484,10 +1484,10 @@ CmdStanMCMC <- R6::R6Class(
 #'   and the \pkg{loo} package [vignettes](https://mc-stan.org/loo/articles/)
 #'   for details.
 #'
-#' @param variables (character vector) The name(s) of the variable(s) in the
-#'   Stan program containing the pointwise log-likelihood. The default is to
-#'   look for `"log_lik"`. This argument is passed to the
-#'   [`$draws()`][fit-method-draws] method.
+#' @param variables (string) The name of the variable in the Stan program
+#'   containing the pointwise log-likelihood. The default is to look for
+#'   `"log_lik"`. This argument is passed to the [`$draws()`][fit-method-draws]
+#'   method.
 #' @param r_eff (multiple options) How to handle the `r_eff` argument for `loo()`:
 #'   * `TRUE` (the default) will automatically call [loo::relative_eff.array()]
 #'   to compute the `r_eff` argument to pass to [loo::loo.array()].
@@ -1526,6 +1526,9 @@ CmdStanMCMC <- R6::R6Class(
 #'
 loo <- function(variables = "log_lik", r_eff = TRUE, moment_match = FALSE, ...) {
   require_suggested_package("loo")
+  if (length(variables) != 1) {
+    stop("Only a single variable name is allowed for the 'variables' argument.", call. = FALSE)
+  }
   LLarray <- self$draws(variables, format = "draws_array")
   if (is.logical(r_eff)) {
     if (isTRUE(r_eff)) {
@@ -1793,6 +1796,11 @@ CmdStanMCMC$set("public", name = "num_chains", value = num_chains)
 #'
 #' @description A `CmdStanMLE` object is the fitted model object returned by the
 #'   [`$optimize()`][model-method-optimize] method of a [`CmdStanModel`] object.
+#'   This object will either contain a penalized maximum likelihood estimate
+#'   (MLE) or a maximum a posteriori estimate (MAP), depending on the value of
+#'   the `jacobian` argument when the model is fit (and whether the model has
+#'   constrained parameters). See [`$optimize()`][model-method-optimize] and the
+#'   CmdStan User's Guide for more details.
 #'
 #' @section Methods: `CmdStanMLE` objects have the following associated methods,
 #'   all of which have their own (linked) documentation pages.
@@ -1862,17 +1870,22 @@ CmdStanMLE <- R6::R6Class(
   )
 )
 
-#' Extract (penalized) maximum likelihood estimate after optimization
+#' Extract point estimate after optimization
 #'
 #' @name fit-method-mle
 #' @aliases mle
-#' @description The `$mle()` method is only available for [`CmdStanMLE`] objects.
-#' It returns the penalized maximum likelihood estimate (posterior mode) as a
-#' numeric vector with one element per variable. The returned vector does *not*
-#' include `lp__`, the total log probability (`target`) accumulated in the
-#' model block of the Stan program, which is available via the
-#' [`$lp()`][fit-method-lp] method and also included in the
-#' [`$draws()`][fit-method-draws] method.
+#' @description The `$mle()` method is only available for [`CmdStanMLE`]
+#'   objects. It returns the point estimate as a numeric vector with one element
+#'   per variable. The returned vector does *not* include `lp__`, the total log
+#'   probability (`target`) accumulated in the model block of the Stan program,
+#'   which is available via the [`$lp()`][fit-method-lp] method and also
+#'   included in the [`$draws()`][fit-method-draws] method.
+#'
+#'   For models with constrained parameters that are fit with `jacobian=TRUE`,
+#'   the `$mle()` method actually returns the maximum a posteriori (MAP)
+#'   estimate (posterior mode) rather than the MLE. See
+#'   [`$optimize()`][model-method-optimize] and the CmdStan User's Guide for
+#'   more details.
 #'
 #' @param variables (character vector) The variables (parameters, transformed
 #'   parameters, and generated quantities) to include. If NULL (the default)
@@ -1897,6 +1910,7 @@ mle <- function(variables = NULL) {
 }
 CmdStanMLE$set("public", name = "mle", value = mle)
 
+
 # CmdStanLaplace ---------------------------------------------------------------
 #' CmdStanLaplace objects
 #'
diff --git a/R/model.R b/R/model.R
index e63c24e5..c09bc74c 100644
--- a/R/model.R
+++ b/R/model.R
@@ -1353,8 +1353,14 @@ CmdStanModel$set("public", name = "sample_mpi", value = sample_mpi)
 #' @family CmdStanModel methods
 #'
 #' @description The `$optimize()` method of a [`CmdStanModel`] object runs
-#'   Stan's optimizer to obtain a (penalized) maximum likelihood estimate or a
-#'   maximum a posteriori estimate (if `jacobian=TRUE`). See the
+#'   Stan's optimizer to obtain a (penalized) maximum likelihood estimate (MLE)
+#'   or a maximum a posteriori estimate (MAP), depending on the value of the
+#'   `jacobian` argument. For models with constrained parameters, when the
+#'   Jacobian adjustment is *not* applied, the point estimate corresponds to a
+#'   penalized MLE, and when the Jacobian adjustment is applied the point
+#'   estimate corresponds to the MAP (posterior mode) of the model we would fit
+#'   if we were instead doing MCMC sampling. The Jacobian adjustment has no
+#'   affect if the model has only unconstrained parameters. See the
 #'   [CmdStan User's Guide](https://mc-stan.org/docs/cmdstan-guide/index.html)
 #'   for more details.
 #'
@@ -1364,6 +1370,7 @@ CmdStanModel$set("public", name = "sample_mpi", value = sample_mpi)
 #'   default arguments. The default values can also be obtained by checking the
 #'   metadata of an example model, e.g.,
 #'   `cmdstanr_example(method="optimize")$metadata()`.
+#'
 #' @template model-common-args
 #' @param threads (positive integer) If the model was
 #'   [compiled][model-method-compile] with threading support, the number of
diff --git a/cmdstanr.Rproj b/cmdstanr.Rproj
index 270314b8..d33dd97c 100644
--- a/cmdstanr.Rproj
+++ b/cmdstanr.Rproj
@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: b2cb035a-82d0-45fa-b984-ec75db1e1feb
 
 RestoreWorkspace: Default
 SaveWorkspace: Default
diff --git a/man/CmdStanMLE.Rd b/man/CmdStanMLE.Rd
index 82c11910..541f9c15 100644
--- a/man/CmdStanMLE.Rd
+++ b/man/CmdStanMLE.Rd
@@ -6,6 +6,11 @@
 \description{
 A \code{CmdStanMLE} object is the fitted model object returned by the
 \code{\link[=model-method-optimize]{$optimize()}} method of a \code{\link{CmdStanModel}} object.
+This object will either contain a penalized maximum likelihood estimate
+(MLE) or a maximum a posteriori estimate (MAP), depending on the value of
+the \code{jacobian} argument when the model is fit (and whether the model has
+constrained parameters). See \code{\link[=model-method-optimize]{$optimize()}} and the
+CmdStan User's Guide for more details.
 }
 \section{Methods}{
  \code{CmdStanMLE} objects have the following associated methods,
diff --git a/man/fit-method-loo.Rd b/man/fit-method-loo.Rd
index c6d98844..8a63681a 100644
--- a/man/fit-method-loo.Rd
+++ b/man/fit-method-loo.Rd
@@ -8,10 +8,10 @@
 loo(variables = "log_lik", r_eff = TRUE, moment_match = FALSE, ...)
 }
 \arguments{
-\item{variables}{(character vector) The name(s) of the variable(s) in the
-Stan program containing the pointwise log-likelihood. The default is to
-look for \code{"log_lik"}. This argument is passed to the
-\code{\link[=fit-method-draws]{$draws()}} method.}
+\item{variables}{(string) The name of the variable in the Stan program
+containing the pointwise log-likelihood. The default is to look for
+\code{"log_lik"}. This argument is passed to the \code{\link[=fit-method-draws]{$draws()}}
+method.}
 
 \item{r_eff}{(multiple options) How to handle the \code{r_eff} argument for \code{loo()}:
 \itemize{
diff --git a/man/fit-method-mle.Rd b/man/fit-method-mle.Rd
index 1e607563..2ece73f6 100644
--- a/man/fit-method-mle.Rd
+++ b/man/fit-method-mle.Rd
@@ -3,7 +3,7 @@
 \name{fit-method-mle}
 \alias{fit-method-mle}
 \alias{mle}
-\title{Extract (penalized) maximum likelihood estimate after optimization}
+\title{Extract point estimate after optimization}
 \usage{
 mle(variables = NULL)
 }
@@ -16,13 +16,18 @@ then all variables are included.}
 A numeric vector. See \strong{Examples}.
 }
 \description{
-The \verb{$mle()} method is only available for \code{\link{CmdStanMLE}} objects.
-It returns the penalized maximum likelihood estimate (posterior mode) as a
-numeric vector with one element per variable. The returned vector does \emph{not}
-include \code{lp__}, the total log probability (\code{target}) accumulated in the
-model block of the Stan program, which is available via the
-\code{\link[=fit-method-lp]{$lp()}} method and also included in the
-\code{\link[=fit-method-draws]{$draws()}} method.
+The \verb{$mle()} method is only available for \code{\link{CmdStanMLE}}
+objects. It returns the point estimate as a numeric vector with one element
+per variable. The returned vector does \emph{not} include \code{lp__}, the total log
+probability (\code{target}) accumulated in the model block of the Stan program,
+which is available via the \code{\link[=fit-method-lp]{$lp()}} method and also
+included in the \code{\link[=fit-method-draws]{$draws()}} method.
+
+For models with constrained parameters that are fit with \code{jacobian=TRUE},
+the \verb{$mle()} method actually returns the maximum a posteriori (MAP)
+estimate (posterior mode) rather than the MLE. See
+\code{\link[=model-method-optimize]{$optimize()}} and the CmdStan User's Guide for
+more details.
 }
 \examples{
 \dontrun{
diff --git a/man/model-method-optimize.Rd b/man/model-method-optimize.Rd
index 7acae8c7..70f79c15 100644
--- a/man/model-method-optimize.Rd
+++ b/man/model-method-optimize.Rd
@@ -210,8 +210,14 @@ A \code{\link{CmdStanMLE}} object.
 }
 \description{
 The \verb{$optimize()} method of a \code{\link{CmdStanModel}} object runs
-Stan's optimizer to obtain a (penalized) maximum likelihood estimate or a
-maximum a posteriori estimate (if \code{jacobian=TRUE}). See the
+Stan's optimizer to obtain a (penalized) maximum likelihood estimate (MLE)
+or a maximum a posteriori estimate (MAP), depending on the value of the
+\code{jacobian} argument. For models with constrained parameters, when the
+Jacobian adjustment is \emph{not} applied, the point estimate corresponds to a
+penalized MLE, and when the Jacobian adjustment is applied the point
+estimate corresponds to the MAP (posterior mode) of the model we would fit
+if we were instead doing MCMC sampling. The Jacobian adjustment has no
+affect if the model has only unconstrained parameters. See the
 \href{https://mc-stan.org/docs/cmdstan-guide/index.html}{CmdStan User's Guide}
 for more details.
 
diff --git a/tests/testthat/test-fit-mcmc.R b/tests/testthat/test-fit-mcmc.R
index 5e81dd71..2d338833 100644
--- a/tests/testthat/test-fit-mcmc.R
+++ b/tests/testthat/test-fit-mcmc.R
@@ -276,6 +276,11 @@ test_that("loo method works if log_lik is available", {
   fit_bernoulli <- testing_fit("bernoulli_log_lik")
   expect_s3_class(suppressWarnings(fit_bernoulli$loo(cores = 1, save_psis = TRUE)), "loo")
   expect_s3_class(suppressWarnings(fit_bernoulli$loo(r_eff = FALSE)), "loo")
+
+  expect_error(
+    fit_bernoulli$loo(variables = c("log_lik", "beta")),
+    "Only a single variable name is allowed"
+  )
 })
 
 test_that("loo method works with moment-matching", {