Bayesian inference in a stochastic predator-prey system Sara - PowerPoint PPT Presentation

Bayesian inference in a stochastic predator-prey system Sara Pasquali and Fabrizio Ruggeri CNR-IMATI Via Bassini 15 - Milano - Italy Gianni Gilioli Universtity of Reggio Calabria Brixen, July 16 - 20, 2007 – p. 1/40

Discrete observations Parameter estimation from time series of predator-prey abundance ⇒ parametric inference of discretely observed diffusion processes. Prakasa-Rao B.L.S. (1999), Sørensen H. (2004). Large number of observations ⇒ consistency and asymptotic normality. Biology ⇒ few observations ⇒ Bayesian approach, MCMC framework. Eraker B. (2001), Elerian O. et al. (2001), Golightly A. and Wilkinson D.J. (2005). Brixen, July 16 - 20, 2007 – p. 2/40

Outline Description of the problem Deterministic model Stochastic model Estimation Simulated and field data Open problems Brixen, July 16 - 20, 2007 – p. 3/40

Biological considerations Study of functional response (i.e. individual predation rate as function of prey abundance) in an acarine predator-prey system Analysis of functional response usually performed in laboratory in a very simple experimental setup A problem is to understand if experiments performed in laboratory may represent the behavior of functional response in natural systems Brixen, July 16 - 20, 2007 – p. 4/40

Biological considerations The problem of extension to natural systems has to take into account scale problems: the behaviour of an individual in a small environment is not the same as in a large environment characteristics of environment organization: the spatial arrangement of the host plant may affect prey and predator behaviors degree of artificiality introduced by the experimental setup It is particularly useful to consider these aspects in the study of the functional response. Important in biological control where the use of predator-prey models may assist decision making. Brixen, July 16 - 20, 2007 – p. 5/40

Ecological assumptions Closed system (no immigration and/or emigration) ⇒ only local dynamics considered Functional response unaffected by variations in abiotic conditions (temperature, humidity, etc.) Interest in a single cycle of the population and not in the long term behavior (i.e. in the prompt control of the prey by the predator and not in the persistence of the system) Knowledge of all biodemographic parameters characterizing prey and predator populations Brixen, July 16 - 20, 2007 – p. 6/40

Deterministic model � dx t = [ rx t G ( x t ) − y t F ( x t , y t ; q )] dt x (0) = x 0 dy t = [ cy t F ( x t , y t ; q ) − uy t ] dt y (0) = y 0 x t , y t normalised biomass of prey and predator r = specific growth rate of the prey c = specific production rate of the predator u = specific loss rate of the predator q = efficiency of the predation process G ( x ) = growth of the prey in absence of predators F ( x, y ; q ) = functional response of the predator to the prey abundance Brixen, July 16 - 20, 2007 – p. 7/40

Lotka-Volterra system G ( x ) = 1 − x ; F ( x, y ; q ) = qx � dx t = [ rx t (1 − x t ) − qx t y t ] dt ⇒ dy t = [ cqx t y t − uy t ] dt Lotka-Volterra system modified by a logistic growth. Main limit: no saturation of the predator when the ingested prey increases; Main advantage: the model is simple and limits the number of parameters to be taken into account. Brixen, July 16 - 20, 2007 – p. 8/40

Demographic stochasticity q subject to noise and dependent on time ⇒ q t = q 0 + σξ t σ positive constant q 0 unknown parameter to be estimated ξ t a Gaussian white noise process ⇒ (early) stochastic model (1)  dx t = [ rx t (1 − x t ) − q 0 x t y t ] dt − σx t y t dw  t  (1)  dy t = [ cq 0 x t y t − uy t ] dt + cσx t y t dw  t w (1) t : Wiener process Brixen, July 16 - 20, 2007 – p. 9/40

