Packages used in present analysis
Bat ecological niche & Leptospira. Bat individual parameters: species, sex, age, stage of annual cycle (Models IS)
Bat ecological niche & Leptospira: Ecology of all sampled bats (Models ES)
Bat ecological niche & Leptospira: ecological parameters only for adult bats (Models EA)

Packages used in present analysis

library(lme4)

## Loading required package: Matrix

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(broom.mixed)
library(ggplot2)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.2     ✔ tidyr     1.3.2
## ✔ readr     2.2.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(svglite)
library(gt)

Bat ecological niche & Leptospira. Bat individual parameters: species, sex, age, stage of annual cycle (Models IS)

lepto_individual <- read.csv("D:/downloads/lepto_stat/leptospira_test.csv")
lepto_individual$Bat.species   <- factor(lepto_individual$Bat.species)
lepto_individual$Bat.sex  <- factor(lepto_individual$Bat.sex)
lepto_individual$Bat.age  <- factor(lepto_individual$Bat.age)
lepto_individual$Physiological.condition  <- factor(lepto_individual$Physiological.condition)
lepto_individual$Male.colonies <- factor(lepto_individual$Male.colonies)
lepto_individual$Annual.cycle <- factor(lepto_individual$Annual.cycle)
lepto_individual$Physiographic.region   <- factor(lepto_individual$Physiographic.region)
lepto_individual$Leptospira.PCR <- factor(lepto_individual$Leptospira.PCR, levels = c("0", "1"))
lepto_individual$Knowledge <- factor(lepto_individual$Knowledge)


lepto_individual_adults <- lepto_individual %>%
  filter(Bat.age == "Adultus")

lepto_myo <- read.csv("D:/downloads/lepto_stat/myo_test.csv")
lepto_myo$Bat.species   <- factor(lepto_myo$Bat.species)
lepto_myo$Bat.age  <- factor(lepto_myo$Bat.age)
lepto_myo$Bat.sex  <- factor(lepto_myo$Bat.sex)
lepto_myo$Leptospira.PCR <- factor(lepto_myo$Leptospira.PCR, levels = c("0", "1"))
lepto_myo$Annual.cycle <- factor(lepto_myo$Annual.cycle)


lepto_sub <- read.csv("D:/downloads/lepto_stat/leptospira_test.csv")
lepto_sub$Bat.species   <- factor(lepto_sub$Bat.species)
lepto_sub$Bat.sex  <- factor(lepto_sub$Bat.sex)
lepto_sub$Bat.age  <- factor(lepto_sub$Bat.age)
lepto_sub$Physiological.condition  <- factor(lepto_sub$Physiological.condition)
lepto_sub$Migration  <- factor(lepto_sub$Migration)
lepto_sub$Foraging.strategy  <- factor(lepto_sub$Foraging.strategy)
lepto_sub$Summer.roost.size  <- factor(lepto_sub$Summer.roost.size)
lepto_sub$Male.colonies <- factor(lepto_sub$Male.colonies)
lepto_sub$Interspecies.colonies  <- factor(lepto_sub$Interspecies.colonies)
lepto_sub$Roost.type <- factor(lepto_sub$Roost.type)
lepto_sub$Annual.cycle <- factor(lepto_sub$Annual.cycle)
lepto_sub$Physiographic.region   <- factor(lepto_sub$Physiographic.region)
lepto_sub$Leptospira.PCR <- factor(lepto_sub$Leptospira.PCR, levels = c("0", "1"))
lepto_sub$Knowledge <- factor(lepto_sub$Knowledge)

lepto_sub <- read.csv("D:/downloads/lepto_stat/leptospira_test.csv")
lepto_sub_ad <- lepto_sub %>%
  filter(Bat.age == "Adultus")

lepto_sub_ad$Bat.species   <- factor(lepto_sub_ad$Bat.species)
lepto_sub_ad$Bat.sex  <- factor(lepto_sub_ad$Bat.sex)
lepto_sub_ad$Bat.age  <- factor(lepto_sub_ad$Bat.age)
lepto_sub_ad$Physiological.condition  <- factor(lepto_sub_ad$Physiological.condition)
lepto_sub_ad$Migration  <- factor(lepto_sub_ad$Migration)
lepto_sub_ad$Foraging.strategy  <- factor(lepto_sub_ad$Foraging.strategy)
lepto_sub_ad$Summer.roost.size  <- factor(lepto_sub_ad$Summer.roost.size)
lepto_sub_ad$Male.colonies <- factor(lepto_sub_ad$Male.colonies)
lepto_sub_ad$Interspecies.colonies  <- factor(lepto_sub_ad$Interspecies.colonies)
lepto_sub_ad$Roost.type <- factor(lepto_sub_ad$Roost.type)
lepto_sub_ad$Annual.cycle <- factor(lepto_sub_ad$Annual.cycle)
lepto_sub_ad$Physiographic.region   <- factor(lepto_sub_ad$Physiographic.region)
lepto_sub_ad$Leptospira.PCR <- factor(lepto_sub_ad$Leptospira.PCR, levels = c("0", "1"))
lepto_sub_ad$Knowledge <- factor(lepto_sub_ad$Knowledge)

Model IS-1. Species (only well-sampled species)

GLMM with a binomial error distribution and logit link, including bat species as a fixed effect and two random intercepts — physiographic region (9 levels) and bat age class (2 levels). The model was based on 2,152 individual bats, including only species where more than 30 individuals were collected. Age influence is known for rodents and bats of other regions. Place of collection may predispose both the species composition and the outdoor conditions for Leptospira storage in the environment.

lepto_filtered <- subset(lepto_individual, Knowledge %in% c("Sufficiently", "Acceptable"))
lepto_filtered <- droplevels(lepto_filtered)

lepto_filtered$Bat.species  <- relevel(lepto_filtered$Bat.species, ref = "Plecotus auritus")

model_species_lepto_filtered  <- glmer(Leptospira.PCR ~ Bat.species + (1|Bat.age) + (1|Physiographic.region),
                                             data = lepto_filtered,
                                             family = binomial(link= "logit")
)
summary(model_species_lepto_filtered)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Leptospira.PCR ~ Bat.species + (1 | Bat.age) + (1 | Physiographic.region)
##    Data: lepto_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2055.4    2174.6   -1006.7    2013.4      2131 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2803 -0.5700 -0.3288 -0.0003 10.6993 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  Physiographic.region (Intercept) 0.1140   0.3376  
##  Bat.age              (Intercept) 0.4171   0.6458  
## Number of obs: 2152, groups:  Physiographic.region, 9; Bat.age, 2
## 
## Fixed effects:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -3.937356   0.767968  -5.127 2.94e-07 ***
## Bat.speciesBarbastella barbastellus -12.831570 550.304714  -0.023 0.981397    
## Bat.speciesEptesicus nilssonii        0.877417   0.687436   1.276 0.201828    
## Bat.speciesMurina hilgendorfi         1.711270   0.750438   2.280 0.022586 *  
## Bat.speciesMyotis blythii             3.415038   0.640589   5.331 9.76e-08 ***
## Bat.speciesMyotis bombinus            1.593563   0.839687   1.898 0.057722 .  
## Bat.speciesMyotis brandtii            0.215317   0.724850   0.297 0.766428    
## Bat.speciesMyotis dasycneme           2.659109   0.609188   4.365 1.27e-05 ***
## Bat.speciesMyotis daubentonii         2.870497   0.602465   4.765 1.89e-06 ***
## Bat.speciesMyotis mystacinus          0.003247   0.934103   0.003 0.997226    
## Bat.speciesMyotis nattereri          -0.005841   0.835686  -0.007 0.994423    
## Bat.speciesMyotis petax               3.270650   0.645400   5.068 4.03e-07 ***
## Bat.speciesMyotis sibiricus           1.170307   0.876296   1.336 0.181707    
## Bat.speciesNyctalus leisleri          2.985154   0.658695   4.532 5.85e-06 ***
## Bat.speciesNyctalus noctula           2.723549   0.638786   4.264 2.01e-05 ***
## Bat.speciesPipistrellus nathusii      2.108929   0.623357   3.383 0.000717 ***
## Bat.speciesPipistrellus pygmaeus      2.370300   0.679124   3.490 0.000483 ***
## Bat.speciesPlecotus ognevi            1.626224   0.697398   2.332 0.019709 *  
## Bat.speciesVespertilio murinus        1.376207   0.655831   2.098 0.035868 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 19 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

