Does brain mass predict how much mammals sleep in a day?

Setup

Load packages

library(dplyr, warn.conflicts = FALSE)
library(ggplot2)
theme_set(theme_minimal())
library(ggrepel)
library(rstan)
library(rstanarm)
library(reshape2)

1 Introduction

This notebook was inspired by a Tristan Mahr's notebook analysing whether brain mass predicts how much mammals sleep in a day. Tristan's original model had the problem that it would predict sleep times over 24h per day.

2 Mammal sleep data

Let’s use the mammal sleep dataset from ggplot2. This dataset contains the number of hours spent sleeping per day for 83 different species of mammals along with each species’ brain mass (kg) and body mass (kg), among other measures. Here’s a first look at the data. Preview sorted by brain/body ratio. The sorting was chosen so that humans show up in the preview.

msleep %>% 
  select(name, sleep_total, brainwt, bodywt, everything()) %>% 
    arrange(desc(brainwt / bodywt))

# A tibble: 83 x 11
   name  sleep_total brainwt bodywt genus vore  order conservation sleep_rem
   <chr>       <dbl>   <dbl>  <dbl> <chr> <chr> <chr> <chr>            <dbl>
 1 Thir…        13.8 4.00e-3  0.101 Sper… herbi Rode… lc                 3.4
 2 Owl …        17   1.55e-2  0.48  Aotus omni  Prim… <NA>               1.8
 3 Less…         9.1 1.40e-4  0.005 Cryp… omni  Sori… lc                 1.4
 4 Squi…         9.6 2.00e-2  0.743 Saim… omni  Prim… <NA>               1.4
 5 Maca…        10.1 1.79e-1  6.8   Maca… omni  Prim… <NA>               1.2
 6 Litt…        19.9 2.50e-4  0.01  Myot… inse… Chir… <NA>               2  
 7 Gala…         9.8 5.00e-3  0.2   Gala… omni  Prim… <NA>               1.1
 8 Mole…        10.6 3.00e-3  0.122 Spal… <NA>  Rode… <NA>               2.4
 9 Tree…         8.9 2.50e-3  0.104 Tupa… omni  Scan… <NA>               2.6
10 Human         8   1.32e+0 62     Homo  omni  Prim… <NA>               1.9
# … with 73 more rows, and 2 more variables: sleep_cycle <dbl>, awake <dbl>

Choose animals with known average brain weight, and add some transformed variables.

msleep <- msleep %>% 
  filter(!is.na(brainwt)) %>% 
  mutate(log_brainwt = log10(brainwt), 
         log_bodywt = log10(bodywt), 
         log_sleep_total = log10(sleep_total),
         logit_sleep_ratio = qlogis(sleep_total/24))

Make a list of examples and give some familiar species shorter names

ex_mammals <- c("Domestic cat", "Human", "Dog", "Cow", "Rabbit",
                "Big brown bat", "House mouse", "Horse", "Golden hamster")
renaming_rules <- c(
  "Domestic cat" = "Cat", 
  "Golden hamster" = "Hamster", 
  "House mouse" = "Mouse")

ex_points <- msleep %>% 
  filter(name %in% ex_mammals) %>% 
  mutate(name = stringr::str_replace_all(name, renaming_rules))

Define these labels only once for all the plots

lab_lines <- list(
  brain_log = "Brain mass (kg., log-scaled)", 
  sleep_raw = "Sleep per day (hours)",
  sleep_log = "Sleep per day (log-hours)"
)

Plot sleep times vs. average brain weights

ggplot(msleep) + 
  aes(x = brainwt, y = sleep_total) + 
  geom_point(color = "grey40") +
  # Circles around highlighted points + labels
  geom_point(size = 3, shape = 1, color = "grey40", data = ex_points) +
  geom_text_repel(aes(label = name), data = ex_points) + 
  # Use log scaling on x-axis
  scale_x_log10(breaks = c(.001, .01, .1, 1)) + 
    labs(x = lab_lines$brain_log, y = lab_lines$sleep_raw)

3 Regression model

Next we use stan_glm from rstanarm package to make a linear model for logit of the sleep ratio given log of the brain weight (Tristan made the model with untransformed variables).

m1 <- stan_glm(
  logit_sleep_ratio ~ log_brainwt, 
  family = gaussian(), 
  data = msleep, 
  prior = normal(0, 3),
  prior_intercept = normal(0, 3))

Check inference

monitor(m1$stanfit)

Inference for the input samples (4 chains: each with iter = 2000; warmup = 0):

                 Q5   Q50   Q95  Mean  SD  Rhat Bulk_ESS Tail_ESS
(Intercept)    -1.5  -1.2  -0.9  -1.2 0.2     1     3726     2756
log_brainwt    -0.6  -0.5  -0.3  -0.5 0.1     1     3894     2578
sigma           0.6   0.7   0.8   0.7 0.1     1     3588     2849
mean_PPD       -0.6  -0.4  -0.1  -0.4 0.1     1     3981     3670
log-posterior -65.8 -62.9 -61.8 -63.2 1.3     1     1564     2564

For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
effective sample size for bulk and tail quantities respectively (an ESS > 100 
per chain is considered good), and Rhat is the potential scale reduction 
factor on rank normalized split chains (at convergence, Rhat <= 1.05).

Prepare \(x\) values for prediction:

x_rng <- range(msleep$log_brainwt) 
x_steps <- seq(x_rng[1], x_rng[2], length.out = 80)
new_data <- data_frame(
  observation = seq_along(x_steps), 
  log_brainwt = x_steps)

Predict expected sleep time at new x values:

preds<-posterior_linpred(m1,newdata=new_data)
preds<-plogis(preds)*24

Plot draws of the expected sleep time lines:

gg<-data.frame(log_brainwt=new_data$log_brainwt,preds=t(preds[1:400,]))
gg<-melt(gg,id=c("log_brainwt"))
names(gg)<-c("log_brainwt","pp","preds")

# aesthetic controllers
alpha_level <- .15
col_draw <- "grey60"
col_median <-  "#3366FF"

ggplot(msleep) + 
  aes(x = log_brainwt, y = sleep_total) + 
  # Plot a random sample of rows as gray semi-transparent lines
  geom_line(aes(x=log_brainwt, y=preds, group=pp), 
              data = gg, color = col_draw, 
              alpha = alpha_level) + 
  geom_point() + 
  scale_x_continuous(labels = function(x) 10 ^ x) +
  labs(x = lab_lines$brain_log, y = lab_lines$sleep_raw)

Predict distribution of sleep times at new x values:

preds_post <- posterior_predict(m1, newdata = new_data)
preds_post<-plogis(preds_post)*24

Plot distribution of sleep times at new x values:

pq<-data.frame(t(apply(t(preds_post), 1, quantile, probs = c(0.025, 0.5, 0.995), na.rm = TRUE)))
names(pq)<-c("lower","median","upper")
pq$log_brainwt<-new_data$log_brainwt

ggplot(msleep) + 
  aes(x = log_brainwt) + 
  geom_ribbon(aes(ymin = lower, ymax = upper), data = pq, 
              alpha = 0.4, fill = "grey60") + 
  geom_line(aes(y = median), data = pq, colour = "#3366FF", size = 1) + 
  geom_point(aes(y = sleep_total)) + 
  scale_x_continuous(labels = function(x) 10 ^ x) +
      labs(x = lab_lines$brain_log, y = lab_lines$sleep_raw)

Modeling the logit of the sleep ratio, and the transforming back to hours of sleep, keeps the distribution of the sleep times restricted to be between 0 and 24 hours.

Licenses

• Code & Text © 2016, Tristan Mahr, 2017-2018, Aki Vehtari, licensed under MIT

Tristan's copyright:

The MIT License (MIT)

Copyright (c) 2016 TJ Mahr

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Original Computing Environment

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.7 LTS

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

locale:
 [1] LC_CTYPE=en_US.utf8       LC_NUMERIC=C             
 [3] LC_TIME=en_US.utf8        LC_COLLATE=en_US.utf8    
 [5] LC_MONETARY=en_US.utf8    LC_MESSAGES=en_US.utf8   
 [7] LC_PAPER=fi_FI.utf8       LC_NAME=C                
 [9] LC_ADDRESS=C              LC_TELEPHONE=C           
[11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C      

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] reshape2_1.4.4       rstanarm_2.21.1      Rcpp_1.0.5          
[4] rstan_2.21.2         StanHeaders_2.21.0-5 ggrepel_0.8.2       
[7] ggplot2_3.3.2        dplyr_1.0.1         

loaded via a namespace (and not attached):
 [1] splines_3.5.1      jsonlite_1.7.0     gtools_3.8.2       RcppParallel_5.0.2
 [5] threejs_0.3.3      shiny_1.5.0        assertthat_0.2.1   statmod_1.4.34    
 [9] stats4_3.5.1       yaml_2.2.1         pillar_1.4.6       lattice_0.20-41   
[13] glue_1.4.1         digest_0.6.25      promises_1.1.1     minqa_1.2.4       
[17] colorspace_1.4-1   Matrix_1.2-18      htmltools_0.5.0    httpuv_1.5.4      
[21] plyr_1.8.6         dygraphs_1.1.1.6   pkgconfig_2.0.3    purrr_0.3.4       
[25] xtable_1.8-4       scales_1.1.1       processx_3.4.3     later_1.1.0.1     
[29] lme4_1.1-23        tibble_3.0.3       farver_2.0.3       bayesplot_1.7.2   
[33] generics_0.0.2     ellipsis_0.3.1     DT_0.15            withr_2.2.0       
[37] shinyjs_1.1        cli_2.0.2          survival_3.2-3     magrittr_1.5      
[41] crayon_1.3.4       mime_0.9           evaluate_0.14      ps_1.3.3          
[45] fansi_0.4.1        nlme_3.1-148       MASS_7.3-51.6      xts_0.12-0        
[49] pkgbuild_1.1.0     colourpicker_1.0   rsconnect_0.8.16   tools_3.5.1       
[53] loo_2.3.1          prettyunits_1.1.1  lifecycle_0.2.0    matrixStats_0.56.0
[57] stringr_1.4.0      V8_3.2.0           munsell_0.5.0      callr_3.4.3       
[61] compiler_3.5.1     rlang_0.4.7        nloptr_1.2.2.1     grid_3.5.1        
[65] ggridges_0.5.2     htmlwidgets_1.5.1  crosstalk_1.1.0.1  igraph_1.2.5      
[69] miniUI_0.1.1.1     labeling_0.3       base64enc_0.1-3    rmarkdown_2.3     
[73] boot_1.3-25        gtable_0.3.0       codetools_0.2-16   inline_0.3.15     
[77] curl_4.3           markdown_1.1       R6_2.4.1           gridExtra_2.3     
[81] rstantools_2.1.1   zoo_1.8-8          knitr_1.29         utf8_1.1.4        
[85] fastmap_1.0.1      shinystan_3.0.0    shinythemes_1.1.2  stringi_1.4.6     
[89] parallel_3.5.1     vctrs_0.3.2        tidyselect_1.1.0   xfun_0.16