Environmental stochasticity Environmental stochasticity affects prey and predator abundance (environmental variables fluctuation) w (2) t : Wiener process independent of w (1) t ε and ρ positive parameters ⇒ stochastic model  dx t = [ rx t (1 − x t ) − q 0 x t y t ] dt − σx t y t dw (1) + εx t dw (2)  t t dy t = [ cq 0 x t y t − uy t ] dt + cσx t y t dw (1) + ρy t dw (2)  t t Additive noise can be interpreted as sampling error affecting population abundance estimates. Solutions not necessarily in the compact [0 , 1] × [0 , 1] Brixen, July 16 - 20, 2007 – p. 10/40

Stochastic model Function χ ( z ) continuously differentiable and Lipschitz equal to 1 in the compact [ θ, 1 − θ ] decreasing in (1 − θ, + ∞ ) increasing in ( −∞ , θ ) lim z →−∞ χ ( z ) = lim z → + ∞ χ ( z ) = 0 χ (0) = χ (1) = θ ⇒ (final) stochastic model  dx t = [ rx t (1 − x t ) − q 0 x t y t ] χ ( x t ) dt − σx t y t χ ( x t ) dw (1) + εx t χ ( x t ) dw (2)  t t dy t = [ cq 0 x t y t − uy t ] χ ( y t ) dt + cσx t y t χ ( y t ) dw (1) + ηy t χ ( y t ) dw (2)  t t Brixen, July 16 - 20, 2007 – p. 11/40

Stochastic model ⇒ (Bivariate) diffusion process dX t = µ ( X t , q 0 ) dt + β ( X t ) dW t , X 0 = x 0 , t ≥ 0 µ ( X t , q 0 ) : drift coefficient � � [ rx t (1 − x t ) − q 0 x t y t ] χ ( x t ) µ ( X t , q 0 ) = [ cq 0 x t y t − uy t ] χ ( y t ) β ( X t ) : diffusion coefficient � � − σx t y t χ ( x t ) εx t χ ( x t ) β ( X t ) = cσx t y t χ ( y t ) ηy t χ ( y t ) q 0 , σ, ε, η unknown parameters Brixen, July 16 - 20, 2007 – p. 12/40

Likelihood � T � − 1 dX t + µ T ( X t , q 0 ) β ( X t ) β T ( X t ) � L ( q 0 ) = exp {− 0 � T � − 1 µ ( X t , q 0 ) dt } µ T ( X t , q 0 ) β ( X t ) β T ( X t ) � +1 / 2 0 Score function (i.e. derivative of log L ( q 0 ) w.r.t. q 0 ) � T � − 1 dX t + µ T ( X t , q 0 ) β ( X t ) β T ( X t ) � S T ( q 0 ) = − ˙ 0 � T � − 1 µ ( X t , q 0 ) dt µ T ( X t , q 0 ) β ( X t ) β T ( X t ) � +1 / 2 ˙ 0 Brixen, July 16 - 20, 2007 – p. 13/40

MLE µ ( X t ; q 0 ) = a ( X t ) q 0 + b ( X t ) ˆ X = ( X 0 , X 1 , ..., X p ) observations at times t 0 , t 1 , ..., t p ∆ i = t i − t i − 1 Discretized score function p − a T ( X i − 1 ) β ( X i − 1 ) β T ( X i − 1 ) � − 1 · [ X i − X i − 1 − µ ( X i − 1 )∆ i ] � � S N ( q 0 ) = i =1 p � − ρ ( x i − x i − 1 ) q 0 ,p = 1 1 � MLE: ˆ x i − 1 χ ( x i − 1 ) + T ρy i − 1 + cεx i − 1 i =1 + rρ (1 − x i − 1 )∆ i + ε ( y i − y i − 1 ) � y i − 1 χ ( y i − 1 ) + εu ∆ i Brixen, July 16 - 20, 2007 – p. 14/40

