Standardized linear model with random effects: fitting
Birthweight model including covariates at municipality level with linear effects.
Load packages, read data and source custom scripts
Paths are defined relative to the git repository location.
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")
path_modelled_data <- file.path(path_modelled, "bw-muni-10-slm-re.rds")
path_modelled_sink <- gsub("\\.rds$", "\\.txt", path_modelled_data)
path_modelled_form <- gsub("(\\.rds)$", "-form\\1", path_modelled_data)
bwdata_model <- fst::read_fst(file.path(path_processed, "bwdata_41_model.fst"))Define formula for our model
Now we define the same models as in the previous study.
form_sigma <- sigma ~ 1
form_mu <- born_weight ~ remoteness_cd + prop_tap_toilet_cd + s(res_muni, bs = "re")
form <- list(form_mu, form_sigma)Run the model of interest and save results
{
sink(path_modelled_sink)
bamlss_model <- bamlss(
form, data = bwdata_model,
n.iter = 1000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4465334. logPost -2454953 logLik -2232662 edf 5.0881 eps 0.0321 iteration 1"
#> [2] "AICc 4458447. logPost -2231450 logLik -2229211 edf 12.679 eps 0.0074 iteration 2"
#> [3] "AICc 4456610. logPost -2228817 logLik -2228272 edf 33.194 eps 0.0038 iteration 3"
#> [4] "AICc 4456383. logPost -2228414 logLik -2228146 edf 44.877 eps 0.0015 iteration 4"
#> [5] "AICc 4456381. logPost -2228395 logLik -2228144 edf 46.357 eps 0.0001 iteration 5"
#> [6] "AICc 4456381. logPost -2228395 logLik -2228144 edf 46.352 eps 0.0000 iteration 6"
#> [7] "AICc 4456381. logPost -2228395 logLik -2228144 edf 46.352 eps 0.0000 iteration 6"
#> [8] "elapsed time: 9.96sec"
#> [9] "Starting the sampler..."
#> [10] "Starting the sampler...Starting the sampler..."
#> [11] ""
#> [12] "Starting the sampler..."
#> [13] ""
#> [14] "| "
#> [15] " "
#> [16] " "
#> [17] " | 0% 12.70min| | | 0% 12.79min | 0% 12.79min| | 0% 12.79min"
#> [18] ""
#> [19] ""
#> [20] ""
#> [21] "|* | 5%|* | 5% 11.52min11.52min 36.37sec36.37sec|* | 5%|* | 5% 11.52min11.52min 36.37sec36.37sec"
#> [22] "|** | 10% 10.88min 1.21min"
#> [23] "|** | 10% 10.88min 1.21min"
#> [24] "|** | 10% 10.89min 1.21min"
#> [25] "|** | 10% 10.91min 1.21min"
#> [26] ""
#> [27] "|*** | 15% 10.22min 1.80min|*** | 15% 10.22min 1.80min"
#> [28] "|*** | 15% 10.23min 1.80min"
#> [29] "|*** | 15% 10.28min 1.81min"
#> [30] ""
#> [31] "|**** | 20%|**** | 20% 9.61min 9.61min 2.40min 2.40min"
#> [32] "|**** | 20% 9.62min 2.41min"
#> [33] "|**** | 20% 9.66min 2.41min"
#> [34] "|***** | 25% 9.03min 3.01min"
#> [35] "|***** | 25% 9.03min 3.01min"
#> [36] "|***** | 25% 9.05min 3.02min"
#> [37] "|***** | 25% 9.07min 3.02min"
#> [38] "|****** | 30% 8.45min 3.62min"
#> [39] "|****** | 30% 8.45min 3.62min"
#> [40] "|****** | 30% 8.47min 3.63min"
#> [41] "|****** | 30% 8.48min 3.63min"
#> [42] ""
#> [43] "|******* | 35%|******* | 35% 7.84min 7.84min 4.22min 4.22min"
#> [44] "|******* | 35% 7.87min 4.24min"
#> [45] "|******* | 35% 7.88min 4.24min"
#> [46] "|******** | 40% 7.25min 4.83min"
#> [47] "|******** | 40% 7.25min 4.83min"
#> [48] "|******** | 40% 7.26min 4.84min"
#> [49] "|******** | 40% 7.28min 4.85min"
#> [50] "|********* | 45% 6.63min 5.43min"
#> [51] "|********* | 45% 6.63min 5.43min"
#> [52] "|********* | 45% 6.65min 5.44min"
#> [53] "|********* | 45% 6.67min 5.46min"
#> [54] "|********** | 50% 6.03min 6.03min"
#> [55] "|********** | 50% 6.03min 6.03min"
#> [56] "|********** | 50% 6.04min 6.04min"
#> [57] "|********** | 50% 6.07min 6.07min"
#> [58] ""
#> [59] "|*********** | 55% |*********** | 55% 5.43min 5.43min 6.63min 6.63min"
#> [60] "|*********** | 55% 5.44min 6.64min"
#> [61] "|*********** | 55% 5.46min 6.67min"
#> [62] ""
#> [63] "|************ | 60%|************ | 60% 4.82min 4.82min 7.23min 7.23min"
#> [64] "|************ | 60% 4.83min 7.25min"
#> [65] "|************ | 60% 4.85min 7.28min"
#> [66] "|************* | 65% 4.22min 7.84min"
#> [67] "|************* | 65% 4.22min 7.84min"
#> [68] "|************* | 65% 4.23min 7.85min"
#> [69] "|************* | 65% 4.25min 7.89min"
#> [70] "|************** | 70% 3.62min 8.45min"
#> [71] "|************** | 70% 3.62min 8.45min"
#> [72] "|************** | 70% 3.63min 8.47min"
#> [73] "|************** | 70% 3.64min 8.50min"
#> [74] "|*************** | 75% 3.02min 9.06min"
#> [75] "|*************** | 75% 3.02min 9.06min"
#> [76] "|*************** | 75% 3.03min 9.08min"
#> [77] "|*************** | 75% 3.04min 9.11min"
#> [78] "|**************** | 80% 2.42min 9.67min"
#> [79] "|**************** | 80% 2.42min 9.67min"
#> [80] "|**************** | 80% 2.42min 9.68min"
#> [81] "|**************** | 80% 2.43min 9.72min"
#> [82] "|***************** | 85% 1.81min 10.28min"
#> [83] "|***************** | 85% 1.81min 10.28min"
#> [84] "|***************** | 85% 1.82min 10.29min"
#> [85] "|***************** | 85% 1.82min 10.33min"
#> [86] "|****************** | 90% 1.21min 10.88min"
#> [87] "|****************** | 90% 1.21min 10.89min"
#> [88] "|****************** | 90% 1.21min 10.90min"
#> [89] "|****************** | 90% 1.22min 10.94min"
#> [90] "|******************* | 95% 36.31sec 11.50min"
#> [91] "|******************* | 95% 36.33sec 11.50min"
#> [92] "|******************* | 95% 36.34sec 11.51min"
#> [93] "|******************* | 95% 36.50sec 11.56min"
#> [94] "|********************| 100% 0.00sec 12.10min"
#> [95] ""
#> [96] "|********************| 100% 0.00sec 12.11min"
#> [97] ""
#> [98] "|********************| 100% 0.00sec 12.11min"
#> [99] ""
#> [100] "|********************| 100% 0.00sec 12.14min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 2.336 0.000 2.340saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 2965.231 8.623 763.231