diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 704b5d2ddb4..d6f7be5c8a8 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -1,4 +1,4 @@
-name: Windows RTools 3.5
+name: Windows RTools 4.0
 
 on:
   pull_request:
@@ -22,19 +22,17 @@ jobs:
       with:
         token: ${{ secrets.GITHUB_TOKEN }}
       if: "!startsWith(github.ref, 'refs/tags/') && github.ref != 'refs/heads/master' && github.ref != 'refs/heads/develop'"
-  
+
     - uses: actions/checkout@v2
     - uses: actions/setup-python@v2
       with:
         python-version: '2.x'
-    - name: Download RTools
-      run: Invoke-WebRequest -Uri https://cran.rstudio.com/bin/windows/Rtools/Rtools35.exe -OutFile ./R35.exe
-    - name: Install RTools
-      shell: powershell
-      run: Start-Process -FilePath ./R35.exe -ArgumentList /VERYSILENT -NoNewWindow -Wait
-    - name: PATH Setup
-      shell: powershell
-      run: echo "C:/Rtools/bin;C:/Rtools/mingw_64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
+    - uses: r-lib/actions/setup-r@v1
+      with:
+        r-version: 'release'
+    - name: Set path for RTools 4.0
+      if: runner.os == 'Windows'
+      run: echo "C:/rtools40/usr/bin;C:/rtools40/mingw64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
     - name: Print g++ & make version and path
       shell: powershell
       run: |
@@ -56,7 +54,7 @@ jobs:
         python.exe runTests.py -j2 test/unit/math/prim
         python.exe runTests.py -j2 test/unit/math/rev
         python.exe runTests.py -j2 test/unit/math/memory
-        
+
     - name: Upload gtest_output xml
       uses: actions/upload-artifact@v2
       if: failure()
@@ -72,14 +70,12 @@ jobs:
     - uses: actions/setup-python@v2
       with:
         python-version: '2.x'
-    - name: Download RTools
-      run: Invoke-WebRequest -Uri https://cran.rstudio.com/bin/windows/Rtools/Rtools35.exe -OutFile ./R35.exe
-    - name: Install RTools
-      shell: powershell
-      run: Start-Process -FilePath ./R35.exe -ArgumentList /VERYSILENT -NoNewWindow -Wait
-    - name: PATH Setup
-      shell: powershell
-      run: echo "C:/Rtools/bin;C:/Rtools/mingw_64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
+    - uses: r-lib/actions/setup-r@v1
+      with:
+        r-version: 'release'
+    - name: Set path for RTools 4.0
+      if: runner.os == 'Windows'
+      run: echo "C:/rtools40/usr/bin;C:/rtools40/mingw64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
     - name: Print g++ & make version and path
       shell: powershell
       run: |
@@ -93,9 +89,6 @@ jobs:
     - name: Add TBB to PATH
       shell: powershell
       run: echo "D:/a/math/math/lib/tbb" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
-    - name: Disable running fwd/mix tests
-      shell: powershell
-      run: echo "CXXFLAGS+= -DSTAN_MATH_TESTS_REV_ONLY" | Out-File -Append -FilePath make/local -Encoding utf8      
     - name: Run fwd unit tests and all the mix tests except those in mix/fun
       shell: powershell
       run: |
@@ -121,14 +114,12 @@ jobs:
     - uses: actions/setup-python@v2
       with:
         python-version: '2.x'
-    - name: Download RTools
-      run: Invoke-WebRequest -Uri https://cran.rstudio.com/bin/windows/Rtools/Rtools35.exe -OutFile ./R35.exe
-    - name: Install RTools
-      shell: powershell
-      run: Start-Process -FilePath ./R35.exe -ArgumentList /VERYSILENT -NoNewWindow -Wait
-    - name: PATH Setup
-      shell: powershell
-      run: echo "C:/Rtools/bin;C:/Rtools/mingw_64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
+    - uses: r-lib/actions/setup-r@v1
+      with:
+        r-version: 'release'
+    - name: Set path for RTools 4.0
+      if: runner.os == 'Windows'
+      run: echo "C:/rtools40/usr/bin;C:/rtools40/mingw64/bin" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
     - name: Print g++ & make version and path
       shell: powershell
       run: |
@@ -142,14 +133,11 @@ jobs:
     - name: Add TBB to PATH
       shell: powershell
       run: echo "D:/a/math/math/lib/tbb" | Out-File -Append -FilePath $env:GITHUB_PATH -Encoding utf8
-    - name: Disable running fwd/mix tests
-      shell: powershell
-      run: echo "CXXFLAGS+= -DSTAN_MATH_TESTS_REV_ONLY" | Out-File -Append -FilePath make/local -Encoding utf8      
     - name: Run mix/fun unit tests
       shell: powershell
       run: |
         python.exe runTests.py test/unit/math/mix/fun
-        
+
     - name: Upload gtest_output xml
       uses: actions/upload-artifact@v2
       if: failure()
