Agent-based modeling in R – habitat diversity and species richness

[This article was first published on Ecology in silico, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

How does habitat diversity affect species richness? Perhaps intuition
suggests that habitat diversity increases
species richness by facilitating niche or resource partitioning among
species. But, for a fixed
area, as habitat heterogeneity increases, the area that can be allocated to each
habitat type decreases. In a recent paper,
Allouche and colleagues (2012)
provide a theoretical and empirical treatment of the habitat
area-heterogeneity trade-off’s consequences for species richness. Both
treatments of the subject indicated that the relationship between
habitat heterogeneity and species richness may be unimodal, rather than
strictly increasing.

Conceptually, this is expected to occur when on the left side of the curve,
increasing habitat heterogeneity opens up new regions in niche space,
facilitating colonization by new species. However, as heterogeneity
continues to increase, each species has fewer habitat patches to utilize,
population sizes decrease, and local extinction risk increases
due to demographic stochasticity. To explore this idea theoretically,
Allouche et al. (2012) developed an individually based model using
a continuous time Markov process. The details of their modeling approach
can be found in
the supplementary material
to their article,
which I recommend. In this post, I’ll
demonstrate how to implement a discrete time version of their model in R.
Thanks to the agent-based modeling
working group at the University of Colorado for providing motivation to
code up model in R.

Model structure

This model is spatially implicit, with A equally connected sites.
Each site falls on an environmental condition axis, receiving
some value E that characterizes local conditions. The environmental
conditions for each site are uniformly distributed between two values that
dictate the range of environmental conditions in a focal area. The
local range of environmental conditions is
a subset of some global range. There are N species
in the regional pool that can colonize habitat patches. Each species has
some environmental optimum , and some niche width ,
which together define a Gaussian function for the probability of
establishment given a colonization attempt and a habitat patch
environmental condition E.

The image above illustrates the probability of establishing for five
species across the global range of environmental conditions possible.
For any focal area, the realized
range of environmental conditions is some subset of this global range.

It is assumed that all individuals that occupy a patch have the same
per-timestep probabilities
of death and reproduction. If an individual reproduces, the number of
offspring it produces is a Poisson distributed random variable, and each
individual offspring attempts to colonize one randomly selected site.
At each time-step, every site has an equal probability of a colonization
attempt by an individual from each species in the regional pool. Every
habitat patch holds only one individual.

Offspring and immigrants from
the regional pool do not displace individuals from habitat patches when
they attempt to colonize. In empty sites, offspring receive
colonization priority,
with regional colonization occurring after breeding. When multiple
offspring or immigrants from the regional pool could establish in
an empty site, one successful individual is randomly chosen to establish
regardless of species identity.

Parameters

The following parameters are supplied to the function
alloucheIBM():

A = number of sites;
N = number of species in the regional pool;
ERmin = global environmental conditions minimum;
ERmax = global environmental conditions maximum;
Emin = local environmental minimum;
Emax = local environmental maximum;
sig = niche width standard deviation for all species;
pM = per timestep probability of mortality;
pR = per timestep probability of reproduction;
R = per capita expected number of offspring; and
I = per timestep probability of attempted colonization by an immigrant
from the regional pool for each patch.

Implementation in R

The function alloucheIBM() does the majority of work for this model:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
<span class="line">alloucheIBM <span class="o"><-</span> <span class="kr">function</span><span class="p">(</span>A<span class="o">=</span><span class="m">100</span><span class="p">,</span> N<span class="o">=</span><span class="m">100</span><span class="p">,</span> ERmin<span class="o">=</span><span class="m">-30</span><span class="p">,</span> ERmax<span class="o">=</span><span class="m">30</span><span class="p">,</span> Emin<span class="o">=</span><span class="m">-30</span><span class="p">,</span> Emax<span class="o">=</span><span class="m">30</span><span class="p">,</span>
</span><span class="line">                    sig<span class="o">=</span><span class="m">5</span><span class="p">,</span> timesteps<span class="o">=</span><span class="m">1000</span><span class="p">,</span> pM<span class="o">=</span><span class="m">.1</span><span class="p">,</span> pR<span class="o">=</span><span class="m">1</span><span class="p">,</span> R<span class="o">=</span><span class="m">2</span><span class="p">,</span> I<span class="o">=</span><span class="m">0.1</span><span class="p">){</span>
</span><span class="line">  E <span class="o"><-</span> runif<span class="p">(</span>A<span class="p">,</span> Emin<span class="p">,</span> Emax<span class="p">)</span> <span class="c1"># habitat patch environment type</span>
</span><span class="line">  mu.i <span class="o"><-</span> runif<span class="p">(</span>N<span class="p">,</span> ERmin<span class="p">,</span> ERmax<span class="p">)</span> <span class="c1"># optimum environment</span>
</span><span class="line">  sigma <span class="o"><-</span> rep<span class="p">(</span>sig<span class="p">,</span> N<span class="p">)</span> <span class="c1"># niche width</span>
</span><span class="line">  Z <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>N<span class="p">))</span> <span class="c1"># normalization constant</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># Z ensures all species have equal Pr(establishment) in regional pool</span>
</span><span class="line">  <span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>N<span class="p">){</span>
</span><span class="line">    integrand <span class="o"><-</span> <span class="kr">function</span><span class="p">(</span>E<span class="p">)</span> <span class="p">{</span>
</span><span class="line">      exp<span class="p">(</span><span class="o">-</span><span class="p">((</span>E <span class="o">-</span> mu.i<span class="p">[</span>i<span class="p">])</span> <span class="o">^</span> <span class="m">2</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="m">2</span> <span class="o">*</span> sigma<span class="p">[</span>i<span class="p">]</span> <span class="o">^</span> <span class="m">2</span><span class="p">))</span>
</span><span class="line">    <span class="p">}</span>
</span><span class="line">    res <span class="o"><-</span> integrate<span class="p">(</span>integrand<span class="p">,</span> lower<span class="o">=</span>ERmin<span class="p">,</span> upper<span class="o">=</span>ERmax<span class="p">)</span>
</span><span class="line">    Z<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> <span class="m">1</span> <span class="o">/</span> res<span class="o">$</span>value
</span><span class="line">  <span class="p">}</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># probability of establishment|colonization attempt</span>
</span><span class="line">  Pcol <span class="o"><-</span> array<span class="p">(</span>dim<span class="o">=</span>c<span class="p">(</span>A<span class="p">,</span> N<span class="p">))</span>
</span><span class="line">  <span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>A<span class="p">){</span>
</span><span class="line">    <span class="kr">for</span> <span class="p">(</span>j <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>N<span class="p">){</span>
</span><span class="line">      Pcol<span class="p">[</span>i<span class="p">,</span> j<span class="p">]</span> <span class="o"><-</span> Z<span class="p">[</span>j<span class="p">]</span> <span class="o">*</span> exp<span class="p">(</span><span class="o">-</span><span class="p">((</span>E<span class="p">[</span>i<span class="p">]</span> <span class="o">-</span> mu.i<span class="p">[</span>j<span class="p">])</span> <span class="o">^</span> <span class="m">2</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="m">2</span><span class="o">*</span>sigma<span class="p">[</span>j<span class="p">]</span> <span class="o">^</span> <span class="m">2</span><span class="p">))</span>
</span><span class="line">    <span class="p">}</span>
</span><span class="line">  <span class="p">}</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># store niche data</span>
</span><span class="line">  species <span class="o"><-</span> rep<span class="p">(</span><span class="m">1</span><span class="o">:</span>N<span class="p">,</span> each<span class="o">=</span>A<span class="p">)</span>
</span><span class="line">  E <span class="o"><-</span> rep<span class="p">(</span>E<span class="p">,</span> N<span class="p">)</span>
</span><span class="line">  Pr.estab <span class="o"><-</span> c<span class="p">(</span>Pcol<span class="p">)</span>
</span><span class="line">  niche.d <span class="o"><-</span> data.frame<span class="p">(</span>species<span class="p">,</span> E<span class="p">,</span> Pr.estab<span class="p">)</span>
</span><span class="line">  niche.d <span class="o"><-</span> niche.d<span class="p">[</span>with<span class="p">(</span>niche.d<span class="p">,</span> order<span class="p">(</span>species<span class="p">,</span> E<span class="p">)),]</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># initialize output</span>
</span><span class="line">  state <span class="o"><-</span> array<span class="p">(</span><span class="m">0</span><span class="p">,</span> dim<span class="o">=</span>c<span class="p">(</span>timesteps<span class="p">,</span> A<span class="p">,</span> N<span class="p">))</span>
</span><span class="line">  richness <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> timesteps<span class="p">)</span>
</span><span class="line">  richness<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o"><-</span> <span class="m">0</span>
</span><span class="line">  p.occ <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> timesteps<span class="p">)</span>
</span><span class="line">  p.occ<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o"><-</span> <span class="m">0</span>
</span><span class="line">
</span><span class="line">  <span class="kr">for</span> <span class="p">(</span>t <span class="kr">in</span> <span class="m">2</span><span class="o">:</span>timesteps<span class="p">){</span>
</span><span class="line">    state<span class="p">[</span>t<span class="p">,,]</span> <span class="o"><-</span> state<span class="p">[</span>t<span class="m">-1</span><span class="p">,,]</span>
</span><span class="line">
</span><span class="line">    <span class="c1">## DEATHS ##</span>
</span><span class="line">    deaths <span class="o"><-</span> array<span class="p">(</span>rbinom<span class="p">(</span>A<span class="o">*</span>N<span class="p">,</span> <span class="m">1</span><span class="p">,</span> c<span class="p">(</span>state<span class="p">[</span>t<span class="p">,,])</span><span class="o">*</span>pM<span class="p">),</span> dim<span class="o">=</span>c<span class="p">(</span>A<span class="p">,</span> N<span class="p">))</span>
</span><span class="line">    state<span class="p">[</span>t<span class="p">,,]</span> <span class="o"><-</span> state<span class="p">[</span>t<span class="p">,,]</span> <span class="o">-</span> deaths
</span><span class="line">
</span><span class="line">    <span class="c1">## BIRTHS ##</span>
</span><span class="line">    pot.fecundity <span class="o"><-</span> array<span class="p">(</span>rpois<span class="p">(</span>A <span class="o">*</span> N<span class="p">,</span> lambda <span class="o">=</span> c<span class="p">(</span>state<span class="p">[</span>t<span class="p">,,]</span> <span class="o">*</span> R<span class="p">)),</span> dim<span class="o">=</span>c<span class="p">(</span>A<span class="p">,</span> N<span class="p">))</span>
</span><span class="line">    repro <span class="o"><-</span> array<span class="p">(</span>rbinom<span class="p">(</span>A<span class="o">*</span>N<span class="p">,</span> <span class="m">1</span><span class="p">,</span> pR<span class="p">),</span> dim<span class="o">=</span>c<span class="p">(</span>A<span class="p">,</span> N<span class="p">))</span>
</span><span class="line">    fecundity <span class="o"><-</span> repro <span class="o">*</span> pot.fecundity
</span><span class="line">    sum.fec <span class="o"><-</span> apply<span class="p">(</span>fecundity<span class="p">,</span> <span class="m">2</span><span class="p">,</span> sum<span class="p">)</span> <span class="c1"># number of offspring per species</span>
</span><span class="line">
</span><span class="line">    <span class="c1">## OFFSPRING COLONIZE ##</span>
</span><span class="line">    occupancy <span class="o"><-</span> apply<span class="p">(</span>state<span class="p">[</span>t<span class="p">,,],</span> <span class="m">1</span><span class="p">,</span> max<span class="p">)</span>
</span><span class="line">    <span class="kr">if</span> <span class="p">(</span>sum<span class="p">(</span>occupancy<span class="p">)</span> <span class="o"><</span> A <span class="o">&</span> sum<span class="p">(</span>sum.fec<span class="p">)</span> <span class="o">></span> <span class="m">0</span><span class="p">){</span> <span class="c1"># if empty sites & new offspring</span>
</span><span class="line">      empty.sites <span class="o"><-</span> which<span class="p">(</span>occupancy <span class="o">==</span> <span class="m">0</span><span class="p">)</span>
</span><span class="line">      occ.sites <span class="o"><-</span> which<span class="p">(</span>occupancy <span class="o">==</span> <span class="m">1</span><span class="p">)</span>
</span><span class="line">      <span class="c1"># randomly assign sites (empty & filled) to each offspring individual</span>
</span><span class="line">      col.sites <span class="o"><-</span> sample<span class="p">(</span><span class="m">1</span><span class="o">:</span>A<span class="p">,</span> sum<span class="p">(</span>sum.fec<span class="p">),</span> replace<span class="o">=</span><span class="k-Variable">T</span><span class="p">)</span>
</span><span class="line">      col.spec <span class="o"><-</span> rep<span class="p">(</span><span class="m">1</span><span class="o">:</span>N<span class="p">,</span> times<span class="o">=</span>sum.fec<span class="p">)</span> <span class="c1"># how many of each species colonizing</span>
</span><span class="line">      colonizing.offspring <span class="o"><-</span> array<span class="p">(</span><span class="m">0</span><span class="p">,</span> dim<span class="o">=</span>c<span class="p">(</span>A<span class="p">,</span> N<span class="p">))</span>
</span><span class="line">      <span class="kr">for</span><span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>length<span class="p">(</span>col.sites<span class="p">)){</span>
</span><span class="line">        colonizing.offspring<span class="p">[</span>col.sites<span class="p">[</span>i<span class="p">],</span> col.spec<span class="p">[</span>i<span class="p">]]</span> <span class="o"><-</span>
</span><span class="line">          colonizing.offspring<span class="p">[</span>col.sites<span class="p">[</span>i<span class="p">],</span> col.spec<span class="p">[</span>i<span class="p">]]</span> <span class="o">+</span> <span class="m">1</span>
</span><span class="line">      <span class="p">}</span>
</span><span class="line">      <span class="c1"># offspring attempting to colonize occupied sites fail to displace</span>
</span><span class="line">      colonizing.offspring<span class="p">[</span>occ.sites<span class="p">,]</span> <span class="o"><-</span> <span class="m">0</span>
</span><span class="line">
</span><span class="line">      <span class="c1"># which colonizing offspring can actually establish?</span>
</span><span class="line">      binom.mat <span class="o"><-</span> ifelse<span class="p">(</span>colonizing.offspring <span class="o">></span> <span class="m">0</span><span class="p">,</span> <span class="m">1</span><span class="p">,</span> <span class="m">0</span><span class="p">)</span>
</span><span class="line">      colonists <span class="o"><-</span> array<span class="p">(</span>rbinom<span class="p">(</span>n <span class="o">=</span> A <span class="o">*</span> N<span class="p">,</span>
</span><span class="line">                                size <span class="o">=</span> c<span class="p">(</span>colonizing.offspring<span class="p">),</span>
</span><span class="line">                                prob <span class="o">=</span> c<span class="p">(</span>binom.mat <span class="o">*</span> Pcol<span class="p">)),</span>
</span><span class="line">                         dim<span class="o">=</span>c<span class="p">(</span>A<span class="p">,</span> N<span class="p">))</span>
</span><span class="line">
</span><span class="line">      <span class="c1"># are there colonization conflicts (> 1 individual trying to colonize each site?)</span>
</span><span class="line">      attempting <span class="o"><-</span> apply<span class="p">(</span>colonists<span class="p">,</span> <span class="m">1</span><span class="p">,</span> sum<span class="p">)</span>
</span><span class="line">      <span class="kr">if</span> <span class="p">(</span>any<span class="p">(</span>attempting <span class="o">></span> <span class="m">1</span><span class="p">)){</span>
</span><span class="line">        <span class="c1"># resolve colonization conflicts</span>
</span><span class="line">        conflicts <span class="o"><-</span> which<span class="p">(</span>attempting <span class="o">></span> <span class="m">1</span><span class="p">)</span> <span class="c1"># which sites have conflicts</span>
</span><span class="line">        <span class="kr">for</span> <span class="p">(</span>k <span class="kr">in</span> conflicts<span class="p">){</span> <span class="c1"># for each conflict</span>
</span><span class="line">          <span class="c1"># how many indiv's of each spp attempting to simultaneously colonize?</span>
</span><span class="line">          n.attempting <span class="o"><-</span> rep<span class="p">(</span><span class="m">1</span><span class="o">:</span>N<span class="p">,</span> times <span class="o">=</span> colonists<span class="p">[</span>k<span class="p">,])</span>
</span><span class="line">          <span class="c1"># randomly select one successful from those attempting</span>
</span><span class="line">          successful <span class="o"><-</span> sample<span class="p">(</span>n.attempting<span class="p">,</span> size<span class="o">=</span><span class="m">1</span><span class="p">)</span>
</span><span class="line">          new.row <span class="o"><-</span> rep<span class="p">(</span><span class="m">0</span><span class="p">,</span> length.out<span class="o">=</span>N<span class="p">)</span>
</span><span class="line">          new.row<span class="p">[</span>successful<span class="p">]</span> <span class="o"><-</span> <span class="m">1</span>
</span><span class="line">          colonists<span class="p">[</span>k<span class="p">,]</span> <span class="o"><-</span> new.row <span class="c1"># individual becomes the only colonist</span>
</span><span class="line">        <span class="p">}</span>
</span><span class="line">      <span class="p">}</span>
</span><span class="line">      <span class="c1"># add successful colonists</span>
</span><span class="line">      state<span class="p">[</span>t<span class="p">,,]</span> <span class="o"><-</span> state<span class="p">[</span>t<span class="p">,,]</span> <span class="o">+</span> colonists
</span><span class="line">    <span class="p">}</span>
</span><span class="line">
</span><span class="line">    <span class="c1">## IMMIGRANTS COLONIZE ##</span>
</span><span class="line">    occupancy <span class="o"><-</span> apply<span class="p">(</span>state<span class="p">[</span>t<span class="p">,,],</span> <span class="m">1</span><span class="p">,</span> max<span class="p">)</span>
</span><span class="line">    <span class="kr">if</span><span class="p">(</span>sum<span class="p">(</span>occupancy<span class="p">)</span> <span class="o"><</span> A<span class="p">){</span>
</span><span class="line">      empty.sites <span class="o"><-</span> which<span class="p">(</span>occupancy <span class="o">==</span> <span class="m">0</span><span class="p">)</span>
</span><span class="line">      <span class="c1"># which species immigrate to each site?</span>
</span><span class="line">      immigration <span class="o"><-</span> array<span class="p">(</span>rbinom<span class="p">(</span>length<span class="p">(</span>empty.sites<span class="p">)</span><span class="o">*</span>N<span class="p">,</span>
</span><span class="line">                                  <span class="m">1</span><span class="p">,</span> I<span class="p">),</span> dim<span class="o">=</span>c<span class="p">(</span>length<span class="p">(</span>empty.sites<span class="p">),</span> N<span class="p">))</span>
</span><span class="line">      <span class="c1"># which immigrants establish?</span>
</span><span class="line">      Pest <span class="o"><-</span> immigration <span class="o">*</span> Pcol<span class="p">[</span>empty.sites<span class="p">,</span> <span class="p">]</span>
</span><span class="line">      establishment <span class="o"><-</span> array<span class="p">(</span>rbinom<span class="p">(</span>length<span class="p">(</span>Pest<span class="p">),</span> <span class="m">1</span><span class="p">,</span> c<span class="p">(</span>Pest<span class="p">)),</span>
</span><span class="line">                             dim<span class="o">=</span>c<span class="p">(</span>length<span class="p">(</span>empty.sites<span class="p">),</span> N<span class="p">))</span>
</span><span class="line">
</span><span class="line">      <span class="c1"># resolve conflicts arising from simultaneous colonization</span>
</span><span class="line">      col.attempts <span class="o"><-</span> apply<span class="p">(</span>establishment<span class="p">,</span> <span class="m">1</span><span class="p">,</span> sum<span class="p">)</span>
</span><span class="line">      <span class="kr">if</span> <span class="p">(</span>any<span class="p">(</span>col.attempts <span class="o">></span> <span class="m">1</span><span class="p">)){</span> <span class="c1"># if individuals are trying to simultaneously colonize</span>
</span><span class="line">        conflicts <span class="o"><-</span> which<span class="p">(</span>col.attempts <span class="o">></span> <span class="m">1</span><span class="p">)</span> <span class="c1"># which empty sites have conflicts</span>
</span><span class="line">        <span class="kr">for</span> <span class="p">(</span>k <span class="kr">in</span> conflicts<span class="p">){</span> <span class="c1"># for each conflict</span>
</span><span class="line">          <span class="c1"># how many of individuals of each species are attempting to simultaneously colonize?</span>
</span><span class="line">          attempting <span class="o"><-</span> rep<span class="p">(</span><span class="m">1</span><span class="o">:</span>N<span class="p">,</span> times <span class="o">=</span> establishment<span class="p">[</span>k<span class="p">,])</span>
</span><span class="line">          <span class="c1"># successful individual randomly selected from those attempting</span>
</span><span class="line">          successful <span class="o"><-</span> sample<span class="p">(</span>attempting<span class="p">,</span> size<span class="o">=</span><span class="m">1</span><span class="p">)</span>
</span><span class="line">          new.row <span class="o"><-</span> rep<span class="p">(</span><span class="m">0</span><span class="p">,</span> length.out<span class="o">=</span>N<span class="p">)</span>
</span><span class="line">          new.row<span class="p">[</span>successful<span class="p">]</span> <span class="o"><-</span> <span class="m">1</span>
</span><span class="line">          establishment<span class="p">[</span>k<span class="p">,]</span> <span class="o"><-</span> new.row
</span><span class="line">        <span class="p">}</span>
</span><span class="line">      <span class="p">}</span>
</span><span class="line">      <span class="c1"># add successful immigrants</span>
</span><span class="line">      state<span class="p">[</span>t<span class="p">,</span> empty.sites<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> state<span class="p">[</span>t<span class="p">,</span> empty.sites<span class="p">,]</span> <span class="o">+</span> establishment
</span><span class="line">    <span class="p">}</span>
</span><span class="line">    richness<span class="p">[</span>t<span class="p">]</span> <span class="o"><-</span> sum<span class="p">(</span>apply<span class="p">(</span>state<span class="p">[</span>t<span class="p">,,],</span> <span class="m">2</span><span class="p">,</span> max<span class="p">))</span>
</span><span class="line">    p.occ<span class="p">[</span>t<span class="p">]</span> <span class="o"><-</span> length<span class="p">(</span>which<span class="p">(</span>state<span class="p">[</span>t<span class="p">,,]</span> <span class="o">==</span> <span class="m">1</span><span class="p">))</span> <span class="o">/</span> A
</span><span class="line">  <span class="p">}</span>
</span><span class="line">  <span class="kr">return</span><span class="p">(</span>list<span class="p">(</span>richness<span class="o">=</span>richness<span class="p">,</span> p.occ <span class="o">=</span> p.occ<span class="p">,</span> state<span class="o">=</span>state<span class="p">,</span> niche.d<span class="o">=</span>niche.d<span class="p">))</span>
</span><span class="line"><span class="p">}</span>
</span>

The function returns a list containing a vector of species richness at
each timestep, the proportion of sites occupied at each timestep,
a state array containing all occupancy information for each patch,
species, and timestep, and lastly
a dataframe containing information on the niches of each species in
the regional pool.

Using this function we can begin to explore the dynamics of the model
through time:

1
2
<span class="line">out <span class="o"><-</span> hostIBM<span class="p">(</span>pM<span class="o">=</span><span class="m">.1</span><span class="p">,</span> Emin<span class="o">=</span><span class="m">-30</span><span class="p">,</span> Emax<span class="o">=</span><span class="m">30</span><span class="p">,</span> sig<span class="o">=</span><span class="m">8</span><span class="p">,</span> R<span class="o">=</span><span class="m">8</span><span class="p">,</span> N<span class="o">=</span><span class="m">100</span><span class="p">,</span> A<span class="o">=</span><span class="m">300</span><span class="p">,</span> I<span class="o">=</span><span class="m">.01</span><span class="p">,</span> timesteps<span class="o">=</span><span class="m">1000</span><span class="p">)</span>
</span><span class="line">plot<span class="p">(</span>out<span class="o">$</span>richness<span class="p">,</span> type<span class="o">=</span><span class="s">"l"</span><span class="p">,</span> xlab<span class="o">=</span><span class="s">"Timestep"</span><span class="p">,</span> ylab<span class="o">=</span><span class="s">"Species richness"</span><span class="p">)</span>
</span>

