In davan690/beech-publication-wr: Small mammal dynamics in NZ beech forests

Methods {#method}

# number of nights of each in [j,t] matrix
#source("./R/r-packages-needed.R")
#source("./R/davidson_2019_theme.r")

# number of nights of each in [j,t] matrix
#source("./R/wrangling/n-nights-code.R")

#saving total number of nights during study
#n.nights.df.1 <- table(filter(n.nights.df,Var2 == "hol R2")$value)

We carried out our study in beech forests in the South Island of New Zealand, in two adjacent valleys: the Eglinton Valley (GPS co-ord, S $44^\circ$ $45'$ S $168^\circ$ $10'$ E $m.a.s.l$) and the Hollyford Valley (GPS co-ord, S $44^\circ$ $2'$ S $168^\circ$ $5'$ E $m.a.s.l$). Stoats were the main predator and mice were the most abundant rodent species [@king1983]. Ship rats (Rattus rattus; @innes2005) and kiore (Rattus exulans; @roberts1991) were also present but at much lower numbers (@king1983; @ruscoe2001). Intensive stoat control has been under taken in the Eglinton Valley since 1983 and no stoat control had been conducted prior our study in the Hollyford Valley (pers comms).

We tested whether reducing predator populations by stoat control influenced mouse population dynamics as shown in Figure \@ref(fig:figure-one-plot1), and whether the outcomes were affected by interactions between mice and rats. Specifically, we tested the following five predictions (@blackwell2001; @blackwell2003; @tompkins2006; and @tompkins2013; Figure \@ref(fig:figure-one-plot1): Prediction A to D). If stoat control affects mouse population dynamics then, relative to sites with stoat control (Figure \@ref(fig:figure-one-plot1); hollow symbols), mouse populations at sites without stoat control should exhibit changes (red arrows in Figure \@ref(fig:figure-one-plot1)).

A) Higher mouse abundance in non-mast years; B) Higher peak mouse abundance in mast years; C) A faster rate of mouse population increase in response to high seed availability in late summer/winter; D) A slower rate of mouse population decline from peak abundance; E) Predictions A-D should hold only when both stoats and rats are controlled.

Data collection

To test Predictions A-E, we use data from four trapping grids, in each valley, with each grid located at least xx kms apart. On each grid, we monitored the abundance of rodent species (mice, rats and kiore) using live capture Elliot traps (Landcare Research Animal Ethics Approval number 02/99). Each grid comprised of $81$ traps, each $20$ metres apart in a $9 \times 9$ grid covering a total area of $2.56$ ha ($25600$ $m^{2}$).

We estimated the abundance of rodents on each grid using capture-mark-recapture techniques (CMR), assuming a closed population during each trip, which typically lasted five nights. Traps were set each night and all rodents captured were marked with ear tags so that individuals could be uniquely identified. Mice had a single unique ear-tag, and each individual was ear notched with a notch position that identified the season of first capture so we could identify animals that had been captured but had lost ear tags. Rats were double ear tagged. We sexed and weighed each captured animal. We aimed to trap for five nights during each trapping session (trip) but due to unfavourable weather, some trips were shorter.

To measure food availability, we recorded the amount of seed reaching the forest floor (seedfall) by placing four standard circular seed-traps on each grid. Each seed trap had a radius of $0.125$ metres ($m$)equating to an area of $0.0491 m^2$ metres and a overall sampled area on each grid of $0.1964 m^2$. These seed traps were in place for the duration of the study, with each trap elevated $1.2$ metres above the ground and covered by wire netting to exclude seed predators. The contents of each seed-trap were collected once during each trip and the contents were sorted to remove unwanted plant matter (e.g. leaf litter). We counted beech seeds, recording only those with a kernel because rodents are known to only consume the kernel of beech seeds [@ruscoe2005].

Rats could influence the population dynamics of mice through overlapping diets [@king2005; @mcqueen2008], or direct predation [@bridgman2013]. During our study we captured rats (ship rats and Kiore combined) less often than mice. Due to the low recaptures in this dataset we regarded it as not complete enough to fit a full capture-recapture analysis. As a less precise approximation we used the number of unique rats captured on each grid during each trip (the minimum number alive) as an index of rat abundance ($R_{j,t}$).

