Growth model with t-distribution and no rainfall: fitting
Birthweight model including covariates at municipality level with linear effects, and controlling for gestational age.
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-25-growth-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 + gestation_weeks +
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 4441698. logPost -2383904 logLik -2220751 edf 97.387 eps 0.0654 iteration 1"
#> [2] "AICc 4413816. logPost -2209738 logLik -2206795 edf 112.20 eps 0.0133 iteration 2"
#> [3] "AICc 4411020. logPost -2206508 logLik -2205372 edf 137.01 eps 0.0036 iteration 3"
#> [4] "AICc 4409641. logPost -2205685 logLik -2204673 edf 147.00 eps 0.0019 iteration 4"
#> [5] "AICc 4408962. logPost -2205207 logLik -2204330 edf 150.99 eps 0.0017 iteration 5"
#> [6] "AICc 4408866. logPost -2205090 logLik -2204277 edf 155.51 eps 0.0007 iteration 6"
#> [7] "AICc 4408851. logPost -2205081 logLik -2204267 edf 158.06 eps 0.0004 iteration 7"
#> [8] "AICc 4408850. logPost -2205082 logLik -2204267 edf 158.16 eps 0.0001 iteration 8"
#> [9] "AICc 4408850. logPost -2205082 logLik -2204267 edf 158.15 eps 0.0001 iteration 9"
#> [10] "AICc 4408850. logPost -2205081 logLik -2204267 edf 158.15 eps 0.0001 iteration 10"
#> [11] "AICc 4408850. logPost -2205081 logLik -2204267 edf 158.15 eps 0.0001 iteration 10"
#> [12] "elapsed time: 4.77min"
#> [13] "Starting the sampler..."
#> [14] "Starting the sampler..."
#> [15] "Starting the sampler...Starting the sampler..."
#> [16] ""
#> [17] ""
#> [18] "| | 0% 871.23min"
#> [19] "| | 0% 876.26min"
#> [20] ""
#> [21] "| | 0% 899.38min| | 0% 899.38min"
#> [22] "|* | 5% 822.03min 43.26min"
#> [23] "|* | 5% 823.87min 43.36min"
#> [24] "|* | 5% 824.16min 43.38min"
#> [25] "|* | 5% 827.88min 43.57min"
#> [26] "|** | 10% 780.23min 86.69min"
#> [27] "|** | 10% 781.40min 86.82min"
#> [28] "|** | 10% 783.63min 87.07min"
#> [29] "|** | 10% 785.47min 87.27min"
#> [30] "|*** | 15% 736.27min 129.93min"
#> [31] "|*** | 15% 736.80min 130.02min"
#> [32] "|*** | 15% 738.25min 130.28min"
#> [33] "|*** | 15% 738.92min 130.40min"
#> [34] "|**** | 20% 692.82min 173.20min"
#> [35] "|**** | 20% 693.19min 173.30min"
#> [36] "|**** | 20% 693.55min 173.39min"
#> [37] "|**** | 20% 694.41min 173.60min"
#> [38] "|***** | 25% 651.29min 217.10min"
#> [39] "|***** | 25% 651.88min 217.29min"
#> [40] "|***** | 25% 651.92min 217.31min"
#> [41] "|***** | 25% 652.85min 217.62min"
#> [42] "|****** | 30% 608.23min 260.67min"
#> [43] "|****** | 30% 608.94min 260.97min"
#> [44] "|****** | 30% 609.28min 261.12min"
#> [45] "|****** | 30% 610.49min 261.64min"
#> [46] "|******* | 35% 564.91min 304.18min"
#> [47] "|******* | 35% 565.55min 304.53min"
#> [48] "|******* | 35% 565.92min 304.73min"
#> [49] "|******* | 35% 566.96min 305.28min"
#> [50] "|******** | 40% 521.40min 347.60min"
#> [51] "|******** | 40% 522.44min 348.29min"
#> [52] "|******** | 40% 522.44min 348.30min"
#> [53] "|******** | 40% 522.97min 348.64min"
#> [54] "|********* | 45% 477.93min 391.03min"
#> [55] "|********* | 45% 478.78min 391.73min"
#> [56] "|********* | 45% 479.10min 391.99min"
#> [57] "|********* | 45% 479.22min 392.09min"
#> [58] "|********** | 50% 434.71min 434.71min"
#> [59] "|********** | 50% 435.50min 435.50min"
#> [60] "|********** | 50% 435.58min 435.58min"
#> [61] "|********** | 50% 435.72min 435.72min"
#> [62] "|*********** | 55% 391.14min 478.06min"
#> [63] "|*********** | 55% 391.84min 478.92min"
#> [64] "|*********** | 55% 391.84min 478.91min"
#> [65] "|*********** | 55% 392.07min 479.20min"
#> [66] "|************ | 60% 347.72min 521.58min"
#> [67] "|************ | 60% 348.31min 522.46min"
#> [68] "|************ | 60% 348.39min 522.59min"
#> [69] "|************ | 60% 348.42min 522.62min"
#> [70] "|************* | 65% 304.29min 565.11min"
#> [71] "|************* | 65% 304.74min 565.95min"
#> [72] "|************* | 65% 304.81min 566.08min"
#> [73] "|************* | 65% 304.83min 566.11min"
#> [74] "|************** | 70% 260.81min 608.55min"
#> [75] "|************** | 70% 261.11min 609.25min"
#> [76] "|************** | 70% 261.22min 609.52min"
#> [77] "|************** | 70% 261.32min 609.75min"
#> [78] "|*************** | 75% 217.17min 651.50min"
#> [79] "|*************** | 75% 217.40min 652.19min"
#> [80] "|*************** | 75% 217.50min 652.50min"
#> [81] "|*************** | 75% 217.58min 652.73min"
#> [82] "|**************** | 80% 173.75min 695.00min"
#> [83] "|**************** | 80% 173.97min 695.88min"
#> [84] "|**************** | 80% 174.09min 696.36min"
#> [85] "|**************** | 80% 174.10min 696.38min"
#> [86] "|***************** | 85% 130.28min 738.28min"
#> [87] "|***************** | 85% 130.43min 739.11min"
#> [88] "|***************** | 85% 130.51min 739.58min"
#> [89] "|***************** | 85% 130.53min 739.66min"
#> [90] "|****************** | 90% 86.91min 782.15min"
#> [91] "|****************** | 90% 86.95min 782.52min"
#> [92] "|****************** | 90% 87.03min 783.28min"
#> [93] "|****************** | 90% 87.04min 783.36min"
#> [94] "|******************* | 95% 43.44min 825.41min"
#> [95] "|******************* | 95% 43.46min 825.82min"
#> [96] "|******************* | 95% 43.50min 826.50min"
#> [97] "|******************* | 95% 43.50min 826.52min"
#> [98] "|********************| 100% 0.00sec 868.59min"
#> [99] ""
#> [100] "|********************| 100% 0.00sec 868.83min"
#> [101] ""
#> [102] "|********************| 100% 0.00sec 869.39min"
#> [103] ""
#> [104] "|********************| 100% 0.00sec 869.48min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 9.556 0.032 9.600saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 209018.551 33.204 52534.708