Statistical inference for stochastic simulation models--theory and application

Statistical models are the traditional choice to test scientific theories when observations, processes or boundary conditions are subject to stochasticity. Many important systems in ecology and biology, however, are difficult to capture with statistical models. Stochastic simulation models offer an alternative, but they were hitherto associated with a major disadvantage: their likelihood functions can usually not be calculated explicitly, and thus it is difficult to couple them to well-established statistical theory such as maximum likelihood and Bayesian statistics. A number of new methods, among them Approximate Bayesian Computing and Pattern-Oriented Modelling, bypass this limitation. These methods share three main principles: aggregation of simulated and observed data via summary statistics, likelihood approximation based on the summary statistics, and efficient sampling. We discuss principles as well as advantages and caveats of these methods, and demonstrate their potential for integrating stochastic simulation models into a unified framework for statistical modelling.     

A likelihood framework for inferring the evolution of geographic range on phylogenetic trees

Abstract At a time when historical biogeography appears to be again expanding its scope after a period of focusing primarily on discerning area relationships using cladograms, new inference methods are needed to bring more kinds of data to bear on questions about the geographic history of lineages. Here we describe a likelihood framework for inferring the evolution of geographic range on phylogenies that models lineage dispersal and local extinction in a set of discrete areas as stochastic events in continuous time. Unlike existing methods for estimating ancestral areas, such as dispersal‐vicariance analysis, this approach incorporates information on the timing of both lineage divergences and the availability of connections between areas (dispersal routes). Monte Carlo methods are used to estimate branch‐specific transition probabilities for geographic ranges, enabling the likelihood of the data (observed species distributions) to be evaluated for a given phylogeny and parameterized paleogeographic model. We demonstrate how the method can be used to address two biogeographic questions: What were the ancestral geographic ranges on a phylogenetic tree? How were those ancestral ranges affected by speciation and inherited by the daughter lineages at cladogenesis events? For illustration we use hypothetical examples and an analysis of a Northern Hemisphere plant clade (Cercis), comparing and contrasting inferences to those obtained from dispersal‐vicariance analysis. Although the particular model we implement is somewhat simplistic, the framework itself is flexible and could readily be modified to incorporate additional sources of information and also be extended to address other aspects of historical biogeography.     

Prediction and error in multi‐stage models for spread of aquatic non‐indigenous species

Aim Predictions of spread of non‐indigenous species allow for greater efficiency in managing invasions by targeting areas for preventative measures. The invasion sequence is a useful concept in predictions of spread, as it allows us to test hypotheses about the transport and establishment of propagules in novel habitats. Our aims are twofold: (1) to develop and validate multi‐stage invasion models for the introduced fishhook waterflea, Cercopagis pengoi, and (2) to assess how variability in the transport patterns of the propagules influences the accuracy and spatial extent for predictions of spread. Location New York State, USA. Methods We developed a two‐stage model for the spread of C. pengoi. First, we developed a stochastic gravity model for dispersal based on surveys of recreational boat traffic in New York State as a proxy for propagule pressure. We then modelled the probability of establishment based on predicted levels of propagule pressure and measures of lakes’ physicochemistry. In addition, we used Monte Carlo simulations based on the gravity model to propagate variability in boater traffic through the establishment model to assess how uncertainty in dispersal influenced predictions of spread. Results The amount recreationalists were willing to spend, lake area and population size of the city nearest to the destination lake were significant factors affecting boater traffic. In turn, boater traffic, lake area, specific conductance and turbidity were significant predictors of establishment. The inclusion of stochastic dispersal reduced the rate of false positives (i.e. incorrect prediction of an invasion) in detecting invasions at the upper 95% prediction interval for the probability of establishment. Main conclusions Combinations of measures of propagule pressure, habitat suitability and stochastic dispersal allow for the most accurate predictions of spread. Further, multi‐stage spread models may overestimate the extent of spread if stochasticity in early stages of the models is not considered.     

Directional asymmetry of long-distance dispersal and colonization could mislead reconstructions of biogeography

