Add option to save log_lik parameter

23 files changed, 111 insertions, 29 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 03dd646..7cd21b7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ubms
 Version: 1.1.9005
-Date: 2022-09-19
+Date: 2022-09-20
 Title: Bayesian Models for Data from Unmarked Animals using 'Stan'
 Authors@R: person("Ken", "Kellner", email="contact@kenkellner.com", 
                   role=c("cre","aut"))
@@ -34,7 +34,7 @@ License: GPL (>=3)
 URL: https://kenkellner.com/ubms/
 BugReports: https://github.com/kenkellner/ubms/issues
 Encoding: UTF-8
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.0
 Biarch: true
 LinkingTo: 
     BH (>= 1.66.0),
diff --git a/R/colext.R b/R/colext.R
index 49fce3a..40cea5b 100644
--- a/R/colext.R
+++ b/R/colext.R
@@ -26,6 +26,9 @@
 #' @param prior_coef_det Prior distribution for the regression coefficients of
 #'  the detection model
 #' @param prior_sigma Prior distribution on random effect standard deviations
+#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values
+#'  in the output. This can greatly increase the size of the model. If
+#'  \code{FALSE}, the values are calculated post-hoc from the posteriors
 #' @param ... Arguments passed to the \code{\link{stan}} call, such as
 #'  number of chains \code{chains} or iterations \code{iter}
 #'
@@ -60,6 +63,7 @@ stan_colext <- function(psiformula = ~1,
                         prior_intercept_det = logistic(0, 1),
                         prior_coef_det = logistic(0, 1),
                         prior_sigma = gamma(1, 1),
+                        log_lik = FALSE,
                         ...){
 
   umf <- process_umf(data)
@@ -80,7 +84,7 @@ stan_colext <- function(psiformula = ~1,
                       prior_intercept_det, prior_coef_det, prior_sigma)
   submodels <- ubmsSubmodelList(state, col, ext, det)
 
-  ubmsFit("colext", match.call(), data, response, submodels, ...)
+  ubmsFit("colext", match.call(), data, response, submodels, log_lik, ...)
 }
 
 
diff --git a/R/fit.R b/R/fit.R
index 5abfc49..dedf6e4 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -19,13 +19,13 @@ setClass("ubmsFitOccu", contains = "ubmsFit")
 # Child class for abundance/N-mixture type models
 setClass("ubmsFitAbun", contains = "ubmsFit")
 
