diff --git a/stan/math/mix/prob.hpp b/stan/math/mix/prob.hpp
index 99babe99e3f..a2fa67c4079 100644
--- a/stan/math/mix/prob.hpp
+++ b/stan/math/mix/prob.hpp
@@ -3,13 +3,11 @@
 
 #include <stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp>
 #include <stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp>
-#include <stan/math/mix/prob/laplace_latent_poisson_log_2_rng.hpp>
 #include <stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp>
 #include <stan/math/mix/prob/laplace_latent_rng.hpp>
 #include <stan/math/mix/prob/laplace_marginal.hpp>
 #include <stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp>
 #include <stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp>
-#include <stan/math/mix/prob/laplace_marginal_poisson_log_2_lpmf.hpp>
 #include <stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp>
 
 #endif
diff --git a/stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp b/stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp
index c30b47e8c96..8d3cde2a36d 100644
--- a/stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp
+++ b/stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp
@@ -19,19 +19,22 @@ namespace math {
  * where the likelihood is a Bernoulli with logit link.
  * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
  * with dynamic sized rows and 1 column.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @tparam RNG A valid boost rng type
  * @param[in] y Vector Vector of total number of trials with a positive outcome.
  * @param[in] n_samples Vector of number of trials.
+ * @param[in] mean the mean of the latent normal variable.
  * \laplace_common_args
  * \laplace_options
  * \rng_arg
  * \msg_arg
  */
-template <typename ThetaVec, typename CovarFun, typename CovarArgs,
-          typename RNG, require_eigen_t<ThetaVec>* = nullptr>
+template <typename ThetaVec, typename Mean, typename CovarFun,
+          typename CovarArgs, typename RNG,
+          require_eigen_vector_t<ThetaVec>* = nullptr>
 inline Eigen::VectorXd laplace_latent_tol_bernoulli_logit_rng(
-    const std::vector<int>& y, const std::vector<int>& n_samples,
+    const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args, ThetaVec&& theta_0,
     const double tolerance, const int max_num_steps,
     const int hessian_block_size, const int solver,
@@ -39,10 +42,11 @@ inline Eigen::VectorXd laplace_latent_tol_bernoulli_logit_rng(
   laplace_options_user_supplied ops{hessian_block_size,    solver,
                                     max_steps_line_search, tolerance,
                                     max_num_steps,         value_of(theta_0)};
-  return laplace_base_rng(bernoulli_logit_likelihood{},
-                          std::forward_as_tuple(to_vector(y), n_samples),
-                          std::forward<CovarFun>(covariance_function),
-                          std::forward<CovarArgs>(covar_args), ops, rng, msgs);
+  return laplace_base_rng(
+      bernoulli_logit_likelihood{},
+      std::forward_as_tuple(to_vector(y), n_samples, std::forward<Mean>(mean)),
+      std::forward<CovarFun>(covariance_function),
+      std::forward<CovarArgs>(covar_args), ops, rng, msgs);
 }
 
 /**
@@ -54,24 +58,27 @@ inline Eigen::VectorXd laplace_latent_tol_bernoulli_logit_rng(
  * return a multivariate normal random variate sampled
  * from the gaussian approximation of p(theta | y, phi),
  * where the likelihood is a Bernoulli with logit link.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @tparam RNG A valid boost rng type
  * @param[in] y Vector Vector of total number of trials with a positive outcome.
  * @param[in] n_samples Vector of number of trials.
+ * @param[in] mean the mean of the latent normal variable.
  * \laplace_common_args
  * \rng_arg
  * \msg_arg
  */
-template <typename CovarFun, typename CovarArgs, typename RNG>
+template <typename Mean, typename CovarFun, typename CovarArgs, typename RNG>
 inline Eigen::VectorXd laplace_latent_bernoulli_logit_rng(
-    const std::vector<int>& y, const std::vector<int>& n_samples,
+    const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args, RNG& rng,
     std::ostream* msgs) {
-  return laplace_base_rng(bernoulli_logit_likelihood{},
-                          std::forward_as_tuple(to_vector(y), n_samples),
-                          std::forward<CovarFun>(covariance_function),
-                          std::forward<CovarArgs>(covar_args),
-                          laplace_options_default{}, rng, msgs);
+  return laplace_base_rng(
+      bernoulli_logit_likelihood{},
+      std::forward_as_tuple(to_vector(y), n_samples, std::forward<Mean>(mean)),
+      std::forward<CovarFun>(covariance_function),
+      std::forward<CovarArgs>(covar_args), laplace_options_default{}, rng,
+      msgs);
 }
 
 }  // namespace math
diff --git a/stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp b/stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp
index 9aae00210f4..ccdf4675c28 100644
--- a/stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp
+++ b/stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp
@@ -23,23 +23,25 @@ namespace math {
  * @tparam Eta A type for the overdispersion parameter.
  * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
  * with dynamic sized rows and 1 column.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @tparam RNG A valid boost rng type
  * @param[in] y Observed counts.
  * @param[in] y_index Index indicating which group each observation belongs to.
  * @param[in] eta Overdisperison parameter.
+ * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
  * \laplace_options
  * \rng_arg
  * \msg_arg
  */
-template <typename Eta, typename ThetaVec, typename CovarFun,
+template <typename Eta, typename ThetaVec, typename Mean, typename CovarFun,
           typename CovarArgs, typename RNG,
-          require_eigen_t<ThetaVec>* = nullptr>
+          require_eigen_vector_t<ThetaVec>* = nullptr>
 inline Eigen::VectorXd laplace_latent_tol_neg_binomial_2_log_rng(
     const std::vector<int>& y, const std::vector<int>& y_index, Eta&& eta,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, ThetaVec&& theta_0,
-    const double tolerance, const int max_num_steps,
+    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
+    ThetaVec&& theta_0, const double tolerance, const int max_num_steps,
     const int hessian_block_size, const int solver,
     const int max_steps_line_search, RNG& rng, std::ostream* msgs) {
   laplace_options_user_supplied ops{hessian_block_size,    solver,
@@ -47,7 +49,8 @@ inline Eigen::VectorXd laplace_latent_tol_neg_binomial_2_log_rng(
                                     max_num_steps,         value_of(theta_0)};
   return laplace_base_rng(
       neg_binomial_2_log_likelihood{},
-      std::forward_as_tuple(std::forward<Eta>(eta), y, y_index),
+      std::forward_as_tuple(std::forward<Eta>(eta), y, y_index,
+                            std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), ops, rng, msgs);
 }
@@ -65,23 +68,27 @@ inline Eigen::VectorXd laplace_latent_tol_neg_binomial_2_log_rng(
  * parameterization of the Negative Binomial.
  *
  * @tparam Eta A type for the overdispersion parameter.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @tparam RNG A valid boost rng type
  * @param[in] y Observed counts.
  * @param[in] y_index Index indicating which group each observation belongs to.
  * @param[in] eta Overdisperison parameter.
+ * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
  * \rng_arg
  * \msg_arg
  */
-template <typename Eta, typename CovarFun, typename CovarArgs, typename RNG>
+template <typename Eta, typename Mean, typename CovarFun, typename CovarArgs,
+          typename RNG>
 inline Eigen::VectorXd laplace_latent_neg_binomial_2_log_rng(
     const std::vector<int>& y, const std::vector<int>& y_index, Eta&& eta,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, RNG& rng,
-    std::ostream* msgs) {
+    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
+    RNG& rng, std::ostream* msgs) {
   return laplace_base_rng(
       neg_binomial_2_log_likelihood{},
-      std::forward_as_tuple(std::forward<Eta>(eta), y, y_index),
+      std::forward_as_tuple(std::forward<Eta>(eta), y, y_index,
+                            std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), laplace_options_default{}, rng,
       msgs);
diff --git a/stan/math/mix/prob/laplace_latent_poisson_log_2_rng.hpp b/stan/math/mix/prob/laplace_latent_poisson_log_2_rng.hpp
deleted file mode 100644
index bcdc0581fc7..00000000000
--- a/stan/math/mix/prob/laplace_latent_poisson_log_2_rng.hpp
+++ /dev/null
@@ -1,85 +0,0 @@
-#ifndef STAN_MATH_MIX_PROB_LAPLACE_LATENT_POISSON_LOG_2_RNG_HPP
-#define STAN_MATH_MIX_PROB_LAPLACE_LATENT_POISSON_LOG_2_RNG_HPP
-
-#include <stan/math/mix/functor/laplace_base_rng.hpp>
-#include <stan/math/mix/functor/laplace_likelihood.hpp>
-#include <stan/math/mix/prob/laplace_marginal_poisson_log_2_lpmf.hpp>
-
-namespace stan {
-namespace math {
-
-/**
- * Wrapper function around the laplace_marginal function for
- * a log poisson likelihood with exposure. Returns the marginal density
- * p(y | phi) by marginalizing out the latent gaussian variable,
- * with a Laplace approximation. See the laplace_marginal function
- * for more details.
- *
- * @tparam YeVec A type inheriting from `Eigen::EigenBase` with dynamic
- *  sized rows and 1 column.
- * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
- * with dynamic sized rows and 1 column.
- * \laplace_common_template_args
- * @tparam RNG A valid boost rng type
- * @param[in] y total counts per group. Second sufficient statistics.
- * @param[in] y_index group to which each observation belongs.
- * @param[in] ye the exposure for each group.
- * \laplace_common_args
- * \laplace_options
- * \rng_arg
- * \msg_arg
- */
-template <typename YeVec, typename ThetaVec, typename CovarFun,
-          typename CovarArgs, typename RNG,
-          require_eigen_t<ThetaVec>* = nullptr>
-inline auto laplace_latent_tol_poisson_2_log_rng(
-    const std::vector<int>& y, const std::vector<int>& y_index, const YeVec& ye,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, ThetaVec&& theta_0,
-    const double tolerance, const int max_num_steps,
-    const int hessian_block_size, const int solver,
-    const int max_steps_line_search, RNG& rng, std::ostream* msgs) {
-  laplace_options_user_supplied ops{hessian_block_size,    solver,
-                                    max_steps_line_search, tolerance,
-                                    max_num_steps,         value_of(theta_0)};
-  return laplace_base_rng(poisson_log_2_likelihood{},
-                          std::forward_as_tuple(y, y_index, ye),
-                          std::forward<CovarFun>(covariance_function),
-                          std::forward<CovarArgs>(covar_args), ops, rng, msgs);
-}
-
-/**
- * Wrapper function around the laplace_marginal function for
- * a log poisson likelihood with exposure. Returns the marginal density
- * p(y | phi) by marginalizing out the latent gaussian variable,
- * with a Laplace approximation. See the laplace_marginal function
- * for more details.
- *
- * @tparam YeVec A type inheriting from `Eigen::EigenBase` with dynamic
- *  sized rows and 1 column.
- * \laplace_common_template_args
- * @tparam RNG A valid boost rng type
- * @param[in] y total counts per group. Second sufficient statistics.
- * @param[in] y_index group to which each observation belongs.
- * @param[in] ye the exposure for each group.
- * \laplace_common_args
- * \rng_arg
- * \msg_arg
- */
-template <typename YeVec, typename CovarFun, typename CovarArgs, typename RNG>
-inline auto laplace_latent_poisson_2_log_rng(const std::vector<int>& y,
-                                             const std::vector<int>& y_index,
-                                             const YeVec& ye,
-                                             CovarFun&& covariance_function,
-                                             CovarArgs&& covar_args, RNG& rng,
-                                             std::ostream* msgs) {
-  return laplace_base_rng(poisson_log_2_likelihood{},
-                          std::forward_as_tuple(y, y_index, ye),
-                          std::forward<CovarFun>(covariance_function),
-                          std::forward<CovarArgs>(covar_args),
-                          laplace_options_default{}, rng, msgs);
-}
-
-}  // namespace math
-}  // namespace stan
-
-#endif
diff --git a/stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp b/stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp
index b02830d3ec2..802f2dde8e4 100644
--- a/stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp
+++ b/stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp
@@ -19,19 +19,22 @@ namespace math {
  * In this specialized function, the likelihood p(y|theta) is a
  * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
  * with dynamic sized rows and 1 column.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @tparam RNG A valid boost rng type
  * @param[in] y Observed counts.
  * @param[in] y_index Index indicating which group each observation belongs to.
+ * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
  * \laplace_options
  * \rng_arg
  * \msg_arg
  */
-template <typename ThetaVec, typename CovarFun, typename CovarArgs,
-          typename RNG, require_eigen_t<ThetaVec>* = nullptr>
+template <typename ThetaVec, typename Mean, typename CovarFun,
+          typename CovarArgs, typename RNG,
+          require_eigen_vector_t<ThetaVec>* = nullptr>
 inline Eigen::VectorXd laplace_latent_tol_poisson_log_rng(
-    const std::vector<int>& y, const std::vector<int>& y_index,
+    const std::vector<int>& y, const std::vector<int>& y_index, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args, ThetaVec&& theta_0,
     const double tolerance, const int max_num_steps,
     const int hessian_block_size, const int solver,
@@ -39,10 +42,11 @@ inline Eigen::VectorXd laplace_latent_tol_poisson_log_rng(
   laplace_options_user_supplied ops{hessian_block_size,    solver,
                                     max_steps_line_search, tolerance,
                                     max_num_steps,         value_of(theta_0)};
-  return laplace_base_rng(poisson_log_likelihood{},
-                          std::forward_as_tuple(y, y_index),
-                          std::forward<CovarFun>(covariance_function),
-                          std::forward<CovarArgs>(covar_args), ops, rng, msgs);
+  return laplace_base_rng(
+      poisson_log_likelihood{},
+      std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
+      std::forward<CovarFun>(covariance_function),
+      std::forward<CovarArgs>(covar_args), ops, rng, msgs);
 }
 
 /**
@@ -55,24 +59,27 @@ inline Eigen::VectorXd laplace_latent_tol_poisson_log_rng(
  * The Laplace approximation is computed using a Newton solver.
  * In this specialized function, the likelihood p(y|theta) is a
  * Poisson with a log link.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @tparam RNG A valid boost rng type
  * @param[in] y Observed counts.
  * @param[in] y_index Index indicating which group each observation belongs to.
+ * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
  * \rng_arg
  * \msg_arg
  */
-template <typename CovarFun, typename CovarArgs, typename RNG>
+template <typename CovarFun, typename CovarArgs, typename RNG, typename Mean>
 inline Eigen::VectorXd laplace_latent_poisson_log_rng(
-    const std::vector<int>& y, const std::vector<int>& y_index,
+    const std::vector<int>& y, const std::vector<int>& y_index, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args, RNG& rng,
     std::ostream* msgs) {
-  return laplace_base_rng(poisson_log_likelihood{},
-                          std::forward_as_tuple(y, y_index),
-                          std::forward<CovarFun>(covariance_function),
-                          std::forward<CovarArgs>(covar_args),
-                          laplace_options_default{}, rng, msgs);
+  return laplace_base_rng(
+      poisson_log_likelihood{},
+      std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
+      std::forward<CovarFun>(covariance_function),
+      std::forward<CovarArgs>(covar_args), laplace_options_default{}, rng,
+      msgs);
 }
 
 }  // namespace math
diff --git a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
index c9a28488009..c0af114da54 100644
--- a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
@@ -20,12 +20,14 @@ namespace stan {
 namespace math {
 
 struct bernoulli_logit_likelihood {
-  template <typename T_theta, typename YVec>
-  inline auto operator()(const T_theta& theta, const YVec& y,
-                         const std::vector<int>& delta_int,
+  template <typename ThetaVec, typename YVec, typename Mean>
+  inline auto operator()(const ThetaVec& theta, const YVec& y,
+                         const std::vector<int>& delta_int, Mean&& mean,
                          std::ostream* pstream) const {
-    return sum(elt_multiply(theta, y)
-               - elt_multiply(to_vector(delta_int), log(add(1.0, exp(theta)))));
+    auto theta_offset = to_ref(add(theta, mean));
+    return sum(
+        elt_multiply(theta_offset, y)
+        - elt_multiply(to_vector(delta_int), log(add(1.0, exp(theta_offset)))));
   }
 };
 
@@ -39,18 +41,21 @@ struct bernoulli_logit_likelihood {
  * @tparam propto boolean ignored
  * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
  * with dynamic sized rows and 1 column.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @param[in] y total counts per group. Second sufficient statistics.
  * @param[in] n_samples number of samples per group. First sufficient
- *            statistics.
+ * statistics.
+ * @param[in] mean the mean of the latent normal variable.
  * \laplace_common_args
  * \laplace_options
  * \msg_arg
  */
-template <bool propto = false, typename ThetaVec, typename CovarFun,
-          typename CovarArgs, require_eigen_t<ThetaVec>* = nullptr>
+template <bool propto = false, typename ThetaVec, typename Mean,
+          typename CovarFun, typename CovarArgs,
+          require_eigen_vector_t<ThetaVec>* = nullptr>
 inline auto laplace_marginal_tol_bernoulli_logit_lpmf(
-    const std::vector<int>& y, const std::vector<int>& n_samples,
+    const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args,
     const ThetaVec& theta_0, double tolerance, int max_num_steps,
     const int hessian_block_size, const int solver,
@@ -60,7 +65,7 @@ inline auto laplace_marginal_tol_bernoulli_logit_lpmf(
                                     max_num_steps,         value_of(theta_0)};
   return laplace_marginal_density(
       bernoulli_logit_likelihood{},
-      std::forward_as_tuple(to_vector(y), n_samples),
+      std::forward_as_tuple(to_vector(y), n_samples, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), ops, msgs);
 }
@@ -73,21 +78,24 @@ inline auto laplace_marginal_tol_bernoulli_logit_lpmf(
  * for more details.
  *
  * @tparam propto boolean ignored
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @param[in] y total counts per group. Second sufficient statistics.
  * @param[in] n_samples number of samples per group. First sufficient
- *            statistics.
+ * statistics.
+ * @param[in] mean the mean of the latent normal variable.
  * \laplace_common_args
  * \msg_arg
  */
-template <bool propto = false, typename CovarFun, typename CovarArgs>
+template <bool propto = false, typename Mean, typename CovarFun,
+          typename CovarArgs>
 inline auto laplace_marginal_bernoulli_logit_lpmf(
-    const std::vector<int>& y, const std::vector<int>& n_samples,
+    const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args,
     std::ostream* msgs) {
   return laplace_marginal_density(
       bernoulli_logit_likelihood{},
-      std::forward_as_tuple(to_vector(y), n_samples),
+      std::forward_as_tuple(to_vector(y), n_samples, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), laplace_options_default{}, msgs);
 }
diff --git a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
index 4ee58edac91..9da39279b52 100644
--- a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
@@ -24,24 +24,29 @@ namespace stan {
 namespace math {
 
 struct neg_binomial_2_log_likelihood {
-  template <typename T_theta, typename T_eta>
-  inline return_type_t<T_theta, T_eta> operator()(
-      const T_theta& theta, const T_eta& eta, const std::vector<int>& y,
-      const std::vector<int>& y_index, std::ostream* pstream) const {
+  template <typename ThetaVec, typename Eta, typename Mean,
+            require_all_eigen_vector_t<ThetaVec>* = nullptr>
+  inline auto operator()(const ThetaVec& theta, const Eta& eta,
+                         const std::vector<int>& y,
+                         const std::vector<int>& y_index, Mean&& mean,
+                         std::ostream* pstream) const {
     Eigen::VectorXi n_per_group = Eigen::VectorXi::Zero(theta.size());
     Eigen::VectorXi counts_per_group = Eigen::VectorXi::Zero(theta.size());
 
     for (int i = 0; i < y.size(); i++) {
-      n_per_group[y_index[i]]++;
-      counts_per_group[y_index[i]] += y[i];
+      n_per_group[y_index[i] - 1]++;
+      counts_per_group[y_index[i] - 1] += y[i];
     }
     Eigen::Map<const Eigen::VectorXi> y_map(y.data(), y.size());
-    auto log_eta_plus_exp_theta = eval(log(add(eta, exp(theta))));
+
+    auto theta_offset = to_ref(add(theta, mean));
+    auto log_eta_plus_exp_theta = eval(log(add(eta, exp(theta_offset))));
     return sum(binomial_coefficient_log(subtract(add(y_map, eta), 1), y_map))
-           + sum(add(elt_multiply(counts_per_group,
-                                  subtract(theta, log_eta_plus_exp_theta)),
-                     elt_multiply(multiply(n_per_group, eta),
-                                  subtract(log(eta), log_eta_plus_exp_theta))));
+           + sum(
+               add(elt_multiply(counts_per_group,
+                                subtract(theta_offset, log_eta_plus_exp_theta)),
+                   elt_multiply(multiply(n_per_group, eta),
+                                subtract(log(eta), log_eta_plus_exp_theta))));
   }
 };
 
@@ -55,21 +60,23 @@ struct neg_binomial_2_log_likelihood {
  * @tparam Eta The type of parameter arguments for the likelihood function.
  * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
  * with dynamic sized rows and 1 column.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @param[in] y observed counts.
  * @param[in] y_index group to which each observation belongs. Each group
  *            is parameterized by one element of theta.
  * @param[in] eta non-marginalized model parameters for the likelihood.
+ * @param[in] mean the mean of the latent normal variable
  * \laplace_common_args
  * \laplace_options
  * \msg_arg
  */
-template <bool propto = false, typename Eta, typename ThetaVec,
+template <bool propto = false, typename Eta, typename ThetaVec, typename Mean,
           typename CovarFun, typename CovarArgs,
-          require_all_eigen_vector_t<ThetaVec>* = nullptr>
+          require_eigen_vector_t<ThetaVec>* = nullptr>
 inline auto laplace_marginal_tol_neg_binomial_2_log_lpmf(
     const std::vector<int>& y, const std::vector<int>& y_index, const Eta& eta,
-    CovarFun&& covariance_function, CovarArgs&& covar_args,
+    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
     const ThetaVec& theta_0, double tolerance, int max_num_steps,
     const int hessian_block_size, const int solver,
     const int max_steps_line_search, std::ostream* msgs) {
@@ -77,7 +84,8 @@ inline auto laplace_marginal_tol_neg_binomial_2_log_lpmf(
                                     max_steps_line_search, tolerance,
                                     max_num_steps,         value_of(theta_0)};
   return laplace_marginal_density(
-      neg_binomial_2_log_likelihood{}, std::forward_as_tuple(eta, y, y_index),
+      neg_binomial_2_log_likelihood{},
+      std::forward_as_tuple(eta, y, y_index, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), ops, msgs);
 }
@@ -91,42 +99,51 @@ inline auto laplace_marginal_tol_neg_binomial_2_log_lpmf(
  *
  * @tparam Eta The type of parameter arguments for the likelihood function.
  * \laplace_common_template_args
+ * @tparam Mean type of the mean of the latent normal distribution
  * @param[in] y observed counts.
  * @param[in] y_index group to which each observation belongs. Each group
  *            is parameterized by one element of theta.
  * @param[in] eta Parameter argument for likelihood function.
+ * @param[in] mean the mean of the latent normal variable
  * \laplace_common_args
  * \msg_arg
  */
-template <bool propto = false, typename Eta, typename CovarFun,
+template <bool propto = false, typename Eta, typename Mean, typename CovarFun,
           typename CovarArgs>
 inline auto laplace_marginal_neg_binomial_2_log_lpmf(
     const std::vector<int>& y, const std::vector<int>& y_index, const Eta& eta,
-    CovarFun&& covariance_function, CovarArgs&& covar_args,
+    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
     std::ostream* msgs) {
   return laplace_marginal_density(
-      neg_binomial_2_log_likelihood{}, std::forward_as_tuple(eta, y, y_index),
+      neg_binomial_2_log_likelihood{},
+      std::forward_as_tuple(eta, y, y_index, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), laplace_options_default{}, msgs);
 }
 
 struct neg_binomial_2_log_likelihood_summary {
-  template <typename T_theta, typename T_eta>
-  inline return_type_t<T_theta, T_eta> operator()(
-      const T_theta& theta, const T_eta& eta, const std::vector<int>& y,
-      const std::vector<int>& n_per_group,
-      const std::vector<int>& counts_per_group, std::ostream* pstream) const {
+  template <typename ThetaVec, typename Eta, typename Mean,
+            require_eigen_vector_t<ThetaVec>* = nullptr>
+  inline auto operator()(const ThetaVec& theta, const Eta& eta,
+                         const std::vector<int>& y,
+                         const std::vector<int>& n_per_group,
+                         const std::vector<int>& counts_per_group, Mean&& mean,
+                         std::ostream* pstream) const {
     Eigen::Map<const Eigen::VectorXi> y_map(y.data(), y.size());
     Eigen::Map<const Eigen::VectorXi> n_per_group_map(n_per_group.data(),
                                                       n_per_group.size());
     Eigen::Map<const Eigen::VectorXi> counts_per_group_map(
         counts_per_group.data(), counts_per_group.size());
-    auto log_eta_plus_exp_theta = eval(log(add(eta, exp(theta))));
+
+    auto theta_offset = to_ref(add(theta, mean));
+    auto log_eta_plus_exp_theta = eval(log(add(eta, exp(theta_offset))));
+
     return sum(binomial_coefficient_log(subtract(add(y_map, eta), 1.0), y_map))
-           + sum(add(elt_multiply(counts_per_group_map,
-                                  subtract(theta, log_eta_plus_exp_theta)),
-                     elt_multiply(multiply(n_per_group_map, eta),
-                                  subtract(log(eta), log_eta_plus_exp_theta))));
+           + sum(
+               add(elt_multiply(counts_per_group_map,
+                                subtract(theta_offset, log_eta_plus_exp_theta)),
+                   elt_multiply(multiply(n_per_group_map, eta),
+                                subtract(log(eta), log_eta_plus_exp_theta))));
   }
 };
 
@@ -140,21 +157,23 @@ struct neg_binomial_2_log_likelihood_summary {
  * @tparam Eta The type of parameter arguments for the likelihood function.
  * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
  * with dynamic sized rows and 1 column.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @param[in] y observations.
  * @param[in] n_per_group number of samples per group
  * @param[in] counts_per_group total counts per group
  * @param[in] eta non-marginalized model parameters for the likelihood.
+ * @param[in] mean the mean of the latent normal variable
  * \laplace_common_args
  * \laplace_options
  * \msg_arg
  */
-template <bool propto = false, typename Eta, typename ThetaVec,
+template <bool propto = false, typename Eta, typename ThetaVec, typename Mean,
           typename CovarFun, typename CovarArgs,
-          require_all_eigen_vector_t<ThetaVec>* = nullptr>
+          require_eigen_vector_t<ThetaVec>* = nullptr>
 inline auto laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
     const std::vector<int>& y, const std::vector<int>& n_per_group,
-    const std::vector<int>& counts_per_group, const Eta& eta,
+    const std::vector<int>& counts_per_group, const Eta& eta, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args,
     const ThetaVec& theta_0, double tolerance, int max_num_steps,
     const int hessian_block_size, const int solver,
@@ -164,7 +183,8 @@ inline auto laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
                                     max_num_steps,         value_of(theta_0)};
   return laplace_marginal_density(
       neg_binomial_2_log_likelihood_summary{},
-      std::forward_as_tuple(eta, y, n_per_group, counts_per_group),
+      std::forward_as_tuple(eta, y, n_per_group, counts_per_group,
+                            std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), ops, msgs);
 }
@@ -177,24 +197,27 @@ inline auto laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
  * See the laplace_marginal function for more details.
  *
  * @tparam Eta The type of parameter arguments for the likelihood function.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @param[in] y observations.
  * @param[in] n_per_group number of samples per group
  * @param[in] counts_per_group total counts per group
  * @param[in] eta non-marginalized model parameters for the likelihood.
+ * @param[in] mean the mean of the latent normal variable
  * \laplace_common_args
  * \msg_arg
  */
-template <bool propto = false, typename Eta, typename CovarFun,
+template <bool propto = false, typename Eta, typename Mean, typename CovarFun,
           typename CovarArgs>
 inline auto laplace_marginal_neg_binomial_2_log_summary_lpmf(
     const std::vector<int>& y, const std::vector<int>& n_per_group,
-    const std::vector<int>& counts_per_group, const Eta& eta,
+    const std::vector<int>& counts_per_group, const Eta& eta, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args,
     std::ostream* msgs) {
   return laplace_marginal_density(
       neg_binomial_2_log_likelihood_summary{},
-      std::forward_as_tuple(eta, y, n_per_group, counts_per_group),
+      std::forward_as_tuple(eta, y, n_per_group, counts_per_group,
+                            std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
       std::forward<CovarArgs>(covar_args), laplace_options_default{}, msgs);
 }
diff --git a/stan/math/mix/prob/laplace_marginal_poisson_log_2_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_poisson_log_2_lpmf.hpp
deleted file mode 100644
index 25f25ae6243..00000000000
--- a/stan/math/mix/prob/laplace_marginal_poisson_log_2_lpmf.hpp
+++ /dev/null
@@ -1,125 +0,0 @@
-#ifndef STAN_MATH_MIX_PROB_LAPLACE_MARGINAL_POISSON_LOG_2_LPMF_HPP
-#define STAN_MATH_MIX_PROB_LAPLACE_MARGINAL_POISSON_LOG_2_LPMF_HPP
-
-#include <stan/math/mix/functor/laplace_marginal_density.hpp>
-#include <stan/math/mix/functor/laplace_likelihood.hpp>
-#include <stan/math/rev/core/operator_addition.hpp>
-#include <stan/math/rev/core/operator_multiplication.hpp>
-#include <stan/math/rev/core/operator_subtraction.hpp>
-#include <stan/math/rev/fun/dot_product.hpp>
-#include <stan/math/rev/fun/exp.hpp>
-#include <stan/math/rev/fun/lgamma.hpp>
-#include <stan/math/rev/fun/log.hpp>
-#include <stan/math/rev/fun/sum.hpp>
-#include <stan/math/fwd/fun/exp.hpp>
-#include <stan/math/fwd/fun/lgamma.hpp>
-#include <stan/math/fwd/fun/sum.hpp>
-
-namespace stan {
-namespace math {
-
-struct poisson_log_2_likelihood {
-  /**
-   * Returns the lpmf for a Poisson with a log link across
-   * multiple groups. No need to compute the log normalizing constant.
-   * Same as above, but includes a exposure term to correct the
-   * log rate for each group.
-   * @tparam Theta Type of the log Poisson rate.
-   * @tparam Eta Type of the auxiliary parameter (not used here).
-   * @param[in] theta log Poisson rate for each group.
-   * @param[in] y First n elements contain the sum of counts in each group
-   * @param[in] y_index group to which each observation belongs.
-   * @param[in] ye next n elements the exposure in each group, where n is the
-   * number of groups.
-   * @param[in, out] pstream msgs that are not used here
-   */
-  template <typename Theta, typename YVec, typename YIndexVec, typename YeVec>
-  inline auto operator()(const Theta& theta, const YVec& y,
-                         const YIndexVec& y_index, const YeVec& ye,
-                         std::ostream* /*pstream*/) const {
-    Eigen::VectorXd y_vec = to_vector(y);
-    Eigen::VectorXd counts_per_group = Eigen::VectorXd::Zero(theta.size());
-    Eigen::VectorXd n_per_group = Eigen::VectorXd::Zero(theta.size());
-    for (int i = 0; i < theta.size(); i++) {
-      counts_per_group(y_index[i]) += y[i];
-      n_per_group(y_index[i]) += 1;
-    }
-    // auto n_samples = to_vector(delta_int);
-    auto shifted_mean = to_ref(add(theta, log(ye)));
-    return -sum(lgamma(add(y_vec, 1))) + dot_product(shifted_mean, y_vec)
-           - dot_product(n_per_group, exp(shifted_mean));
-  }
-};
-
-/**
- * Wrapper function around the laplace_marginal function for
- * a log poisson likelihood. Returns the marginal density
- * p(y | phi) by marginalizing out the latent gaussian variable,
- * with a Laplace approximation. See the laplace_marginal function
- * for more details.
- *
- * @tparam propto boolean ignored
- * @tparam YeVec A type inheriting from `Eigen::EigenBase` with dynamic
- * sized rows and 1 column.
- * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
- * with dynamic sized rows and 1 column.
- * \laplace_common_template_args
- * @param[in] y total counts per group. Second sufficient statistics.
- * @param[in] y_index group to which each observation belongs.
- * @param[in] ye the exposure for each group.
- * \laplace_common_args
- * \laplace_options
- * \msg_arg
- */
-template <bool propto = false, typename YeVec, typename ThetaVec,
-          typename CovarFun, typename CovarArgs,
-          require_all_eigen_vector_t<YeVec, ThetaVec>* = nullptr>
-inline auto laplace_marginal_tol_poisson_2_log_lpmf(
-    const std::vector<int>& y, const std::vector<int>& y_index, const YeVec& ye,
-    CovarFun&& covariance_function, CovarArgs&& covar_args,
-    const ThetaVec& theta_0, double tolerance, int max_num_steps,
-    const int hessian_block_size, const int solver,
-    const int max_steps_line_search, std::ostream* msgs) {
-  laplace_options_user_supplied ops{hessian_block_size,    solver,
-                                    max_steps_line_search, tolerance,
-                                    max_num_steps,         value_of(theta_0)};
-  return laplace_marginal_density(
-      poisson_log_2_likelihood{}, std::forward_as_tuple(y, y_index, ye),
-      std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), ops, msgs);
-}
-
-/**
- * Wrapper function around the laplace_marginal function for
- * a log poisson likelihood. Returns the marginal density
- * p(y | phi) by marginalizing out the latent gaussian variable,
- * with a Laplace approximation. See the laplace_marginal function
- * for more details.
- *
- * @tparam propto boolean ignored
- * @tparam YeVec The type for the global parameter, phi.
- * \laplace_common_template_args
- * @param[in] y total counts per group. Second sufficient statistics.
- * @param[in] y_index group to which each observation belongs.
- * @param[in] ye the exposure for each group.
- * \laplace_common_args
- * \msg_arg
- */
-template <bool propto = false, typename YeVec, typename CovarFun,
-          typename CovarArgs, require_eigen_vector_t<YeVec>* = nullptr>
-inline auto laplace_marginal_poisson_2_log_lpmf(const std::vector<int>& y,
-                                                const std::vector<int>& y_index,
-                                                const YeVec& ye,
-                                                CovarFun&& covariance_function,
-                                                CovarArgs&& covar_args,
-                                                std::ostream* msgs) {
-  return laplace_marginal_density(
-      poisson_log_2_likelihood{}, std::forward_as_tuple(y, y_index, ye),
-      std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), laplace_options_default{}, msgs);
-}
-
-}  // namespace math
-}  // namespace stan
-
-#endif
diff --git a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
index f3ea5942cc8..5652bac7f01 100644
--- a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
@@ -17,30 +17,35 @@ struct poisson_log_likelihood {
   /**
    * Returns the lpmf for a Poisson with a log link across
    * multiple groups. No need to compute the log normalizing constant.
-   * @tparam T_theta Type of the log Poisson rate.
-   * @tparam T_eta Type of the auxiliary parameter (not used here).
+   * @tparam Theta A type inheriting from `Eigen::EigenBase` with dynamic
+   * sized rows and 1 column.
+   * @tparam YVec A vector type containing integers.
+   * @tparam Mean type of the mean of the latent normal distribution
    * @param[in] theta log Poisson rate for each group.
    * @param[in] y observed counts
    * @param[in] y_index group to which each observation belongs
    * return lpmf for a Poisson with a log link.
-   * @param[in] pstream
+   * @param[in] mean the mean of the latent normal variable
+   * \msg_arg
    */
-  template <typename Theta, typename YVec,
+  template <typename Theta, typename YVec, typename Mean,
             require_eigen_vector_t<Theta>* = nullptr>
   inline auto operator()(const Theta& theta, const YVec& y,
-                         const std::vector<int>& y_index,
-                         std::ostream* pstream) const {
+                         const std::vector<int>& y_index, Mean&& mean,
+                         std::ostream* /*pstream*/) const {
     Eigen::VectorXd counts_per_group = Eigen::VectorXd::Zero(theta.size());
     Eigen::VectorXd n_per_group = Eigen::VectorXd::Zero(theta.size());
 
     for (int i = 0; i < theta.size(); i++) {
-      counts_per_group(y_index[i]) += y[i];
-      n_per_group(y_index[i]) += 1;
+      counts_per_group(y_index[i] - 1) += y[i];
+      n_per_group(y_index[i] - 1) += 1;
     }
 
+    auto theta_offset = to_ref(add(theta, mean));
+
     return -sum(lgamma(add(counts_per_group, 1)))
-           + dot_product(theta, counts_per_group)
-           - dot_product(n_per_group, exp(theta));
+           + dot_product(theta_offset, counts_per_group)
+           - dot_product(n_per_group, exp(theta_offset));
   }
 };
 
@@ -54,17 +59,20 @@ struct poisson_log_likelihood {
  * @tparam propto ignored
  * @tparam ThetaVec A type inheriting from `Eigen::EigenBase`
  * with dynamic sized rows and 1 column.
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @param[in] y observed counts
  * @param[in] y_index group to which each observation belongs
+ * @param[in] mean the mean of the latent normal variable
  * \laplace_common_args
  * \laplace_options
  * \msg_arg
  */
-template <bool propto = false, typename ThetaVec, typename CovarFun,
-          typename CovarArgs, require_all_eigen_vector_t<ThetaVec>* = nullptr>
+template <bool propto = false, typename ThetaVec, typename Mean,
+          typename CovarFun, typename CovarArgs,
+          require_eigen_vector_t<ThetaVec>* = nullptr>
 inline auto laplace_marginal_tol_poisson_log_lpmf(
-    const std::vector<int>& y, const std::vector<int>& y_index,
+    const std::vector<int>& y, const std::vector<int>& y_index, Mean&& mean,
     CovarFun&& covariance_function, CovarArgs&& covar_args,
     const ThetaVec& theta_0, double tolerance, int max_num_steps,
     const int hessian_block_size, const int solver,
@@ -73,7 +81,8 @@ inline auto laplace_marginal_tol_poisson_log_lpmf(
                                     max_steps_line_search, tolerance,
                                     max_num_steps,         value_of(theta_0)};
   return laplace_marginal_density(
-      poisson_log_likelihood{}, std::forward_as_tuple(y, y_index),
+      poisson_log_likelihood{},
+      std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
       covariance_function, std::forward<CovarArgs>(covar_args), ops, msgs);
 }
 
@@ -85,20 +94,25 @@ inline auto laplace_marginal_tol_poisson_log_lpmf(
  * for more details.
  *
  * @tparam propto ignored
+ * @tparam Mean type of the mean of the latent normal distribution
  * \laplace_common_template_args
  * @param[in] y observed counts
  * @param[in] y_index group to which each observation belongs
+ * @param[in] mean the mean of the latent normal variable
  * \laplace_common_args
  * \msg_arg
  */
-template <bool propto = false, typename CovarFun, typename CovarArgs>
+template <bool propto = false, typename CovarFun, typename CovarArgs,
+          typename Mean>
 inline auto laplace_marginal_poisson_log_lpmf(const std::vector<int>& y,
                                               const std::vector<int>& y_index,
+                                              Mean&& mean,
                                               CovarFun&& covariance_function,
                                               CovarArgs&& covar_args,
                                               std::ostream* msgs) {
   return laplace_marginal_density(
-      poisson_log_likelihood{}, std::forward_as_tuple(y, y_index),
+      poisson_log_likelihood{},
+      std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
       covariance_function, std::forward<CovarArgs>(covar_args),
       laplace_options_default{}, msgs);
 }
diff --git a/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp b/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
index 3b4af5c63d7..9ad613894c2 100644
--- a/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
@@ -6,8 +6,6 @@
 
 #include <gtest/gtest.h>
 #include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 namespace {
@@ -74,13 +72,15 @@ TEST(laplace_bernoulli_logit_rng, two_dim_diag) {
   std::vector<int> sums = {1, 0};
   Eigen::VectorXd ye(2);
   ye << 1, 1;
+  Eigen::VectorXd mean(2);
+  mean << 0, 0;
   std::vector<double> d0;
   std::vector<int> di0;
   std::vector<Eigen::VectorXd> x_dummy;
   boost::random::mt19937 rng;
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_bernoulli_logit_rng(
-      sums, n_samples, diagonal_kernel_functor{},
+      sums, n_samples, mean, diagonal_kernel_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   // Compute exact mean and covariance
diff --git a/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp
index 9f1b7a178e9..3445a6d1591 100644
--- a/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp
@@ -7,9 +7,6 @@
 #include <test/unit/math/rev/fun/util.hpp>
 
 #include <gtest/gtest.h>
-#include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 TEST(laplace_marginal_bernoulli_logit_lpmf, phi_dim500) {
@@ -33,6 +30,7 @@ TEST(laplace_marginal_bernoulli_logit_lpmf, phi_dim500) {
   }
   std::vector<int> n_samples = stan::math::rep_array(1, dim_theta);
   Eigen::VectorXd theta_0 = Eigen::VectorXd::Zero(dim_theta);
+  Eigen::VectorXd mean = Eigen::VectorXd::Zero(dim_theta);
   std::vector<double> delta;
   std::vector<int> delta_int;
   int dim_phi = 2;
@@ -41,7 +39,7 @@ TEST(laplace_marginal_bernoulli_logit_lpmf, phi_dim500) {
   phi_dbl << 1.6, 1;
   using stan::math::test::sqr_exp_kernel_functor;
   double target = laplace_marginal_bernoulli_logit_lpmf(
-      y, n_samples, sqr_exp_kernel_functor{},
+      y, n_samples, 0, sqr_exp_kernel_functor{},
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), nullptr);
   // Benchmark against gpstuff.
   EXPECT_NEAR(-195.368, target, tol);
@@ -56,7 +54,7 @@ TEST(laplace_marginal_bernoulli_logit_lpmf, phi_dim500) {
           auto&& theta_0) {
         auto f = [&](auto&& alpha, auto&& rho) {
           return laplace_marginal_tol_bernoulli_logit_lpmf(
-              y, n_samples, sqr_exp_kernel_functor{},
+              y, n_samples, mean, sqr_exp_kernel_functor{},
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
               max_steps_line_search, nullptr);
diff --git a/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp b/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp
index b0e193651c7..629d69128b7 100644
--- a/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp
@@ -8,8 +8,6 @@
 #include <test/unit/math/laplace/motorcycle_gp/x_vec.hpp>
 #include <gtest/gtest.h>
 #include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 struct poisson_log_likelihood2 {
@@ -35,7 +33,6 @@ TEST(laplace, poisson_log_phi_dim_2) {
   Eigen::VectorXd theta_0(dim_theta);
   theta_0 << 0, 0;
 
-  int dim_x = 2;
   std::vector<Eigen::VectorXd> x(dim_theta);
   Eigen::VectorXd x_0{{0.05100797, 0.16086164}};
   Eigen::VectorXd x_1{{-0.59823393, 0.98701425}};
@@ -176,7 +173,6 @@ TEST(laplace, bernoulli_logit_phi_dim500) {
   Eigen::Matrix<double, Eigen::Dynamic, 1> phi_dbl(dim_phi);
   phi_dbl << 1.6, 1;
 
-  stan::math::test::sqr_exp_kernel_functor K;
   double target = laplace_marginal<false>(
       bernoulli_logit_likelihood{}, std::forward_as_tuple(y),
       stan::math::test::sqr_exp_kernel_functor{},
diff --git a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp
index 66cf95d5bce..31151c2adf6 100644
--- a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp
@@ -6,9 +6,6 @@
 #include <test/unit/math/rev/fun/util.hpp>
 
 #include <gtest/gtest.h>
-#include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 TEST(laplace_marginal_beg_binomial_log_lpmf, phi_dim_2) {
@@ -18,14 +15,12 @@ TEST(laplace_marginal_beg_binomial_log_lpmf, phi_dim_2) {
   using stan::math::value_of;
   using stan::math::var;
 
-  int dim_phi = 2;
   double alpha_dbl = 1.6;
   double rho_dbl = 0.45;
   int dim_theta = 2;
   Eigen::VectorXd theta_0(dim_theta);
   theta_0 << 0, 0;
 
-  int dim_x = 2;
   std::vector<Eigen::VectorXd> x(dim_theta);
   Eigen::VectorXd x_0(2);
   x_0 << 0.05100797, 0.16086164;
@@ -38,7 +33,7 @@ TEST(laplace_marginal_beg_binomial_log_lpmf, phi_dim_2) {
   std::vector<int> delta_int;
 
   std::vector<int> y = {1, 0};
-  std::vector<int> y_index = {0, 1};
+  std::vector<int> y_index = {1, 2};
   double eta_dbl = 10000;
 
   constexpr double tolerance = 1e-12;
@@ -48,7 +43,7 @@ TEST(laplace_marginal_beg_binomial_log_lpmf, phi_dim_2) {
           auto&& theta_0) {
         auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
           return laplace_marginal_tol_neg_binomial_2_log_lpmf(
-              y, y_index, eta, stan::math::test::squared_kernel_functor{},
+              y, y_index, eta, 0, stan::math::test::squared_kernel_functor{},
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
               max_steps_line_search, nullptr);
@@ -68,7 +63,7 @@ TEST_F(laplace_disease_map_test, laplace_marginal_neg_binomial_2_log_lpmf) {
   double eta = 1;
 
   double marginal_density = laplace_marginal_neg_binomial_2_log_lpmf(
-      y, y_index, eta, stan::math::test::sqr_exp_kernel_functor(),
+      y, y_index, eta, mean, stan::math::test::sqr_exp_kernel_functor(),
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), nullptr);
 
   // ToDo (charlesm93): get benchmark from GPStuff or another software.
@@ -79,7 +74,7 @@ TEST_F(laplace_disease_map_test, laplace_marginal_neg_binomial_2_log_lpmf) {
           auto&& theta_0) {
         auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
           return laplace_marginal_tol_neg_binomial_2_log_lpmf(
-              y, y_index, eta, stan::math::test::sqr_exp_kernel_functor{},
+              y, y_index, eta, mean, stan::math::test::sqr_exp_kernel_functor{},
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
               max_steps_line_search, nullptr);
@@ -92,7 +87,7 @@ TEST_F(laplace_disease_map_test, laplace_marginal_neg_binomial_2_log_lpmf) {
           auto&& theta_0) {
         auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
           return laplace_marginal_tol_neg_binomial_2_log_lpmf(
-              y, y_index, eta, stan::math::test::sqr_exp_kernel_functor{},
+              y, y_index, eta, mean, stan::math::test::sqr_exp_kernel_functor{},
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
               max_steps_line_search, nullptr);
diff --git a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp
index 4420dedcc48..14c22f27767 100644
--- a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp
@@ -6,9 +6,6 @@
 #include <test/unit/math/rev/fun/util.hpp>
 
 #include <gtest/gtest.h>
-#include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 TEST(laplace_marginal_beg_binomial_log_summary_lpmf, phi_dim_2) {
@@ -18,14 +15,12 @@ TEST(laplace_marginal_beg_binomial_log_summary_lpmf, phi_dim_2) {
   using stan::math::value_of;
   using stan::math::var;
 
-  int dim_phi = 2;
   double alpha_dbl = 1.6;
   double rho_dbl = 0.45;
   int dim_theta = 2;
   Eigen::VectorXd theta_0(dim_theta);
   theta_0 << 0, 0;
 
-  int dim_x = 2;
   std::vector<Eigen::VectorXd> x(dim_theta);
   Eigen::VectorXd x_0(2);
   x_0 << 0.05100797, 0.16086164;
@@ -38,13 +33,13 @@ TEST(laplace_marginal_beg_binomial_log_summary_lpmf, phi_dim_2) {
   std::vector<int> delta_int;
 
   std::vector<int> y = {1, 0};
-  std::vector<int> y_index = {0, 1};
+  std::vector<int> y_index = {1, 1};
   double eta_dbl = 10000;
   std::vector<int> n_per_group(theta_0.size(), 0);
   std::vector<int> counts_per_group(theta_0.size(), 0);
   for (int i = 0; i < y.size(); i++) {
-    n_per_group[y_index[i]]++;
-    counts_per_group[y_index[i]] += y[i];
+    n_per_group[y_index[i] - 1]++;
+    counts_per_group[y_index[i] - 1] += y[i];
   }
   constexpr double tolerance = 1e-12;
   constexpr int max_num_steps = 1000;
@@ -53,7 +48,7 @@ TEST(laplace_marginal_beg_binomial_log_summary_lpmf, phi_dim_2) {
           auto&& theta_0) {
         auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
           return laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
-              y, n_per_group, counts_per_group, eta,
+              y, n_per_group, counts_per_group, eta, 0,
               stan::math::test::squared_kernel_functor{},
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
@@ -76,12 +71,12 @@ TEST_F(laplace_disease_map_test,
   std::vector<int> n_per_group(theta_0.size(), 0);
   std::vector<int> counts_per_group(theta_0.size(), 0);
   for (int i = 0; i < y.size(); i++) {
-    n_per_group[y_index[i]]++;
-    counts_per_group[y_index[i]] += y[i];
+    n_per_group[y_index[i] - 1]++;
+    counts_per_group[y_index[i] - 1] += y[i];
   }
 
   double marginal_density = laplace_marginal_neg_binomial_2_log_summary_lpmf(
-      y, n_per_group, counts_per_group, eta,
+      y, n_per_group, counts_per_group, eta, mean,
       stan::math::test::sqr_exp_kernel_functor{},
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), nullptr);
 
@@ -93,7 +88,7 @@ TEST_F(laplace_disease_map_test,
           auto&& theta_0) {
         auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
           return laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
-              y, n_per_group, counts_per_group, eta,
+              y, n_per_group, counts_per_group, eta, 0,
               stan::math::test::sqr_exp_kernel_functor{},
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
@@ -107,7 +102,7 @@ TEST_F(laplace_disease_map_test,
           auto&& theta_0) {
         auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
           return laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
-              y, n_per_group, counts_per_group, eta,
+              y, n_per_group, counts_per_group, eta, mean,
               stan::math::test::sqr_exp_kernel_functor{},
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
diff --git a/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp
index 8604f4b039d..7c9a14d0b9e 100644
--- a/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp
@@ -5,29 +5,23 @@
 #include <test/unit/math/rev/fun/util.hpp>
 
 #include <gtest/gtest.h>
-#include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 TEST(laplace_marginal_poisson_log_lpmf, phi_dim_2) {
-  using stan::math::laplace_marginal_poisson_2_log_lpmf;
   using stan::math::laplace_marginal_poisson_log_lpmf;
-  using stan::math::laplace_marginal_tol_poisson_2_log_lpmf;
   using stan::math::laplace_marginal_tol_poisson_log_lpmf;
 
+  using stan::math::log;
   using stan::math::to_vector;
   using stan::math::value_of;
   using stan::math::var;
 
-  int dim_phi = 2;
   double alpha_dbl = 1.6;
   double rho_dbl = 0.45;
   int dim_theta = 2;
   Eigen::VectorXd theta_0(dim_theta);
   theta_0 << 0, 0;
 
-  int dim_x = 2;
   std::vector<Eigen::VectorXd> x(dim_theta);
   Eigen::VectorXd x_0(2);
   x_0 << 0.05100797, 0.16086164;
@@ -40,7 +34,7 @@ TEST(laplace_marginal_poisson_log_lpmf, phi_dim_2) {
   std::vector<int> delta_int;
 
   std::vector<int> y = {1, 0};
-  std::vector<int> y_index = {0, 1};
+  std::vector<int> y_index = {1, 2};
 
   stan::math::test::squared_kernel_functor sq_kernel;
   constexpr double tolerance = 1e-6;
@@ -57,7 +51,7 @@ TEST(laplace_marginal_poisson_log_lpmf, phi_dim_2) {
       for (int solver_num = 1; solver_num < 4; solver_num++) {
         auto f = [&](auto&& alpha, auto&& rho) {
           return laplace_marginal_tol_poisson_log_lpmf(
-              y, y_index, sq_kernel, std::forward_as_tuple(x, alpha, rho),
+              y, y_index, 0, sq_kernel, std::forward_as_tuple(x, alpha, rho),
               theta_0, tolerance, max_num_steps, hessian_block_size, solver_num,
               max_steps_line_search, nullptr);
         };
@@ -74,9 +68,10 @@ TEST(laplace_marginal_poisson_log_lpmf, phi_dim_2) {
          hessian_block_size++) {
       for (int solver_num = 1; solver_num < 4; solver_num++) {
         auto f = [&](auto&& alpha, auto&& rho) {
-          return laplace_marginal_tol_poisson_2_log_lpmf(
-              y, y_index, ye, sq_kernel, std::forward_as_tuple(x, alpha, rho),
-              theta_0, tolerance, max_num_steps, hessian_block_size, solver_num,
+          return laplace_marginal_tol_poisson_log_lpmf(
+              y, y_index, log(ye), sq_kernel,
+              std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
+              max_num_steps, hessian_block_size, solver_num,
               max_steps_line_search, nullptr);
         };
         stan::test::expect_ad<true>(tols, f, alpha_dbl, rho_dbl);
@@ -86,14 +81,14 @@ TEST(laplace_marginal_poisson_log_lpmf, phi_dim_2) {
 }
 
 TEST_F(laplace_disease_map_test, laplace_marginal_poisson_log_lpmf) {
-  using stan::math::laplace_marginal_poisson_2_log_lpmf;
   using stan::math::laplace_marginal_poisson_log_lpmf;
-  using stan::math::laplace_marginal_tol_poisson_2_log_lpmf;
+  using stan::math::laplace_marginal_tol_poisson_log_lpmf;
+  using stan::math::log;
   using stan::math::value_of;
   using stan::math::var;
 
-  double marginal_density = laplace_marginal_poisson_2_log_lpmf(
-      y, y_index, ye, stan::math::test::sqr_exp_kernel_functor(),
+  double marginal_density = laplace_marginal_poisson_log_lpmf(
+      y, y_index, log(ye), stan::math::test::sqr_exp_kernel_functor(),
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), nullptr);
 
   double tol = 6e-4;
@@ -108,8 +103,8 @@ TEST_F(laplace_disease_map_test, laplace_marginal_poisson_log_lpmf) {
          hessian_block_size++) {
       for (int solver_num = 1; solver_num < 4; solver_num++) {
         auto f = [&](auto&& alpha, auto&& rho) {
-          return laplace_marginal_tol_poisson_2_log_lpmf(
-              y, n_samples, ye, stan::math::test::sqr_exp_kernel_functor(),
+          return laplace_marginal_tol_poisson_log_lpmf(
+              y, y_index, log(ye), stan::math::test::sqr_exp_kernel_functor(),
               std::forward_as_tuple(x, alpha, rho), theta_0, tolerance,
               max_num_steps, hessian_block_size, solver_num,
               max_steps_line_search, nullptr);
@@ -119,3 +114,33 @@ TEST_F(laplace_disease_map_test, laplace_marginal_poisson_log_lpmf) {
     }
   }
 }
+
+struct diag_covariance {
+  template <typename T0__>
+  Eigen::Matrix<stan::return_type_t<T0__>, -1, -1> operator()(
+      const T0__& sigma, const int& N, std::ostream* pstream__) const {
+    return stan::math::diag_matrix(
+        stan::math::rep_vector(stan::math::pow(sigma, 2), N));
+  }
+};
+
+TEST(laplace_marginal_poisson_log_lpmf, mean_argument) {
+  // working example from
+  // https://discourse.mc-stan.org/t/embedded-laplace-numerical-problem/39700
+  using stan::math::laplace_marginal_poisson_log_lpmf;
+
+  const int N = 1;
+  const std::vector<int> y{153};
+  const std::vector<int> y_index{1};
+
+  Eigen::VectorXd mu(1);
+  mu << 4.3;
+
+  const double sigmaz = 2.0;
+
+  double marginal_density = laplace_marginal_poisson_log_lpmf(
+      y, y_index, mu, diag_covariance(), std::tuple<double, int>(sigmaz, N),
+      nullptr);
+
+  EXPECT_FLOAT_EQ(-6.7098737, marginal_density);
+}
diff --git a/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp b/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
index 98bf81c0670..4ace1fc9bdd 100644
--- a/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
@@ -7,8 +7,6 @@
 
 #include <gtest/gtest.h>
 #include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 namespace {
@@ -86,7 +84,7 @@ TEST(laplace_latent_neg_binomial_2_log_rng, count_two_dim_diag) {
   using stan::math::square;
 
   std::vector<int> y = {1, 0};
-  std::vector<int> y_index = {0, 1};
+  std::vector<int> y_index = {1, 2};
   Eigen::VectorXd theta_0{{1, 1}};
   Eigen::VectorXd phi{{3, 2}};
 
@@ -105,7 +103,7 @@ TEST(laplace_latent_neg_binomial_2_log_rng, count_two_dim_diag) {
 
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_neg_binomial_2_log_rng(
-      y, y_index, eta, diagonal_kernel_nb_functor{},
+      y, y_index, eta, 0, diagonal_kernel_nb_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   double tol = 1e-3;
@@ -118,7 +116,7 @@ TEST(laplace_latent_neg_binomial_2_log_rng, count_two_dim_diag) {
   for (int i = 0; i < n_sim; i++) {
     rng.seed(2025 + i);
     Eigen::MatrixXd theta_pred = laplace_latent_neg_binomial_2_log_rng(
-        y, y_index, eta, diagonal_kernel_nb_functor{},
+        y, y_index, eta, 0, diagonal_kernel_nb_functor{},
         std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
     theta_dim0(i) = theta_pred(0);
diff --git a/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp b/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
index 8dc96726c41..c2d358a1586 100644
--- a/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
@@ -6,10 +6,6 @@
 #include <boost/math/distributions.hpp>
 
 #include <gtest/gtest.h>
-#include <iostream>
-#include <istream>
-#include <fstream>
-#include <vector>
 
 TEST_F(laplace_count_two_dim_diag_test, poisson_log_likelihood) {
   using stan::math::laplace_latent_poisson_log_rng;
@@ -27,7 +23,7 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_likelihood) {
   boost::random::mt19937 rng;
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_poisson_log_rng(
-      y, y_index, stan::math::test::diagonal_kernel_functor{},
+      y, y_index, 0, stan::math::test::diagonal_kernel_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   // double tol = 1e-3;
@@ -40,7 +36,7 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_likelihood) {
   for (int i = 0; i < n_sim; i++) {
     rng.seed(2025 + i);
     Eigen::MatrixXd theta_pred = laplace_latent_poisson_log_rng(
-        y, y_index, stan::math::test::diagonal_kernel_functor{},
+        y, y_index, 0, stan::math::test::diagonal_kernel_functor{},
         std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
     theta_dim0(i) = theta_pred(0);
@@ -69,14 +65,15 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_likelihood) {
 }
 
 TEST_F(laplace_count_two_dim_diag_test, poisson_log_exp_likelihood) {
-  using stan::math::laplace_latent_poisson_2_log_rng;
+  using stan::math::laplace_latent_poisson_log_rng;
+  using stan::math::log;
   using stan::math::multi_normal_rng;
   using stan::math::sqrt;
   using stan::math::square;
 
   rng.seed(1954);
-  Eigen::MatrixXd theta_pred_exp = laplace_latent_poisson_2_log_rng(
-      y, y_index, ye, stan::math::test::diagonal_kernel_functor{},
+  Eigen::MatrixXd theta_pred_exp = laplace_latent_poisson_log_rng(
+      y, y_index, log(ye), stan::math::test::diagonal_kernel_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   EXPECT_NEAR(theta_benchmark(0), theta_pred_exp(0), tol);
@@ -86,8 +83,8 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_exp_likelihood) {
   Eigen::VectorXd theta_dim1(n_sim);
   for (int i = 0; i < n_sim; i++) {
     rng.seed(2025 + i);
-    Eigen::MatrixXd theta_pred = laplace_latent_poisson_2_log_rng(
-        y, y_index, ye, stan::math::test::diagonal_kernel_functor{},
+    Eigen::MatrixXd theta_pred = laplace_latent_poisson_log_rng(
+        y, y_index, log(ye), stan::math::test::diagonal_kernel_functor{},
         std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
     theta_dim0(i) = theta_pred(0);
diff --git a/test/unit/math/laplace/laplace_types_test.cpp b/test/unit/math/laplace/laplace_types_test.cpp
index c994121c017..b2f9510d6b6 100644
--- a/test/unit/math/laplace/laplace_types_test.cpp
+++ b/test/unit/math/laplace/laplace_types_test.cpp
@@ -8,8 +8,6 @@
 #include <test/unit/math/laplace/motorcycle_gp/x_vec.hpp>
 #include <gtest/gtest.h>
 #include <iostream>
-#include <istream>
-#include <fstream>
 #include <vector>
 
 namespace {
@@ -49,12 +47,8 @@ struct poisson_re_log_ll {
   stan::return_type_t<stan::base_type_t<T0__>, stan::base_type_t<T2__>>
   operator()(const T0__& theta_arg__, const std::vector<int>& y,
              const T2__& mu_arg__, std::ostream* pstream__) const {
-    using local_scalar_t__
-        = stan::return_type_t<stan::base_type_t<T0__>, stan::base_type_t<T2__>>;
-    // suppress unused var warning
     const auto& theta = stan::math::to_ref(theta_arg__);
     const auto& mu = stan::math::to_ref(mu_arg__);
-    static constexpr bool propto__ = true;
     return stan::math::poisson_log_lpmf<false>(y, stan::math::add(mu, theta));
   }
 };
@@ -66,7 +60,6 @@ struct cov_fun {
                 std::is_floating_point<std::decay_t<T0__>>>>* = nullptr>
   Eigen::Matrix<stan::return_type_t<T0__>, -1, -1> operator()(
       const T0__& sigma, const int& N, std::ostream* pstream__) const {
-    using local_scalar_t__ = stan::return_type_t<T0__>;
     return stan::math::diag_matrix(
         stan::math::rep_vector(stan::math::pow(sigma, 2), N));
   }
@@ -126,7 +119,6 @@ TEST(laplace, poisson_log_phi_dim_2_tuple_extended) {
   Eigen::VectorXd theta_0(dim_theta);
   theta_0 << 0, 0;
 
-  int dim_x = 2;
   std::vector<Eigen::VectorXd> x(dim_theta);
   Eigen::VectorXd x_0{{0.05100797, 0.16086164}};
   Eigen::VectorXd x_1{{-0.59823393, 0.98701425}};
@@ -181,7 +173,6 @@ TEST(laplace, poisson_log_phi_dim_2_tuple) {
   Eigen::VectorXd theta_0(dim_theta);
   theta_0 << 0, 0;
 
-  int dim_x = 2;
   std::vector<Eigen::VectorXd> x(dim_theta);
   Eigen::VectorXd x_0{{0.05100797, 0.16086164}};
   Eigen::VectorXd x_1{{-0.59823393, 0.98701425}};
@@ -259,7 +250,6 @@ TEST(laplace, poisson_log_phi_dim_2_array_tuple) {
   Eigen::VectorXd theta_0(dim_theta);
   theta_0 << 0, 0;
 
-  int dim_x = 2;
   std::vector<Eigen::VectorXd> x(dim_theta);
   Eigen::VectorXd x_0{{0.05100797, 0.16086164}};
   Eigen::VectorXd x_1{{-0.59823393, 0.98701425}};
diff --git a/test/unit/math/laplace/laplace_utility.hpp b/test/unit/math/laplace/laplace_utility.hpp
index 8f995452462..b60494c345d 100644
--- a/test/unit/math/laplace/laplace_utility.hpp
+++ b/test/unit/math/laplace/laplace_utility.hpp
@@ -4,8 +4,6 @@
 #include <test/unit/math/laplace/aki_disease_data/x1.hpp>
 #include <boost/algorithm/string.hpp>
 #include <iostream>
-#include <istream>
-#include <fstream>
 #include <gtest/gtest.h>
 
 namespace stan {
@@ -290,6 +288,7 @@ class laplace_disease_map_test : public ::testing::Test {
       n_samples[i] = 1;
 
     theta_0 = Eigen::VectorXd::Zero(dim_theta);
+    mean = Eigen::VectorXd::Zero(dim_theta);
     dim_phi = 2;
     phi_dbl.resize(dim_phi);
     phi_dbl << 0.3162278, 200;  // variance, length scale
@@ -299,7 +298,7 @@ class laplace_disease_map_test : public ::testing::Test {
     for (int i = 0; i < n_observations; i++) {
       delta_lk(i) = y[i];
       delta_lk(n_observations + i) = ye(i);
-      y_index[i] = i;
+      y_index[i] = i + 1;
     }
   }
 
@@ -317,6 +316,7 @@ class laplace_disease_map_test : public ::testing::Test {
   std::vector<int> delta_int;
 
   Eigen::VectorXd theta_0;
+  Eigen::VectorXd mean;
   int dim_phi;
   Eigen::Matrix<double, -1, 1> phi_dbl;
   Eigen::Matrix<double, -1, 1> eta_dummy_dbl;
@@ -333,7 +333,7 @@ class laplace_count_two_dim_diag_test : public ::testing::Test {
     y.resize(2);
     y = {1, 0};
     y_index.resize(2);
-    y_index = {0, 1};
+    y_index = {1, 2};
 
     theta_root = algebra_solver(stan::math::test::stationary_point(), theta_0,
                                 phi, d0, di0);