Aim Phylogenies are increasingly being used to attempt to answer biogeographical questions. However, a reliance on tree topology alone has emerged without consideration of earth processes or the biology of the organisms in question. Most ancestral‐state optimization methods have inherent problems, including failure to take account of asymmetry, such as unequal probabilities of losses and gains, and the lack of use of independent cost estimates. Here we discuss what we perceive as shortcomings in most current tree‐based biogeography interpretation methods and show that consideration of processes and their likelihoods can turn the conventional biogeographical interpretation on its head. Location Southern hemisphere focus but applicable world‐wide. Methods The logic of existing methods is reviewed with respect to their adequacy in modelling processes such as geographical mode of speciation and likelihood of dispersal, including directional bias. Published reconstructions of dispersal of three plant taxa between Australia and New Zealand were re‐analysed using standard parsimony and maximum likelihood (ML) methods with rate matrices to model expected asymmetry of dispersal. Results Few studies to date incorporate asymmetric dispersal rate matrices or question the simplistic assumption of equal costs. Even when they do, cost matrices typically are not derived independently of tree topology. Asymmetrical dispersal between Australia and New Zealand could be reconstructed using parsimony but not with ML. Main conclusions The inadequacy of current models has important consequences for our interpretation of southern hemisphere biogeography, particularly in relation to dispersal. For example, if repeated directional dispersals and colonization in the direction of prevailing winds have occurred, with intervening periods of speciation, then there is no need to infer dispersals against those winds. Failure to take account of directionality and other biases in reconstruction methods has implications beyond the simple misinterpretation of the biogeography of a taxonomic group, such as calibration of molecular clocks, the dating of vicariance events, and the prioritization of areas for conservation.     

SIMULATION‐BASED LIKELIHOOD APPROACH FOR EVOLUTIONARY MODELS OF PHENOTYPIC TRAITS ON PHYLOGENY

Phylogenetic comparative methods (PCMs) have been used to test evolutionary hypotheses at phenotypic levels. The evolutionary modes commonly included in PCMs are Brownian motion (genetic drift) and the Ornstein–Uhlenbeck process (stabilizing selection), whose likelihood functions are mathematically tractable. More complicated models of evolutionary modes, such as branch‐specific directional selection, have not been used because calculations of likelihood and parameter estimates in the maximum‐likelihood framework are not straightforward. To solve this problem, we introduced a population genetics framework into a PCM, and here, we present a flexible and comprehensive framework for estimating evolutionary parameters through simulation‐based likelihood computations. The method does not require analytic likelihood computations, and evolutionary models can be used as long as simulation is possible. Our approach has many advantages: it incorporates different evolutionary modes for phenotypes into phylogeny, it takes intraspecific variation into account, it evaluates full likelihood instead of using summary statistics, and it can be used to estimate ancestral traits. We present a successful application of the method to the evolution of brain size in primates. Our method can be easily implemented in more computationally effective frameworks such as approximate Bayesian computation (ABC), which will enhance the use of computationally intensive methods in the study of phenotypic evolution.     

Evolutionary processes driving spatial patterns of intraspecific genetic diversity in river ecosystems

Describing, understanding and predicting the spatial distribution of genetic diversity is a central issue in biological sciences. In river landscapes, it is generally predicted that neutral genetic diversity should increase downstream, but there have been few attempts to test and validate this assumption across taxonomic groups. Moreover, it is still unclear what are the evolutionary processes that may generate this apparent spatial pattern of diversity. Here, we quantitatively synthesized published results from diverse taxa living in river ecosystems, and we performed a meta‐analysis to show that a downstream increase in intraspecific genetic diversity (DIGD) actually constitutes a general spatial pattern of biodiversity that is repeatable across taxa. We further demonstrated that DIGD was stronger for strictly waterborne dispersing than for overland dispersing species. However, for a restricted data set focusing on fishes, there was no evidence that DIGD was related to particular species traits. We then searched for general processes underlying DIGD by simulating genetic data in dendritic‐like river systems. Simulations revealed that the three processes we considered (downstream‐biased dispersal, increase in habitat availability downstream and upstream‐directed colonization) might generate DIGD. Using random forest models, we identified from simulations a set of highly informative summary statistics allowing discriminating among the processes causing DIGD. Finally, combining these discriminant statistics and approximate Bayesian computations on a set of twelve empirical case studies, we hypothesized that DIGD were most likely due to the interaction of two of these three processes and that contrary to expectation, they were not solely caused by downstream‐biased dispersal.     

