Discrete binning: effects
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.constructlibrary(gamlss.dist)#> Loading required package: MASSpath_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")
bwdata_file <- file.path(path_processed, "bwdata_51_test.fst")
model_file <- file.path(path_modelled, "bw-08-discrete-bin.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)Compute results
model$results <- results.bamlss.default(model)Fixed effects
summary(model)#>
#> Call:
#> bamlss(formula = form, data = bwdata_model, cores = 4, combine = FALSE,
#> light = TRUE, binning = TRUE, n.iter = 1000, burnin = 0)
#> ---
#> Family: gaussian
#> Link function: mu = identity, sigma = log
#> *---
#> Formula mu:
#> ---
#> born_weight ~ marital_status + race + s(municipality, bs = "re") +
#> s(age_bin)
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 2457.633 2444.368 2457.577 2470.664 2379.54
#> marital_status2 -13.169 -19.804 -13.091 -6.502 -12.96
#> marital_status3 -136.775 -143.568 -136.763 -129.845 -136.87
#> race2 14.589 6.632 14.690 22.598 14.90
#> race3 35.497 27.771 35.436 43.321 35.57
#> race4 57.397 49.406 57.283 65.472 57.85
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9951 0.9757 1.0000 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(municipality).tau21 4.045e+04 2.521e+04 3.926e+04 6.377e+04 3.759e+04
#> s(municipality).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
#> s(municipality).edf 3.992e+01 3.987e+01 3.992e+01 3.995e+01 3.991e+01
#> s(age_bin).tau21 1.224e+06 3.875e+05 9.913e+05 3.413e+06 6.099e+06
#> s(age_bin).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00 NA
#> s(age_bin).edf 8.374e+00 7.846e+00 8.393e+00 8.786e+00 8.874e+00
#> ---
#> Formula sigma:
#> ---
#> sigma ~ 1
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 5.295 5.286 5.296 5.305 5.297
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9974 0.9773 1.0000 1
#> ---
#> Sampler summary:
#> -
#> runtime = 37.748
#> ---
#> Optimizer summary:
#> -
#> AICc = 268762.6 edf = 55.7894 logLik = -134325.4
#> logPost = -134723.6 nobs = 20000 runtime = 5.81Smoothed effects
There seems to be a problem with the labels of the random effects plot.
par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot(model, scale = 0, scheme = 2, spar = FALSE)
Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 10.290 0.119 10.477