species_pvals <- broom.mixed::tidy(model_species_lepto_filtered, effects = "fixed") %>%
  filter(term != "(Intercept)") %>%
  transmute(
    species = str_remove(term, "^Bat.species"),
    estimate,
    std.error,
    statistic,
    p.value
  ) %>%
  arrange(p.value)

species_pvals_adj <- species_pvals %>%
  mutate(p_adj_BH = p.adjust(p.value, method = "BH")) %>%
  arrange(p_adj_BH)

species_pvals_adj

## # A tibble: 18 × 6
##    species                   estimate std.error statistic      p.value  p_adj_BH
##    <chr>                        <dbl>     <dbl>     <dbl>        <dbl>     <dbl>
##  1 Myotis blythii             3.42        0.641   5.33    0.0000000976   1.76e-6
##  2 Myotis petax               3.27        0.645   5.07    0.000000403    3.63e-6
##  3 Myotis daubentonii         2.87        0.602   4.76    0.00000189     1.14e-5
##  4 Nyctalus leisleri          2.99        0.659   4.53    0.00000585     2.63e-5
##  5 Myotis dasycneme           2.66        0.609   4.37    0.0000127      4.58e-5
##  6 Nyctalus noctula           2.72        0.639   4.26    0.0000201      6.03e-5
##  7 Pipistrellus pygmaeus      2.37        0.679   3.49    0.000483       1.24e-3
##  8 Pipistrellus nathusii      2.11        0.623   3.38    0.000717       1.61e-3
##  9 Plecotus ognevi            1.63        0.697   2.33    0.0197         3.94e-2
## 10 Murina hilgendorfi         1.71        0.750   2.28    0.0226         4.07e-2
## 11 Vespertilio murinus        1.38        0.656   2.10    0.0359         5.87e-2
## 12 Myotis bombinus            1.59        0.840   1.90    0.0577         8.66e-2
## 13 Myotis sibiricus           1.17        0.876   1.34    0.182          2.52e-1
## 14 Eptesicus nilssonii        0.877       0.687   1.28    0.202          2.59e-1
## 15 Myotis brandtii            0.215       0.725   0.297   0.766          9.20e-1
## 16 Barbastella barbastellus -12.8       550.     -0.0233  0.981          9.97e-1
## 17 Myotis nattereri          -0.00584     0.836  -0.00699 0.994          9.97e-1
## 18 Myotis mystacinus          0.00325     0.934   0.00348 0.997          9.97e-1

write.csv(species_pvals, "species_pvalues.csv", row.names = FALSE)

Model IS-1-1. Species (all collected bats)

Actually the same as previous model, however based on 2,455 individual bats, adding every collected species, disregarding was its sample size 1, 3, or 300. Bat species as a fixed effect and two random intercepts — physiographic region (9 levels) and bat age class (2 levels).

lepto_individual$Bat.species  <- relevel(lepto_individual$Bat.species, ref = "Plecotus auritus")

model_species  <- glmer(Leptospira.PCR ~ Bat.species + (1|Bat.age) + (1|Physiographic.region),
                                             data = lepto_individual,
                                             family = binomial(link= "logit")
)

## Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
## failure to converge in 10000 evaluations

## Warning in optwrap(optimizer, devfun, start, lower = lower, upper = upper, :
## convergence code 4 from Nelder_Mead: failure to converge in 10000 evaluations

summary(model_species)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Leptospira.PCR ~ Bat.species + (1 | Bat.age) + (1 | Physiographic.region)
##    Data: lepto_individual
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2188.4    2417.9   -1054.2    2108.4      2254 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2425 -0.5280 -0.3275 -0.0001 10.6689 
## 
## Random effects:
##  Groups               Name        Variance Std.Dev.
##  Physiographic.region (Intercept) 0.1045   0.3233  
##  Bat.age              (Intercept) 0.4067   0.6377  
## Number of obs: 2294, groups:  Physiographic.region, 9; Bat.age, 2
## 
## Fixed effects:
##                                         Estimate Std. Error z value Pr(>|z|)
## (Intercept)                           -3.922e+00  7.634e-01  -5.138 2.78e-07
## Bat.speciesBarbastella barbastellus   -1.574e+01  2.348e+03  -0.007 0.994652
## Bat.speciesBarbastella leucomelas     -1.241e+01  2.804e+03  -0.004 0.996467
## Bat.speciesEptesicus nilssonii         8.754e-01  6.873e-01   1.274 0.202775
## Bat.speciesEptesicus serotinus         1.511e+00  1.236e+00   1.222 0.221583
## Bat.speciesHypsugo savii              -1.445e+01  5.501e+03  -0.003 0.997904
## Bat.speciesMiniopterus fuliginosus    -1.337e+01  3.672e+03  -0.004 0.997094
## Bat.speciesMiniopterus schreibersi     4.759e+00  8.231e-01   5.782 7.37e-09
## Bat.speciesMurina hilgendorfi          1.705e+00  7.494e-01   2.275 0.022916
## Bat.speciesMyotis blythii              3.401e+00  6.397e-01   5.316 1.06e-07
## Bat.speciesMyotis bombinus             1.590e+00  8.377e-01   1.898 0.057756
## Bat.speciesMyotis brandtii             2.140e-01  7.247e-01   0.295 0.767789
## Bat.speciesMyotis dasycneme            2.653e+00  6.090e-01   4.357 1.32e-05
## Bat.speciesMyotis daubentonii          2.866e+00  6.023e-01   4.758 1.96e-06
## Bat.speciesMyotis davidii             -1.435e+01  3.191e+03  -0.004 0.996413
## Bat.speciesMyotis emarginatus         -1.389e+01  5.868e+03  -0.002 0.998111
## Bat.speciesMyotis ikonnikovi           4.526e+00  1.261e+00   3.590 0.000331
## Bat.speciesMyotis longicaudatus        2.314e+00  1.325e+00   1.746 0.080726
## Bat.speciesMyotis macrodactylus        4.559e+00  8.851e-01   5.151 2.59e-07
## Bat.speciesMyotis mystacinus           7.970e-04  9.341e-01   0.001 0.999319
## Bat.speciesMyotis nattereri           -7.186e-03  8.356e-01  -0.009 0.993138
## Bat.speciesMyotis petax                3.265e+00  6.441e-01   5.070 3.99e-07
## Bat.speciesMyotis sibiricus            1.163e+00  8.756e-01   1.329 0.183975
## Bat.speciesMyotis tschuliensis        -1.362e+01  1.208e+03  -0.011 0.991003
## Bat.speciesNyctalus lasiopterus       -1.305e+01  5.143e+03  -0.003 0.997975
## Bat.speciesNyctalus leisleri           2.974e+00  6.579e-01   4.520 6.18e-06
## Bat.speciesNyctalus noctula            2.716e+00  6.384e-01   4.254 2.10e-05
## Bat.speciesPipistrellus kuhlii         1.505e+00  1.220e+00   1.233 0.217438
## Bat.speciesPipistrellus nathusii       2.100e+00  6.232e-01   3.370 0.000750
## Bat.speciesPipistrellus pipistrellus  -1.502e+01  4.906e+03  -0.003 0.997557
## Bat.speciesPipistrellus pygmaeus       2.361e+00  6.787e-01   3.479 0.000504
## Bat.speciesPlecotus ognevi             1.619e+00  6.965e-01   2.324 0.020111
## Bat.speciesRhinolophus aff bocharicus -1.390e+01  5.582e+03  -0.002 0.998013
## Bat.speciesRhinolophus euryale         4.865e+00  1.296e+00   3.752 0.000175
## Bat.speciesRhinolophus ferrumequinum   1.318e+00  9.882e-01   1.334 0.182159
## Bat.speciesRhinolophus hipposideros    2.358e+00  9.306e-01   2.534 0.011275
## Bat.speciesVespertilio murinus         1.371e+00  6.556e-01   2.092 0.036473
## Bat.speciesVespertilio sinensis       -1.307e+01  3.160e+03  -0.004 0.996699
##                                          
## (Intercept)                           ***
## Bat.speciesBarbastella barbastellus      
## Bat.speciesBarbastella leucomelas        
## Bat.speciesEptesicus nilssonii           
## Bat.speciesEptesicus serotinus           
## Bat.speciesHypsugo savii                 
## Bat.speciesMiniopterus fuliginosus       
## Bat.speciesMiniopterus schreibersi    ***
## Bat.speciesMurina hilgendorfi         *  
## Bat.speciesMyotis blythii             ***
## Bat.speciesMyotis bombinus            .  
## Bat.speciesMyotis brandtii               
## Bat.speciesMyotis dasycneme           ***
## Bat.speciesMyotis daubentonii         ***
## Bat.speciesMyotis davidii                
## Bat.speciesMyotis emarginatus            
## Bat.speciesMyotis ikonnikovi          ***
## Bat.speciesMyotis longicaudatus       .  
## Bat.speciesMyotis macrodactylus       ***
## Bat.speciesMyotis mystacinus             
## Bat.speciesMyotis nattereri              
## Bat.speciesMyotis petax               ***
## Bat.speciesMyotis sibiricus              
## Bat.speciesMyotis tschuliensis           
## Bat.speciesNyctalus lasiopterus          
## Bat.speciesNyctalus leisleri          ***
## Bat.speciesNyctalus noctula           ***
## Bat.speciesPipistrellus kuhlii           
## Bat.speciesPipistrellus nathusii      ***
## Bat.speciesPipistrellus pipistrellus     
## Bat.speciesPipistrellus pygmaeus      ***
## Bat.speciesPlecotus ognevi            *  
## Bat.speciesRhinolophus aff bocharicus    
## Bat.speciesRhinolophus euryale        ***
## Bat.speciesRhinolophus ferrumequinum     
## Bat.speciesRhinolophus hipposideros   *  
## Bat.speciesVespertilio murinus        *  
## Bat.speciesVespertilio sinensis          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 38 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (Nelder_Mead) convergence code: 4 (failure to converge in 10000 evaluations)
## failure to converge in 10000 evaluations

