Predictions and variance estimation in SDM with random effect

[Phil Bouchet]

Dear list,

In a study of coastal dolphins, I have built a density surface model which includes: (1) a single smooth of depth, (2) a smooth-factor interaction between year (factor) and geographic coordinates (bivariate smooth of x and y), and (3) a random effect term for the transect ID. 

Following the guidance in Pedersen et al. (2019) (https://peerj.com/articles/6876/), my model is specified as:

abundance.est ~ year_f + s(x, y, m = 2) + s(x, y, by = year_f, 

    bs = c("tp", "tp"), m = 1) + s(depth) + s(transect_f, bs = "re")

(Note that the data support a hazard-rate detection function model with group size as a covariate, hence why the response is the estimated abundance per segment, rather than the counts).

The rationale for including a random effect term in this instance is that (1) dolphin abundance is expected to be more similar within a transect than between transects, and (2) a subset of transects were sampled multiple times over the course of the study period. Including the random effect term seems to have a positive effect on model performance, and e.g., propels the %deviance explained from ~5% without it to > 20% with it.

Given the above model structure, what is the correct way of obtaining predictions and estimates of variance from the model (conditional on the random effect)?

As explained for example here (https://stackoverflow.com/questions/61595497/how-to-predict-gam-model-with-random-effect-in-r), I have so far used the "exclude" argument in function "predict" to generate predictions of abundance, and the results look sensible.

I cannot, however, seem to find a similar option in "dsm.var.gam" – which means that repeat runs of the function with a different value for the transect ID in "new.data" give slightly different results. What am I missing?

Any advice much appreciated,

Best wishes,

Phil

[David Lawrence Miller]

Hi Phil, hi listfolks,

I've not implemented exclude as part of dsm.var.gam yet, so that's on
me, sorry.

In the meantime, since you're only using dsm.var.gam (and just adding
the squared CVs of the detection function and GAM) things are a bit
easier. You should be able to use whichever method you'd use for a gam
on the dsm object (e.g., predict.gam(..., exclude="s(transect_f)")) to
get the GAM bit of the prediction or with se.fit=TRUE the uncertainty.
Using the $fit and $se.fit elements you can calculate the CVs which you
can square and add to the average detection function p CV^2 from the
summary output of the detection function. That's then the total CV
assuming independence between the two model components.

It'll be worth trying that in the non-random-effects case first to check
the output from that procedure is the same as dsm.var.gam's.

Feel free to contact me off-list if you have problems and we can post a
full solution to the list if necessary.

cheers,
--dave

> --
> You received this message because you are subscribed to the Google
> Groups "distance-sampling" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to distance-sampl...@googlegroups.com
> <mailto:distance-sampl...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/distance-sampling/a8959a4e-f6d8-439d-8c37-a49c74f1b1a5n%40googlegroups.com
> <https://groups.google.com/d/msgid/distance-sampling/a8959a4e-f6d8-439d-8c37-a49c74f1b1a5n%40googlegroups.com?utm_medium=email&utm_source=footer>.