-ubmsFit <- function(model, call, data, response, submodels, ...){
+ubmsFit <- function(model, call, data, response, submodels, log_lik=TRUE, ...){
   #Find missing
   response <- update_missing(response, submodels)
   submodels <- update_missing(submodels, response)
 
   #Fit model
-  fit <- fit_model(model, response, submodels, ...)
+  fit <- fit_model(model, response, submodels, log_lik, ...)
   # Exit early if just returning Stan inputs
   if(check_return_inputs(...)) return(fit)
 
@@ -75,9 +75,9 @@ fit_class <- function(mod){
 
 #Fit stan model
 #' @include inputs.R
-fit_model <- function(name, response, submodels, ...){
+fit_model <- function(name, response, submodels, log_lik, ...){
   model <- name_to_stanmodel(name, submodels)
-  inp <- build_stan_inputs(name, response, submodels)
+  inp <- build_stan_inputs(name, response, submodels, log_lik)
   # Should just Stan inputs be returned?
   if(check_return_inputs(...)){
     inp$submodels <- submodels
diff --git a/R/inputs.R b/R/inputs.R
index cb0a555..f3c183b 100644
--- a/R/inputs.R
+++ b/R/inputs.R
@@ -1,9 +1,9 @@
- build_stan_inputs <- function(name, response, submodels, ...){
+ build_stan_inputs <- function(name, response, submodels, log_lik, ...){
 
   model_code <- name_to_modelcode(name)
   y_data <- get_stan_data(response)
 
-  pars <- get_pars(submodels, name)
+  pars <- get_pars(submodels, name, log_lik)
   submodels <- unname(submodels@submodels)
   types <- sapply(submodels, function(x) x@type)
   submodel_data <- lapply(submodels, get_stan_data)
@@ -39,18 +39,17 @@ add_placeholder_priors <- function(submodel_data, types){
 
 setGeneric("get_pars", function(object, ...) standardGeneric("get_pars"))
 
-setMethod("get_pars", "ubmsSubmodelList", function(object, model_name, ...){
+setMethod("get_pars", "ubmsSubmodelList", function(object, model_name, log_lik, ...){
   #Remove placeholder submodels - we don't want to save those parameters
   submodels <- object@submodels
   submodels <- submodels[!sapply(submodels, is_placeholder)]
   submodels <- unname(submodels)
   out <- unlist(lapply(submodels, get_pars))
 
-  keep_loglik <- FALSE
-  if(model_name == "distsamp") keep_loglik <- TRUE
-  if(any(sapply(submodels, has_spatial))) keep_loglik <- TRUE
+  if(model_name == "distsamp") log_lik <- TRUE
+  if(any(sapply(submodels, has_spatial))) log_lik <- TRUE
 
-  if(keep_loglik) out <- c(out, "log_lik")
+  if(log_lik) out <- c(out, "log_lik")
   out
 })
 
diff --git a/R/kfold.R b/R/kfold.R
index 9402bdc..3dc0500 100644
--- a/R/kfold.R
+++ b/R/kfold.R
@@ -26,6 +26,9 @@ setMethod("kfold", "ubmsFit", function(x, K=10, folds=NULL, quiet=FALSE, ...){
   if(has_spatial(x)){
     stop("kfold does not work with spatial models", call.=FALSE)
   }
+  if(inherits(x, "ubmsFitDistsamp")){
+    stop("kfold method not yet supported for stan_distsamp models", call.=FALSE)
+  }
 
   op <- pbapply::pboptions()
   if(quiet) pbapply::pboptions(type = "none")
diff --git a/R/multinomPois.R b/R/multinomPois.R
index 9b05ef6..88618d3 100644
--- a/R/multinomPois.R
+++ b/R/multinomPois.R
@@ -15,6 +15,9 @@
 #' @param prior_coef_det Prior distribution for the regression coefficients of
 #'  the detection model
 #' @param prior_sigma Prior distribution on random effect standard deviations
+#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values
+#'  in the output. This can greatly increase the size of the model. If
+#'  \code{FALSE}, the values are calculated post-hoc from the posteriors
 #' @param ... Arguments passed to the \code{\link{stan}} call, such as
 #'  number of chains \code{chains} or iterations \code{iter}
 #'
@@ -39,6 +42,7 @@ stan_multinomPois <- function(formula,
                               prior_intercept_det = logistic(0, 1),
                               prior_coef_det = logistic(0, 1),
                               prior_sigma = gamma(1, 1),
+                              log_lik = FALSE,
                               ...){
 
   forms <- split_formula(formula)
@@ -62,7 +66,7 @@ stan_multinomPois <- function(formula,
                       prior_intercept_det, prior_coef_det, prior_sigma)
   submodels <- ubmsSubmodelList(state, det)
 
-  ubmsFit("multinomPois", match.call(), data, response, submodels, ...)
+  ubmsFit("multinomPois", match.call(), data, response, submodels, log_lik, ...)
 }
 
 get_pifun_type <- function(umf){
diff --git a/R/occu.R b/R/occu.R
index 4846987..8f1e47b 100644
--- a/R/occu.R
+++ b/R/occu.R
@@ -15,6 +15,9 @@
 #' @param prior_coef_det Prior distribution for the regression coefficients of
 #'  the detection model
 #' @param prior_sigma Prior distribution on random effect standard deviations
+#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values
+#'  in the output. This can greatly increase the size of the model. If
+#'  \code{FALSE}, the values are calculated post-hoc from the posteriors
 #' @param ... Arguments passed to the \code{\link{stan}} call, such as
 #'  number of chains \code{chains} or iterations \code{iter}
 #'
@@ -46,6 +49,7 @@ stan_occu <- function(formula,
                       prior_intercept_det = logistic(0, 1),
                       prior_coef_det = logistic(0, 1),
                       prior_sigma = gamma(1, 1),
+                      log_lik = FALSE,
                       ...){
 
   forms <- split_formula(formula)
@@ -68,7 +72,7 @@ stan_occu <- function(formula,
                       prior_intercept_det, prior_coef_det, prior_sigma)
   submodels <- ubmsSubmodelList(state, det)
 
-  ubmsFit("occu", match.call(), data, response, submodels, ...)
+  ubmsFit("occu", match.call(), data, response, submodels, log_lik, ...)
 }
 
 
diff --git a/R/occuRN.R b/R/occuRN.R
index e16bb87..77c8e47 100644
--- a/R/occuRN.R
+++ b/R/occuRN.R
@@ -19,6 +19,9 @@
 #' @param prior_coef_det Prior distribution for the regression coefficients of
 #'  the detection model
 #' @param prior_sigma Prior distribution on random effect standard deviations
+#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values
+#'  in the output. This can greatly increase the size of the model. If
+#'  \code{FALSE}, the values are calculated post-hoc from the posteriors
 #' @param ... Arguments passed to the \code{\link{stan}} call, such as
 #'  number of chains \code{chains} or iterations \code{iter}
 #'
@@ -47,6 +50,7 @@ stan_occuRN <- function(formula,
                         prior_intercept_det = logistic(0, 1),
                         prior_coef_det = logistic(0, 1),
                         prior_sigma = gamma(1, 1),
+                        log_lik = FALSE,
                         ...){
 
   forms <- split_formula(formula)
@@ -69,7 +73,7 @@ stan_occuRN <- function(formula,
                       prior_intercept_det, prior_coef_det, prior_sigma)
   submodels <- ubmsSubmodelList(state, det)
 
-  ubmsFit("occuRN", match.call(), data, response, submodels, ...)
+  ubmsFit("occuRN", match.call(), data, response, submodels, log_lik, ...)
 }
 
 
diff --git a/R/occuTTD.R b/R/occuTTD.R
index 605bec3..83aef23 100644
--- a/R/occuTTD.R
+++ b/R/occuTTD.R
@@ -30,6 +30,9 @@
 #' @param prior_intercept_shape Prior distribution for the intercept of the
 #'  shape parameter (i.e., log(shape)) for Weibull TTD models
 #' @param prior_sigma Prior distribution on random effect standard deviations
+#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values
+#'  in the output. This can greatly increase the size of the model. If
+#'  \code{FALSE}, the values are calculated post-hoc from the posteriors
 #' @param ... Arguments passed to the \code{\link{stan}} call, such as
 #'  number of chains \code{chains} or iterations \code{iter}
 #'
@@ -88,6 +91,7 @@ stan_occuTTD <- function(psiformula=~1,
                          prior_coef_det = normal(0, 2.5),
                          prior_intercept_shape = normal(0,2.5),
                          prior_sigma = gamma(1, 1),
+                         log_lik = FALSE,
                          ...){
 
   if(data@numPrimary > 1) stop("Dynamic models not yet supported", call.=FALSE)
@@ -119,7 +123,7 @@ stan_occuTTD <- function(psiformula=~1,
 
   submodels <- ubmsSubmodelList(state, det, shape)
 
-  ubmsFit("occuTTD", match.call(), data, response, submodels, ...)
+  ubmsFit("occuTTD", match.call(), data, response, submodels, log_lik, ...)
 }
 
 
diff --git a/R/pcount.R b/R/pcount.R
index 219703f..c44abfa 100644
--- a/R/pcount.R
+++ b/R/pcount.R
@@ -19,6 +19,9 @@
 #' @param prior_coef_det Prior distribution for the regression coefficients of
 #'  the detection model
 #' @param prior_sigma Prior distribution on random effect standard deviations
+#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values
+#'  in the output. This can greatly increase the size of the model. If
+#'  \code{FALSE}, the values are calculated post-hoc from the posteriors
 #' @param ... Arguments passed to the \code{\link{stan}} call, such as
 #'  number of chains \code{chains} or iterations \code{iter}
 #'
@@ -47,6 +50,7 @@ stan_pcount <- function(formula,
                         prior_intercept_det = logistic(0, 1),
                         prior_coef_det = logistic(0, 1),
                         prior_sigma = gamma(1, 1),
+                        log_lik = FALSE,
                         ...){
 
   forms <- split_formula(formula)
@@ -69,7 +73,7 @@ stan_pcount <- function(formula,
                       prior_intercept_det, prior_coef_det, prior_sigma)
   submodels <- ubmsSubmodelList(state, det)
 
-  ubmsFit("pcount", match.call(), data, response, submodels, ...)
+  ubmsFit("pcount", match.call(), data, response, submodels, log_lik, ...)
 }
 
 #Output object-----------------------------------------------------------------
diff --git a/R/ubmsFit-methods.R b/R/ubmsFit-methods.R
index ba42790..2f0bac8 100644
--- a/R/ubmsFit-methods.R
+++ b/R/ubmsFit-methods.R
@@ -288,8 +288,8 @@ setMethod("get_stancode", "ubmsFit", function(object, ...){
 
 # Re-create inputs
 rebuild_inputs <- function(object){
-  inps <- build_stan_inputs(object@stanfit@model_name,
-                            object@response, object@submodels)
+  inps <- build_stan_inputs(object@stanfit@model_name, object@response,
+                            object@submodels, loglik_saved(object))
   inps$submodels <- object@submodels
   inps
 }
diff --git a/man/extract_log_lik.Rd b/man/extract_log_lik.Rd
index 371b894..6a2e861 100644
--- a/man/extract_log_lik.Rd
+++ b/man/extract_log_lik.Rd
@@ -3,7 +3,6 @@
 \name{extract_log_lik}
 \alias{extract_log_lik}
 \alias{extract_log_lik,ubmsFit-method}
-\alias{extract_log_lik,ubmsFitOccuTTD-method}
 \alias{extract_log_lik,ubmsFitDistsamp-method}
 \title{Extract Pointwise Log-likelihood From Model}
 \usage{
diff --git a/man/kfold-ubmsFit-method.Rd b/man/kfold-ubmsFit-method.Rd
index 73f4267..c5cd1e7 100644
--- a/man/kfold-ubmsFit-method.Rd
+++ b/man/kfold-ubmsFit-method.Rd
@@ -21,5 +21,8 @@
 An object of class \code{elpd_generic} that is compatible with \code{loo::loo_compare}
 }
 \description{
-K-fold Cross-validation of a ubmsFit Model
+Randomly partition data into K subsets of equal size (by site). Re-fit the model
+K times, each time leaving out one of the subsets. Calculate the log-likelihood
+for each of the sites that was left out. This function is an alternative
+to \code{loo} (leave-one-out cross validation).
 }
diff --git a/man/stan_colext.Rd b/man/stan_colext.Rd
index 8ce5688..9493637 100644
--- a/man/stan_colext.Rd
+++ b/man/stan_colext.Rd
@@ -19,6 +19,7 @@ stan_colext(
   prior_intercept_det = logistic(0, 1),
   prior_coef_det = logistic(0, 1),
   prior_sigma = gamma(1, 1),
+  log_lik = FALSE,
   ...
 )
 }
@@ -60,6 +61,10 @@ the detection model}
 
 \item{prior_sigma}{Prior distribution on random effect standard deviations}
 
+\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values
+in the output. This can greatly increase the size of the model. If
+\code{FALSE}, the values are calculated post-hoc from the posteriors}
+
 \item{...}{Arguments passed to the \code{\link{stan}} call, such as
 number of chains \code{chains} or iterations \code{iter}}
 }
diff --git a/man/stan_multinomPois.Rd b/man/stan_multinomPois.Rd
index ff34467..4237be7 100644
--- a/man/stan_multinomPois.Rd
+++ b/man/stan_multinomPois.Rd
@@ -12,6 +12,7 @@ stan_multinomPois(
   prior_intercept_det = logistic(0, 1),
   prior_coef_det = logistic(0, 1),
   prior_sigma = gamma(1, 1),
+  log_lik = FALSE,
   ...
 )
 }
@@ -35,6 +36,10 @@ the detection model}
 
 \item{prior_sigma}{Prior distribution on random effect standard deviations}
 
+\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values
+in the output. This can greatly increase the size of the model. If
+\code{FALSE}, the values are calculated post-hoc from the posteriors}
+
 \item{...}{Arguments passed to the \code{\link{stan}} call, such as
 number of chains \code{chains} or iterations \code{iter}}
 }
diff --git a/man/stan_occu.Rd b/man/stan_occu.Rd
index 512bc64..86bc7fb 100644
--- a/man/stan_occu.Rd
+++ b/man/stan_occu.Rd
@@ -12,6 +12,7 @@ stan_occu(
   prior_intercept_det = logistic(0, 1),
   prior_coef_det = logistic(0, 1),
   prior_sigma = gamma(1, 1),
+  log_lik = FALSE,
   ...
 )
 }
@@ -35,6 +36,10 @@ the detection model}
 
 \item{prior_sigma}{Prior distribution on random effect standard deviations}
 
+\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values
+in the output. This can greatly increase the size of the model. If
+\code{FALSE}, the values are calculated post-hoc from the posteriors}
+
 \item{...}{Arguments passed to the \code{\link{stan}} call, such as
 number of chains \code{chains} or iterations \code{iter}}
 }
diff --git a/man/stan_occuRN.Rd b/man/stan_occuRN.Rd
index fcf0f02..bdc954e 100644
--- a/man/stan_occuRN.Rd
+++ b/man/stan_occuRN.Rd
@@ -13,6 +13,7 @@ stan_occuRN(
   prior_intercept_det = logistic(0, 1),
   prior_coef_det = logistic(0, 1),
   prior_sigma = gamma(1, 1),
+  log_lik = FALSE,
   ...
 )
 }
@@ -40,6 +41,10 @@ the detection model}
 
 \item{prior_sigma}{Prior distribution on random effect standard deviations}
 
+\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values
+in the output. This can greatly increase the size of the model. If
+\code{FALSE}, the values are calculated post-hoc from the posteriors}
+
 \item{...}{Arguments passed to the \code{\link{stan}} call, such as
 number of chains \code{chains} or iterations \code{iter}}
 }
diff --git a/man/stan_occuTTD.Rd b/man/stan_occuTTD.Rd
index 3583848..0efd163 100644
--- a/man/stan_occuTTD.Rd
+++ b/man/stan_occuTTD.Rd
@@ -18,6 +18,7 @@ stan_occuTTD(
   prior_coef_det = normal(0, 2.5),
   prior_intercept_shape = normal(0, 2.5),
   prior_sigma = gamma(1, 1),
+  log_lik = FALSE,
   ...
 )
 }
