Only covariates: fitting
Birthweight model including covariates at municipality level
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.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")
path_modelled_data <- file.path(path_modelled, "bw-muni-00-covs.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 ~ s(remoteness) + s(prop_tap_toilet)
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 4458377. logPost -2229229 logLik -2229184 edf 4.2021 eps 0.0299 iteration 1"
#> [2] "AICc 4458249. logPost -2229206 logLik -2229119 edf 5.3771 eps 0.0007 iteration 2"
#> [3] "AICc 4458121. logPost -2229167 logLik -2229052 edf 8.1029 eps 0.0007 iteration 3"
#> [4] "AICc 4457879. logPost -2229090 logLik -2228927 edf 12.443 eps 0.0007 iteration 4"
#> [5] "AICc 4457558. logPost -2228948 logLik -2228762 edf 17.218 eps 0.0013 iteration 5"
#> [6] "AICc 4457450. logPost -2228874 logLik -2228705 edf 19.531 eps 0.0008 iteration 6"
#> [7] "AICc 4457423. logPost -2228862 logLik -2228691 edf 19.797 eps 0.0003 iteration 7"
#> [8] "AICc 4457417. logPost -2228859 logLik -2228688 edf 19.770 eps 0.0001 iteration 8"
#> [9] "AICc 4457416. logPost -2228859 logLik -2228688 edf 19.757 eps 0.0000 iteration 9"
#> [10] "AICc 4457416. logPost -2228859 logLik -2228688 edf 19.757 eps 0.0000 iteration 9"
#> [11] "elapsed time: 22.54sec"
#> [12] "Starting the sampler..."
#> [13] "Starting the sampler..."
#> [14] "Starting the sampler...Starting the sampler..."
#> [15] ""
#> [16] ""
#> [17] "| | 0% 10.43min"
#> [18] "| | 0% 10.44min"
#> [19] "| | 0% 11.14min"
#> [20] "| | 0% 10.61min"
#> [21] "|* | 5% 10.11min 31.92sec"
#> [22] "|* | 5% 10.25min 32.35sec"
#> [23] "|* | 5% 10.17min 32.10sec"
#> [24] "|* | 5% 10.31min 32.54sec"
#> [25] ""
#> [26] "|** | 10% 9.54min 1.06min|** | 10% 9.59min 1.07min"
#> [27] "|** | 10% 9.65min 1.07min"
#> [28] "|** | 10% 9.70min 1.08min"
#> [29] "|*** | 15% 9.08min 1.60min"
#> [30] "|*** | 15% 9.04min 1.60min"
#> [31] "|*** | 15% 9.12min 1.61min"
#> [32] "|*** | 15% 9.19min 1.62min"
#> [33] "|**** | 20% 8.55min 2.14min"
#> [34] "|**** | 20% 8.55min 2.14min"
#> [35] "|**** | 20% 8.62min 2.15min"
#> [36] "|**** | 20% 8.64min 2.16min"
#> [37] "|***** | 25% 8.04min 2.68min"
#> [38] "|***** | 25% 8.04min 2.68min"
#> [39] "|***** | 25% 8.08min 2.69min"
#> [40] "|***** | 25% 8.10min 2.70min"
#> [41] "|****** | 30% 7.52min 3.22min"
#> [42] "|****** | 30% 7.53min 3.23min"
#> [43] "|****** | 30% 7.56min 3.24min"
#> [44] "|****** | 30% 7.57min 3.25min"
#> [45] "|******* | 35% 6.99min 3.76min"
#> [46] "|******* | 35% 7.00min 3.77min"
#> [47] "|******* | 35% 7.02min 3.78min"
#> [48] "|******* | 35% 7.04min 3.79min"
#> [49] "|******** | 40% 6.46min 4.31min"
#> [50] "|******** | 40% 6.47min 4.31min"
#> [51] "|******** | 40% 6.49min 4.32min"
#> [52] "|******** | 40% 6.50min 4.33min"
#> [53] "|********* | 45% 5.91min 4.84min"
#> [54] "|********* | 45% 5.92min 4.84min"
#> [55] "|********* | 45% 5.93min 4.85min"
#> [56] "|********* | 45% 5.94min 4.86min"
#> [57] "|********** | 50% 5.36min 5.36min"
#> [58] "|********** | 50% 5.37min 5.37min"
#> [59] "|********** | 50% 5.38min 5.38min"
#> [60] "|********** | 50% 5.39min 5.39min"
#> [61] "|*********** | 55% 4.82min 5.89min"
#> [62] "|*********** | 55% 4.82min 5.89min"
#> [63] "|*********** | 55% 4.83min 5.91min"
#> [64] "|*********** | 55% 4.84min 5.92min"
#> [65] "|************ | 60% 4.28min 6.41min"
#> [66] "|************ | 60% 4.28min 6.42min"
#> [67] "|************ | 60% 4.29min 6.43min"
#> [68] "|************ | 60% 4.29min 6.44min"
#> [69] "|************* | 65% 3.74min 6.94min"
#> [70] "|************* | 65% 3.74min 6.95min"
#> [71] "|************* | 65% 3.75min 6.96min"
#> [72] "|************* | 65% 3.75min 6.97min"
#> [73] "|************** | 70% 3.20min 7.46min"
#> [74] "|************** | 70% 3.20min 7.47min"
#> [75] "|************** | 70% 3.21min 7.48min"
#> [76] "|************** | 70% 3.21min 7.49min"
#> [77] "|*************** | 75% 2.66min 7.99min"
#> [78] "|*************** | 75% 2.67min 8.00min"
#> [79] "|*************** | 75% 2.67min 8.01min"
#> [80] "|*************** | 75% 2.67min 8.02min"
#> [81] "|**************** | 80% 2.13min 8.51min"
#> [82] "|**************** | 80% 2.13min 8.52min"
#> [83] "|**************** | 80% 2.13min 8.53min"
#> [84] "|**************** | 80% 2.13min 8.54min"
#> [85] "|***************** | 85% 1.60min 9.04min"
#> [86] "|***************** | 85% 1.60min 9.06min"
#> [87] "|***************** | 85% 1.60min 9.07min"
#> [88] "|***************** | 85% 1.60min 9.07min"
#> [89] "|****************** | 90% 1.06min 9.57min"
#> [90] "|****************** | 90% 1.07min 9.59min"
#> [91] "|****************** | 90% 1.07min 9.59min"
#> [92] "|****************** | 90% 1.07min 9.60min"
#> [93] "|******************* | 95% 31.86sec 10.09min"
#> [94] "|******************* | 95% 31.92sec 10.11min"
#> [95] "|******************* | 95% 31.95sec 10.12min"
#> [96] "|******************* | 95% 31.96sec 10.12min"
#> [97] "|********************| 100% 0.00sec 10.62min"
#> [98] ""
#> [99] "|********************| 100% 0.00sec 10.63min"
#> [100] ""
#> [101] "|********************| 100% 0.00sec 10.64min"
#> [102] ""
#> [103] "|********************| 100% 0.00sec 10.64min"
system.time(saveRDS(bamlss_model, file = path_modelled_data))
#> user system elapsed
#> 0.116 0.000 0.116
saveRDS(form, file = path_modelled_form)
Time to execute the task
Only useful when executed with Rscript
.
proc.time()
#> user system elapsed
#> 2626.596 2.224 684.595