Modelling with group size and estimating abundance with group detections

[IJCorreia]

Good afternoon,

First of all, thank you for your work on this package and for this forum, both are very useful! I have been exploring them and, along with the online documentation provided, I have been able to dive into inlabru!

I would like to ask for your help in the analysis I am trying to do with some data using inlabru. The context is the following. The data were collected during an aerial survey in three coastal regions in Greenland meant for distance sampling analysis. The GPS locations of each detection were also recorded. Also, some spatial covariates are provided in a raster object

We thought about fitting a point process model that accounts for both detectability and location of objects of interest following your example available here.

However, since each detection corresponds to one group of objects of interest (instead of a single object of interest), I also have the group size variable and that must be accounted for too. I saw that this question was already raised here and thus, if I understood it correctly, group size needs to be accounted in a separate process and then a joint model must be estimated in a similar way as shown here.

After considering all this, I tried to fit a joint model for the objects’ location (and detectability) and group size. Below is an example of the model components, formula and likelihoods for what I call “jointfit2” which is the full model, i.e., the model with all covariates provided:

cmp.def2 <- ~ Intercept(1) + size.Intercept(1) + sigma_theta(1, prec.linear = 1) +
          mySPDE(main = geometry, model = matern) + size.SPDE(main = geometry, model = matern) +
          regiRW(main = f.regi(data.trunc), model = "rw2") +
          elev(main = f.elev(data.trunc), model = "linear") + slop(main = f.slop(data.trunc), model = "linear") +
          aspe(main = f.aspe(data.trunc), model = "linear") + roug(main = f.roug(data.trunc), model = "linear") +
          flow(main = f.flow(data.trunc), model = "linear") + dist(main = f.dist(data.trunc), model = "linear")

fml.defcoord2 <- geometry + distance ~ Intercept + log(hn(distance, sigma.bru(sigma_theta))) +
              mySPDE + regiRW + elev + slop + aspe + roug + flow + dist

#just for group size
fml.defsize <- size ~ size.Intercept + size.SPDE

#for geometry + distance`
lik.coord2 <- inlabru::like(family = "cp", formula = fml.defcoord2,
                            data = data.trunc, samplers = tran.inlabru,
                            domain = list(geometry = meshUTM,
                                          distance = fmesher::fm_mesh_1d(seq(0, max(bins), length.out = 30))))

#for group size
lik.size <- inlabru::like(family = "poisson", formula = fml.defsize, 
                          data = data.trunc)

#full
jointfit2 <- inlabru::bru(cmp.def2, lik.coord2, lik.size,
                          options = list(control.inla = list(int.strategy = "eb"),
                                         bru_max_iter = 1,
                                         control.predictor = list(compute = TRUE)))

The same was done for other models (both model components and formulas were modified accordingly):

• jointfit0 - a null model (no covariates),
• jointfit1 - similar to the full model (jointfit2) but without the region covariate
  Below is a print of the summary of the jointfit2 object (full model):

When calculating and plotting the model predictions for the whole study region, some issues arise:

1. when the model includes the region covariate, we get “bands” where estimates are really low (print below);
2. when calculating the abundance estimate, it “explodes” (print and code below)

#prediction points
pxl.size <- fmesher::fm_pixels(mesh = meshUTM, dims = c(150, 150), mask = sf::st_as_sfc(border.utm))
predpts.size <- fmesher::fm_int(domain = meshUTM, samplers = border.utm) 

#prediction plots
Lambda.pred0 <- predict(object = jointfit0, newdata = pxl.size,
                        ~ exp(mySPDE + Intercept) * exp(size.Intercept + size.SPDE))
Lambda.pred1 <- predict(object = jointfit1, newdata = pxl.size,
                        ~ exp(mySPDE + Intercept + regiRW) * exp(size.Intercept + size.SPDE))