Model IS-2. Age (sad/ad). Species adjusted

To test whether Leptospira PCR positivity differed between age classes, we fitted a GLMM with a binomial error distribution and logit link, with bat age (subadults vs. adult) as a fixed effect and bat species included as a random intercept. The model was based on 2,294 individual bats across 37 species.

lepto_individual$Bat.age  <- relevel(lepto_individual$Bat.age, ref = "Subadultus")


model_individual_age  <- glmer(Leptospira.PCR ~ Bat.age + (1|Bat.species),
                                             data = lepto_individual,
                                             family = binomial(link= "logit")
)
summary(model_individual_age)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.age + (1 | Bat.species)
##    Data: lepto_individual
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2221.5    2238.7   -1107.8    2215.5      2291 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6192 -0.5017 -0.3355 -0.1791  7.9423 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.845    1.358   
## Number of obs: 2294, groups:  Bat.species, 38
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.7838     0.3020  -9.217  < 2e-16 ***
## Bat.ageAdultus   1.2559     0.1541   8.150 3.63e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Bat.agAdlts -0.460

Model IS-2f. Age (sad/ad) only for females of well sampled species. Species adjusted

Specifying the same age model in one sex (females) and species collected up to 30 individuals. The model was based on 1,188 individual bats across 19 species.

lepto_females <- subset(lepto_filtered, Bat.sex %in% ("Female"))
lepto_females <- droplevels(lepto_females)

lepto_females$Bat.age  <- relevel(lepto_females$Bat.age, ref = "Subadultus")

model_females_age  <- glmer(Leptospira.PCR ~ Bat.age + (1|Bat.species),
                                             data = lepto_females,
                                             family = binomial(link= "logit")
)
summary(model_females_age)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.age + (1 | Bat.species)
##    Data: lepto_females
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1214.2    1229.4    -604.1    1208.2      1185 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1154 -0.5156 -0.3217  0.9857  6.6173 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.198    1.095   
## Number of obs: 1188, groups:  Bat.species, 19
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -3.2913     0.3608  -9.121   <2e-16 ***
## Bat.ageAdultus   2.0035     0.2431   8.240   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Bat.agAdlts -0.632

Model IS-2m. Age (sad/ad) only for males of well-sampled species. Species adjusted

Specifying the same age model in one sex (males) and species collected up to 30 individuals. The model was based on 955 individual bats across 19 species.

lepto_males <- subset(lepto_filtered, Bat.sex %in% ("Male"))
lepto_males <- droplevels(lepto_males)

lepto_males$Bat.age  <- relevel(lepto_males$Bat.age, ref = "Subadultus")


model_males_age  <- glmer(Leptospira.PCR ~ Bat.age + (1|Bat.species),
                                             data = lepto_males,
                                             family = binomial(link= "logit")
)
summary(model_males_age)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.age + (1 | Bat.species)
##    Data: lepto_males
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     842.8     857.3    -418.4     836.8       952 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.9908 -0.6037 -0.3034 -0.1685  6.2539 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.315    1.147   
## Number of obs: 955, groups:  Bat.species, 19
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.1704     0.3548  -6.116 9.57e-10 ***
## Bat.ageAdultus   0.1047     0.2358   0.444    0.657    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Bat.agAdlts -0.494

Model IS-2myo. Age for three Myotis species

To examine whether Leptospira PCR positivity varied across age categories within three Myotis species, where the age could be estimated more precisely (Myotis daubentonii, Myotis dasycneme, and Myotis petax), we fitted a GLMM with a binomial error distribution and logit link, with bat age (four levels: Subadultus [reference], Adultus of 1st year, Adultus of 2nd year, Adultus older than 2 years) as a fixed effect, and bat species and bat sex included as random intercepts. The model was based on 309 individuals.

lepto_myo$Bat.age  <- relevel(lepto_myo$Bat.age, ref = "Subadultus")

model_myo_age  <- glmer(Leptospira.PCR ~ Bat.age + (1|Bat.species) + (1|Bat.sex),
                                             data = lepto_myo,
                                             family = binomial(link= "logit")
)
summary(model_myo_age)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.age + (1 | Bat.species) + (1 | Bat.sex)
##    Data: lepto_myo
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     390.0     412.4    -189.0     378.0       303 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0039 -0.7149 -0.5097  1.0030  2.0757 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  Bat.species (Intercept) 4.939e-08 0.0002222
##  Bat.sex     (Intercept) 1.161e-02 0.1077382
## Number of obs: 309, groups:  Bat.species, 3; Bat.sex, 2
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -1.4044     0.2616  -5.369  7.9e-08 ***
## Bat.ageAdultus, 1st year   0.6766     0.3489   1.940  0.05243 .  
## Bat.ageAdultus, 2nd year   1.3556     0.4293   3.158  0.00159 ** 
## Bat.ageAdultus, older 2    1.3418     0.3226   4.159  3.2e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) B.A,1y B.A,2y
## Bt.gAdlt,1y -0.684              
## Bt.gAdlt,2y -0.552  0.415       
## Bt.gAdlt,o2 -0.742  0.554  0.447

Model IS-3-1. Sex (all collected bats). Species adjusted

To test whether Leptospira PCR positivity differed by sex, we fitted GLMM (binomial family, logit link) including bat sex as a fixed effect, and bat species as a random intercept. The model was based on 2,294 individuals.

lepto_individual$Bat.sex  <- relevel(lepto_individual$Bat.sex, ref = "Male")

model_individual_AnC  <- glmer(Leptospira.PCR ~ Annual.cycle + (1|Bat.species) + (1|Bat.age),
                                             data = lepto_individual,
                                             family = binomial(link= "logit")
)
summary(model_individual_AnC)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Annual.cycle + (1 | Bat.species) + (1 | Bat.age)
##    Data: lepto_individual
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2208.9    2237.6   -1099.4    2198.9      2289 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6081 -0.5600 -0.3044 -0.1764  6.1678 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.6502   1.2846  
##  Bat.age     (Intercept) 0.5549   0.7449  
## Number of obs: 2294, groups:  Bat.species, 38; Bat.age, 2
## 
## Fixed effects:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -2.1693     0.6002  -3.614 0.000301 ***
## Annual.cycleHibernation      -0.3655     0.1755  -2.083 0.037282 *  
## Annual.cycleSummer activity   0.3921     0.1718   2.282 0.022508 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Annl.H
## Annl.cyclHb -0.169       
## Annl.cyclSa -0.178  0.600

Model IS-3-2. Sex (all collected bats). Species, presence of male colonies, age adjusted

To test whether Leptospira PCR positivity differed by sex, we fitted GLMM (binomial family, logit link) including bat sex as a fixed effect, and bat species, presence of male colonies, and bat as random intercepts. The model was based on 2,285 individuals.

lepto_individual$Bat.sex  <- relevel(lepto_individual$Bat.sex, ref = "Male")

