Use saved log_lik for spatial models

3 files changed, 35 insertions, 27 deletions

diff --git a/R/inputs.R b/R/inputs.R
index a3845d6..cb0a555 100644
--- a/R/inputs.R
+++ b/R/inputs.R
@@ -3,7 +3,7 @@
   model_code <- name_to_modelcode(name)
   y_data <- get_stan_data(response)
 
-  pars <- get_pars(submodels, name %in% c("distsamp"))
+  pars <- get_pars(submodels, name)
   submodels <- unname(submodels@submodels)
   types <- sapply(submodels, function(x) x@type)
   submodel_data <- lapply(submodels, get_stan_data)
@@ -39,12 +39,17 @@ add_placeholder_priors <- function(submodel_data, types){
 
 setGeneric("get_pars", function(object, ...) standardGeneric("get_pars"))
 
-setMethod("get_pars", "ubmsSubmodelList", function(object, keep_loglik=FALSE, ...){
+setMethod("get_pars", "ubmsSubmodelList", function(object, model_name, ...){
   #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(keep_loglik) out <- c(out, "log_lik")
   out
 })
diff --git a/R/loglik.R b/R/loglik.R
index ef4dc71..c6ab712 100644
--- a/R/loglik.R
+++ b/R/loglik.R
@@ -128,6 +128,10 @@ setGeneric("extract_log_lik",
 
 setMethod("extract_log_lik", "ubmsFit",
           function(object, parameter_name = "log_lik", merge_chains=TRUE){
+  #if log_lik was saved as an output parameter, just return that
+  if(loglik_saved(object)){
+    return(loo::extract_log_lik(object@stanfit, parameter_name, merge_chains))
+  }
   inps <- rebuild_inputs(object)
   ll <- get_loglik(object, inps, merge_chains)
   if(!merge_chains){
@@ -138,8 +142,6 @@ setMethod("extract_log_lik", "ubmsFit",
   ll
 })
 
-#' @include distsamp.R
-setMethod("extract_log_lik", "ubmsFitDistsamp",
-          function(object, parameter_name = "log_lik", merge_chains=TRUE){
-  loo::extract_log_lik(object@stanfit, parameter_name, merge_chains)
-})
+loglik_saved <- function(fit){
+  "log_lik" %in% fit@stanfit@sim$pars_oi
+}
diff --git a/vignettes/spatial-models.Rmd b/vignettes/spatial-models.Rmd
index 8e6a780..63cc2f7 100644
--- a/vignettes/spatial-models.Rmd
+++ b/vignettes/spatial-models.Rmd
@@ -250,21 +250,22 @@ fit_ubms
 ## 
 ## Occupancy (logit-scale):
 ##             Estimate     SD    2.5%  97.5% n_eff Rhat
-## (Intercept)   -1.038 0.2286 -1.5060 -0.601  1269 1.00
-## habCov12      -0.212 0.3231 -0.8503  0.427  2404 1.00
-## habCov13      -0.258 0.3273 -0.9035  0.374  2058 1.00
-## habCov2        1.080 0.1708  0.7683  1.439  1429 1.00
-## RSR [tau]      0.155 0.0589  0.0704  0.293   253 1.01
+## (Intercept)   -1.054 0.2324 -1.5121 -0.611  1177 1.01
+## habCov12      -0.216 0.3201 -0.8831  0.383  2499 1.00
+## habCov13      -0.261 0.3223 -0.8950  0.356  2262 1.00
+## habCov2        1.089 0.1801  0.7416  1.461  1295 1.01
+## RSR [tau]      0.144 0.0579  0.0675  0.283   264 1.02
 ## 
 ## Detection (logit-scale):
 ##             Estimate    SD   2.5% 97.5% n_eff Rhat
-## (Intercept)    1.727 0.370  1.004 2.481  1610    1
-## detCov1        0.113 0.210 -0.310 0.533  4962    1
-## detCov21       1.185 0.555  0.122 2.320  2189    1
-## detCov22       1.372 0.557  0.356 2.496  2455    1
-## detCov23      -0.242 0.484 -1.197 0.719  2013    1
+## (Intercept)    1.734 0.371  1.035 2.488  1756    1
+## detCov1        0.108 0.203 -0.291 0.497  4281    1
+## detCov21       1.183 0.556  0.102 2.270  2319    1
+## detCov22       1.379 0.561  0.307 2.533  2471    1
+## detCov23      -0.245 0.486 -1.204 0.704  2218    1
 ## 
-## LOOIC: 668.02
+## LOOIC: 668.475
+## Runtime: 88.788 sec
 ```
 
 The precision term (tau, $\tau$) associated with the spatial random effect is the final row of the occupancy model estimates.
@@ -291,12 +292,12 @@ head(psi)
 
 ```
 ##   Predicted         SD       2.5%     97.5%
-## 1 0.1441475 0.04413909 0.07266996 0.2439409
-## 2 0.2360008 0.08352206 0.10501601 0.4271224
-## 3 0.3578073 0.11260628 0.17170766 0.6029410
-## 4 0.4527547 0.12485601 0.23020833 0.7192525
-## 5 0.4522272 0.11820159 0.23762047 0.6955722
-## 6 0.3847294 0.11279317 0.18753530 0.6292943
+## 1 0.1441804 0.04624230 0.07020541 0.2517342
+## 2 0.2390063 0.08739192 0.10297331 0.4388321
+## 3 0.3640516 0.11781411 0.16512449 0.6239752
+## 4 0.4593940 0.13070577 0.22866630 0.7265760
+## 5 0.4570641 0.12368794 0.23700396 0.7078097
+## 6 0.3865731 0.11876489 0.18351350 0.6434938
 ```
 
 ## Plotting
@@ -349,8 +350,8 @@ round(modSel(fl),2)
 
 ```
 ##                   elpd nparam elpd_diff se_diff weight
-## fit_ubms       -334.01  44.80      0.00     0.0      1
-## fit_nonspatial -380.54   8.43    -46.53     8.1      0
+## fit_ubms       -334.24  46.25      0.00    0.00      1
+## fit_nonspatial -380.46   8.36    -46.22    8.26      0
 ```
 
 A measure of predictive accuracy (`elpd`) is shown for each model, and the model with the highest predictive accuracy is listed first.
@@ -376,7 +377,7 @@ names <- list(
 fit_stocc <- spatial.occupancy(
   detection.model = ~ detCov1 + detCov2,
   occupancy.model = ~ habCov1 + habCov2,
-  spatial.model = list(model="rsr", threshold=1, moran.cut=0.1*nrow(habData)),
+  spatial.model = list(model="rsr", threshold=1.01, moran.cut=0.1*nrow(habData)),
   so.data = make.so.data(visitData, habData, names),
   prior = list(a.tau=0.5, b.tau=0.00005, Q.b=0.1, Q.g=0.1),
   control = list(burnin=1000/5, iter=4000/5, thin=5)