A hierarchical Bayesian state space model for fluctuating rodent populations, with application to the populations at Finse, Norway

fluctuating rodent populations, with application to the populations at Finse, Norway

by

Magnus Nygård Osnes

Thesis

for the degree of Master of Science Main field

in Ecology and Evolution

60 credits

Department of Biosciences

Faculty of Mathematics and Natural Sciences UNIVERSITY OF OSLO

2018

© Magnus Nygård Osnes 2018

A hierarchical Bayesian state space model for fluctuating rodent populations, with application to the populations at Finse, Norway

http://www.duo.uio.no/

Printed: Reprosentralen, University of Oslo

I want to thank my principal supervisor, Torbjørn Håkan Ergon. To even begin and digest the comprehensive amount of work that has been done on fluctuating rodent populations, is overwhelming... Thank you so much for all the encouragement, discussions and time you took with me. It helped me a whole lot with keeping structure in a confused mind. In helping me define this project, you have allowed me to combine the fields in which I love to work, and for that I am thankful. I have been learning new exciting things till the last moment in the production of this thesis.

Erik Framstad, thank you for the warm welcome in the fieldwork at Finse, and also for your insightful suggestions. If you are ever in need of an field assistant, I am standing by.

Geir Storvik, thank you kindly for your suggestions on the statistical methods. The model development would not have been possible without your help.

For the last six years at Blindern I have been torn between the pure biology studies and the dirty statistical ways. I want to thank all my friends in the SDATA-crew. I already miss the time at the 8th floor of Abel. I also want to thank the BIO-crew. Ingvild, for helping me keeping motivation up. Lisa for being my partner in crime throughout my six years at Blindern. Silje, Malin, Harald, and Julie for making lunchtime a great priority.

I want to thank my brother for proofreading, in general being very supportive. You are of one of the smartest people I know, and always very dependable. I also want to thank the rest of my family for allowing me to go on about rodents in every family gathering.

A special thanks goes to Katinka Lund Bergskaug who kept me sane in times of frus- tration and confusion, and was enormously patient with me throughout the entire process.

Thanks you for your boundless kindness.

I am grateful to Martyn Plummer, for providing the modified version of JAGS that was used in this thesis. I also want to acknowledge the help received from the Department for Research Computing at USIT, the University of Oslo IT-department. The computa- tions and simulations were performed on resources provided by UNINETT Sigma2 - the National Infrastructure for High Performance Computing and Data Storage in Norway.

January 2018, Oslo Magnus Nygård Osnes

i

Population abundance fluctuations of rodents have been observed to dampen in a syn- chronized fashion in multiple locations across Europe. Hypotheses for the causes of the dampening exist, but have not been confirmed. In need to determine the causes, precise descriptions of the population cycles in each location is essential. To do this and produce results that are comparable between studies with differences in sampling procedures, one must incorporate known sources of variation in the sampling procedures.

In this work, I analyze one of the longest running time series on rodents in Scandinavia.

The data contains snap-trap captures from two sites on seven small mammals species (5 rodent species and 2 species of shrews). Previous analyses have used the number of cap- tures per effort (numbers of captures divided by number of traps multiplied by the number of hundred days) when making inferences about the population abundance dynamics. This relies on the assumption that this index is proportional to the rodent abundance and ignores important problems caused by the sampling procedure. In this thesis, I develop a new ob- servation model, that incorporates sources of variance that are known to follow with the removal sampling method, such as trap-saturation, population depletion and competition for traps in a multi-species setting. Moreover, the model produces estimates of abundance itself, not an index of it. Through simulations, I show that this model produces reliable estimates of abundance under combinations of parameter values that are likely to be ob- served in the field. The observation model is incorporated into a hierarchical Bayesian state space model, which is used to estimate population abundance dynamics of three of the rodent species: Norway lemmings (Lemmus lemmus), root vole (Microtus oeconomus) and field vole (M. agrestis). In light of the Europe wide dampening of population cycles, which includes these rodents, I test a hypothesis of shifted dynamics between the intervals [1970,1994]and[1994,2017]. I do this by fitting the state space model with independent sets of coefficients in each half of the time series. The resulting posterior distributions indicate changes in direct and delayed density dependence and underlying dynamics (pe- riodicity and amplitudes), for root vole and field vole. Until 1994, the estimated dynamics are cycles with periods of 3 5 years. After this, the estimated dynamics are dampened cycles with longer periods, and the credibility intervals of the posterior distributions cover

iii

section of the process model to seasonal components indicates that for root vole there has been a decrease in a positive effect of spring density on population growth rate during summer, at one of the sites. For the other species, sites and time intervals, the seasonal dissections are inconclusive.

iv

Notation xi

1 Introduction 1

2 Methods 8

2.1 Locations . . . 8

2.2 Data gathering procedure . . . 9

2.3 Data . . . 11

2.4 Formatting of the data . . . 11

2.5 Problems with the data and handling . . . 12

2.6 Study species and environment . . . 13

2.7 Modeling abundance . . . 14

2.8 State space model . . . 18

2.9 Fitting the models (Bayesian paradigm) . . . 18

2.10 The observation model . . . 21

2.10.1 Model specification . . . 22

2.10.1.1 The observations . . . 22

2.10.1.2 Inter-arrival times . . . 23

2.10.1.3 Order of the captured species . . . 25

2.10.1.4 The prior distribution ofb_s . . . 27

2.10.2 The effect of sampling on subsequent trap-nights . . . 27

2.11 Process model . . . 29

2.11.1 Fitting the AR(4) model . . . 30 vi

2.11.2.2 Density dependent coefficients . . . 32

2.11.2.3 Process errors . . . 33

2.12 Investigating differences in dynamics . . . 33

2.13 Interpretation of the density dependent coefficients from the process model 36 2.14 State space model evaluation . . . 38

2.14.1 Convergence diagnostics of the MCMC chains . . . 38

2.14.2 Performance of the observation model . . . 39

2.14.3 Residual process component . . . 40

3 Results 42 3.1 Performance on simulated data . . . 42

3.2 Relation between SSM estimates and the CE100 index . . . 46

3.3 Goodness of fit of the SSM . . . 48

3.4 The hazard coefficients . . . 50

3.5 Estimates of abundance . . . 51

3.6 Annual density dependence and underlying dynamics . . . 53

3.6.1 Estimated dynamical properties . . . 60

3.7 Dissection to seasonal changes . . . 64

3.7.1 The remaining seasonal posterior distributions . . . 71

4 Discussion 76 4.1 Observation model properties . . . 77

4.2 Weaknesses of the model . . . 77

4.3 Estimated abundances . . . 79

4.3.1 Reappearance of the cycles of Norway lemming and root vole? . . 79

4.3.2 Shapes of the cycles . . . 79

4.4 Annual dynamics . . . 80

4.5 Seasonal dynamics and dissections . . . 81

4.6 Changes in dynamics and relation to the hypotheses . . . 82 vii

4.6.3 The overgrazing hypothesis and nutrient cycling . . . 84

4.6.4 Climate and snow condition hypothesis . . . 85

4.6.5 Multi factor hypothesis . . . 85

5 Summary and conclusions 87 6 Further studies 89 Bibliography 91 A Theory and additional detail 101 A.1 The bernoulli trick . . . 101

A.2 Hazard function plotsb . . . 101

A.3 Derivation of the distribution of the latent variable t . . . 102

A.4 Removal model . . . 103

A.5 RRMSE . . . 105

A.6 Coverage rate . . . 105

A.7 The annual-level and direct and delayed density dependent effects . . . . 105

A.7.1 Annual fall to fall effects (AR4) . . . 105

A.7.2 Annual spring to spring effects, AR(4) . . . 107

A.8 The distribution of theF_kandS_k . . . 107

B Model diagnostics and summaries 109 B.1 Abundance estimates . . . 109

B.2 Full coefficient tables for the seasonal density dependencies . . . 113

B.2.1 Residual process components of the AR4 model . . . 116

B.3 Correlations of the seasonal density dependent parameters . . . 117

B.4 Correlations of the annual density dependent parameters . . . 118

B.5 Tables of decomposition and posterior distributions . . . 119

B.6 Trace and density plots of the MCMC . . . 125 viii

C Environmental variables 134 C.1 Data from SeNorge . . . 134 C.2 NAO-index . . . 136

D Supplementary 137

D.1 Observation model in JAGS . . . 144 D.2 State space model in JAGS with interval specific coefficients . . . 145

ix

The notation is used as follows. Symbols that use subscripts are used in the order: i, k, s, r, and p as a superscript. If one of these indexes are left out, I am considering a simplified scenario, where that index is not of interest. If the index is replaced by a dot, then the symbol is summed over that index. For vectors I use the same letter as the original symbol, but written inboldface.

List of indices and constants

j Index running over events throughout a trap-night. The total number of events within a trap-night is the sum of all species that are caught.

J The total number of trap events in a trap-night. This is always equal to the total number of individual rodents caught, that is,J=Â^S_s=1Ys. (See list of symbols) i index running over trap-nights (i as short for interval).

I The total number of trap-nights in a survey.

p Index running over the halves of the seasons series. p = S means spring and p = F means fall.

k Index running over the years 1970-2017

K The total number of years. This index is only used in the simulation study.

s Index running over the number of species in the study.s=1 means species 1. The total number of species denoted byS.

xi

r index running over the halves of the time series. r=1 means 1970-1994 while r=2 means 1995-2017.

F The total number of traps in a survey.

List of symbols

g N_0,k|N_{0,k 1},yyy System process. This describes how the system changes through time given its parametersyyy (in this study density dependence and noise)

f Y_k|N_0,k,qqq The observation model. This describes how the observations (the number of captured rodents) have been generated from the state (the unknown true abun- dance) by the sampling procedure.

y y

y The set of parameters that is used in the system process qqq The set of parameters that is used in the observation model aj,s The capture hazard rate of speciessat event j

aaa_s The vector(a₁,a₂, . . . ,a_Ys), containing the piecewise constant hazard rates of species jthroughout one trap-night.

bs Constant of proportionality that is used in the hazard function. This represents the efficiency of the traps and behavior of the species.

hj The time the j’th capture was made

t_j The j’th interarrival time. That is, the time difference in time between subsequent captures:h_j h_{j 1}.

t Is the vector:(t₁,t₂, . . . ,t_Y.), that contains the interarrival timest_j

oj a categorical variable the denotes the species that is caught at capture event jwithin a trap-night.

xii

Y_i,k,s^p The observed number of captures in trap-night i, in year k, for species s, and in season p. In cases where some of these subscripts have been left out, I am referring to a simplified situation where those indexes are not of interest.

Y Vector of the observed number of captures on consecutive trap-nigths (i) N_k^S The abundance in spring in yeark

Nc_j,s Counting variable for number of captures of species s at event j in a trap-night.

This variable contains the number caught of each species after the first jcaptures.

N_k^F The abundance in fall in yeark S_k The logarithm ofN_k^S

F_k The logarithm ofN_k^F

R_Wi,k The population growth rate over winter in year k.

R_Su,k The population growth rate over summer in year k.

e_Wi,k System process error over winter. That is, the residual process component from the fitted AR(4) model for the winter equation

eSu,k System process error over summer. hat is, the residual process component from the fitted AR(4) model for the summer equation

µ_k^S The expected mean of of the spring log abundance, predicted from the AR(4) model.

xiii

model.

sWi,s The system process error standard deviation for the winter equation s_Su,s The system process error standard deviation for the summer equation

a_k,r Density dependent coefficient withk years of lag acting on the growth rate over winter. The r index indicates if it belongs to the interval 1970-1994 (r=1) or 1995-2017 (r=2).

b_k,r Density dependent coefficient withk years of lag acting on the growth rate over summer. Ther index indicates if it belongs to the interval 1970-1994 (r=1) or 1995-2017 (r=2).

W1 Annual direct density dependence. With the AR(4) structure, this made up by a combination of the seasonal coefficients in the following way:

W1= (a1+b1+b2+a2+a1b1),

W2 Annual delayed density dependence. With the AR(4) structure, this made up by a combination of the seasonal coefficients in the following way:

W₂= (a₃+b₃+b₄+a₄+a₁b₃+a₃b₁ a₂b₂), N_0,H High abundance species in the simulation study.

N_0,M Medium abundance species in the simulation study. Always set to half of the high abundance speciesN_0,H

N_0,L Low abundance species in the simulation study. Always set to one tenth of the high abundance speciesN_0,H

List of distributions

⇠dexp(·) Drawn from the exponential distribution

⇠dcat(·) Drawn from the catergorical distribution xiv

⇠dnorm(·) Drawn from the normal distribution List of abbreviations

K1 The productive grid with the longer productive season of the two grids. Almost all captures of root vole and field vole were obtained in this grid

T1 The less productive grid located. In general there was lower numbers of other species than Norway lemming at this grid

Trap-night The trapping interval from the traps are set, until they are checked the next day.

CE100 Captures per 100 Trap-nights of effort. This index has been used as an abundance measure in a number of studies. In this work, I estimate the distribution of the latent unobserved abundance by state-space modeling.

JAGS Just Another Gibbs Sampler. A software designed for fitting MCMC’s through the Gibbs sampling algorithm

MCMC Markov Chain Monte Carlo

Seasonal equations Referring to the fitted equtions of the AR(4) model, where the winter equation is:

S_k=a₀+ (1+a₁)F_{k 1}+a₂S_{k 1}+a₃F_{k 2}+a₄S_{k 2}+eWi,k, and the summer equations is:

F_k=b₀+ (1+b₁)S_k+b₂F_{k 1}+b₃S_{k 1}+b₄F_{k 2}+eSu,k

AR(2) Second order log-linear autoreggressive model.

AR(4) Fourth order log-linear autoreggressive model.

burn-in Period of iterations in MCMC chains that are affected by initial values, and have not yet converged. This must be discarded to have a representible sample of the true distribution

xv

MA Short forMicrotus agrestis (field vole). This abbreviation is only used in some figures and tables.

MO Short forMicrotus oeconomus (root vole).This abbreviation is only used in some figures and tables.

NAO The North Atlantic Oscillation

SSM State space model. A decomposition of the dynamics into a system processg N_0,k|N_{0,k 1},yyy and observation model, f Y_k|N_0,k,qqq .

xvi

Chapter 1 Introduction

Since Elton described population cycles to the scientific community (Elton, 1924, 1942), they have gained widespread attention and have been studied extensively. The fluctuations of rodent populations in alpine environments were what originally sparked Elton’s interest, and much effort have been devoted to determine the causes of their cycles. However, many aspects of the cycles are still poorly understood and there are many questions left to answer. The problem is the high complexity of ecosystems that manifest through a network of interacting factors.

Accurate descriptions of this network through data is difficult, and generation of such data can be very time-consuming and difficult to produce. Acknowledging these challenges motivates the development of models that make best use of data. This thesis is motivated by such a problem, and tries to make better use of a very long time series (still running) in an analysis of a rodent populations in Scandinavia.

When considering population cycles in ecology, one usually mean cycles that have a reg- ular period of three or more years (Berryman, 2002). The population fluctuations of rodents in alpine ecosystems are examples of such cycles, and have been described for many species (Stenseth and Ims, 1993b). They usually exhibit periods of 3-5 years between successive lows, with highly variable peak densities (Hansson and Henttonen, 1985). A well known example is the 3-5 year cyclicity of Norway Lemmings, (Lemmus lemmus) where the population den- sities change from very low abundance, to high densities over large areas. At high densities they are often seen migrating in large numbers (see Henttonen and Kaikusalo (1993)). In peak years, Norway lemmings and other voles dominate in biomass in the alpine ecosystem. As a consequence, many species have become well adapted to their fluctuations and the presence of some predators are dependent on them (for example, snowy owls (Bubo scandicus) and arctic

foxes (Vulpes lagopus)) (Gilg et al., 2006). Fluctuations in the population size of the rodents are known to propagate to numerical responses in predators and other species (Hörnfeldt, 1978;

Ims and Fuglei, 2005; Ims et al., 2013). Another characteristic of rodent cycles is that sympatric cyclic vole populations are usually synchronous. By this I mean that their populations tend to peak and decline through the same periods of time. This suggests that their cycles are some- how linked by their causes, or that they experience the same perturbations which force them into phase (Moran, 1953). Although cyclicity seems to characterize many small rodent popula- tions, there are populations that only fluctuate annually, and thus are not cyclic in the ecological sense. In fact, a gradient between cyclic populations in northern Scandinavia and non-cyclic populations at lower latitudes, has been observed (Hansson and Henttonen, 1985; Hanski et al., 1991; Bjørnstad et al., 1995). This suggests that important clues to the causes lie in differences between these ecosystems.

The data used in this thesis was collected at Hardangervidda at Finse, close to the Finse Alpine Research Center. At this location, rather stable lemming cycles occurred regularly until the mid 90’s. After this, there was a period with dampened cycles, where the large peaks were more or less absent until 2014 (described until 2007 in Kausrud et al. (2008)). Notably, this period conforms to an observed pattern of dampening amplitudes in rodent cycles across Eu- rope (Hörnfeldt, 2004; Hörnfeldt et al., 2005; Bierman et al., 2006; Ims et al., 2008; Cornulier et al., 2013; Gouveia et al., 2015). Small mammal herbivores form an essential link between pri- mary production and higher trophic levels and are keystone herbivores (Ims and Fuglei, 2005).

Phenomenological shifts in the cycles of their populations is therefore expected to cause large changes in the species composition at these locations. Efforts should therefore be dedicated to understand the causes of this pattern. In 2014 there was a reappearance of a peak year for Norway lemming at Finse, breaking the dampened pattern. This was followed by another one in 2017, suggesting that the regular high amplitude peaks might have come back (data to be presented in this thesis).

What needs to be investigated is the cause of the dampened cycles. Dampening can occur if the strength of the factors governing the population dynamics have changed, but it does not have to be so. It is well known that population cycles can seem to fade out, even when the factors governing them have not changed. This can happen if the growth rate of a population is dominated by noise terms for some time (e.g. Ims et al. (2008); Royama (1981)). It can also happen if different initial population sizes have stable limit cycles of different periods, and disturbances moves the system to the proximity of different points of attraction (Turchin, 2003;

Frisman et al., 2016). These complications can result in temporal subsets of the time series

looking more similar than the time series as a whole, and the resulting dynamics appear differ- ent even when the causal factors are unaltered. However, the synchronized dampening pattern across Europe implies that the changes are less likely to have arisen as a result of coincidental simultaneous regime shifts. Care must be taken not to misinterpret this, as large scale synchrony in cycles can be caused by correlated environmental variabilities, even if the cause of the cycles themselves is a completely different mechanism. That is, the cause of the synchrony may be entirely different than the cause(s) of the cycles (Moran, 1953; Royama, 1992), what is referred to a Moran’s theorem in population ecology.

For my thesis I have available a time series table containing data of small mammal captures from 1970 till 2017 at Finse. This contains captures of five rodents species: Norway lemming (Lemmus lemmus), root vole (Microtus oeconomus), field vole (Microtus agrestis), bank vole (Myodes glareolus, formerly Clethrionomys glareolus), grey-sided vole (Myodes rufocanus, formerly Clethrionomys rufocanus), and two species of shrews: Common shrew (Sorex aran- eus) and Eurasian pygme shrew (Sorex minutus). In this work I analyze the dynamics of the three most common species in the data and describe their dynamics during the cyclic phase 1970-1994, and compare it to the dynamics during the period 1994-2017.

An approach to modeling that has been successful in describing the dynamics of many ro- dent populations is the description of factors affecting the population growth rates in terms of direct density dependence and delayed density dependence (Bjørnstad et al., 1995; Bierman et al., 2006; Cornulier et al., 2013). Direct density dependence affects the population growth rate without delay. Such factors can be intrinsic, e.g. competition among individuals, maternal effects, genetics, food restrictions, or extrinsic factors such as generalist predators, parasites, epidemic diseases and environmental effects. It is generally recognized that direct density de- pendence has a stabilizing effect on population fluctuations (Hanski et al., 1991), but in theory it can also produce population cycles (Turchin, 2003; Barraquand et al., 2014). For northern eu- ropean rodents, the importance of delayed density dependent factors have been empasized (e.g.

Turchin (1993); Hornfeldt (1994)), which allows the population to grow for some time, but then responds by decreasing the population growth rate. Examples of delayed density dependent fac- tors include numerical responses of specialist predators, delayed regrowth of grazen vegetation, or time delays in maternal effects.

Although the reigning hypotheses and focus of research on population cycles in rodents have changed over time (see Stenseth and Ims (1993a)), some hypotheses have gained more credit than others. These hypotheses can have similar net effects on density dependence, and causal factors may be difficult to distinguish between. Nevertheless, the widespread and synchronized

changes in cycles indicate that factors that can change globally, and might do so in a synchro- nized fashion, may be a plausible explanation. In the following paragraphs, I present an overview of some of the hypotheses of the causes of the small rodent population cycles. I also present the- ories of how the dampening of the cyclic amplitudes or loss of cyclicity could be related to these hypotheses.

1. The predator hypothesis. In the mid 1920’s, Lotka (1925) and Volterra (1926) indepen- dently showed that a simplistic consumer-resource relationships of predators and prey can produce population oscillations. Although the ideas of Lotka and Volterra lacked realism (such as population self regulation), this important result spiked interest in the importance of predator and prey relationship. Predators are capable of producing both direct density dependence (Korpimäki and Norrdahl, 1991a,b) and delayed density dependence (Kor- pimäki, 1993) in prey populations. The direct density dependence is usually explained by generalist predators or nomadic specialists (see Hanski et al. (1991)), since these are known to respond almost immediately to prey densities (Korpimäki, 1994). Examples at Finse includes the avian predators rough-legged buzzard (Buteo lagopus), kestrel(Falco tinnunculus) and hawk owl (Surnia ulula). The delayed part of the hypothesis is explained by specialistic predators that respond numerically in a delayed fashion to rodent densities.

At Finse there are two mammalian specialist rodent predators: the least weasel (Mustela nivalis) and more abundantly the stoat (Mustela erminea). These predators are known to respond numerically to the density of the rodent community (Henttonen et al., 1987; Gilg et al., 2006), and absence are observed to coincide with prolonged peaks (Henttonen et al., 1987). The fact that cycles of sympatric rodent species are often synchronous, has often been taken a support of the predator hypothesis (e.g. Lack (1954)). This is because one can easily imagine that predators syncronize the rodent community by alternating prey when their preferred prey decreases profoundly, or by inducing proportional mortalities among all prey species (Korpimäki et al., 2005). See e.g. Hanski et al. (2001) for a review of the predator hypothesis.

2. Intrinsic factor hypothesis. Under this hypothesis some factor within the population is assumed to be the cause of the decline phase of the cycles. Proposed mechanisms include genetic changes (see Krebs (1978) and Myers and Krebs (1974)), stress levels (Christian, 1950; Boonstra, 1994), maternal effects (Ginzburg and Taneyhill, 1994), and sociality through mechanism such as dispersal or interspecies aggression (Charnov and Finerty, 1980; Andreassen et al., 2013). Theoretical models have shown that maternal effects can produce delayed density dependent effects capable of generating population cycles with

periods similar to those observed empirically (Inchausti and Ginzburg, 1998). In fact, the maternal hypothesis is one of the most supported hypotheses of intrinsic mechanisms that can cause cycles (Turchin and Hanski, 2001; Inchausti and Ginzburg, 2009).

3. The overgrazing hypothesis and nutrient cycling. It is clear that a range of theoretical models of plants and their interaction with rodents can produce both direct and delayed effects that produced cycles with a range of amplitudes and periods, similar to those ob- served empirically (Turchin and Batzli, 2001). The overgrazing hypothesis, have been given less credit of explaining the cycles of voles than for Norway lemming (Oksanen and Oksanen, 1992; Turchin and Batzli, 2001; Turchin et al., 2000), since experiments have shown that voles residing in habitats that are hypothesized to be of low quality do not experience a negative effect of this (Klemola et al., 2000). Further supporting this is the observation that Norway lemmings cause significant decreases in the vegetation of graminoids and mosses (Moen et al., 1993; Virtanen et al., 2002), while the decrease caused by voles might not be nearly as significant (Krebs and Myers, 1974, table XIII).

This could result in reduced quantity and quality of food in the following seasons (review in Batzli 1992). The delayed part of this hypothesis is the time is takes for plants to com- pensate by regrowth and fixation of nutrients. At Finse the growth season is very short which may further enhance the scarcity of resources and time delay. Another explanation for the delayed part of this hypothesis is the induction of secondary defense compounds in plants.

4. Climate and snow condition hypothesis. Climate fluctuations are known to affect the populations of a number of animals (Stenseth et al., 2002, 2003a). The rodents at Finse experience especially strong seasonality, and spend about two thirds of the year in the subnivean habitat. Changes in winter conditions can be expected to have strong effects on their survival. One pathway could be through increased frequencies of melt and freeze events, due to warmer climate. Such events are known to have detrimental effects on ro- dent survival (Aars and Ims, 2002). In 2008 Kausrud et al. showed that a combination of humidity, snow hardness, and the previous years rodent abundance, could predict the spring catch rates of rodents at Finse, emphasizing the correlation between these events.

Another pathway could be if changed conditions or length of the subnivean habitat medi- ates different interactions with predators (Saitoh et al., 2006). If for example, predators obtain access to the prey at an earlier date of the year, this could have detrimental ef- fects on the increase phase before the peaks. It has also been shown that season length itself is capable of generating population cycles, a phenomenon termed “seasonal forc-

ing” (Stenseth et al., 2003b). If this is the case, changed lengths of seasons due to climate changes would result in different dynamics.

5. Multi factor hypothesis. This idea has arisen from studies showing that experimental exlusion of single factors in cyclic populations results in cycles where the different phases are different fromt the emperically observed cycles. Instead, a combination of factors is hypothesized to have important effects through different parts of the cycles. For example, Radchuk et al. (2016) showed that by excluding predation, the cycles did not produce lows that last more than one year, although the cycles of many species have lows lasting for two years. Many scientists acknowledge that a multitude of factors are important for explaining the cycles, but research on rodent cycles have definitely been coloured by single hypotheses receiving enormous attention. As Chitty (1996) noted,

Unfortunately, in spite of 70-odd years of learned effort there’s still no answer to this ecological equivalent to the riddle of the sphinx. Plausible hypotheses do indeed abound, each championed by one of several quarrelsome cults; but in spite of evidence in their favor, all lack the hallmark of scientific proof (so called), namely, success in giving repeatable results. For until corroborated by risky predictions, favorable evidence must always be mistrusted, especially if it’s your own.

With these hypothesis in mind, I want to describe the dynamical properties of the populations.

For this I need estimates of abundance through the time series. Previous works have used a capture index (numbers of captures divided by number of traps multiplied by the number of hundred days) when making inferences about the population dynamics (Framstad et al., 1997;

Kausrud et al., 2008). This relies on the assumption that this index is proportional to the rodent abundance, and ignores important problems caused by the sampling procedure. I expect the es- timates to be biased because sampling removes animals from the local population, traps become saturated, and the presence of multiple species saturates the traps even faster. To handle these problems there is a need for a better observation model.

In this work I develop a hierarchically formulated Bayesian removal model, that uses the daily number of captures, data on sympatric species, and latent variables that impose structures that are hidden in the data, to obtain estimates of the abundance. This model is then extended to a fully integrated state space model, which models the dynamics through the time series. The hypothesis of changed dynamics from the interval 1970-1994 to 1995-2017 is tested by allowing

the model to have different sets of coefficients in each interval. The resulting changes in direct and delayed density dependence and estimated and underlying dynamics is discussed.

Chapter 2 Methods

The data used in this work was supplied by Erik Framstad from Norwegian Institute for Na- ture Research (NINA). It is a continuation of long-term snap-trap capture series of rodents at Hardangervidda at Finse, started by Eivind Østbye at UiO. Originally, this was as part of the studies in the Norwegian section of the International Biological Program (IBP) (1964-1974, in Norway 1968-1975). The focus of the IBP program was to study biological communities and changes in the natural environment, and on the conservation and growth of natural resources for human benefit. Since the mid 1960’s various studies have been conducted on several grids in Hardangervidda, as an effort to obtain various information on small mammal populations.

As such, the sampling protocol was not optimized for this work, and sources of uncertainty are included in the data.

2.1 Locations

The location of the sample grids are indicated in Figure 2.1, which shows a satelite photo of the Finse area. Figure 2.2 shows of the vegetation in of each area. Clearly, the two areas are qualitatively different, as discussed in more detail in the following section.

K1 60.36°’N, 7°32’E

Finse Alpine Research Center

T1 60°35’N, 7°31'E Finse train station

Figure 2.1: The sample sites K1 and T1 along with their gps-coordinates (latitude, longitude).

Figure 2.2: Grid locations: K1 on the left and T1 on the right.

2.2 Data gathering procedure

Rodents have been captured at two geographically separated grids “K1” and and “T1” (Figure 2.1, 2.2). The two areas are different in their primary productivity, and therefore food availabil- ity. The snow melts earlier at K1 than at T1, due to more sun exposure (south facing hill), and the growth season is longer in this location. This is reflected in the number of captures, as K1

1 2 3 4 5 6 7 8 9 10

A B C D E F G H I J

Figure 2.3: Trap setup. There are two traps (denoted by blue dots) at each trap station. The locations are given by the intersection of letters and numbers.

should be able to support many more nests than T1, although the captures are not always higher in this field. This suggests that the differences should be accounted for in the analysis, and the grids may not be considered as replications.

Each grid is the size of one hectare (100m·100m), and there are 100 trap stations evenly distributed along 10 transection lines (Figure 2.3). At each trap station, 2 traps were placed at spots where capture was deemed most probable within approximately 2 meters, such as spots with signs of activities or obvious nest spots (e.g. holes with piles of secrements, such as shown in photo to the right in Figure 2.4). The traps were baited with a piece of wick that was coated in cooking oil. The traps were set for a full day, and checked and emptied in the morning. In some years, and especially in the spring, a great portion of the trap-stations were covered with snow. No traps were placed at these stations, resulting in a lower sampling effort. Additionally, traps were sometimes lost to other causes, which could happen if for example a trap failed to kill a large Norway lemming, and was dragged with it into tunnels. Other causes could be that larger predators (e.g birds or mustelids) feeds on animals that are caught in the traps, and drags the animal along with the traps out of the grid. Missing traps were replaced on the following day. In roughly half of the years there were fewer than the normal 200 active traps, and in the worst years as much as half of the traps were covered in snow. This was accounted for in in the analysis by including the total number of active traps in the models.

Figure 2.4: Trap preparation and placement. Middle: applying cooking oil on wick as bait.

Right: Trap placed in front of an obvious nest spot (with piles of rodent excrements).

2.3 Data

The data has been collected into a excel spreadsheet with data on individual captures and their characteristics in columns. The time series table contains capture data on seven species:Lemmus lemmus, Microtus oeconomus, Microtus agrestis, Myodes glareolus (formerly Clethrionomys glareolus), Myodes rufocanus (formerly Clethrionomys rufocanus), Sorex araneus and Sorex minutus,along with variables for sex, weight, reproductive status, and age. The capturing pro- cedure has been consistent since 1980. The data from the 1970’s is less regular, as there were performed many field experiments and trapping proceeded over grids of varying size. However, for the analysis in this thesis, only data that was captured with the same procedure and on the same grids was used.

2.4 Formatting of the data

To fit the models introduced in Section 2.8, the data had to be formatted. This meant aggregating species specific captures on the day they were captured. Table 2.1 shows an excerpt of the data after formatting for Norway lemming. The tables contain the captures, with columns indicating the day the captures were made and rows indicating the years (spring and autumn, respectively).

Similar tables were made for all species. The formatted full dataset is included in Appendix A on page 101.

Table 2.1: Excerpt from the formatted data containing number of captures of Norway lemming on the productive grid K1. Table elements show the number of animals caught. The year and season is given as rows and the capture-night is given as columns

Day 1 Day 2 Day 3 Day 4 Day 5 Day 6

S 1986 9 0 2 0 0 1

A 1986 1 0 0 0 3 0

S 1987 1 0 0 1 0 0

A 1987 0 0 0 1

S 1988 15 8 6 6 10 6

A 1988 141 72 46 50 58 37

S 1989 0 0 0 0 0 0

A 1989 0 0 0 0 0 0

S 1990 0 0 0 0

A 1990 0 0 1

S 1991 6 5 2 3 1 0

A 1991 27 13 11 17 17 15

S 1992 1 0 0 0 1 0

A 1992 0 0 0 0 0 0

S 1993 0 1 0 0 0 1

A 1993 4 2 4 0 0 0

S 1994 8 11 7 16 12 6

A 1994 60 57 51 48 61 60

2.5 Problems with the data and handling

During the 1970’s the capture procedure was performed in a inconsistent way compared to the remainder of the years. In some of the peak years, traps were checked twice a day without clear notes on when it was done so, and when it was not. As a result the datafile contains captures on one more day than what is expected from the sampling protocol (traps that were emptied the same day as they were placed). This was handled by simply removing these captures from the data, and adding the number of individuals removed to the estimates of abundance. Other than this, it was not possible to account for between year differences in the procedure.

In two of the years, 1978 and 1979, it was noted very low presence of animals from field observations. As a result the sampling was not carried out as normally, with the reasoning that zero captures was expected. Instead they placed two transect lines at the study grids, with half the number of traps as were usually placed (100 instead of 200), and checked them daily for two days. In 1979, T1 was not surveyed at all. This is of course problematic, in terms of fitting an

model, since the sampling procedure is not the same. For K1, I handled this as two days with very little traps, but ignored that it was not on the exact grid. For T1, 1978 was entered as normal year with half the number of traps and only two days of survey. 1979 was entered as NA in the data (which results in large uncertainty in the estimates of abundance for this year).

2.6 Study species and environment

Table 2.2 shows the total number of captures of each species during the years 1970-2017. The most common species in the captures and the focal species of this thesis were Norway lemmings (Lemmus lemmus), root vole (Microtus oeconomus) and field vole (Microtus agrestis). These are rodents belonging to the familyCricetidaeand subfamilyArvicolinaecharacterized by a herbi- vorous diet and large population fluctuations. The Norway lemming is endemic to Scandinavia, and the Kola peninsula in Russia (Henttonen, 2016). Root vole has a holoarctic distribution (Linzey et al., 2016) and field vole has a paleartic distribution (Kryštufek et al., 2016). They both commonly occur alongside the Norway lemming. The Norway lemming feed on grasses, mosses and sedges (Tast, 1991; Saetnan et al., 2009; Soininen et al., 2013), while root and field vole uses a wider range of resources (Hansson, 1971; Saetnan et al., 2009).

The rodents have high reproductive rates under favourable conditions. The field vole often has several litters of 4 7 pups (Myllymäki, 1977) that can become sexually mature in 3 weeks.

The root vole often has several litters of 4 8 pups, and the pups can become sexually mature in three weeks (Ims, 1997). In laboratory conditions Norway lemming is known to have up to 16 pups and can become pregnant in 2 weeks (Semb-Johansson et al., 1993), but in the data used in this work the average number of fetuses found in pregnant females was 5.6, and the maximum was 13.

The climate at Finse is very harsh with long winters and short summers. The ground is usu- ally free of snow from July until the onset of snowfall, which begins early in October or (late) in November, when the formation of the subnivean habitat begins. As such, the inhabitant species at Finse are highly adapted to the strong seasonality. Under favourable conditions, the three ro- dent species can reproduce under the snow-cover throughout winter (Tast and Kaikusalo, 1976;

Myllymäki, 1977; Hansson, 1984). This overwintering reproduction seem to be especially im- portant for initiation of the increase phase of the Norway lemming cycles. Norway lemmings are famous for their migrations which occur in two phases. In early spring when extensive snowmelt begins, they migrate from their winter snow-bed habitats, often down from the mountain highs,

Table 2.2: Numbers caught of the different species on the capture grid K1 and T1 in the period 1970-2017

Species Captures on K1 Captures on T1

Lemmus lemmus 1751 1755

Microtus oeconomus 768 34

Microtus agrestis 132 26

Myodes glareolus 47 52

Sorex sp. (unidentified) 92 2

Sorex araneus 50 2

Myodes rufocanus 10 33

Sorex minitus 3 0

on to moist summer habitats (Henttonen and Kaikusalo, 1993). When the population reaches some threshold density, migration from their summer habitats begins, probably because good quality habitats are sparse. The spring migration is short, lasting two to three weeks, while the autumn migration is longer and can last up to three months (Henttonen and Kaikusalo, 1993).

Studies of genetic diversity that show small genetic differences also indicate that there is much dispersal of individuals between populations (Andersen and Wiig, 1982). These are testimonies of the large dispersal capacity of this species.

2.7 Modeling abundance

In this work I want to model the dynamics of Norway lemming, root vole and field vole. To do this I need some estimate of abundance throughout the time series. For animals where accurate observations are difficult, as is the case with rodents, reliable measures of abundance can be difficult to obtain. A major problem is how counts of the species are related to the number of animals present. How abundance should be estimated will depend on the sampling protocol. In this study, rodents were caught by sequential removal sampling.

Previous work on the same data used the number of rodents captured per 100 trap-nights (Kausrud et al., 2008; Framstad et al., 1997) as an index of the abundance, hereafter referred to as CE100 (captures per 100 trap-nights of effort). This index is commonly used for this type

of data (see for example Henttonen et al. (1987); Framstad et al. (1997); Hörnfeldt (2004)).

Here I refer to the trapping interval, from the moment the traps are set until the moment they are checked the next day, as a trap-night. This term was first used by Grinnell (1914), in his description of a capture procedure for animals with nocturnal habits, but has since become the convention for describing trapping intervals for small mammals. Although the CE100 index ought to be positively correlated with abundance, there is no obvious relation to the true number at the beginning of sampling, which is the quantity I am interested in. There are also factors in the data that are known to cause variation in this index, which motivates the use of different measure of abundance. In the following paragraph I list some of the problems associated with the index which I attempt to handle in this work:

1. Population depletion and unequal length of sampling schedule

An important problem with using CE100 as an index of abundance is that sampling depletes the population of animals. Each successful capture will reduce the population by one. Intensive sampling over a restricted area could substantially reduce the population, and the number of captures at later points in time is expected to be lower than what it was at beginning of sampling.

This could bias the CE100-index, as illustrated in the following example.

Say all animals have the same capture probabilityr, and the abundance before sampling is N₀. With the parameter values N₀=10, r =0.5 and F =5 traps, a sampling period could yield a vectorY= (2,1,1,3,3)of captures on consequtive trap-nights. The sum of captures is Â^I_i=1Y_i=10=N₀, which means that all animals have been caught. If captures continues for 3 more nights the resulting sequence is Y= (2,1,1,3,3,0,0,0), as there are no more animals available for capture. For the first case the CE100-index would be ¹⁰₂₅trap-nights^animals =0.4. For the second case it would be ¹⁰₄₀trap-nights^animals =0.25, effectively a downscaling just by increasing the number of trap-nights sampled. The problem could occur with the rodent data, since the captures were obtained over periods ranging from 3 to 6 trap-nights. Without accounting for the animals that have been removed the capture index could be substantially reduced, especially in years with low abundance and many trap-nights.

2. Trap saturation

Within trap-nights (between trap-checks), the availability of traps is depleted in the same way as the population of animals, an effect called “trap saturation”. As traps are sprung throughout the trap-night, the remaining number of traps decreases. Further captures are therefore produced

at a lower effort. I denote the number of active traps by F. In the extreme one can imagineF individuals being caught at some time throughout the trap-night, which leaves no active traps.

The probability of capturing animals has then decreased to zero. As this point is reached, the traps are fully saturated, and the precision of the estimates at very high abundances may be poor.

This is a clear limitation of the measurements of the density, and it needs to be implemented in the model. This problem has been noted by (Nelson and Clark, 1973; Xia and Boonstra, 1992; Beauvais and Buskirk, 1999) and ways of accounting for the effect have been suggested (Nelson and Clark, 1973; Beauvais and Buskirk, 1999). Hanski et al. (1994) argued that many studies are robust to this effect with the capture index being sufficiently linearly related to the abundance under moderate trap-saturation. However, trap saturation at Finse was high in many years. On each study site there were 200 traps, or less in years when snow covered trap stations.

The highest number of captures observed in one trap-night was 142, yielding considerable trap saturation. Figure 2.5 shows the mean trap saturation (over consecutive trap-nights) on the grid K1. The peak years can be identified in the figure as spikes in the trap saturation rates. During the peak years, trap saturation ranges from 10-40%, and for these values the CE100 can not safely be assumed to be linearly related to the abundance. For an indication an assessment of the magnitude of this problem with these levels of saturation, see Beauvais and Buskirk (1999).

In these years captures were not produced at equal effort, biasing the CE100-index downwards.

●

●

●●●●●

●

●

●

●●●●●

●

●

●●●●●

●

●

●

●●●●

●●

●●

●●●

●

●

●●●●

●

●

●

●●●

●

●

●●●●●●●

●●

●

●●●●

●●●

●

●

●

●

●●●●●●

●●

●

●

●

●

●●●●

●●

●

●●●●

●

●

1970 1975 1980 1985 1990 1995 2000 2005 2010 2015

0.000.050.100.150.200.250.300.35

Year

Mean (# captures / # traps)

Figure 2.5: Mean trap saturation over trap nights within seasons on the grid K1. Spring values are shown with red dots and autumn values in blue. Note how the peak years appear as peaks in the saturation values.

3. Multiple species saturates traps even faster

There are multiple species at the study site. If one catches an animal of species A, a trap becomes occupied, which reduces the effort of capturing species B. If many animals of species A are captured, there will be few traps left to catch species B, in effect lowering the effort substantially.

Using the CE100-index independently across species can therefore cause negative bias. This should be accounted for by incorporate the number of captures of all species in a model.

4. Migration

The capture index relies on the assumption of a closed population. The number of trap nights in a survey ranged from 3-8. During this time it was possible to move in and out of the sampling grid, which makes the assumption of a closed population questionable. Moreover, numerous reports have been made of the large-scale migrations of Norway lemmings during years with high density (Henttonen and Kaikusalo, 1993). Migration clearly violates the assumption of a

closed population, and I expect estimates to be strongly affected by this. However, attempts to incorporate migration proved unsuccessfull in this work (probably because the data does not contain enough information to estimate this with the structure I tried), and remained a source of uncertainty.

5. Handling zeros

Surveys with zero captures occurred often in the data set, and is common in capture series of rodents in general. When modeling dynamics it is common to model the rate of change on log scale, which causes trouble whenever there is a zero. The usual way of dealing with this is by adding an arbitrary constant which makes it possible to calculate the log of the estimate. The chosen constant can bias estimates, and therefore many authors advice against log transforming count data (e.g. O’hara and Kotze (2010)).

In the following sections, I propose a model which incorporates the effects that could result in problems when using the capture index.

2.8 State space model

To model the dynamics of the rodent populations, I use a hierarchical state space model. A state space model (SSM) decomposes the dynamics into a process model: g N_0,k|N_{0,k 1},yyy and observation model: f Y_k|N_0,k,qqq . Figure 2.6 shows a schematic diagram of a state space model. The process model g N_0,k|N_{0,k 1},yyy describes how the system behaves, given its parameters yyy. This process and its parameters is often what is of ecological interest. The observation model f Y_k|N_0,k,qqq describes how the observationsY_k have been retrieved from the state of the system at yeark. The SSM framework decomposes the variation into ecological process variation (the process model) and variation that is due to the data gathering procedure (the observation model). This framework is also attractive in that it allows us to break down the system into simpler parts which might be easier to model.

2.9 Fitting the models (Bayesian paradigm)

To fit the models I adopt a hierarchical Bayesian perspective. In the Bayesian framework pa- rameters are summarized by their posterior distributions. The posterior distribution is the distri-

N_0,k-1

Y_1,k-1 N_0,k

Y_1,k

N_0,k+1

Y_1,k+1 N_0,k+2 N_0,k-2

Y_1,k+2 Y_1,k-2

Process model

Observation model

Figure 2.6: Graphical representation of a state space model. N_0,k is the unobserved abundance of the population in yearkat the beginning of sampling.Y_i,kis the observed number of captures the end of trap-nightiin yeark. Note that season (p) is left out of this Figure. becuase of this, the actual transitions is more complicated than more complicated than in this figure.

bution of the parameters given the data. It can be obtained from Bayes theorem:

P(qqq |data)µP(data|qqq)P(qqq).

This says that the product of the likelihood of the dataP(data|qqq)given the set of parameters qqq, multiplied by the prior distribution of these parametersP(qqq), forms the posterior distribution P(qqq |data). From the posterior distributions one can observe the variation of each parameter through the quantiles of the posterior distribution, giving us credibility intervals. This is useful for example when a survey yields zero captures, as more trap-nights with zero captures will shrink the posterior distribution closer towards zero, reflecting that there is more evidence for lower abundance. In this way the posterior distribution “learns“ from the data, and make better use of it. In the same situation the CE100-index is zero, and does not distinguish between surveys where you catch e.g. zero animals in one trap-night or zero animals over six trap-nights, even though the latter yields much more information.

Estimation of posterior distributions was done by Markov Chain Monte Carlo (MCMC). The model was fitted in JAGS software (Plummer, 2003) by using the the rjags-package (Plummer, 2016) in the R computing environment (RDevelopment Core R Core Team, 2016). JAGS is designed for implementing MCMCs through the Gibbs sampling algorithm. A Gibbs sampler works by starting at some initial values in the set of parametersqqq, and then sampling from the conditional distribution of the other parameters. It runs sequentially through the entire param- eter space qqq given all previously sampled parameters. It is therefore sufficient to define the conditional distributions of the parameters in order to define a model in JAGS.

Due to the extensive number of random variables in the model, model runs were performed on the Abel Cluster ("The Abel Computer Cluster", 2012), owned by the University of Oslo and the Norwegian metacenter for High Performance Computing (NOTUR), and operated by the Department for Research Computing at USIT, the University of Oslo IT-department (IT- department, ). To fit the model in JAGS, a modified version of JAGS was needed where the sum sampler in the bugs module was slightly modified. This version was provided upon request by Martyn Plummer.

In the next section I explain the observation model by specifying the conditional distribution of the parameters, moving through Figure 2.7 from the bottom to the top.

2.10 The observation model

N_0,s β_s

α_s

t

Y_1,s

Latent Variables

Observations F

Constant

Top level parameter

o Process model

Figure 2.7: Schematic view of the relations between the variables in the observation model.

Further details of the figure is explained in the introduction of Section 2.10. I label the variable b “top level parameter”, because this is the only variable for which one must specify a prior distribution. The figure illustrates the relations on the first trap-night of sampling, since the process model produces the initial abundance (N_0,s). The following trap-nights of sampling, e.g.

i=2,3. . .I, could be represented by similar figures by replacingN_0,s withNi,s, andY_1,s byYi,s, and removing the arrow from the process modelN_i,s.

A critical part of any state space model, is the observation model, f(Y_k|N_k,qqq). This describes how the observations, in this case the number of captures, are related to the unobserved abun- dance of animalsN_i,s. The primary goal in formulating this observation model was to account for the problems associated with the capture index (see Section 2.7). To do this I used a hierar- chically formulated removal model. Figure 2.7 shows a schematic diagram of this model. Here I assume an process model (which is specified in Section 2.11), that yields unobserved abundances N_0,s, as shown in the upper middle of Figure 2.7. Additionally, I have the constant of propor- tionalityb_s which represents the efficiency of the traps and behavior of the species. Given the

unobserved true abundances (N_0,s), the number of trapsFand a constant of proportionalitybs, it is possible to define the hazard rateaj,s, here represented by its vectoraaas= (a1,s,a2,s, . . . ,aYs,s).

The hazard rateaj,s, describes the rate that captures of speciessare arriving at, in the interval of time after event j 1, and the variables on the levels below the hazard vector are given by these events. These are the inter-arrival times (t_j) and species identification variable (o_j), which are represented by their vectorst= (t₁,t₂, . . . ,t_J)ando= (o₁,o₂, . . . ,o_J)in Figure 2.7. The final level contains the observationsY_i,s (in the Figure represented by the captures on the first trap- night (Y₁)), which depends on the vector of inter-arrival timest, and the order vectorooo. To fit the model in JAGS, one only needs to define the conditional distribution of all lower level variables, given the variables at higher levels of the hierarchy. In the following sections I describe how this model accounts for the problems that was mentioned in Section 2.7, and specify the conditional distributions of each of the variables.

2.10.1 Model specification

To avoid messy notation (multiple indices on each variable), the observation model is explained in terms of a single capture survey within a single year. The data does, however, consider two surveys, one in spring and one in autumn, over 48 years. The generalization of the model is done by adding an indexkfor years, and pfor season. These indices are included in the model (see Appendix D.1), but are omitted here for clarity.

2.10.1.1 The observations

The bottom level in Figure 2.7 the hierarchy is the number of captures of each species (Ys) (Figure 2.7). The total number of captures in trap-night is

J=

Â

Y_s,

whereSthe number of species used in the model (in this case 3). J is assumed to depend on the unobserved times of the captures. I denote the time of the j’th capture byhj, with h0denoting the time sampling begins. I assume that no captures arrive at the exact same time, ensuring that h_j is unique for each capture. The difference between the times of capturet_j=h_j h_{j 1} is

called the inter-arrival times. I denote the vector of inter-arrival times byt, so that t= (h₁ h₀,h₂ h₁, . . . ,h_J h_{J 1}) = (t₁,t₂, . . . ,t_J).

With this setup the sum over this vector Â^J

j=1t_j, is the time when the last rodent is caught during this trap-night. The number of captures therefore only depends on howmanyof these times that sum to less than one trap-night. If for example ²⁰Â

j=1tj=0.97 and ²¹Â

j=1tj=1.04, then there would 20 captures that trap-night, as the 21’th would arrive in the next trap-night. In other words, the times must satisfy the conditions Â^J

j=1t_j<1 and^J+1Â

j=1t_j>1. It might seem odd to introduce these times at first, especially since they are not observed, but the motivation for including them will explained in the next subsubsection. In JAGS the constraint on the times is specified by:

sum.t <- sum(t[1:(J+1)]) tlast <- t[J+1]

#Condition 1 - sum of times 1:J should be within 1 trap-night

C1 <- step(1-(sum.t-tlast))

#Condition 2 - sum of times 1:J+1 should exceed 1 trap-night

C2 <- step(sum.t-1)

C <- C1*C2 #Combined conditions

ones ~ dbern(C) #Bernoulli trick.

This utilizes what is called a “Bernoulli trick” which is explained in detail in Appendix A.1. It ensures that proposal values of times which do not satisfy the condition are rejected, by having a likelihood contribution that is zero. They will therefore be resampled until they satisfy the conditions.

2.10.1.2 Inter-arrival times

The sampling procedure was treated as a single process giving rise to “events”, which in this setting means that a rodent is captured. The time until an event depends only on the rate that events are arriving at. This is called the “hazard rate” and is defined as,

a(t) = lim

Dt!0P(tT <t+Dt |T t).

I assume that the hazard rate is constant in the time between captures, but depends on species considered. To capture how the intensity of the interaction between traps and animals changes over times, I consider the hazard rate to be a product of the abundance of the speciesN0,s, the number of active trapsF and a constant of proportionalitybs, which is related to the efficiency of the trap and behavior of the species. The initial hazard rate is then:

a_1,s=N_0,s·F·b_s.

When a trap is triggered there will be one less animal and one less trap, meaning that an event decreases the hazard rate until the next capture. Let Nc_j,s (“Number captured”) denote the number of captures of speciessbefore event j. The total number of captures after capture event

j 1 and before capture event jis then:

Nc_j.=

Â

Nc_j,s.

It follows that Nc_j. is also the number of traps that have become incactive (sprung). At this point, the remaining number of animals of speciessisN_0,s Nc_j,s, and the remaining number of active traps isF Nc_j.. This yields the general expression for hazard rate for speciess:

aj,s= N_0,s Nc_j,s · F Nc_j. ·bs.

Plots of this hazard function is included in Figure A.1 in Appendix A. Since I consider a multi- species setting, the rate that events are arriving at depends on the total hazard, which is obtained by summing over species.

aj.=

Â

s=1aj,s.

Under these assumptions the process producing the captures is a piecewise constant hazard that is monotonously decreasing by each capture. This is also the motivation for including the inter-arrival times as unobserved random variables. By their inclusion, I model the process of a decreasing capture effort, and hence trap saturation, and in addition the lower probability of captures that is associated with population depletion. In Appendix A.3, I show that the given N0,sandbs, the times are exponentially distributed with rate equal the total hazard: