The lme4 R package can be used to fit generalized linear mixed models (GLMMs), which extend the class of linear mixed models (LMMs). The two main extensions provided by GLMMs are (1) allowing for the conditional distribution of the response given the random effects to be non-Gaussian (e.g. binomial, Poisson) and (2) allowing the conditional mean to be a nonlinear function of a linear combination of the fixed and random effect coefficients, via an inverse link function. The conditional mode of the random effects given the observed data, the variance-covariance matrix of the random effects, and the fixed effect parameters are determined using penalized iteratively reweighted least squares. We compute an approximation of the integral over the distributions of the conditional modes to compute the maximum likelihood estimate for a given set of parameters (by default we use the Laplace approximation or, alternatively, the more computationally expensive adaptive Gauss-Hermite quadrature). The package provides all the standard features available for GLMs in base R, including the standard set of accessor functions as well as the possibility of user-specified distributions (within the exponential dispersion family) and link functions.

Multivariate random effects with unstructured variance-covariance matrices of large dimensions, q, can be a major challenge to estimate. In this paper, we introduce a new implementation of a reduced-rank approach to fit large dimensional multivariate random effects by writing them as a linear combination of d < q latent variables. By adding reduced-rank functionality to the package glmmTMB, we enhance the mixed models available to include random effects of dimensions that were previously not possible. We apply the reduced-rank random effect to two examples, estimating a generalized latent variable model for multivariate abundance data and a random-slopes model.

Invasive species are a global problem with large ecological and economic costs. A better understanding of how invasive species populations change over time, how these species become integrated into ecosystems, and how their population demographics vary across different environments could help inform management priorities and shape control strategies. For 20 years (2002–2022), we have monitored round goby (Neogobius melanostomus) in Hamilton Harbour, Canada, an industrial harbour and an Area of Concern with high levels of contaminants. We sampled round goby across six sites that vary in contamination levels. We first quantified changes in round goby population demographics and morphology over a twenty-year period and second, we compared how abundance and other life history trajectories differ between sites of high and low contamination. Round goby abundance and body length both decreased over the study period. In contrast, body condition, gonadosomatic index (GSI), and the proportion of guarding parental males in the population increased over time. Over the many years of monitoring, there was no clear difference in round goby abundance between sites of high and low contamination, but individuals from sites of high contamination were smaller, had larger gonad investment, and higher hepatosomatic index compared to round goby from sites of low contamination. We also found there were fewer guarding parental males at sites of high contamination. Our results are valuable because they provide insights into how invasive species interact with different invaded habitats over the long-term. This information can help researchers and managers understand the effects of invasive species and develop strategies to predict, prevent, and manage them.

Phylogenetic generalised linear mixed models (PGLMMs) help ecologists to distinguish ecological drivers from other processes shaping evolutionary patterns, yet existing implementations are often limited in distributional scope or computational speed. We compare five R packages for fitting PGLMMs and highlight the new covariance structure \`propto\` in the general-purpose GLMM package glmmTMB. Simulations show that glmmTMB fits PGLMMs faster overall than brms, MCMCglmm, INLA, and phyr, while producing similar model estimates. We present the first practical application of glmmTMB for fitting phylogenetic random effects using likelihood-based models that accommodate repeated measures, demonstrated through case studies of evolutionary trait data. By improving both speed and flexibility, glmmTMB broadens access to PGLMM and supports deeper insights into trait evolution and diversification. ### Competing Interest Statement BB and MM are contributors and maintainers of the glmmTMB R package. Australian Research Council (ARC) Discovery Grant, DP230101248 Canada Excellence Research Chair, CERC-2022-00074

Reproductive accessory glands are organs involved in reproduction that do not directly produce or release gametes but can play crucial roles in securing reproductive success. In fishes, the 2 leading hypotheses about why accessory glands evolved are (1) in response to sperm competition, or (2) to facilitate parental care activities. Here, we investigate the evolutionary history of accessory glands and test these hypotheses by estimating quantitative differences in evolutionary rates. We found that accessory glands are present in 116 of the 607 sampled species of ray-finned fishes, representing 26/267 families. We estimated that accessory glands have arisen independently ~20 times and that these glands were gained 5.8 times faster in lineages with male parental care, compared to those without male care, supporting the hypothesis that they evolved to facilitate care. In contrast, group spawning, used as a proxy for sperm competition risk, seemed to select against the evolution of accessory glands, as lineages exhibiting group spawning gained accessory glands 3.9 times slower than those with pair spawning (though this failed to reach statistical significance). This study provides new insights into the evolutionary history of accessory glands in fishes and highlights the importance of parental care in shaping reproductive anatomy.

Yersinia pestis has spilled over from wild rodent reservoirs to commensal rodents and humans, causing three historically recorded pandemics. Depletion in the copy number of the plasmid-encoded virulence gene pla occurred in later-dated strains of the first and second pandemics, yet the biological relevance of the pla deletion has been difficult to test. We identified modern Y. pestis strains that independently acquired the same pla depletion as ancient strains and herein show that excision of pla from the multicopy pPCP1 plasmid is accompanied by the integration of a separate full pPCP1 harboring pla into the single-copy pCD1 plasmid, reducing pla dosage. Moreover, we demonstrate that this depletion decreases the mortality of mice in models of bubonic plague but not in the pneumonic and septicemic forms of the disease. We hypothesize that pla depletion may have been selectively advantageous in bubonic plague, owing to rodent fragmentation after pandemic-induced mortality.

Rabies spread by domestic dogs continues to cause tens of thousands of human deaths every year in low- and middle-income countries. Nevertheless rabies is often neglected, perhaps because it has already been eliminated from high-income countries through dog vaccination. Estimates of canine rabies’s intrinsic reproductive number ( ℛ ), a metric of disease spread, from a wide range of times and locations are relatively low (values < 2), with narrow confidence intervals. Given rabies’s persistence, this consistently low and narrow range of estimates is surprising. We combined incidence data from historical outbreaks of canine rabies from around the world with in-depth contact-tracing data from Tanzania to investigate initial growth rates ( r ), generation-interval distributions ( G ), and reproductive numbers ( ℛ ). We improved on earlier estimates by choosing outbreak windows algorithmically; fitting r using a more appropriate statistical method that accounts for decreases through time; and incorporating uncertainty from both r and G in our confidence intervals on ℛ . Our ℛ estimates are larger than previous estimates, with wider confidence intervals. These revised ℛ estimates suggest that a greater level of vaccination effort will be required to eliminate rabies than previously thought, but that the level of coverage required remains feasible. Our hybrid approach for estimating ℛ and its uncertainty is applicable to other disease systems where researchers estimate ℛ by combining data-based estimates of r and G .### Competing Interest StatementThe authors have declared no competing interest.

Information-theoretic (IT) and multi-model averaging (MMA) statistical approaches are widely used but suboptimal tools for pursuing a multifactorial approach (also known as the method of multiple working hypotheses) in ecology. (1) Conceptually, IT encourages ecologists to perform tests on sets of artificially simplified models. (2) MMA improves on IT model selection by implementing a simple form of shrinkage estimation (a way to make accurate predictions from a model with many parameters relative to the amount of data, by “shrinking” parameter estimates toward zero). However, other shrinkage estimators such as penalized regression or Bayesian hierarchical models with regularizing priors are more computationally efficient and better supported theoretically. (3) In general, the procedures for extracting confidence intervals from MMA are overconfident, providing overly narrow intervals. If researchers want to use limited data sets to accurately estimate the strength of multiple competing ecological processes along with reliable confidence intervals, the current best approach is to use full (maximal) statistical models (possibly with Bayesian priors) after making principled, a priori decisions about model complexity.

Background: Canadian notifiable disease surveillance provides data on the incidence of communicable diseases, dating back to the late 19th century. The Public Health Agency of Canada offers summaries of these data (from 1924-2022) through an online portal. These summaries provide historical context for Canadian health researchers, but lack information on intra-annual and inter-provincial patterns. Sub-annual and sub-national data can be found in published documents, or requested from archives of government agencies, but are only available in typewritten or handwritten hard copies. We digitized and collated these data sources to create a resource for epidemiology and public health. Methods: We manually entered data from scans of hard copies into spreadsheets resembling the originals, facilitating accurate transcription through easier cross-checking. We developed open-source pipelines to harmonize these spreadsheets into CSV files that blend data across sources. Results: We assembled and processed 1,631,380 incidence values from 1903-2021. Focusing on sub-annual and sub-national data and removing redundancy yielded 934,009 weekly, monthly, or quarterly incidence values broken down by province/territory, containing 139 diseases. We give two examples of sub-annual and sub-national patterns: strong annual cycles of poliomyelitis that peaked simultaneously across provinces, and spatially heterogeneous resurgence of whooping cough in the 1990s. Interpretation: Canada's history of infectious disease surveillance has produced a detailed record of sub-annual and sub-national disease incidence patterns that remains largely unexplored. This important record is now available as the Canadian Disease Incidence Dataset (CANDID), hosted on a publicly accessible website along with the pipelines used to create it and scans of the original sources. ### Competing Interest Statement The authors have declared no competing interest. ### Funding Statement This study was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) via an Emerging Infectious Disease Modelling (EIDM) grant to the Canadian Network for Modelling Infectious Diseases (CANMOD). ### Author Declarations I confirm all relevant ethical guidelines have been followed, and any necessary IRB and/or ethics committee approvals have been obtained. Yes The details of the IRB/oversight body that provided approval or exemption for the research described are given below: The Research Ethics Board of McMaster University gave ethical approval for this work. I confirm that all necessary patient/participant consent has been obtained and the appropriate institutional forms have been archived, and that any patient/participant/sample identifiers included were not known to anyone (e.g., hospital staff, patients or participants themselves) outside the research group so cannot be used to identify individuals. Yes I understand that all clinical trials and any other prospective interventional studies must be registered with an ICMJE-approved registry, such as ClinicalTrials.gov. I confirm that any such study reported in the manuscript has been registered and the trial registration ID is provided (note: if posting a prospective study registered retrospectively, please provide a statement in the trial ID field explaining why the study was not registered in advance). Yes I have followed all appropriate research reporting guidelines, such as any relevant EQUATOR Network research reporting checklist(s) and other pertinent material, if applicable. Yes All data produced are available online at https://github.com/canmod/iidda

State variables such as abundance and occurrence of species are central to many questions in ecology and conservation, but our ability to detect and enumerate species is imperfect and often varies across space and time. Accounting for imperfect and variable detection is important for obtaining unbiased estimates of state variables. Here, I investigate whether closed spatial capture-recapture (SCR) and single season occupancy models are robust to ignoring temporal variation in detection probability. Ignoring temporal variation allows collapsing detection data across repeated sampling occasions, speeding up computations, which can be important when analyzing large datasets with complex models. I simulated data under different scenarios of temporal and spatio-temporal variation in detection, analyzed data with the data-generating model and an alternative model ignoring temporal variation in detection, and compared estimates between these two models with respect to relative bias, coefficient of variation (CV) and relative root mean squared error (RMSE). SCR model estimates of abundance, the density-covariate coefficient β and the movement-related scale parameter of the detection function σ were robust to ignoring temporal variation in detection, with relative bias, CV and RMSE of the two models generally being within 4% of each other. An SCR case study for brown tree snakes showed identical estimates of density and σ under models accounting for or ignoring temporal variation in detection. Occupancy model estimates of the occupancy-covariate coefficient β and average occupancy were also largely robust to ignoring temporal variation in detection, and differences in occupancy predictions were mostly <<0.1. But there was a slight tendency for bias in β under the alternative model to increase when detection varied more strongly over time. Thus, when temporal variation in detection is extreme, it may be necessary to model that variation to avoid bias in parameter estimates in occupancy models. An occupancy case study for ten bird species with a more complex model structure showed considerable differences in occupancy parameter estimates under models accounting for or ignoring temporal variation in detection; but estimates and predictions from the latter were always within 95% confidence intervals of the former. There are cases where we cannot or may not want to ignore temporal variation in detection: a behavioral response to detection and certain SCR observation models do not allow collapsing data across sampling occasions; and temporal variation in detection may be informative of species phenology/behavior or for future study planning. But this study shows that it can be safely ignored under a range of conditions when analyzing SCR or occupancy data.

Fred Brauer was an eminent mathematician who studied dynamical systems, especially differential equations. He made many contributions to mathematical epidemiology, a field that is strongly connected to data, but he always chose to avoid data analysis. Nevertheless, he recognized that fitting models to data is usually necessary when attempting to apply infectious disease transmission models to real public health problems. He was curious to know how one goes about fitting dynamical models to data, and why it can be hard. Initially in response to Fred’s questions, we developed a user-friendly R package, fitode, that facilitates fitting ordinary differential equations to observed time series. Here, we use this package to provide a brief tutorial introduction to fitting compartmental epidemic models to a single observed time series. We assume that, like Fred, the reader is familiar with dynamical systems from a mathematical perspective, but has limited experience with statistical methodology or optimization techniques.

Taking a representative sample to determine prevalence of variables such as disease or vaccination in a population presents challenges, especially when little is known about the population. Several methods have been proposed for second stage cluster sampling. They include random sampling in small areas (the approach used in several international surveys), random walks within a specified geographic area, and using a grid superimposed on a map. We constructed 50 virtual populations with varying characteristics, such as overall prevalence of disease and variability of population density across towns. Each population comprised about a million people spread over 300 towns. We applied ten sampling methods to each. In 1,000 simulations, with different sample sizes per cluster, we estimated the prevalence of disease and the relative risk of disease given an exposure and calculated the Root Mean Squared Error (RMSE) of these estimates. We compared the sampling methods using the RMSEs. In our simulations a grid method was the best statistically in the great majority of circumstances. It showed less susceptibility to clustering effects, likely because it sampled over a much wider area than the other methods. We discuss the findings in relation to practical sampling issues.

We present an approach to computing the probability of epidemic “burnout,” i.e., the probability that a newly emergent pathogen will go extinct after a major epidemic. Our analysis is based on the standard stochastic formulation of the Susceptible-Infectious-Removed (SIR) epidemic model including host demography (births and deaths) and corresponds to the standard SIR ordinary differential equations (ODEs) in the infinite population limit. Exploiting a boundary layer approximation to the ODEs and a birth-death process approximation to the stochastic dynamics within the boundary layer, we derive convenient, fully analytical approximations for the burnout probability. We demonstrate—by comparing with computationally demanding individual-based stochastic simulations and with semi-analytical approximations derived previously—that our fully analytical approximations are highly accurate for biologically plausible parameters. We show that the probability of burnout always decreases with increased mean infectious period. However, for typical biological parameters, there is a relevant local minimum in the probability of persistence as a function of the basic reproduction number R 0 . For the shortest infectious periods, persistence is least likely if R 0 ≈ 2.57 ; for longer infectious periods, the minimum point decreases to R 0 ≈ 2 . For typical acute immunizing infections in human populations of realistic size, our analysis of the SIR model shows that burnout is almost certain in a well-mixed population, implying that susceptible recruitment through births is insufficient on its own to explain disease persistence.

Productivity is strongly associated with terrestrial species richness patterns, although the mechanisms underpinning such patterns have long been debated. Despite considerable consumption of primary productivity by fire, its influence on global diversity has received relatively little study. Here we examine the sensitivity of terrestrial vertebrate biodiversity (amphibians, birds and mammals) to fire, while accounting for other drivers. We analyse global data on terrestrial vertebrate richness, net primary productivity, fire occurrence (fraction of productivity consumed) and additional influences unrelated to productivity (i.e., historical phylogenetic and area effects) on species richness. For birds, fire is associated with higher diversity, rivalling the effects of productivity on richness, and for mammals, fire's positive association with diversity is even stronger than productivity; for amphibians, in contrast, there are few clear associations. Our findings suggest an underappreciated role for fire in the generation of animal species richness and the conservation of global biodiversity.

Compartmental models are valuable tools for investigating infectious diseases. Researchers building such models typically begin with a simple structure where compartments correspond to individuals with different epidemiological statuses, e.g., the classic SIR model which splits the population into susceptible, infected, and recovered compartments. However, as more information about a specific pathogen is discovered, or as a means to investigate the effects of heterogeneities, it becomes useful to stratify models further -- for example by age, geographic location, or pathogen strain. The operation of constructing stratified compartmental models from a pair of simpler models resembles the Cartesian product used in graph theory, but several key differences complicate matters. In this article we give explicit mathematical definitions for several so-called ``model products'' and provide examples where each is suitable. We also provide examples of model stratification where no existing model product will generate the desired result.

The management of forest pests relies on an accurate understanding of the species’ phenology. Thermal performance curves (TPCs) have traditionally been used to model insect phenology. Many such models have been proposed and fitted to data from both wild and laboratory-reared populations. Using Hamiltonian Monte Carlo for estimation, we implement and fit an individual-level, Bayesian hierarchical model of insect development to the observed larval stage durations of a population reared in a laboratory at constant temperatures. This hierarchical model handles interval censoring and temperature transfers between two constant temperatures during rearing. It also incorporates individual variation, quadratic variation in development rates across insects’ larval stages, and “flexibility” parameters that allow for deviations from a parametric TPC. Using a Bayesian method ensures a proper propagation of parameter uncertainty into predictions and provides insights into the model at hand. The model is applied to a population of eastern spruce budworm ( Choristoneura fumiferana ) reared at 7 constant temperatures. Resulting posterior distributions can be incorporated into a workflow that provides prediction intervals for the timing of life stages under different temperature regimes. We provide a basic example for the spruce budworm using a year of hourly temperature data from Timmins, Ontario, Canada. Supplementary materials accompanying this paper appear on-line.

The Cox proportional hazards model is commonly used in evaluating risk factors in cancer survival data. The model assumes an additive, linear relationship between the risk factors and the log hazard. However, this assumption may be too simplistic. Further, failure to take time-varying covariates into account, if present, may lower prediction accuracy. In this retrospective, population-based, prognostic study of data from patients diagnosed with cancer from 2008 to 2015 in Ontario, Canada, we applied machine learning-based time-to-event prediction methods and compared their predictive performance in two sets of analyses: (1) yearly-cohort-based time-invariant and (2) fully time-varying covariates analysis. Machine learning-based methods-gradient boosting model (gbm), random survival forest (rsf), elastic net (enet), lasso and ridge-were compared to the traditional Cox proportional hazards (coxph) model and the prior study which used the yearly-cohort-based time-invariant analysis. Using Harrell's C index as our primary measure, we found that using both machine learning techniques and incorporating time-dependent covariates can improve predictive performance. Gradient boosting machine showed the best performance on test data in both time-invariant and time-varying covariates analysis.

Testing individuals for pathogens can affect the spread of epidemics. Understanding how individual-level processes of sampling and reporting test results can affect community- or population-level spread is a dynamical modeling question. The effect of testing processes on epidemic dynamics depends on factors underlying implementation, particularly testing intensity and on whom testing is focused. Here, we use a simple model to explore how the individual-level effects of testing might directly impact population-level spread. Our model development was motivated by the COVID-19 epidemic, but has generic epidemiological and testing structures. To the classic SIR framework we have added a per capita testing intensity, and compartment-specific testing weights, which can be adjusted to reflect different testing emphases -- surveillance, diagnosis, or control. We derive an analytic expression for the relative reduction in the basic reproductive number due to testing, test-reporting and related isolation behaviours. Intensive testing and fast test reporting are expected to be beneficial at the community level because they can provide a rapid assessment of the situation, identify hot spots, and may enable rapid contact-tracing. Direct effects of fast testing at the individual level are less clear, and may depend on how individuals' behaviour is affected by testing information. Our simple model shows that under some circumstances both increased testing intensity and faster test reporting can reduce the effectiveness of control, and allows us to explore the conditions under which this occurs. Conversely, we find that focusing testing on infected individuals always acts to increase effectiveness of control.