Issues with calculating WAIC using waicAbund()

[Duston Duffie]

Hello Jeff and the spOccupancy group,

A colleague and I are using msNMix() function in the spAbundance package to run multi-species n-mixture models on an amphibian community dataset in a stacked data format for each survey year. The goal of these analyses is to evaluate impacts of different vegetation communities on the amphibian community, so we are interested in the community-level responses to various vegetation covariates. When running the waicAbund() function on a null model one of the species is returning Inf values in the elpd and WAIC columns (see table below). We are hoping for some advice on this issue and whether we should change our priors or attempt to use a different model to analyze these data.

Additional information:

The data is from herpetofauna capture data from 6 sites over a 4-year period, and we are modeling data from 11 amphibian species. We selected a negative binominal distribution after fitting simple models with single covariate in both the abundance and detection formula with a negative binomial and Poisson distribution. We compared the fit of these models using the ppcAbund() function, but we did not initially run WAIC test on these models. After realizing the issues with WAIC, we did attempt to calculate WAIC for these simplified models, and they similarly calculated Inf values.

The species with Inf values returned by WAIC does not appear to have converged for the species-level abundance or detection beta coefficients. I think one potential issue in this dataset may be that during the second trapping year, we captured many more individuals of the species with convergence issues than in the previous year or in the following two years, so capture numbers/ abundance is greatly inflated during that one year (0-10 captures at sites in the first year, 100-600 captures at sites in the second year, 0-50 captures at sites in the third year, and 0-10 captures at sites in the fourth year).

R code used for running null model:

n.batch = 5000

batch.length = 25

n.thin = 50

n.chains = 3

n.samples <- n.batch * batch.length

n.burn <- 0.6*(n.batch*batch.length)

n.sp <- dim(data_list$y)[1]

ms.inits <- list(alpha.comm = 0, # community-level detection coefficients

                   beta.comm = 0, # community-level abundance coefficients

                   beta = 0, # species-level detection coefficients

                   alpha = 0, # species-level abundance coefficients

                   tau.sq.beta = 1, # community-level abundance variance parameters

                   kappa = 1, # species-level negative binomial overdispersion parameters

                   tau.sq.alpha = 1, # community-level detection variance parameters

                   sigma.sq.mu = 0.5, # random effect variances for abundance

                   N = apply(data_list$y, c(1, 2), max, na.rm = TRUE)) #latent abundance variables for all species

 

ms.priors <- list(beta.comm.normal = list(mean = 0, var = 100), #normal prior on the community-level abundance mean effects

                    alpha.comm.normal = list(mean = 0, var = 2.72), #normal prior on the community-level detection mean effects

                    tau.sq.beta.ig = list(a = 0.1, b = 0.1), #inverse-Gamma prior on the community-level abundance variance parameters

                    tau.sq.alpha.ig = list(a = 0.1, b = 0.1), #inverse-Gamma prior on the community-level detection variance parameters

                    sigma.sq.mu.ig = list(a = 0.1, b = 0.1), #inverse-Gamma prior on the abundance random effect variances

                    kappa.unif = list(a = 0, b = 100)) # uniform prior on the species-specific overdispersion parameters

# Specify initial tuning values

ms.tuning <- list(beta = 0.3, alpha = 0.3, beta.star = 0.5, alpha.star = 0.5, kappa = 0.5)

 

# null model

null_mod <- msNMix(abund.formula = ~ 1,

                              det.formula = ~ 1 ,

                              data = data_list,

                              inits = ms.inits,

                              tuning = ms.tuning,

                              family = 'NB',

                              priors = ms.priors,

                              n.omp.threads = 7,

                              verbose = TRUE,

                              n.report = 400,

                              n.burn = n.burn,

                              n.thin = n.thin,

                              n.chains = n.chains,

                              n.batch = n.batch,

                              batch.length = batch.length)

 

waicAbund(null_mod, by.species = TRUE)

elpd        pD                 WAIC
-90.3544 2.954032 186.6168
-41.5423 1.155127 85.39485
-75.6755 1.920855 155.1927
-104.83         3.767612 217.1948
Inf                 1.197839 Inf
-348.324 2.062787 700.7727
-18.0473 2.653417 41.40146
-180.792 1.924675 365.4332
-8.68403 1.078886 19.52584
-8.67433 1.034205 19.41708
-8.66752 1.055561 19.44617

Any advice would be greatly appreciated.

Thank you,

Duston 

[Jeffrey Doser]

Hi Dustin, 

Apologies for the delay, I’ve been in and out of the office over the last couple weeks. It sounds like the problem here is likely related to the use of N-mixture models rather than a problem with WAIC itself. In your case, for the species that isn’t converging, I’m guessing based on your description that the estimated abundance values are extremely large as a result of the large variation over time that’s not going to be accounted for in the null model. Since your null model does not account for any variation in detection probability across the visits or abundance across sites/years, it’s probably going to estimate detection probability to be very small and subsequently generate massive abundance estimates. When estimates of abundance are estimated to be extremely large, N-mixture models can have difficulty converging, which is likely what’s going on in your case. When that happens, the way WAIC is calculated will result in numerical challenges and lead to a value of Inf. So, I don’t think there’s much you can do in this case to fix the WAIC value without putting an extremely informative prior on the beta coefficients, which I would not recommend doing without any justification beyond trying to get it to converge. In general, based on the small number of sites you have and the large variation in abundance values, an N-mixture model may have a lot of trouble converging. In that case, a relative abundance model (i.e. a GLMM akin to the one fit by msAbund) might be a good way to go. 

Jeff

Jeffrey W. Doser, Ph.D.

Assistant Professor

Department of Forestry and Environmental Resources

North Carolina State University

Statistical Ecology and Forest Science Lab

Pronouns: he/him/his

--
You received this message because you are subscribed to the Google Groups "spOccupancy and spAbundance users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spocc-spabund-u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/spocc-spabund-users/a86be893-2fef-4446-b7d1-5f3e76b4f992n%40googlegroups.com.

[Duston Duffie]

Hello Jeff,

Thank you for your response. I think the relative abundance models may be what we need to use for this dataset. However, one suggestion from my advisor is to remove the species that is not converging from these analyses given that the increase in abundance for that species is likely due to a factor that is not the main focus of this project. This species is an explosive breeder, so the increase in abundance may be more of a response to heavy rainfall in that year. I am relatively new to multispecies/ community models. Given that the focus of this project is the community-level response, would it be a valid approach to remove the most commonly captured species in this community? I am assuming removing this species would impact any community-level estimates. 

Thanks,

Duston