Bayesian estimation Posterior distribution p � � q 0 | ˆ � ∝ π ( q 0 ) f ( X i | X i − 1 , q 0 ) π X i =1 prior distribution π ( q 0 ) 1 2 · � � − 1 � β ( X i − 1 ) β T ( X i − 1 ) � f ( X i | X i − 1 , q 0 ) ∝ � � � � � 2 [ X i − X i − 1 − µ ( X i − 1 , q 0 )∆ i ] T − 1 · exp � − 1 [ X i − X i − 1 − µ ( X i − 1 , q 0 )∆ i ] � ∆ i β ( X i − 1 ) β T ( X i − 1 ) � · Brixen, July 16 - 20, 2007 – p. 15/40

Bayesian estimation � σ 2 I µ X + µ I σ 2 , σ 2 X σ 2 � µ I , σ 2 � � Prior N ⇒ posterior N , X I σ 2 I + σ 2 σ 2 I + σ 2 I X X X = σ 2 q 0 ,p and σ 2 with µ X = ˆ h � � i e t p σ 2 I µ X + µ I σ 2 ⇒ Bayesian estimator , X σ 2 I + σ 2 X i.e. weighted average of MLE and prior mean Improper prior π ( q 0 ) ∝ I (0 , ∞ ) ( q 0 ) ( q 0 − µ X ) 2 1 − 1 2 σ 2 ⇒ posterior √ I (0 , ∞ ) ( q 0 ) X µX 2 πσ 2 1 − Φ − X σX Gamma prior (later) Brixen, July 16 - 20, 2007 – p. 16/40

MCMC M latent data generated between two consecutive observations Matrix of all data � � x ( t 0 ) x ∗ ( t 1 ) x ∗ ( t M ) x ( t M +1 ) x ∗ ( t n − 1 ) x ( t n ) ... ... Y = y ( t 0 ) y ∗ ( t 1 ) y ∗ ( t M ) y ( t M +1 ) y ∗ ( t n − 1 ) y ( t n ) ... ... n = ( p − 1) M + p total number of observations X i = ( x ( t i ) , y ( t i )) the datum sampled at time t i i = ( x ∗ ( t i ) , y ∗ ( t i )) the latent datum at t i X ∗ Y i both a real and a latent datum at time t i Brixen, July 16 - 20, 2007 – p. 17/40

Generation of latent data 1st step: generate latent data through linear interpolation s -th iteration: we generate the latent datum X ∗ i from the conditional distribution π ( X ∗ i | Y i − 1 , Y i +1 ; q 0 ) where Y i − 1 is obtained at iteration s and Y i +1 at iteration s − 1 . π ( X ∗ i | Y i − 1 , Y i +1 ; q 0 ) ∝ P ( Y i +1 | X ∗ i ; q 0 ) P ( X ∗ i | Y i − 1 ; q 0 ) ∝ � − 1 ∆ i β ( Y i − 1 ) β T ( Y i − 1 ) � − 1 i − Y i − 1 − µ ( Y i − 1 ; q 0 )∆ i ] T � 2 [ X ∗ exp � · [ X ∗ i − Y i − 1 − µ ( Y i − 1 ; q 0 )∆ i ] � − 1 i ) β T ( X ∗ � − 1 i ; q 0 )∆ i +1 ] T � · exp 2 [ Y i +1 − X ∗ i − µ ( X ∗ ∆ i +1 β ( X ∗ i ) · [ Y i +1 − X ∗ i − µ ( X ∗ i ; q 0 )∆ i +1 ] } Brixen, July 16 - 20, 2007 – p. 18/40

Block size sampling Update latent observations in block of length m (Elerian et al., 2001) � � X ∗ ( k,m ) = X ∗ k , X ∗ k +1 , ..., X ∗ k + m − 1 k + m − 1 � � � � � ( k,m ) | Y k − 1 , Y k + m ; q 0 ∝ j | Y j − 1 , Y k + m ; q 0 X ∗ X ∗ f π j = k � � where π j | Y j − 1 , Y k + m ; q 0 depends on 2 observations: X ∗ one before the datum and the other after the block. Brixen, July 16 - 20, 2007 – p. 19/40