Replicating table 1 from Rouder and Morey (2012)

Bailey_Geary_2009.csv

table1.Rmd

---
title: "Replicate Table 1 of Rouder and Morey (2012)"
author: "Richard D. Morey"
date: "24/08/2021"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)

library(dplyr)
library(gt)
library(ggplot2)
library(colorspace)
```


```{r functions}
## Regression functions as they were in 2012
dinvgamma = function (x, shape, scale = 1) 
{
  if (shape <= 0 | scale <= 0) {
    stop("Shape or scale parameter negative in dinvgamma().\n")
  }
  alpha <- shape
  beta <- scale
  log.density <- alpha * log(beta) - lgamma(alpha) - (alpha + 
                                                        1) * log(x) - (beta/x)
  return(exp(log.density))
}

integrand.regression=function(g,N,p,R2)
{
  a=.5*((N-p-1)*log(1+g)-(N-1)*log(1+g*(1-R2)))
  exp(a)*dinvgamma(g,shape=.5,scale=N/2)
}

linearReg.Quad=function(N,p,R2) {
  h=integrate(integrand.regression,lower=0,upper=Inf,N=N,p=p,R2=R2)
  return(h$value)
}

```

```{r readdata}

## Read from public gist
data_files = gistr::gist('https://gist.github.com/richarddmorey/1bca5c91a5a16a100a07644f9be6d7bf')

# I extracted these values using mathpix because I do not 
# have the original .tex file
data_files$files[['table1.tsv']]$content %>%
  textConnection() %>%
  read.delim(.,
           sep = '\t',
           header = TRUE) -> tab1

data_files$files[['Bailey_Geary_2009.csv']]$content %>%
  textConnection() %>%
  read.csv(.,
           header = TRUE) %>%
  rename(
    Density = dpop_30,
    Parasites = tpar,
    Local = diftemp,
    Global = Isosd,
    Capacity = CC
  ) %>%
  select(Density, Parasites, Local, Global, Capacity ) -> b_g

