Seasonal temperature acclimatization in a semi-fossorial mammal and the role of burrows as thermal refugia
Milling CR, Rachlow JL, Chappell MA, Camp MJ, Johnson TR, Shipley LA, Paul DR, Forbey JS
Peer J

Abbreviations and descriptions of data:
Milling_etal_EnviroTemps.csv

DateTime		Observation timestamp
Season			winter or summer
Tbrw			Mean burrow temperature (degrees C)
Tbrw_Cost		Predicted RMR (mL O2/min) for ave. size animal at rest in burrow (predicted from nlme regression)
Mean_Te			Mean above-ground operative temperature (degrees C)
Te1_Cost		Predicted RMR (mL O2/min) for ave. size animal at rest above ground at Mean_Te (predicted from nlme regression for Mean_Te =<35 degrees C;
			When Mean_Te>35, Te1_Cost = 0.209286Te - 2.54501 -- see manuscript for justification)
Min_Te			Minimum above-ground operative temperature (degrees C)
Te2_Cost		Predicted RMR (mL O2/min) for ave. size animal at rest above ground, assuming animal is capable of behavioral thermoregulation. 
			When Mean_Te > 35, RMR recalculated using TeMin (predicted from nlme regression for Min_Te =< 35 degrees C; when Min_Te>35,
			Te2_Cost = 0.209286Te - 2.54501)
BurrowIsRefuge		TRUE if the burrow satisfies the definition of thermal refuge (winter: Te < Tbrw; summer: Te > 35)
Opt_Cost		Predicted RMR (mL O2/min) for ave. size animal that uses a burrow according to definition of refuge (if BurrowIsRefuge == TRUE, cost = Tbrw_Cost;
			If BurrowIsRefuge == blank, cost = Te1_Cost [ winter] or Te2_cost [summer])
mass			Mass of average size animal used in this study (462 g)

Milling_etal_Temps_ABBR.csv
Abbreviated file of the same information contained in Milling_etal_EnviroTemps.csv for calculating SE of predicted RMR
Using methods outlined in Mandel (2013).

season			winter or summer
Microhabitat		Burrow, Above ground, or both (refuge according to riles in Milling_etal_EnviroTemps.csv)
Tc			Mean temperature (degrees C) of the microhabitat for the observation
mass			Mass of average size animal used in this study (462 g)

Milling_etal_nlmePYRA.R
Data and R-code for performing non-linear mixed effects segmented regression on respirometry data. 
Annotated code appears below and can be copied/pasted into R.

