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.construct

library(gamlss.dist)

#> Loading required package: MASS

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")

bwdata_file <- file.path(path_processed, "bwdata_51_test.fst")
model_file <- file.path(path_modelled, "bw-04-inter-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)

# plot(bwdata_model$age, bwdata_model$age_bin)
# mean(bwdata_model$age)
# mean(bwdata_model$age_bin)
# length(unique(bwdata_model$age))
# length(unique(bwdata_model$age_bin))

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)
#> -
#> Parametric coefficients:
#>                     Mean     2.5%      50%    97.5% parameters
#> (Intercept)     2828.838 2822.160 2828.826 2835.651    2796.59
#> marital_status2  -13.163  -19.911  -13.171   -6.572     -13.21
#> marital_status3 -136.723 -143.201 -136.741 -130.116    -136.83
#> race2             14.944    7.095   14.877   23.087      14.93
#> race3             35.840   28.093   35.811   43.724      35.87
#> race4             57.353   49.512   57.326   65.188      57.43
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.9945 0.9738 0.9999     1
#> -
#> Smooth terms:
#>                            Mean      2.5%       50%     97.5% parameters
#> s(municipality).tau21 4.039e+04 2.550e+04 3.885e+04 6.306e+04  3.846e+04
#> s(municipality).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
#> s(municipality).edf   3.992e+01 3.988e+01 3.992e+01 3.995e+01  3.992e+01
#> s(age).tau21          7.308e+05 2.337e+05 5.902e+05 2.081e+06  2.287e+05
#> s(age).alpha          1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
#> s(age).edf            8.607e+00 8.195e+00 8.631e+00 8.888e+00  8.181e+00
#> ---
#> Formula sigma:
#> ---
#> sigma ~ 1
#> -
#> Parametric coefficients:
#>              Mean  2.5%   50% 97.5% parameters
#> (Intercept) 5.295 5.285 5.295 5.304      5.293
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.9972 0.9774 1.0000     1
#> ---
#> Sampler summary:
#> -
#> runtime = 40.429
#> ---
#> Optimizer summary:
#> -
#> AICc = 268602.2 edf = 55.0982 logLik = -134245.9
#> logPost = -134635.2 nobs = 20000 runtime = 1.443

Smoothed 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.975   0.116  11.113