diff --git a/Jenkinsfile b/Jenkinsfile
index 1c10966c46a..6ce297c951f 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -256,7 +256,8 @@ pipeline {
             }
             steps {
                 unstash 'MathSetup'
-	            sh "echo CXXFLAGS += -fsanitize=address >> make/local"
+                sh "echo CXX=${CLANG_CXX} > make/local"
+	              sh "echo CXXFLAGS += -fsanitize=address >> make/local"
                 script {
                     if (isUnix()) {
                         runTests("test/unit", false)
diff --git a/make/compiler_flags b/make/compiler_flags
index b59b75b7f82..830bb6b22e5 100644
--- a/make/compiler_flags
+++ b/make/compiler_flags
@@ -116,7 +116,7 @@ CPPFLAGS_SUNDIALS ?= -DNO_FPRINTF_OUTPUT $(CPPFLAGS_OPTIM_SUNDIALS) $(CXXFLAGS_F
 
 
 ## setup compiler flags
-CXXFLAGS_LANG ?= -std=c++1y
+CXXFLAGS_LANG ?= -std=c++17
 #CXXFLAGS_BOOST ?=
 CXXFLAGS_SUNDIALS ?= -pipe $(CXXFLAGS_OPTIM_SUNDIALS) $(CPPFLAGS_FLTO_SUNDIALS)
 #CXXFLAGS_GTEST
diff --git a/stan/math/fwd/functor/operands_and_partials.hpp b/stan/math/fwd/functor/operands_and_partials.hpp
index bf8e087d521..e852c29f0ef 100644
--- a/stan/math/fwd/functor/operands_and_partials.hpp
+++ b/stan/math/fwd/functor/operands_and_partials.hpp
@@ -111,10 +111,10 @@ class operands_and_partials<Op1, Op2, Op3, Op4, Op5, fvar<Dx>> {
 namespace internal {
 
 // Vectorized Univariate
-template <typename Dx>
-class ops_partials_edge<Dx, std::vector<fvar<Dx>>> {
+// std::vector<fvar<Dx>>
+template <typename Dx, typename Op>
+class ops_partials_edge<Dx, Op, require_std_vector_vt<is_fvar, Op>> {
  public:
-  using Op = std::vector<fvar<Dx>>;
   using partials_t = Eigen::Matrix<Dx, -1, 1>;
   partials_t partials_;                       // For univariate use-cases
   broadcast_array<partials_t> partials_vec_;  // For multivariate
@@ -137,11 +137,12 @@ class ops_partials_edge<Dx, std::vector<fvar<Dx>>> {
   }
 };
 
-template <typename Dx, int R, int C>
-class ops_partials_edge<Dx, Eigen::Matrix<fvar<Dx>, R, C>> {
+// Eigen::Matrix<fvar<Dx>, R, C>
+template <typename Dx, typename Op>
+class ops_partials_edge<Dx, Op, require_eigen_vt<is_fvar, Op>> {
  public:
-  using partials_t = Eigen::Matrix<Dx, R, C>;
-  using Op = Eigen::Matrix<fvar<Dx>, R, C>;
+  using partials_t
+      = Eigen::Matrix<Dx, Op::RowsAtCompileTime, Op::ColsAtCompileTime>;
   partials_t partials_;                       // For univariate use-cases
   broadcast_array<partials_t> partials_vec_;  // For multivariate
   explicit ops_partials_edge(const Op& ops)
@@ -164,11 +165,15 @@ class ops_partials_edge<Dx, Eigen::Matrix<fvar<Dx>, R, C>> {
 };
 
 // Multivariate; vectors of eigen types
-template <typename Dx, int R, int C>
-class ops_partials_edge<Dx, std::vector<Eigen::Matrix<fvar<Dx>, R, C>>> {
+// std::vector<Eigen::Matrix<fvar<Dx>, R, C>>
+template <typename Dx, typename Op>
+class ops_partials_edge<
+    Dx, Op,
+    require_all_t<is_std_vector<Op>, is_eigen<value_type_t<Op>>,
+                  is_fvar<value_type_t<value_type_t<Op>>>>> {
  public:
-  using Op = std::vector<Eigen::Matrix<fvar<Dx>, R, C>>;
-  using partial_t = Eigen::Matrix<Dx, R, C>;
+  using partial_t = Eigen::Matrix<Dx, value_type_t<Op>::RowsAtCompileTime,
+                                  value_type_t<Op>::ColsAtCompileTime>;
   std::vector<partial_t> partials_vec_;
   explicit ops_partials_edge(const Op& ops)
       : partials_vec_(ops.size()), operands_(ops) {
@@ -193,10 +198,12 @@ class ops_partials_edge<Dx, std::vector<Eigen::Matrix<fvar<Dx>, R, C>>> {
   }
 };
 
-template <typename Dx>
-class ops_partials_edge<Dx, std::vector<std::vector<fvar<Dx>>>> {
+template <typename Dx, typename Op>
+class ops_partials_edge<
+    Dx, Op,
+    require_all_t<is_std_vector<Op>, is_std_vector<value_type_t<Op>>,
+                  is_fvar<value_type_t<value_type_t<Op>>>>> {
  public:
-  using Op = std::vector<std::vector<fvar<Dx>>>;
   using partial_t = std::vector<Dx>;
   std::vector<partial_t> partials_vec_;
   explicit ops_partials_edge(const Op& ops)
diff --git a/stan/math/opencl/kernel_generator/operation_cl.hpp b/stan/math/opencl/kernel_generator/operation_cl.hpp
index ef2d4977c97..5f0a4f6a918 100644
--- a/stan/math/opencl/kernel_generator/operation_cl.hpp
+++ b/stan/math/opencl/kernel_generator/operation_cl.hpp
@@ -372,7 +372,9 @@ class operation_cl : public operation_cl_base {
    * expression.
    * @return number of elements
    */
-  inline int size() const { return derived().rows() * derived().cols(); }
+  inline Eigen::Index size() const {
+    return derived().rows() * derived().cols();
+  }
 
   /**
    * Number of rows threads need to be launched for. For most expressions this
diff --git a/stan/math/opencl/matrix_cl.hpp b/stan/math/opencl/matrix_cl.hpp
index 6f92f6f62b6..caac68bd01d 100644
--- a/stan/math/opencl/matrix_cl.hpp
+++ b/stan/math/opencl/matrix_cl.hpp
@@ -64,7 +64,7 @@ class matrix_cl : public matrix_cl_base {
 
   int cols() const { return cols_; }
 
-  int size() const { return rows_ * cols_; }
+  Eigen::Index size() const { return rows_ * cols_; }
 
   const matrix_cl_view& view() const { return view_; }
 
@@ -283,7 +283,7 @@ class matrix_cl : public matrix_cl_base {
   matrix_cl(const int rows, const int cols,
             matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(rows), cols_(cols), view_(partial_view) {
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     cl::Context& ctx = opencl_context.context();
@@ -321,7 +321,7 @@ class matrix_cl : public matrix_cl_base {
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(A.rows()), cols_(A.cols()), view_(partial_view) {
     using Mat_type = std::decay_t<ref_type_for_opencl_t<Mat>>;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer_no_heap_if<
@@ -457,7 +457,7 @@ class matrix_cl : public matrix_cl_base {
       return *this;
     }
     this->wait_for_read_write_events();
-    if (size() != a.size()) {
+    if (this->size() != a.size()) {
       buffer_cl_ = cl::Buffer(opencl_context.context(), CL_MEM_READ_WRITE,
                               sizeof(T) * a.size());
     }
@@ -518,7 +518,7 @@ class matrix_cl : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(const T* A) {
     cl::Event transfer_event;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -538,7 +538,7 @@ class matrix_cl : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(T* A) {
     cl::Event transfer_event;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -577,7 +577,7 @@ class matrix_cl : public matrix_cl_base {
    */
   template <bool No_heap, typename U, std::enable_if_t<No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer(obj.data());
@@ -587,7 +587,7 @@ class matrix_cl : public matrix_cl_base {
   template <bool No_heap, typename U, std::enable_if_t<!No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
     using U_val = std::decay_t<ref_type_for_opencl_t<U>>;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     auto* obj_heap = new U_val(std::move(obj));
diff --git a/stan/math/opencl/prim/bernoulli_cdf.hpp b/stan/math/opencl/prim/bernoulli_cdf.hpp
index 219e5785050..38e987f3b4d 100644
--- a/stan/math/opencl/prim/bernoulli_cdf.hpp
+++ b/stan/math/opencl/prim/bernoulli_cdf.hpp
@@ -37,7 +37,7 @@ return_type_t<T_prob_cl> bernoulli_cdf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 1.0;
   }
diff --git a/stan/math/opencl/prim/bernoulli_lccdf.hpp b/stan/math/opencl/prim/bernoulli_lccdf.hpp
index 44530cbde89..d2fd582b775 100644
--- a/stan/math/opencl/prim/bernoulli_lccdf.hpp
+++ b/stan/math/opencl/prim/bernoulli_lccdf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_prob_cl> bernoulli_lccdf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/bernoulli_lcdf.hpp b/stan/math/opencl/prim/bernoulli_lcdf.hpp
index 2b3dafc778e..30e42538f43 100644
--- a/stan/math/opencl/prim/bernoulli_lcdf.hpp
+++ b/stan/math/opencl/prim/bernoulli_lcdf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_prob_cl> bernoulli_lcdf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
index 038e024e7b1..9de5422b293 100644
--- a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
@@ -66,13 +66,14 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
 
   if (N == 0) {
diff --git a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
index 8db546a1853..accc947ab17 100644
--- a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
@@ -36,7 +36,7 @@ return_type_t<T_prob_cl> bernoulli_logit_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/bernoulli_lpmf.hpp b/stan/math/opencl/prim/bernoulli_lpmf.hpp
index b6335842b60..ddbfb5e9074 100644
--- a/stan/math/opencl/prim/bernoulli_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_lpmf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_prob_cl> bernoulli_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
index 337daf326dd..baf4974f60d 100644
--- a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
@@ -61,10 +61,10 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
   static const char* function = "categorical_logit_glm_lpmf";
   if (is_y_vector) {
     check_size_match(function, "Rows of ", "x", N_instances, "size of ", "y",
-                     size(y));
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "beta", N_classes, "size of ",
-                   "alpha", size(alpha));
+                   "alpha", math::size(alpha));
   check_size_match(function, "Columns of ", "x", N_attributes, "Rows of",
                    "beta", beta.rows());
 
diff --git a/stan/math/opencl/prim/cumulative_sum.hpp b/stan/math/opencl/prim/cumulative_sum.hpp
index 6027052c0f5..b31a7ee3e04 100644
--- a/stan/math/opencl/prim/cumulative_sum.hpp
+++ b/stan/math/opencl/prim/cumulative_sum.hpp
@@ -38,8 +38,8 @@ inline auto cumulative_sum(T_vec&& v) {
       = opencl_kernels::cumulative_sum<T_scal>::kernel1.get_option(
           "LOCAL_SIZE_");
   const int work_groups = std::min(
-      (v.size() + local_size - 1) / local_size,
-      static_cast<int>(
+      static_cast<Eigen::Index>((v.size() + local_size - 1) / local_size),
+      static_cast<Eigen::Index>(
           opencl_context.device()[0].getInfo<CL_DEVICE_MAX_COMPUTE_UNITS>())
           * 16);
   const int local_size2
diff --git a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
index 4f0ac0b2aec..1bfec8258b6 100644
--- a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
@@ -82,17 +82,18 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_phi_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "phi",
-                     size(phi));
+                     math::size(phi));
   }
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (N == 0) {
     return 0;
diff --git a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
index b2210f94683..b999c25e535 100644
--- a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
+++ b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
@@ -74,17 +74,18 @@ normal_id_glm_lpdf(const T_y_cl& y, const T_x_cl& x, const T_alpha_cl& alpha,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_sigma_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "sigma",
-                     size(sigma));
+                     math::size(sigma));
   }
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (!include_summand<propto, T_y_cl, T_x_cl, T_alpha_cl, T_beta_cl,
                        T_sigma_cl>::value) {
diff --git a/stan/math/opencl/prim/num_elements.hpp b/stan/math/opencl/prim/num_elements.hpp
index a5d965290d3..546f5758e3f 100644
--- a/stan/math/opencl/prim/num_elements.hpp
+++ b/stan/math/opencl/prim/num_elements.hpp
@@ -3,6 +3,8 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/opencl/prim/size.hpp>
 
 namespace stan {
 namespace math {
@@ -16,7 +18,7 @@ namespace math {
 template <typename T,
           require_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr>
 size_t num_elements(const T& m) {
-  return size(m);
+  return math::size(m);
 }
 
 }  // namespace math
diff --git a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
index b09e48c990c..5f5912b2328 100644
--- a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
@@ -63,19 +63,19 @@ return_type_t<T_x, T_beta, T_cuts> ordered_logistic_glm_lpmf(
 
   const size_t N_instances = x.rows();
   const size_t N_attributes = x.cols();
-  const size_t N_classes = size(cuts) + 1;
+  const size_t N_classes = math::size(cuts) + 1;
 
   if (is_y_vector) {
     check_size_match(function, "Rows of ", "x", N_instances, "rows of ", "y",
-                     size(y));
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x", N_attributes, "Size of",
-                   "beta", size(beta));
+                   "beta", math::size(beta));
 
   const auto& cuts_val = eval(value_of(cuts));
   if (N_classes >= 2) {
-    auto cuts_head = block_zero_based(cuts_val, 0, 0, size(cuts) - 1, 1);
-    auto cuts_tail = block_zero_based(cuts_val, 1, 0, size(cuts) - 1, 1);
+    auto cuts_head = block_zero_based(cuts_val, 0, 0, math::size(cuts) - 1, 1);
+    auto cuts_tail = block_zero_based(cuts_val, 1, 0, math::size(cuts) - 1, 1);
     check_cl(function, "Cuts", cuts_head, "ordered and finite")
         = cuts_head < cuts_tail && isfinite(cuts_head) && isfinite(cuts_tail);
   } else {
diff --git a/stan/math/opencl/prim/ordered_logistic_lpmf.hpp b/stan/math/opencl/prim/ordered_logistic_lpmf.hpp
index e72276b05aa..4d01c155ef9 100644
--- a/stan/math/opencl/prim/ordered_logistic_lpmf.hpp
+++ b/stan/math/opencl/prim/ordered_logistic_lpmf.hpp
@@ -72,9 +72,9 @@ inline return_type_t<T_y_cl, T_loc_cl, T_cuts_cl> ordered_logistic_lpmf(
   constexpr bool is_y_vector = !is_stan_scalar<T_y_cl>::value;
   static const char* function = "ordered_logistic_lpmf(OpenCL)";
 
-  if (size(y) != 1) {
-    check_size_match(function, "Size of ", "y", size(y), "Size of", "lambda",
-                     size(lambda));
+  if (stan::math::size(y) != 1) {
+    check_size_match(function, "Size of ", "y", math::size(y), "Size of",
+                     "lambda", math::size(lambda));
   }
 
   int N_instances = max_size(y, lambda);
diff --git a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
index 541fe38d5be..c6b730e7644 100644
--- a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
@@ -66,13 +66,14 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (N == 0) {
     return 0;
diff --git a/stan/math/opencl/prim/poisson_log_lpmf.hpp b/stan/math/opencl/prim/poisson_log_lpmf.hpp
index 4d01d570441..31ed3dfaecd 100644
--- a/stan/math/opencl/prim/poisson_log_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_lpmf.hpp
@@ -39,7 +39,7 @@ return_type_t<T_log_rate_cl> poisson_log_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n, "Log rate parameter",
                          alpha);
-  const size_t N = is_n_vector ? size(n) : size(alpha);
+  const size_t N = is_n_vector ? math::size(n) : math::size(alpha);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/poisson_lpmf.hpp b/stan/math/opencl/prim/poisson_lpmf.hpp
index 6e7ef9dd2aa..2878e9d2849 100644
--- a/stan/math/opencl/prim/poisson_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_lpmf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_rate_cl> poisson_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n, "Rate parameter",
                          lambda);
-  const size_t N = is_n_vector ? size(n) : size(lambda);
+  const size_t N = is_n_vector ? math::size(n) : math::size(lambda);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/std_normal_cdf.hpp b/stan/math/opencl/prim/std_normal_cdf.hpp
index 23448b45747..739d496780b 100644
--- a/stan/math/opencl/prim/std_normal_cdf.hpp
+++ b/stan/math/opencl/prim/std_normal_cdf.hpp
@@ -29,7 +29,7 @@ return_type_t<T_y_cl> std_normal_cdf(const T_y_cl& y) {
   using std::isfinite;
   using std::isnan;
 
-  const size_t N = size(y);
+  const size_t N = math::size(y);
   if (N == 0) {
     return 1.0;
   }
diff --git a/stan/math/opencl/prim/std_normal_lccdf.hpp b/stan/math/opencl/prim/std_normal_lccdf.hpp
index c16138d944c..1e531fdaf0f 100644
--- a/stan/math/opencl/prim/std_normal_lccdf.hpp
+++ b/stan/math/opencl/prim/std_normal_lccdf.hpp
@@ -30,7 +30,7 @@ return_type_t<T_y_cl> std_normal_lccdf(const T_y_cl& y) {
   using std::isfinite;
   using std::isnan;
 
-  const size_t N = size(y);
+  const size_t N = math::size(y);
   if (N == 0) {
     return 1.0;
   }
diff --git a/stan/math/opencl/prim/std_normal_lcdf.hpp b/stan/math/opencl/prim/std_normal_lcdf.hpp
index ffee5704877..86966cf7a60 100644
--- a/stan/math/opencl/prim/std_normal_lcdf.hpp
+++ b/stan/math/opencl/prim/std_normal_lcdf.hpp
@@ -183,7 +183,7 @@ return_type_t<T_y_cl> std_normal_lcdf(const T_y_cl& y) {
   using std::isfinite;
   using std::isnan;
 
-  const size_t N = size(y);
+  const size_t N = math::size(y);
   if (N == 0) {
     return 1.0;
   }
diff --git a/stan/math/opencl/prim/std_normal_lpdf.hpp b/stan/math/opencl/prim/std_normal_lpdf.hpp
index b9782db683e..2bf4de436c1 100644
--- a/stan/math/opencl/prim/std_normal_lpdf.hpp
+++ b/stan/math/opencl/prim/std_normal_lpdf.hpp
@@ -35,7 +35,7 @@ inline return_type_t<T_y_cl> std_normal_lpdf(const T_y_cl& y) {
   using std::isfinite;
   using std::isnan;
 
-  const size_t N = size(y);
+  const size_t N = math::size(y);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/rev/operands_and_partials.hpp b/stan/math/opencl/rev/operands_and_partials.hpp
index d25c2a157a8..1102be4a265 100644
--- a/stan/math/opencl/rev/operands_and_partials.hpp
+++ b/stan/math/opencl/rev/operands_and_partials.hpp
@@ -30,7 +30,7 @@ class ops_partials_edge<double, var_value<Op>,
   template <typename, typename, typename, typename, typename, typename>
   friend class stan::math::operands_and_partials;
   var_value<Op> operands_;
-  static constexpr int size() noexcept { return 0; }
+  static constexpr Eigen::Index size() noexcept { return 0; }
   inline auto& operand() noexcept { return this->operands_; }
   inline auto& partial() noexcept { return this->partials_; }
 };
diff --git a/stan/math/opencl/rev/softmax.hpp b/stan/math/opencl/rev/softmax.hpp
index f61ba86963d..e031050f1ba 100644
--- a/stan/math/opencl/rev/softmax.hpp
+++ b/stan/math/opencl/rev/softmax.hpp
@@ -3,6 +3,7 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/opencl/prim/log_sum_exp.hpp>
+#include <stan/math/opencl/prim/dot_product.hpp>
 #include <stan/math/opencl/kernel_generator.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
diff --git a/stan/math/opencl/zeros_strict_tri.hpp b/stan/math/opencl/zeros_strict_tri.hpp
index 1c916c0e90c..3ff0b92842a 100644
--- a/stan/math/opencl/zeros_strict_tri.hpp
+++ b/stan/math/opencl/zeros_strict_tri.hpp
@@ -40,7 +40,7 @@ inline void matrix_cl<T>::zeros_strict_tri() try {
         "zeros_strict_tri", "matrix_view",
         "matrix_cl_view::Diagonal is not a valid template parameter value", "");
   }
-  if (size() == 0) {
+  if (this->size() == 0) {
     return;
   }
   this->view_ = both(this->view_, invert(matrix_view));
diff --git a/stan/math/prim/fun/max_size.hpp b/stan/math/prim/fun/max_size.hpp
index 9b0c8f2af8f..b23cb33c21e 100644
--- a/stan/math/prim/fun/max_size.hpp
+++ b/stan/math/prim/fun/max_size.hpp
@@ -16,8 +16,9 @@ namespace math {
  * @return the size of the largest input
  */
 template <typename T1, typename... Ts>
-inline size_t max_size(const T1& x1, const Ts&... xs) {
-  return std::max({stan::math::size(x1), stan::math::size(xs)...});
+inline constexpr Eigen::Index max_size(const T1& x1, const Ts&... xs) {
+  return std::max({static_cast<Eigen::Index>(stan::math::size(x1)),
+                   static_cast<Eigen::Index>(stan::math::size(xs))...});
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/offset_multiplier_constrain.hpp b/stan/math/prim/fun/offset_multiplier_constrain.hpp
index a1fd18ee228..189341b45c9 100644
--- a/stan/math/prim/fun/offset_multiplier_constrain.hpp
+++ b/stan/math/prim/fun/offset_multiplier_constrain.hpp
@@ -106,8 +106,8 @@ inline auto offset_multiplier_constrain(const T& x, const M& mu, const S& sigma,
   check_finite("offset_multiplier_constrain", "offset", value_of_rec(mu_ref));
   check_positive_finite("offset_multiplier_constrain", "multiplier",
                         value_of_rec(sigma_ref));
-  if (size(sigma_ref) == 1) {
-    lp += sum(multiply_log(size(x), sigma_ref));
+  if (stan::math::size(sigma_ref) == 1) {
+    lp += sum(multiply_log(stan::math::size(x), sigma_ref));
   } else {
     lp += sum(log(sigma_ref));
   }
diff --git a/stan/math/prim/fun/scalar_seq_view.hpp b/stan/math/prim/fun/scalar_seq_view.hpp
index 5e7a176b8fc..66547912fb5 100644
--- a/stan/math/prim/fun/scalar_seq_view.hpp
+++ b/stan/math/prim/fun/scalar_seq_view.hpp
@@ -31,19 +31,19 @@ class scalar_seq_view<C, require_eigen_vector_t<C>> {
    */
   inline auto operator[](size_t i) const { return c_.coeff(i); }
 
-  inline auto size() const noexcept { return c_.size(); }
+  inline Eigen::Index size() const noexcept { return c_.size(); }
 
   inline const value_type_t<C>* data() const noexcept { return c_.data(); }
   inline value_type_t<C>* data() noexcept { return c_.data(); }
 
   template <typename T = C, require_st_arithmetic<T>* = nullptr>
   inline decltype(auto) val(size_t i) const {
-    return c_.coeffRef(i);
+    return c_.coeff(i);
   }
 
   template <typename T = C, require_st_autodiff<T>* = nullptr>
   inline decltype(auto) val(size_t i) const {
-    return c_.coeffRef(i).val();
+    return c_.coeff(i).val();
   }
 
  private:
@@ -66,7 +66,7 @@ class scalar_seq_view<C, require_var_matrix_t<C>> {
   inline const auto* data() const noexcept { return c_.vi_; }
   inline auto* data() noexcept { return c_.vi_; }
 
-  inline auto size() const noexcept { return c_.size(); }
+  inline Eigen::Index size() const noexcept { return c_.size(); }
 
   template <typename T = C, require_st_autodiff<T>* = nullptr>
   inline auto val(size_t i) const {
@@ -95,7 +95,7 @@ class scalar_seq_view<C, require_std_vector_t<C>> {
    * @return the element at the specified position in the container
    */
   inline auto operator[](size_t i) const { return c_[i]; }
-  inline auto size() const noexcept { return c_.size(); }
+  inline Eigen::Index size() const noexcept { return c_.size(); }
   inline const auto* data() const noexcept { return c_.data(); }
 
   template <typename T = C, require_st_arithmetic<T>* = nullptr>
@@ -125,8 +125,9 @@ class scalar_seq_view<C, require_t<std::is_pointer<C>>> {
    * @return the element at the specified position in the container
    */
   inline auto operator[](size_t i) const { return c_[i]; }
-  inline auto size() const noexcept {
+  static constexpr Eigen::Index size() {
     static_assert(1, "Cannot Return Size of scalar_seq_view with pointer type");
+    return 1;
   }
   inline const auto* data() const noexcept { return &c_[0]; }
 
diff --git a/stan/math/prim/fun/size.hpp b/stan/math/prim/fun/size.hpp
index b5540638290..fb8bac5d9d0 100644
--- a/stan/math/prim/fun/size.hpp
+++ b/stan/math/prim/fun/size.hpp
@@ -14,8 +14,8 @@ namespace math {
  * that are always of length 1.
  */
 template <typename T, require_stan_scalar_t<T>* = nullptr>
-inline size_t size(const T& /*x*/) {
-  return 1U;
+inline constexpr Eigen::Index size(const T& /*x*/) noexcept {
+  return static_cast<Eigen::Index>(1);
 }
 
 /** \ingroup type_trait
@@ -25,12 +25,12 @@ inline size_t size(const T& /*x*/) {
  * @tparam T type of m
  */
 template <typename T, require_container_t<T>* = nullptr>
-inline size_t size(const T& m) {
+inline constexpr Eigen::Index size(const T& m) noexcept {
   return m.size();
 }
 
 template <typename T, require_var_matrix_t<T>* = nullptr>
-inline size_t size(const T& m) {
+inline Eigen::Index size(const T& m) noexcept {
   return m.size();
 }
 
diff --git a/stan/math/prim/fun/size_mvt.hpp b/stan/math/prim/fun/size_mvt.hpp
index 6543412964f..1c225047de0 100644
--- a/stan/math/prim/fun/size_mvt.hpp
+++ b/stan/math/prim/fun/size_mvt.hpp
@@ -21,17 +21,18 @@ namespace math {
  * @throw std::invalid_argument since the type is a scalar.
  */
 template <typename ScalarT, require_stan_scalar_t<ScalarT>* = nullptr>
-size_t size_mvt(const ScalarT& /* unused */) {
+inline Eigen::Index size_mvt(const ScalarT& /* unused */) {
   throw std::invalid_argument("size_mvt passed to an unrecognized type.");
+  return 1;
 }
 
 template <typename MatrixT, require_matrix_t<MatrixT>* = nullptr>
-size_t size_mvt(const MatrixT& /* unused */) {
-  return 1U;
+inline Eigen::Index constexpr size_mvt(const MatrixT& /* unused */) {
+  return static_cast<Eigen::Index>(1);
 }
 
 template <typename MatrixT, require_matrix_t<MatrixT>* = nullptr>
-size_t size_mvt(const std::vector<MatrixT>& x) {
+inline Eigen::Index size_mvt(const std::vector<MatrixT>& x) {
   return x.size();
 }
 
diff --git a/stan/math/prim/fun/size_zero.hpp b/stan/math/prim/fun/size_zero.hpp
index 8af469eee9f..11d128cd751 100644
--- a/stan/math/prim/fun/size_zero.hpp
+++ b/stan/math/prim/fun/size_zero.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_PRIM_FUN_SIZE_ZERO_HPP
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/size.hpp>
 #include <utility>
 
 namespace stan {
@@ -15,8 +16,8 @@ namespace math {
  * @return 0 or 1
  */
 template <typename T>
-inline bool size_zero(const T& x) {
-  return !size(x);
+inline constexpr bool size_zero(const T& x) noexcept {
+  return !stan::math::size(x);
 }
 
 /**
@@ -28,7 +29,7 @@ inline bool size_zero(const T& x) {
  * @return 0 or 1
  */
 template <typename T, typename... Ts>
-inline bool size_zero(const T& x, const Ts&... xs) {
+inline constexpr bool size_zero(const T& x, const Ts&... xs) noexcept {
   return (size_zero(x) || size_zero(xs...));
 }
 }  // namespace math
diff --git a/stan/math/prim/fun/vector_seq_view.hpp b/stan/math/prim/fun/vector_seq_view.hpp
index 44903ca3db5..23d368ed367 100644
--- a/stan/math/prim/fun/vector_seq_view.hpp
+++ b/stan/math/prim/fun/vector_seq_view.hpp
@@ -36,7 +36,7 @@ template <typename T>
 class vector_seq_view<T, require_matrix_t<T>> {
  public:
   explicit vector_seq_view(const T& m) : m_(m) {}
-  static constexpr auto size() { return 1; }
+  static constexpr Eigen::Index size() { return 1; }
   inline const auto& operator[](size_t /* i */) const noexcept { return m_; }
 
   template <typename C = T, require_st_arithmetic<C>* = nullptr>
@@ -75,7 +75,7 @@ class vector_seq_view<
     T, require_std_vector_vt<internal::is_matrix_or_std_vector, T>> {
  public:
   explicit vector_seq_view(const T& v) noexcept : v_(v) {}
-  inline auto size() const noexcept { return v_.size(); }
+  inline Eigen::Index size() const noexcept { return v_.size(); }
 
   inline decltype(auto) operator[](size_t i) const { return v_[i]; }
 
diff --git a/stan/math/prim/functor/coupled_ode_system.hpp b/stan/math/prim/functor/coupled_ode_system.hpp
index 90df071de8b..67cf38f5f13 100644
--- a/stan/math/prim/functor/coupled_ode_system.hpp
+++ b/stan/math/prim/functor/coupled_ode_system.hpp
@@ -67,9 +67,9 @@ struct coupled_ode_system_impl<true, F, T_y0, Args...> {
 
     dz_dt.resize(y.size());
 
-    Eigen::VectorXd f_y_t
-        = apply([&](const Args&... args) { return f_(t, y, msgs_, args...); },
-                args_tuple_);
+    Eigen::VectorXd f_y_t = math::apply(
+        [&](const Args&... args) { return f_(t, y, msgs_, args...); },
+        args_tuple_);
 
     check_size_match("coupled_ode_system", "dy_dt", f_y_t.size(), "states",
                      y.size());
@@ -82,7 +82,7 @@ struct coupled_ode_system_impl<true, F, T_y0, Args...> {
    *
    * @return size of the coupled system.
    */
-  size_t size() const { return N_; }
+  inline Eigen::Index size() const noexcept { return N_; }
 
   /**
    * Returns the initial state of the coupled system. For arithmetic types
@@ -91,7 +91,7 @@ struct coupled_ode_system_impl<true, F, T_y0, Args...> {
    * @return the initial condition of the coupled system
    */
   std::vector<double> initial_state() const {
-    std::vector<double> initial(size(), 0.0);
+    std::vector<double> initial(this->size(), 0.0);
 
     for (size_t i = 0; i < N_; i++) {
       initial[i] = value_of(y0_(i));
diff --git a/stan/math/prim/functor/ode_ckrk.hpp b/stan/math/prim/functor/ode_ckrk.hpp
index 3f8fd0a0972..4f093093683 100644
--- a/stan/math/prim/functor/ode_ckrk.hpp
+++ b/stan/math/prim/functor/ode_ckrk.hpp
@@ -1,11 +1,11 @@
 #ifndef STAN_MATH_PRIM_FUNCTOR_ODE_CKRK_HPP
 #define STAN_MATH_PRIM_FUNCTOR_ODE_CKRK_HPP
 
+#include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/coupled_ode_system.hpp>
 #include <stan/math/prim/functor/ode_store_sensitivities.hpp>
-#include <stan/math/prim/fun/Eigen.hpp>
 #include <boost/numeric/odeint.hpp>
 #include <ostream>
 #include <tuple>
@@ -78,7 +78,7 @@ ode_ckrk_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   std::tuple<ref_type_t<Args>...> args_ref_tuple(args...);
 
-  apply(
+  math::apply(
       [&](const auto&... args_ref) {
         // Code from https://stackoverflow.com/a/17340003
         std::vector<int> unused_temp{
@@ -101,7 +101,7 @@ ode_ckrk_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   using return_t = return_type_t<T_y0, T_t0, T_ts, Args...>;
   // creates basic or coupled system by template specializations
-  auto&& coupled_system = apply(
+  auto&& coupled_system = math::apply(
       [&](const auto&... args_ref) {
         return coupled_ode_system<F, T_y0_t0, ref_type_t<Args>...>(f, y0, msgs,
                                                                    args_ref...);
@@ -127,7 +127,7 @@ ode_ckrk_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
       observer_initial_recorded = true;
       return;
     }
-    apply(
+    math::apply(
         [&](const auto&... args_ref) {
           y.emplace_back(ode_store_sensitivities(
               f, coupled_state, y0, t0, ts[time_index], msgs, args_ref...));
diff --git a/stan/math/prim/functor/ode_rk45.hpp b/stan/math/prim/functor/ode_rk45.hpp
index 8a90bd1e2c7..a03a4aebd8d 100644
--- a/stan/math/prim/functor/ode_rk45.hpp
+++ b/stan/math/prim/functor/ode_rk45.hpp
@@ -1,11 +1,11 @@
 #ifndef STAN_MATH_PRIM_FUNCTOR_ODE_RK45_HPP
 #define STAN_MATH_PRIM_FUNCTOR_ODE_RK45_HPP
 
+#include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/coupled_ode_system.hpp>
 #include <stan/math/prim/functor/ode_store_sensitivities.hpp>
-#include <stan/math/prim/fun/Eigen.hpp>
 #include <boost/numeric/odeint.hpp>
 #include <ostream>
 #include <tuple>
@@ -79,7 +79,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   std::tuple<ref_type_t<Args>...> args_ref_tuple(args...);
 
-  apply(
+  math::apply(
       [&](const auto&... args_ref) {
         // Code from https://stackoverflow.com/a/17340003
         std::vector<int> unused_temp{
@@ -102,7 +102,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   using return_t = return_type_t<T_y0, T_t0, T_ts, Args...>;
   // creates basic or coupled system by template specializations
-  auto&& coupled_system = apply(
+  auto&& coupled_system = math::apply(
       [&](const auto&... args_ref) {
         return coupled_ode_system<F, T_y0_t0, ref_type_t<Args>...>(f, y0, msgs,
                                                                    args_ref...);
@@ -128,7 +128,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
       observer_initial_recorded = true;
       return;
     }
-    apply(
+    math::apply(
         [&](const auto&... args_ref) {
           y.emplace_back(ode_store_sensitivities(
               f, coupled_state, y0, t0, ts[time_index], msgs, args_ref...));
diff --git a/stan/math/prim/functor/operands_and_partials.hpp b/stan/math/prim/functor/operands_and_partials.hpp
index 3d2f70b6fb4..c9817419071 100644
--- a/stan/math/prim/functor/operands_and_partials.hpp
+++ b/stan/math/prim/functor/operands_and_partials.hpp
@@ -66,7 +66,7 @@ struct ops_partials_edge<ViewElt, Op, require_st_arithmetic<Op>> {
    * Return the size of the operand for the edge. For doubles this is a compile
    * time expression returning zero.
    */
-  static constexpr int size() noexcept { return 0; }  // reverse mode
+  static constexpr Eigen::Index size() noexcept { return 0; }  // reverse mode
 
  private:
   template <typename, typename, typename, typename, typename, typename>
diff --git a/stan/math/prim/functor/reduce_sum.hpp b/stan/math/prim/functor/reduce_sum.hpp
index 868de060f32..c9b70a461da 100644
--- a/stan/math/prim/functor/reduce_sum.hpp
+++ b/stan/math/prim/functor/reduce_sum.hpp
@@ -3,7 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-
+#include <stan/math/prim/functor/apply.hpp>
 #include <tbb/task_arena.h>
 #include <tbb/parallel_reduce.h>
 #include <tbb/blocked_range.h>
@@ -80,7 +80,7 @@ struct reduce_sum_impl<ReduceFunction, require_arithmetic_t<ReturnType>,
         sub_slice.emplace_back(vmapped_[i]);
       }
 
-      sum_ += apply(
+      sum_ += math::apply(
           [&](auto&&... args) {
             return ReduceFunction()(sub_slice, r.begin(), r.end() - 1, &msgs_,
                                     args...);
diff --git a/stan/math/prim/meta/seq_view.hpp b/stan/math/prim/meta/seq_view.hpp
index 6f3decda7a6..ac469cf3e56 100644
--- a/stan/math/prim/meta/seq_view.hpp
+++ b/stan/math/prim/meta/seq_view.hpp
@@ -42,7 +42,7 @@ class seq_view {
  public:
   explicit seq_view(typename pass_type<S>::type x) : x_(x) {}
   inline typename pass_type<T>::type operator[](int n) const { return x_; }
-  int size() const { return 1; }
+  Eigen::Index size() const { return 1; }
 };
 
 template <typename T, typename S>
@@ -55,7 +55,7 @@ class seq_view<T, Eigen::Matrix<S, Eigen::Dynamic, 1> > {
       typename pass_type<Eigen::Matrix<S, Eigen::Dynamic, 1> >::type x)
       : x_(x) {}
   inline typename pass_type<T>::type operator[](int n) const { return x_(n); }
-  int size() const { return x_.size(); }
+  Eigen::Index size() const { return x_.size(); }
 };
 
 template <typename T, typename S>
@@ -68,7 +68,7 @@ class seq_view<T, Eigen::Matrix<S, 1, Eigen::Dynamic> > {
       typename pass_type<Eigen::Matrix<S, 1, Eigen::Dynamic> >::type x)
       : x_(x) {}
   inline typename pass_type<T>::type operator[](int n) const { return x_(n); }
-  int size() const { return x_.size(); }
+  Eigen::Index size() const { return x_.size(); }
 };
 
 // row-major order of returns to match std::vector
@@ -85,7 +85,7 @@ class seq_view<T, Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> > {
   inline typename pass_type<T>::type operator[](int n) const {
     return x_(n / x_.cols(), n % x_.cols());
   }
-  int size() const { return x_.size(); }
+  Eigen::Index size() const { return x_.size(); }
 };
 
 // question is how expensive the ctor is
@@ -101,7 +101,7 @@ class seq_view<T, std::vector<S> > {
   inline typename pass_type<T>::type operator[](int n) const {
     return seq_view<T, S>(x_[n / elt_size_])[n % elt_size_];
   }
-  int size() const { return x_.size() * elt_size_; }
+  Eigen::Index size() const { return x_.size() * elt_size_; }
 };
 
 // BELOW HERE JUST FOR EFFICIENCY
@@ -114,7 +114,7 @@ class seq_view<T, std::vector<T> > {
  public:
   explicit seq_view(typename pass_type<std::vector<T> >::type x) : x_(x) {}
   inline typename pass_type<T>::type operator[](int n) const { return x_[n]; }
-  int size() const { return x_.size(); }
+  Eigen::Index size() const { return x_.size(); }
 };
 
 // if vector of S with S assignable to T, also works
@@ -131,7 +131,7 @@ class seq_view<T, std::vector<std::vector<T> > > {
   inline typename pass_type<T>::type operator[](int n) const {
     return x_[n / cols_][n % cols_];
   }
-  int size() const { return x_.size() * cols_; }
+  Eigen::Index size() const { return x_.size() * cols_; }
 };
 
 template <>
@@ -142,7 +142,7 @@ class seq_view<double, std::vector<int> > {
  public:
   explicit seq_view(pass_type<std::vector<int> >::type x) : x_(x) {}
   inline pass_type<double>::type operator[](int n) const { return x_[n]; }
-  int size() const { return x_.size(); }
+  Eigen::Index size() const { return x_.size(); }
 };
 
 }  // namespace math
diff --git a/stan/math/prim/prob/bernoulli_lpmf.hpp b/stan/math/prim/prob/bernoulli_lpmf.hpp
index 29fc79d1bf8..0a0e0ed39aa 100644
--- a/stan/math/prim/prob/bernoulli_lpmf.hpp
+++ b/stan/math/prim/prob/bernoulli_lpmf.hpp
@@ -59,7 +59,7 @@ return_type_t<T_prob> bernoulli_lpmf(const T_n& n, const T_prob& theta) {
   scalar_seq_view<T_theta_ref> theta_vec(theta_ref);
   size_t N = max_size(n, theta);
 
-  if (size(theta) == 1) {
+  if (stan::math::size(theta) == 1) {
     size_t sum = 0;
     for (size_t n = 0; n < N; n++) {
       sum += n_vec.val(n);
diff --git a/stan/math/prim/prob/beta_proportion_lpdf.hpp b/stan/math/prim/prob/beta_proportion_lpdf.hpp
index c82ba7d40b5..523fdc1c138 100644
--- a/stan/math/prim/prob/beta_proportion_lpdf.hpp
+++ b/stan/math/prim/prob/beta_proportion_lpdf.hpp
@@ -89,7 +89,7 @@ return_type_t<T_y, T_loc, T_prec> beta_proportion_lpdf(const T_y& y,
   size_t N = max_size(y, mu, kappa);
   T_partials_return logp(0);
   if (include_summand<propto, T_prec>::value) {
-    logp += sum(lgamma(kappa_val)) * N / size(kappa);
+    logp += sum(lgamma(kappa_val)) * N / math::size(kappa);
   }
   if (include_summand<propto, T_loc, T_prec>::value) {
     logp -= sum(lgamma(mukappa) + lgamma(kappa_val - mukappa)) * N
diff --git a/stan/math/prim/prob/binomial_logit_lpmf.hpp b/stan/math/prim/prob/binomial_logit_lpmf.hpp
index d530ee15a4b..08e42e2873a 100644
--- a/stan/math/prim/prob/binomial_logit_lpmf.hpp
+++ b/stan/math/prim/prob/binomial_logit_lpmf.hpp
@@ -87,10 +87,11 @@ return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
       ops_partials.edge1_.partials_
           = n_val * inv_logit_neg_alpha - (N_val - n_val) * inv_logit_alpha;
     } else {
-      T_partials_return sum_n = sum(n_val) * maximum_size / size(n);
+      T_partials_return sum_n = sum(n_val) * maximum_size / math::size(n);
       ops_partials.edge1_.partials_[0] = forward_as<T_partials_return>(
           sum_n * inv_logit_neg_alpha
-          - (sum(N_val) * maximum_size / size(N) - sum_n) * inv_logit_alpha);
+          - (sum(N_val) * maximum_size / math::size(N) - sum_n)
+                * inv_logit_alpha);
     }
   }
 
diff --git a/stan/math/prim/prob/cauchy_lpdf.hpp b/stan/math/prim/prob/cauchy_lpdf.hpp
index 59dba34dcfc..6e16581e786 100644
--- a/stan/math/prim/prob/cauchy_lpdf.hpp
+++ b/stan/math/prim/prob/cauchy_lpdf.hpp
@@ -84,7 +84,7 @@ return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
     logp -= N * LOG_PI;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
diff --git a/stan/math/prim/prob/chi_square_cdf.hpp b/stan/math/prim/prob/chi_square_cdf.hpp
index ea439f86bf6..6f8c1a00a9e 100644
--- a/stan/math/prim/prob/chi_square_cdf.hpp
+++ b/stan/math/prim/prob/chi_square_cdf.hpp
@@ -69,9 +69,9 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lccdf.hpp b/stan/math/prim/prob/chi_square_lccdf.hpp
index 8347406a373..3330a8913e9 100644
--- a/stan/math/prim/prob/chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/chi_square_lccdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> chi_square_lccdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lcdf.hpp b/stan/math/prim/prob/chi_square_lcdf.hpp
index dc7c196a454..3ccbc065a3d 100644
--- a/stan/math/prim/prob/chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/chi_square_lcdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> chi_square_lcdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lpdf.hpp b/stan/math/prim/prob/chi_square_lpdf.hpp
index d2044d75986..885d9a3c988 100644
--- a/stan/math/prim/prob/chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/chi_square_lpdf.hpp
@@ -74,12 +74,12 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {
 
   T_partials_return logp(0);
   if (include_summand<propto, T_dof>::value) {
-    logp -= sum(nu_val * HALF_LOG_TWO + lgamma(half_nu)) * N / size(nu);
+    logp -= sum(nu_val * HALF_LOG_TWO + lgamma(half_nu)) * N / math::size(nu);
   }
   logp += sum((half_nu - 1.0) * log_y);
 
   if (include_summand<propto, T_y>::value) {
-    logp -= 0.5 * sum(y_val) * N / size(y);
+    logp -= 0.5 * sum(y_val) * N / math::size(y);
   }
 
   operands_and_partials<T_y_ref, T_nu_ref> ops_partials(y_ref, nu_ref);
diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp
index e15dfeceebe..3db919362d2 100644
--- a/stan/math/prim/prob/double_exponential_lpdf.hpp
+++ b/stan/math/prim/prob/double_exponential_lpdf.hpp
@@ -83,7 +83,7 @@ return_type_t<T_y, T_loc, T_scale> double_exponential_lpdf(
     logp -= N * LOG_TWO;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   logp -= sum(scaled_diff);
 
diff --git a/stan/math/prim/prob/exp_mod_normal_lccdf.hpp b/stan/math/prim/prob/exp_mod_normal_lccdf.hpp
index 115eb02ee81..6ec01d2f6dd 100644
--- a/stan/math/prim/prob/exp_mod_normal_lccdf.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_lccdf.hpp
@@ -66,8 +66,8 @@ return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_lccdf(
 
   scalar_seq_view<decltype(y_val)> y_vec(y_val);
   for (size_t n = 0, size_y = stan::math::size(y); n < size_y; n++) {
-    if (is_inf(y_vec[n])) {
-      return ops_partials.build(y_vec[n] > 0 ? negative_infinity() : 0);
+    if (is_inf(y_vec.val(n))) {
+      return ops_partials.build(y_vec.val(n) > 0 ? negative_infinity() : 0);
     }
   }
 
diff --git a/stan/math/prim/prob/exp_mod_normal_lcdf.hpp b/stan/math/prim/prob/exp_mod_normal_lcdf.hpp
index ee2f0bd3438..940cde7f6a8 100644
--- a/stan/math/prim/prob/exp_mod_normal_lcdf.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_lcdf.hpp
@@ -66,8 +66,8 @@ return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_lcdf(
 
   scalar_seq_view<decltype(y_val)> y_vec(y_val);
   for (size_t n = 0, size_y = stan::math::size(y); n < size_y; n++) {
-    if (is_inf(y_vec[n])) {
-      return ops_partials.build(y_vec[n] < 0 ? negative_infinity() : 0);
+    if (is_inf(y_vec.val(n))) {
+      return ops_partials.build(y_vec.val(n) < 0 ? negative_infinity() : 0);
     }
   }
 
diff --git a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
index b2db8f387c9..ba3eddbcbde 100644
--- a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
@@ -79,7 +79,7 @@ return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_lpdf(
     logp -= LOG_TWO * N;
   }
   if (include_summand<propto, T_inv_scale>::value) {
-    logp += sum(log(lambda_val)) * N / size(lambda);
+    logp += sum(log(lambda_val)) * N / math::size(lambda);
   }
   const auto& log_erfc_calc = log(erfc_calc);
   logp
diff --git a/stan/math/prim/prob/exponential_lccdf.hpp b/stan/math/prim/prob/exponential_lccdf.hpp
index d3e9d09deb4..c48c90ab9d3 100644
--- a/stan/math/prim/prob/exponential_lccdf.hpp
+++ b/stan/math/prim/prob/exponential_lccdf.hpp
@@ -47,7 +47,7 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     using beta_val_array = Eigen::Array<beta_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_y>::value && !is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = T_partials_array::Constant(
-          size(y), -forward_as<beta_val_scalar>(beta_val));
+          math::size(y), -forward_as<beta_val_scalar>(beta_val));
     } else if (is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = -forward_as<beta_val_array>(beta_val);
     } else {
@@ -59,7 +59,7 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     using y_val_array = Eigen::Array<y_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_inv_scale>::value && !is_vector<T_y>::value) {
       ops_partials.edge2_.partials_ = T_partials_array::Constant(
-          size(beta), -forward_as<y_val_scalar>(y_val));
+          math::size(beta), -forward_as<y_val_scalar>(y_val));
     } else if (is_vector<T_y>::value) {
       ops_partials.edge2_.partials_ = -forward_as<y_val_array>(y_val);
     } else {
diff --git a/stan/math/prim/prob/exponential_lpdf.hpp b/stan/math/prim/prob/exponential_lpdf.hpp
index 4328c1af2d6..18f3bbfd0da 100644
--- a/stan/math/prim/prob/exponential_lpdf.hpp
+++ b/stan/math/prim/prob/exponential_lpdf.hpp
@@ -76,7 +76,7 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
 
   T_partials_return logp(0.0);
   if (include_summand<propto, T_inv_scale>::value) {
-    logp = sum(log(beta_val)) * max_size(y, beta) / size(beta);
+    logp = sum(log(beta_val)) * max_size(y, beta) / math::size(beta);
   }
   if (include_summand<propto, T_y, T_inv_scale>::value) {
     logp -= sum(beta_val * y_val);
@@ -87,7 +87,7 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
     using beta_val_array = Eigen::Array<beta_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_y>::value && !is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = T_partials_array::Constant(
-          size(y), -forward_as<beta_val_scalar>(beta_val));
+          math::size(y), -forward_as<beta_val_scalar>(beta_val));
     } else if (is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = -forward_as<beta_val_array>(beta_val);
     } else {
diff --git a/stan/math/prim/prob/frechet_lpdf.hpp b/stan/math/prim/prob/frechet_lpdf.hpp
index fe0f6783bef..b305f0a8f20 100644
--- a/stan/math/prim/prob/frechet_lpdf.hpp
+++ b/stan/math/prim/prob/frechet_lpdf.hpp
@@ -72,7 +72,7 @@ return_type_t<T_y, T_shape, T_scale> frechet_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, sigma);
   T_partials_return logp = -sum(sigma_div_y_pow_alpha);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp -= sum((alpha_val + 1.0) * log_y) * N / max_size(y, alpha);
diff --git a/stan/math/prim/prob/gamma_cdf.hpp b/stan/math/prim/prob/gamma_cdf.hpp
index 2918ebb9f22..c3768682032 100644
--- a/stan/math/prim/prob/gamma_cdf.hpp
+++ b/stan/math/prim/prob/gamma_cdf.hpp
@@ -81,9 +81,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_cdf(const T_y& y,
   using std::pow;
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(stan::math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(stan::math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp
index 12fbd24a9f2..39bfd3e5ed7 100644
--- a/stan/math/prim/prob/gamma_lccdf.hpp
+++ b/stan/math/prim/prob/gamma_lccdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(stan::math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(stan::math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lcdf.hpp b/stan/math/prim/prob/gamma_lcdf.hpp
index b40c75bc7f3..3ffb09aa114 100644
--- a/stan/math/prim/prob/gamma_lcdf.hpp
+++ b/stan/math/prim/prob/gamma_lcdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(stan::math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(stan::math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lpdf.hpp b/stan/math/prim/prob/gamma_lpdf.hpp
index ba289902c9d..d7a033c70eb 100644
--- a/stan/math/prim/prob/gamma_lpdf.hpp
+++ b/stan/math/prim/prob/gamma_lpdf.hpp
@@ -84,7 +84,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lpdf(const T_y& y,
 
   scalar_seq_view<decltype(y_val)> y_vec(y_val);
   for (size_t n = 0; n < stan::math::size(y); n++) {
-    if (y_vec[n] < 0) {
+    if (y_vec.val(n) < 0) {
       return LOG_ZERO;
     }
   }
@@ -92,7 +92,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, beta);
   T_partials_return logp(0.0);
   if (include_summand<propto, T_shape>::value) {
-    logp = -sum(lgamma(alpha_val)) * N / size(alpha);
+    logp = -sum(lgamma(alpha_val)) * N / math::size(alpha);
   }
   const auto& log_y = to_ref_if<is_constant_all<T_shape>::value>(log(y_val));
   if (include_summand<propto, T_shape, T_inv_scale>::value) {
diff --git a/stan/math/prim/prob/gumbel_lpdf.hpp b/stan/math/prim/prob/gumbel_lpdf.hpp
index 998e4140f37..22814b8982b 100644
--- a/stan/math/prim/prob/gumbel_lpdf.hpp
+++ b/stan/math/prim/prob/gumbel_lpdf.hpp
@@ -78,7 +78,7 @@ return_type_t<T_y, T_loc, T_scale> gumbel_lpdf(const T_y& y, const T_loc& mu,
   size_t N = max_size(y, mu, beta);
   T_partials_return logp = -sum(y_minus_mu_over_beta + exp_y_m_mu_over_beta);
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(beta_val)) * N / size(beta);
+    logp -= sum(log(beta_val)) * N / math::size(beta);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
diff --git a/stan/math/prim/prob/inv_chi_square_cdf.hpp b/stan/math/prim/prob/inv_chi_square_cdf.hpp
index 1f1ff29c386..8865ff1c734 100644
--- a/stan/math/prim/prob/inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_cdf.hpp
@@ -69,9 +69,9 @@ return_type_t<T_y, T_dof> inv_chi_square_cdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/inv_chi_square_lccdf.hpp b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
index 394e7956f97..7e0a7031ab3 100644
--- a/stan/math/prim/prob/inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> inv_chi_square_lccdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/inv_chi_square_lcdf.hpp b/stan/math/prim/prob/inv_chi_square_lcdf.hpp
index abd879cb2c6..cb8894b507b 100644
--- a/stan/math/prim/prob/inv_chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lcdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> inv_chi_square_lcdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/inv_chi_square_lpdf.hpp b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
index aff4dece7f8..d47c4355017 100644
--- a/stan/math/prim/prob/inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
@@ -82,11 +82,12 @@ return_type_t<T_y, T_dof> inv_chi_square_lpdf(const T_y& y, const T_dof& nu) {
   size_t N = max_size(y, nu);
   T_partials_return logp = -sum((half_nu + 1.0) * log_y);
   if (include_summand<propto, T_dof>::value) {
-    logp -= (sum(nu_val) * HALF_LOG_TWO + sum(lgamma(half_nu))) * N / size(nu);
+    logp -= (sum(nu_val) * HALF_LOG_TWO + sum(lgamma(half_nu))) * N
+            / math::size(nu);
   }
   if (include_summand<propto, T_y>::value) {
     const auto& inv_y = to_ref_if<!is_constant_all<T_y>::value>(inv(y_val));
-    logp -= 0.5 * sum(inv_y) * N / size(y);
+    logp -= 0.5 * sum(inv_y) * N / math::size(y);
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_ = (0.5 * inv_y - half_nu - 1.0) * inv_y;
     }
diff --git a/stan/math/prim/prob/inv_gamma_cdf.hpp b/stan/math/prim/prob/inv_gamma_cdf.hpp
index 89271caa42d..12278edefe7 100644
--- a/stan/math/prim/prob/inv_gamma_cdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_cdf.hpp
@@ -80,9 +80,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_cdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(stan::math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(stan::math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lccdf.hpp b/stan/math/prim/prob/inv_gamma_lccdf.hpp
index 22212e5adff..52b1a367fac 100644
--- a/stan/math/prim/prob/inv_gamma_lccdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lccdf.hpp
@@ -66,9 +66,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lccdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(stan::math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(stan::math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lcdf.hpp b/stan/math/prim/prob/inv_gamma_lcdf.hpp
index 27e6800b444..e65ebe262d3 100644
--- a/stan/math/prim/prob/inv_gamma_lcdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lcdf.hpp
@@ -66,9 +66,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lcdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(stan::math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(stan::math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lpdf.hpp b/stan/math/prim/prob/inv_gamma_lpdf.hpp
index ceee125329f..0184a6ba572 100644
--- a/stan/math/prim/prob/inv_gamma_lpdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lpdf.hpp
@@ -83,7 +83,7 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lpdf(const T_y& y,
 
   size_t N = max_size(y, alpha, beta);
   if (include_summand<propto, T_shape>::value) {
-    logp -= sum(lgamma(alpha_val)) * N / size(alpha);
+    logp -= sum(lgamma(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_shape, T_scale>::value) {
     const auto& log_beta
diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp
index 5ed1b4982ba..f6359a50142 100644
--- a/stan/math/prim/prob/logistic_lpdf.hpp
+++ b/stan/math/prim/prob/logistic_lpdf.hpp
@@ -66,7 +66,7 @@ return_type_t<T_y, T_loc, T_scale> logistic_lpdf(const T_y& y, const T_loc& mu,
   T_partials_return logp = -sum(y_minus_mu_div_sigma)
                            - 2.0 * sum(log1p(exp(-y_minus_mu_div_sigma)));
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_scale>::value) {
diff --git a/stan/math/prim/prob/lognormal_lpdf.hpp b/stan/math/prim/prob/lognormal_lpdf.hpp
index 4afdbae67b5..422b19d65d6 100644
--- a/stan/math/prim/prob/lognormal_lpdf.hpp
+++ b/stan/math/prim/prob/lognormal_lpdf.hpp
@@ -74,10 +74,10 @@ return_type_t<T_y, T_loc, T_scale> lognormal_lpdf(const T_y& y, const T_loc& mu,
   T_partials_return logp
       = N * NEG_LOG_SQRT_TWO_PI - 0.5 * sum(square(logy_m_mu) * inv_sigma_sq);
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y>::value) {
-    logp -= sum(log_y) * N / size(y);
+    logp -= sum(log_y) * N / math::size(y);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
diff --git a/stan/math/prim/prob/neg_binomial_2_lccdf.hpp b/stan/math/prim/prob/neg_binomial_2_lccdf.hpp
index d7575ba64b8..f98276e1851 100644
--- a/stan/math/prim/prob/neg_binomial_2_lccdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_lccdf.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/divide.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
@@ -31,18 +32,9 @@ return_type_t<T_location, T_precision> neg_binomial_2_lccdf(
     return 0;
   }
 
-  scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
-  scalar_seq_view<T_phi_ref> phi_vec(phi_ref);
-  size_t size_beta = max_size(mu, phi);
-
-  VectorBuilder<true, return_type_t<T_location, T_precision>, T_location,
-                T_precision>
-      beta_vec(size_beta);
-  for (size_t i = 0; i < size_beta; ++i) {
-    beta_vec[i] = phi_vec[i] / mu_vec[i];
-  }
-
-  return neg_binomial_lccdf(n, phi_ref, beta_vec.data());
+  auto beta_vec = stan::math::divide(as_column_vector_or_scalar(phi_ref),
+                                     as_column_vector_or_scalar(mu_ref));
+  return neg_binomial_lccdf(n, phi_ref, beta_vec);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp
index fa0a1302994..9cf7f28f4ae 100644
--- a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp
@@ -3,7 +3,9 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/add.hpp>
 #include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/divide.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
@@ -38,30 +40,19 @@ return_type_t<T_location, T_precision> neg_binomial_2_lcdf(
   }
 
   scalar_seq_view<T_n_ref> n_vec(n_ref);
-  scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
-  scalar_seq_view<T_phi_ref> phi_vec(phi_ref);
   size_t size_n = stan::math::size(n);
-  size_t size_phi_mu = max_size(mu, phi);
 
   for (size_t i = 0; i < size_n; i++) {
     if (n_vec[i] < 0) {
       return LOG_ZERO;
     }
   }
+  auto&& phi_vec = as_column_vector_or_scalar(phi_ref);
+  auto phi_mu = stan::math::divide(
+      phi_vec, stan::math::add(phi_vec, as_column_vector_or_scalar(mu_ref)));
 
-  VectorBuilder<true, return_type_t<T_location, T_precision>, T_location,
-                T_precision>
-      phi_mu(size_phi_mu);
-  for (size_t i = 0; i < size_phi_mu; i++) {
-    phi_mu[i] = phi_vec[i] / (phi_vec[i] + mu_vec[i]);
-  }
-
-  VectorBuilder<true, return_type_t<T_n>, T_n> np1(size_n);
-  for (size_t i = 0; i < size_n; i++) {
-    np1[i] = n_vec[i] + 1.0;
-  }
-
-  return beta_cdf_log(phi_mu.data(), phi_ref, np1.data());
+  auto np1 = as_array_or_scalar(n_ref) + 1.0;
+  return beta_cdf_log(phi_mu, phi_ref, np1);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/prob/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_log_glm_lpmf.hpp
index 42694d69f6c..73a0cf7d05e 100644
--- a/stan/math/prim/prob/neg_binomial_2_log_glm_lpmf.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_log_glm_lpmf.hpp
@@ -168,7 +168,8 @@ return_type_t<T_x, T_alpha, T_beta, T_precision> neg_binomial_2_log_glm_lpmf(
     if (is_vector<T_precision>::value) {
       scalar_seq_view<decltype(phi_val_vec)> phi_vec(phi_val_vec);
       for (size_t n = 0; n < N_instances; ++n) {
-        logp += multiply_log(phi_vec[n], phi_vec[n]) - lgamma(phi_vec[n]);
+        logp += multiply_log(phi_vec.val(n), phi_vec.val(n))
+                - lgamma(phi_vec.val(n));
       }
     } else {
       logp += N_instances
diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp
index e5a1be8265a..5d688964549 100644
--- a/stan/math/prim/prob/normal_lpdf.hpp
+++ b/stan/math/prim/prob/normal_lpdf.hpp
@@ -85,7 +85,7 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
     logp += NEG_LOG_SQRT_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_scale, T_loc>::value) {
diff --git a/stan/math/prim/prob/normal_sufficient_lpdf.hpp b/stan/math/prim/prob/normal_sufficient_lpdf.hpp
index 475a13aefd1..abf5924b75d 100644
--- a/stan/math/prim/prob/normal_sufficient_lpdf.hpp
+++ b/stan/math/prim/prob/normal_sufficient_lpdf.hpp
@@ -107,7 +107,7 @@ return_type_t<T_y, T_s, T_loc, T_scale> normal_sufficient_lpdf(
   size_t N = max_size(y_bar, s_squared, n_obs, mu, sigma);
   T_partials_return logp = -sum(cons_expr / (2 * sigma_squared));
   if (include_summand<propto>::value) {
-    logp += NEG_LOG_SQRT_TWO_PI * sum(n_obs_val) * N / size(n_obs);
+    logp += NEG_LOG_SQRT_TWO_PI * sum(n_obs_val) * N / math::size(n_obs);
   }
   if (include_summand<propto, T_scale>::value) {
     logp -= sum(n_obs_val * log(sigma_val)) * N / max_size(n_obs, sigma);
@@ -141,7 +141,7 @@ return_type_t<T_y, T_s, T_loc, T_scale> normal_sufficient_lpdf(
       } else {
         forward_as<internal::broadcast_array<T_partials_return>>(
             ops_partials.edge2_.partials_)
-            = -0.5 / sigma_squared * N / size(sigma);
+            = -0.5 / sigma_squared * N / math::size(sigma);
       }
     }
   }
diff --git a/stan/math/prim/prob/ordered_logistic_lpmf.hpp b/stan/math/prim/prob/ordered_logistic_lpmf.hpp
index ccc92be4bd5..858c6c720ce 100644
--- a/stan/math/prim/prob/ordered_logistic_lpmf.hpp
+++ b/stan/math/prim/prob/ordered_logistic_lpmf.hpp
@@ -86,7 +86,7 @@ return_type_t<T_loc, T_cut> ordered_logistic_lpmf(const T_y& y,
 
   T_cut_ref c_ref = c;
   vector_seq_view<T_cut_ref> c_vec(c_ref);
-  int N = size(lambda);
+  int N = math::size(lambda);
   int C_l = size_mvt(c);
 
   check_consistent_sizes(function, "Integers", y, "Locations", lambda);
diff --git a/stan/math/prim/prob/ordered_probit_lpmf.hpp b/stan/math/prim/prob/ordered_probit_lpmf.hpp
index b7a89401737..18bbfdd0174 100644
--- a/stan/math/prim/prob/ordered_probit_lpmf.hpp
+++ b/stan/math/prim/prob/ordered_probit_lpmf.hpp
@@ -54,7 +54,7 @@ return_type_t<T_loc, T_cut> ordered_probit_lpmf(const T_y& y,
   static const char* function = "ordered_probit";
 
   check_nonzero_size(function, "Cut-points", c);
-  int N = size(lambda);
+  int N = math::size(lambda);
   int C_l = size_mvt(c);
   vector_seq_view<T_cut> c_vec(c);
   int K = c_vec[0].size() + 1;
diff --git a/stan/math/prim/prob/pareto_lpdf.hpp b/stan/math/prim/prob/pareto_lpdf.hpp
index 65f7d0e646d..e003450d7d7 100644
--- a/stan/math/prim/prob/pareto_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_lpdf.hpp
@@ -63,7 +63,7 @@ return_type_t<T_y, T_scale, T_shape> pareto_lpdf(const T_y& y,
   size_t N = max_size(y, y_min, alpha);
   T_partials_return logp(0);
   if (include_summand<propto, T_shape>::value) {
-    logp = sum(log(alpha_val)) * N / size(alpha);
+    logp = sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp -= sum(alpha_val * log_y + log_y) * N / max_size(alpha, y);
diff --git a/stan/math/prim/prob/pareto_type_2_lpdf.hpp b/stan/math/prim/prob/pareto_type_2_lpdf.hpp
index 3946f0a29d8..37500c11fd8 100644
--- a/stan/math/prim/prob/pareto_type_2_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_type_2_lpdf.hpp
@@ -66,10 +66,10 @@ return_type_t<T_y, T_loc, T_scale, T_shape> pareto_type_2_lpdf(
   size_t N = max_size(y, mu, lambda, alpha);
   T_partials_return logp(0.0);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(lambda_val)) * N / size(lambda);
+    logp -= sum(log(lambda_val)) * N / math::size(lambda);
   }
   if (include_summand<propto, T_y, T_loc, T_scale, T_shape>::value) {
     logp -= sum((alpha_val + 1.0) * log1p_scaled_diff);
diff --git a/stan/math/prim/prob/poisson_binomial_cdf.hpp b/stan/math/prim/prob/poisson_binomial_cdf.hpp
index 47ef2eca983..8b6cd988eb8 100644
--- a/stan/math/prim/prob/poisson_binomial_cdf.hpp
+++ b/stan/math/prim/prob/poisson_binomial_cdf.hpp
@@ -37,13 +37,13 @@ return_type_t<T_theta> poisson_binomial_cdf(const T_y& y,
                                             const T_theta& theta) {
   static const char* function = "poisson_binomial_cdf";
 
-  size_t size_theta = size_mvt(theta);
+  Eigen::Index size_theta = size_mvt(theta);
   if (size_theta > 1) {
     check_consistent_sizes(function, "Successes variables", y,
                            "Probability parameters", theta);
   }
 
-  size_t max_sz = std::max(stan::math::size(y), size_theta);
+  Eigen::Index max_sz = std::max(stan::math::size(y), size_theta);
   scalar_seq_view<T_y> y_vec(y);
   vector_seq_view<T_theta> theta_vec(theta);
 
diff --git a/stan/math/prim/prob/poisson_binomial_lccdf.hpp b/stan/math/prim/prob/poisson_binomial_lccdf.hpp
index f700c4795ef..b982cfb1ae3 100644
--- a/stan/math/prim/prob/poisson_binomial_lccdf.hpp
+++ b/stan/math/prim/prob/poisson_binomial_lccdf.hpp
@@ -38,17 +38,17 @@ return_type_t<T_theta> poisson_binomial_lccdf(const T_y& y,
                                               const T_theta& theta) {
   static const char* function = "poisson_binomial_lccdf";
 
-  size_t size_theta = size_mvt(theta);
+  const Eigen::Index size_theta = size_mvt(theta);
   if (size_theta > 1) {
     check_consistent_sizes(function, "Successes variables", y,
                            "Probability parameters", theta);
   }
 
-  size_t max_sz = std::max(stan::math::size(y), size_mvt(theta));
+  const Eigen::Index max_sz = std::max(stan::math::size(y), size_mvt(theta));
   scalar_seq_view<T_y> y_vec(y);
   vector_seq_view<T_theta> theta_vec(theta);
 
-  for (size_t i = 0; i < max_sz; ++i) {
+  for (Eigen::Index i = 0; i < max_sz; ++i) {
     check_bounded(function, "Successes variable", y_vec[i], 0,
                   theta_vec[i].size());
     check_finite(function, "Probability parameters", theta_vec.val(i));
@@ -57,7 +57,7 @@ return_type_t<T_theta> poisson_binomial_lccdf(const T_y& y,
   }
 
   return_type_t<T_theta> lccdf = 0.0;
-  for (size_t i = 0; i < max_sz; ++i) {
+  for (Eigen::Index i = 0; i < max_sz; ++i) {
     if (stan::math::size(theta_vec[i]) == 1) {
       if (y_vec[i] == 0) {
         lccdf += log(theta_vec[i][0]);
diff --git a/stan/math/prim/prob/poisson_binomial_lcdf.hpp b/stan/math/prim/prob/poisson_binomial_lcdf.hpp
index d3006865707..c7d3148cff4 100644
--- a/stan/math/prim/prob/poisson_binomial_lcdf.hpp
+++ b/stan/math/prim/prob/poisson_binomial_lcdf.hpp
@@ -37,13 +37,13 @@ return_type_t<T_theta> poisson_binomial_lcdf(const T_y& y,
                                              const T_theta& theta) {
   static const char* function = "poisson_binomial_lcdf";
 
-  size_t size_theta = size_mvt(theta);
+  const Eigen::Index size_theta = size_mvt(theta);
   if (size_theta > 1) {
     check_consistent_sizes(function, "Successes variables", y,
                            "Probability parameters", theta);
   }
 
-  size_t max_sz = std::max(stan::math::size(y), size_theta);
+  const Eigen::Index max_sz = std::max(stan::math::size(y), size_theta);
   scalar_seq_view<T_y> y_vec(y);
   vector_seq_view<T_theta> theta_vec(theta);
 
diff --git a/stan/math/prim/prob/poisson_binomial_lpmf.hpp b/stan/math/prim/prob/poisson_binomial_lpmf.hpp
index bca1f7680de..d6cba790d73 100644
--- a/stan/math/prim/prob/poisson_binomial_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_binomial_lpmf.hpp
@@ -29,13 +29,13 @@ return_type_t<T_theta> poisson_binomial_lpmf(const T_y& y,
                                              const T_theta& theta) {
   static const char* function = "poisson_binomial_lpmf";
 
-  size_t size_theta = size_mvt(theta);
+  const Eigen::Index size_theta = size_mvt(theta);
   if (size_theta > 1) {
     check_consistent_sizes(function, "Successes variables", y,
                            "Probability parameters", theta);
   }
 
-  size_t max_sz = std::max(stan::math::size(y), size_theta);
+  const Eigen::Index max_sz = std::max(stan::math::size(y), size_theta);
   scalar_seq_view<T_y> y_vec(y);
   vector_seq_view<T_theta> theta_vec(theta);
 
diff --git a/stan/math/prim/prob/poisson_log_lpmf.hpp b/stan/math/prim/prob/poisson_log_lpmf.hpp
index 44c884d8839..85be278124c 100644
--- a/stan/math/prim/prob/poisson_log_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_log_lpmf.hpp
@@ -73,10 +73,10 @@ return_type_t<T_log_rate> poisson_log_lpmf(const T_n& n,
 
   T_partials_return logp = sum(n_val * alpha_val);
   if (include_summand<propto, T_log_rate>::value) {
-    logp -= sum(exp_alpha) * N / size(alpha);
+    logp -= sum(exp_alpha) * N / math::size(alpha);
   }
   if (include_summand<propto>::value) {
-    logp -= sum(lgamma(n_val + 1.0)) * N / size(n);
+    logp -= sum(lgamma(n_val + 1.0)) * N / math::size(n);
   }
 
   if (!is_constant_all<T_log_rate>::value) {
diff --git a/stan/math/prim/prob/poisson_lpmf.hpp b/stan/math/prim/prob/poisson_lpmf.hpp
index 85732a158ed..01ab3832a95 100644
--- a/stan/math/prim/prob/poisson_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_lpmf.hpp
@@ -68,10 +68,10 @@ return_type_t<T_rate> poisson_lpmf(const T_n& n, const T_rate& lambda) {
 
   T_partials_return logp = stan::math::sum(multiply_log(n_val, lambda_val));
   if (include_summand<propto, T_rate>::value) {
-    logp -= sum(lambda_val) * N / size(lambda);
+    logp -= sum(lambda_val) * N / math::size(lambda);
   }
   if (include_summand<propto>::value) {
-    logp -= sum(lgamma(n_val + 1.0)) * N / size(n);
+    logp -= sum(lgamma(n_val + 1.0)) * N / math::size(n);
   }
 
   if (!is_constant_all<T_rate>::value) {
diff --git a/stan/math/prim/prob/rayleigh_lpdf.hpp b/stan/math/prim/prob/rayleigh_lpdf.hpp
index aebd139c317..1e447b33956 100644
--- a/stan/math/prim/prob/rayleigh_lpdf.hpp
+++ b/stan/math/prim/prob/rayleigh_lpdf.hpp
@@ -56,10 +56,10 @@ return_type_t<T_y, T_scale> rayleigh_lpdf(const T_y& y, const T_scale& sigma) {
   size_t N = max_size(y, sigma);
   T_partials_return logp = -0.5 * sum(square(y_over_sigma));
   if (include_summand<propto, T_scale>::value) {
-    logp -= 2.0 * sum(log(sigma_val)) * N / size(sigma);
+    logp -= 2.0 * sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y>::value) {
-    logp += sum(log(y_val)) * N / size(y);
+    logp += sum(log(y_val)) * N / math::size(y);
   }
 
   if (!is_constant_all<T_y, T_scale>::value) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
index f441852d649..acae28cafbb 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
@@ -78,9 +78,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_cdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
index 83b31549806..9d6bafc3ca5 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lccdf(
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
index c25413b2831..ac56d58af1c 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lcdf(
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
index 5187e4b03f5..db4f84bdb71 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
@@ -86,7 +86,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
                 T_partials_return, T_dof>
-      half_nu(size(nu));
+      half_nu(stan::math::size(nu));
   for (size_t i = 0; i < stan::math::size(nu); i++) {
     if (include_summand<propto, T_dof, T_y, T_scale>::value) {
       half_nu[i] = 0.5 * nu_vec.val(i);
@@ -95,7 +95,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y>::value, T_partials_return,
                 T_y>
-      log_y(size(y));
+      log_y(stan::math::size(y));
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (include_summand<propto, T_dof, T_y>::value) {
       log_y[i] = log(y_vec.val(i));
@@ -104,7 +104,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
                 T_partials_return, T_y>
-      inv_y(size(y));
+      inv_y(stan::math::size(y));
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (include_summand<propto, T_dof, T_y, T_scale>::value) {
       inv_y[i] = 1.0 / y_vec.val(i);
@@ -113,7 +113,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_scale>::value,
                 T_partials_return, T_scale>
-      log_s(size(s));
+      log_s(stan::math::size(s));
   for (size_t i = 0; i < stan::math::size(s); i++) {
     if (include_summand<propto, T_dof, T_scale>::value) {
       log_s[i] = log(s_vec.val(i));
@@ -121,11 +121,11 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
   }
 
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      log_half_nu(size(nu));
+      log_half_nu(stan::math::size(nu));
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      lgamma_half_nu(size(nu));
+      lgamma_half_nu(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_half_nu_over_two(size(nu));
+      digamma_half_nu_over_two(stan::math::size(nu));
   for (size_t i = 0; i < stan::math::size(nu); i++) {
     if (include_summand<propto, T_dof>::value) {
       lgamma_half_nu[i] = lgamma(half_nu[i]);
diff --git a/stan/math/prim/prob/skew_double_exponential_lpdf.hpp b/stan/math/prim/prob/skew_double_exponential_lpdf.hpp
index c14ba95f0c1..d3b6d5cea3e 100644
--- a/stan/math/prim/prob/skew_double_exponential_lpdf.hpp
+++ b/stan/math/prim/prob/skew_double_exponential_lpdf.hpp
@@ -97,10 +97,10 @@ return_type_t<T_y, T_loc, T_scale, T_skewness> skew_double_exponential_lpdf(
     logp += N * LOG_TWO;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_skewness>::value) {
-    logp += sum(log(tau_val) + log1m(tau_val)) * N / size(tau);
+    logp += sum(log(tau_val) + log1m(tau_val)) * N / math::size(tau);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/skew_normal_lpdf.hpp b/stan/math/prim/prob/skew_normal_lpdf.hpp
index c8481abf6c2..d811f507ae5 100644
--- a/stan/math/prim/prob/skew_normal_lpdf.hpp
+++ b/stan/math/prim/prob/skew_normal_lpdf.hpp
@@ -76,7 +76,7 @@ return_type_t<T_y, T_loc, T_scale, T_shape> skew_normal_lpdf(
     logp -= HALF_LOG_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y, T_loc, T_scale>::value) {
     logp -= sum(square(y_minus_mu_over_sigma)) * 0.5 * N
diff --git a/stan/math/prim/prob/student_t_cdf.hpp b/stan/math/prim/prob/student_t_cdf.hpp
index 38232b2ea0a..c781a818b5e 100644
--- a/stan/math/prim/prob/student_t_cdf.hpp
+++ b/stan/math/prim/prob/student_t_cdf.hpp
@@ -65,11 +65,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_cdf(const T_y& y,
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lccdf.hpp b/stan/math/prim/prob/student_t_lccdf.hpp
index 5c9e39865d2..72794e414a2 100644
--- a/stan/math/prim/prob/student_t_lccdf.hpp
+++ b/stan/math/prim/prob/student_t_lccdf.hpp
@@ -65,11 +65,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lccdf(
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lcdf.hpp b/stan/math/prim/prob/student_t_lcdf.hpp
index 401e938cc5b..ee821f77ce4 100644
--- a/stan/math/prim/prob/student_t_lcdf.hpp
+++ b/stan/math/prim/prob/student_t_lcdf.hpp
@@ -67,11 +67,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lcdf(const T_y& y,
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(stan::math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(stan::math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lpdf.hpp b/stan/math/prim/prob/student_t_lpdf.hpp
index 7df5ea97a22..b6ee5f3f176 100644
--- a/stan/math/prim/prob/student_t_lpdf.hpp
+++ b/stan/math/prim/prob/student_t_lpdf.hpp
@@ -111,10 +111,10 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lpdf(const T_y& y,
   if (include_summand<propto, T_dof>::value) {
     logp += (sum(lgamma(half_nu + 0.5)) - sum(lgamma(half_nu))
              - 0.5 * sum(log(nu_val)))
-            * N / size(nu);
+            * N / math::size(nu);
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/uniform_lcdf.hpp b/stan/math/prim/prob/uniform_lcdf.hpp
index 174a847b0db..83e8b22e0d8 100644
--- a/stan/math/prim/prob/uniform_lcdf.hpp
+++ b/stan/math/prim/prob/uniform_lcdf.hpp
@@ -68,7 +68,7 @@ return_type_t<T_y, T_low, T_high> uniform_lcdf(const T_y& y, const T_low& alpha,
   if (!is_constant_all<T_y>::value) {
     if (!is_vector<T_y>::value && is_vector<T_high>::value
         && !is_vector<T_low>::value) {
-      ops_partials.edge1_.partials_ = size(beta) * inv(y_minus_alpha);
+      ops_partials.edge1_.partials_ = math::size(beta) * inv(y_minus_alpha);
     } else {
       ops_partials.edge1_.partials_ = inv(y_minus_alpha);
     }
@@ -80,7 +80,7 @@ return_type_t<T_y, T_low, T_high> uniform_lcdf(const T_y& y, const T_low& alpha,
   if (!is_constant_all<T_high>::value) {
     if (is_vector<T_y>::value && !is_vector<T_low>::value
         && !is_vector<T_high>::value) {
-      ops_partials.edge3_.partials_ = inv(-b_minus_a) * size(y);
+      ops_partials.edge3_.partials_ = inv(-b_minus_a) * math::size(y);
     } else {
       ops_partials.edge3_.partials_ = inv(-b_minus_a);
     }
diff --git a/stan/math/prim/prob/uniform_lpdf.hpp b/stan/math/prim/prob/uniform_lpdf.hpp
index cbb2d50261a..e19290bad1c 100644
--- a/stan/math/prim/prob/uniform_lpdf.hpp
+++ b/stan/math/prim/prob/uniform_lpdf.hpp
@@ -95,7 +95,7 @@ return_type_t<T_y, T_low, T_high> uniform_lpdf(const T_y& y, const T_low& alpha,
     if (!is_constant_all<T_high>::value) {
       if (is_vector<T_y>::value && !is_vector<T_low>::value
           && !is_vector<T_high>::value) {
-        ops_partials.edge3_.partials_ = -inv_beta_minus_alpha * size(y);
+        ops_partials.edge3_.partials_ = -inv_beta_minus_alpha * math::size(y);
       } else {
         ops_partials.edge3_.partials_ = -inv_beta_minus_alpha;
       }
@@ -103,7 +103,7 @@ return_type_t<T_y, T_low, T_high> uniform_lpdf(const T_y& y, const T_low& alpha,
     if (!is_constant_all<T_low>::value) {
       if (is_vector<T_y>::value && !is_vector<T_low>::value
           && !is_vector<T_high>::value) {
-        ops_partials.edge2_.partials_ = inv_beta_minus_alpha * size(y);
+        ops_partials.edge2_.partials_ = inv_beta_minus_alpha * math::size(y);
       } else {
         ops_partials.edge2_.partials_ = std::move(inv_beta_minus_alpha);
       }
diff --git a/stan/math/prim/prob/von_mises_lpdf.hpp b/stan/math/prim/prob/von_mises_lpdf.hpp
index dec4681dc87..dcaef3036b1 100644
--- a/stan/math/prim/prob/von_mises_lpdf.hpp
+++ b/stan/math/prim/prob/von_mises_lpdf.hpp
@@ -63,7 +63,8 @@ return_type_t<T_y, T_loc, T_scale> von_mises_lpdf(T_y const& y, T_loc const& mu,
     logp -= LOG_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log_modified_bessel_first_kind(0, kappa_val)) * N / size(kappa);
+    logp -= sum(log_modified_bessel_first_kind(0, kappa_val)) * N
+            / math::size(kappa);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/weibull_lpdf.hpp b/stan/math/prim/prob/weibull_lpdf.hpp
index 5f69eff760d..b0eeeb43ea7 100644
--- a/stan/math/prim/prob/weibull_lpdf.hpp
+++ b/stan/math/prim/prob/weibull_lpdf.hpp
@@ -89,7 +89,7 @@ return_type_t<T_y, T_shape, T_scale> weibull_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, sigma);
   T_partials_return logp = -sum(y_div_sigma_pow_alpha);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp += sum((alpha_val - 1.0) * log_y) * N / max_size(alpha, y);
diff --git a/stan/math/rev/core/matrix_vari.hpp b/stan/math/rev/core/matrix_vari.hpp
index 3ed15c86b3f..1adad4ae2cd 100644
--- a/stan/math/rev/core/matrix_vari.hpp
+++ b/stan/math/rev/core/matrix_vari.hpp
@@ -25,7 +25,7 @@ class op_matrix_vari : public vari {
         = vs.vi();
   }
   vari* operator[](size_t n) const { return vis_[n]; }
-  size_t size() { return size_; }
+  inline Eigen::Index size() { return size_; }
 };
 
 }  // namespace math
diff --git a/stan/math/rev/core/profiling.hpp b/stan/math/rev/core/profiling.hpp
index ea51aa79499..7870f6eaf5f 100644
--- a/stan/math/rev/core/profiling.hpp
+++ b/stan/math/rev/core/profiling.hpp
@@ -19,22 +19,22 @@ namespace math {
  */
 class profile_info {
  private:
-  bool active_;
-
-  double fwd_pass_time_;
-  double rev_pass_time_;
-  size_t n_fwd_AD_passes_;
-  size_t n_fwd_no_AD_passes_;
-  size_t n_rev_passes_;
-  size_t chain_stack_size_sum_;
-  size_t nochain_stack_size_sum_;
+  bool active_{false};
+
+  double fwd_pass_time_{0.0};
+  double rev_pass_time_{0.0};
+  size_t n_fwd_AD_passes_{0};
+  size_t n_fwd_no_AD_passes_{0};
+  size_t n_rev_passes_{0};
+  size_t chain_stack_size_sum_{0};
+  size_t nochain_stack_size_sum_{0};
   std::chrono::time_point<std::chrono::steady_clock> fwd_pass_tp_;
   std::chrono::time_point<std::chrono::steady_clock> rev_pass_tp_;
-  size_t start_chain_stack_size_;
-  size_t start_nochain_stack_size_;
+  size_t start_chain_stack_size_{0};
+  size_t start_nochain_stack_size_{0};
 
  public:
-  profile_info()
+  profile_info() noexcept
       : active_(false),
         fwd_pass_time_(0.0),
         rev_pass_time_(0.0),
@@ -79,9 +79,11 @@ class profile_info {
     active_ = false;
   }
 
-  void rev_pass_start() { rev_pass_tp_ = std::chrono::steady_clock::now(); }
+  void rev_pass_start() noexcept {
+    rev_pass_tp_ = std::chrono::steady_clock::now();
+  }
 
-  void rev_pass_stop() {
+  void rev_pass_stop() noexcept {
     rev_pass_time_ += std::chrono::duration<double>(
                           std::chrono::steady_clock::now() - rev_pass_tp_)
                           .count();
@@ -135,8 +137,9 @@ class profile {
   profile_info* profile_;
 
  public:
-  profile(std::string name, profile_map& profiles)
-      : key_({name, std::this_thread::get_id()}) {
+  template <typename String>
+  profile(String&& name, profile_map& profiles)
+      : key_({std::forward<String>(name), std::this_thread::get_id()}) {
     profile_map::iterator p = profiles.find(key_);
     if (p == profiles.end()) {
       profiles[key_] = profile_info();
diff --git a/stan/math/rev/core/vector_vari.hpp b/stan/math/rev/core/vector_vari.hpp
index 57f7ebed6c6..2f2d2e5ce04 100644
--- a/stan/math/rev/core/vector_vari.hpp
+++ b/stan/math/rev/core/vector_vari.hpp
@@ -25,7 +25,7 @@ class op_vector_vari : public vari {
     }
   }
   vari* operator[](size_t n) const { return vis_[n]; }
-  size_t size() { return size_; }
+  inline Eigen::Index size() { return size_; }
 };
 
 }  // namespace math
diff --git a/stan/math/rev/functor/algebra_solver_newton.hpp b/stan/math/rev/functor/algebra_solver_newton.hpp
index d77921cd90c..04b426b212d 100644
--- a/stan/math/rev/functor/algebra_solver_newton.hpp
+++ b/stan/math/rev/functor/algebra_solver_newton.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/eval.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/algebra_solver_adapter.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <unsupported/Eigen/NonLinearOptimization>
 #include <iostream>
 #include <string>
@@ -142,7 +143,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
     const int64_t max_num_steps, const T_Args&... args) {
   const auto& x_ref = to_ref(value_of(x));
   auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...));
-  auto args_vals_tuple = apply(
+  auto args_vals_tuple = math::apply(
       [&](const auto&... args) {
         return std::make_tuple(to_ref(value_of(args))...);
       },
@@ -157,7 +158,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
   check_positive("algebra_solver_newton", "max_num_steps", max_num_steps);
 
   // Solve the system
-  Eigen::VectorXd theta_dbl = apply(
+  Eigen::VectorXd theta_dbl = math::apply(
       [&](const auto&... vals) {
         return kinsol_solve(f, x_ref, scaling_step_size, function_tolerance,
                             max_num_steps, 1, 10, KIN_LINESEARCH, msgs,
@@ -166,8 +167,8 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
       args_vals_tuple);
 
   auto f_wrt_x = [&](const auto& x) {
-    return apply([&](const auto&... args) { return f(x, msgs, args...); },
-                 args_vals_tuple);
+    return math::apply([&](const auto&... args) { return f(x, msgs, args...); },
+                       args_vals_tuple);
   };
 
   Eigen::MatrixXd Jf_x;
@@ -190,7 +191,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
           nested_rev_autodiff rev;
 
           Eigen::VectorXd ret_val = ret.val();
-          auto x_nrad_ = apply(
+          auto x_nrad_ = math::apply(
               [&ret_val, &f, msgs](const auto&... args) {
                 return eval(f(ret_val, msgs, args...));
               },
diff --git a/stan/math/rev/functor/algebra_solver_powell.hpp b/stan/math/rev/functor/algebra_solver_powell.hpp
index 5415f361cad..2d09f97ea8c 100644
--- a/stan/math/rev/functor/algebra_solver_powell.hpp
+++ b/stan/math/rev/functor/algebra_solver_powell.hpp
@@ -136,7 +136,7 @@ Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
   auto args_vals_tuple = std::make_tuple(to_ref(args)...);
 
   auto f_wrt_x = [&args_vals_tuple, &f, msgs](const auto& x) {
-    return apply(
+    return math::apply(
         [&x, &f, msgs](const auto&... args) { return f(x, msgs, args...); },
         args_vals_tuple);
   };
@@ -338,14 +338,14 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_powell_impl(
     const int64_t max_num_steps, const T_Args&... args) {
   auto x_ref = eval(value_of(x));
   auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...));
-  auto args_vals_tuple = apply(
+  auto args_vals_tuple = math::apply(
       [&](const auto&... args) {
         return std::make_tuple(to_ref(value_of(args))...);
       },
       *arena_args_tuple);
 
   auto f_wrt_x = [&args_vals_tuple, &f, msgs](const auto& x) {
-    return apply(
+    return math::apply(
         [&x, &f, msgs](const auto&... args) { return f(x, msgs, args...); },
         args_vals_tuple);
   };
@@ -384,7 +384,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_powell_impl(
     {
       nested_rev_autodiff rev;
       Eigen::VectorXd ret_val = ret.val();
-      auto x_nrad_ = apply(
+      auto x_nrad_ = math::apply(
           [&](const auto&... args) { return eval(f(ret_val, msgs, args...)); },
           *arena_args_tuple);
       x_nrad_.adj() = eta;
diff --git a/stan/math/rev/functor/coupled_ode_system.hpp b/stan/math/rev/functor/coupled_ode_system.hpp
index 1f043fe0e9a..5bb151650cf 100644
--- a/stan/math/rev/functor/coupled_ode_system.hpp
+++ b/stan/math/rev/functor/coupled_ode_system.hpp
@@ -7,6 +7,7 @@
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stdexcept>
 #include <ostream>
@@ -116,7 +117,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
                   double t) {
     using std::vector;
 
-    dz_dt.resize(size());
+    dz_dt.resize(this->size());
 
     // Run nested autodiff in this scope
     nested_rev_autodiff nested;
@@ -125,9 +126,9 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
     for (size_t n = 0; n < N_; ++n)
       y_vars.coeffRef(n) = z[n];
 
-    Eigen::Matrix<var, Eigen::Dynamic, 1> f_y_t_vars
-        = apply([&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
-                local_args_tuple_);
+    Eigen::Matrix<var, Eigen::Dynamic, 1> f_y_t_vars = math::apply(
+        [&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
+        local_args_tuple_);
 
     check_size_match("coupled_ode_system", "dy_dt", f_y_t_vars.size(), "states",
                      N_);
@@ -143,7 +144,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
                sizeof(double) * args_adjoints_.size());
       }
 
-      apply(
+      math::apply(
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },
@@ -189,7 +190,9 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
    *
    * @return size of the coupled system.
    */
-  size_t size() const { return N_ + N_ * num_y0_vars_ + N_ * num_args_vars; }
+  inline Eigen::Index size() const {
+    return N_ + N_ * num_y0_vars_ + N_ * num_args_vars;
+  }
 
   /**
    * Returns the initial state of the coupled system.
@@ -210,7 +213,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
    *   parameters at the initial time-point, which is zero.
    */
   std::vector<double> initial_state() const {
-    std::vector<double> initial(size(), 0.0);
+    std::vector<double> initial(this->size(), 0.0);
     for (size_t i = 0; i < N_; i++) {
       initial[i] = value_of(y0_(i));
     }
diff --git a/stan/math/rev/functor/cvodes_integrator.hpp b/stan/math/rev/functor/cvodes_integrator.hpp
index 637526938e7..473dfad1a3d 100644
--- a/stan/math/rev/functor/cvodes_integrator.hpp
+++ b/stan/math/rev/functor/cvodes_integrator.hpp
@@ -7,6 +7,7 @@
 #include <stan/math/rev/functor/ode_store_sensitivities.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <sundials/sundials_context.h>
 #include <cvodes/cvodes.h>
 #include <nvector/nvector_serial.h>
@@ -104,9 +105,9 @@ class cvodes_integrator {
   inline void rhs(double t, const double y[], double dy_dt[]) const {
     const Eigen::VectorXd y_vec = Eigen::Map<const Eigen::VectorXd>(y, N_);
 
-    Eigen::VectorXd dy_dt_vec
-        = apply([&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
-                value_of_args_tuple_);
+    Eigen::VectorXd dy_dt_vec = math::apply(
+        [&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
+        value_of_args_tuple_);
 
     check_size_match("cvodes_integrator", "dy_dt", dy_dt_vec.size(), "states",
                      N_);
@@ -123,8 +124,9 @@ class cvodes_integrator {
     Eigen::MatrixXd Jfy;
 
     auto f_wrapped = [&](const Eigen::Matrix<var, Eigen::Dynamic, 1>& y) {
-      return apply([&](auto&&... args) { return f_(t, y, msgs_, args...); },
-                   value_of_args_tuple_);
+      return math::apply(
+          [&](auto&&... args) { return f_(t, y, msgs_, args...); },
+          value_of_args_tuple_);
     };
 
     jacobian(f_wrapped, Eigen::Map<const Eigen::VectorXd>(y, N_), fy, Jfy);
@@ -214,7 +216,7 @@ class cvodes_integrator {
 
     // Code from: https://stackoverflow.com/a/17340003 . Should probably do
     // something better
-    apply(
+    math::apply(
         [&](auto&&... args) {
           std::vector<int> unused_temp{
               0, (check_finite(function_name, "ode parameters and data", args),
@@ -313,7 +315,7 @@ class cvodes_integrator {
           }
         }
 
-        y.emplace_back(apply(
+        y.emplace_back(math::apply(
             [&](auto&&... args) {
               return ode_store_sensitivities(f_, coupled_state_, y0_, t0_,
                                              ts_[n], msgs_, args...);
diff --git a/stan/math/rev/functor/cvodes_integrator_adjoint.hpp b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
index 117b869e89f..d698226110e 100644
--- a/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
+++ b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <sundials/sundials_context.h>
 #include <cvodes/cvodes.h>
 #include <nvector/nvector_serial.h>
@@ -544,7 +545,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
   template <typename yT, typename... ArgsT>
   constexpr auto rhs(double t, const yT& y,
                      const std::tuple<ArgsT...>& args_tuple) const {
-    return apply(
+    return math::apply(
         [&](auto&&... args) { return solver_->f_(t, y, msgs_, args...); },
         args_tuple);
   }
@@ -641,7 +642,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
                      f_y_t_vars.size(), "states", N_);
     f_y_t_vars.adj() = -Eigen::Map<Eigen::VectorXd>(NV_DATA_S(yB), N_);
     grad();
-    apply(
+    math::apply(
         [&qBdot](auto&&... args) {
           accumulate_adjoints(NV_DATA_S(qBdot), args...);
         },
diff --git a/stan/math/rev/functor/dae.hpp b/stan/math/rev/functor/dae.hpp
index b3ca6d2c6dc..2d9f40017ad 100644
--- a/stan/math/rev/functor/dae.hpp
+++ b/stan/math/rev/functor/dae.hpp
@@ -66,14 +66,14 @@ dae_tol_impl(const char* func, const F& f, const T_yy& yy0, const T_yp& yp0,
   check_positive(func, "max_num_steps", max_num_steps);
 
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  apply(
+  stan::math::apply(
       [&](auto&&... args) {
         std::vector<int> unused_temp{
             0, (check_finite("dae", "DAE parameters and data", args), 0)...};
       },
       args_ref_tuple);
 
-  return apply(
+  return stan::math::apply(
       [&](const auto&... args_refs) {
         dae_system<F, T_yy, T_yp, ref_type_t<T_Args>...> dae(f, yy0, yp0, msgs,
                                                              args_refs...);
diff --git a/stan/math/rev/functor/dae_system.hpp b/stan/math/rev/functor/dae_system.hpp
index 4e694524b7f..5b05671ab4b 100644
--- a/stan/math/rev/functor/dae_system.hpp
+++ b/stan/math/rev/functor/dae_system.hpp
@@ -99,7 +99,7 @@ class dae_system {
   }
 
   void eval_residual(double t) {
-    dbl_rr = apply(
+    dbl_rr = stan::math::apply(
         [&](auto&&... args) { return f_(t, dbl_yy, dbl_yp, msgs_, args...); },
         dbl_args_tuple_);
   }
@@ -144,7 +144,7 @@ class dae_system {
     }
 
     Eigen::VectorXd g(m);
-    Eigen::Matrix<var, -1, 1> fy = apply(
+    Eigen::Matrix<var, -1, 1> fy = stan::math::apply(
         [&](auto&&... args) {
           return dae->f_(t, yy_var, yp_var, dae->msgs_, args...);
         },
diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp
index be8ff106393..ac4f6edc6d8 100644
--- a/stan/math/rev/functor/integrate_1d.hpp
+++ b/stan/math/rev/functor/integrate_1d.hpp
@@ -9,6 +9,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/functor/integrate_1d.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <cmath>
 #include <functional>
 #include <ostream>
@@ -57,7 +58,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
 
     double integral = integrate(
         [&](const auto &x, const auto &xc) {
-          return apply(
+          return math::apply(
               [&](auto &&... val_args) { return f(x, xc, msgs, val_args...); },
               args_val_tuple);
         },
@@ -79,7 +80,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
     }
 
     if (is_var<T_a>::value && !is_inf(a)) {
-      *partials_ptr = apply(
+      *partials_ptr = math::apply(
           [&f, a_val, msgs](auto &&... val_args) {
             return -f(a_val, 0.0, msgs, val_args...);
           },
@@ -88,7 +89,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
     }
 
     if (!is_inf(b) && is_var<T_b>::value) {
-      *partials_ptr = apply(
+      *partials_ptr = math::apply(
           [&f, b_val, msgs](auto &&... val_args) {
             return f(b_val, 0.0, msgs, val_args...);
           },
@@ -104,7 +105,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
 
       // Save the varis so it's easy to efficiently access the nth adjoint
       std::vector<vari *> local_varis(num_vars_args);
-      apply(
+      math::apply(
           [&](const auto &... args) {
             save_varis(local_varis.data(), args...);
           },
@@ -118,7 +119,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
               argument_nest.set_zero_all_adjoints();
 
               nested_rev_autodiff gradient_nest;
-              var fx = apply(
+              var fx = math::apply(
                   [&f, &x, &xc, msgs](auto &&... local_args) {
                     return f(x, xc, msgs, local_args...);
                   },
diff --git a/stan/math/rev/functor/kinsol_data.hpp b/stan/math/rev/functor/kinsol_data.hpp
index 73f74d9e684..a0847cf9fbd 100644
--- a/stan/math/rev/functor/kinsol_data.hpp
+++ b/stan/math/rev/functor/kinsol_data.hpp
@@ -74,7 +74,7 @@ class kinsol_system_data {
 
     Eigen::Map<Eigen::VectorXd> f_eval_map(N_VGetArrayPointer(f_eval),
                                            explicit_system->N_);
-    auto result = apply(
+    auto result = math::apply(
         [&](const auto&... args) {
           return explicit_system->f_(x_eigen, explicit_system->msgs_, args...);
         },
@@ -109,7 +109,7 @@ class kinsol_system_data {
                                            explicit_system->N_);
 
     auto f_wrt_x = [&](const auto& x) {
-      return apply(
+      return math::apply(
           [&](const auto&... args) {
             return explicit_system->f_(x, explicit_system->msgs_, args...);
           },
diff --git a/stan/math/rev/functor/ode_adams.hpp b/stan/math/rev/functor/ode_adams.hpp
index ee0bdafbbd5..08f6955d767 100644
--- a/stan/math/rev/functor/ode_adams.hpp
+++ b/stan/math/rev/functor/ode_adams.hpp
@@ -1,9 +1,10 @@
 #ifndef STAN_MATH_REV_FUNCTOR_ODE_ADAMS_HPP
 #define STAN_MATH_REV_FUNCTOR_ODE_ADAMS_HPP
 
+#include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/functor/cvodes_integrator.hpp>
-#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <ostream>
 #include <vector>
 
@@ -54,7 +55,7 @@ ode_adams_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                    long int max_num_steps,  // NOLINT(runtime/int)
                    std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return apply(
+  return math::apply(
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_ADAMS, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/ode_bdf.hpp b/stan/math/rev/functor/ode_bdf.hpp
index a07af2e3339..e817ecd9210 100644
--- a/stan/math/rev/functor/ode_bdf.hpp
+++ b/stan/math/rev/functor/ode_bdf.hpp
@@ -1,9 +1,10 @@
 #ifndef STAN_MATH_REV_FUNCTOR_ODE_BDF_HPP
 #define STAN_MATH_REV_FUNCTOR_ODE_BDF_HPP
 
+#include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/functor/cvodes_integrator.hpp>
-#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <ostream>
 #include <vector>
 
@@ -55,7 +56,7 @@ ode_bdf_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                  long int max_num_steps,  // NOLINT(runtime/int)
                  std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return apply(
+  return math::apply(
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_BDF, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/operands_and_partials.hpp b/stan/math/rev/functor/operands_and_partials.hpp
index 908d5803174..d2917ca7e97 100644
--- a/stan/math/rev/functor/operands_and_partials.hpp
+++ b/stan/math/rev/functor/operands_and_partials.hpp
@@ -37,7 +37,7 @@ class ops_partials_edge<double, var> {
   template <typename, typename, typename, typename, typename, typename>
   friend class stan::math::operands_and_partials;
   var operand_;
-  static constexpr int size() noexcept { return 1; }
+  static constexpr Eigen::Index size() noexcept { return 1; }
   inline auto operand() const noexcept { return this->operand_; }
   inline auto partial() const noexcept { return this->partial_; }
 };
@@ -190,7 +190,7 @@ class ops_partials_edge<double, std::vector<var>> {
   friend class stan::math::operands_and_partials;
   Op operands_;
 
-  inline int size() const noexcept { return this->operands_.size(); }
+  inline Eigen::Index size() const noexcept { return this->operands_.size(); }
   inline auto&& operand() noexcept { return std::move(this->operands_); }
   inline auto& partial() noexcept { return this->partials_; }
 };
@@ -210,7 +210,7 @@ class ops_partials_edge<double, Op, require_eigen_st<is_var, Op>> {
   template <typename, typename, typename, typename, typename, typename>
   friend class stan::math::operands_and_partials;
   arena_t<Op> operands_;
-  inline int size() const noexcept { return this->operands_.size(); }
+  inline Eigen::Index size() const noexcept { return this->operands_.size(); }
   inline auto operand() noexcept { return this->operands_; }
   inline auto partial() noexcept { return this->partials_; }
 };
@@ -234,7 +234,7 @@ class ops_partials_edge<double, var_value<Op>, require_eigen_t<Op>> {
 
   static constexpr void dump_operands(vari** varis) {}
   static constexpr void dump_partials(double* partials) {}
-  static constexpr int size() noexcept { return 0; }
+  static constexpr Eigen::Index size() noexcept { return 0; }
   inline auto operand() { return this->operands_; }
   inline auto partial() { return this->partials_; }
 };
@@ -260,7 +260,7 @@ class ops_partials_edge<double, std::vector<Eigen::Matrix<var, R, C>>> {
   friend class stan::math::operands_and_partials;
   Op operands_;
 
-  inline int size() const noexcept {
+  inline Eigen::Index size() const noexcept {
     if (unlikely(this->operands_.size() == 0)) {
       return 0;
     }
@@ -289,7 +289,7 @@ class ops_partials_edge<double, std::vector<std::vector<var>>> {
   template <typename, typename, typename, typename, typename, typename>
   friend class stan::math::operands_and_partials;
   Op operands_;
-  inline int size() const noexcept {
+  inline Eigen::Index size() const noexcept {
     return this->operands_.size() * this->operands_[0].size();
   }
   inline auto&& operand() noexcept { return std::move(this->operands_); }
@@ -315,7 +315,7 @@ class ops_partials_edge<double, std::vector<var_value<Op>>,
   friend class stan::math::operands_and_partials;
   std::vector<var_value<Op>, arena_allocator<var_value<Op>>> operands_;
 
-  static constexpr int size() noexcept { return 0; }
+  static constexpr Eigen::Index size() noexcept { return 0; }
   inline auto&& operand() noexcept { return std::move(this->operands_); }
   inline auto&& partial() noexcept { return std::move(this->partials_vec_); }
 };
diff --git a/stan/math/rev/functor/reduce_sum.hpp b/stan/math/rev/functor/reduce_sum.hpp
index 8002a320117..ba9a626ccd4 100644
--- a/stan/math/rev/functor/reduce_sum.hpp
+++ b/stan/math/rev/functor/reduce_sum.hpp
@@ -1,9 +1,9 @@
 #ifndef STAN_MATH_REV_FUNCTOR_REDUCE_SUM_HPP
 #define STAN_MATH_REV_FUNCTOR_REDUCE_SUM_HPP
 
+#include <stan/math/rev/core.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/functor.hpp>
-#include <stan/math/rev/core.hpp>
 
 #include <tbb/task_arena.h>
 #include <tbb/parallel_reduce.h>
@@ -114,7 +114,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
         // scope. In this case no need for zeroing adjoints, since the
         // fresh copy has all adjoints set to zero.
         local_args_tuple_scope_.stack_.execute([&]() {
-          apply(
+          math::apply(
               [&](auto&&... args) {
                 local_args_tuple_scope_.args_tuple_holder_ = std::make_unique<
                     typename scoped_args_tuple::args_tuple_t>(
@@ -141,7 +141,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
       }
 
       // Perform calculation
-      var sub_sum_v = apply(
+      var sub_sum_v = math::apply(
           [&](auto&&... args) {
             return ReduceFunction()(local_sub_slice, r.begin(), r.end() - 1,
                                     &msgs_, args...);
@@ -159,7 +159,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
                           std::move(local_sub_slice));
 
       // Accumulate adjoints of shared_arguments
-      apply(
+      math::apply(
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },
diff --git a/test/expressions/expression_test_helpers.hpp b/test/expressions/expression_test_helpers.hpp
index d9c373d2e09..5968c552129 100644
--- a/test/expressions/expression_test_helpers.hpp
+++ b/test/expressions/expression_test_helpers.hpp
@@ -1,7 +1,5 @@
+#include <stan/math/mix.hpp>
 #include <gtest/gtest.h>
-#include <stan/math/prim.hpp>
-#include <stan/math/rev.hpp>
-#include <stan/math/fwd.hpp>
 #include <vector>
 #include <random>
 
diff --git a/test/unit/math/opencl/prim/size_test.cpp b/test/unit/math/opencl/prim/size_test.cpp
index 8066203a331..77e15e86c08 100644
--- a/test/unit/math/opencl/prim/size_test.cpp
+++ b/test/unit/math/opencl/prim/size_test.cpp
@@ -7,19 +7,19 @@ TEST(MathMatrixCL, size_prim) {
   using stan::math::matrix_cl;
   using stan::math::size;
   matrix_cl<double> m0;
-  EXPECT_EQ(0, size(m0));
+  EXPECT_EQ(0, stan::math::size(m0));
 
   matrix_cl<double> m1(0, 5);
-  EXPECT_EQ(0, size(m1));
+  EXPECT_EQ(0, stan::math::size(m1));
 
   matrix_cl<double> m2(5, 0);
-  EXPECT_EQ(0, size(m2));
+  EXPECT_EQ(0, stan::math::size(m2));
 
   matrix_cl<double> m3(5, 4);
-  EXPECT_EQ(20, size(m3));
+  EXPECT_EQ(20, stan::math::size(m3));
 
   matrix_cl<double> m4(3, 2);
-  EXPECT_EQ(6, size(m4));
+  EXPECT_EQ(6, stan::math::size(m4));
 }
 
 #endif
diff --git a/test/unit/math/opencl/rev/size_test.cpp b/test/unit/math/opencl/rev/size_test.cpp
index f831360d291..11ceaebd10c 100644
--- a/test/unit/math/opencl/rev/size_test.cpp
+++ b/test/unit/math/opencl/rev/size_test.cpp
@@ -10,23 +10,23 @@ TEST(MathMatrixCL, size_rev) {
   using stan::math::var_value;
   Eigen::Matrix<double, -1, -1> m0(0, 0);
   var_value<matrix_cl<double>> m0_cl = to_matrix_cl(m0);
-  EXPECT_EQ(0, size(m0_cl));
+  EXPECT_EQ(0, stan::math::size(m0_cl));
 
   Eigen::Matrix<double, -1, -1> m1(0, 5);
   var_value<matrix_cl<double>> m1_cl = to_matrix_cl(m1);
-  EXPECT_EQ(0, size(m1_cl));
+  EXPECT_EQ(0, stan::math::size(m1_cl));
 
   Eigen::Matrix<double, -1, -1> m2(5, 0);
   var_value<matrix_cl<double>> m2_cl = to_matrix_cl(m2);
-  EXPECT_EQ(0, size(m2_cl));
+  EXPECT_EQ(0, stan::math::size(m2_cl));
 
   Eigen::Matrix<double, -1, -1> m3(5, 4);
   var_value<matrix_cl<double>> m3_cl = to_matrix_cl(m3);
-  EXPECT_EQ(20, size(m3_cl));
+  EXPECT_EQ(20, stan::math::size(m3_cl));
 
   Eigen::Matrix<double, -1, -1> m4(3, 2);
   var_value<matrix_cl<double>> m4_cl = to_matrix_cl(m4);
-  EXPECT_EQ(6, size(m4_cl));
+  EXPECT_EQ(6, stan::math::size(m4_cl));
 }
 
 #endif
diff --git a/test/unit/math/prim/meta/StdVectorBuilder_test.cpp b/test/unit/math/prim/meta/StdVectorBuilder_test.cpp
index 45707477ec7..cf9d082560d 100644
--- a/test/unit/math/prim/meta/StdVectorBuilder_test.cpp
+++ b/test/unit/math/prim/meta/StdVectorBuilder_test.cpp
@@ -9,7 +9,7 @@ TEST(MathMetaPrim, StdVectorBuilder_false_false_scalar) {
 
   double a_double(1);
 
-  StdVectorBuilder<false, double, double> dvv1(size(a_double));
+  StdVectorBuilder<false, double, double> dvv1(stan::math::size(a_double));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -20,7 +20,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_scalar) {
 
   double a_double(1);
 
-  StdVectorBuilder<true, double, double> dvv1(size(a_double));
+  StdVectorBuilder<true, double, double> dvv1(stan::math::size(a_double));
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   EXPECT_FLOAT_EQ(0.0, dvv1[1]);
   EXPECT_FLOAT_EQ(0.0, dvv1[100]);
@@ -47,7 +47,7 @@ TEST(MathMetaPrim, StdVectorBuilder_false_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  StdVectorBuilder<false, double, double> dvv2(size(a_std_vector));
+  StdVectorBuilder<false, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -59,7 +59,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  StdVectorBuilder<true, double, double> dvv2(size(a_std_vector));
+  StdVectorBuilder<true, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_FLOAT_EQ(0.0, dvv2[0]);
   EXPECT_FLOAT_EQ(0.0, dvv2[1]);
   EXPECT_FLOAT_EQ(0.0, dvv2[2]);
@@ -89,11 +89,11 @@ TEST(MathMetaPrim, StdVectorBuilder_false_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  StdVectorBuilder<false, double, double> dvv3(size(a_vector));
+  StdVectorBuilder<false, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  StdVectorBuilder<false, double, double> dvv4(size(a_row_vector));
+  StdVectorBuilder<false, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv4.data(), std::logic_error);
 }
@@ -107,7 +107,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  StdVectorBuilder<true, double, double> dvv3(size(a_vector));
+  StdVectorBuilder<true, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_FLOAT_EQ(0.0, dvv3[0]);
   EXPECT_FLOAT_EQ(0.0, dvv3[1]);
   EXPECT_FLOAT_EQ(0.0, dvv3[2]);
@@ -115,7 +115,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_matrix) {
   EXPECT_NO_THROW(data3 = dvv3.data());
   EXPECT_FLOAT_EQ(0.0, data3);
 
-  StdVectorBuilder<true, double, double> dvv4(size(a_row_vector));
+  StdVectorBuilder<true, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_FLOAT_EQ(0.0, dvv4[0]);
   EXPECT_FLOAT_EQ(0.0, dvv4[1]);
   EXPECT_FLOAT_EQ(0.0, dvv4[2]);
diff --git a/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp b/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp
index 8c0a105f4fe..a8be33072b4 100644
--- a/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp
+++ b/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp
@@ -9,7 +9,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_false_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilderHelper<double, false, false> dvv1(size(a_double));
+  VectorBuilderHelper<double, false, false> dvv1(stan::math::size(a_double));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -20,7 +20,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilderHelper<double, true, false> dvv1(size(a_double));
+  VectorBuilderHelper<double, true, false> dvv1(stan::math::size(a_double));
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   EXPECT_FLOAT_EQ(0.0, dvv1[1]);
   EXPECT_FLOAT_EQ(0.0, dvv1[100]);
@@ -36,7 +36,8 @@ TEST(MathMetaPrim, VectorBuilderHelper_false_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilderHelper<double, false, false> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, false, false> dvv2(
+      stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -48,7 +49,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilderHelper<double, true, false> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, true, false> dvv2(stan::math::size(a_std_vector));
   EXPECT_FLOAT_EQ(0.0, dvv2[0]);
   EXPECT_FLOAT_EQ(0.0, dvv2[1]);
   EXPECT_FLOAT_EQ(0.0, dvv2[2]);
@@ -66,11 +67,12 @@ TEST(MathMetaPrim, VectorBuilderHelper_false_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, false, false> dvv3(size(a_vector));
+  VectorBuilderHelper<double, false, false> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  VectorBuilderHelper<double, false, false> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, false, false> dvv4(
+      stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv4.data(), std::logic_error);
 }
@@ -84,7 +86,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, true, false> dvv3(size(a_vector));
+  VectorBuilderHelper<double, true, false> dvv3(stan::math::size(a_vector));
   EXPECT_FLOAT_EQ(0.0, dvv3[0]);
   EXPECT_FLOAT_EQ(0.0, dvv3[1]);
   EXPECT_FLOAT_EQ(0.0, dvv3[2]);
@@ -92,7 +94,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_matrix) {
   EXPECT_NO_THROW(data3 = dvv3.data());
   EXPECT_FLOAT_EQ(0.0, data3);
 
-  VectorBuilderHelper<double, true, false> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, true, false> dvv4(stan::math::size(a_row_vector));
   EXPECT_FLOAT_EQ(0.0, dvv4[0]);
   EXPECT_FLOAT_EQ(0.0, dvv4[1]);
   EXPECT_FLOAT_EQ(0.0, dvv4[2]);
diff --git a/test/unit/math/prim/meta/VectorBuilder_test.cpp b/test/unit/math/prim/meta/VectorBuilder_test.cpp
index a1dca9ae971..e48eb3dbc74 100644
--- a/test/unit/math/prim/meta/VectorBuilder_test.cpp
+++ b/test/unit/math/prim/meta/VectorBuilder_test.cpp
@@ -9,7 +9,7 @@ TEST(MathMetaPrim, VectorBuilder_false_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilder<false, double, double> dvv1(size(a_double));
+  VectorBuilder<false, double, double> dvv1(stan::math::size(a_double));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -20,7 +20,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilder<true, double, double> dvv1(size(a_double));
+  VectorBuilder<true, double, double> dvv1(stan::math::size(a_double));
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   EXPECT_FLOAT_EQ(0.0, dvv1[1]);
   EXPECT_FLOAT_EQ(0.0, dvv1[100]);
@@ -36,7 +36,7 @@ TEST(MathMetaPrim, VectorBuilder_false_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilder<false, double, double> dvv2(size(a_std_vector));
+  VectorBuilder<false, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -48,7 +48,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilder<true, double, double> dvv2(size(a_std_vector));
+  VectorBuilder<true, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_FLOAT_EQ(0.0, dvv2[0]);
   EXPECT_FLOAT_EQ(0.0, dvv2[1]);
   EXPECT_FLOAT_EQ(0.0, dvv2[2]);
@@ -66,11 +66,11 @@ TEST(MathMetaPrim, VectorBuilder_false_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<false, double, double> dvv3(size(a_vector));
+  VectorBuilder<false, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  VectorBuilder<false, double, double> dvv4(size(a_row_vector));
+  VectorBuilder<false, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv4.data(), std::logic_error);
 }
@@ -84,7 +84,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<true, double, double> dvv3(size(a_vector));
+  VectorBuilder<true, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_FLOAT_EQ(0.0, dvv3[0]);
   EXPECT_FLOAT_EQ(0.0, dvv3[1]);
   EXPECT_FLOAT_EQ(0.0, dvv3[2]);
@@ -92,7 +92,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_matrix) {
   EXPECT_NO_THROW(data3 = dvv3.data());
   EXPECT_FLOAT_EQ(0.0, data3);
 
-  VectorBuilder<true, double, double> dvv4(size(a_row_vector));
+  VectorBuilder<true, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_FLOAT_EQ(0.0, dvv4[0]);
   EXPECT_FLOAT_EQ(0.0, dvv4[1]);
   EXPECT_FLOAT_EQ(0.0, dvv4[2]);
diff --git a/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp b/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp
index ac845bbfe35..af0d0234e78 100644
--- a/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp
+++ b/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp
@@ -11,7 +11,7 @@ TEST(MetaTraitsRevScal, VectorBuilderHelper_false_true) {
 
   var a_var(1);
 
-  VectorBuilderHelper<double, false, true> dvv1(size(a_var));
+  VectorBuilderHelper<double, false, true> dvv1(stan::math::size(a_var));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -24,7 +24,7 @@ TEST(MetaTraitsRevArr, VectorBuilderHelper_false_true) {
 
   std::vector<var> a_std_vector(3);
 
-  VectorBuilderHelper<double, false, true> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, false, true> dvv2(stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -37,14 +37,14 @@ TEST(MetaTraitsRevArr, VectorBuilderHelper_true_true) {
 
   var a_var(1);
   std::vector<var> a_std_vector(3);
-  VectorBuilderHelper<double, true, true> dvv1(size(a_var));
+  VectorBuilderHelper<double, true, true> dvv1(stan::math::size(a_var));
   dvv1[0] = 0.0;
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   std::vector<double> data1;
   EXPECT_NO_THROW(data1 = dvv1.data());
-  EXPECT_EQ(size(a_var), data1.size());
+  EXPECT_EQ(stan::math::size(a_var), data1.size());
 
-  VectorBuilderHelper<double, true, true> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, true, true> dvv2(stan::math::size(a_std_vector));
   dvv2[0] = 0.0;
   dvv2[1] = 1.0;
   dvv2[2] = 2.0;
@@ -53,7 +53,7 @@ TEST(MetaTraitsRevArr, VectorBuilderHelper_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv2[2]);
   std::vector<double> data2;
   EXPECT_NO_THROW(data2 = dvv2.data());
-  EXPECT_EQ(size(a_std_vector), data2.size());
+  EXPECT_EQ(stan::math::size(a_std_vector), data2.size());
 }
 
 TEST(MetaTraitsRevMat, VectorBuilderHelper_false_true) {
@@ -66,11 +66,11 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_false_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, false, true> dvv3(size(a_vector));
+  VectorBuilderHelper<double, false, true> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  VectorBuilderHelper<double, false, true> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, false, true> dvv4(stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 }
@@ -85,7 +85,7 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_true_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, true, true> dvv3(size(a_vector));
+  VectorBuilderHelper<double, true, true> dvv3(stan::math::size(a_vector));
   dvv3[0] = 0.0;
   dvv3[1] = 1.0;
   dvv3[2] = 2.0;
@@ -94,9 +94,9 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv3[2]);
   std::vector<double> data3;
   EXPECT_NO_THROW(data3 = dvv3.data());
-  EXPECT_EQ(size(a_vector), data3.size());
+  EXPECT_EQ(stan::math::size(a_vector), data3.size());
 
-  VectorBuilderHelper<double, true, true> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, true, true> dvv4(stan::math::size(a_row_vector));
   dvv4[0] = 0.0;
   dvv4[1] = 1.0;
   dvv4[2] = 2.0;
@@ -105,5 +105,5 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv4[2]);
   std::vector<double> data4;
   EXPECT_NO_THROW(data4 = dvv4.data());
-  EXPECT_EQ(size(a_row_vector), data4.size());
+  EXPECT_EQ(stan::math::size(a_row_vector), data4.size());
 }
diff --git a/test/unit/math/rev/meta/VectorBuilder_test.cpp b/test/unit/math/rev/meta/VectorBuilder_test.cpp
index 481870f24a8..edd3df52fdc 100644
--- a/test/unit/math/rev/meta/VectorBuilder_test.cpp
+++ b/test/unit/math/rev/meta/VectorBuilder_test.cpp
@@ -13,11 +13,12 @@ TEST(MetaTraitsRevArr, VectorBuilder_false_true) {
   var a_var(1);
   std::vector<var> a_std_vector(3);
 
-  VectorBuilder<false, double, std::vector<var> > dvv1(size(a_var));
+  VectorBuilder<false, double, std::vector<var> > dvv1(stan::math::size(a_var));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 
-  VectorBuilder<false, double, std::vector<var> > dvv2(size(a_std_vector));
+  VectorBuilder<false, double, std::vector<var> > dvv2(
+      stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -31,14 +32,15 @@ TEST(MetaTraitsRevArr, VectorBuilder_true_true) {
   var a_var(1);
   std::vector<var> a_std_vector(3);
 
-  VectorBuilder<true, double, std::vector<var> > dvv1(size(a_var));
+  VectorBuilder<true, double, std::vector<var> > dvv1(stan::math::size(a_var));
   dvv1[0] = 0.0;
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   std::vector<double> data1;
   EXPECT_NO_THROW(data1 = dvv1.data());
-  EXPECT_EQ(size(a_var), data1.size());
+  EXPECT_EQ(stan::math::size(a_var), data1.size());
 
-  VectorBuilder<true, double, std::vector<var> > dvv2(size(a_std_vector));
+  VectorBuilder<true, double, std::vector<var> > dvv2(
+      stan::math::size(a_std_vector));
   dvv2[0] = 0.0;
   dvv2[1] = 1.0;
   dvv2[2] = 2.0;
@@ -47,7 +49,7 @@ TEST(MetaTraitsRevArr, VectorBuilder_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv2[2]);
   std::vector<double> data2;
   EXPECT_NO_THROW(data2 = dvv2.data());
-  EXPECT_EQ(size(a_std_vector), data2.size());
+  EXPECT_EQ(stan::math::size(a_std_vector), data2.size());
 }
 
 TEST(MetaTraitsRevMat, VectorBuilder_false_true) {
@@ -60,7 +62,8 @@ TEST(MetaTraitsRevMat, VectorBuilder_false_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<false, double, Matrix<var, Dynamic, 1> > dvv3(size(a_vector));
+  VectorBuilder<false, double, Matrix<var, Dynamic, 1> > dvv3(
+      stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
@@ -80,7 +83,8 @@ TEST(MetaTraitsRevMat, VectorBuilder_true_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<true, double, Matrix<var, Dynamic, 1> > dvv3(size(a_vector));
+  VectorBuilder<true, double, Matrix<var, Dynamic, 1> > dvv3(
+      stan::math::size(a_vector));
   dvv3[0] = 0.0;
   dvv3[1] = 1.0;
   dvv3[2] = 2.0;
@@ -89,7 +93,7 @@ TEST(MetaTraitsRevMat, VectorBuilder_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv3[2]);
   std::vector<double> data3;
   EXPECT_NO_THROW(data3 = dvv3.data());
-  EXPECT_EQ(size(a_vector), data3.size());
+  EXPECT_EQ(stan::math::size(a_vector), data3.size());
 
   VectorBuilder<true, double, Matrix<var, 1, Dynamic> > dvv4(
       stan::math::size(a_row_vector));
@@ -101,5 +105,5 @@ TEST(MetaTraitsRevMat, VectorBuilder_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv4[2]);
   std::vector<double> data4;
   EXPECT_NO_THROW(data4 = dvv4.data());
-  EXPECT_EQ(size(a_row_vector), data4.size());
+  EXPECT_EQ(stan::math::size(a_row_vector), data4.size());
 }