@@ -61,6 +62,10 @@ shape parameter (i.e., log(shape)) for Weibull TTD models}
 
 \item{prior_sigma}{Prior distribution on random effect standard deviations}
 
+\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values
+in the output. This can greatly increase the size of the model. If
+\code{FALSE}, the values are calculated post-hoc from the posteriors}
+
 \item{...}{Arguments passed to the \code{\link{stan}} call, such as
 number of chains \code{chains} or iterations \code{iter}}
 }
diff --git a/man/stan_pcount.Rd b/man/stan_pcount.Rd
index 2f5a635..66c0bec 100644
--- a/man/stan_pcount.Rd
+++ b/man/stan_pcount.Rd
@@ -14,6 +14,7 @@ stan_pcount(
   prior_intercept_det = logistic(0, 1),
   prior_coef_det = logistic(0, 1),
   prior_sigma = gamma(1, 1),
+  log_lik = FALSE,
   ...
 )
 }
@@ -43,6 +44,10 @@ the detection model}
 
 \item{prior_sigma}{Prior distribution on random effect standard deviations}
 
+\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values
+in the output. This can greatly increase the size of the model. If
+\code{FALSE}, the values are calculated post-hoc from the posteriors}
+
 \item{...}{Arguments passed to the \code{\link{stan}} call, such as
 number of chains \code{chains} or iterations \code{iter}}
 }
diff --git a/tests/testthat/test_distsamp.R b/tests/testthat/test_distsamp.R
index b875be3..c26afa0 100644
--- a/tests/testthat/test_distsamp.R
+++ b/tests/testthat/test_distsamp.R
@@ -124,8 +124,8 @@ test_that("stan_distsamp handles NA values",{
 test_that("extract_log_lik method works",{
   ll <- extract_log_lik(fit_pt_hn)
   expect_is(ll, "matrix")
-  expect_equal(dim(ll), c(200/2 * 2, numSites(fit_pt_hn@data))) 
-  expect_between(sum(ll), -39000, -38000) 
+  expect_equal(dim(ll), c(200/2 * 2, numSites(fit_pt_hn@data)))
+  expect_between(sum(ll), -39000, -38000)
 })
 
 test_that("ubmsFitDistsamp gof method works",{
@@ -318,3 +318,8 @@ test_that("distsamp spatial works", {
   ps <- plot_spatial(fit_spat)
   expect_is(ps, "gg")
 })
+
+test_that("kfold errors when used on a stan_distsamp model", {
+  expect_error(kfold(fit_line_hn),
+               "kfold method not yet supported for stan_distsamp models")
+})
diff --git a/tests/testthat/test_fit.R b/tests/testthat/test_fit.R
index 6ff842d..bfd2497 100644
--- a/tests/testthat/test_fit.R
+++ b/tests/testthat/test_fit.R
@@ -23,7 +23,7 @@ resp <- ubmsResponse(umf@y,"binomial","binomial",max_primary=1)
 
 #Build a stanfit object
 set.seed(123)
-inp <- build_stan_inputs("occu", resp, sl)
+inp <- build_stan_inputs("occu", resp, sl, log_lik=FALSE)
 
 good_fit <- TRUE
 tryCatch({
@@ -66,7 +66,7 @@ test_that("remove_placeholders removes placeholder submodels from list",{
 
 test_that("fit_model builds model correctly",{
   ufit <- suppressWarnings(
-    fit_model("occu", resp, sl, chains=2, iter=20, refresh=0))
+    fit_model("occu", resp, sl, log_lik=FALSE, chains=2, iter=20, refresh=0))
   expect_true(inherits(ufit, "stanfit"))
   nms <- stanfit_names(sl)
   expect_equal(ufit@sim$fnames_oi[1:length(nms)], nms)
@@ -81,7 +81,7 @@ test_that("fit_model builds model correctly",{
 test_that("specific model name is shown in console output",{
   # e.g. 'occu' instead of 'single_season'
   out <- capture.output(ufit <- suppressWarnings(
-    fit_model("occu", resp, sl, chains=2, iter=20)))
+    fit_model("occu", resp, sl, log_lik=FALSE, chains=2, iter=20)))
   expect_true(any(grepl("occu", out)))
   expect_false(any(grepl("single_season", out)))
 })
diff --git a/tests/testthat/test_inputs.R b/tests/testthat/test_inputs.R
index 2a221a7..9bf5e9e 100644
--- a/tests/testthat/test_inputs.R
+++ b/tests/testthat/test_inputs.R
@@ -10,7 +10,7 @@ test_that("stan inputs are built correctly", {
   y <- matrix(c(1,0,0,1,1,1,0,0,1), nrow=3, byrow=T)
   resp <- ubmsResponse(y, "binomial", "P")
 
-  inp <- build_stan_inputs("occu", resp, sl)
+  inp <- build_stan_inputs("occu", resp, sl, log_lik=FALSE)
   expect_is(inp, "list")
   expect_equal(names(inp), c("stan_data", "pars"))
 
@@ -33,6 +33,10 @@ test_that("parameter list for stan is generated correctly",{
   sl <- ubmsSubmodelList(state, det)
   expect_equal(get_pars(det), "beta_det")
   expect_equal(get_pars(state), c("beta_state", "b_state", "sigma_state"))
+  expect_equal(get_pars(sl, "occu", log_lik=FALSE),
+               c("beta_state", "b_state", "sigma_state", "beta_det"))
+  expect_equal(get_pars(sl, "occu", log_lik=TRUE),
+               c("beta_state", "b_state", "sigma_state", "beta_det", "log_lik"))
 })
 
 test_that("get_stan_data pulls necessary info from response object",{
diff --git a/tests/testthat/test_occu.R b/tests/testthat/test_occu.R
index 87577ac..29b0c9f 100644
--- a/tests/testthat/test_occu.R
+++ b/tests/testthat/test_occu.R
@@ -88,6 +88,21 @@ test_that("extract_log_lik works when there are missing values and random effect
   expect_equal(dim(ll), c(50,9))
 })
 
+test_that("log_lik argument controls saving log_lik parameter", {
+  skip_on_cran()
+  set.seed(123)
+  fit <- suppressWarnings(stan_occu(~x2~x1, umf[1:10,], chains=2,
+                                  iter=100, refresh=0))
+  set.seed(123)
+  fit2 <- suppressWarnings(stan_occu(~x2~x1, umf[1:10,], chains=2,
+                                  iter=100, refresh=0, log_lik=TRUE))
+
+  expect_equal(fit@loo$estimates, fit2@loo$estimates)
+  expect_false("log_lik" %in% fit@stanfit@sim$pars_oi)
+  expect_true("log_lik" %in% fit2@stanfit@sim$pars_oi)
+
+})
+
 test_that("ubmsFitOccu gof method works",{
   set.seed(123)
   g <- gof(fit, draws=5, quiet=TRUE)