Full model with t-distribution: dic


Load packages, read data and source custom scripts

rm(list = ls())
library(bamlss)
#> Loading required package: coda
#> Loading required package: colorspace
#> Loading required package: mgcv
#> Loading required package: nlme
#> This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
#> 
#> Attaching package: 'bamlss'
#> The following object is masked from 'package:mgcv':
#> 
#>     smooth.construct
library(gamlss.dist)
#> Loading required package: MASS
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(ggplot2)

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")

source(file.path(path_proj, "src", "51-bamlss.R"))

bwdata_file <- file.path(path_processed, "bwdata_41_model.fst")
model_file <- file.path(path_modelled, "bw-10-full-re-t.rds")
form_file <- gsub("(\\.rds)$", "-form\\1", model_file)
model_file_burned <- gsub("(\\.rds)$", "-burned\\1", model_file)

bwdata_model <- fst::read_fst(bwdata_file)
form <- readRDS(form_file)
model <- readRDS(model_file_burned)

Diagnostics

frame <- bamlss.frame(form, family = TF, data = bwdata_model, model.matrix = FALSE,
                      smooth.construct = FALSE)
model$y <- frame$y

par(mar = c(4, 4, 3, 1), mfrow = c(1, 2))
plot(model, which = "qq-resid", spar = FALSE)
plot(model, which = "hist-resid", spar = FALSE)

DIC

system.time(dic <- DIC(model))
#>     user   system  elapsed 
#> 1563.380  219.187 1782.851
dic
#>      DIC       pd
#>  4420308 234.9712

Time to execute the task

Only useful when executed with Rscript.

proc.time()
#>     user   system  elapsed 
#> 3338.252  540.038 3879.996