DOI: 10.5061/dryad.0p2ngf29s
title: Data and code from: A timid choice: Risk-taking behavior predicts individualized niche in a varying landscape of safety

description: Individual niche specialization predicts a match between an individual’s
 phenotype and environment. Yet, whether animals achieve this match through
 phenotypic change (niche conformance) or by selecting the environment
 (niche choice), remains unexplored. Individual variation in risk-taking
 behavior should contribute to the process of realizing individualized
 niches. Using wild populations of common voles (Microtus arvalis), we
 provide first evidence of how animals in the wild realize individualized
 niches and match their risk-taking behavior to microhabitats of varying
 safety. Under natural conditions, risk-averse individuals used safer
 microhabitats than risk-prone conspecifics. This correlation strengthened
 when we in situ experimentally made the environment riskier. A change in
 microhabitat use did not result in a change of risk-taking behavior. Our
 results demonstrate that animals choose environments matching their
 risk-taking phenotype and support the hypothesis that individual-level
 selection of environments of varying safety can be an extended phenotypic
 trait.
description: Study site (trapping & animals) The
 study was conducted in a meadow landscape with regular mowing practices
 between September and November 2019, and in October 2022 in the Havelland
 district, Brandenburg, Germany (52°44’00” N 12°12’44” E). The study
 site<em> </em>is on an island surrounded by the river Havel,
 and the smaller rivers Gülper Havel and Pirre. Each year we set up two
 grid systems (50 m x 50 m) for capturing voles with six by six Ugglan
 traps (Grahnab Sweden, special no. 2, with added shrew exit) evenly spaced
 out 10 m apart. One grid in 2022 had 4 x 6 traps due to space
 restrictions. We conducted capture-mark-recapture (CMR) for at least two
 weeks in each grid. Grids in 2019 were surrounded by an electric fence (90
 cm high; mesh size 17.5 cm wide and 9 to 15 cm high) to restrict
 disturbance of traps by raccoons. Traps were activated around sunset and
 checked after sunrise (maximum of 8 hours between checking active traps).
 In 2019 traps were additionally activated and checked during the day.
 Traps were baited with oat flakes and apples or cucumber. In colder
 temperatures (&lt;5°C) traps were provided with wood shavings as
 nesting material to prevent hypothermia. We weighed, measured head width,
 and determined sex and age class based on body mass (≥17 g = mature/adult,
 &lt;17 g = juvenile/immature) of all common voles. Individuals were
 fur-marked with a unique pattern for individual identification. After this
 handling procedure, animals were released next to the trap they came
 from. Based on head width (mm) and body mass (g), we
 calculated body condition using the scaled body mass index (Peig &
 Green 2009). This index was calculated as an individual's body mass
 (in g) standardized to the mean head width to retain a body size measure
 of all individuals in the population with a population-specific allometric
 relationship. Vegetation We surveyed
 the natural microhabitat in an area of one meter around each trap by
 measuring the height (average maximum height of five measuring points) and
 estimating cover at ten centimeters above ground by counting the number of
 grasses/branches touching a 1-meter ruler. To
 experimentally make the habitat less safe, we mowed and removed the grass
 10 m x 10 m around every other trap in the grid (except the first grid in
 2019 where half of the home range of VHF-tracked individuals was mowed).
 Thereby, we intended to affect the home ranges of all individuals to some
 degree. The size of the mowed patches was chosen to resemble home range
 sizes of male common voles estimated from VHF-tracking (mean: 202
 m<sup>2</sup> Jacob & Hempel 2003); 3,830 to 5,950
 m<sup>2</sup>: A. Schirmer, J. Eccard & M. Dammhahn,
 unpublished data). We decided against a more individual-tailored
 manipulation since this required work-intensive efforts of individual
 tracking and excluded smaller individuals. Based on a pilot VHF-telemetry
 tracking study in the study population, individuals' home ranges
 covered areas larger than 10 m x 10 m around each trap (F. Erixon, M.
 Gilmour, A. Schirmer, unpublished data). After mowing, the mowed patches
 had significantly lower vegetation height (&lt; 10 cm) and, thus,
 lower vegetation cover (Fig. SI 1). Behavioral
 assays Upon first capture, we assess among-individual
 differences in behaviors related to risk-taking of all adult common voles
 (n = 82) using a combined <em>emergence</em> and
 <em>open-field</em> test (Herde & Eccard 2013;
 Schirmer <em>et al.</em> 2019). Tests were conducted indoors
 and all handling and sample collection was done after testing. We
 reassessed recaptured animals after 5-17 days (median = 7).
 The setup consists of an opaque tube (11 cm diameter, 30 cm
 length) with a swing door connected to a round arena (120 cm diameter, 60
 cm height) divided into 16 equal-sized sections: eight at the edge and
 eight at the exposed central region (84 cm diameter). First, the animal is
 taken from the trap into a can that is placed in the tub and left to
 habituate for one minute. The experimenter then opens the swing door to
 start the <em>emergence</em> test. Once the vole exits the
 tube the swing door is closed and the open-field test starts. The
 experiments were scored in real time in 2019 and recorded on camera and
 scored from videos in 2022. Between tests, the setup was cleaned with 70%
 ethanol. In the emergence test, we measured the latency
 (seconds) to exit the tube with its head (including the ears, ‘latency
 head’) and body (whole body excluding the tail, ‘latency body’). The
 experimenter gently pushed the vole out of the tube if it did not exit
 after five minutes (maximum latency 300 s). In the open-field test, we
 observe the vole for five minutes and record the latency to enter the
 center section (‘latency center’ in s), the number of crossings into the
 center (with the whole body excluding the tail, ‘crossings’), the number
 of sections entered (with the whole body excluding the tail, ‘sections’),
 and the activity (running, jumping, scanning, grooming) instantaneously
 every 10 seconds (‘activity’). See supplements for details on the
 setup. Animal trapping, handling, and experiments were
 conducted according to German legislation on animal research under
 permission of the Landesamt für Umwelt, Verbraucherschutz und Gesundheit,
 Brandenburg (LUGV_RW7-4744/41+5#243052/2015 and
 LFU-N4-4730/11+10#120786/2021) and the Landesamt für Arbeitschutz,
 Verbraucherschutz und Gesundheit (2347-46-2018).
 Statistical analyses We first ran simple
 (generalized) linear mixed effects models ((G)LMMs) on single variables,
 using the ‘lmer’ and ‘glmmTMB’ functions (R-packages lme4: Bates
 <em>et al.</em> 2015; glmmTMB: Brooks <em>et
 al.</em> 2017, respectively). Latencies were log-transformed and
 proportions (‘activity’ and ‘sections’) were arcsin square root
 transformed. For behavioral variables, we controlled for the year, site,
 test occasion, if the animal was tested under natural or altered
 conditions, and their two-way interactions as fixed effects, and identity
 as a random effect. For variables related to microhabitat use, we
 controlled for the year and (only for models on the total dataset of
 microhabitat use) if the animal was trapped under natural or altered
 conditions as fixed effects, and individual identity and site as random
 effects. We included the site as a random rather than a fixed effect to
 control for heterogeneity due to differences across grids. The models on
 the use of microhabitats under both natural and altered conditions (all)
 and only under altered conditions were adjusted for zero inflation. Model
 selection was determined based on the lowest AIC values (&lt; Δ2 AIC
 considered equally good); interactions were excluded when they did not
 reduce the AIC value &gt; Δ2. Variables with a good model fit
 (assessed on normality of residuals, non-significant dispersion,
 non-significant deviation from uniformity, and homogeneity of variance)
 were kept for subsequent analyses. Secondly, to
 estimate repeatability and adjusted repeatability of single variables, we
 used the ´rpt´ function (rptR package: Nakagawa & Schielzeth 2013;
 Stoffel <em>et al.</em> 2017) using 1,000 simulations and
 1,000 permutations. We estimated raw repeatability and repeatability
 adjusted for predictors that had an effect in the simple univariate
 models. We report raw repeatability; see supplements for adjusted
 repeatabilities. Repeatable variables were used in subsequent bivariate
 mixed effects models. To evaluate the among-individual
 correlations between microhabitat use and personality, and the use of
 microhabitat under natural and altered conditions, we assessed the
 covariance between repeatable behavioral variables and microhabitat use
 (H1 and H2) by running six Bayesian bivariate mixed effects models (BMMs),
 using the ‘mcmcglmm’ function (´MCMCglmm´ package: Hadfield 2010, as in
 Dingemanse & Dochtermann 2013; Houslay & Wilson 2017), and
 calculated the pairwise among-individual correlation from the posterior
 distribution of covariances. All models included individual identity as a
 random effect and variable-specific fixed effects determined from simple
 univariate models. Models with microhabitat use included site as a random
 effect. We fixed the within-individual variance to zero because we did not
 measure all variables at the same time. We ran all models with two
 different weakly informative priors (as suggested by Hadfield 2010). Model
 results were robust against different priors (results not shown). The
 reported model results are based on a weak prior (V=diag(<em>number
 dependent variables</em>), nu=1.002). We used 250,000 iterations, a
 thinning interval of 100, and a burn-in of 50,000. Using the posterior
 distributions from the BMMs, we calculated repeatability for each
 dependent variable, pairwise among-individual correlations, and their
 credibility intervals. To assess whether there is a
 change in behavior with a change in the use of safe microhabitats under
 natural versus altered conditions (H3) we performed Spearman’s rank
 correlations due to restrictions by small sample sizes (n = 11 or 12
 depending on analysis). We correlated changes (i.e., delta) in repeatable
 behavior with changes in microhabitat use from natural to altered
 conditions. To assess the association between body
 condition and repeatable behavioral variables and microhabitat use under
 natural conditions, we ran (G)LMMs with the repeatable behavioral
 variables and microhabitat use as dependent variable, body condition and
 year as fixed effects, and individual identity and site (only for
 microhabitat use) as random effects. We included two-way interactions
 between body condition and year when they reduced the AIC value.
 We assessed population-level avoidance of mowed traps with a
 χ²-test, from the observed (use of traps under altered conditions) and
 expected (use of traps under natural conditions). Additionally, we ran a
 GLMM to assess if behavior affected reappearance (yes/no) under altered
 conditions with repeatable behavioral variables as fixed effects, and
 individual identity as a random effect. Lastly, we assessed if the
 percentage of traps used under natural conditions that were affected by
 mowing changed the use of vegetation cover after mowing by running a GLMM
 with year and percentage mowed traps as fixed effects and individual
 identity and site as random effects as well as adjusting for
 zero-inflation. In keeping with the ‘language of
 evidence’ (Muff <em>et al.</em> 2022) and Dingemanse
 <em>et al.</em> (2020) we interpreted an effect as “strongly
 supported” when zero was not included in the 95% CIs, “moderately or
 weakly supported” when the point estimate was skewed away from zero whilst
 the 95% CIs overlapped zero, and strong support for the absence of an
 effect when estimates were centered at zero. All analyses were carried out
 in R, Version 2023.06.0+421 (R Core Team 2023).
description: # Data and code from: A timid choice: Risk-taking behavior predicts
 individualized niche in a varying landscape of safety
 [https://doi.org/10.5061/dryad.0p2ngf29s](https://doi.org/10.5061/dryad.0p2ngf29s) ## Description of the data and file structure ### Study site (trapping & animals) The study was conducted in a meadow landscape with regular mowing practices between September and November 2019, and in October 2022 in the Havelland district, Brandenburg, Germany (52°44’00” N 12°12’44” E). The study site is on an island surrounded by the river Havel, and the smaller rivers Gülper Havel and Pirre. Each year we set up two grid systems (50 m x 50 m) for capturing voles with six by six Ugglan traps (Grahnab Sweden, special no. 2, with added shrew exit) evenly spaced out 10 m apart. One grid in 2022 had 4 x 6 traps due to space restrictions. We conducted capture-mark-recapture (CMR) for at least two weeks in each grid. Grids in 2019 were surrounded by an electric fence (90 cm high; mesh size 17.5 cm wide and 9 to 15 cm high) to restrict disturbance of traps by raccoons. Traps were activated around sunset and checked after sunrise (maximum of 8 hours between checking active traps). In 2019 traps were additionally activated and checked during the day. Traps were baited with oat flakes and apples or cucumber. In colder temperatures (&lt;5°C) traps were provided with wood shavings as nesting material to prevent hypothermia. We weighed, measured head width, and determined sex and age class based on body mass (≥17 g = mature/adult, &lt;17 g = juvenile/immature) of all common voles. Individuals were fur-marked with a unique pattern for individual identification. After this handling procedure, animals were released next to the trap they came from. Based on head width (mm) and body mass (g), we calculated body condition using the scaled body mass index (Peig & Green 2009). This index was calculated as an individual's body mass (in g) standardized to the mean head width to retain a body size measure of all individuals in the population with a population-specific allometric relationship. ### Vegetation We surveyed the natural microhabitat in an area of one meter around each trap by measuring the height (average maximum height of five measuring points) and estimating cover at ten centimeters above ground by counting the number of grasses/branches touching a 1-meter ruler. To experimentally make the habitat less safe, we mowed and removed the grass 10 m x 10 m around every other trap in the grid (except the first grid in 2019 where half of the home range of VHF-tracked individuals was mowed). Thereby, we intended to affect the home ranges of all individuals to some degree. The size of the mowed patches was chosen to resemble home range sizes of male common voles estimated from VHF-tracking (mean: 202 m2 Jacob & Hempel 2003); 3,830 to 5,950 m2: A. Schirmer, J. Eccard & M. Dammhahn, unpublished data). We decided against a more individual-tailored manipulation since this required work-intensive efforts of individual tracking and excluded smaller individuals. Based on a pilot VHF-telemetry tracking study in the study population, individuals' home ranges covered areas larger than 10 m x 10 m around each trap (F. Erixon, M. Gilmour, A. Schirmer, unpublished data). After mowing, the mowed patches had significantly lower vegetation height (&lt; 10 cm) and, thus, lower vegetation cover (Fig. SI 1). ### Behavioral assays Upon first capture, we assess among-individual differences in behaviors related to risk-taking of all adult common voles (n = 82) using a combined *emergence* and *open-field* test (Herde & Eccard 2013; Schirmer *et al.* 2019). Tests were conducted indoors and all handling and sample collection was done after testing. We reassessed recaptured animals after 5-17 days (median = 7). The setup consists of an opaque tube (11 cm diameter, 30 cm length) with a swing door connected to a round arena (120 cm diameter, 60 cm height) divided into 16 equal-sized sections: eight at the edge and eight at the exposed central region (84 cm diameter). First, the animal is taken from the trap into a can that is placed in the tub and left to habituate for one minute. The experimenter then opens the swing door to start the *emergence* test. Once the vole exits the tube the swing door is closed and the open-field test starts. The experiments were scored in real time in 2019 and recorded on camera and scored from videos in 2022. Between tests, the setup was cleaned with 70% ethanol. In the emergence test, we measured the latency (seconds) to exit the tube with its head (including the ears, ‘latency head’) and body (whole body excluding the tail, ‘latency body’). The experimenter gently pushed the vole out of the tube if it did not exit after five minutes (maximum latency 300 s). In the open-field test, we observe the vole for five minutes and record the latency to enter the center section (‘latency center’ in s), the number of crossings into the center (with the whole body excluding the tail, ‘crossings’), the number of sections entered (with the whole body excluding the tail, ‘sections’), and the activity (running, jumping, scanning, grooming) instantaneously every 10 seconds (‘activity’). See supplements for details on the setup. Animal trapping, handling, and experiments were conducted according to German legislation on animal research under permission of the Landesamt für Umwelt, Verbraucherschutz und Gesundheit, Brandenburg (LUGV_RW7-4744/41+5#243052/2015 and LFU-N4-4730/11+10#120786/2021) and the Landesamt für Arbeitschutz, Verbraucherschutz und Gesundheit (2347-46-2018). ### References in this README file Bates, D., Maechler, M., Bolker, B. & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. *J. Stat. Softw.*, 67, 1–48. Brooks, M.E., Kristensen, K., van Benthem, K.J., Magnusson, A., Berg, C.W., Nielsen, A., *et al.* (2017). glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. *R J.*, 9, 378--400. Dingemanse, N.J. & Dochtermann, N.A. (2013). Quantifying individual variation in behaviour: mixed-effect modelling approaches. *J. Anim. Ecol.*, 82, 39–54. Dingemanse, N.J., Moiron, M., Araya‐Ajoy, Y.G., Mouchet, A. & Abbey‐Lee, R.N. (2020). Individual variation in age‐dependent reproduction: Fast explorers live fast but senesce young? *J. Anim. Ecol.*, 89, 601–613. Hadfield, J.D. (2010). MCMC Methods for Multi-Response Generalized Linear Mixed Models: The **MCMCglmm** *R* Package. *J. Stat. Softw.*, 33, 1–22. Herde, A. & Eccard, J.A. (2013). Consistency in boldness, activity and exploration at different stages of life. *BMC Ecol.*, 13, 49. Houslay, T.M. & Wilson, A.J. (2017). Avoiding the misuse of BLUP in behavioural ecology. *Behav. Ecol.*, 28, 948–952. Jacob, J. & Hempel, N. (2003). Effects of farming practices on spatial behaviour of common voles. *J. Ethol.*, 21, 45–50. Muff, S., Nilsen, E.B., O’Hara, R.B. & Nater, C.R. (2022). Rewriting results sections in the language of evidence. *Trends Ecol. Evol.*, 37, 203–210. Nakagawa, S. & Schielzeth, H. (2013). A general and simple method for obtaining *R*2 from generalized linear mixed‐effects models. *Methods Ecol. Evol.*, 4, 133–142. Peig, J. & Green, A.J. (2009). New perspectives for estimating body condition from mass/length data: the scaled mass index as an alternative method. *Oikos*, 118, 1883–1891. R Core Team. (2023). R: A Language and Environment for Statistical Computing. Schirmer, A., Herde, A., Eccard, J.A. & Dammhahn, M. (2019). Individuals in space: personality-dependent space use, movement and microhabitat use facilitate individual spatial niche specialization. *Oecologia*, 189, 647–660. Stoffel, M.A., Nakagawa, S. & Schielzeth, H. (2017). rptR: repeatability estimation and variance decomposition by generalized linear mixed‐effects models. *Methods Ecol. Evol.*, 8, 1639–1644. ### Files and variables #### File: behavior.csv **Description:** Contains data on behavior assessed in the standardized behavioral assay combining an emergence and open-field test. Missing data value: NA ##### Variables * identity: identification of individual voles * year: the year of testing (2019 or 2022) * site: name of site where an individual was live-trapped and tested (Grid1.1, Grid2, GS1, GS3) * date: the date of when the combined emergence and open-field test was conducted * tested_after_mowing: defining whether or not the current test was conducted before or after alteration of the environment through mowing (0 = tested under natural conditions, 1 = tested under altered conditions) * sex: the sex of the individual tested (F = female, M = male, NA = missing data/sex not possible to assess) * test_occasion: defines what test occasion, i.e., first, second, third, etc., time tested (1 = first time tested, 2 = second time tested, 3 = third time tested, 4 = fourth time tested) * lat_head: latency head, the latency (in seconds) to emerge from the emergence test with the head (including the ears) (NA = missing data) * lat_body: latency body, the latency (in seconds) to emerge from the emergence test with the whole body (excluding the tail) (NA = missing data) * lat_center: latency center, the latency (in seconds) to cross into the center section of the open-field test with the whole body (excluding the tail) for the first time (NA = missing data) * no_sections: number of sections entered in the open field test with the whole body (excluding the tail) (NA = missing data) * active: number of 10 second instances spent active (running, jumping, scanning, grooming), (between 0-35) (NA = missing data) * inactive: number of 10 second instances inactive, (between 0-35) (NA = missing data) * crossings: number of crossings into the center section of the open field with the whole body (excluding the tail) (NA = missing data) #### File: body\_condition.csv **Description:** Contains head width (mm) and body mass (g) information, one measurement per individual. If multiple head width/body mass measurements have been taken for individuals, the first head width and closest body mass measurement are included in this data set. ##### Variables * identity: identification of individual voles * date: date of when the head width was measured * bm: the closest body mass measurement, typically taken the same day (grams) * hw: head width (mm) #### File: vegetation.csv **Description:** Contains data on vegetation height and cover measured at all traps in all sites under natural conditions, and who took the measurements.  ##### Variables * site: name of site where the vegetation was assessed (Grid1.1, Grid2, GS1, GS3) * trap: name of trap at which the vegetation was assessed * height: mean of five measurements of max vegetation height at 1 m x 1m around the trap (cm) * cover: the vegetation cover at the trap (number of branches touching a meter) * who: initials of the person who assessed the vegetation #### File: cmr\_altered.csv **Description:** Contains information gathered during capture-mark-recapture under altered conditions (after manipulating the microhabitat of predation risk through mowing). Missing data value: NA ##### Variables * identity: identification of individual voles * year: the year of testing (2019 or 2022) * site: name of site where an individual was live-trapped and tested (Grid1.1, Grid2, GS1, GS3) * date: date of when the trapping occasion was conducted * time: time of day of when traps were checked * trap: name of trap at which the individual was caught * capture_recapture: whether or not an individual was trapped for the first time ( = new) or had already been trapped before (= recapture) * sex: the sex of the individual tested (F = female, M = male, NA = missing data/sex not possible to assess/information not taken) * age: age of individual, juvenile (&lt;17 g) or adult (&gt;= 17 g) (NA = information not taken,* typically when the individual is recaptured and this information has already been taken*)  * bm: body mass measurement (g) taken (NA = information not taken,* typically when the individual is recaptured and this information has already been taken*)  * hw: head width measurement (mm) taken (NA = information not taken,* typically when the individual is recaptured and this information has already been taken*)  * trap_mowed: defines whether or not the particular trap captured had been mowed, i.e., altered the vegetation higher to create a risky microhabitat, or not (0 = not mowed, 1 = mowed) * height: mean of five measurements of max vegetation height 1 m x 1 m around the trap (cm) * cover: the vegetation cover at the trap (number of branches touching a meter) #### File: cmr\_natural.csv **Description:** Contains information gathered during capture-mark-recapture under natural conditions (before manipulating the microhabitat of predation risk through mowing). Missing data value: NA ##### Variables * identity: identification of individual voles * year: the year of testing (2019 or 2022) * site: name of site where an individual was live-trapped and tested (Grid1.1, Grid2, GS1, GS3) * date: date of when the trapping occasion was conducted * time: time of day of when the trapping occasion was conducted * trap: name of trap at which the individual was caught * capture_recapture: whether or not an individual was trapped for the first time ( = new) or had already been trapped before (= recapture) * sex: the sex of the individual tested (F = female, M = male, NA = missing data/not possible to assess sex/information not taken) * age: age of individual, juvenile (&lt;17 g) or adult (&gt;= 17 g) (NA = information not taken,* typically when the individual is recaptured and this information has already been taken*)  * bm: body mass measurement (g) taken (NA = information not taken,*typically when the individual is recaptured and this information has already been taken*)  * hw: head width measurement (mm) taken (NA = information not taken,*typically when the individual is recaptured and this information has already been taken*)  * trap_mowed: defines whether or not the particular trap captured had been mowed, i.e. altered the vegetation higher to create a risky microhabitat, or not (0 = not mowed, 1 = mowed) * height: mean of five measurements of max vegetation height 1 m x 1 m around the trap (cm) * cover: the vegetation cover at the trap (number of branches touching a meter) ## Code/software All analyses were carried out in Rstudio, Version 2023.06.0+421 (R Core Team 2023). All included data sets (csv files) are used for the analysis in the code. The following packages have been used to view and analyse the data, and plot the results and are required to run the code: tidyverse;  Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” \_Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 &lt;[https://doi.org/10.21105/joss.01686&gt;](https://doi.org/10.21105/joss.01686&gt;). here; Müller K (2020). \_here: A Simpler Way to Find Your Files_. R package version 1.0.1, [https://CRAN.R-project.org/package=here](https://CRAN.R-project.org/package=here). rptR; Stoffel, M. A., Nakagawa, S. and Schielzeth, H. (2017), rptR: repeatability estimation and variance decomposition by generalized linear mixed-effects models. Methods Ecol Evol, 8: 1639???1644. doi:10.1111/2041-210X.12797 lme4; Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01. glmmTMB; Mollie E. Brooks, Kasper Kristensen, Koen J. van Benthem, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Maechler and Benjamin M. Bolker (2017). glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal, 9(2), 378-400. doi: 10.32614/RJ-2017-066. lmerTest; Kuznetsova A, Brockhoff PB, Christensen RHB (2017). “lmerTest Package: Tests in Linear Mixed Effects Models.” \_Journal of Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13 [https://doi.org/10.18637/jss.v082.i13](https://doi.org/10.18637/jss.v082.i13). car; Fox J, Weisberg S (2019). \_An R Companion to Applied Regression_, Third edition. Sage, Thousand Oaks CA. &lt;[https://socialsciences.mcmaster.ca/jfox/Books/Companion/&gt;](https://socialsciences.mcmaster.ca/jfox/Books/Companion/&gt;). DHARMa; Hartig F (2022). \_DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. R package version 0.4.6, &lt;[https://CRAN.R-project.org/package=DHARMa&gt;](https://CRAN.R-project.org/package=DHARMa&gt;). ggplot2; H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016. lmtest;  Achim Zeileis, Torsten Hothorn (2002). Diagnostic Checking in Regression Relationships. R News 2(3), 7-10. URL [https://CRAN.R-project.org/doc/Rnews/](https://CRAN.R-project.org/doc/Rnews/) sjPlot; Lüdecke D (2023). \_sjPlot: Data Visualization for Statistics in Social Science_. R package version 2.8.15, &lt;[https://CRAN.R-project.org/package=sjPlot&gt;](https://CRAN.R-project.org/package=sjPlot&gt;). MCMCglmm; Jarrod D Hadfield (2010). MCMC Methods for Multi-Response Generalized Linear Mixed Models: The MCMCglmm R Package. Journal of Statistical Software, 33(2), 1-22. URL [https://www.jstatsoft.org/v33/i02/](https://www.jstatsoft.org/v33/i02/). bayestestR; Makowski, D., Ben-Shachar, M., & Lüdecke, D. (2019). bayestestR: Describing Effects and their Uncertainty, Existence and Significance within the Bayesian Framework. Journal of Open Source Software, 4(40), 1541. doi:10.21105/joss.01541