A diagnostic for Cox regression with discrete failure-time models

Changes in maximum likelihood parameter estimates due to deletion of individual observations are useful statistics, both for regression diagnostics and for computing robust estimates of covariance. For many likelihoods, including those in the exponential family, these delete-one statistics can be approximated analytically from a one-step Newton-Raphson iteration on the full maximum likelihood solution. But for general conditional likelihoods and the related Cox partial likelihood, the one-step method does not reduce to an analytic solution. For these likelihoods, an alternative analytic approximation that relies on an appropriately augmented design matrix has been proposed. In this paper, we extend the augmentation approach to explicitly deal with discrete failure-time models. In these models, an individual subject may contribute information at several time points, thereby appearing in multiple risk sets before eventually experiencing a failure or being censored. Our extension also allows the covariates to be time dependent. The new augmentation requires no additional computational resources while improving results.     

Metapopulation dynamics and biological invasions: a spatially explicit model applied to a freshwater snail

The spatial spread of invading organisms is a major contemporary concern. We focus here on invasions in inherently fragmented habitats, such as freshwater systems, and explore the usefulness of metapopulation models in this context. Maximum-likelihood methods allow the estimation of colonization and extinction rates, as functions of habitat patch sizes and positions, from time series of presence/absence data. This framework also provides confidence intervals of these estimates and hypotheses tests. We analyze a previously unpublished 12-year survey of the spread of the introduced snail Tarebia granifera in 47 Martinican rivers. Simple metapopulation models reproduce with reasonable accuracy several quantitative aspects of the invasion, including regional abundance, spatiotemporal structure, and site-by-site colonization dates. Sensitivity analysis reveals that the invasion sequence depended strongly on metapopulation size (number of sites) and spatial structure (distances among sites). The invasion history has also been accelerated by stochastic events, as illustrated by a large, central river that happened to be colonized very early and served as an invasion pool. Finally, we discuss the benefits of this approach for the understanding of invasions in fragmented landscapes.     

Analytical methods for predicting the behaviour of population models with general spatial interactions

Many biologists use population models that are spatial, stochastic and individual based. Analytical methods that describe the behaviour of these models approximately are attracting increasing interest as an alternative to expensive computer simulation. The methods can be employed for both prediction and fitting models to data. Recent work has extended existing (mean field) methods with the aim of accounting for the development of spatial correlations. A common feature is the use of closure approximations for truncating the set of evolution equations for summary statistics. We investigate an analytical approach for spatial and stochastic models where individuals interact according to a generic function of their distance; this extends previous methods for lattice models with interactions between close neighbours, such as the pair approximation. Our study also complements work by Bolker and Pacala (BP) [Theor. Pop. Biol. 52 (1997) 179; Am. Naturalist 153 (1999) 575]: it treats individuals as being spatially discrete (defined on a lattice) rather than as a continuous mass distribution; it tests the accuracy of different closure approximations over parameter space, including the additive moment closure (MC) used by BP and the Kirkwood approximation. The study is done in the context of an susceptible-infected-susceptible epidemic model with primary infection and with secondary infection represented by power-law interactions. MC is numerically unstable or inaccurate in parameter regions with low primary infection (or density-independent birth rates). A modified Kirkwood approximation gives stable and generally accurate transient and long-term solutions; we argue it can be applied to lattice and to continuous-space models as a substitute for MC. We derive a generalisation of the basic reproduction ratio, R(0), for spatial models.     

Parameter inference for stochastic biochemical models from perturbation experiments parallelised at the single cell level