Repeating this process a few times, we can get a picture of the
expected species richness and it’s variance for some set of parameters.

Finally, we can address the issue of habitat heterogeneity and its
effect on species richness. There are many ways to approach this issue,
and many parameter combinations to consider. Allouche et al. (2012)
provides a thorough treatment of the subject; I’ll
demonstrate just one result: that under certain conditions, species
richness
peaks at intermediate levels of habitat heterogeneity.

To construct a range of habitat heterogeneity values, let’s construct
an interval and take subsequently narrower intervals centered around the
middle of the original interval.

1
2
3
4
5
6
7
<span class="line">ERmin <span class="o"><-</span> <span class="m">-50</span>
</span><span class="line">ERmax <span class="o"><-</span> <span class="m">50</span>
</span><span class="line">global.median <span class="o"><-</span> median<span class="p">(</span>c<span class="p">(</span>ERmin<span class="p">,</span> ERmax<span class="p">))</span>
</span><span class="line">n.intervals <span class="o"><-</span> <span class="m">30</span>
</span><span class="line">lower.limits <span class="o"><-</span> seq<span class="p">(</span>ERmin<span class="p">,</span> global.median<span class="m">-.5</span><span class="p">,</span> length.out<span class="o">=</span>n.intervals<span class="p">)</span>
</span><span class="line">upper.limits <span class="o"><-</span> seq<span class="p">(</span>ERmax<span class="p">,</span> global.median<span class="p">,</span> length.out<span class="o">=</span>n.intervals<span class="p">)</span>
</span><span class="line">hab.het <span class="o"><-</span> upper.limits <span class="o">-</span> lower.limits <span class="c1"># habitat heterogeneity</span>
</span>

Here are the intervals:

Now, for each interval, we can iteratively run the model and track
species richness. Because species richness tends to vary through time,
let’s take the mean of the final 100 timesteps as a measure of species
richness for each model run, and record the standard deviation to track
variability.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
<span class="line">n.iter <span class="o"><-</span> <span class="m">3</span> <span class="c1"># number of iterations per interval</span>
</span><span class="line">sd.rich <span class="o"><-</span> array<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> dim<span class="o">=</span>c<span class="p">(</span>n.intervals<span class="p">,</span> n.iter<span class="p">))</span>
</span><span class="line">end.rich <span class="o"><-</span> array<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> dim<span class="o">=</span>c<span class="p">(</span>n.intervals<span class="p">,</span> n.iter<span class="p">))</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>n.intervals<span class="p">){</span>
</span><span class="line">  <span class="kr">for</span> <span class="p">(</span>iter <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>n.iter<span class="p">){</span>
</span><span class="line">    Emin <span class="o"><-</span> lower.limits<span class="p">[</span>i<span class="p">]</span>
</span><span class="line">    Emax <span class="o"><-</span> upper.limits<span class="p">[</span>i<span class="p">]</span>
</span><span class="line">    out <span class="o"><-</span> hostIBM<span class="p">(</span>pM<span class="o">=</span><span class="m">.9</span><span class="p">,</span> Emin<span class="o">=</span>Emin<span class="p">,</span> Emax<span class="o">=</span>Emax<span class="p">,</span> sig<span class="o">=</span><span class="m">10</span><span class="p">,</span> pR<span class="o">=</span><span class="m">1</span><span class="p">,</span> R<span class="o">=</span><span class="m">1</span><span class="p">,</span>
</span><span class="line">                   N<span class="o">=</span><span class="m">100</span><span class="p">,</span> A<span class="o">=</span><span class="m">100</span><span class="p">,</span> I<span class="o">=</span><span class="m">.9</span><span class="p">,</span> timesteps<span class="o">=</span><span class="m">300</span><span class="p">)</span>
</span><span class="line">    timesteps <span class="o"><-</span> length<span class="p">(</span>out<span class="o">$</span>richness<span class="p">)</span>
</span><span class="line">    end.rich<span class="p">[</span>i<span class="p">,</span> iter<span class="p">]</span> <span class="o"><-</span> mean<span class="p">(</span>out<span class="o">$</span>richness<span class="p">[</span>timesteps<span class="m">-100</span><span class="o">:</span>timesteps<span class="p">])</span>
</span><span class="line">    sd.rich<span class="p">[</span>i<span class="p">,</span> iter<span class="p">]</span> <span class="o"><-</span> sd<span class="p">(</span>out<span class="o">$</span>richness<span class="p">[</span>timesteps<span class="m">-100</span><span class="o">:</span>timesteps<span class="p">])</span>
</span><span class="line">  <span class="p">}</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">end.richness <span class="o"><-</span> c<span class="p">(</span>t<span class="p">(</span>end.rich<span class="p">))</span>
</span><span class="line">end.sd <span class="o"><-</span> c<span class="p">(</span>t<span class="p">(</span>sd.rich<span class="p">))</span>
</span><span class="line">interval <span class="o"><-</span> rep<span class="p">(</span>hab.het<span class="p">,</span> each<span class="o">=</span>n.iter<span class="p">)</span>
</span><span class="line">plot.d <span class="o"><-</span> data.frame<span class="p">(</span>interval<span class="p">,</span> end.richness<span class="p">,</span> end.sd<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># visualize the results</span>
</span><span class="line">require<span class="p">(</span>ggplot2<span class="p">)</span>
</span><span class="line">range <span class="o"><-</span> aes<span class="p">(</span>ymax<span class="o">=</span>end.richness <span class="o">+</span> end.sd<span class="p">,</span> ymin<span class="o">=</span>end.richness <span class="o">-</span> end.sd<span class="p">)</span>
</span><span class="line">ggplot<span class="p">(</span>plot.d<span class="p">,</span> aes<span class="p">(</span>x<span class="o">=</span>interval<span class="p">,</span> y<span class="o">=</span>end.richness<span class="p">))</span> <span class="o">+</span>
</span><span class="line">  geom_point<span class="p">(</span>shape<span class="o">=</span><span class="m">1</span><span class="p">,</span> size<span class="o">=</span><span class="m">2</span><span class="p">)</span> <span class="o">+</span> geom_errorbar<span class="p">(</span>range<span class="p">)</span> <span class="o">+</span>
</span><span class="line">  theme_classic<span class="p">()</span> <span class="o">+</span>
</span><span class="line">  xlab<span class="p">(</span><span class="s">"Habitat heterogeneity"</span><span class="p">)</span> <span class="o">+</span> ylab<span class="p">(</span><span class="s">"Species richness"</span><span class="p">)</span>
</span>

Of course, the shape of this relationship is sensitive to the parameters.
As an example, changing niche width to increase or decrease niche overlap
will mediate the strength of interspecific competition for space.
Also, increasing reproductive rates may buffer each species from
stochastic extinction so that the relationship between environmental
heterogeneity and richness is monotonically increasing. Furthermore,
here I centered all intervals around the same value, but the exact
position of the environmental heterogeneity interval will affect the net
establishment probability for each site, depending on how the interval
relates to species niches. The parameter space is yours to explore.

These types of stochastic simulation models are fairly
straightforward to implement in R. Indeed there’s a package
dedicated to facilitating the implementation of such models:
simecol. There’s even a book:
A Practical Guide to Ecological Modelling: Using R as a Simulation Platform.

Full reference

Allouche O, Kalyuzhny M, Moreno-Rueda G, Pizarro M, & Kadmon R. (2012)
Area-heterogeneity tradeoff and the diversity of ecological communities.
Proceedings of the National Academy of Sciences of the United States of
America, 109 (43): 17495-17500.

To leave a comment for the author, please follow the link and comment on their blog: Ecology in silico.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)