inst/scripts/table_relaxed.R

In joshwlambert/relaxedDAISIE: Reproduce Relaxed-rate DAISIE Paper

# script to produce tables that can be copied and modified in latex

data("Canaries")
data("Comoros")
data("Galapagos")
data("Hawaii")
data("Marquesas")
data("New_Caledonia")
data("SaoTome_Principe")

tab_list <- list()
models <- c("cr_dd", "rr_lac_dd", "rr_mu_dd", "rr_k", "rr_laa_dd")
islands <- c("Canaries", "Comoros", "Galapagos", "Hawaii", "Marquesas",
             "New_Caledonia", "SaoTome_Principe")

for (i in seq_along(islands)) {
  tab <- data.frame()
  for (model in models) {
    best_model <- choose_best_model(data_name = islands[i], model = model)
    best_model <- cbind(best_model, aic = NaN, aicc = NaN)
    best_model[["aic"]] <- calc_aic(
      results = best_model,
      daisie_data = eval(parse(text = islands[i]))
    )
    best_model[["aicc"]] <- calc_aicc(
      results = best_model,
      daisie_data = eval(parse(text = islands[i]))
    )
    if (model == "cr_dd") {
      best_model <- cbind(best_model, sd = NA_real_)
      best_model <-
        best_model[, c(
          "lambda_c", "mu", "K", "gamma", "lambda_a", "sd", "loglik",
          "bic", "aic", "aicc"
        )]
    } else {
      rm_index <- which(colnames(best_model) %in% c("df", "conv"))
      best_model <- best_model[, -rm_index]
    }
    best_model <- cbind(model = model, best_model)
    tab <- rbind(tab, best_model)
  }
  tab_list[[i]] <- tab
}

names(tab_list) <- islands

tab_list <- lapply(tab_list, \(x) {
  x$model <- toupper(x$model)
  x$model <- gsub(pattern = "CR", replacement = "HR", x = x$model)
  x$model <- gsub(pattern = "_", replacement = " ", x = x$model)
  x$model <- gsub(pattern = "LAC", replacement = "$\\\\lambda^c$", x = x$model)
  x$model <- gsub(pattern = "MU", replacement = "$\\\\mu$", x = x$model)
  x$model <- gsub(pattern = "LAA", replacement = "$\\\\lambda^a$", x = x$model)
  x$model <- gsub(pattern = "RR K", replacement = "RR $K'$", x = x$model)
  x$model <- gsub(pattern = "RR ", replacement = "RR", x = x$model)
  x
})

tab_list <- lapply(tab_list, \(x) {
  colnames(x) <- c("Model", "$\\lambda^c$", "$\\mu$", "$K'$", "$\\gamma$",
                   "$\\lambda^a$", "$\\sigma$", "loglik", "BIC", "AIC", "AICc")
  x
})

# add archipelago names to data frame
tab_list <- Map(
  f = function(x, arch) {
    cbind(Archipelago = replicate(expr = arch, n = nrow(x)), x)
  },
  tab_list,
  names(tab_list)
)

tab_list <- lapply(tab_list, \(x) {
  x$Archipelago[1] <- paste0(
    "\\multirow{", nrow(x), "}{*}{", x$Archipelago[1], "}"
  )
  x$Archipelago[2:nrow(x)] <- rep("", length(2:nrow(x)))
  x
})

# remove DD from model name
tab_list <- lapply(tab_list, \(x) {
  x$Model <- gsub(pattern = " DD", replacement = "", x = x$Model)
  x
})

tab <- do.call(rbind, tab_list)
rownames(tab) <- NULL

print(
  xtable::xtable(
    tab,
    digits = c(
      3, # rownames
      3, # archipelago name
      3, # model name
      4, # cladogenesis
      4, # extincton
      5, # carrying capacity
      3, # colonisation
      4, # anagenesis
      6, # sd
      6, # loglik
      6, # BIC
      6, # AIC
      6  #AICc
    ),
    display = c("s", "s", rep("g", 11)),
    align = "ccccccccccccc",
    caption = paste0(
      "Maximum likelihood parameter estimates and model diagnostics for ",
      "each of the seven archipelagos. Results from fitting homogeneous-rate ",
      "(HR) and relaxed-rate (RR) models. Parameters estimated are: ",
      "cladogenesis ($\\lambda^c$), extinction ",
      "($\\mu$), carrying capacity (\\textit{K}), colonisation ($\\gamma$), ",
      "anagenesis ($\\lambda^a$), standard deviation of relaxed parameter ",
      "($\\sigma$), as well as the maximum log likelihood (loglik), and ",
      "Bayesian Information Criterion (BIC), Akaike Information Criterion ",
      "(AIC) and bias-corrected AIC (AICc). The homogeneous-rate model has ",
      "five free parameters (degrees of freedom), because $\\sigma$ is fixed ",
      "to zero, while the relaxed-rate models have six free parameters."),
    label = paste0("tab:archipelagos_ml")
  ),
  size = "footnotesize",
  math.style.exponents = TRUE,
  NA.string = "NA",
  include.rownames = FALSE,
  hline.after = c(0, seq(0, 35, 5)),
  sanitize.colnames.function = identity,
  sanitize.text.function = force
)