In davan690/beech-publication-wr: Small mammal dynamics in NZ beech forests
Summary {#summary}
[still needs work but last]
The conflicting outcomes from previous studies provides the motivation for this more comprehensive field study assessing the effects of stoat control on mouse dynamics. Here we clarify the discrepancy between research suggesting that mesopredator release of rodents will occur in New Zealand forest systems and field studies that have presented limited but conflicting support for increases in mouse abundance following predator control. We used an experimental design to test for differences between areas with and without stoat control after accounting for known effects of resources, density and competition on mouse dynamics in NZ beech forests. We found no evidence to suggest mice will become more abundant after predator removal.
# libraries needed
source("./R/r-packages-needed.R", echo = FALSE)
require("dplyr")
require("knitr")
require("docxtools")
require("citr")
#document global rules
knitr::opts_chunk$set(comment=NA,
# fig.path = "../figs/",
echo=FALSE,
fig.height=5,
fig.width=7,
message=FALSE,
warning=FALSE,
tidy = TRUE
# , cache.extra = packageVersion('tufte')
)
# options(knitr.table.format = 'markdown')
# options(htmltools.dir.version = FALSE)
# how do I do this??
# ,eval = FALSE,include = FALSE
# export .r code only
# knitr::purl("./Davidson_2019_BeechForest.Rmd")
# render draft to webpage
# rmarkdown::render(input = "Davidson_2019_BeechForest.Rmd")
# ,
# output_format = "html_document",
# output_file = "Davidson_2019_t.html")
## ----figure themes------------
source("./R/theme_raw_fig3s.r", echo = FALSE)
source("./R/davidson_2019_theme.r", echo = FALSE)
## ----sorting-references--------------------------------------------------
# citr::md_cite("./Beech-forest.bib")
# bib(file = "Beech-forests.bib")
# automatically create a bib database for R packages
knitr::write_bib(c(.packages(), 'bookdown', 'knitr', 'rmarkdown','tidyverse', "compareGroups", "jagsUI"
), 'packages.bib')
#data and plotting import
# source("./R/r-packages-needed.R")
# source("./R/theme_raw_fig3s.r")
# source("./R/davidson_2019_theme.r")
# # source("./R/Davidson_2019_Simple_models.R")
# source("./R/Davidson_2019_Model_wrangling.R")
# source("./R/Davidson_2019_Data_wrangling.R")
# source("./R/figures/pb-source-code.R")
# source("./R/figures/pa-plot-code.R")
# source("./R/figures/pb-plot-code.R")
#
# source("./R/data_input.R")
# source("./R/figures/pd-plot-code.R")
#
# source("./R/figures/pd-plot-code.R")
# source("./R/figures/rate-density-plot.R")
# source("./R/figures/pa-plot-code.R")
bibliography: ["Beech-forests.bib", "references.bib"]
Introduction {#intro}
Worldwide, but particularly on islands, introduced predators have significant impacts on native species [@gurevitch2004], and in the worst cases extinctions [@doherty2016; @towns2011]. For example, on the Island of Guam, the invasion of the Brown Tree Snake (Boiga irregularis) was the primary driver of multiple extinction events [@fritts1998]. In New Zealand, mammalian predators contribute to a disprorpoitatly large group of invaders impacting native populations nationally [cite]. However, removing, or reducing the abundance of the top predator can lead to increasing numbers of mesopredators (for a review see @prugh2009 and for examples see @rayner2009; @robles2002), which in turn, can lead to unintended and suprising outcomes [cite suprise island papers]. For example foxes predate on native and evolutionaryly distinct marsupials throughout Australia [cite margaritas paper??], heavyily impacting the full range of many native species [@kinnear1998].
Introduced species also have destructive impacts in New Zealand’s (NZ) remaining/remnant native forests. In these systems direct predation is commonly proposed as the reason for native species decline [@fitts1998; @kinnear1998]. These contemporary ecological communities commonly contain of at least one of the following four introduced mammalian predators;
- Stoats (Mustela erminea; Veale (2015));
- Brushtail possums (Trichosurus vulpecula Kerr; Clout and Gaze (1984));
- Ship rats (Rattus rattus; Innes (2005) referred to as rats);
- House mice (Mus mus musculus; @allen2006 referred to as mice).
Stoats are the generally regarded as the top predator when they are present in New Zealand ecosystems [@king1983], following their deliberate introduction in the late nineteenth century [@king2017]. Stoat control is commonly undertaken in these systems to protect native birds that are vulnerable to mammalian predation [@white2006], in particular, hole-nesting species like mohua (Mohoua ochrocephala; @Odonnell2006]. Stoat control is now commonly undertaken to protect native birds that are vulnerable to predation [@white2006], in particular, hole-nesting species like mohua (Mohoua ochrocephala; [@odonnell2006]). However, the primary food source for stoats in NZ forests are rats and mice [@king2006], and consequently there is a concern that reducing stoat populations to protect native birds may allow rodent populations to increase (e.g @rayner2009). An increase in the number of rats and mice could offset the benefits of stoat control because rodents are known to consume the eggs and chicks of native birds [@russell2016; @towns2006; @latham2018], directly compete with native species for food resources such as flowers and seeds [@mcqueen2008] and predate on invertebrates [@ruscoe2012]. In this paper we address the question: does stoat control lead to increased abundance of rodents, particularly mice, in NZ beech forests?
####UP TO HERE>>>>
We carried out our study in beech forests in the South Island of New Zealand where stoats are the main predator and mice are the most abundant rodent species [@king1983]. Ship rats (Rattus rattus; @innes2005) and kiore (Rattus exulans; @roberts1991) are also present but at much lower numbers [@jones2011]. We tested whether reducing stoat populations by predator control influenced mouse population dynamics as shown in Figure 1, 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 1: Prediction A to D). If stoat control affects mouse population dynamics then, relative to sites with stoat control (Figure 1; hollow symbols), mouse populations at sites without stoat control should exhibit changes in abundance in the direction of red arrow in Figure 1:
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.
Studies elsewhere in the world have shown that removing or reducing the abundance of a top predator often leads to an increase in the numbers of predators at lower trophic levels (termed mesopredator release), which in turn, can lead to unintended and often negative outcomes for native species (for a review see @prugh2009 and for examples see @rayner2007; @robles2002). While mesopredator release has been widely documented elsewhere, it is unclear if stoat control in NZ forests will cause rodent populations to increase. Rodent populations in NZ forests are known to respond strongly to variation in food supply [@choquenot2000; @ruscoe2001; @blackwell2001; @blackwell2003; @ruscoe2005; @tompkins2006; @tompkins2013; @holland2015; @latham2017], primarily seed availability (Figure 1).
# #data file
# source("./R/ecosystem-simulation/sim-raw-data.R", echo = FALSE)
#file code
source("./R/figure_1_script.R")
result.plot
cap1.stes <- c("*Expected changes in mouse populations over time in New Zealand forests (bottom panel) in response to changes in seed availability (top panel) during non-mast and mast years (grey blocks). In the bottom panel, differences between the hollow and filled symbols represent the proposed differences between areas with and without stoat control respectively. The arrows labelled A-D represent the four predicted outcomes of stoat removal that we tested; A) during non-mast years when little seed is available we expect that mice populations will be similar, B) at the peak of mouse abundance (during winter and spring in mast years) when we would expect larger mouse populations in areas where stoats are no longer predating on mice, C) during increased seed availability in areas with stoat control would expect to be able to support larger mouse populations; D) when mouse populations are declining from peak abundance, stoats would additionally reduce mouse abundance through predation in addition to other density dependent processes.*")
This is particularly pronounced in beech forests (spp. nothofagus) throughout New Zealand, where, between years beech seed production is highly variable. Little seed is produced in most years (Figure 1a: non-mast years) with occasional years of high seed production (Figure 1a: mast years; grey boxes). Mouse populations are low in the non-mast years, due to low food availability [@choquenot2000; @king1983]. In mast years, when seed becomes abundant, mouse populations can increase quickly following a predictable seasonal cycle. Seed begins to fall and accumulate on the forest floor in late summer (February) allowing mouse populations to increase, with mouse populations typically remaining high through winter (August) and into the following spring (November). Beech seed that is not consumed by mice and other seed predators germinates in spring to early summer, meaning this food resource disappears and mouse populations begin to decline in the subsequent seasons. If the following year is a non-mast year with little seed available, mouse populations fall to low levels. It is unclear whether stoat populations can increase at a great enough rate to exert sufficiently strong predation pressure to alter these food-driven population eruptions.
Previous studies have investigated the likely response of mouse populations to stoat control by modelling the outcome of interactions between stoats, mice and seed availability. @blackwell2001 made four predictions regarding the likely effects of stoat predation on mouse dynamics (see Figure 1) with a subsequent field study concluding that stoat predation should have minimal effects on the population dynamics of mice. The authors identified three different phases in the eruption cycle where stoats could have an effect (Figure 1) and the subsequent field study [@blackwell2003] indicated that stoat control has little detectable effect on mouse populations during the peak (Prediction B), decline (Prediction C) and low (Prediction A) phases of the beech eruption cycle.
Subsequent modelling work reached similar conclusions [@tompkins2006; @tompkins2013] but identified that the response of mice to stoat control should depend on interactions with rats. Specifically, [@tompkins2013] concluded that, where rats were present, stoat control alone should allow rats to increase, which would have a suppressive effect on mouse populations through either predation or competition. In contrast, when both stoats and rats were controlled, mouse populations could increase to higher levels than in the absence of control (see Figure 1).
Our aim was to test the predictions outlined in both @blackwell2003 and @tompkins2006 using data from a large-scale field study. Specifically, we measured the abundance of mice and rats on trapping grids over six years in beech forest in two adjacent valleys, one with intensive stoat trapping and one without. In each valley we also manipulated rat densities by including trapping grids where rats were removed and compared these to grids without rat removal. This allowed us to examine if predicted responses of mouse populations to stoat control (see Prediction A-D) was influenced by interactions with rats (Prediction E).
Methods {#method}
Data collection
# 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)
Our study was carried out in two adjacent valleys: the Hollyford Valley (GPS co-ord, S $44^\circ$ $45'$ S $168^\circ$ $10'$ E) which has had intensive stoat control since 1983 to protect a vulnerable population of mohua [@odonnellPredictingIncidenceMohua1996] and the Eglinton Valley (GPS co-ord, S $44^\circ$ $2'$ S $168^\circ$ $5'$ E, 450 $m.a.s.l$) which had no stoat control conducted prior to and during the first part of our study (pers comms).
To test Predictions A-E, we set up 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 Elliot traps (permit numbers?. Each grid comprised $81$ traps with each trap located $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 5 nights. Traps were set prior to 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, rats were double-ear tagged, and each individual was ear notched with a notch position that identified the night of first capture so we could identify animals that had been captured but had lost ear tags. We sexed and weighed each captured animal. We aimed to trap for five nights on each trip but due to unfavourable weather, fewer nights were trapped on some trips.
To measure food availability, we recorded the amount of seed reaching the forest floor by placing four standard circular seed-traps at random locations on each grid. Each trap had a radius of $0.125$ equating to an area of $m = 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 preclude 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.
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 ref ethics?). 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).
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. 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 stopped collecting data from one rat control grid in each valley.
Statistical modelling
Mouse abundance
We used the capture-mark-recapture (CMR) data to estimate the number of mice $(N_{j,k})$ 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 augumentation [@royle2009; @royle2012]. This involved augmenting data for the individual mice captured on each grid during each trip with pseudo-individuals. These psuedo-individuals represent mice 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 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}}$$
Prediction A and B
We tested if stoat control would resulted in more mice during non-mast years when mouse abundance was low (Prediction A) by we comparing the average mouse abundance ($N_{j,t}$) produced by the CMR model above. We filtered the full model output into years when all grids have low seed availability (Figure 1). We used these non-mast years to compare the average differences in areas with and without stoat control using a two-way analysis of variance (ANOVA).
We accounted for heterogeneity between valleys (Eglinton and Hollyford Valleys) and our rat removal conditions where, on some grids, rats were artificially reduced). We assessed normality by visual examination of histograms, quantile–quantile (Q–Q) plots. The were violated for mouse abundance was log transformed ($log_{base}$) to account for the general parametric assumptions of this simple test. 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 food availability [@king1983]. We therefore compared the number of mice per seed (${N_{j,t+1}}/{S_{j,t}}$) among grids with and without stoat control using the same ANOVA test as Prediction A ($log(\frac{N_{j,t+1}}{S_{j,t}})$ $\sim$ $valley + control + rat + rat:control$).
Prediction C and D
To test Prediction C and D we modelled the rate of increase of mice($r_{j,t}$) during the relevant seasons (Figure 1) with rate of increase calculated as:
$$r_{j,t} = log\left (\frac {N_{j,t + 1}}{N_{j,t}} \right )$$
To estimate ($r_{j,t}$) and propagate the uncertainty in our abundance estimates we used the complete data likelihood approach described by @schofield2014 (for details see supplementary material). To test predictions C and D we compared rates of increase with and without stoat control. To do this, we needed to account for other factors known to affect rates of increase in mouse populations, including seed availability, density dependence, and potential interactions with rats.
Seed availability
With respect to seed availability, our study aimed to determine the relationship between the rate of increase of mice populations (a numerical response), however, within the population, each mouse can only consume a maximum amount of seed at any given time. To incorporate this "individual consumption rate ($S_{j,t}$)" we fitted three previously proposed functions to our field data and selected the best fitting model. The three functions were: 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 this model were collected in the Eglinton Valley [@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. We used the transformed variable from the best fitting model above to incorporate individual "seed consumption" $S_{j,t}$.
Numerical response
After accounting for the relationship between seed fall and mouse abundance using an intake rate ($S_{j,t}$), there are still other demographic and environmental processes that are known to regulated mouse populations in New Zealand beech forests. Population abundance has been identified as a key driver of future mouse abundance in multiple studies under similar conditions to this research ($N_{j,t}$; @choquenot2000; @holland2015; @ruscoe2005). It is still not clear what these density dependent processes are but it is still regarded that they are an important component of mice dynamics [@holland2015]. In our model we included the estimate of mouse abundance at the previous time step as a measure of density dependence ($N_{j,t-1}$).
Rats may also influence mouse populations via competition, because rats and mice have overlapping diets [@king2005; @mcqueen2008], or possibly through direct predation [@bridgman2013]. We were interested in whether there was a difference between areas with and without rats in relation to the outcome of mice dynamics. To do this we using the minimum number of rats in the season before ($R_{j,t-1}$).
We identified the most appropriate numerical response was:
$$r_{j,t} = \beta_{0,s,v,c} + \beta_{1,s,v,c} (S_{j,t})+ \beta_{2,s,v,c} (N_{j,t-1}) + \beta_{3,s,v,c} (R_{j,t-1})$$
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 mesopredator release only occurs when rat populations are reduced to low densities. (Prediction E) by comparing differences between population model coefficient for rats ($\beta_3$) grouped into the grids where we reduced rat densities and the others that where not. 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 showing differing population dynamics.
We captured rats (ship rats and Kiore combined) less often on grids than mice. Rather than a full capture-recapture analysis, we used the number of individual rats captured on each grid during each trip (the minimum number present) as an index of rat abundance. For each trip and grid we estimated seed availability as the number of seeds per $m^2$, averaged over all seed traps on a grid.
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].
Results {#results}
We captured at total of 2370 individual mice, 219 ship rats, and Kiore during a total of 94 trips. The stoat control operation in the Eglinton Valley removed a total of 792 individual stoats between January 1999 and June 2001 [@ruscoe2003]. During this operation, low numbers of rats were also removed as by-catch in stoat traps but no other independent predator control was conducted in either study Valley during our data collection. Incidental stoat captures during our routine rodent trapping sessions were uncommon (proportion of total captures that were stoats < 0.01). Lower numbers of stoats were captured in the Eglinton Valley ($N_{stoats} = 8$) where the long-term predator control program was undertaken than the Hollyford Valley ($N_{stoats} = 20$).
Seedfall
Our study captured two years of high seedfall ($2000$ and $2003$). During all other years except the first year ($2001$, $2002$ and $2004$) we observed low seedfall on all grids (Figure 2). As expected, we observed high variability both between and within each trip. For example, in $2000$, the difference in seedfall between grids varied from the highest largest estimate of available seed ($Seed_{7,2} = 3387 m^2$; Eglinton Valley) and two grids recording zero seedfall for the same period ($Seed_{4,2}$ and $Seed_{4,5}= 0$ $m^2$; Hollyford Valley) on the same trip.
# data
source("./R/r-packages-needed.R")
source("./R/wrangling/Data_CRinput_mice_jan2019.R", echo = FALSE)
source("./R/Davidson_2019_Data_wrangling2.R")
source("./R/figures/seed_black_white.R", echo = FALSE)
# fig.2.plot.design
plot.final.f3.1
cap.seed <- c("The average avaliable beech seed collected each season during the study period (Autumn 1999 to Winter 2004). The shape distinguishes the two valleys apart (triangle = Eglinton Valley, circle = Hollyford Valley). Hollow symobols represent cases where stoats are controlled and solid fill in valleys with stoat were uncontrolled.")
#blown up
knitr::include_graphics(c("./figs/seeds.png"))
Mice
We estimated mouse abundance ($N_{j,t}$), where $j$ represents each unique grid and $t$ for each trip from the overall community dynamics model. We found that mouse abundance was greatest in the Eglinton Valley in grids and trips where seed availability was also high. In any given year mouse abundance was on average greater in the Eglinton Valley (Figure 3; triangles) than the Hollyford Valley (Figure 3; circles).
# # Database.plot
# data
source("./R/wrangling/Data_CRinput_mice_jan2019.R", echo = FALSE)
source("./R/Davidson_2019_Data_wrangling2.R", echo = FALSE)
# # glimpse(meanM)
# source("./R/figures/mouse-plot-may2019.R", echo = FALSE)
fig.3.N
cap.mice <- "The abundance of mice ($N_{j,t}$; point estimates) collected each trip during the study period (Autumn 1999 to Winter 2004). The shape distinguishes the two valleys apart (triangle = Eglinton Valley, circle = Hollyford Valley) and ) and solid symbols represent cases where stoats are uncontrolled and hollow in valleys where stoat control was undertaken."
#blown up
knitr::include_graphics(c("./figs/mice.png"))
Rats
Rats and mice displayed a similar response to beech seed (Figure 4 and 5 respectively). However, the overall number of rats captured ($R_{j,t}$) remained relatively low in both valleys throughout our six-year study, except for a single trip in Spring 2004 (n = 25). The highest rat abundance was recorded where stoat control was undertaken (Hollyford Valley).
# plot rat CIs
rat.plot <-
ggplot(p.design2, aes(y = N,col = Rats,shape = Valley,fill = Control,x = Date)) +
geom_rect(aes(xmin=ymd('2000-01-01'),xmax = ymd('2000-12-31'),ymin = -Inf,ymax = Inf), colour = "grey90", fill = "grey90") +
geom_rect(aes(xmin=ymd('2002-01-01'),xmax = ymd('2002-12-31'),ymin = -Inf,ymax = Inf), colour = "grey90", fill = "grey90") +
geom_rect(aes(xmin=ymd('2004-01-01'),xmax = ymd('2004-12-31'),ymin = -Inf,ymax = Inf), colour = "grey90", fill = "grey90") +
geom_errorbar(aes(ymin = lcl.slog, ymax = ucl.slog, width = 0), lwd = 0.7, col = "black", position = location.move) +
geom_line(data = rat.mean,
aes(y = mean.rat,
x = Date),
size = 0.95,
col = "grey50", position = location.move) +
geom_point(data = rat.mean,aes( y = mean.rat, col = Rats,
shape = Valley,
fill = Control,
x = Date
),size = 6 , alpha = 0.7, position = location.move) +
xlab(expression(paste("Time", "(", italic(t), ")"))) +
ylab(expression(atop(paste("Minimum "," ", " number"," "),
paste(" ", "of"," ",
"rats"," ","(",italic(R[jt]),")"))) ) +
scale_color_manual(name = "Stoat Control",
values = c("black", "white", "white")) +
scale_shape_manual(name = "Ecosystem",
values = c(24, 21)) +
scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
# manually define the fill colours
scale_fill_manual(name = "Stoat Control",
values = c("white", "black", "white")) +
xlab(expression(paste("Time", "(", italic(t), ")"))) +
ylab(expression(atop(paste("Minimum "," ", " number"," "),
paste(" ", "of"," ",
"rats"," ","(",italic(R[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
# Remove fill legend and replace the fill legend using the newly created size
guides(
col = "none",
size = guide_legend(override.aes = list(
shape = c(15,0),alpha = 1
)),
shape = guide_legend(override.aes = list(
shape = c(24, 21), size = 4
)),
fill = guide_legend(override.aes = list(
col = c("white", "black"),shape = c("square"),
size = 4
))
) +
theme_tufte() +
theme_bw() +
theme(strip.background = element_blank(),
strip.text.y = element_blank(),
plot.title = element_text(hjust = 0, size=24, family = "Times", color="black", margin = margin(t = 10, b = 10)),
plot.subtitle=element_text(size=16, face="italic", color="black"),
legend.position = "none",
legend.key = element_blank(),
legend.background = element_rect(fill="white", size=1),
legend.key.size=unit(1,"cm"),
legend.text = element_text(colour = "black", size =16, family = "Times"),
legend.title = element_text(colour = "black", size =16, family = "Times"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.spacing = unit(2, "lines"),
panel.border = element_blank(),
axis.title.y = element_text(colour = "black",size =20, family = "Times", angle = 90),
axis.title.x = element_text(colour = "black", size =20, family = "Times"),
axis.text.y=element_text(colour = "black",size = 20, family = "Times"),
axis.text.x = element_text(colour = "black", size =20, family = "Times"),
axis.ticks.x = element_line(size = 1),
axis.ticks.y = element_line(size = 1),
axis.line.x = element_line(size = 1),
axis.line.y = element_line(size = 1),
strip.text = element_text(face="bold",colour = "black",size =14, family = "Times"))
rat.plot
# plot
#rat.plot
#caption
# , fig.cap=cap.rat
cap.rat <- c("The avarage number of rat captures ($R_{j,t}$; point estimates) collected each trip during the study period (Autumn 1999 to Winter 2004). The shape of each point distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley) and solid symbols represent cases where stoats are uncontrolled and hollow in valleys with stoat removal.")
#blown up
knitr::include_graphics(c("./figs/rats.png"))
All other trips with high rat records were in areas with stoat control. In the Hollyford Valley during both mast-years, rats increased irrespective of stoat control. No rats were captured after May 2002 in the Eglinton Valley.
Numerical response
We found that the intake rate we transformed from the seedfall data $(S_{j,t})$ always had the greatest impact on the rate of increase of mice $(r_{j,k})$. Greater than both mouse density ($N_{j,t}$) and rat presence ($R_{j,t}$; see supplementary material). Our model accounted for the confounding effects of numerical processes such as density ($N_{j,t}$) and food availability ($S_{j,t}$) using our Bayesian hierarchical model and then tested for evidence of mesopredator release of mice during our study. We used the fitted model parameters to compare the differences between the key predictions proposed for areas with and without control after accounting for both the observation error (data collection and population estimates) and process error (population dynamics; @ahrestani2013).
# s.final.model$coefficients[2,4]
# s.final.model$coefficients[3,4]
Prediction A
Stoat control does not increase mouse populations at times of low abundance.
We compared the abundance of mice during trips when both mouse abundance and food availability were low ($2002$ and $2003$ using the model. We found no visually observable or statistically significant differences in the estimates between areas with and without stoat control ($1.00$ SE = $1.99$, p-value = $0.62$) during any of the trips during our study.
# summary tables
# if rendering in word
# , ft.split = TRUE, ft.align = "left"
# summary table
# flextable::flextable(data.frame(co.effs, summary(low.mod3)$coefficients))
# # # data
# source("./R/figures/pa-plot-code.R", echo = FALSE)
# #
# # # Plot
# low.plot.time
# +
# theme(legend.position = "right")
# high.plot.time
# gridExtra::grid.arrange(low.plot.time, output.dens, ncol = 2)
# Figure Caption
cap.Pa <- c("The unmodified abundance of mice in non-mast years (2001, 2003). The shape difference distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley). The grids are divided into the two groups on the x-axis (areas with and without stoat control). Solid symbols represent cases where rats are present and removed at points with hollow symbols.")
# knitr::include_graphics(c("./figs/fig41.png"))
# Figure Caption
cap.Pa <- c("The unmodified abundance of mice in non-mast years (2001, 2003). The shape difference distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley). The grids are divided into the two groups on the x-axis (areas with and without stoat control). Solid symbols represent cases where rats are present and removed at points with hollow symbols.")
source("./R/data_input.R")
source("./R/figures/pa-plot-code.R")
low.plot.out
For Prediction A we did not need to account for varying seed input because seed fall was biological equivalent to zero for all grids and trips during low abundance seasons [@choquenot2000].
Prediction B
Stoat control does not result in higher peak abundance of mice
# run source script
# source("./R/figures/pb-plot-code.R", echo = FALSE)
# summary
# s.final.model <- summary(high.mod3)
# # s.final.model$coefficients[2,4]
#
# # output from models in workable format
# co.effs.m2 <- c(row.names(as.data.frame(summary(high.mod2)$coefficients)))
# co.effs.m3 <- c(row.names(as.data.frame(summary(high.mod3)$coefficients)))
#
# # results data
# pb.results.m2 <- data.frame(co.effs.m2, summary(high.mod2)$coefficients)
# pb.results.m3 <- data.frame(co.effs.m3, summary(high.mod3)$coefficients)
Our simple ANOVA model (see here for analysis) accounted for the differences in sample size between the areas with and without stoat control and un-equal variances within groups and found no significant difference between the stoat controlled and uncontrolled areas (What is best here and I can repeat for all?).
# simple anova comparison
#flextable::flextable(pb.results.m3)
During the years of peak abundance in our study we did observe lower mouse abundance in the Hollyford Valley compared to the Eglinton Valley ($-0.17$, SE = $0.09$, p-value = $ 0.06 $). Although this is not below the 0.05 threshold. It should be noted that when the rat removal treatment was removed the p-value did drop to 0.05.
# #data
# source("./R/figures/pb-plot-code.R", echo = FALSE)
#plot
# high.plot.time
# +
# theme(legend.position = "right")
# Figure Caption
# low.plot.time
# gridExtra::grid.arrange(high.plot.time, out.dens.high, ncol = 2)
cap.Pb <- c("Difference between areas with and without stoat control at peak mouse abundance. The shape difference distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley). The grids are divided into the two groups on the x-axis (areas with and without stoat control). Solid symbols represent cases where rats are present and removed at points with hollow symbols.")
source("./R/figures/pa-topb-test.R")
low.plot.out
Prediction C
Stoat control did not increase mouse populations as food becomes available.
We tested the differences in the rate of increase of mouse populations between the four seasons, two valleys and stoat control areas (excluding the rat treatment) during seasons when mouse population were increasing (August; Autumn - Winter). We found was no significant differences between mice populations ($ r_{j,t} $) during the increasing seasons of mouse dynamics. All 95% credible intervals for the differences between these trips included zero (Figure 10).
#rate vs seed data with lines for average prediction
source("./R/pc-plot-code.R", echo = FALSE)
# # just most interesting variable
# # caption
cap.Pc = c("Visual representation of the effect of stoat control on the rate of increase of mouse populations during the increasing seedfall (Autumn to Winter; **Prediction C**). The different shapes distinguish the two valley systems (triangles = Eglinton Valley; circles = Hollyford Valley) and all solid symbols represent cases where rats are present and removed at points with hollow symbols.")
pc.plot
# knitr::include_graphics(c("./figs/seeded.png"))
# Save plot
# export plot for example vignette
# png("C://Code/beech-publication-wr/figs/seeded.png")
# pc.plot
# dev.off()
Prediction D
Stoat control did not hasten the decline from peak abundance.
# knitr::purl("Davidson_2019_Simple_models.Rmd")
# source("./R/Davidson_2019_Model_wrangling.R")
# source("./R/figures/pd-plot-code.R")
# pd.plot
This prediction specifically focused on the Spring and Summer seasons, when rodent populations are crashing (Figure 3 vs. Figure 4). Overall density had a weak negative effect on mice dynamics under all conditions (Figure 3) but was greatest during the declining seasons of mice dynamics. After accounting for seed and density in the population model there was a faster rate of decline in the Hollyford Valley (large solid circles) compared to the Eglinton Valley (large hollow circles) but this was not statistically significant.
#data for plot
# source("./R/figures/rate-density-plot.R", echo = FALSE)
#plot
# pd.plot.feb +
# theme(legend.position = "right")
#caption
cap.Pd = c("Visual representation of the effect of stoat control during decreasing mouse abundance (Spring to Summer; *Prediction D*). The different shapes distinguish the two valley systems (triangles = Eglinton Valley; circles = Hollyford Valley) and all solid symbols represent cases where rats are present and removed at points with hollow symbols.")
# source("./R/data_input.R")
source("./R/pd-plot-code.R", echo = FALSE)
pD.plot
# export plot for example vignette
# png("C://Code/beech-publication-wr/figs/mouse_dens.png")
# pD.plot
# dev.off()
# knitr::include_graphics(c("./figs/mouse_dens.png"))
Prediction E
Does the presence of rats impact the population dynamics of mice at each of the seasons tested in Predictions A-D.
# knitr::purl("Davidson_2019_Simple_models.Rmd")
# source("./R/Davidson_2019_Model_wrangling.R")
# source("./R/figures/pe-source-code.R")
# pe.plot
We found that rat numbers $R_{j,t}$ had the smallest overall effect on mice dynamics ($ \mu(\beta_3 $ $|$ $All$ $Seasons$$ )=$ $0.007$) and was the most variable parameter in our community model ($ CV_{rats} = 5.74 $). We recorded the highest number of rats in February and May, at the same times of high mouse abundance (Figure 4).
#plot code
source("./R/pe-plot-code.R", echo = FALSE)
#plot
# rat.plot.all
# theme(legend.position = "right")
#caption
pE.cap <- c("Differences in estimates between rats during each season of mice dynamics.")
# Save plot
# export plot for example vignette
# png("C://Code/beech-publication-wr/figs/fig_rat.png")
# pe.plot
# dev.off()
# knitr::include_graphics(c("./figs/fig_rat.png"))
The relationship between rats and mice was examined by comparing the estimated differences between treatments (Figure 8; lines represent mean relationships in the three treatment groups). We also statistically testing the difference between each different treatment and found no differences across all parameters (Figure 9).
`r
pe.plot
We did not have enough data to estimate a statistical interaction between rats and mice (non-linear relationship) or other more complex non-linear models for rat dynamics.
Discussion {#discussion}
Overall we found no evidence of mesopredator release in the mouse populations we monitored. We used a bayesian mixed model (bmm; Appendix 1) to model mice dynamics in New Zealand Beech forests. We tested five predictions in a relation to stoat control and mouse abundance. We experimentally manipulating predator control in two independent beech forests and compared the differences between the rate of increase and abundance of mice. The model replicated mouse dynamics in our data and produced similar results to other NZ mammal models (@choquenot2000; @ruscoe2001; @ruscoe2005; @tompkins2006; @tompkins2013; @holland2015; summarised in Figure 1). Generally in NZ Beech forests, high food availability results in high mouse abundance in the subsequent months until natural germation of the seed and is commonly observed in both our data and a range of other NZ forest types [@innes; @sweetapple; @ruscoe2011; @ruscoe2012]. Modelling has incorporated these parameter estimates to represent invasive species dynamics in NZ forests [@tompkins2006; @tompkins2013; @holland2015] A naïve investigation of our summarised data (Figure 3.1, 3.2, and 3.3) might have falsely lead to similarly low powered conclusions, that rats may have an effect on mice dynamics as well as many generally accepted relationships. By using a BHM framework we were able to test the theoretical relationships proposed using experimental data.
Seed fall
The intake rate had the greatest impact on mice populations in all seasons. We found that $r_{j,t}$ was strongest on grid and trips when beech seed was increasing. After comparing a selection of functional responses to our data (Appendix 1.3) we found the best type II response and converted our estimated $Seed_{j,k}$ per $m^2$ to the "intake rate" $S_{j,t}$. We found that the effect of mouse density ($N_{j,t}$) and rats ($R_{j,t}$) was lower than the "intake rate" ($S_{j,t}$). We observed that the effects of density were greatest during the declining seasons. These differences in effect sizes are similar to other studies.
Mouse density
The effect of mouse density and rats on was much lower and were also similar to other published results [@choquenot2000]. We and many other studies observed high heterogeneity in abundance estimates. This type of data can be difficult to model correctly in many frequentist frameworks. To account for this heterogeneity in the biological data we used experimental manipulation and Bayesian modelling framework to account for the confounding effects of food availability, mouse density, and rats to simultaneously assess the overall effect of stoat control on mice dynamics. Importantly, our model captured both the observational and process variability within NZ beech forests. We also conclude that the underlying heterogeneity in small mammal populations was sufficiently large to make estimating true response of mice to rats difficult in biological experiments (Figures 4-9).
Rat presence
Our results suggest that it is unlikely that rats have strong impacts on the rate of increase of mice in our large scale study system. Contrary to this, previous laboratory experiments suggested that interactions between mice and rats were negatively related to predation or competition effects and limited evidence of competition release has been observed in mixed podocarp-tawa forests in NZ. We found that rats had the smallest impact on mice dynamics (Figure 8). In part, this is due to the high variability in the number of rats and kiore between grids and sampling trips in our study (??), which is often the case in non-laboratory studies. Instead of rats/kiore being present but at low numbers in all grids (low effect size and low uncertainty), there were only a few grids that had comparability high numbers. This result suggests that refugee areas may exist and can support larger rodents in beech forests that current averages suggest. Overall, it is likely that beech forests do not have enough resources to support larger spatially uniform populations of rats.
Bayesian modelling
The application of BHM's in ecology has been increasing over the past decade [@king2012]. BHM modelling can provide challenges when selecting models and assessing model fit [@auger-methe2016]. To account for this, we used the general associations above as well as AIC and r-squared. We then used these parameters to fit the proposed theoretical population model. We investigated the patterns in our model to insure we adequately encapsulated the patterns in the data with our model to assess the effects of our parametrisation above. Large-scale ecological experiments are very informative but often need to be interpreted with caution due to often unique problems. We address these below and add additional direction to address these uncertainties in the results we found.
Our population model representation is relatively simplistic compared to some of the previous models that incorporated more interacting species and subsequently complex interactions. Our experimental design was just a big taking advantage of the situation to built a simple population model from CR data. A benefit to a simple BHM model is the ability to use bayes theorem to re-assess additional data in context of the proposed predictions. In the context of the data we have used in this experiment and uninformative priors. We choose to focus on high quality data collection methods and the key dynamics already identified with invasive species dynamics in NZ beech forests. There are however several key aspects that may be confounding our results, to the observational and population model choices.
Stoat control
An alternative explanation for the lack of any observed effects of stoat control on mouse abundance was that the removal of animals was potentially limited by the effectiveness of a particular stoat control program, management or experimental design.
If trapping did not remove a sufficient number of individual stoats to reduce predation pressure we will not observe the theoretical responses proposed by @tompkins2013 and others.
High individual heterogeneity in the capture rates of stoats at both high and low mouse densities makes reductions to zero individuals virtually impossible targets. This was the case in our study where by-catch of stoats in controlled areas was observed. As an indication of removal, $4.75$ times more stoats were removed from the stoat control program (Hollyford Valley 2001-2004) than caught in by-catch in areas already being controlled prior to the experiment. Our before/after treatment tested this and verified that we did not find differences within the grid when experimental treatments were changed. We may also be missing the impact of by-catch of other pests in stoat traps. We believe that any effect this may have on mice dynamics would be limited due to the limited impacts of rat removal in our data. Biologically, the differences in home range size of rodent species (10's ha; [@bramley2014; @innes1983; @pryde2005]) relative to the home range of stoats and the scale of trap line (100's ha; [@miller2001; @murphy1995]) which suggests that population level impacts from these removals is unlikely due to the limited number of animals removed and the re-invasion biology of rodents [@bramley2014].
<- is this too critical given ?... ->
[@blackwell2001; @blackwell2003]: However, the limited geographical scale of the study, low statistical power and replication between different forest types makes it difficult to compare the effect of different treatment types during the four different seasons.
Repleccate and extend on larger scale.
<- is this too critical given ?... ->
Resource pulse systems
During our study we found stoat control had minimal impact on mice populations compared to beech seed. Beech forests have previously defined as resource pulse systems [@wardle]. Strong resource pulses cause strong effects throughout an ecosystem [@yang2008]. The process behind for this phenomenon is that even in the absence of control, stoats cannot attain high enough densities to limit mouse populations under conditions of high food availability [@king1983]. Furthermore, the intrinsic differences in the reproductive and growth rates (life history traits) of each species reduces the ability of stoats to regulate mouse populations at the beginning of a masting event[@king2011]. By taking the proposed processes that effect mice dynamics and applying a hierarchical model framework allows us to test and estimate effect sizes that further reduce uncertainty in the outcomes of predator control in NZ beech forests.
The effects of climate change and other human lead impacts will have lasting impacts on on future rodent dynamics [@holland2015]. Population models of climate change indicate a shift to prolonged but lower overall mast events where it may be possible to observe mesopredator release of mice with the removal of predators [@tompkins2013]. Mice populations monitored in more complex forest ecosystems suggest complex interactions between top-down and bottom-up effects [@ruscoe2011; @ruscoe2012]. These studies suggest that rodent populations would increase in other NZ forest ecosystems after the removal of stoats [@ruscoe2011] or other invasive predators [@rayner2007].
Challanges
Population level responses such as abundance and rate of increase are often measured using indices (e.g. MNA). This can lead to difficulty in estimating population level parameters such as density. We aimed to avoid this bias as we estimated mice abundance using capture-recapture (CR) data and an integrated population model to correctly account for uncertainty in abundance estimation when fitting the population model of mice dynamics. However due to limited population numbers and data, only indices could be used to represent rat interactions. Independent research in beech forests has shown a high correlation between indices of rats in this valley and CR data. Nevertheless, increasing the quality of the rat/kiore data would most likely reduce uncertainty on the estimates of these parameter estimates. By calculating the maximum increase from the fitted population model we were able to compare the maximum population growth rate to other studies ($R.max_{mice} = actual number here$; @hone2010), however, analogous models are not directly comparable because they do not include a resource component (e.g. beech seed production) or use indices (e.g. $R_{j,t}$). Field studies are also difficult to directly compare because may studies incorporate structurally different ecosystems (e.g @ruscoe2011), species indices instead of true abundance estimation using CR models or did not report effect sizes. Although this limitations exist our comparisons to general trends were all comparable to our model estimates and our estimates are all biologically viable.
Management implications
Our model has subtle but important differences compared previous research when in comes to management application. We have incorporated a reproducible workflow[@britishecologicalsociety2018] for the future development of model testing for different and new datasets [@wickham2014]. In doing this we have incorporated leading reproducible science techniques to help address any aspects of our study may need for further support. Often studies of such large scale can address key population-level questions however these studies struggle to find replication because of the scale at which they are conducted [@oksanen2001]. Many issues that many confront a PFNZ2050 will be reduced by using simple but well fitted BHMs for predicting and allocation management resources.
Predator-free NZ
Our integrated population modelling approach and the large, high-quality CR data used to parametrise the model allowed us to predict and report the biologically relevant effect sizes. We used this to address the direct question of stoat control in NZ Beech Forests. Future studies can now test and assess the beech forest ecosystem compared to the data they collect. By using the unified modelling framework presented in this paper (Appendix 3; [@davidson2019), reproducible research practises and "tidy-data" [@wickham2014]. Many of the processes we have verified were previously hypothesized using indices of mice abundance [@blackwell2001; @blackwell2003], laboratory experiments [@bridgman2013], novel systems such as remote islands [@mulder2009] and small patches of mainland forest [@blackwell2001; @blackwell2003]. We have integrated these results into a framework that can now test these results in other systems that are hypothesized to be different. It will be important for future work to continue monitoring this system for anomalies and unexpected patterns such as increases in other invasive predators, particularly rats.
`r if (knitr:::is_html_output()) '
References {-}
'`
Sorting references has been much hard than I expected when there are multiple types of software to use and different ways to link them.
The core problems I have needed to solve are:
-
Storing pdfs locally and working on them.
-
Mendeley works well for looking and searching through pdfs and was the reason I switched from EndNote (that my supervisors were using at the time) to Mendeley.
-
Convert references attached to pdfs in reference manager into a bibtex file so that citr and refmanageR can read these.
-
Zotero because it has good extentions that can turn bib references into bibtex structure for R packages that use a variation on the .bib structure in R.
-
Convert stable referencing in .Rmd to other document types and structures.
-
This has now become hard to understand how this is all working so I have had to layout a software path (and maybe future Docker image) to ensure that the refences converge and the errors are solvable (or at least I can work out where they come from)
Editor??
To work with .bib files (to find errors to begin with) it is much nice to have some code snippit support (highlighting syntax)
davan690/beech-publication-wr documentation built on March 29, 2020, 11:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Summary {#summary}
[still needs work but last]
The conflicting outcomes from previous studies provides the motivation for this more comprehensive field study assessing the effects of stoat control on mouse dynamics. Here we clarify the discrepancy between research suggesting that mesopredator release of rodents will occur in New Zealand forest systems and field studies that have presented limited but conflicting support for increases in mouse abundance following predator control. We used an experimental design to test for differences between areas with and without stoat control after accounting for known effects of resources, density and competition on mouse dynamics in NZ beech forests. We found no evidence to suggest mice will become more abundant after predator removal.
# libraries needed source("./R/r-packages-needed.R", echo = FALSE) require("dplyr") require("knitr") require("docxtools") require("citr") #document global rules knitr::opts_chunk$set(comment=NA, # fig.path = "../figs/", echo=FALSE, fig.height=5, fig.width=7, message=FALSE, warning=FALSE, tidy = TRUE # , cache.extra = packageVersion('tufte') ) # options(knitr.table.format = 'markdown') # options(htmltools.dir.version = FALSE) # how do I do this?? # ,eval = FALSE,include = FALSE # export .r code only # knitr::purl("./Davidson_2019_BeechForest.Rmd") # render draft to webpage # rmarkdown::render(input = "Davidson_2019_BeechForest.Rmd") # , # output_format = "html_document", # output_file = "Davidson_2019_t.html") ## ----figure themes------------ source("./R/theme_raw_fig3s.r", echo = FALSE) source("./R/davidson_2019_theme.r", echo = FALSE)
## ----sorting-references-------------------------------------------------- # citr::md_cite("./Beech-forest.bib") # bib(file = "Beech-forests.bib") # automatically create a bib database for R packages knitr::write_bib(c(.packages(), 'bookdown', 'knitr', 'rmarkdown','tidyverse', "compareGroups", "jagsUI" ), 'packages.bib')
#data and plotting import # source("./R/r-packages-needed.R") # source("./R/theme_raw_fig3s.r") # source("./R/davidson_2019_theme.r") # # source("./R/Davidson_2019_Simple_models.R") # source("./R/Davidson_2019_Model_wrangling.R") # source("./R/Davidson_2019_Data_wrangling.R") # source("./R/figures/pb-source-code.R") # source("./R/figures/pa-plot-code.R") # source("./R/figures/pb-plot-code.R") # # source("./R/data_input.R") # source("./R/figures/pd-plot-code.R") # # source("./R/figures/pd-plot-code.R") # source("./R/figures/rate-density-plot.R") # source("./R/figures/pa-plot-code.R")
bibliography: ["Beech-forests.bib", "references.bib"]
Introduction {#intro}
Worldwide, but particularly on islands, introduced predators have significant impacts on native species [@gurevitch2004], and in the worst cases extinctions [@doherty2016; @towns2011]. For example, on the Island of Guam, the invasion of the Brown Tree Snake (Boiga irregularis) was the primary driver of multiple extinction events [@fritts1998]. In New Zealand, mammalian predators contribute to a disprorpoitatly large group of invaders impacting native populations nationally [cite]. However, removing, or reducing the abundance of the top predator can lead to increasing numbers of mesopredators (for a review see @prugh2009 and for examples see @rayner2009; @robles2002), which in turn, can lead to unintended and suprising outcomes [cite suprise island papers]. For example foxes predate on native and evolutionaryly distinct marsupials throughout Australia [cite margaritas paper??], heavyily impacting the full range of many native species [@kinnear1998].
Introduced species also have destructive impacts in New Zealand’s (NZ) remaining/remnant native forests. In these systems direct predation is commonly proposed as the reason for native species decline [@fitts1998; @kinnear1998]. These contemporary ecological communities commonly contain of at least one of the following four introduced mammalian predators;
- Stoats (Mustela erminea; Veale (2015));
- Brushtail possums (Trichosurus vulpecula Kerr; Clout and Gaze (1984));
- Ship rats (Rattus rattus; Innes (2005) referred to as rats);
- House mice (Mus mus musculus; @allen2006 referred to as mice).
Stoats are the generally regarded as the top predator when they are present in New Zealand ecosystems [@king1983], following their deliberate introduction in the late nineteenth century [@king2017]. Stoat control is commonly undertaken in these systems to protect native birds that are vulnerable to mammalian predation [@white2006], in particular, hole-nesting species like mohua (Mohoua ochrocephala; @Odonnell2006]. Stoat control is now commonly undertaken to protect native birds that are vulnerable to predation [@white2006], in particular, hole-nesting species like mohua (Mohoua ochrocephala; [@odonnell2006]). However, the primary food source for stoats in NZ forests are rats and mice [@king2006], and consequently there is a concern that reducing stoat populations to protect native birds may allow rodent populations to increase (e.g @rayner2009). An increase in the number of rats and mice could offset the benefits of stoat control because rodents are known to consume the eggs and chicks of native birds [@russell2016; @towns2006; @latham2018], directly compete with native species for food resources such as flowers and seeds [@mcqueen2008] and predate on invertebrates [@ruscoe2012]. In this paper we address the question: does stoat control lead to increased abundance of rodents, particularly mice, in NZ beech forests?
####UP TO HERE>>>>
We carried out our study in beech forests in the South Island of New Zealand where stoats are the main predator and mice are the most abundant rodent species [@king1983]. Ship rats (Rattus rattus; @innes2005) and kiore (Rattus exulans; @roberts1991) are also present but at much lower numbers [@jones2011]. We tested whether reducing stoat populations by predator control influenced mouse population dynamics as shown in Figure 1, 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 1: Prediction A to D). If stoat control affects mouse population dynamics then, relative to sites with stoat control (Figure 1; hollow symbols), mouse populations at sites without stoat control should exhibit changes in abundance in the direction of red arrow in Figure 1:
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.
Studies elsewhere in the world have shown that removing or reducing the abundance of a top predator often leads to an increase in the numbers of predators at lower trophic levels (termed mesopredator release), which in turn, can lead to unintended and often negative outcomes for native species (for a review see @prugh2009 and for examples see @rayner2007; @robles2002). While mesopredator release has been widely documented elsewhere, it is unclear if stoat control in NZ forests will cause rodent populations to increase. Rodent populations in NZ forests are known to respond strongly to variation in food supply [@choquenot2000; @ruscoe2001; @blackwell2001; @blackwell2003; @ruscoe2005; @tompkins2006; @tompkins2013; @holland2015; @latham2017], primarily seed availability (Figure 1).
# #data file # source("./R/ecosystem-simulation/sim-raw-data.R", echo = FALSE) #file code source("./R/figure_1_script.R")
result.plot cap1.stes <- c("*Expected changes in mouse populations over time in New Zealand forests (bottom panel) in response to changes in seed availability (top panel) during non-mast and mast years (grey blocks). In the bottom panel, differences between the hollow and filled symbols represent the proposed differences between areas with and without stoat control respectively. The arrows labelled A-D represent the four predicted outcomes of stoat removal that we tested; A) during non-mast years when little seed is available we expect that mice populations will be similar, B) at the peak of mouse abundance (during winter and spring in mast years) when we would expect larger mouse populations in areas where stoats are no longer predating on mice, C) during increased seed availability in areas with stoat control would expect to be able to support larger mouse populations; D) when mouse populations are declining from peak abundance, stoats would additionally reduce mouse abundance through predation in addition to other density dependent processes.*")
This is particularly pronounced in beech forests (spp. nothofagus) throughout New Zealand, where, between years beech seed production is highly variable. Little seed is produced in most years (Figure 1a: non-mast years) with occasional years of high seed production (Figure 1a: mast years; grey boxes). Mouse populations are low in the non-mast years, due to low food availability [@choquenot2000; @king1983]. In mast years, when seed becomes abundant, mouse populations can increase quickly following a predictable seasonal cycle. Seed begins to fall and accumulate on the forest floor in late summer (February) allowing mouse populations to increase, with mouse populations typically remaining high through winter (August) and into the following spring (November). Beech seed that is not consumed by mice and other seed predators germinates in spring to early summer, meaning this food resource disappears and mouse populations begin to decline in the subsequent seasons. If the following year is a non-mast year with little seed available, mouse populations fall to low levels. It is unclear whether stoat populations can increase at a great enough rate to exert sufficiently strong predation pressure to alter these food-driven population eruptions.
Previous studies have investigated the likely response of mouse populations to stoat control by modelling the outcome of interactions between stoats, mice and seed availability. @blackwell2001 made four predictions regarding the likely effects of stoat predation on mouse dynamics (see Figure 1) with a subsequent field study concluding that stoat predation should have minimal effects on the population dynamics of mice. The authors identified three different phases in the eruption cycle where stoats could have an effect (Figure 1) and the subsequent field study [@blackwell2003] indicated that stoat control has little detectable effect on mouse populations during the peak (Prediction B), decline (Prediction C) and low (Prediction A) phases of the beech eruption cycle.
Subsequent modelling work reached similar conclusions [@tompkins2006; @tompkins2013] but identified that the response of mice to stoat control should depend on interactions with rats. Specifically, [@tompkins2013] concluded that, where rats were present, stoat control alone should allow rats to increase, which would have a suppressive effect on mouse populations through either predation or competition. In contrast, when both stoats and rats were controlled, mouse populations could increase to higher levels than in the absence of control (see Figure 1).
Our aim was to test the predictions outlined in both @blackwell2003 and @tompkins2006 using data from a large-scale field study. Specifically, we measured the abundance of mice and rats on trapping grids over six years in beech forest in two adjacent valleys, one with intensive stoat trapping and one without. In each valley we also manipulated rat densities by including trapping grids where rats were removed and compared these to grids without rat removal. This allowed us to examine if predicted responses of mouse populations to stoat control (see Prediction A-D) was influenced by interactions with rats (Prediction E).
Methods {#method}
Data collection
# 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)
Our study was carried out in two adjacent valleys: the Hollyford Valley (GPS co-ord, S $44^\circ$ $45'$ S $168^\circ$ $10'$ E) which has had intensive stoat control since 1983 to protect a vulnerable population of mohua [@odonnellPredictingIncidenceMohua1996] and the Eglinton Valley (GPS co-ord, S $44^\circ$ $2'$ S $168^\circ$ $5'$ E, 450 $m.a.s.l$) which had no stoat control conducted prior to and during the first part of our study (pers comms). To test Predictions A-E, we set up 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 Elliot traps (permit numbers?. Each grid comprised $81$ traps with each trap located $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 5 nights. Traps were set prior to 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, rats were double-ear tagged, and each individual was ear notched with a notch position that identified the night of first capture so we could identify animals that had been captured but had lost ear tags. We sexed and weighed each captured animal. We aimed to trap for five nights on each trip but due to unfavourable weather, fewer nights were trapped on some trips.
To measure food availability, we recorded the amount of seed reaching the forest floor by placing four standard circular seed-traps at random locations on each grid. Each trap had a radius of $0.125$ equating to an area of $m = 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 preclude 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.
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 ref ethics?). 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).
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. 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 stopped collecting data from one rat control grid in each valley.
Statistical modelling
Mouse abundance
We used the capture-mark-recapture (CMR) data to estimate the number of mice $(N_{j,k})$ 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 augumentation [@royle2009; @royle2012]. This involved augmenting data for the individual mice captured on each grid during each trip with pseudo-individuals. These psuedo-individuals represent mice 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 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}}$$
Prediction A and B
We tested if stoat control would resulted in more mice during non-mast years when mouse abundance was low (Prediction A) by we comparing the average mouse abundance ($N_{j,t}$) produced by the CMR model above. We filtered the full model output into years when all grids have low seed availability (Figure 1). We used these non-mast years to compare the average differences in areas with and without stoat control using a two-way analysis of variance (ANOVA).
We accounted for heterogeneity between valleys (Eglinton and Hollyford Valleys) and our rat removal conditions where, on some grids, rats were artificially reduced). We assessed normality by visual examination of histograms, quantile–quantile (Q–Q) plots. The were violated for mouse abundance was log transformed ($log_{base}$) to account for the general parametric assumptions of this simple test. 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 food availability [@king1983]. We therefore compared the number of mice per seed (${N_{j,t+1}}/{S_{j,t}}$) among grids with and without stoat control using the same ANOVA test as Prediction A ($log(\frac{N_{j,t+1}}{S_{j,t}})$ $\sim$ $valley + control + rat + rat:control$).
Prediction C and D
To test Prediction C and D we modelled the rate of increase of mice($r_{j,t}$) during the relevant seasons (Figure 1) with rate of increase calculated as:
$$r_{j,t} = log\left (\frac {N_{j,t + 1}}{N_{j,t}} \right )$$
To estimate ($r_{j,t}$) and propagate the uncertainty in our abundance estimates we used the complete data likelihood approach described by @schofield2014 (for details see supplementary material). To test predictions C and D we compared rates of increase with and without stoat control. To do this, we needed to account for other factors known to affect rates of increase in mouse populations, including seed availability, density dependence, and potential interactions with rats.
Seed availability
With respect to seed availability, our study aimed to determine the relationship between the rate of increase of mice populations (a numerical response), however, within the population, each mouse can only consume a maximum amount of seed at any given time. To incorporate this "individual consumption rate ($S_{j,t}$)" we fitted three previously proposed functions to our field data and selected the best fitting model. The three functions were: 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 this model were collected in the Eglinton Valley [@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. We used the transformed variable from the best fitting model above to incorporate individual "seed consumption" $S_{j,t}$.
Numerical response
After accounting for the relationship between seed fall and mouse abundance using an intake rate ($S_{j,t}$), there are still other demographic and environmental processes that are known to regulated mouse populations in New Zealand beech forests. Population abundance has been identified as a key driver of future mouse abundance in multiple studies under similar conditions to this research ($N_{j,t}$; @choquenot2000; @holland2015; @ruscoe2005). It is still not clear what these density dependent processes are but it is still regarded that they are an important component of mice dynamics [@holland2015]. In our model we included the estimate of mouse abundance at the previous time step as a measure of density dependence ($N_{j,t-1}$).
Rats may also influence mouse populations via competition, because rats and mice have overlapping diets [@king2005; @mcqueen2008], or possibly through direct predation [@bridgman2013]. We were interested in whether there was a difference between areas with and without rats in relation to the outcome of mice dynamics. To do this we using the minimum number of rats in the season before ($R_{j,t-1}$).
We identified the most appropriate numerical response was:
$$r_{j,t} = \beta_{0,s,v,c} + \beta_{1,s,v,c} (S_{j,t})+ \beta_{2,s,v,c} (N_{j,t-1}) + \beta_{3,s,v,c} (R_{j,t-1})$$
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 mesopredator release only occurs when rat populations are reduced to low densities. (Prediction E) by comparing differences between population model coefficient for rats ($\beta_3$) grouped into the grids where we reduced rat densities and the others that where not. 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 showing differing population dynamics.
We captured rats (ship rats and Kiore combined) less often on grids than mice. Rather than a full capture-recapture analysis, we used the number of individual rats captured on each grid during each trip (the minimum number present) as an index of rat abundance. For each trip and grid we estimated seed availability as the number of seeds per $m^2$, averaged over all seed traps on a grid.
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].
Results {#results}
We captured at total of 2370 individual mice, 219 ship rats, and Kiore during a total of 94 trips. The stoat control operation in the Eglinton Valley removed a total of 792 individual stoats between January 1999 and June 2001 [@ruscoe2003]. During this operation, low numbers of rats were also removed as by-catch in stoat traps but no other independent predator control was conducted in either study Valley during our data collection. Incidental stoat captures during our routine rodent trapping sessions were uncommon (proportion of total captures that were stoats < 0.01). Lower numbers of stoats were captured in the Eglinton Valley ($N_{stoats} = 8$) where the long-term predator control program was undertaken than the Hollyford Valley ($N_{stoats} = 20$).
Seedfall
Our study captured two years of high seedfall ($2000$ and $2003$). During all other years except the first year ($2001$, $2002$ and $2004$) we observed low seedfall on all grids (Figure 2). As expected, we observed high variability both between and within each trip. For example, in $2000$, the difference in seedfall between grids varied from the highest largest estimate of available seed ($Seed_{7,2} = 3387 m^2$; Eglinton Valley) and two grids recording zero seedfall for the same period ($Seed_{4,2}$ and $Seed_{4,5}= 0$ $m^2$; Hollyford Valley) on the same trip.
# data source("./R/r-packages-needed.R") source("./R/wrangling/Data_CRinput_mice_jan2019.R", echo = FALSE) source("./R/Davidson_2019_Data_wrangling2.R") source("./R/figures/seed_black_white.R", echo = FALSE) # fig.2.plot.design plot.final.f3.1 cap.seed <- c("The average avaliable beech seed collected each season during the study period (Autumn 1999 to Winter 2004). The shape distinguishes the two valleys apart (triangle = Eglinton Valley, circle = Hollyford Valley). Hollow symobols represent cases where stoats are controlled and solid fill in valleys with stoat were uncontrolled.")
#blown up knitr::include_graphics(c("./figs/seeds.png"))
Mice
We estimated mouse abundance ($N_{j,t}$), where $j$ represents each unique grid and $t$ for each trip from the overall community dynamics model. We found that mouse abundance was greatest in the Eglinton Valley in grids and trips where seed availability was also high. In any given year mouse abundance was on average greater in the Eglinton Valley (Figure 3; triangles) than the Hollyford Valley (Figure 3; circles).
# # Database.plot # data source("./R/wrangling/Data_CRinput_mice_jan2019.R", echo = FALSE) source("./R/Davidson_2019_Data_wrangling2.R", echo = FALSE) # # glimpse(meanM) # source("./R/figures/mouse-plot-may2019.R", echo = FALSE) fig.3.N cap.mice <- "The abundance of mice ($N_{j,t}$; point estimates) collected each trip during the study period (Autumn 1999 to Winter 2004). The shape distinguishes the two valleys apart (triangle = Eglinton Valley, circle = Hollyford Valley) and ) and solid symbols represent cases where stoats are uncontrolled and hollow in valleys where stoat control was undertaken."
#blown up knitr::include_graphics(c("./figs/mice.png"))
Rats
Rats and mice displayed a similar response to beech seed (Figure 4 and 5 respectively). However, the overall number of rats captured ($R_{j,t}$) remained relatively low in both valleys throughout our six-year study, except for a single trip in Spring 2004 (n = 25). The highest rat abundance was recorded where stoat control was undertaken (Hollyford Valley).
# plot rat CIs rat.plot <- ggplot(p.design2, aes(y = N,col = Rats,shape = Valley,fill = Control,x = Date)) + geom_rect(aes(xmin=ymd('2000-01-01'),xmax = ymd('2000-12-31'),ymin = -Inf,ymax = Inf), colour = "grey90", fill = "grey90") + geom_rect(aes(xmin=ymd('2002-01-01'),xmax = ymd('2002-12-31'),ymin = -Inf,ymax = Inf), colour = "grey90", fill = "grey90") + geom_rect(aes(xmin=ymd('2004-01-01'),xmax = ymd('2004-12-31'),ymin = -Inf,ymax = Inf), colour = "grey90", fill = "grey90") + geom_errorbar(aes(ymin = lcl.slog, ymax = ucl.slog, width = 0), lwd = 0.7, col = "black", position = location.move) + geom_line(data = rat.mean, aes(y = mean.rat, x = Date), size = 0.95, col = "grey50", position = location.move) + geom_point(data = rat.mean,aes( y = mean.rat, col = Rats, shape = Valley, fill = Control, x = Date ),size = 6 , alpha = 0.7, position = location.move) + xlab(expression(paste("Time", "(", italic(t), ")"))) + ylab(expression(atop(paste("Minimum "," ", " number"," "), paste(" ", "of"," ", "rats"," ","(",italic(R[jt]),")"))) ) + scale_color_manual(name = "Stoat Control", values = c("black", "white", "white")) + scale_shape_manual(name = "Ecosystem", values = c(24, 21)) + scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) + # manually define the fill colours scale_fill_manual(name = "Stoat Control", values = c("white", "black", "white")) + xlab(expression(paste("Time", "(", italic(t), ")"))) + ylab(expression(atop(paste("Minimum "," ", " number"," "), paste(" ", "of"," ", "rats"," ","(",italic(R[jt]),")"))) ) + # scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) + # Remove fill legend and replace the fill legend using the newly created size guides( col = "none", size = guide_legend(override.aes = list( shape = c(15,0),alpha = 1 )), shape = guide_legend(override.aes = list( shape = c(24, 21), size = 4 )), fill = guide_legend(override.aes = list( col = c("white", "black"),shape = c("square"), size = 4 )) ) + theme_tufte() + theme_bw() + theme(strip.background = element_blank(), strip.text.y = element_blank(), plot.title = element_text(hjust = 0, size=24, family = "Times", color="black", margin = margin(t = 10, b = 10)), plot.subtitle=element_text(size=16, face="italic", color="black"), legend.position = "none", legend.key = element_blank(), legend.background = element_rect(fill="white", size=1), legend.key.size=unit(1,"cm"), legend.text = element_text(colour = "black", size =16, family = "Times"), legend.title = element_text(colour = "black", size =16, family = "Times"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.spacing = unit(2, "lines"), panel.border = element_blank(), axis.title.y = element_text(colour = "black",size =20, family = "Times", angle = 90), axis.title.x = element_text(colour = "black", size =20, family = "Times"), axis.text.y=element_text(colour = "black",size = 20, family = "Times"), axis.text.x = element_text(colour = "black", size =20, family = "Times"), axis.ticks.x = element_line(size = 1), axis.ticks.y = element_line(size = 1), axis.line.x = element_line(size = 1), axis.line.y = element_line(size = 1), strip.text = element_text(face="bold",colour = "black",size =14, family = "Times")) rat.plot # plot #rat.plot #caption # , fig.cap=cap.rat cap.rat <- c("The avarage number of rat captures ($R_{j,t}$; point estimates) collected each trip during the study period (Autumn 1999 to Winter 2004). The shape of each point distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley) and solid symbols represent cases where stoats are uncontrolled and hollow in valleys with stoat removal.")
#blown up knitr::include_graphics(c("./figs/rats.png"))
All other trips with high rat records were in areas with stoat control. In the Hollyford Valley during both mast-years, rats increased irrespective of stoat control. No rats were captured after May 2002 in the Eglinton Valley.
Numerical response
We found that the intake rate we transformed from the seedfall data $(S_{j,t})$ always had the greatest impact on the rate of increase of mice $(r_{j,k})$. Greater than both mouse density ($N_{j,t}$) and rat presence ($R_{j,t}$; see supplementary material). Our model accounted for the confounding effects of numerical processes such as density ($N_{j,t}$) and food availability ($S_{j,t}$) using our Bayesian hierarchical model and then tested for evidence of mesopredator release of mice during our study. We used the fitted model parameters to compare the differences between the key predictions proposed for areas with and without control after accounting for both the observation error (data collection and population estimates) and process error (population dynamics; @ahrestani2013).
# s.final.model$coefficients[2,4] # s.final.model$coefficients[3,4]
Prediction A
Stoat control does not increase mouse populations at times of low abundance.
We compared the abundance of mice during trips when both mouse abundance and food availability were low ($2002$ and $2003$ using the model. We found no visually observable or statistically significant differences in the estimates between areas with and without stoat control ($1.00$ SE = $1.99$, p-value = $0.62$) during any of the trips during our study.
# summary tables # if rendering in word # , ft.split = TRUE, ft.align = "left" # summary table # flextable::flextable(data.frame(co.effs, summary(low.mod3)$coefficients)) # # # data # source("./R/figures/pa-plot-code.R", echo = FALSE) # # # # # Plot # low.plot.time # + # theme(legend.position = "right") # high.plot.time # gridExtra::grid.arrange(low.plot.time, output.dens, ncol = 2) # Figure Caption cap.Pa <- c("The unmodified abundance of mice in non-mast years (2001, 2003). The shape difference distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley). The grids are divided into the two groups on the x-axis (areas with and without stoat control). Solid symbols represent cases where rats are present and removed at points with hollow symbols.")
# knitr::include_graphics(c("./figs/fig41.png")) # Figure Caption cap.Pa <- c("The unmodified abundance of mice in non-mast years (2001, 2003). The shape difference distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley). The grids are divided into the two groups on the x-axis (areas with and without stoat control). Solid symbols represent cases where rats are present and removed at points with hollow symbols.")
source("./R/data_input.R") source("./R/figures/pa-plot-code.R") low.plot.out
For Prediction A we did not need to account for varying seed input because seed fall was biological equivalent to zero for all grids and trips during low abundance seasons [@choquenot2000].
Prediction B
Stoat control does not result in higher peak abundance of mice
# run source script # source("./R/figures/pb-plot-code.R", echo = FALSE) # summary # s.final.model <- summary(high.mod3) # # s.final.model$coefficients[2,4] # # # output from models in workable format # co.effs.m2 <- c(row.names(as.data.frame(summary(high.mod2)$coefficients))) # co.effs.m3 <- c(row.names(as.data.frame(summary(high.mod3)$coefficients))) # # # results data # pb.results.m2 <- data.frame(co.effs.m2, summary(high.mod2)$coefficients) # pb.results.m3 <- data.frame(co.effs.m3, summary(high.mod3)$coefficients)
Our simple ANOVA model (see here for analysis) accounted for the differences in sample size between the areas with and without stoat control and un-equal variances within groups and found no significant difference between the stoat controlled and uncontrolled areas (What is best here and I can repeat for all?).
# simple anova comparison #flextable::flextable(pb.results.m3)
During the years of peak abundance in our study we did observe lower mouse abundance in the Hollyford Valley compared to the Eglinton Valley ($-0.17$, SE = $0.09$, p-value = $ 0.06 $). Although this is not below the 0.05 threshold. It should be noted that when the rat removal treatment was removed the p-value did drop to 0.05.
# #data # source("./R/figures/pb-plot-code.R", echo = FALSE) #plot # high.plot.time # + # theme(legend.position = "right") # Figure Caption # low.plot.time # gridExtra::grid.arrange(high.plot.time, out.dens.high, ncol = 2) cap.Pb <- c("Difference between areas with and without stoat control at peak mouse abundance. The shape difference distinguishes the two valleys (triangle = Eglinton Valley, circle = Hollyford Valley). The grids are divided into the two groups on the x-axis (areas with and without stoat control). Solid symbols represent cases where rats are present and removed at points with hollow symbols.")
source("./R/figures/pa-topb-test.R") low.plot.out
Prediction C
Stoat control did not increase mouse populations as food becomes available.
We tested the differences in the rate of increase of mouse populations between the four seasons, two valleys and stoat control areas (excluding the rat treatment) during seasons when mouse population were increasing (August; Autumn - Winter). We found was no significant differences between mice populations ($ r_{j,t} $) during the increasing seasons of mouse dynamics. All 95% credible intervals for the differences between these trips included zero (Figure 10).
#rate vs seed data with lines for average prediction source("./R/pc-plot-code.R", echo = FALSE) # # just most interesting variable # # caption cap.Pc = c("Visual representation of the effect of stoat control on the rate of increase of mouse populations during the increasing seedfall (Autumn to Winter; **Prediction C**). The different shapes distinguish the two valley systems (triangles = Eglinton Valley; circles = Hollyford Valley) and all solid symbols represent cases where rats are present and removed at points with hollow symbols.")
pc.plot
# knitr::include_graphics(c("./figs/seeded.png"))
# Save plot # export plot for example vignette # png("C://Code/beech-publication-wr/figs/seeded.png") # pc.plot # dev.off()
Prediction D
Stoat control did not hasten the decline from peak abundance.
# knitr::purl("Davidson_2019_Simple_models.Rmd") # source("./R/Davidson_2019_Model_wrangling.R") # source("./R/figures/pd-plot-code.R") # pd.plot
This prediction specifically focused on the Spring and Summer seasons, when rodent populations are crashing (Figure 3 vs. Figure 4). Overall density had a weak negative effect on mice dynamics under all conditions (Figure 3) but was greatest during the declining seasons of mice dynamics. After accounting for seed and density in the population model there was a faster rate of decline in the Hollyford Valley (large solid circles) compared to the Eglinton Valley (large hollow circles) but this was not statistically significant.
#data for plot # source("./R/figures/rate-density-plot.R", echo = FALSE) #plot # pd.plot.feb + # theme(legend.position = "right") #caption cap.Pd = c("Visual representation of the effect of stoat control during decreasing mouse abundance (Spring to Summer; *Prediction D*). The different shapes distinguish the two valley systems (triangles = Eglinton Valley; circles = Hollyford Valley) and all solid symbols represent cases where rats are present and removed at points with hollow symbols.")
# source("./R/data_input.R") source("./R/pd-plot-code.R", echo = FALSE) pD.plot # export plot for example vignette # png("C://Code/beech-publication-wr/figs/mouse_dens.png") # pD.plot # dev.off()
# knitr::include_graphics(c("./figs/mouse_dens.png"))
Prediction E
Does the presence of rats impact the population dynamics of mice at each of the seasons tested in Predictions A-D.
# knitr::purl("Davidson_2019_Simple_models.Rmd") # source("./R/Davidson_2019_Model_wrangling.R") # source("./R/figures/pe-source-code.R") # pe.plot
We found that rat numbers $R_{j,t}$ had the smallest overall effect on mice dynamics ($ \mu(\beta_3 $ $|$ $All$ $Seasons$$ )=$ $0.007$) and was the most variable parameter in our community model ($ CV_{rats} = 5.74 $). We recorded the highest number of rats in February and May, at the same times of high mouse abundance (Figure 4).
#plot code source("./R/pe-plot-code.R", echo = FALSE) #plot # rat.plot.all # theme(legend.position = "right") #caption pE.cap <- c("Differences in estimates between rats during each season of mice dynamics.") # Save plot # export plot for example vignette # png("C://Code/beech-publication-wr/figs/fig_rat.png") # pe.plot # dev.off()
# knitr::include_graphics(c("./figs/fig_rat.png"))
The relationship between rats and mice was examined by comparing the estimated differences between treatments (Figure 8; lines represent mean relationships in the three treatment groups). We also statistically testing the difference between each different treatment and found no differences across all parameters (Figure 9).
`r
pe.plot
We did not have enough data to estimate a statistical interaction between rats and mice (non-linear relationship) or other more complex non-linear models for rat dynamics.
Discussion {#discussion}
Overall we found no evidence of mesopredator release in the mouse populations we monitored. We used a bayesian mixed model (bmm; Appendix 1) to model mice dynamics in New Zealand Beech forests. We tested five predictions in a relation to stoat control and mouse abundance. We experimentally manipulating predator control in two independent beech forests and compared the differences between the rate of increase and abundance of mice. The model replicated mouse dynamics in our data and produced similar results to other NZ mammal models (@choquenot2000; @ruscoe2001; @ruscoe2005; @tompkins2006; @tompkins2013; @holland2015; summarised in Figure 1). Generally in NZ Beech forests, high food availability results in high mouse abundance in the subsequent months until natural germation of the seed and is commonly observed in both our data and a range of other NZ forest types [@innes; @sweetapple; @ruscoe2011; @ruscoe2012]. Modelling has incorporated these parameter estimates to represent invasive species dynamics in NZ forests [@tompkins2006; @tompkins2013; @holland2015] A naïve investigation of our summarised data (Figure 3.1, 3.2, and 3.3) might have falsely lead to similarly low powered conclusions, that rats may have an effect on mice dynamics as well as many generally accepted relationships. By using a BHM framework we were able to test the theoretical relationships proposed using experimental data.
Seed fall
The intake rate had the greatest impact on mice populations in all seasons. We found that $r_{j,t}$ was strongest on grid and trips when beech seed was increasing. After comparing a selection of functional responses to our data (Appendix 1.3) we found the best type II response and converted our estimated $Seed_{j,k}$ per $m^2$ to the "intake rate" $S_{j,t}$. We found that the effect of mouse density ($N_{j,t}$) and rats ($R_{j,t}$) was lower than the "intake rate" ($S_{j,t}$). We observed that the effects of density were greatest during the declining seasons. These differences in effect sizes are similar to other studies.
Mouse density
The effect of mouse density and rats on was much lower and were also similar to other published results [@choquenot2000]. We and many other studies observed high heterogeneity in abundance estimates. This type of data can be difficult to model correctly in many frequentist frameworks. To account for this heterogeneity in the biological data we used experimental manipulation and Bayesian modelling framework to account for the confounding effects of food availability, mouse density, and rats to simultaneously assess the overall effect of stoat control on mice dynamics. Importantly, our model captured both the observational and process variability within NZ beech forests. We also conclude that the underlying heterogeneity in small mammal populations was sufficiently large to make estimating true response of mice to rats difficult in biological experiments (Figures 4-9).
Rat presence
Our results suggest that it is unlikely that rats have strong impacts on the rate of increase of mice in our large scale study system. Contrary to this, previous laboratory experiments suggested that interactions between mice and rats were negatively related to predation or competition effects and limited evidence of competition release has been observed in mixed podocarp-tawa forests in NZ. We found that rats had the smallest impact on mice dynamics (Figure 8). In part, this is due to the high variability in the number of rats and kiore between grids and sampling trips in our study (??), which is often the case in non-laboratory studies. Instead of rats/kiore being present but at low numbers in all grids (low effect size and low uncertainty), there were only a few grids that had comparability high numbers. This result suggests that refugee areas may exist and can support larger rodents in beech forests that current averages suggest. Overall, it is likely that beech forests do not have enough resources to support larger spatially uniform populations of rats.
Bayesian modelling
The application of BHM's in ecology has been increasing over the past decade [@king2012]. BHM modelling can provide challenges when selecting models and assessing model fit [@auger-methe2016]. To account for this, we used the general associations above as well as AIC and r-squared. We then used these parameters to fit the proposed theoretical population model. We investigated the patterns in our model to insure we adequately encapsulated the patterns in the data with our model to assess the effects of our parametrisation above. Large-scale ecological experiments are very informative but often need to be interpreted with caution due to often unique problems. We address these below and add additional direction to address these uncertainties in the results we found.
Our population model representation is relatively simplistic compared to some of the previous models that incorporated more interacting species and subsequently complex interactions. Our experimental design was just a big taking advantage of the situation to built a simple population model from CR data. A benefit to a simple BHM model is the ability to use bayes theorem to re-assess additional data in context of the proposed predictions. In the context of the data we have used in this experiment and uninformative priors. We choose to focus on high quality data collection methods and the key dynamics already identified with invasive species dynamics in NZ beech forests. There are however several key aspects that may be confounding our results, to the observational and population model choices.
Stoat control
An alternative explanation for the lack of any observed effects of stoat control on mouse abundance was that the removal of animals was potentially limited by the effectiveness of a particular stoat control program, management or experimental design.
If trapping did not remove a sufficient number of individual stoats to reduce predation pressure we will not observe the theoretical responses proposed by @tompkins2013 and others.
High individual heterogeneity in the capture rates of stoats at both high and low mouse densities makes reductions to zero individuals virtually impossible targets. This was the case in our study where by-catch of stoats in controlled areas was observed. As an indication of removal, $4.75$ times more stoats were removed from the stoat control program (Hollyford Valley 2001-2004) than caught in by-catch in areas already being controlled prior to the experiment. Our before/after treatment tested this and verified that we did not find differences within the grid when experimental treatments were changed. We may also be missing the impact of by-catch of other pests in stoat traps. We believe that any effect this may have on mice dynamics would be limited due to the limited impacts of rat removal in our data. Biologically, the differences in home range size of rodent species (10's ha; [@bramley2014; @innes1983; @pryde2005]) relative to the home range of stoats and the scale of trap line (100's ha; [@miller2001; @murphy1995]) which suggests that population level impacts from these removals is unlikely due to the limited number of animals removed and the re-invasion biology of rodents [@bramley2014].
<- is this too critical given ?... ->
[@blackwell2001; @blackwell2003]: However, the limited geographical scale of the study, low statistical power and replication between different forest types makes it difficult to compare the effect of different treatment types during the four different seasons.
Repleccate and extend on larger scale. <- is this too critical given ?... ->
Resource pulse systems
During our study we found stoat control had minimal impact on mice populations compared to beech seed. Beech forests have previously defined as resource pulse systems [@wardle]. Strong resource pulses cause strong effects throughout an ecosystem [@yang2008]. The process behind for this phenomenon is that even in the absence of control, stoats cannot attain high enough densities to limit mouse populations under conditions of high food availability [@king1983]. Furthermore, the intrinsic differences in the reproductive and growth rates (life history traits) of each species reduces the ability of stoats to regulate mouse populations at the beginning of a masting event[@king2011]. By taking the proposed processes that effect mice dynamics and applying a hierarchical model framework allows us to test and estimate effect sizes that further reduce uncertainty in the outcomes of predator control in NZ beech forests.
The effects of climate change and other human lead impacts will have lasting impacts on on future rodent dynamics [@holland2015]. Population models of climate change indicate a shift to prolonged but lower overall mast events where it may be possible to observe mesopredator release of mice with the removal of predators [@tompkins2013]. Mice populations monitored in more complex forest ecosystems suggest complex interactions between top-down and bottom-up effects [@ruscoe2011; @ruscoe2012]. These studies suggest that rodent populations would increase in other NZ forest ecosystems after the removal of stoats [@ruscoe2011] or other invasive predators [@rayner2007].
Challanges
Population level responses such as abundance and rate of increase are often measured using indices (e.g. MNA). This can lead to difficulty in estimating population level parameters such as density. We aimed to avoid this bias as we estimated mice abundance using capture-recapture (CR) data and an integrated population model to correctly account for uncertainty in abundance estimation when fitting the population model of mice dynamics. However due to limited population numbers and data, only indices could be used to represent rat interactions. Independent research in beech forests has shown a high correlation between indices of rats in this valley and CR data. Nevertheless, increasing the quality of the rat/kiore data would most likely reduce uncertainty on the estimates of these parameter estimates. By calculating the maximum increase from the fitted population model we were able to compare the maximum population growth rate to other studies ($R.max_{mice} = actual number here$; @hone2010), however, analogous models are not directly comparable because they do not include a resource component (e.g. beech seed production) or use indices (e.g. $R_{j,t}$). Field studies are also difficult to directly compare because may studies incorporate structurally different ecosystems (e.g @ruscoe2011), species indices instead of true abundance estimation using CR models or did not report effect sizes. Although this limitations exist our comparisons to general trends were all comparable to our model estimates and our estimates are all biologically viable.
Management implications
Our model has subtle but important differences compared previous research when in comes to management application. We have incorporated a reproducible workflow[@britishecologicalsociety2018] for the future development of model testing for different and new datasets [@wickham2014]. In doing this we have incorporated leading reproducible science techniques to help address any aspects of our study may need for further support. Often studies of such large scale can address key population-level questions however these studies struggle to find replication because of the scale at which they are conducted [@oksanen2001]. Many issues that many confront a PFNZ2050 will be reduced by using simple but well fitted BHMs for predicting and allocation management resources.
Predator-free NZ
Our integrated population modelling approach and the large, high-quality CR data used to parametrise the model allowed us to predict and report the biologically relevant effect sizes. We used this to address the direct question of stoat control in NZ Beech Forests. Future studies can now test and assess the beech forest ecosystem compared to the data they collect. By using the unified modelling framework presented in this paper (Appendix 3; [@davidson2019), reproducible research practises and "tidy-data" [@wickham2014]. Many of the processes we have verified were previously hypothesized using indices of mice abundance [@blackwell2001; @blackwell2003], laboratory experiments [@bridgman2013], novel systems such as remote islands [@mulder2009] and small patches of mainland forest [@blackwell2001; @blackwell2003]. We have integrated these results into a framework that can now test these results in other systems that are hypothesized to be different. It will be important for future work to continue monitoring this system for anomalies and unexpected patterns such as increases in other invasive predators, particularly rats.
`r if (knitr:::is_html_output()) '
References {-}
'`
Sorting references has been much hard than I expected when there are multiple types of software to use and different ways to link them.
The core problems I have needed to solve are:
-
Storing pdfs locally and working on them.
-
Mendeley works well for looking and searching through pdfs and was the reason I switched from EndNote (that my supervisors were using at the time) to Mendeley.
-
Convert references attached to pdfs in reference manager into a
bibtexfile so thatcitrandrefmanageRcan read these. -
Zotero because it has good extentions that can turn bib references into bibtex structure for R packages that use a variation on the
.bibstructure in R. -
Convert stable referencing in
.Rmdto other document types and structures. -
This has now become hard to understand how this is all working so I have had to layout a software path (and maybe future Docker image) to ensure that the refences converge and the errors are solvable (or at least I can work out where they come from)
Editor??
To work with .bib files (to find errors to begin with) it is much nice to have some code snippit support (highlighting syntax)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.