Standardized linear model with random effects: 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")

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-muni-10-slm-re.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)

Summary

summary(model)

#> 
#> Call:
#> bamlss(formula = form, data = bwdata_model, cores = 4, combine = FALSE, 
#>     light = TRUE, n.iter = 1000, burnin = 0)
#> ---
#> Family: gaussian 
#> Link function: mu = identity, sigma = log
#> *---
#> Formula mu:
#> ---
#> born_weight ~ remoteness_cd + prop_tap_toilet_cd + s(res_muni, 
#>     bs = "re")
#> -
#> Parametric coefficients:
#>                        Mean     2.5%      50%    97.5% parameters
#> (Intercept)        3226.918 3223.352 3227.019 3230.216    3191.57
#> remoteness_cd       -60.545 -122.232  -59.422    5.945     -57.49
#> prop_tap_toilet_cd  124.857  -36.923  131.311  249.893     -51.00
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.9981 0.9913 1.0000     1
#> -
#> Smooth terms:
#>                      Mean    2.5%     50%   97.5% parameters
#> s(res_muni).tau21 2191.78 1347.33 2097.70 3478.19    3823.09
#> s(res_muni).alpha    1.00    1.00    1.00    1.00         NA
#> s(res_muni).edf     41.84   41.27   41.86   42.28      42.35
#> ---
#> Formula sigma:
#> ---
#> sigma ~ 1
#> -
#> Parametric coefficients:
#>              Mean  2.5%   50% 97.5% parameters
#> (Intercept) 6.225 6.223 6.225 6.228      6.225
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.9992 0.9940 1.0000     1
#> ---
#> Sampler summary:
#> -
#> runtime = 739.441
#> ---
#> Optimizer summary:
#> -
#> AICc = 4456381 edf = 46.3523 logLik = -2228144
#> logPost = -2228395 nobs = 291479 runtime = 9.96

Parametric effects

par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot2d_bamlss(model, bwdata_model, model = "mu", term = "remoteness_cd", grid = 50,
              FUN = c95)
plot2d_bamlss(model, bwdata_model, model = "mu", term = "prop_tap_toilet_cd", grid = 50,
              FUN = c95)

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  22.030   0.323  22.371