Understanding and characterising biochemical processes inside single cells requires experimental platforms that allow one to perturb and observe the dynamics of such processes as well as computational methods to build and parameterise models from the collected data. Recent progress with experimental platforms and optogenetics has made it possible to expose each cell in an experiment to an individualised input and automatically record cellular responses over days with fine time resolution. However, methods to infer parameters of stochastic kinetic models from single-cell longitudinal data have generally been developed under the assumption that experimental data is sparse and that responses of cells to at most a few different input perturbations can be observed. Here, we investigate and compare different approaches for calculating parameter likelihoods of single-cell longitudinal data based on approximations of the chemical master equation (CME) with a particular focus on coupling the linear noise approximation (LNA) or moment closure methods to a Kalman filter. We show that, as long as cells are measured sufficiently frequently, coupling the LNA to a Kalman filter allows one to accurately approximate likelihoods and to infer model parameters from data even in cases where the LNA provides poor approximations of the CME. Furthermore, the computational cost of filtering-based iterative likelihood evaluation scales advantageously in the number of measurement times and different input perturbations and is thus ideally suited for data obtained from modern experimental platforms. To demonstrate the practical usefulness of these results, we perform an experiment in which single cells, equipped with an optogenetic gene expression system, are exposed to various different light-input sequences and measured at several hundred time points and use parameter inference based on iterative likelihood evaluation to parameterise a stochastic model of the system.     

Evidence for non-uniform dispersal in the biological invasions of two naturalized North American bird species

Most previous attempts to model the geographical range expansion of an invading species assume random dispersal of organisms through a homogeneous environment. These models result in a series of uniformly increasing circles radiating out from the centre of origin over time. Although these models often give reasonable fits to available data, they do not typically include mechanisms of dispersal. Alternatively, models that include assumptions of non‐random dispersal and a heterogeneous environment inevitably result in an anisotropic or jagged invasion front. This front will include propagules of pioneer individuals for the expanding species. Existing data from biological invasions reveal that the spatial structure of an invading species usually exhibits these propagules. Using population data gathered from the past century, we investigated the propagules of two North American invading bird species: the European starling (Sturnus vulgaris Linnaeus), and the house finch (Carpodacus mexicanus Müller), and found a correlation between propagule location and habitat quality. These results suggest that dispersing individuals seek out favourable habitat and remain there, thus introducing a possible mechanism for explaining non‐uniform dispersal during invasions. When combined with results from other studies, our results suggest that propagules provide starting points for future population expansion of an invading species.     

Accounting for spatial autocorrelation in null models of tree species association

A commonly used null model for species association among forest trees is a well‐mixed community (WMC). A WMC represents a non‐spatial, or spatially implicit, model, in which species form nearest‐neighbor pairs at a rate equal to the product of their community proportions. WMC models assume that the outcome of random dispersal and demographic processes is complete spatial randomness (CSR) in the species’ spatial distributions. Yet, stochastic dispersal processes often lead to spatial autocorrelation (SAC) in tree species densities, giving rise to clustering, segregation, and other nonrandom patterns. Although methods exist to account for SAC in spatially‐explicit models, its impact on non‐spatial models often remains unaccounted for. To investigate the potential for SAC to bias tests based upon non‐spatial models, we developed a spatially‐heterogeneous (SH) modelling approach that incorporates measured levels of SAC. Using the mapped locations of individuals in a tropical tree community, we tested the hypothesis that the identity of nearest‐neighbors represents a random draw from neighborhood species pools. Correlograms of Moran's I confirmed that, for 50 of 51 dominant species, stem density was significantly autocorrelated over distances ranging from 50 to 200 m. The observed patterns of SAC were consistent with dispersal limitation, with most species occurring in distinct patches. For nearly all of the 106 species in the community, the frequency of pairwise association was statistically indistinguishable from that projected by the null models. However, model comparisons revealed that non‐spatial models more strongly underestimated observed species‐pair frequencies, particularly for conspecific pairs. Overall, the CSR models projected more significant facilitative interactions than did SH models, yielding a more liberal test of niche differences. Our results underscore the importance of accounting for stochastic spatial processes in tests of association, regardless of whether spatial or non‐spatial models are employed.     

Using ABC and microsatellite data to detect multiple introductions of invasive species from a single source

