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

▶
bayesian inference in a stochastic predator prey system
SMART_READER_LITE
LIVE PREVIEW

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

Jun 08, 2023 350 likes •770 views

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

slide-1
SLIDE 1

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

slide-2
SLIDE 2

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

slide-3
SLIDE 3

Outline

Description of the problem Deterministic model Stochastic model Estimation Simulated and field data Open problems

Brixen, July 16 - 20, 2007 – p. 3/40

slide-4
SLIDE 4

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

slide-5
SLIDE 5

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

slide-6
SLIDE 6

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

slide-7
SLIDE 7

Deterministic model

  • dxt = [rxtG (xt) − ytF (xt, yt; q)] dt

x(0) = x0 dyt = [cytF (xt, yt; q) − uyt] dt y(0) = y0 xt, yt 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

slide-8
SLIDE 8

Lotka-Volterra system

G(x) = 1 − x

;

F(x, y; q) = qx ⇒

  • dxt = [rxt (1 − xt) − qxtyt] dt

dyt = [cqxtyt − uyt] 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

slide-9
SLIDE 9

Demographic stochasticity

q subject to noise and dependent on time ⇒ qt = q0 + σξt σ positive constant q0 unknown parameter to be estimated ξt a Gaussian white noise process ⇒ (early) stochastic model      dxt = [rxt (1 − xt) − q0xtyt] dt − σxtytdw

(1) t

dyt = [cq0xtyt − uyt] dt + cσxtytdw

(1) t

w(1)

t : Wiener process

Brixen, July 16 - 20, 2007 – p. 9/40

slide-10
SLIDE 10

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    dxt = [rxt(1 − xt) − q0xtyt]dt − σxtytdw(1)

t

+ εxtdw(2)

t

dyt = [cq0xtyt − uyt]dt + cσxtytdw(1)

t

+ ρytdw(2)

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

slide-11
SLIDE 11

Stochastic model

Function χ(z) continuously differentiable and Lipschitz equal to 1 in the compact [θ, 1 − θ] decreasing in (1 − θ, +∞) increasing in (−∞, θ)

limz→−∞ χ(z) = limz→+∞ χ(z) = 0 χ(0) = χ(1) = θ ⇒ (final) stochastic model

   dxt = [rxt(1 − xt) − q0xtyt]χ(xt)dt − σxtytχ(xt)dw(1)

t

+ εxtχ(xt)dw(2)

t

dyt = [cq0xtyt − uyt]χ(yt)dt + cσxtytχ(yt)dw(1)

t

+ ηytχ(yt)dw(2)

t

Brixen, July 16 - 20, 2007 – p. 11/40

slide-12
SLIDE 12

Stochastic model

⇒ (Bivariate) diffusion process dXt = µ (Xt, q0) dt + β (Xt) dWt, X0 = x0, t ≥ 0 µ (Xt, q0): drift coefficient µ (Xt, q0) =

  • [rxt(1 − xt) − q0xtyt]χ(xt)

[cq0xtyt − uyt] χ(yt)

  • β (Xt): diffusion coefficient

β (Xt) =

  • −σxtytχ(xt)

εxtχ(xt) cσxtytχ(yt) ηytχ(yt)

  • q0, σ, ε, η unknown parameters

Brixen, July 16 - 20, 2007 – p. 12/40

slide-13
SLIDE 13

Likelihood

L(q0) = exp{− T µT (Xt, q0)

  • β(Xt)βT (Xt)

−1 dXt+ +1/2 T µT (Xt, q0)

  • β(Xt)βT (Xt)

−1 µ (Xt, q0) dt}

Score function (i.e. derivative of log L(q0) w.r.t. q0)

ST (q0) = − T ˙ µT (Xt, q0)

  • β(Xt)βT (Xt)

−1 dXt+ +1/2 T ˙ µT (Xt, q0)

  • β(Xt)βT (Xt)

−1 µ (Xt, q0) dt

Brixen, July 16 - 20, 2007 – p. 13/40

slide-14
SLIDE 14

MLE

µ(Xt; q0) = a(Xt)q0 + b(Xt) ˆ X = (X0, X1, ..., Xp) observations at times t0, t1, ..., tp ∆i = ti − ti−1

Discretized score function

SN(q0) =

p

  • i=1

−aT (Xi−1)

  • β(Xi−1)βT (Xi−1)

−1·[Xi − Xi−1 − µ(Xi−1)∆i] MLE: ˆ q0,p = 1 T

p

  • i=1

1 ρyi−1 + cεxi−1

  • −ρ (xi − xi−1)

xi−1χ(xi−1) + +rρ(1 − xi−1)∆i + ε(yi − yi−1) yi−1χ(yi−1) + εu∆i

  • Brixen, July 16 - 20, 2007 – p. 14/40
slide-15
SLIDE 15

Bayesian estimation

Posterior distribution

π

  • q0| ˆ

X

  • ∝ π(q0)

p

  • i=1

f (Xi|Xi−1, q0)

prior distribution π(q0)

f (Xi|Xi−1, q0) ∝

  • β (Xi−1) βT (Xi−1)

−1

  • 1

2 ·

· exp

  • −1

2 [Xi − Xi−1 − µ(Xi−1, q0)∆i]T

·

  • ∆i β (Xi−1) βT (Xi−1)

−1 [Xi − Xi−1 − µ(Xi−1, q0)∆i]

  • Brixen, July 16 - 20, 2007 – p. 15/40
slide-16
SLIDE 16

Bayesian estimation

Prior N

  • µI, σ2

I

  • ⇒ posterior

N

  • σ2

IµX+µIσ2 X

σ2

I+σ2 X

, σ2

Xσ2 I

σ2

I+σ2 X

  • ,

with µX = ˆ

q0,p and σ2

X = σ2 tp

⇒ Bayesian estimator

σ2

IµX+µIσ2 X

σ2

I+σ2 X

, i.e. weighted average of MLE and prior mean Improper prior π (q0) ∝ I(0,∞)(q0)

⇒ posterior

1

√

2πσ2

X

h

1−Φ

  • −

µX σX

  • ie

−

1 2σ2 X

(q0−µX)2

I(0,∞)(q0)

Gamma prior (later)

Brixen, July 16 - 20, 2007 – p. 16/40

slide-17
SLIDE 17

MCMC

M latent data generated between two consecutive

  • bservations

Matrix of all data

Y =

  • x(t0)

x∗(t1) ... x∗(tM) x(tM+1) ... y(t0) y∗(t1) ... y∗(tM) y(tM+1) ... x∗(tn−1) x(tn) y∗(tn−1) y(tn)

  • n = (p − 1)M + p total number of observations

Xi = (x(ti), y(ti)) the datum sampled at time ti X∗

i = (x∗(ti), y∗(ti)) the latent datum at ti

Yi both a real and a latent datum at time ti

Brixen, July 16 - 20, 2007 – p. 17/40

slide-18
SLIDE 18

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 |Yi−1, Yi+1; q0) where Yi−1 is

  • btained at iteration s and Yi+1 at iteration s − 1.

π (X∗

i |Yi−1, Yi+1; q0) ∝ P (Yi+1|X∗ i ; q0) P (X∗ i |Yi−1; q0) ∝

exp

  • −1

2 [X∗

i − Yi−1 − µ(Yi−1; q0)∆i]T

∆iβ(Yi−1)βT (Yi−1) −1 · [X∗

i − Yi−1 − µ(Yi−1; q0)∆i]

  • ·exp
  • −1

2 [Yi+1 − X∗

i − µ(X∗ i ; q0)∆i+1]T

∆i+1β(X∗

i )βT (X∗ i )

−1 · [Yi+1 − X∗

i − µ(X∗ i ; q0)∆i+1]}

Brixen, July 16 - 20, 2007 – p. 18/40

slide-19
SLIDE 19

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

  • f
  • X∗

(k,m)|Yk−1, Yk+m; q0

  • ∝

k+m−1

  • j=k

π

  • X∗

j |Yj−1, Yk+m; q0

  • where π
  • X∗

j |Yj−1, Yk+m; q0

  • depends on 2 observations:
  • ne before the datum and the other after the block.

Brixen, July 16 - 20, 2007 – p. 19/40

slide-20
SLIDE 20

Block size sampling

π

  • X∗

j |Yj−1, Yk+m; q0

  • ∝ P
  • Yk+m|X∗

j ; q0

  • P
  • X∗

j |Yj−1; q0

  • ∝
  • β (Yj−1) βT (Yj−1)

−1

  • 1

2

  • β
  • X∗

j

  • βT

X∗

j

−1

  • 1

2

· exp

  • −1

2

  • X∗

j − Yj−1 − µ(Yj−1; q0)∆j

T ∆j β(Yj−1)βT (Yj−1) −1 ·

  • X∗

j − Yj−1 − µ(Yj−1; q0)∆j

  • · exp
  • −1

2

  • Yk+m − X∗

j − µ(X∗ j ; q0) (tk+m − tj)

T ·

  • (tk+m − tj) β(X∗

j )βT (X∗ j )

−1 ·

  • Yk+m − X∗

j − µ(X∗ j ; q0) (tk+m − tj)

  • Brixen, July 16 - 20, 2007 – p. 20/40
slide-21
SLIDE 21

M-H algorithm

f

  • X∗

(k,m)|Yk−1, Yk+m; q0

  • is not of a standard form,

M-H algorithm to simulate from f

  • X∗

(k,m)|Yk−1, Yk+m; q0

  • Proposal:

g

  • X∗

(k,m)|Yk−1, Yk+m; q0

  • ∝

k+m−1

  • j=k

N 1 2 (Yj−1 + Yk+m) , 1 2∆j · β(Yj−1)βT (Yj−1)

  • ·I(0,∞)×(0,∞)(X∗

j )

truncated bivariate normal distribution (truncation of the normal distribution proposed by Eraker (2001) and then used by Golightly and Wilkinson (2005)). Alternative choice based on the modified Brownian bridge construction of Durham and Gallant (2002).

Brixen, July 16 - 20, 2007 – p. 21/40

slide-22
SLIDE 22

M-H algorithm

X∗(s)

(k,m): value of X∗ (k,m) at the sth iteration.

Calculate X∗

(k,m) at the (s + 1)th iteration.

w: candidate vector generated from the proposal density.

Acceptance probability:

α(X∗

(k,m), w|Yk−1, Yk+m; q0) =

min   1, f (w|Yk−1, Yk+m; q0) g

  • X∗(s)

(k,m)|Yk−1, Yk+m; q0

  • f
  • X∗(s)

(k,m)|Yk−1, Yk+m; q0

  • g (w|Yk−1, Yk+m; q0)

  

We generate a uniform random number a – if a ≤ α(X∗

(k,m), w|Yk−1, Yk+m; q0) ⇒ X∗(s+1) (k,m)

= w

– otherwise X∗(s+1)

(k,m)

= X∗(s)

(k,m)

Brixen, July 16 - 20, 2007 – p. 22/40

slide-23
SLIDE 23

Simulation of q0

After simulation of Y ⇒ simulation of q0 Gamma prior on q0 ⇒ posterior has not a standard form

⇒ M-H algorithm to generate q0 from π(q0|Y ) ∝ π(q0)

n

  • i=1

π(Yi|Yi−1, q0) ⇒ gamma proposal

Improper prior on q0 ⇒ posterior is truncated normal

Brixen, July 16 - 20, 2007 – p. 23/40

slide-24
SLIDE 24

Application

Functional response parameter estimation in one of the most important system relevant to biological control problems: prey mite Tetranychus urticae and predator mite Phytoseiulus persimilis Population dynamics

   dxt = [rxt(1 − xt) − q0xtyt]χ(xt)dt − σxtytχ(xt)dw(1)

t

+ εxtχ(xt)dw(2)

t

dyt = [cq0xtyt − uyt]χ(yt)dt + cσxtytχ(yt)dw(1)

t

+ ηytχ(yt)dw(2)

t

x(0) = 0.1 y(0) = 0.007 r = 0.11 c = 0.35 u = 0.09

(Buffoni and Gilioli, 2003)

Brixen, July 16 - 20, 2007 – p. 24/40

slide-25
SLIDE 25

Application

Two datasets extensive survey of 8 separated local predator-prey dynamics (used to obtain estimates for σ, ε, η, q0 through LS method)

σ = 0.321 ε = 0.079 η = 0.106 q0 = 2.4767

intensive sampling of 12 separated local predator-prey dynamics (used to improve, through MCMC methods, the estimate of q0) Two cases simulated data (to test method performance) field data

Brixen, July 16 - 20, 2007 – p. 25/40

slide-26
SLIDE 26

Simulated data

Initial conditions: x0 = 0.1

y0 = 0.007

Total time: T = 100 days

q0 = 1.5

10 data simulated from the system MLE: 1.7357 Posterior median without latent data: 1.7347 No differences between MLE and Bayesian estimate.

Brixen, July 16 - 20, 2007 – p. 26/40

slide-27
SLIDE 27

No latent data

1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 2 4 6 8 10 12 14 q0 posterior density

Posterior density with median 1.7347, obtained using a gamma prior

Brixen, July 16 - 20, 2007 – p. 27/40

slide-28
SLIDE 28

Latent data

1 1.5 2 2.5 5 10 15 density (a) M = 1 1 1.5 2 2.5 5 10 15 (b) M = 5 1 1.5 2 2.5 5 10 15 density (c) M = 10 1 1.5 2 2.5 5 10 15 (d) M = 15 1 1.5 2 2.5 5 10 15 q0 density (e) M = 20 1 1.5 2 2.5 5 10 15 q0 (f) M = 30

MCMC with 100000 simulations, burn-in of 10000 iterations. Case of gamma prior for q0.

Brixen, July 16 - 20, 2007 – p. 28/40

slide-29
SLIDE 29

MCMC estimates

Latent data 1 5 10 15 20 30 posterior median 1.6358 1.5828 1.5729 1.5689 1.5678 1.5659

Improper prior leads to similar posterior MLE underestimates maximum values of prey and predators and display advanced cycles

Brixen, July 16 - 20, 2007 – p. 29/40

slide-30
SLIDE 30

Prey and predator dynamics

10 20 30 40 50 60 70 80 90 100 0.1 0.2 0.3 0.4 0.5 prey 10 20 30 40 50 60 70 80 90 100 0.05 0.1 time predator classical estimate Bayesian estimate simulated data

Initial conditions x0 = 0.1, y0 = 0.007. Classical trajectories: mean over 500 simulations. Bayesian setting: mean over 500 simulations and 100 values of q0.

Brixen, July 16 - 20, 2007 – p. 30/40

slide-31
SLIDE 31

Field data

10 20 30 40 50 60 70 80 90 100 0.2 0.4 0.6 0.8 prey 10 20 30 40 50 60 70 80 90 100 0.05 0.1 0.15 0.2 time predator

Data collected at 13 different time in a one-hectare strawberry crop in Ispica (Ragusa, Italy).

Brixen, July 16 - 20, 2007 – p. 31/40

slide-32
SLIDE 32

Field data

Predator and prey average abundance over 12 separated local dynamics - intensive sampling Only first cycle of the prey (48 days) no interest in predator persistence

⇒ reject data 8-13

rarity conditions strongly affect functional response

⇒ reject data 1

MLE: 2.6218 Posterior median: 2.6204

Brixen, July 16 - 20, 2007 – p. 32/40

slide-33
SLIDE 33

Latent data

1 1.5 2 2.5 3 5 10 density (a) M = 1 1 1.5 2 2.5 3 5 10 (b) M = 5 1 1.5 2 2.5 3 5 10 density (c) M = 10 1 1.5 2 2.5 3 5 10 (d) M = 15 1 1.5 2 2.5 3 5 10 q0 density (e) M = 20 1 1.5 2 2.5 3 5 10 q0 (f) M = 30

MCMC with 100000 simulations, burn-in of 10000 iterations. Case of gamma prior for q0.

Brixen, July 16 - 20, 2007 – p. 33/40

slide-34
SLIDE 34

Indexes

Smax: maximum size of the population ⇒ related to prey population impact (damage) on plants Tmax: time to reach maximum size ⇒ related to capability of the predator to restrain prey

population exponential growth

T0.5 (or T0.1): time to halve (or reduce to one tenth) Smax ⇒ related to capability of the predator to cause a rapid

decrease of the prey population

I: integral of the population up to T0.1 ⇒ measurement of population pressure on the resource

Similar definitions for predator population dynamics

Brixen, July 16 - 20, 2007 – p. 34/40

slide-35
SLIDE 35

Indexes: prey

Latent median of Smax tmax t0.5 t0.1 I Residual data

  • post. distr.

(days) (days) (days) 1 2.0386 0.6059 22.93 30.61 39.17 14.0505 0.3693 2 1.8732 0.6155 23.72 31.69 41.49 14.8856 0.3428 5 1.7332 0.6213 24.07 32.80 44.23 15.7425 0.3213 6 1.7161 0.6221 24.07 32.96 44.62 15.8585 0.3190 7 1.7038 0.6227 24.40 33.06 44.88 15.9436 0.3174 8 1.6942 0.6232 24.40 33.16 45.11 16.0110 0.3163 9 1.6869 0.6236 24.40 33.22 45.31 16.0628 0.3154 10 1.6809 0.6239 24.40 33.29 45.44 16.1058 0.3147 15 1.6639 0.6248 24.40 33.45 45.86 16.2296 0.3129 20 1.6554 0.6252 24.40 33.55 46.09 16.2926 0.3120 30 1.6467 0.6256 24.40 33.61 46.32 16.3579 0.3112

  • Class. est.

2.6218 0.6007 21.89 28.03 33.87 12.0858 0.4475 Field data 0.6393 21 37.8 48.58 16.0472

Brixen, July 16 - 20, 2007 – p. 35/40

slide-36
SLIDE 36

Indexes: predator

Latent median of Smax tmax t0.5 t0.1 I Residual data

  • post. distr.

(days) (days) (days) 1 2.0386 0.1685 34.69 47.27 70.13 3.4652 0.2031 2 1.8732 0.1673 34.57 49.03 74.74 3.5383 0.1998 5 1.7332 0.1665 36.75 50.83 82.00 3.5974 0.1984 6 1.7161 0.1664 36.85 51.06 82.00 3.6042 0.1983 7 1.7038 0.1663 37.24 51.22 82.00 3.6090 0.1983 8 1.6942 0.1663 37.24 51.38 82.00 3.6127 0.1983 9 1.6869 0.1662 37.24 51.48 82.00 3.6154 0.1983 10 1.6809 0.1662 37.24 51.55 82.00 3.6176 0.1983 15 1.6639 0.1661 37.24 51.74 82.00 3.6239 0.1983 20 1.6554 0.1660 37.24 51.87 82.00 3.6269 0.1984 30 1.6467 0.1659 36.57 51.97 82.00 3.6300 0.1984

  • Class. est.

2.6218 0.1738 31.36 43.02 62.72 3.2338 0.2217 Field data 0.1763 42 50.84 56.84 3.0145

Brixen, July 16 - 20, 2007 – p. 36/40

slide-37
SLIDE 37

Prey and predator dynamics

10 20 30 40 50 60 70 80 90 100 0.2 0.4 0.6 0.8 prey 10 20 30 40 50 60 70 80 90 100 0.05 0.1 0.15 0.2 time predator classical estimate Bayesian estimate simulated data

Initial conditions x0 = 0.1, y0 = 0.007. Classical trajectories: mean over 500 simulations. Bayesian setting: mean over 500 simulations and 100 values of q0.

Brixen, July 16 - 20, 2007 – p. 37/40

slide-38
SLIDE 38

Open problems

Estimate also parameter σ, η, ε with the proposed method Consider different forms of functional response (e.g. Ivlev) More data

Brixen, July 16 - 20, 2007 – p. 38/40

slide-39
SLIDE 39

Bibliography - biology

  • G. Buffoni, G. Di Cola, J. Baumgartener, V. Maurer (1995),

A mathematical model of trophic interactions in an acarine predator-prey system, Journal of Biological System, 3 (2), 303-312.

  • G. Buffoni, G. Gilioli (2003), A lumped parameter model for

acarine predator-prey population interactions, Ecological Modelling, 170, 155-171.

Brixen, July 16 - 20, 2007 – p. 39/40

slide-40
SLIDE 40

Bibliography - MCMC

  • O. Elerian, S. Chib, N. Shephard (2001), Likelihood

inference for discretely observed nonlinear diffusions, Econometrica, 69(4), 959-993.

  • B. Eraker (2001), MCMC analysis of diffusion models with

application to finance, Journal of Business & Economic Statistics, 19(2), 177-191.

  • A. Golightly, D.J. Wilkinson (2005), On Bayesian inference

for stochastic kinetic models using diffusion approximations, Biometrics, 61(3), 781-788. G.B. Durham, A.R. Gallant (2002), Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes, Journal of Business & Economic Statistics, 20(3), 297-316.

Brixen, July 16 - 20, 2007 – p. 40/40