```

## Table as published

```{r tab1.1}
tab1 %>%
  gt(rowname_col = 'abbr') %>%
  cols_label(
    R2 = '\\(R^2\\)',
    Bm0 = '\\(B_{m0}\\)',
    Bmf = '\\(B_{mf}\\)'
  ) %>%
  fmt_scientific(c('Bm0','Bmf')) %>%
  tab_header(title = 'TABLE 1',
             subtitle = 'Bayes Factor Analysis of Hominid Cranial Capacity (Data From Bailey & Geary, 2009)') %>%
  tab_source_note(
    source_note = html('<i>Note</i>. Local = local climate; Global = global temperature; Parasites = parasite load; Density
= population density.')
  ) 
```

## Some checking

I wanted to compare the results to the 2012 code, to see whether the differences could be accounted for by integration tweaks.

```{r}

## Do some checks
tab1 %>%
  pull(Model) %>%
  paste('Capacity ~',.) %>%
  sapply(X = ., FUN = function(fmla){
    lm(as.formula(fmla), data = b_g) %>%
      summary() %>%
      `$`('r.squared')
  }) -> r2s

npars = stringr::str_count(tab1$Model, pattern = '\\+') + 1 

tab1 %>%
  pull(Model) %>%
  paste('Capacity ~',.) %>%
  sapply(X = ., FUN = function(fmla){
    BayesFactor::lmBF(as.formula(fmla), data = b_g, rscaleCont = 1) %>%
      as.vector()
  }) -> bfs

names(bfs) = names(r2s)

bfs2012 = mapply(p = npars,
                 R2 = r2s,
                 MoreArgs = list(N = nrow(b_g)), 
                 FUN = linearReg.Quad)

## These are almost exactly the same
tibble(
  idx = 1:length(bfs),
  bfs = bfs,
  bfs2012 = bfs2012,
  perc_diff = 100*(bfs - bfs2012)/bfs
  ) %>%
  ggplot(aes(x = idx, y = perc_diff)) +
  geom_point() + 
  scale_x_continuous(
    name = 'Model index',
    breaks = 1:length(bfs)
  ) +
  scale_y_continuous(
    name = 'Difference (% of BayesFactor)',
    limits = c(-.001,.001)
  ) +
  ggtitle(label = 'Difference between 2012 BF and BayesFactor') +
  theme_minimal()
  
```

Apparently not. The differences are very small.

## Table as it should be

In the table below, I have computed the values using the `BayesFactor`
 package. I have highlighted the background color of the cells based on their original table's deviation from this table; redder cells are "worse". As you can see, the $\cal{M}_13$ row seems to be the main error, but there are still minor deviations elsewhere.
 
```{r}
tab1 %>%
  mutate(
    Bm0_2021 = bfs,
    Bmf_2021 = bfs / bfs[1],
    err_bm0 = abs(100 * (Bm0 - Bm0_2021)/Bm0_2021),
    err_bmf = abs(100 * (Bmf - Bmf_2021)/Bmf_2021)
    ) -> tab1_corrected

tab1_corrected %>%  
  gt(rowname_col = 'abbr') %>%
  cols_hide(c('Bm0','Bmf','err_bm0','err_bmf')) %>%
  cols_label(
    R2 = '\\(R^2\\)',
    Bm0_2021 = '\\(B_{m0}\\)',
    Bmf_2021 = '\\(B_{mf}\\)'
  ) %>%
  fmt_scientific(c('Bm0_2021','Bmf_2021')) %>%
  tab_header(title = 'TABLE 1',
             subtitle = 'Bayes Factor Analysis of Hominid Cranial Capacity (Data From Bailey & Geary, 2009)') %>%
  tab_source_note(
    source_note = html('<i>Note</i>. Local = local climate; Global = global temperature; Parasites = parasite load; Density
= population density.')
  ) -> gt_obj

# Coloring of cells
# adapted from https://stackoverflow.com/a/63945239/1129889
heat_palette <- leaflet::colorNumeric(
  palette = colorRamp(c("#FFFFFF", "#FF4444"), interpolate = "spline"),
  domain = c(tab1_corrected$err_bm0,
             tab1_corrected$err_bmf)
  )

ht_bm0 <- heat_palette(tab1_corrected$err_bm0)
ht_bmf <- heat_palette(tab1_corrected$err_bmf)

for(i in seq_along(tab1_corrected$err_bm0)) {
   gt_obj <- gt_obj %>%
         tab_style(
           style = cell_fill(color = ht_bm0[i]),
           locations = cells_body(columns = "Bm0_2021", rows = i)
           )      
   gt_obj <- gt_obj %>%
         tab_style(
           style = cell_fill(color = ht_bmf[i]),
           locations = cells_body(columns = "Bmf_2021", rows = i)
         )
   }

gt_obj

```

## Rounded \(R^2\)?

As a check, I wanted to see if the smaller differences could be due to using the rounded \(R^2\) values in the table. Maybe we wanted people to be able to type the value from the table into a calculator and get the Bayes factor from the table.

```{r}
bfs_r2tab = mapply(p = npars,
                 R2 = tab1$R2,
                 MoreArgs = list(N = nrow(b_g)), 
                 FUN = linearReg.Quad)

tab1 %>%
  mutate(
    Bm0_rounded = bfs_r2tab,
    Bmf_rounded = bfs_r2tab / bfs_r2tab[1],
    err_bm0 = abs(100 * (Bm0 - Bm0_rounded)/Bm0_rounded),
    err_bmf = abs(100 * (Bmf - Bmf_rounded)/Bmf_rounded)
    ) -> tab1_rounded

tab1_rounded %>%  
  gt(rowname_col = 'abbr') %>%
  cols_hide(c('Bm0','Bmf','err_bm0','err_bmf')) %>%
  cols_label(
    R2 = '\\(R^2\\)',
    Bm0_rounded = '\\(B_{m0}\\)',
    Bmf_rounded = '\\(B_{mf}\\)'
  ) %>%
  fmt_scientific(c('Bm0_rounded','Bmf_rounded')) %>%
  tab_header(title = 'TABLE 1',
             subtitle = 'Bayes Factor Analysis of Hominid Cranial Capacity (Data From Bailey & Geary, 2009)') %>%
  tab_source_note(
    source_note = html('<i>Note</i>. Local = local climate; Global = global temperature; Parasites = parasite load; Density
= population density.')
  ) -> gt_obj

# Coloring of cells
# adapted from https://stackoverflow.com/a/63945239/1129889
heat_palette <- leaflet::colorNumeric(
  palette = colorRamp(c("#FFFFFF", "#FF4444"), interpolate = "spline"),
  domain = c(tab1_rounded$err_bm0,
             tab1_rounded$err_bmf)
  )

ht_bm0 <- heat_palette(tab1_rounded$err_bm0)
ht_bmf <- heat_palette(tab1_rounded$err_bmf)

for(i in seq_along(tab1_rounded$err_bm0)) {
   gt_obj <- gt_obj %>%
         tab_style(
           style = cell_fill(color = ht_bm0[i]),
           locations = cells_body(columns = "Bm0_rounded", rows = i)
           )      
   gt_obj <- gt_obj %>%
         tab_style(
           style = cell_fill(color = ht_bmf[i]),
           locations = cells_body(columns = "Bmf_rounded", rows = i)
         )
   }

gt_obj


```

Using rounded $R^2$ seems to account for the smaller differences in the table. 

table1.tsv