Analysis of the capture-mark-recapture data

We used a capture-mark-recapture model (CMR) to estimate the number of mice $(N_{j,t})$ on each grid ($j$) during each trip ($t$), assuming the numbers of mice on each grid were independent of each other and that the populations on each grid were closed for the duration of each trip. We had data on the number of times each mouse, $i$, was captured on each grid, during each trip ($y_{i,j,t}$). We used the method described in @royle2007 to estimate the true abundance of mice from the CMR data using data augmentation [@royle2009; @royle2012]. The benefit of this approach allowed us to combine datasets from different sites and to build study specific predictions/tests. To do this we augmented our data for each individual mice captured on each grid during each trip with pseudo-individuals. Each pseudo-individual represents a mouse that were potentially present on the grid but were not captured during our study. For these pseudo-individuals, the number of captures $y_{i,j,t}$ was set to zero. Each individual mouse, including the pseudo-individuals, was assigned an indicator variable $z_{i,j,t}$ taking the value 1 if the individual was present on the grid and 0 otherwise. For mice that were captured, $z_{i,j,t}$ was known and takes the value 1. For the pseudo-individuals, the value of $z_{i,j,t}$ is unknown and was estimated by the model as being drawn from a Bernoulli distribution with overall probability of an individual mouse being present on a grid during a trip as $psi_{j,t}$:

$$z_{i,j,t} \sim Bernoulli(psi_{j,t})$$

The probability that an individual mouse (including pseudo-individuals) was captured $y_{i,j,t}$ times was drawn from a binomial distribution with probability, $\mu_{i,j,t}$, conditional on the number of nights trapping on each grid during each trip $J_{j,t}$.

$$y_{i,j,t} \sim Binomial \left (\mu_{i,j,t},J_{j,t} \right )$$

$\mu_{i,j,t}$ is the probability an individual was captured during one night of trapping on a grid. This probability depends on whether the individual was present on the grid or not ($z_{i,j,t}$) and the probability of capture conditional on presence ($p_{i,j,t}$).

$$\mu_{i,j,t} = z_{i,j,t} \times p_{i,j,t}$$

We specified the number of pseudo-individuals ($m_{j,t}$) to be $4 \times$ the total number of mice captured on each grid during each trip ($N_{j,t}$), which should be sufficient to ensure a non-informative prior for the number of mice on each grid during each trip [@royle2012a; @ruscoe2011]. We anticipated that the capture probabilities could vary between grids and trips, and there is often heterogeneity among individuals in capture probability for small mammals [@krebs2011]. We therefore modelled variation in individual capture probability on the logit scale, allowing for heterogeneity by assuming values were drawn from a normal distribution with a different mean for each trip and grid ($\mu.ind_{i,j,t}$). and a variance estimated from the data, which reflected unobserved among-individual heterogeneity. The mean capture probabilities for each trip and grid ($\mu.ind_{j,t}$) were modelled hierarchically, treating them as draws from a normal distribution with overall mean and variance estimated from the data.

$$logit(p_{i,j,t}) = lp_{i,j,t}$$

$$lp_{i,j,t} \sim Normal(\mu.ind_{j,t},\tau.ind)$$

$$\mu.ind_{j,t} \sim Normal(\mu_p,\tau_p)$$

We estimated the number of mice on each grid during each trip $N_{j,t}$ by summing the values of $z_{i,j,k}$.

$$N_{j,t} = \sum{z_{i,j,t}}$$

Seedfall

Our study aimed to determine the relationship between stoat control and mice populations. Differences in food availability between populations (grids) was collected using a standard seed trap methodology for beech forests [@choquenot2000; @ruscoe2001; @blackwell2001; @blackwell2003; @ruscoe2005]. Within each population, an individual mouse can only consume a maximum amount of seed in any given time period between t and t+1. To incorporate this individual "intake rate"($S_{j,t}$) we fitted three relationships previously proposed in the literature to our observed field data (seedfall);

1) a linear relationship between seed availability and $r_{j,t} = \beta_0 + \beta_1(Seed_{j,t})$;

2) a decelerating relationship similar to a Type II response $r_{j,t}=\beta_0+\beta_1(log(Seed_{j,t}))$;

