SPIFA Urban Ipixuna: Model comparison with DIC

Load required libraries and data

rm(list = ls())
library(day2day)
library(spifa)
library(xtable)
library(tibble)

path_main <- git_path()
path_data <- file.path(path_main, "data")
path_raw <- file.path(path_data, "raw")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")
path_fig <- file.path(path_data, "figures")

samples <- readRDS(file.path(path_modelled, "cifa-ipixuna-urban.rds"))
samples1 <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-1gp.rds"))
samples2 <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-2gp.rds"))
samples3 <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-3gp.rds"))
iter <- nrow(samples1[["theta"]])

DIC computation

dic_values <- lapply(
    list("CIFA" = samples, "SPIFA I" = samples1, "SPIFA II" = samples2, "SPIFA III" = samples3),
    function(x) dic(as_tibble(x, burnin = iter/5))
)

dic_df <- do.call(cbind, lapply(dic_values, function(x) as.numeric(x))) |>
    t() |> data.frame() |>
    setNames(c("Posterior Mean Deviance", "Effective Number of Parameters", "DIC"))
dic_df

#>           Posterior Mean Deviance Effective Number of Parameters      DIC
#> CIFA                     1894.307                       329.8208 2224.128
#> SPIFA I                  1866.667                       335.1415 2201.808
#> SPIFA II                 1861.217                       338.6363 2199.853
#> SPIFA III                1857.537                       341.2626 2198.799

Save Table

# convert to latex
dic_table <- xtable(
    dic_df,
    align = c("l", "c", "c", "c"),
    caption = paste0(
        "Deviance Information Criteria (DIC) for the Proposed Models: ",
        "without spatial correlation (CIFA), with spatial correlation in factor 1 (SPIFA I), ",
        "with spatial correlation in factor 1 and 3 (SPIFA II) and ",
        "with spatial correlation in all factors (SPIFA III)."),
    label = "tab:dics"
)

# improve row names and add lines
addtorow <- list()
addtorow$pos <- list(0, 0, 0, nrow(dic_table))
addtorow$command <- c(
    "\\toprule[.1em] \\multirow{2}[3]{*}{\\textbf{Model}} & \\multicolumn{3}{c}{\\textbf{Diagnostics}} \\\\\n",
    "\\cmidrule(rl){2-4} & Posterior Mean Deviance & Effective Number of Parameters & DIC \\\\\n",
    "\\midrule[.05em] ",
    "\\bottomrule[.1em] "
)

# print latex with additional formating
print(dic_table,
      file = file.path(path_fig, "results-dic-comparison.tex"),
      caption.placement = "top",
      size = "\\footnotesize",
      add.to.row = addtorow,
      include.colnames = FALSE,
      hline.after = NULL,
      )

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  26.208   4.169 117.756