rabbits <- read.table(header = TRUE, text = "
animal	season	mass	Ta	Tc	VO2
eya	summer	479	5	6.43	6.36
eya	summer	481	10	12.86	8.52
eya	summer	480	15	17.06	7.41
eya	summer	481	20	22.27	4.95
eya	summer	479	25	25	5.05
eya	summer	485	30	30.47	5.68
olive	summer	447	20	21.48	5.03
olive	summer	453	5	5.79	9.34
olive	summer	448	10	12.24	7.67
olive	summer	448	15	16.93	5.62
olive	summer	444	25	24.59	3.61
olive	summer	461	30	30.43	6.13
scout	summer	450	5	7.06	9.86
scout	summer	469	10	12.15	5.74
scout	summer	476	15	17	8.66
scout	summer	446	20	21.14	7.06
scout	summer	466	25	24.78	6.16
scout	summer	468	30	30.14	5.73
taylor	summer	501	5	6.43 NA
taylor	summer	504	10	12.54	8.23
taylor	summer	494	15	17.25	7.02
taylor	summer	506	20	21.49	5.22
taylor	summer	525	25	24.75	4.8
taylor	summer	496	30	30.26	4.76
wicket	summer	495	5	6.43	8.05
wicket	summer	493	10	12.86	9.56
wicket	summer	490	15	17.06	7.58
wicket	summer	496	20	21.73	5.33
wicket	summer	495	25	25.5	3.83
wicket	summer	495	30	29.99	3.34
winston	summer	430	5	6.43	6.58
winston	summer	425	10	12.88	8.37
winston	summer	420	15	17.06	5.23
winston	summer	433	20	21.47	4.1
winston	summer	424	25	24.76	2.54
winston	summer	419	30	29.94	4.2
bella	winter	437	-5	-3.04	10.18
bella	winter	430	0	3.37	9.37
bella	winter	437	5	6.88	7.17
bella	winter	447	10	10.87	8.35
bella	winter	440	15	15.73	8.24
bella	winter	436	20	19.66	8.03
bella	winter	437	25	25	7.14
cinna	winter	384	0	1.86	5.97
cinna	winter	382	5	6.55	6.03
cinna	winter	391	15	15.61	5.36
cinna	winter	392	20	19.5	4.6
cinna	winter	384	25	24.35	4.31
prim	winter	404	-5	-4.77	7.73
prim	winter	407	0	2.25	8.01
prim	winter	420	5	7.9	8.14
prim	winter	402	10	10.95	4.58
prim	winter	420	15	15.61	8.17
prim	winter	415	20	19.61	7.61
prim	winter	421	25	25.73	5.18
Scout	winter	519	-5	-5.02	7
Scout	winter	513	0	1.37	9.55
Scout	winter	559	5	6.71	8.36
Scout	winter	558	10	11.07	7.62
Scout	winter	567	15	15.39	8.37
Scout	winter	567	20	18.5	7.05
Scout	winter	514	25	24.67	4.96
taylor	winter	485	-5	-4.52	9.5
taylor	winter	468	0	0.92	6.82
taylor	winter	466	5	5.55	7.12
taylor	winter	485	10	10.06	4.34
taylor	winter	481	15	15.27	4.51
taylor	winter	485	20	19.1	4.14
taylor	winter	477	25	24.87	6.65
winston	winter	459	-5	-4.32	9.12
winston	winter	457	0	2.62	8.89
winston	winter	465	5	5.73 NA
winston	winter	454	10	10.87	7.91
winston	winter	475	15	15.84	7.19
winston	winter	455	20	18.56	4.91
winston	winter	462	25	24.78	3.92
")

library(ggplot2)
library(nlme)

rabbits <- subset(rabbits, !is.na(VO2)) # remove the two cases with missing values

# Model allowing mass to affect the function height (b0) and season to affect height and slope, but not knot.
# The model has three parameters: the height of the flat part after the knot (b0), the slope of the part before the knot (b1),
# and the knot itself (knot). These parameters could vary over other variables like season, mass, and/or animal.
mymodel <- nlme(VO2 ~ b0 + b1*(Tc < knot)*(Tc - knot),
                data = rabbits,
                fixed = list(b0 ~ mass + season, b1 ~ season, knot ~ 1), # b0 varies by mass & season, and b1 by season
                random = b0 ~ 1,
                groups = ~ animal,
                start = c(5, 0, 0, -0.4, 0, 25)) 
# Note that if you look at summary(mymodel) you'll see inferences for b1.seasonwinter, which represents the 
# effect of season on the line slope (which is the b1 parameter). So, for example, b1.(Intercept) is the slope for summer
# and b1.(Intercept) + b1.seasonwinter is the slope for winter.

#To obtain inferences for the slope during winter we
# reparameterize the model by changing the reference category for season as shown below.
rabbits$season <- relevel(rabbits$season, ref = "winter")
mymodel <- nlme(VO2 ~ b0 + b1*(Tc < knot)*(Tc - knot),
                data = rabbits,
                fixed = list(b0 ~ mass + season, b1 ~ season, knot ~ 1), # b0 varies by mass & season, and b1 by season
                random = b0 ~ 1,
                groups = ~ animal,
                start = c(5, 0, 0, -0.4, 0, 25)) 
summary(mymodel)
# Here inferences concerning the slope during winter are listed by b1.(Intercept). 

# To get the height of the line for summer past the knot (i.e., the flat line) for a given but arbitrary value of mass, 
#we can center mass at that value (here 462g which was the mean mass of the animals used). 
#Note switching back to original parameterization. 
rabbits$season <- relevel(rabbits$season, ref = "summer")
mymodel <- nlme(VO2 ~ b0 + b1*(Tc < knot)*(Tc - knot),
                data = rabbits,
                fixed = list(b0 ~ I(mass - 462) + season, b1 ~ season, knot ~ 1), # b0 varies by mass & season, and b1 by season
                random = b0 ~ 1,
                groups = ~ animal,
                start = c(5, 0, 0, -0.4, 0, 25)) 
summary(mymodel)
# Here b0.(Intercept) is the estimated height of the line for summer past the knot for an average animal of 462g.  

### To evaluate model fit:###
# The following gives standardized residuals.
resid(mymodel, level = 0:1, type = "pearson")
# As a rough assessment of fit, examining the proportion of residuals at each level that are greater than two in absolute value.
apply(resid(mymodel, level = 0:1, type = "pearson"), 2, function(x) mean(abs(x) > 2))
# Assuming these residuals have approximately standard normal distributions, only about 5% should exceed 2 in absolute value 
#(this is only a rough approximation, in part due to the dependence of the residuals.


### ANCOVA with random effect of individual to evaluate effect of season on minimum thermal conductance ###

conductance<- read.table(header = TRUE, text = "
animal	mass	season	C
eya	479	summer	11.9736429
olive	453	summer	17.23777302
scout	450	summer	18.93725994
wicket	495	summer	15.15531846
winston	430	summer	12.38782554
bella	437	winter	14.7750363
prim	404	winter	10.76851638
scout	519	winter	9.69529086
taylor	485	winter	13.31153664
winston	459	winter	12.8390427
")

library(lme4)
mymod<- lmer(C~ mass + season + (1|animal), data = conductance) # Linear form of ANCOVA, with random effect for individual (repeated obs on 2 animals)
summary(mymod)
confint(mymod, level =0.95) #lmer does not compute p-values -- use 95% CI overlapping 0

### To estimate thermoregulatory costs above ground, in burrows, both (refuge) ###
# Bootstrapping to obtain standard errors.

d <- read.csv("Milling_etal_Temps_ABBR.csv", header = TRUE) # set path or working directory as needed


library(MASS)

B <- 1000 # number of bootstrap samples

mb <- mymodel$coefficients$fixed # mean vector for bootstrap samples
sb <- vcov(mymodel) # covariance matrix for bootstrap samples

m <- mymodel # copy model for changing parameter values for prediction
y <- rep(NA, B) # hold simulated predicted values

# create data frame to hold results
out <- expand.grid(season = c("summer","winter"), habitat = c("Above","Burrow","Refuge"), estimate = NA, se = NA)

for (i in c("summer","winter")) {
  for (j in c("Above","Burrow","Refuge")) {
    estimate <- sum(1.206 * predict(mymodel, newdata = subset(d, season == i & Microhabitat == j), level = 0))
    for (b in 1:B) {
      beta <- mvrnorm(1, mb, sb) # generate bootstrap sample
      m$coefficients$fixed <- beta # replace parameter estimates
      y[b] <- sum(1.206 * predict(m, newdata = subset(d, season == i & Microhabitat == j), level = 0)) # predicted cost
    }    
    out$estimate[out$season == i & out$habitat == j] <- estimate # point estimate
    out$se[out$season == i & out$habitat == j]<- sqrt(sum((y - estimate)^2)/B) # estimated standard error
  }
}

print(out)