The introduction of invasive species to new locations (that is, biological invasions) can have major impact on biodiversity, agriculture and public health. As such, determining the routes and modality of introductions with genetic data has become a fundamental goal in molecular ecology. To assist with this goal, new statistical methods and frameworks have been developed, such as approximate Bayesian computation (ABC) for inferring invasion history. Here, we present a model of invasion accounting for multiple introductions from a single source (MISS), a heretofore largely unexplored model. We simulate microsatellite data to evaluate the power of ABC to distinguish between single and multiple introductions from the same source, under a range of demographic parameters. We also apply ABC to microsatellite data from three invasions of bumblebee in New Zealand. In addition, we assess the performance of several methods of summary statistics selection. Our simulated results suggested good ability to distinguish between one- and two-wave models over much but not all of the parameter space tested, independent of summary statistics used. Globally, parameter estimation was good except for bottleneck timing. For one of the bumblebee species, we clearly rejected the MISS model, while for the other two we found inconclusive results. Since a second wave may provide genetic reinforcement to initial colonists, help relieve inbreeding among founders, or increase the hazard of the invasion, its detection may be crucial for managing invasions; we suggest that the MISS model could be considered as a potential model in future theoretical and empirical studies of invasions.     

Invasion speeds for structured populations in fluctuating environments

We live in a time where climate models predict future increases in environmental variability and biological invasions are becoming increasingly frequent. A key to developing effective responses to biological invasions in increasingly variable environments will be estimates of their rates of spatial spread and the associated uncertainty of these estimates. Using stochastic, stage-structured, integrodifference equation models, we show analytically that invasion speeds are asymptotically normally distributed with a variance that decreases in time. We apply our methods to a simple juvenile–adult model with stochastic variation in reproduction and an illustrative example with published data for the perennial herb, Calathea ovandensis. These examples buttressed by additional analysis reveal that increased variability in vital rates simultaneously slow down invasions yet generate greater uncertainty about rates of spatial spread. Moreover, while temporal autocorrelations in vital rates inflate variability in invasion speeds, the effect of these autocorrelations on the average invasion speed can be positive or negative depending on life history traits and how well vital rates “remember” the past.     

Saddle-Point Approximations,Integrodifference Equations,and Invasions

Invasion, the growth in numbers and spatial spread of a population over time, is a fundamental process in ecology. Governments and businesses expend vast sums to prevent and control invasions of pests and pestilences and to promote invasions of endangered species and biological control agents. Many mathematical models of biological invasions use nonlinear integrodifference equations to describe the growth and dispersal processes and to predict the speed of invasion fronts. Linear models have received less attention, perhaps because they are difficult to simulate for large times. In this paper, we use the saddle-point method, alias the method of steepest descent, to derive asymptotic approximations for the solutions of linear integrodifference equations. We work through five examples, for Gaussian, Laplace, and uniform dispersal kernels in one dimension and for asymmetric Gaussian and radially symmetric Laplace kernels in two dimensions. Our approximations are extremely close to the exact solutions, even for intermediate times. We also employ an empirical saddle-point approximation to predict densities using dispersal data. We use our approximations to examine the effects of censored dispersal data on estimates of invasion speed and population density.     

Testing environmental and genetic effects in the presence of spatial autocorrelation

Spatial autocorrelation is a well‐recognized concern for observational data in general, and more specifically for spatial data in ecology. Generalized linear mixed models (GLMMs) with spatially autocorrelated random effects are a potential general framework for handling these spatial correlations. However, as the result of statistical and practical issues, such GLMMs have been fitted through the undocumented use of procedures based on penalized quasi‐likelihood approximations (PQL), and under restrictive models of spatial correlation. Alternatively, they are often neglected in favor of simpler but more questionable approaches. In this work we aim to provide practical and validated means of inference under spatial GLMMs, that overcome these limitations. For this purpose, a new software is developed to fit spatial GLMMs. We use it to assess the performance of likelihood ratio tests for fixed effects under spatial autocorrelation, based on Laplace or PQL approximations of the likelihood. Expectedly, the Laplace approximation performs generally slightly better, although a variant of PQL was better in the binary case. We show that a previous implementation of PQL methods in the R language, glmmPQL, is not appropriate for such applications. Finally, we illustrate the efficiency of a bootstrap procedure for correcting the small sample bias of the tests, which applies also to non‐spatial models.     

Predicted range expansion of Chinese tallow tree (Triadica sebifera) in forestlands of the southern United States