model_individual_sex  <- glmer(Leptospira.PCR ~ Bat.sex + (1|Bat.species) + (1|Male.colonies) + (1|Bat.age),
                                             data = lepto_individual,
                                             family = binomial(link= "logit")
)
summary(model_individual_sex)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.sex + (1 | Bat.species) + (1 | Male.colonies) +  
##     (1 | Bat.age)
##    Data: lepto_individual
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2210.5    2239.2   -1100.3    2200.5      2280 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8439 -0.5125 -0.3421 -0.1662  8.2563 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Bat.species   (Intercept) 1.5848   1.2589  
##  Male.colonies (Intercept) 0.4164   0.6453  
##  Bat.age       (Intercept) 0.4720   0.6870  
## Number of obs: 2285, groups:  Bat.species, 38; Male.colonies, 2; Bat.age, 2
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    -2.0511     0.7328  -2.799  0.00512 **
## Bat.sexFemale   0.3466     0.1134   3.056  0.00224 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Bat.sexFeml -0.075

Model IA-3. Sex (adult bats). Species, presence of male colonies adjusted

Almost the same model as previous one, but restricted the sample to adults only.It included the same fixed effect of sex, with bat species and presence of male colonies as random intercepts (age was necessarily dropped, since the subset is age-restricted by definition). The model was based on 1,696 individuals.

lepto_individual_adults$Bat.sex  <- relevel(lepto_individual_adults$Bat.sex, ref = "Male")

model_individualAD_sex  <- glmer(Leptospira.PCR ~ Bat.sex + (1|Bat.species) + (1|Male.colonies),
                                             data = lepto_individual_adults,
                                             family = binomial(link= "logit")
)
summary(model_individualAD_sex)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.sex + (1 | Bat.species) + (1 | Male.colonies)
##    Data: lepto_individual_adults
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1751.6    1773.3    -871.8    1743.6      1692 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0038 -0.6812 -0.3028  0.8906  5.7270 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Bat.species   (Intercept) 1.6957   1.3022  
##  Male.colonies (Intercept) 0.3221   0.5676  
## Number of obs: 1696, groups:  Bat.species, 37; Male.colonies, 2
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.6427     0.5149  -3.190  0.00142 ** 
## Bat.sexFemale   0.5961     0.1296   4.598 4.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Bat.sexFeml -0.152

Model IS-3-4. Sex (species with male colonies). Species adjusted

Most of the Palaearctic bat species form nursery colonies, containing females and juveniles. Males usually roost solitary. However, for some species — this was group with “male colony: present” status — both sexes are known to roost in colonies during warmer period. Gregariousness in winter was not taken into account. To test whether the sex effect on PCR positivity held specifically among bats sampled where male colonies were present, we fitted a GLMM (binomial family, logit link) with bat sex as a fixed effect and bat species as a random intercept, restricted to this subset. The model was based on 576 individuals of 8 species (Vespertilio murinus, Nyctalus noctula, Myotis dasycneme, M. petax, M. macrodactylus, Plecotus ognevi, Miniopterus fuliginosus, and M. schreibersi).

lepto_individual_male_present <- lepto_individual_adults %>%
  filter(Male.colonies == "Present")

lepto_individual_male_present$Bat.sex <- relevel(lepto_individual_male_present$Bat.sex, ref = "Male")

model_individual_sex_male_col_p  <- glmer(Leptospira.PCR ~ Bat.sex + (1|Bat.species),
                                             data = lepto_individual_male_present,
                                             family = binomial(link= "logit")
)
summary(model_individual_sex_male_col_p)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.sex + (1 | Bat.species)
##    Data: lepto_individual_male_present
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     719.7     732.8    -356.9     713.7       573 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7518 -0.8637 -0.4500  1.1263  2.5730 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 0.9843   0.9921  
## Number of obs: 576, groups:  Bat.species, 8
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -0.6175     0.4049  -1.525    0.127
## Bat.sexFemale   0.2932     0.1967   1.491    0.136
## 
## Correlation of Fixed Effects:
##             (Intr)
## Bat.sexFeml -0.282

Model IS-3-5. Sex (species without male colonies). Species adjusted

Most of the Palaearctic bat species form nursery colonies, containing females and juveniles. Males usually roost solitary.

lepto_individual_male_Unknown <- lepto_individual_adults %>%
  filter(Male.colonies == "Unknown")

lepto_individual_male_Unknown$Bat.sex <- relevel(lepto_individual_male_Unknown$Bat.sex, ref = "Male")

model_individual_sex_nomale_col_p  <- glmer(Leptospira.PCR ~ Bat.sex + (1|Bat.species),
                                             data = lepto_individual_male_Unknown,
                                             family = binomial(link= "logit")
)
summary(model_individual_sex_nomale_col_p)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Bat.sex + (1 | Bat.species)
##    Data: lepto_individual_male_Unknown
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1026.8    1041.9    -510.4    1020.8      1117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0753 -0.5042 -0.2806 -0.1553  6.3341 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.786    1.336   
## Number of obs: 1120, groups:  Bat.species, 29
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -2.4002     0.3370  -7.123 1.05e-12 ***
## Bat.sexFemale   0.8186     0.1740   4.705 2.54e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Bat.sexFeml -0.338

Model IS-4. Annual cycle (all collected bats) bats. Species adjusted

To test whether PCR positivity varied across stages of the annual activity cycle, we fitted a GLMM with annual cycle stage as a fixed effect (three levels, with the baseline/reference stage “Flux” when reporting). Bat annual life cycle defined as: hibernation, migratory flux, and summer activity. For each individual, the life stage was determined not only by the date of capture, but also by the circumstances of trapping such as collecting in a cave or on a flyway, physiological state. Bat species and bat age were included as random intercepts. The model was based on 2,294 individuals of 37 species.

lepto_individual$Annual.cycle  <- relevel(lepto_individual$Annual.cycle, ref = "Hibernation")

model_individual_AnC  <- glmer(Leptospira.PCR ~ Annual.cycle + (1|Bat.species) + (1|Bat.age),
                                             data = lepto_individual,
                                             family = binomial(link= "logit")
)
summary(model_individual_AnC)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Annual.cycle + (1 | Bat.species) + (1 | Bat.age)
##    Data: lepto_individual
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2208.9    2237.6   -1099.4    2198.9      2289 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6081 -0.5600 -0.3044 -0.1764  6.1678 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.6502   1.2846  
##  Bat.age     (Intercept) 0.5549   0.7449  
## Number of obs: 2294, groups:  Bat.species, 38; Bat.age, 2
## 
## Fixed effects:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -2.5349     0.5960  -4.253 2.11e-05 ***
## Annual.cycleFlux              0.3655     0.1755   2.083   0.0373 *  
## Annual.cycleSummer activity   0.7576     0.1555   4.874 1.10e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Annl.F
## Annl.cyclFl -0.124       
## Annl.cyclSa -0.143  0.466

Model IS-4myo. Annual cycle (three Myotis species). Species adjusted

To test whether the seasonal pattern in PCR positivity observed in the full dataset held within a more homogeneous subset of three closely ecological related well-sampled Myotis species, we fitted a GLMM (binomial family, logit link) with annual cycle stage as a fixed effect and bat species and bat age included as random intercepts. The model was based on 698 individuals, 3 species — Myotis daubentonii, Myotis dasycneme, and Myotis petax — and 2 age groups.

lepto_three_Myotis <- lepto_individual %>%
filter(Bat.species %in% c("Myotis daubentonii", "Myotis dasycneme", "Myotis petax"))

lepto_three_Myotis$Annual.cycle  <- relevel(lepto_three_Myotis$Annual.cycle, ref = "Hibernation")

model_individual_AnC  <- glmer(Leptospira.PCR ~ Annual.cycle + (1|Bat.species) + (1|Bat.age),
                                             data = lepto_three_Myotis,
                                             family = binomial(link= "logit")
)
summary(model_individual_AnC)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Annual.cycle + (1 | Bat.species) + (1 | Bat.age)
##    Data: lepto_three_Myotis
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     929.5     952.2    -459.7     919.5       693 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0405 -0.8241 -0.5372  1.1551  2.2243 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 0.02089  0.1445  
##  Bat.age     (Intercept) 0.41099  0.6411  
## Number of obs: 698, groups:  Bat.species, 3; Bat.age, 2
## 
## Fixed effects:
##                             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                  -1.0161     0.4975  -2.042   0.0411 *
## Annual.cycleFlux              0.1180     0.2207   0.535   0.5929  
## Annual.cycleSummer activity   0.3757     0.2047   1.835   0.0665 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Annl.F
## Annl.cyclFl -0.146       
## Annl.cyclSa -0.237  0.280

