diff --git a/R/fit.R b/R/fit.R
index 4e15e167..5f07c5ee 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -326,13 +326,12 @@ CmdStanFit$set("public", name = "init", value = init)
 #' @examples
 #' \dontrun{
 #' fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-#' fit_mcmc$init_model_methods()
 #' }
 #' @seealso [log_prob()], [grad_log_prob()], [constrain_variables()],
 #'   [unconstrain_variables()], [unconstrain_draws()], [variable_skeleton()],
 #'   [hessian()]
 #'
-init_model_methods <- function(seed = 0, verbose = FALSE, hessian = FALSE) {
+init_model_methods <- function(seed = 1, verbose = FALSE, hessian = FALSE) {
   if (os_is_wsl()) {
     stop("Additional model methods are not currently available with ",
           "WSL CmdStan and will not be compiled",
@@ -348,11 +347,12 @@ init_model_methods <- function(seed = 0, verbose = FALSE, hessian = FALSE) {
             "which is still experimental. Please report any compilation ",
             "errors that you encounter")
   }
-  message("Compiling additional model methods...")
   if (is.null(private$model_methods_env_$model_ptr)) {
     expose_model_methods(private$model_methods_env_, verbose, hessian)
   }
-  initialize_model_pointer(private$model_methods_env_, self$data_file(), seed)
+  if (!("model_ptr_" %in% ls(private$model_methods_env_))) {
+    initialize_model_pointer(private$model_methods_env_, self$data_file(), seed)
+  }
   invisible(NULL)
 }
 CmdStanFit$set("public", name = "init_model_methods", value = init_model_methods)
@@ -372,7 +372,6 @@ CmdStanFit$set("public", name = "init_model_methods", value = init_model_methods
 #' @examples
 #' \dontrun{
 #' fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-#' fit_mcmc$init_model_methods()
 #' fit_mcmc$log_prob(unconstrained_variables = c(0.5, 1.2, 1.1, 2.2))
 #' }
 #'
@@ -385,10 +384,7 @@ log_prob <- function(unconstrained_variables, jacobian = TRUE, jacobian_adjustme
     warning("'jacobian_adjustment' is deprecated. Please use 'jacobian' instead.", call. = FALSE)
     jacobian <- jacobian_adjustment
   }
-  if (is.null(private$model_methods_env_$model_ptr)) {
-    stop("The method has not been compiled, please call `init_model_methods()` first",
-        call. = FALSE)
-  }
+  self$init_model_methods()
   if (length(unconstrained_variables) != private$model_methods_env_$num_upars_) {
     stop("Model has ", private$model_methods_env_$num_upars_, " unconstrained parameter(s), but ",
           length(unconstrained_variables), " were provided!", call. = FALSE)
@@ -410,7 +406,6 @@ CmdStanFit$set("public", name = "log_prob", value = log_prob)
 #' @examples
 #' \dontrun{
 #' fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-#' fit_mcmc$init_model_methods()
 #' fit_mcmc$grad_log_prob(unconstrained_variables = c(0.5, 1.2, 1.1, 2.2))
 #' }
 #'
@@ -423,10 +418,7 @@ grad_log_prob <- function(unconstrained_variables, jacobian = TRUE, jacobian_adj
     warning("'jacobian_adjustment' is deprecated. Please use 'jacobian' instead.", call. = FALSE)
     jacobian <- jacobian_adjustment
   }
-  if (is.null(private$model_methods_env_$model_ptr)) {
-    stop("The method has not been compiled, please call `init_model_methods()` first",
-        call. = FALSE)
-  }
+  self$init_model_methods()
   if (length(unconstrained_variables) != private$model_methods_env_$num_upars_) {
     stop("Model has ", private$model_methods_env_$num_upars_, " unconstrained parameter(s), but ",
           length(unconstrained_variables), " were provided!", call. = FALSE)
@@ -461,10 +453,7 @@ hessian <- function(unconstrained_variables, jacobian = TRUE, jacobian_adjustmen
     warning("'jacobian_adjustment' is deprecated. Please use 'jacobian' instead.", call. = FALSE)
     jacobian <- jacobian_adjustment
   }
-  if (is.null(private$model_methods_env_$model_ptr)) {
-    stop("The method has not been compiled, please call `init_model_methods()` first",
-        call. = FALSE)
-  }
+  self$init_model_methods()
   if (length(unconstrained_variables) != private$model_methods_env_$num_upars_) {
     stop("Model has ", private$model_methods_env_$num_upars_, " unconstrained parameter(s), but ",
           length(unconstrained_variables), " were provided!", call. = FALSE)
@@ -487,7 +476,6 @@ CmdStanFit$set("public", name = "hessian", value = hessian)
 #' @examples
 #' \dontrun{
 #' fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-#' fit_mcmc$init_model_methods()
 #' fit_mcmc$unconstrain_variables(list(alpha = 0.5, beta = c(0.7, 1.1, 0.2)))
 #' }
 #'
@@ -496,10 +484,7 @@ CmdStanFit$set("public", name = "hessian", value = hessian)
 #'   [hessian()]
 #'
 unconstrain_variables <- function(variables) {
-  if (is.null(private$model_methods_env_$model_ptr)) {
-    stop("The method has not been compiled, please call `init_model_methods()` first",
-        call. = FALSE)
-  }
+  self$init_model_methods()
   model_par_names <- self$metadata()$stan_variables[self$metadata()$stan_variables != "lp__"]
   prov_par_names <- names(variables)
 
@@ -539,11 +524,12 @@ CmdStanFit$set("public", name = "unconstrain_variables", value = unconstrain_var
 #' @param draws A `posterior::draws_*` object.
 #' @param format (string) The format of the returned draws. Must be a valid
 #'   format from the \pkg{posterior} package.
+#' @param inc_warmup (logical) Should warmup draws be included? Defaults to
+#'  `FALSE`.
 #'
 #' @examples
 #' \dontrun{
 #' fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-#' fit_mcmc$init_model_methods()
 #'
 #' # Unconstrain all internal draws
 #' unconstrained_internal_draws <- fit_mcmc$unconstrain_draws()
@@ -560,7 +546,9 @@ CmdStanFit$set("public", name = "unconstrain_variables", value = unconstrain_var
 #'   [hessian()]
 #'
 unconstrain_draws <- function(files = NULL, draws = NULL,
-                              format = getOption("cmdstanr_draws_format", "draws_array")) {
+                              format = getOption("cmdstanr_draws_format", "draws_array"),
+                              inc_warmup = FALSE) {
+  self$init_model_methods()
   if (!(format %in% valid_draws_formats())) {
     stop("Invalid draws format requested!", call. = FALSE)
   }
@@ -570,22 +558,25 @@ unconstrain_draws <- function(files = NULL, draws = NULL,
           call. = FALSE)
     }
     if (!is.null(files)) {
-      read_csv <- read_cmdstan_csv(files = files, format = "draws_matrix")
-      draws <- read_csv$post_warmup_draws
-    }
-    if (!is.null(draws)) {
-      draws <- maybe_convert_draws_format(draws, "draws_matrix")
-    }
-  } else {
-    if (is.null(private$draws_)) {
-      if (!length(self$output_files(include_failed = FALSE))) {
-        stop("Fitting failed. Unable to retrieve the draws.", call. = FALSE)
+      read_csv <- read_cmdstan_csv(files = files)
+      if (inc_warmup) {
+        draws <- posterior::bind_draws(read_csv$warmup_draws,
+                                        read_csv$post_warmup_draws,
+                                        along = "iteration")
+      } else {
+        draws <- read_csv$post_warmup_draws
+      }
+    } else if (!is.null(draws)) {
+      if (inc_warmup) {
+        message("'inc_warmup' cannot be used with a draws object. Ignoring.")
       }
-      private$read_csv_(format = "draws_df")
     }
-    draws <- maybe_convert_draws_format(private$draws_, "draws_matrix")
+  } else {
+    draws <- self$draws(inc_warmup = inc_warmup)
   }
   
+  draws <- maybe_convert_draws_format(draws, "draws_matrix")
+  
   chains <- posterior::nchains(draws)
 
   model_par_names <- self$metadata()$stan_variables[self$metadata()$stan_variables != "lp__"]
@@ -624,7 +615,6 @@ CmdStanFit$set("public", name = "unconstrain_draws", value = unconstrain_draws)
 #' @examples
 #' \dontrun{
 #' fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-#' fit_mcmc$init_model_methods()
 #' fit_mcmc$variable_skeleton()
 #' }
 #'
@@ -633,11 +623,7 @@ CmdStanFit$set("public", name = "unconstrain_draws", value = unconstrain_draws)
 #'   [hessian()]
 #'
 variable_skeleton <- function(transformed_parameters = TRUE, generated_quantities = TRUE) {
-  if (is.null(private$model_methods_env_$model_ptr)) {
-    stop("The method has not been compiled, please call `init_model_methods()` first",
-        call. = FALSE)
-  }
-
+  self$init_model_methods()
   create_skeleton(private$model_methods_env_$param_metadata_,
                   self$runset$args$model_variables,
                   transformed_parameters,
@@ -662,7 +648,6 @@ CmdStanFit$set("public", name = "variable_skeleton", value = variable_skeleton)
 #' @examples
 #' \dontrun{
 #' fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-#' fit_mcmc$init_model_methods()
 #' fit_mcmc$constrain_variables(unconstrained_variables = c(0.5, 1.2, 1.1, 2.2))
 #' }
 #'
@@ -671,12 +656,8 @@ CmdStanFit$set("public", name = "variable_skeleton", value = variable_skeleton)
 #'   [hessian()]
 #'
 constrain_variables <- function(unconstrained_variables, transformed_parameters = TRUE,
-                            generated_quantities = TRUE) {
-  if (is.null(private$model_methods_env_$model_ptr)) {
-    stop("The method has not been compiled, please call `init_model_methods()` first",
-        call. = FALSE)
-  }
-
+                                generated_quantities = TRUE) {
+  self$init_model_methods()
   skeleton <- self$variable_skeleton(transformed_parameters, generated_quantities)
 
   if (length(unconstrained_variables) != private$model_methods_env_$num_upars_) {
diff --git a/R/utils.R b/R/utils.R
index 0fb77245..7485db14 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -786,6 +786,9 @@ rcpp_source_stan <- function(code, env, verbose = FALSE, ...) {
 }
 
 expose_model_methods <- function(env, verbose = FALSE, hessian = FALSE) {
+  if (rlang::is_interactive()) {
+    message("Compiling additional model methods...")
+  }
   code <- c(env$hpp_code_,
             readLines(system.file("include", "model_methods.cpp",
                                   package = "cmdstanr", mustWork = TRUE)))
@@ -1034,7 +1037,9 @@ expose_stan_functions <- function(function_env, global = FALSE, verbose = FALSE)
       })
     }
   } else {
-    message("Compiling standalone functions...")
+    if (rlang::is_interactive()) {
+      message("Compiling standalone functions...")
+    }
     compile_functions(function_env, verbose, global)
   }
   invisible(NULL)
diff --git a/man/fit-method-constrain_variables.Rd b/man/fit-method-constrain_variables.Rd
index 56e0386f..d46b3712 100644
--- a/man/fit-method-constrain_variables.Rd
+++ b/man/fit-method-constrain_variables.Rd
@@ -28,7 +28,6 @@ to the constrained scale.
 \examples{
 \dontrun{
 fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-fit_mcmc$init_model_methods()
 fit_mcmc$constrain_variables(unconstrained_variables = c(0.5, 1.2, 1.1, 2.2))
 }
 
diff --git a/man/fit-method-grad_log_prob.Rd b/man/fit-method-grad_log_prob.Rd
index b74f19ba..5012d79c 100644
--- a/man/fit-method-grad_log_prob.Rd
+++ b/man/fit-method-grad_log_prob.Rd
@@ -27,7 +27,6 @@ model's \code{log_prob} function and its derivative.
 \examples{
 \dontrun{
 fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-fit_mcmc$init_model_methods()
 fit_mcmc$grad_log_prob(unconstrained_variables = c(0.5, 1.2, 1.1, 2.2))
 }
 
diff --git a/man/fit-method-init_model_methods.Rd b/man/fit-method-init_model_methods.Rd
index 1a96cd59..64fa8a76 100644
--- a/man/fit-method-init_model_methods.Rd
+++ b/man/fit-method-init_model_methods.Rd
@@ -6,7 +6,7 @@
 \title{Compile additional methods for accessing the model log-probability function
 and parameter constraining and unconstraining.}
 \usage{
-init_model_methods(seed = 0, verbose = FALSE, hessian = FALSE)
+init_model_methods(seed = 1, verbose = FALSE, hessian = FALSE)
 }
 \arguments{
 \item{seed}{(integer) The random seed to use when initializing the model.}
@@ -29,7 +29,6 @@ these can be ignored so long as they are warnings and not errors.
 \examples{
 \dontrun{
 fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-fit_mcmc$init_model_methods()
 }
 }
 \seealso{
diff --git a/man/fit-method-log_prob.Rd b/man/fit-method-log_prob.Rd
index 41edbc3c..456fb82c 100644
--- a/man/fit-method-log_prob.Rd
+++ b/man/fit-method-log_prob.Rd
@@ -22,7 +22,6 @@ The \verb{$log_prob()} method provides access to the Stan model's
 \examples{
 \dontrun{
 fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-fit_mcmc$init_model_methods()
 fit_mcmc$log_prob(unconstrained_variables = c(0.5, 1.2, 1.1, 2.2))
 }
 
diff --git a/man/fit-method-unconstrain_draws.Rd b/man/fit-method-unconstrain_draws.Rd
index 999f3339..9f4d0e6f 100644
--- a/man/fit-method-unconstrain_draws.Rd
+++ b/man/fit-method-unconstrain_draws.Rd
@@ -8,7 +8,8 @@
 unconstrain_draws(
   files = NULL,
   draws = NULL,
-  format = getOption("cmdstanr_draws_format", "draws_array")
+  format = getOption("cmdstanr_draws_format", "draws_array"),
+  inc_warmup = FALSE
 )
 }
 \arguments{
@@ -19,6 +20,9 @@ be files generated by running CmdStanR or running CmdStan directly.}
 
 \item{format}{(string) The format of the returned draws. Must be a valid
 format from the \pkg{posterior} package.}
+
+\item{inc_warmup}{(logical) Should warmup draws be included? Defaults to
+\code{FALSE}.}
 }
 \description{
 The \verb{$unconstrain_draws()} method transforms all parameter draws
@@ -31,7 +35,6 @@ character vector of paths to CSV files can be passed.
 \examples{
 \dontrun{
 fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-fit_mcmc$init_model_methods()
 
 # Unconstrain all internal draws
 unconstrained_internal_draws <- fit_mcmc$unconstrain_draws()
diff --git a/man/fit-method-unconstrain_variables.Rd b/man/fit-method-unconstrain_variables.Rd
index cd546526..8ac1eef2 100644
--- a/man/fit-method-unconstrain_variables.Rd
+++ b/man/fit-method-unconstrain_variables.Rd
@@ -18,7 +18,6 @@ parameters to the unconstrained scale.
 \examples{
 \dontrun{
 fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-fit_mcmc$init_model_methods()
 fit_mcmc$unconstrain_variables(list(alpha = 0.5, beta = c(0.7, 1.1, 0.2)))
 }
 
diff --git a/man/fit-method-variable_skeleton.Rd b/man/fit-method-variable_skeleton.Rd
index 3fc24b12..df7801c1 100644
--- a/man/fit-method-variable_skeleton.Rd
+++ b/man/fit-method-variable_skeleton.Rd
@@ -22,7 +22,6 @@ parameter values to a named list.
 \examples{
 \dontrun{
 fit_mcmc <- cmdstanr_example("logistic", method = "sample", force_recompile = TRUE)
-fit_mcmc$init_model_methods()
 fit_mcmc$variable_skeleton()
 }
 
diff --git a/tests/testthat/test-model-methods.R b/tests/testthat/test-model-methods.R
index 276d0b26..4ad98ef0 100644
--- a/tests/testthat/test-model-methods.R
+++ b/tests/testthat/test-model-methods.R
@@ -1,43 +1,16 @@
 context("model-methods")
+skip_if(os_is_wsl())
 
 set_cmdstan_path()
 mod <- cmdstan_model(testing_stan_file("bernoulli_log_lik"), force_recompile = TRUE)
 data_list <- testing_data("bernoulli")
 fit <- mod$sample(data = data_list, chains = 1, refresh = 0)
 
-test_that("Methods error if not compiled", {
-  skip_if(os_is_wsl())
-  expect_error(
-    fit$log_prob(NULL),
-    "The method has not been compiled, please call `init_model_methods()` first",
-    fixed = TRUE
-  )
-  expect_error(
-    fit$grad_log_prob(NULL),
-    "The method has not been compiled, please call `init_model_methods()` first",
-    fixed = TRUE
-  )
-  expect_error(
-    fit$hessian(NULL),
-    "The method has not been compiled, please call `init_model_methods()` first",
-    fixed = TRUE
-  )
-  expect_error(
-    fit$unconstrain_variables(NULL),
-    "The method has not been compiled, please call `init_model_methods()` first",
-    fixed = TRUE
-  )
-  expect_error(
-    fit$constrain_variables(NULL),
-    "The method has not been compiled, please call `init_model_methods()` first",
-    fixed = TRUE
-  )
+test_that("Model methods automatically initialise when needed", {
+  expect_no_error(fit$log_prob(unconstrained_variables=c(0.1)))
 })
 
-
 test_that("Methods return correct values", {
-  skip_if(os_is_wsl())
-  fit$init_model_methods(verbose = TRUE)
   lp <- fit$log_prob(unconstrained_variables=c(0.1))
   expect_equal(lp, -8.6327599208828509347)
 
@@ -83,7 +56,6 @@ test_that("Methods return correct values", {
 })
 
 test_that("Model methods environments are independent", {
-  skip_if(os_is_wsl())
   data_list_2 <- data_list
   data_list_2$N <- 20
   data_list_2$y <- c(data_list$y, data_list$y)
@@ -95,7 +67,6 @@ test_that("Model methods environments are independent", {
 })
 
 test_that("methods error for incorrect inputs", {
-  skip_if(os_is_wsl())
   expect_error(
     fit$log_prob(c(1,2)),
     "Model has 1 unconstrained parameter(s), but 2 were provided!",
@@ -120,8 +91,6 @@ test_that("methods error for incorrect inputs", {
   logistic_mod <- cmdstan_model(testing_stan_file("logistic"), force_recompile = TRUE)
   logistic_data_list <- testing_data("logistic")
   logistic_fit <- logistic_mod$sample(data = logistic_data_list, chains = 1)
-  # Init without Hessian, as bernoulli_logit_glm currently not fully fvar<var>
-  # compatible
   logistic_fit$init_model_methods(verbose = TRUE)
 
   expect_error(
@@ -132,7 +101,6 @@ test_that("methods error for incorrect inputs", {
 })
 
 test_that("Methods error with already-compiled model", {
-  skip_if(os_is_wsl())
   precompile_mod <- testing_model("bernoulli")
   mod <- testing_model("bernoulli")
   data_list <- testing_data("bernoulli")
@@ -145,11 +113,9 @@ test_that("Methods error with already-compiled model", {
 })
 
 test_that("Methods can be compiled with model", {
-  skip_if(os_is_wsl())
   mod <- cmdstan_model(testing_stan_file("bernoulli"),
                        force_recompile = TRUE,
-                       compile_model_methods = TRUE,
-                       compile_hessian_method = TRUE)
+                       compile_model_methods = TRUE)
   fit <- mod$sample(data = data_list, chains = 1)
 
   lp <- fit$log_prob(unconstrained_variables=c(0.6))
@@ -174,7 +140,6 @@ test_that("Methods can be compiled with model", {
 })
 
 test_that("unconstrain_variables correctly handles zero-length containers", {
-  skip_if(os_is_wsl())
   model_code <- "
   data {
     int N;
@@ -189,6 +154,7 @@ test_that("unconstrain_variables correctly handles zero-length containers", {
   }
   "
   mod <- cmdstan_model(write_stan_file(model_code),
+                       force_recompile = TRUE,
                        compile_model_methods = TRUE)
   fit <- mod$sample(data = list(N = 0), chains = 1)
   unconstrained <- fit$unconstrain_variables(variables = list(x = 5))
@@ -196,7 +162,6 @@ test_that("unconstrain_variables correctly handles zero-length containers", {
 })
 
 test_that("unconstrain_draws returns correct values", {
-  skip_if(os_is_wsl())
 
   # With no constraints, the parameter draws should be the same as the
   # unconstrained draws
@@ -214,21 +179,34 @@ test_that("unconstrain_draws returns correct values", {
   mod <- cmdstan_model(write_stan_file(model_code),
                        compile_model_methods = TRUE,
                        force_recompile = TRUE)
-  fit <- mod$sample(data = list(N = 0), chains = 2)
+  fit <- mod$sample(data = list(N = 0), chains = 2, save_warmup = TRUE)
+  fit_no_warmup <- mod$sample(data = list(N = 0), chains = 2)
 
   x_draws <- fit$draws(format = "draws_df")$x
-
+  x_draws_warmup <- fit$draws(format = "draws_df", inc_warmup = TRUE)$x
+  
   # Unconstrain all internal draws
   unconstrained_internal_draws <- fit$unconstrain_draws()
+  unconstrained_internal_draws_warmup <- fit$unconstrain_draws(inc_warmup = TRUE)
   expect_equal(as.numeric(x_draws), as.numeric(unconstrained_internal_draws))
-
+  expect_equal(as.numeric(x_draws_warmup), as.numeric(unconstrained_internal_draws_warmup))
+  
+  expect_error({unconstrained_internal_draws <- fit_no_warmup$unconstrain_draws(inc_warmup = TRUE)},
+               "Warmup draws were requested from a fit object without them! Please rerun the model with save_warmup = TRUE.")
+  
   # Unconstrain external CmdStan CSV files
   unconstrained_csv <- fit$unconstrain_draws(files = fit$output_files())
+  unconstrained_csv_warmup <- fit$unconstrain_draws(files = fit$output_files(),
+                                                    inc_warmup = TRUE)
   expect_equal(as.numeric(x_draws), as.numeric(unconstrained_csv))
+  expect_equal(as.numeric(x_draws_warmup), as.numeric(unconstrained_csv_warmup))
 
   # Unconstrain existing draws object
   unconstrained_draws <- fit$unconstrain_draws(draws = fit$draws())
   expect_equal(as.numeric(x_draws), as.numeric(unconstrained_draws))
+  
+  expect_message(fit$unconstrain_draws(draws = fit$draws(), inc_warmup = TRUE),
+                 "'inc_warmup' cannot be used with a draws object. Ignoring.")
 
   # With a lower-bounded constraint, the parameter draws should be the
   # exponentiation of the unconstrained draws
@@ -263,7 +241,6 @@ test_that("unconstrain_draws returns correct values", {
 })
 
 test_that("Model methods can be initialised for models with no data", {
-  skip_if(os_is_wsl())
 
   stan_file <- write_stan_file("parameters { real x; } model { x ~ std_normal(); }")
   mod <- cmdstan_model(stan_file, compile_model_methods = TRUE, force_recompile = TRUE)
@@ -272,7 +249,6 @@ test_that("Model methods can be initialised for models with no data", {
 })
 
 test_that("Variable skeleton returns correct dimensions for matrices", {
-  skip_if(os_is_wsl())
 
   stan_file <- write_stan_file("
   data {