3) a relationship based on a Type II functional response $r_{j,t}=1042.1\left(1-e^{-(Seed_{j,t} \times 0.00139)}\right)$.

The data used to fit our models were collected in the same valley (Eglinton Valley) as @ruscoe2005. @ruscoe2005 fitted a wide range of functional responses to laboratory data harvested from the same valley as rates of increase in mice populations ($r_{j,t}$). They concluded that the IR (type II functional response) proposed by @choquenot2000 and extended in @ruscoe2005 was the best fitting functional response from experimental data. We used the transformed variable from the best fitting model above to incorporate individual "intake rate" $S_{j,t}$. After fitting the curve between seedfall and mouse abundance in our data we selected the best fitting model using AIC [@akalie1964]. The intake rate we used was $S_{j,t}=log(Seed_{j,t})$.

Other demographic and environmental processes are also known to regulated mouse populations in New Zealand beech forests [@ruscoe2005]. Population density has been identified as a key driver of mouse abundance in multiple studies in similar forest systems to this research (@choquenot2000; @holland2015; @ruscoe2005). In previous research it was still not clear what exact density dependent processes were acting on mouse populations, however, it is still regarded that they are an important component of mice dynamics [@holland2015]. In our model we included the estimate of mouse abundance ($N_{j,t}$) at the beginning of each time step as a crude measure of density dependence.

Treatments

We aimed to test our predictions about the effects of stoat control by comparing mouse population dynamics in the two valleys. However, because we lacked replication at the valley-level, we could not be certain that differences in mouse population dynamics between valleys was due only to the presence or absence of stoat control. To overcome this, we undertook stoat control in the Hollyford Valley commencing in May 2002, allowing us to examine mouse dynamics in the presence or absence of stoat control in the same valley in addition to the between-valley comparison.

Stoat treatment

Stoats were removed in the Hollyford Valley using 13 Fenn Traps (are they following the road wendy? On map?). The original treatments applied to the eight grids (two valleys; one with stoat control, one without, and; four grids in each valley, two with rat control and two without) were maintained from May 1999 to May 2002. From November 2002 until the end of the study (August 2004), we reduced the number of trapping grids to three in each valley.

Rat treatment

We experimentally reduced rat (kiore and ship rats collectively) densities on two randomly selected grids in each valley by removing and humanely killing all captured rats on those grids. On the remaining two grids in each valley, we marked and released rats, as we did for mice on all grids. We undertook a total of $20$ trips to the two valleys, with each trip spaced three months apart from May 1999 to February 2004. The timing of trips during the year corresponded to each of the four seasons (February = Summer, May = Autumn, August = Winter, and November = Spring).

Statistical tests

Prediction A and B were tested using a mixed-model structure (GLM) to account for heterogeneity in estimates. We accounted for the spatial heterogeneity between valleys (Eglinton and Hollyford Valleys) and our rat removal conditions where, on some grids, rat numbers were artificially reduced).

Mice are known to be under different processes during each of the four seasons of the year. We added season as a random co-variate to test prediction A to allow for the effect of varying seed by taking the relative "intake rate". We assessed model fit using AIC and assessed model fit using visual assessment of residual diagnostics to check the model assumptions were meet.

Prediction A and B

#results from bayesian model
# source("./R/Davidson_2019_Data_wrangling2.R")
plot.dat.all1 <- read_csv("./data/plot-all-data1.csv")
# names(dat_all_output)

# head(abund.dat5)
low.abund.dat <- plot.dat.all1 %>%
  # mutate(N = log(N) %>%  # for when using out.N
  filter(trip == 10 | trip == 11 | trip == 12 | trip == 13 | trip == 14 | trip == 15)
# filter(low.abund.dat, trip.no == 15)$true.date

We tested if stoat control would resulted in more mice during non-mast years when mouse abundance was low (Prediction A) using equation below. We did not account for varying seed input because seed fall was biological equivalent to zero for all grids and trips during low abundance seasons [@choquenot2000].

$$ N_{j} \sim \beta_0 + \beta_{1}(Valley_j) + \beta_{2}(Control_j) + \beta_{3}(Rat_j) $$

