Model main effects and posterior comparisons

[Julien Piquet]

Hi all,

I am runnning an N-mixture model in spAbundance, which includes both continuous and categorical factors. As I am interested in reporting the effect of each factor on species density/detection, I would like to know if there is a way of obtaining model main effects or anything similar. Model summary provides level-specific estimates, which is not exactly what I need. Also, is there any way to run posterior pairwise comparisons among levels of the categorical predictors? Thank you in advance.

Best

Julien

[Bruce Wayne]

Hello Julien,

I have never responded in this forum before but long time reader.  I think you are looking for species level comparisons but my experience might lead you in the direction you need.  I used multi-species spOccupancy for community level effects (the effects were categorical).  In my case I extracted the `beta.comm.samples` for each level, formatted them with a column for .chain and .draw and used that input into the package `tidybayes`.  Tidybayes has a really cool function called `compare_levels` which generates posterior pairwise comparisons for all levels of your categorical predictor. Then you input those results into `ggplot`.

I can share some code if you need, you just have to direct it to extract the appropriate samples from your model output.

Alan

[Julien Piquet]

Hi Alan,

Thank you for your answer. I would definitely appreciate if you could share some code. However, I think your answer only solves one of my problems. Right now I am running the model with a categorical predictor that includes 7 land cover uses. When running the model sumary I obtain what frequentists call simple effects (I have no idea whether these receive the same name in Bayesian stats). That is, I get how much abundance increases or decreases when land cover uses change from the reference level (let's say urban) to the remaining levels (agricultural landscape, forests and so on). However, I have no idea of the the effect land cover use has as a factor (main effect, I have also seen it called omnibus effect). With conventional GLMs this can be obtained through the function Anova() but I have no idea about how to get it here. If I understood corrrectly, the approach you propose would not allow me to get these main effects, but rather to execute pairwise comparisons among land conver use levels. This is also great, as model summary does not inform on the differences among non-reference levels. It would just be perfect if I could also obtain the other thing. Anyway, thank you in advance.

Best

Julien

[Bruce Wayne]

Julien,

It's not clear what you mean between main effects and simple effects.  How you extract your results depends on how you parameterized your model.  For the example that I attached code for, my model is:  occupancy ~ use.  The use covariate has three categories.  The intercept is the reference level, so the actual means for the other two levels need to be added to the intercept.

It might help if you post your model designation for abundance to help us with a diagnosis. Again, my code is for occupancy but the structure should be the same from `spOccupancy` to `spAbundance`.

Thanks

Alan

[Julien Piquet]

Hi Alan,

First of all, thank you for providing example code. So, with main effects I mean the overall effect of a factor, which is not the same as the abundance estimate for the factor levels. In my case (again, I have a categorial predictor with multiple levels, see model below), I would like to know if the factor land cover use has a significant effect on species abundance (in addition to the specific abundance for each factor level). However, after reading a bit more about Bayesia outputs, I am getting the idea that such a procedure might not work for Bayesian models and could be exclusive of frequentist approaches. Anyway, here would be my model:

NMix(abund.formula = ~ scale(tmed) + landcover,
       det.formula = ~ scale(t50) + (1|observers),
       data = sp,
       n.batch = 1600,
       batch.length = 25,
       n.burn = 20000,
       thin = 20,
       n.chains = 3,
       n.report = 500,
       family = "Poisson",
       verbose = T)

I have 11 species in 7 islands, so putting the entire dataset would be a bit too much, but the data for one species are attached in case you would want to a look.

Best and again thank you,

P.S. for this example I just set n.batch to a standard value but in my actual analyses, this is higher (e.g. 4000). Also, I just specified a reduced number of det.covs but the original model would include other predictors.

Julien

[Jeffrey Doser]

Hi Julien,

I'll chime in here in regards to your interest in an overall effect of a factor via an ANOVA style output. I think it is technically possible to get something along the lines of what you're imagining by generating a variety of derived quantities from the model output (see e.g., the Bayesian section in this document related to stan). However, if the goal of getting an "overall effect" is to test the hypothesis that landcover is important, then perhaps the easiest thing to do (with the tools in spAbundance) would be to simply compare WAIC by species for a model with and without landcover. Alternatively, a more standard approach to Bayesian hypothesis testing would be to use what are called Bayes Factors, which should be possible to derive with output from spOccupancy/spAbundance, although I admittedly have not attempted to do that. While neither gives you the same output as the ANOVA style table, I think they would both give you an overall sense of the importance of landcover and then using Alan's nice code you could get the pairwise comparisons for more detailed information.

Kind regards,

Jeff