In cReedHranac/batwintor: Tool for bat hibernation energetics

library(knitr)
opts_chunk$set(tidy.opts=list(width.cutoff=80),tidy=TRUE)
knitr::opts_chunk$set(echo = TRUE)

This vignette is designed to show how the batwintor package works from start to finish using real data. Metabolic, and morphometric data was collected over 2 winters in British Columbia, Canada.

From Raw

The raw data for this analysis is available through (via)

library(batwintor);library(data.table);library(dplyr)
dat.dirty <- read.csv("../data-raw/RawData.csv")
dat <- tbl_df(dat.dirty)
names(dat) <- c("date", "bat", "species", "sex", "band", "mass.in", "mass.out", "Ta", "Ta.p",
                "flow", "FiO2", "FeO2", "VO2.ml.min", "VO2.ml.h", "VO2.ml.h.g")
dat$Ta <- as.factor(dat$Ta)

Now that the data it in it's always a good idea to visualize some of it and make sure that what you think is going in, actually is. Models are only ever as good as the data that you put into them.

myca <- filter(dat, dat$species == "MYCA")

myyu <- filter(dat, dat$species == "MYYU")

hist(myca$mass.in);hist(myyu$mass.in)

Those look close enough to normal, and indeed if you run normality tests shapiro.test() they are. Another very important parameter is the actual metabolic rate (MR) which is labeled in dat as VO2.ml.h.g.

myca.6 <- filter(myca, myca$Ta == 6) # For trials conducted at 6 degrees
myca.2 <- filter(myca, myca$Ta == 2) # For trials conducted at 2 degrees
hist(myca.6$VO2.ml.h.g)
hist(myca.2$VO2.ml.h.g)

As you can see the trials at 6 degrees look pretty normal, however the trials at 2 degrees seem to have two clusters suggesting that not all of the bats fully thermoconformed implying that there may be some outliers in our data. To handle this we use the MRFromRaw function which fits semi-parametric models using the mixtools package and from that we can visualize all of our temperature conditions, as well as letting the function pull out the resting metabolic rate (RMR) estimate, the minimum torpor metabolic rate (TMRmin) estimate, what temperature that rate was reached, and the mass estimate. These parameters will be used throughout the rest of the modeling process.

myca.mr <- MRFromRaw(dat, "MYCA")
myyu.mr <- MRFromRaw(dat, "MYYU")

Once we have the estimates we can update the rest of the parameters that are stored internally, and are generally conserved between runs.

data("bat.params")
myca.params <- MRUpdate(myca.mr, "M.californicus", params = bat.params)
myyu.params <- MRUpdate(myyu.mr, "M.yumanensis", params = bat.params)

Now our bat metabolic parameters are updated with our estimates drawn from the raw data and we're ready to begin the modeling portion of the analysis.

Fungal Options

An important factor of this mechanistic model is the inclusion of fungal growth parameters drawn from the literature (see ?LoadFung and `data("fung.params")' for more information). Since the fungal growth parameters are generally updated much less frequently than our metabolic estimates they are provided with the package, and the calls are simple. In this case we'll chose the growth rate documented by Chaturveti et al. which is the more aggressive of the two growth options available.

fung.ch <- FungSelect("Chaturvedi")

Model Environment

In order to run the model we need to create an environmental space across which to model. This requires either vecotrs or rasters the describe the the range of humidity and temperature values to run the model across, as well as the resolution of that range. This step can be a major determinate for how long a model run will take. Temperature and humidity ranges extracted from environmental rasters will generally be quite large, and they may extend past conditions that would be biologically relivant (ie. outside of a hibernaculum), however if only range vectors are passed extreme conditions may be excluded.

env <- BuildEnv(temp = c(-5,20), hum = c(60,100), range.res = 1)

Model Run

The final thing before running the model is to determine the maximal length of winter to run the model accross. In this example we will say 9 month, and create an hourly vecor so we can have good resolution of for the differential equations. This step can be another that changes the the calculation length but it is generally suggested that your twinter vecter is at least daily. For the final exampel we will use the built in data set for Myotis lucifigus one of the species most heavilly effected by the WNS epidemic.

twinter <- seq(from = 0, to = 9*24*30*1, by = 24) #Maximum time of winter is 9 months, broken at daily intervals 
data("mylu.params")
# system.time(
# mylu.mod <- DynamicEnergyPd(env, mylu.params, fung.ch)  
# )

This model with the parameters described above benchmarks in at ~359 seconds on run-of-the-mill laptop but will run much faster on most desktops and is generally all the resolution that you'd need however if you want to explore parameter space you may want to keep computational time in mind.

cReedHranac/batwintor documentation built on Jan. 27, 2020, 7:39 p.m.