Extracting odds ratio and confidence intervals

models_list_individual <- list(
  "Age" = model_individual_age,
  "Sex" = model_individual_sex,
  "Sex adults" = model_individualAD_sex,
  "Sex male colonies present" = model_individual_sex_male_col_p,
  "Anual cycle" = model_individual_AnC
)

extract_or_ci_ind <- function(model, model_name) {
  tidy_res <- broom.mixed::tidy(model, exponentiate = TRUE, conf.int = TRUE)
  tidy_res %>%
    filter(effect == "fixed", term != "(Intercept)") %>%   # drops BOTH intercept types
    mutate(model = model_name) %>%
    rename(OR = estimate, OR_low = conf.low, OR_high = conf.high)
}

or_df_ind <- bind_rows(lapply(names(models_list_individual), function(nm) {
  extract_or_ci_ind(models_list_individual[[nm]], nm)
}))

write.csv(or_df_ind, "odds_ratios_indmodels.csv", row.names = FALSE)

# ---- Strip the variable-name prefix from each term, longest-name first ----
predictor_prefixes_models_list_individual <- c("Age", "Sex",
                         "Sex male colonies present", "Sex & age in male colonies",
                         "Sex adults", "Anual cycle")
predictor_models_list_individual <- predictor_prefixes_models_list_individual[order(-nchar(models_list_individual))]
strip_pattern_models_list_individual <- paste0("^(", paste(predictor_prefixes_models_list_individual, collapse = "|"), ")")

or_df_ind <- or_df_ind %>%
  mutate(
    clean_term = str_remove(term, strip_pattern_models_list_individual),
    clean_term = sub("^(.)", "\\L\\1", clean_term, perl = TRUE)
  )


or_df_ind$p.adj <- p.adjust(or_df_ind$p.value, method = "BH")

Bat ecological niche & Leptospira: Ecology of all sampled bats (Models ES)

lepto_sub <- read.csv("D:/downloads/lepto_stat/leptospira_test.csv")
lepto_sub$Bat.species   <- factor(lepto_sub$Bat.species)
lepto_sub$Bat.sex  <- factor(lepto_sub$Bat.sex)
lepto_sub$Bat.age  <- factor(lepto_sub$Bat.age)
lepto_sub$Physiological.condition  <- factor(lepto_sub$Physiological.condition)
lepto_sub$Migration  <- factor(lepto_sub$Migration)
lepto_sub$Foraging.strategy  <- factor(lepto_sub$Foraging.strategy)
lepto_sub$Summer.roost.size  <- factor(lepto_sub$Summer.roost.size)
lepto_sub$Male.colonies <- factor(lepto_sub$Male.colonies)
lepto_sub$Interspecies.colonies  <- factor(lepto_sub$Interspecies.colonies)
lepto_sub$Roost.type <- factor(lepto_sub$Roost.type)
lepto_sub$Annual.cycle <- factor(lepto_sub$Annual.cycle)
lepto_sub$Physiographic.region   <- factor(lepto_sub$Physiographic.region)
lepto_sub$Leptospira.PCR <- factor(lepto_sub$Leptospira.PCR, levels = c("0", "1"))
lepto_sub$Knowledge <- factor(lepto_sub$Knowledge)

Model ES-1. Summer roost size

Summer roost number was estimated as small for bats forming nursery colonies for up to 10-30 individuals, moderate for up to 100, and large for 100 up to thousands of individuals). Ranging these categories we reckon average colony size, not maximum registered (for example, Myotis dasycneme and Vespertilio murinus were considered species with moderate maternity colony size, notwithstanding colonies up to 500 individuals were occasionally recorded). Gregariousness in winter was not taken into account.

lepto_sub$Summer.roost.size  <- relevel(lepto_sub$Summer.roost.size, ref = "Small")

model_Summer.roost.size_sp  <- glmer(Leptospira.PCR ~ Summer.roost.size + (1|Bat.species),
                                             data = lepto_sub,
                                             family = binomial(link= "logit")
)
summary(model_Summer.roost.size_sp)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Summer.roost.size + (1 | Bat.species)
##    Data: lepto_sub
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2444.9    2468.2   -1218.5    2436.9      2451 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6380 -0.6013 -0.3742 -0.1886  4.4424 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.047    1.023   
## Number of obs: 2455, groups:  Bat.species, 38
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -2.1846     0.2910  -7.507 6.04e-14 ***
## Summer.roost.sizeLarge      1.8678     0.6212   3.007  0.00264 ** 
## Summer.roost.sizeModerate   0.8471     0.4554   1.860  0.06284 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Smm..L
## Smmr.rst.sL -0.434       
## Smmr.rst.sM -0.631  0.283

Model ES-2. Male colonies

According to the likelihood of forming male colonies bats fall into two categories: such colonies they have been registered in the literature or by authors’ unpublished observations, or unknown for the present moment.

lepto_sub$Male.colonies  <- relevel(lepto_sub$Male.colonies, ref = "Unknown")

model_Male.colonies_sp  <- glmer(Leptospira.PCR ~ Male.colonies + (1|Bat.species),
                                     data = lepto_sub,
                                     family = binomial(link= "logit")
)
summary(model_Male.colonies_sp)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Male.colonies + (1 | Bat.species)
##    Data: lepto_sub
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2443.8    2461.3   -1218.9    2437.8      2452 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6414 -0.6062 -0.3792 -0.1816  4.4744 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.335    1.156   
## Number of obs: 2455, groups:  Bat.species, 38
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -2.0478     0.2726  -7.512 5.81e-14 ***
## Male.coloniesPresent   1.4019     0.5167   2.713  0.00666 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Ml.clnsPrsn -0.520

Model ES-3. Interspecies colonies

According to the likelihood of forming interspecies colonies colonies bats fall into two categories: such colonies they have been registered in the literature or by authors’ unpublished observations, or unknown for the present moment.

lepto_sub$Interspecies.colonies  <- relevel(lepto_sub$Interspecies.colonies, ref = "Unknown")

model_Interspecies.colonies_sp  <- glmer(Leptospira.PCR ~ Interspecies.colonies + (1|Bat.species),
                                 data = lepto_sub,
                                 family = binomial(link= "logit")
)
summary(model_Interspecies.colonies_sp)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Interspecies.colonies + (1 | Bat.species)
##    Data: lepto_sub
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2446.8    2464.3   -1220.4    2440.8      2452 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5888 -0.6136 -0.3712 -0.1580  4.6930 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.527    1.236   
## Number of obs: 2455, groups:  Bat.species, 38
## 
## Fixed effects:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -2.5002     0.4754  -5.259 1.45e-07 ***
## Interspecies.coloniesPresent   1.0960     0.5485   1.998   0.0457 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Intrspcs.cP -0.855

Model ES-4. Roost type

Preferred roosting types were classified according to Bat Roosts in Trees: A Guide to Identification and Assessment for Tree-Care and Ecology Professionals by Rotherham (2019). This classification comprises void-roosting and crevice-roosting types, as well as a comprehensive type for species that follow both roosting strategies. For seven bat species their preferred roost type was unknown (Supplementary Table S1), so these bats (34 individuals) were excluded from roosting model sample.

lepto_sub$Roost.type  <- relevel(lepto_sub$Roost.type, ref = "Crevice-roosting")

model_Roost.type_sp  <- glmer(Leptospira.PCR ~ Roost.type + (1|Bat.species),
                                         data = lepto_sub,
                                         family = binomial(link= "logit")
)
summary(model_Roost.type_sp)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Roost.type + (1 | Bat.species)
##    Data: lepto_sub
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2429.1    2452.3   -1210.6    2421.1      2417 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6272 -0.6085 -0.3643 -0.1877  4.2307 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.077    1.038   
## Number of obs: 2421, groups:  Bat.species, 30
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -2.1738     0.3428  -6.342 2.27e-10 ***
## Roost.typeComprehensive   0.5031     0.5267   0.955 0.339457    
## Roost.typeVoid-roosting   1.7417     0.5245   3.321 0.000897 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rst.tC
## Rst.typCmpr -0.644       
## Rst.typVd-r -0.652  0.420

Model ES-5. Foraging strategy

