Bat Foraging Calls
M. Roberts
2025-09-26
For my thesis project, I set out acoustics across Fort Campbell Military Base and collected bat calls over 8 days. My field technicians and I completed habitat assessments for each acoustic placed and recorded characteristics such as vegetation cover, mean canopy cover, basal area, snags, water presence, etc. For this mixed-models example, I am going to examine the three vegetation covers I recorded, hardwood, softwood, and mixed, at 19 acoustic locations and see if that fixed effect has an influence on foraging calls. A poisson distribution was used throughout my example due to calls being counts throughout the night.
Let’s read in our data set and look at vegetation cover as our fixed effect while looking at vegetation cover and site ID as our random effects.
ad <- read.csv("acoustics.data.csv")
ad1 <- glmer(num_calls ~ veg_cover + (veg_cover | site_id), family=poisson(), data=ad)
summary(ad1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: num_calls ~ veg_cover + (veg_cover | site_id)
## Data: ad
##
## AIC BIC logLik -2*log(L) df.resid
## 5382.7 5407.1 -2682.4 5364.7 102
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -18.0988 -1.4681 -0.3253 1.3374 31.5268
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## site_id (Intercept) 2.716 1.6480
## veg_covermixed 0.292 0.5404 -0.74
## veg_coversoftwood 2.716 1.6480 -1.00 0.74
## Number of obs: 111, groups: site_id, 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.1856 0.4461 4.899 9.62e-07 ***
## veg_covermixed 1.1892 0.7991 1.488 0.136694
## veg_coversoftwood 1.7339 0.4489 3.862 0.000112 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vg_cvrm
## veg_covrmxd -0.558
## vg_cvrsftwd -0.994 0.555
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')flexplot::visualize(ad1, plot="model", sample=19) +
ggtitle("Random intercepts-random slopes") +
xlab("Vegetation Cover") +
ylab("Number of Calls")
This summary concluded that the hardwood and softwood vegetation cover has significance in foraging calls (p < 0.05) . It’s worth taking a look at the same fixed effect, vegetation cover, but only using site ID as a random effect.
ad2 <- glmer(num_calls ~ veg_cover + (1|site_id), family=poisson(), data=ad)
summary(ad2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: num_calls ~ veg_cover + (1 | site_id)
## Data: ad
##
## AIC BIC logLik -2*log(L) df.resid
## 5379.9 5390.7 -2686.0 5371.9 107
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -18.0982 -1.4685 -0.3279 1.3304 31.5290
##
## Random effects:
## Groups Name Variance Std.Dev.
## site_id (Intercept) 2.365 1.538
## Number of obs: 111, groups: site_id, 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.1912 0.4169 5.256 1.47e-07 ***
## veg_covermixed 1.1613 0.8832 1.315 0.189
## veg_coversoftwood 1.7270 1.5939 1.084 0.279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vg_cvrm
## veg_covrmxd -0.471
## vg_cvrsftwd -0.262 0.123flexplot::visualize(ad2, plot="model", sample=19) +
ggtitle("Random intercepts-same slopes") +
xlab("Vegetation Cover") +
ylab("Number of Calls")
This summary concluded that the hardwood vegetation cover is the only vegetation cover that significantly effects foraging calls (p < 0.05). Now that we have compared both random slope and random intercepts with random slopes and same intercepts, lets complete a comparative analysis.
results <- parameters::compare_parameters(ad1, ad2, select = "{estimate}<br>({se})|{p}")
print_html(results)| Parameter |
ad1
|
ad2
|
||
|---|---|---|---|---|
| Log-Mean(SE) | p | Log-Mean(SE) | p | |
| (Intercept) | 2.19 (0.45) |
<0.001 | 2.19 (0.42) |
<0.001 |
| veg cover (mixed) | 1.19 (0.80) |
0.137 | 1.16 (0.88) |
0.189 |
| veg cover (softwood) | 1.73 (0.45) |
<0.001 | 1.73 (1.59) |
0.279 |
| Observations | 111 | 111 | ||
arm::display(ad1)## glmer(formula = num_calls ~ veg_cover + (veg_cover | site_id),
## data = ad, family = poisson())
## coef.est coef.se
## (Intercept) 2.19 0.45
## veg_covermixed 1.19 0.80
## veg_coversoftwood 1.73 0.45
##
## Error terms:
## Groups Name Std.Dev. Corr
## site_id (Intercept) 1.65
## veg_covermixed 0.54 -0.74
## veg_coversoftwood 1.65 -1.00 0.74
## Residual 1.00
## ---
## number of obs: 111, groups: site_id, 19
## AIC = 5382.7, DIC = 4131.4
## deviance = 4748.1arm::display(ad2)## glmer(formula = num_calls ~ veg_cover + (1 | site_id), data = ad,
## family = poisson())
## coef.est coef.se
## (Intercept) 2.19 0.42
## veg_covermixed 1.16 0.88
## veg_coversoftwood 1.73 1.59
##
## Error terms:
## Groups Name Std.Dev.
## site_id (Intercept) 1.54
## Residual 1.00
## ---
## number of obs: 111, groups: site_id, 19
## AIC = 5379.9, DIC = 4124.1
## deviance = 4748.0performance::compare_performance(ad1, ad2)## # Comparison of Model Performance Indices
##
## Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | RMSE
## ----------------------------------------------------------------------------
## ad1 | glmerMod | 5382.7 (0.195) | 5384.5 (0.107) | 5407.1 (<.001) | 104.160
## ad2 | glmerMod | 5379.9 (0.805) | 5380.3 (0.893) | 5390.7 (>.999) | 104.160
##
## Name | Sigma | Score_log | Score_spherical | R2 (cond.) | R2 (marg.) | ICC
## ----------------------------------------------------------------------------
## ad1 | 1 | -23.684 | 0.060 | | |
## ad2 | 1 | -23.684 | 0.060 | 0.993 | 0.138 | 0.991From this comparison model, my second model (ad2) is better to use. As of now, my sample size is small, so this conclusion might change, but for now the only random effect that is worth noting would be site ID. Ad2 has a lower AIC than ad1 indicating that ad2 is a better fit model.