Full model with t-distribution: 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-10-full-re-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_aux <- 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_mu <- update(form_mu_aux, . ~ . +
s(neg_mbsi_mean_1wk, pos_mbsi_mean_1wk, k = 70) +
s(neg_mckee_mean_8wk, pos_mckee_mean_8wk, k = 70) +
s(neg_ext_mbsi_mean_8wk, pos_ext_mbsi_mean_8wk, k = 70))
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 = 1000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4448681. logPost -2386947 logLik -2224240 edf 100.38 eps 0.0663 iteration 1"
#> [2] "AICc 4440900. logPost -2223494 logLik -2220332 edf 117.54 eps 0.0202 iteration 2"
#> [3] "AICc 4423473. logPost -2212970 logLik -2211597 edf 138.75 eps 0.0151 iteration 3"
#> [4] "AICc 4421401. logPost -2212141 logLik -2210533 edf 166.59 eps 0.0064 iteration 4"
#> [5] "AICc 4420583. logPost -2211746 logLik -2210087 edf 203.97 eps 0.0049 iteration 5"
#> [6] "AICc 4420401. logPost -2211676 logLik -2209959 edf 240.68 eps 0.0032 iteration 6"
#> [7] "AICc 4420322. logPost -2211685 logLik -2209896 edf 264.82 eps 0.0025 iteration 7"
#> [8] "AICc 4420307. logPost -2211683 logLik -2209886 edf 266.38 eps 0.0018 iteration 8"
#> [9] "AICc 4420301. logPost -2211678 logLik -2209884 edf 265.69 eps 0.0013 iteration 9"
#> [10] "AICc 4420297. logPost -2211673 logLik -2209883 edf 264.94 eps 0.0010 iteration 10"
#> [11] "AICc 4420295. logPost -2211670 logLik -2209883 edf 264.20 eps 0.0008 iteration 11"
#> [12] "AICc 4420294. logPost -2211666 logLik -2209883 edf 263.51 eps 0.0007 iteration 12"
#> [13] "AICc 4420293. logPost -2211663 logLik -2209883 edf 262.85 eps 0.0005 iteration 13"
#> [14] "AICc 4420293. logPost -2211660 logLik -2209884 edf 262.22 eps 0.0004 iteration 14"
#> [15] "AICc 4420293. logPost -2211658 logLik -2209884 edf 261.62 eps 0.0003 iteration 15"
#> [16] "AICc 4420293. logPost -2211655 logLik -2209885 edf 261.07 eps 0.0002 iteration 16"
#> [17] "AICc 4420292. logPost -2211653 logLik -2209885 edf 260.57 eps 0.0002 iteration 17"
#> [18] "AICc 4420293. logPost -2211651 logLik -2209886 edf 260.13 eps 0.0001 iteration 18"
#> [19] "AICc 4420293. logPost -2211649 logLik -2209886 edf 259.75 eps 0.0001 iteration 19"
#> [20] "AICc 4420293. logPost -2211647 logLik -2209886 edf 259.41 eps 0.0001 iteration 20"
#> [21] "AICc 4420293. logPost -2211647 logLik -2209886 edf 259.41 eps 0.0001 iteration 20"
#> [22] "elapsed time: 92.69min"
#> [23] "Starting the sampler..."
#> [24] "Starting the sampler..."
#> [25] "Starting the sampler..."
#> [26] "Starting the sampler..."
#> [27] ""
#> [28] "| | 0% 455.90min"
#> [29] "| | 0% 456.27min"
#> [30] "| | 0% 457.15min"
#> [31] "| | 0% 457.97min"
#> [32] "|* | 5% 435.45min 22.92min"
#> [33] "|* | 5% 435.41min 22.92min"
#> [34] "|* | 5% 436.16min 22.96min"
#> [35] "|* | 5% 437.41min 23.02min"
#> [36] "|** | 10% 412.22min 45.80min"
#> [37] "|** | 10% 413.00min 45.89min"
#> [38] "|** | 10% 413.27min 45.92min"
#> [39] "|** | 10% 414.97min 46.11min"
#> [40] "|*** | 15% 389.19min 68.68min"
#> [41] "|*** | 15% 389.97min 68.82min"
#> [42] "|*** | 15% 390.04min 68.83min"
#> [43] "|*** | 15% 391.59min 69.10min"
#> [44] "|**** | 20% 366.23min 91.56min"
#> [45] "|**** | 20% 366.94min 91.73min"
#> [46] "|**** | 20% 367.01min 91.75min"
#> [47] "|**** | 20% 368.44min 92.11min"
#> [48] "|***** | 25% 343.28min 114.43min"
#> [49] "|***** | 25% 343.93min 114.64min"
#> [50] "|***** | 25% 344.05min 114.68min"
#> [51] "|***** | 25% 345.31min 115.10min"
#> [52] "|****** | 30% 320.36min 137.30min"
#> [53] "|****** | 30% 320.92min 137.54min"
#> [54] "|****** | 30% 321.10min 137.61min"
#> [55] "|****** | 30% 322.25min 138.11min"
#> [56] "|******* | 35% 297.47min 160.18min"
#> [57] "|******* | 35% 297.95min 160.43min"
#> [58] "|******* | 35% 298.17min 160.56min"
#> [59] "|******* | 35% 299.19min 161.10min"
#> [60] "|******** | 40% 274.64min 183.09min"
#> [61] "|******** | 40% 275.12min 183.41min"
#> [62] "|******** | 40% 275.24min 183.50min"
#> [63] "|******** | 40% 276.04min 184.02min"
#> [64] "|********* | 45% 251.76min 205.99min"
#> [65] "|********* | 45% 252.25min 206.39min"
#> [66] "|********* | 45% 252.30min 206.43min"
#> [67] "|********* | 45% 252.93min 206.94min"
#> [68] "|********** | 50% 228.88min 228.88min"
#> [69] "|********** | 50% 229.36min 229.36min"
#> [70] "|********** | 50% 229.37min 229.37min"
#> [71] "|********** | 50% 229.87min 229.87min"
#> [72] "|*********** | 55% 206.00min 251.77min"
#> [73] "|*********** | 55% 206.44min 252.31min"
#> [74] "|*********** | 55% 206.47min 252.35min"
#> [75] "|*********** | 55% 206.83min 252.80min"
#> [76] "|************ | 60% 183.10min 274.65min"
#> [77] "|************ | 60% 183.50min 275.25min"
#> [78] "|************ | 60% 183.53min 275.30min"
#> [79] "|************ | 60% 183.82min 275.73min"
#> [80] "|************* | 65% 160.20min 297.52min"
#> [81] "|************* | 65% 160.56min 298.18min"
#> [82] "|************* | 65% 160.61min 298.27min"
#> [83] "|************* | 65% 160.82min 298.67min"
#> [84] "|************** | 70% 137.31min 320.39min"
#> [85] "|************** | 70% 137.61min 321.10min"
#> [86] "|************** | 70% 137.66min 321.20min"
#> [87] "|************** | 70% 137.84min 321.62min"
#> [88] "|*************** | 75% 114.42min 343.27min"
#> [89] "|*************** | 75% 114.67min 344.02min"
#> [90] "|*************** | 75% 114.70min 344.10min"
#> [91] "|*************** | 75% 114.87min 344.62min"
#> [92] "|**************** | 80% 91.55min 366.21min"
#> [93] "|**************** | 80% 91.72min 366.90min"
#> [94] "|**************** | 80% 91.75min 367.01min"
#> [95] "|**************** | 80% 91.90min 367.60min"
#> [96] "|***************** | 85% 68.68min 389.17min"
#> [97] "|***************** | 85% 68.79min 389.79min"
#> [98] "|***************** | 85% 68.81min 389.91min"
#> [99] "|***************** | 85% 68.92min 390.56min"
#> [100] "|****************** | 90% 45.79min 412.12min"
#> [101] "|****************** | 90% 45.85min 412.67min"
#> [102] "|****************** | 90% 45.87min 412.85min"
#> [103] "|****************** | 90% 45.94min 413.48min"
#> [104] "|******************* | 95% 22.90min 435.08min"
#> [105] "|******************* | 95% 22.93min 435.61min"
#> [106] "|******************* | 95% 22.93min 435.76min"
#> [107] "|******************* | 95% 22.97min 436.40min"
#> [108] "|********************| 100% 0.00sec 458.01min"
#> [109] ""
#> [110] "|********************| 100% 0.00sec 458.54min"
#> [111] ""
#> [112] "|********************| 100% 0.00sec 458.66min"
#> [113] ""
#> [114] "|********************| 100% 0.00sec 459.35min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 39.288 0.088 39.380saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 116556.549 52.634 33672.227