Unrealistically low density estimates

[horchard]

Hi there,

I previously fitted a model with homogeneous density to select my detection parameters, the top performing was the null model which produced realistic density estimates for my study species (0.00044/ha, 44/1000km2). I am now trying to fit a density model, testing a few different spatial covariates (rasters at 1000m resolution). I scaled and centered them all when adding them to my mask with the following script:

conifercover7 = rast("path to file/conifercover7_1000m.tif")
conifercover7 = scale(conifercover7, center = T, scale = T)
conifercover7_poly = terra::as.polygons(conifercover7, trunc = F, dissolve = T, values = T)
conifercover7_poly = st_as_sf(conifercover7_poly)
mask = addCovariates(mask, conifercover7_poly)

I started by running a separate model for each covariate, and most of them produced reasonable outputs. However, my conifer forest cover covariate, which is the top fitting, is producing unrealistically low density estimates. When plotting the density surface, it predicts density to be zero almost everywhere aside from two small regions in the study area where conifer cover has the highest values. After scaling, values ranged from -1 to 3 (average of 0.4) and it appears that density is 0 aside from where conifer cover = 3. The raw raster data seems to have a good variation of values between 0 and 1. I didn't have this issue with any of the other covariates, is there a different way that the raster should be scaled?

Here is the model script and summary:

model.3 <- secr.fit(capthist='caps', model = list(D~conifercover7, sigma~1, g0~1), details = list(maxdistance = 50000), buffer = 30000, trace = TRUE, method = 'BFGS', detectfn = 'HEX', CL = FALSE, verify = FALSE, mask = mask, ncores = 4)

> summary(model.3)
$versiontime
[1] "5.1.0, run 15:12:01 17 Apr 2025, elapsed 8653.23 s"

$traps
Detector Number Spacing UsagePct
count 2379 2999.188 300

$capthist
Occasions Detections Animals Detectors Moves Animals2
1 280 245 2379 35 31

$mask
Cells Spacing Area
45524 1000 4552400

$modeldetails
CL fixed distribution hcov
FALSE none poisson

$AICtable
model detectfn npar logLik AIC AICc
D~conifercover7 lambda0~1 sigma~1 hazard exponential 4 -1121.04 2250.079 2250.246

$coef
beta SE.beta lcl ucl
D -12.201445 0.50195906 -13.185267 -11.217623
D.conifercover7 3.001695 0.21379646 2.582662 3.420729
lambda0 -5.800216 0.21901588 -6.229479 -5.370953
sigma 8.858707 0.07829967 8.705242 9.012171

$predicted
link estimate SE.estimate lcl ucl
D log 5.023192e-06 2.688925e-06 1.878068e-06 1.343532e-05
lambda0 log 3.026901e-03 6.709693e-04 1.970478e-03 4.649699e-03
sigma log 7.035379e+03 5.517133e+02 6.034464e+03 8.202312e+03

Any advice is greatly appreciated!

Thank you,

Hailey

[murray.efford]

It's very hard to diagnose this issue, if it is one, without seeing the detail. Maybe send me the data offline (capthist object and mask bundled as an RData file with save(..., file =). That's a lot of detectors!

[murray.efford]

Comments after checking the data that Hailey sent offline -

I think this is doing what it oughta. The covariate conifercover7 takes the integer values -1:3, and very few habitat cells are in class 3. The covariate is highly predictive of density (possibly no animals centred in classes -1,0,1 that are most common). So the loglinear model is struggling to contrast class 2 and class <2. It does this quite well, but as a result class 3 is modelled as extreme density.

I think you are thrown by the predicted density from the raw call to predict() that (from memory) uses the first level of conifercover7 (0). Better to specify conifercover7 in newdata:

predict(model.3, newdata=data.frame(conifercover7=-1:3))

or look at the whole predicted surface
predsurf <- predictDsurface(model.3)
mean(covariates(predsurf)$D.0) * 1e5
[1] 35.89872 mean density per 1000 km^2 across mask2 - but is that a meaningful number?

This looks reasonable to me:
plot(predsurf, breaks = c(seq(0,300,20), 5000), scale = 1e5)

By the way - I sped up fitting with a coarser mask,
mask2 <- make.mask(traps, buffer = 30000, type='trapbuffer', spacing = 2000)
mask2 <- addCovariates(mask2, mask)
and you can drop the argument 'buffer' as it is overridden by 'mask'.

Murray

[horchard]

Thanks so much for taking a look at this! It makes much more sense now.
I had been incorrectly interpreting the summary() density estimate as average per ha across the mask area, rather than per ha when covariate = 0.
35.89/1000km^2 is a very reasonable number. I'm relieved to hear this all looks as it should! Again Murray this has been incredibly helpful and I really appreciate your time!

Hailey