Lambda.pred2 <- predict(object = jointfit2, newdata = pxl.size,
                        ~ exp(mySPDE + Intercept + regiRW + elev + slop + aspe + roug + flow + dist) * exp(size.Intercept + size.SPDE))

#estimating abundance
Lambda.pred0sum <- predict(object = jointfit0, newdata = predpts.size,
                           ~ sum(weight * exp(mySPDE + Intercept) * exp(size.Intercept + size.SPDE)), n.samples = 1000)
Lambda.pred1sum <- predict(object = jointfit1, newdata = predpts.size,
                           ~ sum(weight * exp(mySPDE + Intercept + regiRW) * exp(size.Intercept + size.SPDE)), n.samples = 1000)
Lambda.pred2sum <- predict(object = jointfit2, newdata = predpts.size,
                           ~ sum(weight * exp(mySPDE + Intercept + regiRW + elev + slop + aspe + roug + flow + dist) * exp(size.Intercept + size.SPDE)), n.samples = 1000)

Can you please help me understand why this happens when I add the region covariate to the model?
I also tried modifying the model components and formula by adding the variable "size" containing the group size in my data.trunc object but it did not work... Or maybe I am not understanding how to do it...

If there is any more information needed, let me know and I will provide it asap.

Thank you in advance for your attention!
Iúri J. F. Correia

[finnlindgren]

Sorry, I missed this due to the holidays!

I think the problem is that you're hardcoding the input data set in the component definitions, so that the internally constructed data required for the "cp" likelihood isn't used for covariate evaluations. Change to this instead:

cmp.def2 <- ~ Intercept(1) + size.Intercept(1) + sigma_theta(1, prec.linear = 1) +
          mySPDE(main = geometry, model = matern) + size.SPDE(main = geometry, model = matern) +
          regiRW(main = f.regi(.data.), model = "rw2") +
          elev(main = f.elev(.data.), model = "linear") + slop(main = f.slop(.data.), model = "linear") +
          aspe(main = f.aspe(.data.), model = "linear") + roug(main = f.roug(.data.), model = "linear") +
          flow(main = f.flow(.data.), model = "linear") + dist(main = f.dist(.data.), model = "linear")```

Depending on how your covariates are stored, you likely don't need the f.name() functions at all; e.g. if they are SpatRaster objects you can just give the whole object as component input; inlabru will then call eval_spatial() that has a method for SpatRaster objects.

[IJCorreia]

I understand, it's ok. Thank you for the suggestion but unfortunately it does not work... I tried setting ".data." but the model estimates keep skyrocketing (update: I had to change the region model to "factor_full")... This is applicable to both models 0 and 2. Model 1 crashes R and session aborts... Also, for model 2, having the region covariate, those "bands" do not appear anymore but rather random pixels get very high/low values:

[finnlindgren]

If you use "factor.full" you shouldn't also have an intercept in the model, as that gives an unidentifiable model.

But the main problem seems to be that the covariate values aren't extracted correctly. This is likely due to some bug in your f.covar() functions, or that your data storage is inconsistent/incomplete (e.g. incorrectly specified crs information). Try using the output from integrationpoints<-fm_int(mesh, mypolygon) as the second argument of values = eval_spatial(mycovariateobject, integrationpoints) and plot with ggplot()+geom_sf(data = integrationpoints, aes(color=values)) and see if that looks correct. Also, values_f=f.mycovariatobject(integrationpoints) should match values if your f.*() function is correct.

[IJCorreia]

You are definitely right! Reverted to "factor_contrast". About the second part, I tested your suggestions, indeed it looks correct and values and values_f match, any other ideas? I also noted that the last plot I sent was not correct, below is what I obtained

[finnlindgren]

Do you have consistent crs information on the object's? I may need the covariate object and mesh to debug it.

[IJCorreia]

Yes, indeed, the CRS seem to be consistent (print below). I'll send you an .Rdata file via e-mail, there you can find the mesh (meshUTM), covariates (dem.new.inlabru) and polygons (border.utm)