Aim We present an integrated approach for predicting future range expansion of an invasive species (Chinese tallow tree) that incorporates statistical forecasting and analytical techniques within a spatially explicit, agent‐based, simulation framework. Location East Texas and Louisiana, USA. Methods We drew upon extensive field data from the US Forest Service and the US Geological Survey to calculate spread rate from 2003 to 2008 and to parameterize logistic regression models estimating habitat quality for Chinese tallow within individual habitat cells. We applied the regression analyses to represent population spread rate as a function of habitat quality, integrated this function into a logistic model representing local spread, and coupled this model with a dispersal model based on a lognormal kernel within the simulation framework. We simulated invasions beginning in 2003 based on several different dispersal velocities and compared the resulting spatial patterns to those observed in 2008 using cross Mantel’s tests. We then used the best dispersal velocity to predict range expansion to the year 2023. Results Chinese tallow invasion is more likely in low and flat areas adjacent to water bodies and roads, and less likely in mature forest stands and in pine plantations where artificial regeneration by planting seedlings is used. Forecasted invasions resulted in a distribution that extended from the Gulf Coast of Texas and Louisiana northward and westward as much as 300 km, representing approximately 1.58 million ha. Main conclusions Our new approach of calculating time series projections of annual range expansion should assist land managers and restoration practitioners plan proactive management strategies and treatments. Also, as field sampling continues on the national array of FIA plots, these new data can be incorporated easily into the present model, as well as being used to develop and/or improve models of other invasive plant species.     

Experiments with the Site Frequency Spectrum

Evaluating the likelihood function of parameters in highly-structured population genetic models from extant deoxyribonucleic acid (DNA) sequences is computationally prohibitive. In such cases, one may approximately infer the parameters from summary statistics of the data such as the site-frequency-spectrum (SFS) or its linear combinations. Such methods are known as approximate likelihood or Bayesian computations. Using a controlled lumped Markov chain and computational commutative algebraic methods, we compute the exact likelihood of the SFS and many classical linear combinations of it at a non-recombining locus that is neutrally evolving under the infinitely-many-sites mutation model. Using a partially ordered graph of coalescent experiments around the SFS, we provide a decision-theoretic framework for approximate sufficiency. We also extend a family of classical hypothesis tests of standard neutrality at a non-recombining locus based on the SFS to a more powerful version that conditions on the topological information provided by the SFS.     

Accounting for multi‐scale spatial autocorrelation improves performance of invasive species distribution modelling (iSDM)

Aim Analyses of species distributions are complicated by various origins of spatial autocorrelation (SAC) in biogeographical data. SAC may be particularly important for invasive species distribution models (iSDMs) because biological invasions are strongly influenced by dispersal and colonization processes that typically create highly structured distribution patterns. We examined the efficacy of using a multi‐scale framework to account for different origins of SAC, and compared non‐spatial models with models that accounted for SAC at multiple levels. Location We modelled the spatial distribution of an invasive forest pathogen, Phytophthora ramorum, in western USA. Methods We applied one conventional statistical method (generalized linear model, GLM) and one nonparametric technique (maximum entropy, Maxent) to a large dataset on P. ramorum occurrence (n = 3787) to develop four types of model that included environmental variables and that either ignored spatial context or incorporated it at a broad scale using trend surface analysis, a local scale using autocovariates, or multiple scales using spatial eigenvector mapping. We evaluated model accuracies and amounts of explained spatial structure, and examined the changes in predictive power of the environmental and spatial variables. Results Accounting for different scales of SAC significantly enhanced the predictive capability of iSDMs. Dramatic improvements were observed when fine‐scale SAC was included, suggesting that local range‐confining processes are important in P. ramorum spread. The importance of environmental variables was relatively consistent across all models, but the explanatory power decreased in spatial models for factors with strong spatial structure. While accounting for SAC reduced the amount of residual autocorrelation for GLM but not for Maxent, it still improved the performance of both approaches, supporting our hypothesis that dispersal and colonization processes are important factors to consider in distribution models of biological invasions. Main conclusions Spatial autocorrelation has become a paradigm in biogeography and ecological modelling. In addition to avoiding the violation of statistical assumptions, accounting for spatial patterns at multiple scales can enhance our understanding of dynamic processes that explain ecological mechanisms of invasion and improve the predictive performance of static iSDMs.