Full model with t-distribution and no rainfall: fitting
Birthweight model including covariates at municipality level with linear 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")
path_modelled_data <- file.path(path_modelled, "bw-15-full-no-rain-t.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
form_mu <- born_weight ~ sex + born_race + birth_place +
marital_status + study_years + consult_num + s(age) +
s(wk_ini) + s(rivwk_conception, bs = "cc") +
remoteness + prop_tap_toilet + s(res_muni, bs = "re")
form_sigma <- sigma ~ sex + born_race + birth_place +
study_years + consult_num + s(age) +
s(wk_ini) + 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, family = TF, data = bwdata_model,
n.iter = 8000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4448776. logPost -2387552 logLik -2224293 edf 94.412 eps 0.0662 iteration 1"
#> [2] "AICc 4440890. logPost -2223450 logLik -2220334 edf 111.21 eps 0.0202 iteration 2"
#> [3] "AICc 4423476. logPost -2212706 logLik -2211608 edf 129.51 eps 0.0151 iteration 3"
#> [4] "AICc 4421423. logPost -2211649 logLik -2210566 edf 144.99 eps 0.0065 iteration 4"
#> [5] "AICc 4420664. logPost -2211060 logLik -2210183 edf 148.40 eps 0.0049 iteration 5"
#> [6] "AICc 4420573. logPost -2210920 logLik -2210133 edf 152.79 eps 0.0031 iteration 6"
#> [7] "AICc 4420557. logPost -2210911 logLik -2210123 edf 155.18 eps 0.0023 iteration 7"
#> [8] "AICc 4420551. logPost -2210908 logLik -2210120 edf 155.17 eps 0.0017 iteration 8"
#> [9] "AICc 4420547. logPost -2210906 logLik -2210118 edf 155.16 eps 0.0013 iteration 9"
#> [10] "AICc 4420544. logPost -2210905 logLik -2210117 edf 155.14 eps 0.0011 iteration 10"
#> [11] "AICc 4420543. logPost -2210904 logLik -2210116 edf 155.13 eps 0.0008 iteration 11"
#> [12] "AICc 4420542. logPost -2210903 logLik -2210115 edf 155.12 eps 0.0007 iteration 12"
#> [13] "AICc 4420541. logPost -2210902 logLik -2210115 edf 155.11 eps 0.0005 iteration 13"
#> [14] "AICc 4420541. logPost -2210902 logLik -2210115 edf 155.10 eps 0.0004 iteration 14"
#> [15] "AICc 4420540. logPost -2210902 logLik -2210115 edf 155.09 eps 0.0003 iteration 15"
#> [16] "AICc 4420540. logPost -2210902 logLik -2210115 edf 155.08 eps 0.0002 iteration 16"
#> [17] "AICc 4420540. logPost -2210901 logLik -2210115 edf 155.08 eps 0.0002 iteration 17"
#> [18] "AICc 4420540. logPost -2210901 logLik -2210115 edf 155.07 eps 0.0001 iteration 18"
#> [19] "AICc 4420540. logPost -2210901 logLik -2210115 edf 155.06 eps 0.0001 iteration 19"
#> [20] "AICc 4420540. logPost -2210901 logLik -2210115 edf 155.06 eps 0.0001 iteration 20"
#> [21] "AICc 4420540. logPost -2210901 logLik -2210115 edf 155.06 eps 0.0001 iteration 20"
#> [22] "elapsed time: 9.83min"
#> [23] "Starting the sampler..."
#> [24] "Starting the sampler...Starting the sampler..."
#> [25] ""
#> [26] "Starting the sampler..."
#> [27] ""
#> [28] "| | 0% 865.28min"
#> [29] "| | 0% 873.81min"
#> [30] "| | 0% 910.33min"
#> [31] "| | 0% 910.47min"
#> [32] "|* | 5% 818.03min 43.05min"
#> [33] "|* | 5% 819.78min 43.15min"
#> [34] "|* | 5% 819.85min 43.15min"
#> [35] "|* | 5% 821.04min 43.21min"
#> [36] "|** | 10% 778.44min 86.49min"
#> [37] "|** | 10% 780.25min 86.69min"
#> [38] "|** | 10% 780.50min 86.72min"
#> [39] "|** | 10% 781.17min 86.80min"
#> [40] "|*** | 15% 734.58min 129.63min"
#> [41] "|*** | 15% 735.61min 129.81min"
#> [42] "|*** | 15% 735.62min 129.81min"
#> [43] "|*** | 15% 736.09min 129.90min"
#> [44] "|**** | 20% 691.55min 172.89min"
#> [45] "|**** | 20% 693.05min 173.26min"
#> [46] "|**** | 20% 693.09min 173.27min"
#> [47] "|**** | 20% 693.59min 173.40min"
#> [48] "|***** | 25% 649.76min 216.59min"
#> [49] "|***** | 25% 650.55min 216.85min"
#> [50] "|***** | 25% 650.73min 216.91min"
#> [51] "|***** | 25% 651.15min 217.05min"
#> [52] "|****** | 30% 606.35min 259.87min"
#> [53] "|****** | 30% 606.75min 260.03min"
#> [54] "|****** | 30% 606.97min 260.13min"
#> [55] "|****** | 30% 607.65min 260.42min"
#> [56] "|******* | 35% 562.54min 302.91min"
#> [57] "|******* | 35% 562.87min 303.08min"
#> [58] "|******* | 35% 563.08min 303.20min"
#> [59] "|******* | 35% 563.57min 303.46min"
#> [60] "|******** | 40% 519.25min 346.17min"
#> [61] "|******** | 40% 519.64min 346.43min"
#> [62] "|******** | 40% 519.75min 346.50min"
#> [63] "|******** | 40% 520.54min 347.03min"
#> [64] "|********* | 45% 476.06min 389.50min"
#> [65] "|********* | 45% 476.36min 389.75min"
#> [66] "|********* | 45% 476.49min 389.86min"
#> [67] "|********* | 45% 477.13min 390.38min"
#> [68] "|********** | 50% 432.67min 432.67min"
#> [69] "|********** | 50% 432.91min 432.91min"
#> [70] "|********** | 50% 433.06min 433.06min"
#> [71] "|********** | 50% 433.57min 433.57min"
#> [72] "|*********** | 55% 389.53min 476.09min"
#> [73] "|*********** | 55% 389.73min 476.34min"
#> [74] "|*********** | 55% 390.00min 476.67min"
#> [75] "|*********** | 55% 390.30min 477.04min"
#> [76] "|************ | 60% 346.31min 519.47min"
#> [77] "|************ | 60% 346.78min 520.18min"
#> [78] "|************ | 60% 346.82min 520.23min"
#> [79] "|************ | 60% 347.05min 520.57min"
#> [80] "|************* | 65% 302.95min 562.63min"
#> [81] "|************* | 65% 303.36min 563.38min"
#> [82] "|************* | 65% 303.41min 563.47min"
#> [83] "|************* | 65% 303.59min 563.81min"
#> [84] "|************** | 70% 259.80min 606.19min"
#> [85] "|************** | 70% 260.12min 606.94min"
#> [86] "|************** | 70% 260.16min 607.05min"
#> [87] "|************** | 70% 260.26min 607.28min"
#> [88] "|*************** | 75% 216.76min 650.29min"
#> [89] "|*************** | 75% 216.96min 650.89min"
#> [90] "|*************** | 75% 216.99min 650.98min"
#> [91] "|*************** | 75% 217.05min 651.14min"
#> [92] "|**************** | 80% 173.43min 693.71min"
#> [93] "|**************** | 80% 173.56min 694.25min"
#> [94] "|**************** | 80% 173.58min 694.33min"
#> [95] "|**************** | 80% 173.65min 694.59min"
#> [96] "|***************** | 85% 130.08min 737.11min"
#> [97] "|***************** | 85% 130.15min 737.52min"
#> [98] "|***************** | 85% 130.19min 737.75min"
#> [99] "|***************** | 85% 130.27min 738.22min"
#> [100] "|****************** | 90% 86.68min 780.08min"
#> [101] "|****************** | 90% 86.72min 780.47min"
#> [102] "|****************** | 90% 86.75min 780.71min"
#> [103] "|****************** | 90% 86.80min 781.21min"
#> [104] "|******************* | 95% 43.35min 823.60min"
#> [105] "|******************* | 95% 43.36min 823.93min"
#> [106] "|******************* | 95% 43.39min 824.34min"
#> [107] "|******************* | 95% 43.41min 824.77min"
#> [108] "|********************| 100% 0.00sec 866.70min"
#> [109] ""
#> [110] "|********************| 100% 0.00sec 867.02min"
#> [111] ""
#> [112] "|********************| 100% 0.00sec 867.27min"
#> [113] ""
#> [114] "|********************| 100% 0.00sec 867.63min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 9.272 0.032 9.305saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 208859.318 40.262 52728.021