Foraging strategies were classified as synthesis of findings made by Kruskop (1998) and Fenton & Bogdanowicz (2002). Kruskop on the basis of wing morphology highlighted three foraging types — gleaning, cluttered-area foraging, and open-space foraging — ranging each category by the remoteness from substrate (soil or water). Fenton & Bogdanowicz (2002) based on literature data and authors’ bat measuring suggest dividing bats into four groups: gleaning, trawling, foraging over water, and aerial feeding. Taking about Leptospira, which can maintain both in water and damp soil (Bierque et al., 2020), we assume both the a type of substrate and a degree of attachment to it may be worth noting, and divide our sample as: aerial feeders with cluttered-area foraging (Barbastella, Eptesicus, Hypsugo, most Myotis, Pipistrellus, and Vespertilio), foraging over water (Myotis dasycneme, M. daubentonii, M. macrodactylus, M. petax), aerial feeders with open-space foraging (Miniopterus and Nyctalus), and gleaning (Murina, Plecotus, Rhinolophus, Myotis blythii, M. bombinus, M. emarginatus, M. nattereri, and M. tschuliensis). Vespertilio murinus according to Kruskop (1998) was considered a cluttered-area forager, however, as mentioned in his work and was observed during authors research, it can also be treated as an open-space foraging species.

lepto_sub$Foraging.strategy  <- relevel(lepto_sub$Foraging.strategy, ref = "Aerial feeding, cluttered area foraging")

model_Foraging.strategy_sp  <- glmer(Leptospira.PCR ~ Foraging.strategy  + (1|Bat.species),
                             data = lepto_sub,
                             family = binomial(link= "logit")
)
summary(model_Foraging.strategy_sp)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Foraging.strategy + (1 | Bat.species)
##    Data: lepto_sub
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2442.1    2471.2   -1216.1    2432.1      2450 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5462 -0.6210 -0.3625 -0.1924  4.2533 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 0.9103   0.9541  
## Number of obs: 2455, groups:  Bat.species, 38
## 
## Fixed effects:
##                                                      Estimate Std. Error
## (Intercept)                                           -2.2556     0.3174
## Foraging.strategyAerial feeding, open-space foraging   1.5041     0.6189
## Foraging.strategyForaging over water                   2.1279     0.5933
## Foraging.strategyGleaner                               0.3487     0.4683
##                                                      z value Pr(>|z|)    
## (Intercept)                                           -7.107 1.18e-12 ***
## Foraging.strategyAerial feeding, open-space foraging   2.430 0.015094 *  
## Foraging.strategyForaging over water                   3.586 0.000336 ***
## Foraging.strategyGleaner                               0.745 0.456524    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) F.Afof Fr.Fow
## Frgn.Af,o-f -0.496              
## Frgng.stFow -0.541  0.265       
## Frgng.strtG -0.667  0.337  0.358

Model ES-6. Migration

Migration activity types were divided into three categories according to average flight distances. Short-range migrants — bats that do not typically undertake long flights, wintering and breeding within 30–50 km, up to a maximum 150 — genera Plecotus, Eptesicus, Barbastella, Rhinolophus, Murina, Hypsugo, Miniopterus, and most Myotis; medium-range migrants — migration between summer and winter roosts within one to three regions, usually from 100 to 500 km — Myotis daubentonii, M. petax, M. dasycneme; and long-range migrations — seasonal flux may occur over considerable distances, up to 1,000–2,500 km or more — genera Vespertilio, Nyctalus, Pipistrellus. Within the migration types two subtypes were identified: partial, where part of the population migrates and the remainder stays within the territory of summer roosts; and differential, where only bats of a particular sex, usually females, undertake migration. Yet, this parameters were not included in the analysis and can be found only in Supplementary Table S1.

lepto_sub$Migration <- relevel(lepto_sub$Migration, ref = "Short-range")


model_Migration_sp  <- glmer(Leptospira.PCR ~ Migration + (1|Bat.species),
                                             data = lepto_sub,
                                             family = binomial(link= "logit")
)
summary(model_Migration_sp)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Migration + (1 | Bat.species)
##    Data: lepto_sub
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2448.6    2471.8   -1220.3    2440.6      2451 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5250 -0.6114 -0.3695 -0.1812  4.4654 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.453    1.205   
## Number of obs: 2455, groups:  Bat.species, 38
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            -1.9394     0.2897  -6.695 2.16e-11 ***
## MigrationLong-range     0.3259     0.5652   0.577   0.5642    
## MigrationMiddle-range   1.5723     0.7582   2.074   0.0381 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MgrtL-
## MgrtnLng-rn -0.465       
## MgrtnMddl-r -0.382  0.178

Extracting odds ratio and confidence intervals and Visualisation

models_list <- list(
  "migration" = model_Migration_sp,
  "foraging" = model_Foraging.strategy_sp,
  "summer roost size" = model_Summer.roost.size_sp,
  "male colonies" = model_Male.colonies_sp,
  "interspecies colonies" = model_Interspecies.colonies_sp,
  "roost type" = model_Roost.type_sp
)

extract_or_ci <- function(model, model_name) {
  tidy_res <- broom.mixed::tidy(model, exponentiate = TRUE, conf.int = TRUE)
  tidy_res %>%
    filter(effect == "fixed", term != "(Intercept)") %>%   # drops fixed AND sd__(Intercept)
    mutate(model = model_name) %>%
    rename(OR = estimate, OR_low = conf.low, OR_high = conf.high)
}

or_df <- bind_rows(lapply(names(models_list), function(nm) {
  extract_or_ci(models_list[[nm]], nm)
}))

write.csv(or_df, "odds_ratios_models.csv", row.names = FALSE)

# ---- Strip the variable-name prefix from each term ----
predictor_prefixes <- c("Interspecies.colonies", "Summer.roost.size",
                         "Foraging.strategy", "Male.colonies",
                         "Roost.type", "Migration")
predictor_prefixes <- predictor_prefixes[order(-nchar(predictor_prefixes))]
strip_pattern <- paste0("^(", paste(predictor_prefixes, collapse = "|"), ")")

or_df <- or_df %>%
  mutate(
    clean_term = str_remove(term, strip_pattern),
    clean_term = sub("^(.)", "\\L\\1", clean_term, perl = TRUE) 
  )

or_df$p.adj <- p.adjust(or_df$p.value, method = "BH")

or_plot <- or_df %>%
  mutate(y_label = paste0(model, ": ", clean_term)) %>%
  ggplot(aes(x = OR, y = fct_inorder(y_label))) +
  geom_point(color = "#325153", size = 2) +
  geom_errorbarh(aes(xmin = OR_low, xmax = OR_high), height = 0.15, color = "#325153") +
  geom_vline(xintercept = 1, linetype = "dashed", color = "#651e14", linewidth = 0.8) +
  scale_x_log10(breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10)) +
  labs(
    x = "Odds Ratio (95% CI)",
    y = "Model & Predictor",
    title = "Odds Ratios for ecologicy based models for all 2,455 bats (GLMM, binomial)"
  ) +
  theme_minimal(base_size = 11, base_family = "Helvetica Light") +
  theme(axis.text.y = element_text(size = 10))

## Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggsave(
  filename = "odds_ratios_all.svg",
  plot = or_plot,
  device = "svg",
  width = 8, height = 6, units = "in"
)

## `height` was translated to `width`.

The same as a table:

or_df %>%
  select(-group, -effect, -term) %>%
  mutate(
    p.value = ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value)),
    p.adj   = ifelse(p.adj   < 0.001, "<0.001", sprintf("%.3f", p.adj))
  ) %>%
  gt() %>%
  fmt_number(columns = c(OR, std.error, statistic, OR_low, OR_high), decimals = 2) %>%
  cols_label(
    OR = "Odds Ratio",
    std.error = "SE",
    statistic = "z",
    OR_low = "95% CI (lower)",
    OR_high = "95% CI (upper)",
    p.value = "p",
    p.adj = "p (BH-adjusted)",
    model = "Model",
    clean_term = "Predictor"
  ) %>%
  tab_header(title = "Odds ratios for ecological predictors of PCR positivity")

