In vullioud/torpor: Assigning metabolic rate measurements to torpor and euthermia in heterothermic endotherms

How to use torpor ? An example

This vignette presents a suggested workflow of the torpor package. This package is firstly aimed at users not familiar with rjags, yet interested in investigating the thermoregulatory pattern of their model species, especially when the latter is a heterotherm and when metabolic rate measurements are to be assigned to either torpor or euthermia.

This workflow is based on one of the provided datasets. We go through a complete analysis to illustrate how to use the package functions. At each step, we also present some options for the users and the possible solution to expected problems. Note, however, that this example does not present all aspects of the package.

Installing the package

As the package is still not available on CRAN, users can download it by calling the following command.

remotes::install_github("vullioud/torpor", build_vignettes = TRUE ,force=TRUE)

The program JAGS can be downloaded from the website: https://mcmc-jags.sourceforge.io/

Fitting a model

The first step in any analysis is to import the data. In order to fit a model using the function tor_fit(), we need a set of measurements of Metabolic rate $M$ and of the ambient temperature $T_a$ at which the measurements have been recorded. These values can be vectors of the same length or two columns of a data frame. In the chosen example on the Tasmanian pygmy possum, Cercartetus lepidus, we have 103 metabolic rate values at various $T_a$ by digitalization of a figure (Geiser, 1987).

The data are accessible with the following command:

library(torpor)
data(test_data2)
str(test_data2)

The lower critical temperature ($T_{lc}$) of the thermoneutral zone ($TNZ$) and the metabolic rate within $TNZ$ ($M_{TNZ}$) can be estimated by or provided to the tor_fit() function. These two values were originally estimated by the author and can be found in the documentation of the dataset (test_data2) for the present example, we will estimate $T_{lc}$ and provide $M_{TNZ}$. tor_fit() represents the core function of the model and should be the first step in any analysis using the torpor package.

The model is fitted with the following call:

model <-tor_fit(Ta = test_data2$Ta, 
                M = test_data2$VO2ms,
                Mtnz = 1.8,
                fitting_options = list(ni = 500,
                                            nt = 2,
                                            nb = 200,
                                            nc = 2,
                                            parallel =FALSE), 
                confidence = 0.5)

Note that for the purpose of this vignette we only run two chains with 500 iterations and 200 burned-in. The results have to be considered carefully (i.e. Checking the adequacy of the model, below). Default settings are more suitable but require more time for the model to run.

The output of the tor_fit() unction is a list containing information on the posterior distributions of the investigated parameters as well as on the convergence of the chains and on the prior-posterior distributions overlap of some parameters. These information can also be called by the function tor_summarise(). The rest of the package offers several functions to make more sense of the model output, to check its consistency and to represent it graphically. Researchers familiar with jagsUI can develop their analysis from that point.

Checking the adequacy of the model

Once the model is fitted, it is recommended to verify the adequacy of the results. This can be achieved via the function tor_summarise(). The latter will return a list with the essentials: A dataframe with the mean and median of parameter estimates, the 95% credible interval and the $\hat{R}$ value (i.e. chain convergence estimation). It also reports the parameters’ identifiability.

summary <- tor_summarise(model)

The parameter estimates are accessible with:

summary$params

Checking the convergence

If the $\hat{R}$ given in the summary dataframe are larger than 1.1, we recommend to refit the model and increase the number of iterations and burn-ins respectively. This can be done by setting the parameter fitting_options = list(ni = , nb = ). It is also possible to look graphicaly at the convergence using the jagsUI function traceplot(). Let’s have a look at the convergence for the parameters betat and betac. The necessary object to run the jagsUI::traceplot() function can be accessed through model$mod_parameter.

jagsUI::traceplot(model$mod_parameter, c("betat", "betac"))

checking the parameters identifiability

In addition to a verification of the convergence it is advised to control the identifiability of some parameters. This is done by comparing the prior and the posterior distributions. A warning will be sent to the user for estimated parameters whose PPO are higher than 75% (Fasel et al. in prep.). Parameter identifiability is provided for $T_{lc}$, $M_r$, and $T_{be}$, $TMR$ and $T_{bt}$. The ovelap is also given by the tor_summarise() function. Let’s continue with our analysis of Cercartetus lepidus:

summary$ppo

Making sense of the model

Once the basic checks have been done, we can go on with the evaluation of the output. There are two main functions that deal with predictions. tor_assign() firstly gives the classification of the raw data based on the predictions of the model. tor_predict() further gives the predicted M in torpor and euthermia for a given $T_a$. Finally, we can plot the result using the function tor_plot().

Getteing the assignation

To look at the assignments of the data the corresponding predicted values we use the function tor_assign(), , which will assign the measured metabolic values to either torpor or euthermia and returns a dataframe with the measured $M$, the measured $T_a$, the predicted $M$ and the assignment. (predicted_state).

assignment <- tor_assign(tor_obj = model)

head(assignment)

Prediction

The tor_predict() function is slightly different as it takes a vector of $T_a$ as input and returns the predicted $M$ both in torpor and in euthermia and the 95% credible interval. For example, let’s see what are the predicted $M$ at $T_a$ 22°!

prediction <- tor_predict(tor_obj = model, Ta = 22, CI = TRUE)
head(prediction)

The computation of the credible interval is slow for large jobs and running tor_predict() with the parameter CI = FALSE drastically speed up the computation of the prediction.

plotting the data and predicted values

Finally a built-in function allows plot the results. The user can modulate labels xlab and ylab and the colors with col_torp, col_eut and col_Mtnz and can save the plot using the arguments pdf = TRUE. The plot_type argument also allows to choose between base (default) and ggplot.