improved nmixture posterior generation

nicholasjclark · Feb 13, 2024 · 5a1a61e · 5a1a61e
1 parent ee10b43
commit 5a1a61e
Show file tree

Hide file tree

Showing 5 changed files with 21 additions and 26 deletions.
diff --git a/R/add_nmixture.R b/R/add_nmixture.R
@@ -626,7 +626,7 @@ add_nmix_posterior = function(model_output,
   latentypreds_vec[whichobs] <- condpreds_vec
   latentypreds <- matrix(latentypreds_vec, nrow = NROW(betas))
 
-  # TEST ####
+  # Update parameter names to match expected order
   expand.grid(time = 1:model_output@sim$dims_oi$trend[1],
               series = 1:model_output@sim$dims_oi$trend[2]) %>%
     dplyr::arrange(series, time) %>%
@@ -637,13 +637,9 @@ add_nmix_posterior = function(model_output,
                                 ']')) %>%
     dplyr::pull(name) -> parnames
 
-  # old_parnames <- dimnames(trend)[[2]]
-  # latentypreds <- latentypreds[, order(match(old_parnames,parnames))]
-  parnames <- gsub('trend', 'latent_ypred', parnames)
-
   # Add latent_ypreds to the posterior samples
   model_output <- add_samples(model_output = model_output,
-                        names = parnames,
+                        names = gsub('trend', 'latent_ypred', parnames),
                         samples = latentypreds,
                         nsamples = NROW(latentypreds) / nchains,
                         nchains = nchains,
@@ -670,8 +666,7 @@ add_nmix_posterior = function(model_output,
                        prob = as.vector(detprob))
   ypreds <- matrix(ypreds_vec, nrow = NROW(betas))
   model_output <- add_samples(model_output = model_output,
-                        names = gsub('trend', 'ypred',
-                                     dimnames(trend)[[2]]),
+                        names = gsub('trend', 'ypred', parnames),
                         samples = ypreds,
                         nsamples = NROW(ypreds) / nchains,
                         nchains = nchains,
@@ -683,8 +678,7 @@ add_nmix_posterior = function(model_output,
   mus_vec <- as.vector(detprob) * latentypreds_vec
   mus <- matrix(mus_vec, nrow = NROW(betas))
   model_output <- add_samples(model_output = model_output,
-                        names = gsub('trend', 'mus',
-                                     dimnames(trend)[[2]]),
+                        names = gsub('trend', 'mus', parnames),
                         samples = mus,
                         nsamples = NROW(mus) / nchains,
                         nchains = nchains,

diff --git a/docs/articles/nmixtures.html b/docs/articles/nmixtures.html
diff --git a/src/mvgam.dll b/src/mvgam.dll
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
diff --git a/vignettes/nmixtures.Rmd b/vignettes/nmixtures.Rmd
@@ -136,7 +136,7 @@ head(testdat, 12)
 
 ### Setting up the `trend_map`
 
-Finally, we need to set up the `trend_map` object. This is crucial for allowing multiple observations to be linked to the same latent process model (see more information about this argument in the [Shared latent states vignette](https://nicholasjclark.github.io/mvgam/articles/shared_states.html){target="_blank"}. In this case, the mapping operates by species and site to state that each set of replicate observations from the same time point should all share the exact same latent abundance model:
+Finally, we need to set up the `trend_map` object. This is crucial for allowing multiple observations to be linked to the same latent process model (see more information about this argument in the [Shared latent states vignette](https://nicholasjclark.github.io/mvgam/articles/shared_states.html){target="_blank"}). In this case, the mapping operates by species and site to state that each set of replicate observations from the same time point should all share the exact same latent abundance model:
 ```{r}
 testdat %>%
   # each unique combination of site*species is a separate process
@@ -208,7 +208,7 @@ The posterior summary of this model shows that it has converged nicely
 summary(mod)
 ```
 
-`loo()` functionality works just as it does for all `mvgam` models to aid in model comparison / selection
+`loo()` functionality works just as it does for all `mvgam` models to aid in model comparison / selection (though note that Pareto K values often give warnings for mixture models so these may not be too helpful)
 ```{r}
 loo(mod)
 ```
@@ -218,7 +218,7 @@ Plot the estimated smooths of time from each species' latent abundance process (
 plot(mod, type = 'smooths', trend_effects = TRUE)
 ```
 
-`marginaleffects` support allows for more useful prediction-based interrogations on different scales. Objects that use family `nmix()` have a few additional prediction scales that can be used (i.e. `link`, `response`, `detection` or `latent_N`). For example, here are the estimated detection probabilities per species, which shows that the model has over-estimated detection probability for species 2 (originally simulated to be 0.45):
+`marginaleffects` support allows for more useful prediction-based interrogations on different scales. Objects that use family `nmix()` have a few additional prediction scales that can be used (i.e. `link`, `response`, `detection` or `latent_N`). For example, here are the estimated detection probabilities per species, which show that the model has done a nice job of estimating these parameters:
 ```{r}
 plot_predictions(mod, condition = 'species',
                  type = 'detection') +