Odds Ratio	SE	z	p	95% CI (lower)	95% CI (upper)	Model	Predictor	p (BH-adjusted)
Odds ratios for ecological predictors of PCR positivity
1.39	0.78	0.58	0.564	0.46	4.19	migration	long-range	0.564
4.82	3.65	2.07	0.038	1.09	21.29	migration	middle-range	0.070
4.50	2.79	2.43	0.015	1.34	15.14	foraging	aerial feeding, open-space foraging	0.033
8.40	4.98	3.59	<0.001	2.62	26.86	foraging	foraging over water	0.004
1.42	0.66	0.74	0.457	0.57	3.55	foraging	gleaner	0.502
6.47	4.02	3.01	0.003	1.92	21.87	summer roost size	large	0.010
2.33	1.06	1.86	0.063	0.96	5.70	summer roost size	moderate	0.086
4.06	2.10	2.71	0.007	1.48	11.19	male colonies	present	0.018
2.99	1.64	2.00	0.046	1.02	8.77	interspecies colonies	present	0.072
1.65	0.87	0.96	0.339	0.59	4.64	roost type	comprehensive	0.415
5.71	2.99	3.32	<0.001	2.04	15.95	roost type	void-roosting	0.005

Bat ecological niche & Leptospira: ecological parameters only for adult bats (Models EA)

Because infection prevalence and immune responses may differ between age groups (Aguillon et al., 2025), ecological analyses were performed both in the full dataset (n = 2,455) and after exclusion of subadult individuals (n = 1,726).

Model EA-1. Summer roost size

lepto_sub_ad$Summer.roost.size  <- relevel(lepto_sub_ad$Summer.roost.size, ref = "Small")

model_Summer.roost.size_spAd  <- glmer(Leptospira.PCR ~ Summer.roost.size + (1|Bat.species),
                                             data = lepto_sub_ad,
                                             family = binomial(link= "logit")
)
summary(model_Summer.roost.size_spAd)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Summer.roost.size + (1 | Bat.species)
##    Data: lepto_sub_ad
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1774.2    1795.9    -883.1    1766.2      1700 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6657 -0.8130 -0.3348  0.9972  4.8588 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.258    1.122   
## Number of obs: 1704, groups:  Bat.species, 37
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -2.0988     0.3166  -6.628  3.4e-11 ***
## Summer.roost.sizeLarge      1.7433     0.6745   2.585  0.00975 ** 
## Summer.roost.sizeModerate   1.2048     0.5006   2.407  0.01609 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Smm..L
## Smmr.rst.sL -0.435       
## Smmr.rst.sM -0.619  0.284

Model EA-2. Male colonies

lepto_sub_ad$Male.colonies  <- relevel(lepto_sub_ad$Male.colonies, ref = "Unknown")

model_Male.colonies_spAd  <- glmer(Leptospira.PCR ~ Male.colonies + (1|Bat.species),
                                     data = lepto_sub_ad,
                                     family = binomial(link= "logit")
)
summary(model_Male.colonies_spAd)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Male.colonies + (1 | Bat.species)
##    Data: lepto_sub_ad
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1774.0    1790.4    -884.0    1768.0      1701 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6941 -0.8162 -0.3371  0.9553  4.8939 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.611    1.269   
## Number of obs: 1704, groups:  Bat.species, 37
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -1.8783     0.3008  -6.244 4.27e-10 ***
## Male.coloniesPresent   1.4133     0.5688   2.485    0.013 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Ml.clnsPrsn -0.522

Model EA-3. Interspecies colonies

According to the likelihood of forming interspecies colonies bats fall into two categories: such colonies they have been registered in the literature or by authors’ unpublished observations, or unknown for the present moment.

lepto_sub_ad$Interspecies.colonies  <- relevel(lepto_sub_ad$Interspecies.colonies, ref = "Unknown")

model_Interspecies.colonies_spAd  <- glmer(Leptospira.PCR ~ Interspecies.colonies + (1|Bat.species),
                                 data = lepto_sub_ad,
                                 family = binomial(link= "logit")
)
summary(model_Interspecies.colonies_spAd)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Interspecies.colonies + (1 | Bat.species)
##    Data: lepto_sub_ad
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1774.3    1790.6    -884.2    1768.3      1701 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6325 -0.8126 -0.3457  0.9300  5.2104 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.648    1.284   
## Number of obs: 1704, groups:  Bat.species, 37
## 
## Fixed effects:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -2.4990     0.4935  -5.064  4.1e-07 ***
## Interspecies.coloniesPresent   1.3763     0.5716   2.408   0.0161 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Intrspcs.cP -0.850

Model EA-4. Roost type

lepto_sub_ad$Roost.type  <- relevel(lepto_sub_ad$Roost.type, ref = "Crevice-roosting")

model_Roost.type_spAd  <- glmer(Leptospira.PCR ~ Roost.type + (1|Bat.species),
                                         data = lepto_sub_ad,
                                         family = binomial(link= "logit")
)
summary(model_Roost.type_spAd)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Roost.type + (1 | Bat.species)
##    Data: lepto_sub_ad
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1758.9    1780.5    -875.4    1750.9      1668 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6875 -0.8090 -0.3359  1.0445  4.5474 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.273    1.128   
## Number of obs: 1672, groups:  Bat.species, 30
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -2.0208     0.3792  -5.330 9.84e-08 ***
## Roost.typeComprehensive   0.5558     0.5761   0.965  0.33471    
## Roost.typeVoid-roosting   1.8233     0.5711   3.193  0.00141 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rst.tC
## Rst.typCmpr -0.651       
## Rst.typVd-r -0.659  0.430

Model EA-5. Foraging strategy

lepto_sub_ad$Foraging.strategy  <- relevel(lepto_sub_ad$Foraging.strategy, ref = "Aerial feeding, cluttered area foraging")

model_Foraging.strategy_spAd  <- glmer(Leptospira.PCR ~ Foraging.strategy  + (1|Bat.species),
                             data = lepto_sub_ad,
                             family = binomial(link= "logit")
)
summary(model_Foraging.strategy_spAd)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Foraging.strategy + (1 | Bat.species)
##    Data: lepto_sub_ad
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1769.4    1796.6    -879.7    1759.4      1699 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6879 -0.8198 -0.3359  1.0368  4.5572 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 0.9612   0.9804  
## Number of obs: 1704, groups:  Bat.species, 37
## 
## Fixed effects:
##                                                      Estimate Std. Error
## (Intercept)                                           -2.0605     0.3349
## Foraging.strategyAerial feeding, open-space foraging   2.1698     0.6629
## Foraging.strategyForaging over water                   2.0133     0.6145
## Foraging.strategyGleaner                               0.1956     0.4864
##                                                      z value Pr(>|z|)    
## (Intercept)                                           -6.153 7.62e-10 ***
## Foraging.strategyAerial feeding, open-space foraging   3.273  0.00106 ** 
## Foraging.strategyForaging over water                   3.277  0.00105 ** 
## Foraging.strategyGleaner                               0.402  0.68752    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) F.Afof Fr.Fow
## Frgn.Af,o-f -0.490              
## Frgng.stFow -0.551  0.268       
## Frgng.strtG -0.675  0.335  0.370

Model EA-6. Migration

lepto_sub_ad$Migration  <- relevel(lepto_sub_ad$Migration , ref = "Short-range")


model_Migration_spAd  <- glmer(Leptospira.PCR ~ Migration + (1|Bat.species),
                                             data = lepto_sub_ad,
                                             family = binomial(link= "logit")
)
summary(model_Migration_spAd)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Leptospira.PCR ~ Migration + (1 | Bat.species)
##    Data: lepto_sub_ad
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1776.5    1798.3    -884.2    1768.5      1700 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5569 -0.8172 -0.3369  0.9224  4.9041 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Bat.species (Intercept) 1.619    1.272   
## Number of obs: 1704, groups:  Bat.species, 37
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            -1.8942     0.3099  -6.113 9.77e-10 ***
## MigrationLong-range     1.0374     0.6423   1.615   0.1063    
## MigrationMiddle-range   1.6232     0.8018   2.024   0.0429 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MgrtL-
## MgrtnLng-rn -0.471       
## MgrtnMddl-r -0.386  0.182

Extracting odds ratio and confidence intervals and Visualisation

models_listAd <- list(
  "Migration" = model_Migration_spAd,
  "Foraging" = model_Foraging.strategy_spAd,
  "Summer roost size" = model_Summer.roost.size_spAd,
  "Male colonies" = model_Male.colonies_spAd,
  "Interspecies colonies" = model_Interspecies.colonies_spAd,
  "Roost type" = model_Roost.type_spAd
)

extract_or_ciAd <- function(model, model_name) {
  tidy_res <- broom.mixed::tidy(model, exponentiate = TRUE, conf.int = TRUE)
  tidy_res %>%
    filter(effect == "fixed", term != "(Intercept)") %>%   # drops fixed AND sd__(Intercept)
    mutate(model = model_name) %>%
    rename(OR = estimate, OR_low = conf.low, OR_high = conf.high)
}