We compared the average mouse abundance from our CMR model in each grid, grouped by our three treatments ($Valley_j, Control_j, Rat_j$. We only selected grids and trips in low seedfall years (Figure \@ref(fig:figure-one-plot1; label A) for this test. We used the replicates in these non-mast years (n = r length(low.abund.dat$N)) to compare the mean differences in areas with and without stoat control using a two-way analysis of variance (ANOVA). We log transformed mouse abundance for the general parametric assumption of normality of this simple test and to bound the results to above one (count data).

We tested if stoat control resulted in higher peak mouse abundance (Prediction B) after accounting for the fact that peak abundance on grids would be partly driven by seedfall [@king1983]. We therefore compared the number of mice per seed ($\frac{S_{j}}{N_{j}}$) among grids with and without stoat control using the same ANOVA test as Prediction A.

$log(\frac{N_{j}}{S_{j}}) \sim \beta_0 + \beta_{1}(Valley_j) + \beta_{2}(Control_j) + \beta_{3}(Rat_j)$

Prediction C and D

To test prediction C and D (increasing and decreasing phases of mouse dynamics) we needed to incorporate the other factors known to affect mouse populations. We estimated ($r_{j,t}$) by grouping replicates by Stoat Control(S), Valley(V) and Rat removal conditions(C). we used the complete data likelihood approach described by @schofield2014 (for details see supplementary material) to propagate the uncertainty in our abundance estimates. We then tested predictions C and D we compared rates of increase with and without stoat control using the associated model coefficients.

To test Prediction C and D we modelled the rate of increase of mice ($r_{j,t}$) between successive seasons (Figure \@ref(fig:figure-one-plot1)) with rate of increase calculated as:

$$r_{j,t} = log\left (\frac {N_{j,t + 1}}{N_{j,t}} \right )$$

We identified the most appropriate numerical response to test with our predictions and model system.

$$r_{j,t} = \beta_{0,S,V,C} + \beta_{1,S,V,C} (S_{j,t})+ \beta_{2,S,V,C} (N_{j,t}) + \beta_{3,S,V,C} (R_{j,t})$$

To test Prediction C and D respectively (does stoat control increase the rate at which mouse populations initially increasing/decreasing in response to beech seed intake) we compared the population model intercept ($\beta_{i}$) between grids, with and without stoat control. Where $\beta_{0,S,V,C}$ (model intercept) is equivalent to $r_{j,t}$ when all other parameters are set to mean for each group (treatment); $\beta_{1,S,V,C}$ is the adjusted^[1] effect of a intake rate ($S_{j,t}$) on $r_{j,t}$; $\beta_{2,S,V,C}$ is the adjusted effect of a single unit change in mouse density at the beginning of the previous season $N_{j,t-1}$ and $\beta_{3,S,V,C}$ is the adjusted effect of a single unit change of the minimum number of rats alive ($R_{j,t}$) at the beginning of the previous season $R_{j,t}$. If mesopredator release occurs in the increasing seasons; the intercept ($\beta_0$) for the areas with stoat control will be greater than the uncontrolled areas.

Prediction E

We tested if $r_{j,t}$ only occurs when rat populations are reduced to low densities. (Prediction E). We compared differences between the model coefficient for rats ($\beta_3$) grouped into the grids where we reduced rat densities and the others that where not. We expected that if rat interactions are sufficiently high to impact mice populations but only when stoats are at low levels would result in differences between these parameter estimates with replicates (trips and grids) showing differing population dynamics.

Overall priors

Non-informative prior distributions were selected to allow the data to drive parameter estimation [@gelman2006]. The prior distributions of capture histories were assigned binomial distributions and variance terms were assigned broad uniform prior (0-100) as suggested by @gelman2006. The remaining parameters were assigned normally distributed prior with mean 0 and variance 0.00001. Each model was run for 100 000 iterations with a burn-in of 50 000 iterations, which was sufficient to achieve convergence as judged by visual inspection of the chain histories [@bronder]. For a full model description see Appendix.

Software

Models were fitted in a Bayesian framework using Markov Chain Monte Carlo (MCMC) methods and implemented in JAGS [@plummer2011], called from R v.4.3.4 using the jagsUI package [@kellner2018]. Software pipeline used to reproduce this manuscript was generated using workflowr, rmarkdown using a tidyverse approach and the following additional packages (cite all other packages?)

davan690/beech-publication-wr documentation built on March 29, 2020, 11:09 a.m.