Print bookPrint book

Model Simulation - Part I

Site: LifeWatch ERIC Training Platform
Course: Spatially explicit metapopulation model - Incidence Function Model
Book: Model Simulation - Part I
Printed by: Guest user
Date: Monday, 27 July 2026, 12:18 PM

Description

Model Simulation (I)

Table of contents

1. Model Simulation

Before discussing model parameters, we introduce a few examples showing how to simulate metapopulation dynamics following an IFM.

Regardless of the specific IFM model (the basic IFM, the “time dependent” or any of the two models with rescue effect), dynamics simulation (i.e., presence/absence at subsequent times) can be obtained in the following way. For a given parameters set, start the simulation from an initial vector of states (either 0 or 1) of the n local populations. For each time t calculate Si(t), Ci(t) and Ei(t) according to the model-specific formulas, and generate a random number U according to a uniform distribution on (0,1). Finally, determine the population state of patch i at time t+1, xi(t+1),  by comparing Ei(t) and Ci(t) to U :

\(x_i(t+1)= \begin{cases} 1,& if\: x_i(t) = 0 \: \: \& \:\:U \le C_i(t)\\ 0,& if\: x_i(t) = 0 \: \: \& \:\:U \gt C_i(t)\\ 0,& if\: x_i(t) = 1 \: \: \& \:\:U \le E_i(t)\\ 1,& if\: x_i(t) = \: \: \& \:\:U \gt E_i(t)\\ \end{cases} \)

Example 7.1. BASIC IFM.

We will consider, as usual, the reference metapopulation of Example 4.1.

We show how to simulate the dynamics of a basic IFM (see Unit 4). Vectors p and p1 are here used as “historical data” to determine the (constant) colonization probability, according to Eqs. (4.5) and (5.1).

 

DATA<-read.table("REFMETAPOP.txt",header=TRUE)

attach(DATA)

 

# parameter values

alpha<-1

A0<-0.0004

xhat<-0.05

yhat<-0.005

 

# extinction probabilities

E<-pmin(A0^xhat/A^xhat,1)

 

# connectivity and colonization probabilities, based on the mean          
# observed occupancy time, pm

d<-dist(cbind(x.crd,y.crd))/1000

edis<-as.matrix(exp(-alpha*d))

edis<-sweep(edis,2,A,"*")

pm<-(p+p1)/2

S<-rowSums(edis[,pm>0])

C<-S^2/(S^2+yhat^2)

 

# We introduce the function sim1: the argument of sim1 is any            
# presence/absence vector of length equal to the metapopulation size    
# (100).

# The key element of sim1 is the function runif.  runif(1) generates a  
# number between 0 and 1 (U)according to the uniform distribution.      
# runif(k) generates k numbers.

 

sim1<-function(o){</p> <p class="formule">o<- o > 0 </p> <p class="formule">if(any(o)){</p> <p class="formule">cond<-ifelse(o,E,C)</p> <p class="formule">o<-ifelse(runif(length(o))<cond,!o,o)</p> <p class="formule">}

as.numeric(o)

}

# the main instruction in the definition of sim1 is                     
# ifelse(runif(length(o))<cond,!o,o) that compares each of the 100      
# values “U” generated by runif to E if the corresponding element of o  
# is equal to 1 and to C otherwise.

 

# One simulation of 10 subsequent years of metapopulation dynamics,     
# starting from p1 (see Remark 3.2, Unit 3): we start with a null matrix
# sim; the initial state p1 is assigned to the first column, then by the
# command for a new column, i.e. metapopulation state, is generated from
# the previous one.

#

nt=10;

sim<-matrix(0,nrow=length(p1),ncol=nt+1)

sim[,1]<-p1

sim[1:10,]

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]

 [1,]    1    0    0    0    0    0    0    0    0     0     0

 [2,]    0    0    0    0    0    0    0    0    0     0     0

 [3,]    1    0    0    0    0    0    0    0    0     0     0

 [4,]    0    0    0    0    0    0    0    0    0     0     0

 [5,]    0    0    0    0    0    0    0    0    0     0     0

 [6,]    0    0    0    0    0    0    0    0    0     0     0

 [7,]    0    0    0    0    0    0    0    0    0     0     0

 [8,]    1    0    0    0    0    0    0    0    0     0     0

 [9,]    1    0    0    0    0    0    0    0    0     0     0

[10,]    1    0    0    0    0    0    0    0    0     0     0

 

for(t in 1:nt) sim[ ,t+1]<-sim1(sim[ ,t])

sim[1:10,]

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]

 [1,]    1    0    1    1    0    1    0    1    0     1     0

 [2,]    0    1    0    1    0    1    0    1    0     1     0

 [3,]    1    0    1    0    1    0    1    0    1     0     1

 [4,]    0    1    0    1    0    1    0    1    0     1     0

 [5,]    0    1    1    0    1    0    1    0    1     1     1

 [6,]    0    1    0    1    1    1    0    1    0     1     0

 [7,]    0    1    0    1    0    0    0    1    1     0     1

 [8,]    1    0    1    0    1    0    1    1    1     0     1

 [9,]    1    0    1    1    1    1    1    0    1     1     0

[10,]    1    0    1    0    1    0    1    0    1     0     1

 

# Note that if you run these commands on your computer, this

# output and the following picture as well will differ as runif      

# generates different values at each run.

 

# number of occupied sites in each year (i.e., column of sim)

plot(colSums(sim),type="o",lwd=2,xlab="t",ylab="n. occupied sites")

Figure 7.1 shows the effect of considering rescue effect in the case of Example 7.1 by comparing the yearly number of occupied sites, under the two models. This figure has been obtained with the following instructions:

 

plot(colSums(sim),type="o",lwd=2,xlab="t",ylab="n. occupied sites", ylim=c(0,100))

# ylim sets the vertical plot range

# simRE is the matrix of simulations obtained by including rescue     # effect (see Exercise 7.1)

points(colSums(simRE),col="red", type="o", lwd=2,xlab="t",ylab="n. occupied sites", ylim=c(0,100))

# points works like plot, however it allows to add the new graphic to # the previous one

text(3,98,labels="RESCUE EFFECT")

text(3,48,labels="NO RESCUE EFFECT")

Figure 7.9 – Number of occupied sites in two simulations of the basic IFM with and without rescue effect.

 

and clearly shows that the number of occupied sites strongly increases by the rescue effect.

Example 7.2. “Time-dependent” IFM.
We now consider the model of Unit 6, without rescue effect (Eq. (6.2)). As both connectivity and colonization probability are not constant, the computation of S and C cannot be carried out in advance, but should be included in the function sim1.

 

sim1TD<-function(o){</p> <p class="formule">o<- o > 0 </p> <p class="formule">if(any(o)){</p> <p class="formule">S<-rowSums(edis[,o,drop=FALSE])</p> <p class="formule">C<-S^2/(S^2+yhat^2)</p> <p class="formule">cond<-ifelse(o,E,C)</p> <p class="formule">o<-ifelse(runif(length(o))<cond,!o,o)</p> <p class="formule">}

as.numeric(o)

}

 

The simulated numbers of occupied sites are compared in the following picture. The chain starts at p1. Instructions are as before.

Pictures like Figure 7.1 or the one in Example 7.2 are not sufficient to clearly show the differences among models, as the visual effect is determined by the result of the specific simulation. For instance, the picture of Example 7.2 could suggest that the time-dependent model adds more variability to the yearly number of occupied sites than the basic IFM (excluding rescue effect in both cases). However, another run of sim1TD yields Figure 7.2, where the black and the blue lines are very similar. 

Figure 7.2 - Number of occupied sites in one simulation of the basic IFM with and without rescue effect compared to one simulation of the time-dependent IFM without rescue effect.

 

More precise information is obtained by averaging over a large number of simulations. By repeating simulations for several times, information on the distributions of the vectors Xt  at different times t can also be obtained, as shown in Example 7.3. We discuss this latter issue in the following Section 7.1 and go back to model comparison in Unit 8.

2. A simple approximation of the distribution of Xt

From  Unit 3 we know that P(Xt+1  = x) can be obtained as a sum of all the possible paths X1 = xk,1 , … , Xt = xk,t, Xt+1 = x,  where xk,j  is a vector of 0s and 1s.  Let assume that the initial metapopulation state is certain, that is π assigns the unit mass to a given vector (see Remark 3.2, Unit 3). Then, what the simulation function “sim1” does is just to provide one of such a paths, the matrix “sim”, thanks to the transition probabilities “E” and “C” (the j-th column of “sim” is xk,j ).  Therefore, it can be intuitively understood that by repeating this sampling procedure many times to obtain several “likely” paths and then summing up over all the simulations we can obtain an approximation of P(Xt+1  = x).

Example 7.3. Basic IFM, whitout rescue effect.
As far as the basic IFM is concerned, simulations are obtained by applying the function “sim1” of Example 7.1, where “C” is computed a part based on “pm”.

 

sim1<-function(o){

o<- o > 0 

if(any(o)){

cond<-ifelse(o,E,C)

o<-ifelse(runif(length(o))<cond,!o,o)

}

as.numeric(o)

}

 

# To carry out ns ≥ 2 simulations we repeat the procedure of Example  # 7.1 and we sum all together the ns matrices sim to obtain the matrix  # simT. Then, the [i,j]-th cell of simT says how many of the ns       # simulations assign 1 to the patch i. All the simulations start from # p1.

 

ns=10000

nt=10

 

simT<-matrix(0,nrow=length(p1),ncol=nt+1)

 

for(k in 1:ns) {

 sim<-matrix(0,nrow=length(p1),ncol=nt+1)

 sim[,1]<-p1

 for(t in 1:nt) sim[ ,t+1]<-sim1(sim[ ,t]) 

 simT=simT+sim

}

 

simT[1:13,]

      [,1]  [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]

 [1,] 10000  2374 7763 4049 6561 4943 5977 5271 5654  5453  5646

 [2,]     0  9996 3316 7811 4854 6745 5464 6380 5713  6253  5853

 [3,] 10000  2429 8197 3806 7093 4638 6489 5083 6175  5363  5898

 [4,]     0  9998 2860 7964 4416 6847 5206 6296 5550  6097  5675

 [5,]     0  9984 3044 7848 4518 6775 5262 6303 5597  6041  5813

 [6,]     0 10000 3414 7691 4840 6811 5497 6358 5852  6072  5911

 [7,]     0  7682 4392 5753 5212 5404 5335 5353 5371  5307  5458

 [8,] 10000  2879 7947 4244 6973 5004 6442 5385 6179  5631  5970

 [9,] 10000  3249 7817 4760 6756 5512 6289 5725 6087  5896  6016

[10,] 10000  3467 7741 4984 6774 5599 6368 5939 6138  5996  6105

[11,] 10000  3261 7780 4713 6845 5389 6385 5688 6125  5864  6075

[12,] 10000  2845 7783 4387 6684 5048 6271 5444 5961  5585  5945

[13,]     0  9990 3371 7765 4833 6810 5417 6375 5783  6268  5839

 

simT[1:13,]/ns

      [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]  [,11]

 [1,]    1 0.2374 0.7763 0.4049 0.6561 0.4943 0.5977 0.5271 0.5654 0.5453 0.5646

 [2,]    0 0.9996 0.3316 0.7811 0.4854 0.6745 0.5464 0.6380 0.5713 0.6253 0.5853

 [3,]    1 0.2429 0.8197 0.3806 0.7093 0.4638 0.6489 0.5083 0.6175 0.5363 0.5898

 [4,]    0 0.9998 0.2860 0.7964 0.4416 0.6847 0.5206 0.6296 0.5550 0.6097 0.5675

 [5,]    0 0.9984 0.3044 0.7848 0.4518 0.6775 0.5262 0.6303 0.5597 0.6041 0.5813

 [6,]    0 1.0000 0.3414 0.7691 0.4840 0.6811 0.5497 0.6358 0.5852 0.6072 0.5911

 [7,]    0 0.7682 0.4392 0.5753 0.5212 0.5404 0.5335 0.5353 0.5371 0.5307 0.5458

 [8,]    1 0.2879 0.7947 0.4244 0.6973 0.5004 0.6442 0.5385 0.6179 0.5631 0.5970

 [9,]    1 0.3249 0.7817 0.4760 0.6756 0.5512 0.6289 0.5725 0.6087 0.5896 0.6016

[10,]    1 0.3467 0.7741 0.4984 0.6774 0.5599 0.6368 0.5939 0.6138 0.5996 0.6105

[11,]    1 0.3261 0.7780 0.4713 0.6845 0.5389 0.6385 0.5688 0.6125 0.5864 0.6075

[12,]    1 0.2845 0.7783 0.4387 0.6684 0.5048 0.6271 0.5444 0.5961 0.5585 0.5945

[13,]    0 0.9990 0.3371 0.7765 0.4833 0.6810 0.5417 0.6375 0.5783 0.6268 0.5839