or_dfAd <- bind_rows(lapply(names(models_listAd), function(nm) {
  extract_or_ciAd(models_listAd[[nm]], nm)
}))

write.csv(or_dfAd, "odds_ratios_models.csv", row.names = FALSE)


or_dfAd$p.adj <- p.adjust(or_dfAd$p.value, method = "BH")

# ---- Strip the variable-name prefix from each term ----
predictor_prefixesAd <- c("Interspecies.colonies", "Summer.roost.size",
                         "Foraging.strategy", "Male.colonies",
                         "Roost.type", "Migration")
predictor_prefixesAd <- predictor_prefixesAd[order(-nchar(predictor_prefixesAd))]
strip_pattern <- paste0("^(", paste(predictor_prefixesAd, collapse = "|"), ")")

or_dfAd <- or_dfAd %>%
  mutate(
    clean_term = str_remove(term, strip_pattern),
    clean_term = sub("^(.)", "\\L\\1", clean_term, perl = TRUE)  # lowercase first letter
  )

or_Ad_plot <- or_dfAd %>%
  mutate(y_label = paste0(model, ": ", clean_term)) %>%
  ggplot(aes(x = OR, y = fct_inorder(y_label))) +
  geom_point(color = "#325153", size = 2) +
  geom_errorbarh(aes(xmin = OR_low, xmax = OR_high), height = 0.15, color = "#325153") +
  geom_vline(xintercept = 1, linetype = "dashed", color = "#651e14", linewidth = 0.8) +
  scale_x_log10(breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10)) +
  labs(
    x = "Odds Ratio (95% CI)",
    y = "Model & Predictor",
    title = "Odds Ratios for ecology based models for adult bats (GLMM, binomial)"
  ) +
  theme_minimal(base_size = 11, base_family = "Helvetica Light") +
  theme(axis.text.y = element_text(size = 10))

ggsave(
  filename = "odds_ratios_adult.svg",
  plot = or_Ad_plot,
  device = "svg",
  width = 8, height = 6, units = "in"
)

## `height` was translated to `width`.

The same as table

library(gt)

or_dfAd %>%
  select(-group, -effect, -term) %>%
  mutate(
    p.value = ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value)),
    p.adj   = ifelse(p.adj   < 0.001, "<0.001", sprintf("%.3f", p.adj))
  ) %>%
  gt() %>%
  fmt_number(columns = c(OR, std.error, statistic, OR_low, OR_high), decimals = 2) %>%
  cols_label(
    OR = "Odds Ratio",
    std.error = "SE",
    statistic = "z",
    OR_low = "95% CI (lower)",
    OR_high = "95% CI (upper)",
    p.value = "p",
    p.adj = "p (BH-adjusted)",
    model = "Model",
    clean_term = "Predictor"
  ) %>%
  tab_header(title = "Odds ratios for ecological predictors of PCR positivity for adult bats")

Odds Ratio	SE	z	p	95% CI (lower)	95% CI (upper)	Model	p (BH-adjusted)	Predictor
Odds ratios for ecological predictors of PCR positivity for adult bats
2.82	1.81	1.62	0.106	0.80	9.94	Migration	0.130	long-range
5.07	4.06	2.02	0.043	1.05	24.40	Migration	0.059	middle-range
8.76	5.80	3.27	0.001	2.39	32.10	Foraging	0.005	aerial feeding, open-space foraging
7.49	4.60	3.28	0.001	2.25	24.97	Foraging	0.005	foraging over water
1.22	0.59	0.40	0.688	0.47	3.15	Foraging	0.688	gleaner
5.72	3.86	2.58	0.010	1.52	21.44	Summer roost size	0.025	large
3.34	1.67	2.41	0.016	1.25	8.90	Summer roost size	0.025	moderate
4.11	2.34	2.48	0.013	1.35	12.53	Male colonies	0.025	present
3.96	2.26	2.41	0.016	1.29	12.14	Interspecies colonies	0.025	present
1.74	1.00	0.96	0.335	0.56	5.39	Roost type	0.368	comprehensive
6.19	3.54	3.19	0.001	2.02	18.97	Roost type	0.005	void-roosting

Session Info:

sessionInfo()

## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Moscow
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gt_1.3.0            svglite_2.2.2       lubridate_1.9.5    
##  [4] forcats_1.0.1       stringr_1.6.0       purrr_1.2.2        
##  [7] readr_2.2.0         tidyr_1.3.2         tibble_3.3.1       
## [10] tidyverse_2.0.0     ggplot2_4.0.3       broom.mixed_0.2.9.7
## [13] dplyr_1.2.1         lme4_2.0-6          Matrix_1.7-5       
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.60          bslib_0.11.0       lattice_0.22-9    
##  [5] tzdb_0.5.0         vctrs_0.7.3        tools_4.6.1        Rdpack_2.6.6      
##  [9] generics_0.1.4     parallel_4.6.1     pkgconfig_2.0.3    RColorBrewer_1.1-3
## [13] S7_0.2.2           lifecycle_1.0.5    compiler_4.6.1     farver_2.1.2      
## [17] textshaping_1.0.5  codetools_0.2-20   htmltools_0.5.9    sass_0.4.10       
## [21] yaml_2.3.12        pillar_1.11.1      furrr_0.4.0        nloptr_2.2.1      
## [25] jquerylib_0.1.4    MASS_7.3-66        cachem_1.1.0       reformulas_0.4.4  
## [29] boot_1.3-32        nlme_3.1-169       parallelly_1.48.0  tidyselect_1.2.1  
## [33] digest_0.6.39      stringi_1.8.7      future_1.75.0      listenv_1.0.0     
## [37] splines_4.6.1      fastmap_1.2.0      grid_4.6.1         cli_3.6.6         
## [41] magrittr_2.0.5     utf8_1.2.6         broom_1.0.13       withr_3.0.3       
## [45] scales_1.4.0       backports_1.5.1    timechange_0.4.0   rmarkdown_2.31    
## [49] globals_0.19.1     otel_0.2.0         ragg_1.5.2         hms_1.1.4         
## [53] evaluate_1.0.5     knitr_1.51         rbibutils_2.4.1    rlang_1.3.0       
## [57] Rcpp_1.1.2         glue_1.8.1         xml2_1.6.0         rstudioapi_0.19.0 
## [61] minqa_1.2.8        jsonlite_2.0.0     R6_2.6.1           fs_2.1.0          
## [65] systemfonts_1.3.2

Bat ecological traits behind Leptospira distribution and genetics in the north-eastern Palaearctic. Supplementary Material S2

Statistical analysis was conducted by Alex Lintu Viskontene & Maria S. Malyavina

2026

Packages used in present analysis

Bat ecological niche & Leptospira. Bat individual parameters: species, sex, age, stage of annual cycle (Models IS)

Model IS-1. Species (only well-sampled species)

Model IS-1-1. Species (all collected bats)

Model IS-2. Age (sad/ad). Species adjusted

Model IS-2f. Age (sad/ad) only for females of well sampled species. Species adjusted

Model IS-2m. Age (sad/ad) only for males of well-sampled species. Species adjusted

Model IS-2myo. Age for three Myotis species

Model IS-3-1. Sex (all collected bats). Species adjusted

Model IS-3-2. Sex (all collected bats). Species, presence of male colonies, age adjusted

Model IA-3. Sex (adult bats). Species, presence of male colonies adjusted

Model IS-3-4. Sex (species with male colonies). Species adjusted

Model IS-3-5. Sex (species without male colonies). Species adjusted

Model IS-4. Annual cycle (all collected bats) bats. Species adjusted

Model IS-4myo. Annual cycle (three Myotis species). Species adjusted

Extracting odds ratio and confidence intervals

Bat ecological niche & Leptospira: Ecology of all sampled bats (Models ES)

Model ES-1. Summer roost size

Model ES-2. Male colonies

Model ES-3. Interspecies colonies

Model ES-4. Roost type

Model ES-5. Foraging strategy

Model ES-6. Migration

Extracting odds ratio and confidence intervals and Visualisation

Bat ecological niche & Leptospira: ecological parameters only for adult bats (Models EA)

Model EA-1. Summer roost size

Model EA-2. Male colonies

Model EA-3. Interspecies colonies

Model EA-4. Roost type

Model EA-5. Foraging strategy

Model EA-6. Migration

Extracting odds ratio and confidence intervals and Visualisation

Session Info: