Growth model with t-distribution: 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-20-growth-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 + 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_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 = 30, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4441542. logPost -2383263 logLik -2220668 edf 103.35 eps 0.0655 iteration 1"
#> [2] "AICc 4413785. logPost -2209756 logLik -2206773 edf 118.61 eps 0.0133 iteration 2"
#> [3] "AICc 4411011. logPost -2206609 logLik -2205360 edf 144.98 eps 0.0036 iteration 3"
#> [4] "AICc 4409619. logPost -2206041 logLik -2204646 edf 163.55 eps 0.0019 iteration 4"
#> [5] "AICc 4408925. logPost -2205778 logLik -2204269 edf 193.24 eps 0.0017 iteration 5"
#> [6] "AICc 4408809. logPost -2205768 logLik -2204188 edf 216.23 eps 0.0007 iteration 6"
#> [7] "AICc 4408766. logPost -2205814 logLik -2204152 edf 230.33 eps 0.0005 iteration 7"
#> [8] "AICc 4408744. logPost -2205864 logLik -2204120 edf 251.50 eps 0.0004 iteration 8"
#> [9] "AICc 4408739. logPost -2205867 logLik -2204115 edf 253.51 eps 0.0002 iteration 9"
#> [10] "AICc 4408738. logPost -2205868 logLik -2204115 edf 253.65 eps 0.0001 iteration 10"
#> [11] "AICc 4408737. logPost -2205867 logLik -2204115 edf 253.57 eps 0.0001 iteration 11"
#> [12] "AICc 4408737. logPost -2205867 logLik -2204115 edf 253.57 eps 0.0001 iteration 11"
#> [13] "elapsed time: 14.11min"
#> [14] "Starting the sampler..."
#> [15] "Starting the sampler..."
#> [16] "Starting the sampler..."
#> [17] "Starting the sampler..."
#> [18] ""
#> [19] ""
#> [20] ""
#> [21] ""
#> [22] "|* | 3% 5.38min 11.14sec|* | 3% 5.39min 11.14sec|* | 3%|* | 3% 5.38min 5.39min 11.14sec11.16sec"
#> [23] "|* | 7% 5.14min 22.02sec"
#> [24] "|* | 7% 5.14min 22.04sec"
#> [25] "|* | 7% 5.14min 22.03sec"
#> [26] "|* | 7% 5.14min 22.04sec"
#> [27] "|** | 10% 4.91min 32.76sec"
#> [28] "|** | 10% 4.93min 32.84sec"
#> [29] ""
#> [30] "|** | 10%|** | 10% 4.93min 4.93min 32.87sec32.86sec"
#> [31] "|*** | 13% 4.60min 42.44sec"
#> [32] "|*** | 13% 4.67min 43.15sec"
#> [33] ""
#> [34] "|*** | 13% |*** | 13% 4.68min 4.68min43.16sec 43.18sec"
#> [35] "|*** | 17% 4.35min 52.25sec"
#> [36] "|*** | 17% 4.45min 53.35sec"
#> [37] ""
#> [38] "|*** | 17% 4.45min 53.37sec|*** | 17% 4.45min 53.35sec"
#> [39] "|**** | 20% 4.13min 1.03min"
#> [40] "|**** | 20% 4.25min 1.06min"
#> [41] ""
#> [42] "|**** | 20% 4.25min 1.06min|**** | 20% 4.25min 1.06min"
#> [43] "|***** | 23% 3.92min 1.19min"
#> [44] "|***** | 23% 4.05min 1.23min"
#> [45] "|***** | 23% 4.05min 1.23min"
#> [46] "|***** | 23% 4.05min 1.23min"
#> [47] "|***** | 27% 3.73min 1.36min"
#> [48] "|***** | 27% 3.84min 1.40min"
#> [49] ""
#> [50] "|***** | 27% |***** | 27% 3.85min 3.85min 1.40min 1.40min"
#> [51] "|****** | 30% 3.54min 1.52min"
#> [52] "|****** | 30% 3.63min 1.56min"
#> [53] "|****** | 30% 3.66min 1.57min"
#> [54] "|****** | 30% 3.66min 1.57min"
#> [55] "| | 0% 3.35min"
#> [56] "|******* | 33% 3.35min 1.68min"
#> [57] "| | 0% 3.43min"
#> [58] "|******* | 33% 3.43min 1.71min"
#> [59] "| | 0% 3.47min"
#> [60] "|******* | 33% 3.47min 1.73min"
#> [61] "| | 0% 3.47min"
#> [62] "|******* | 33% 3.47min 1.73min"
#> [63] "|******* | 37% 3.17min 1.84min"
#> [64] "|******* | 37% 3.24min 1.87min"
#> [65] "|******* | 37% 3.28min 1.90min"
#> [66] "|******* | 37% 3.28min 1.90min"
#> [67] "|******** | 40% 3.00min 2.00min"
#> [68] "|******** | 40% 3.05min 2.04min"
#> [69] "|******** | 40% 3.10min 2.07min"
#> [70] "|******** | 40% 3.11min 2.07min"
#> [71] "|********* | 43% 2.82min 2.16min"
#> [72] "|********* | 43% 2.87min 2.20min"
#> [73] "|********* | 43% 2.92min 2.23min"
#> [74] "|********* | 43% 2.92min 2.23min"
#> [75] "|********* | 47% 2.65min 2.32min"
#> [76] "|********* | 47% 2.70min 2.36min"
#> [77] "|********* | 47% 2.74min 2.40min"
#> [78] "|********* | 47% 2.74min 2.40min"
#> [79] "|********** | 50% 2.49min 2.49min"
#> [80] "|********** | 50% 2.53min 2.53min"
#> [81] ""
#> [82] "|********** | 50% 2.58min 2.58min|********** | 50% 2.59min 2.59min"
#> [83] "|*********** | 53% 2.44min 2.79min"
#> [84] ""
#> [85] ""
#> [86] "|*********** | 53%|*********** | 53%|*********** | 53% 2.44min 2.45min 2.45min 2.79min 2.79min 2.79min"
#> [87] "|*********** | 57% 2.27min 2.97min"
#> [88] "|*********** | 57% 2.27min 2.97min"
#> [89] "|*********** | 57% 2.27min 2.97min"
#> [90] "|*********** | 57% 2.27min 2.97min"
#> [91] "|************ | 60% 2.09min 3.14min"
#> [92] "|************ | 60% 2.09min 3.14min"
#> [93] "|************ | 60% 2.09min 3.14min"
#> [94] "|************ | 60% 2.09min 3.14min"
#> [95] ""
#> [96] "|************* | 63%|************* | 63% 1.91min 1.91min 3.30min 3.30min"
#> [97] "|************* | 63% 1.91min 3.31min"
#> [98] "|************* | 63% 1.91min 3.31min"
#> [99] "|************* | 67% 1.73min 3.47min"
#> [100] "|************* | 67% 1.73min 3.47min"
#> [101] "|************* | 67% 1.74min 3.47min"
#> [102] "|************* | 67% 1.74min 3.47min"
#> [103] "|************** | 70% 1.56min 3.63min"
#> [104] "|************** | 70% 1.56min 3.63min"
#> [105] "|************** | 70% 1.56min 3.64min"
#> [106] "|************** | 70% 1.56min 3.64min"
#> [107] "|*************** | 73% 1.38min 3.80min"
#> [108] "|*************** | 73% 1.38min 3.80min"
#> [109] "|*************** | 73% 1.38min 3.81min"
#> [110] "|*************** | 73% 1.38min 3.81min"
#> [111] "|*************** | 77% 1.21min 3.97min"
#> [112] "|*************** | 77% 1.21min 3.97min"
#> [113] "|*************** | 77% 1.21min 3.97min"
#> [114] "|*************** | 77% 1.21min 3.97min"
#> [115] "|**************** | 80% 1.03min 4.13min"
#> [116] "|**************** | 80% 1.03min 4.13min"
#> [117] "|**************** | 80% 1.04min 4.14min"
#> [118] "|**************** | 80% 1.04min 4.14min"
#> [119] ""
#> [120] "|***************** | 83% 51.64sec 4.30min|***************** | 83% 51.64sec 4.30min"
#> [121] "|***************** | 83% 51.75sec "
#> [122] " 4.31min|***************** | 83% 51.76sec 4.31min"
#> [123] "|***************** | 87% 41.25sec 4.47min"
#> [124] "|***************** | 87% 41.25sec 4.47min"
#> [125] "|***************** | 87% 41.34sec 4.48min"
#> [126] "|***************** | 87% 41.34sec 4.48min"
#> [127] "|****************** | 90% 30.90sec 4.63min"
#> [128] "|****************** | 90% 30.94sec 4.64min"
#> [129] "|****************** | 90% 30.96sec 4.64min"
#> [130] "|****************** | 90% 30.99sec 4.65min"
#> [131] "|******************* | 93% 20.54sec 4.79min"
#> [132] "|******************* | 93% 20.56sec 4.80min"
#> [133] "|******************* | 93% 20.59sec 4.80min"
#> [134] "|******************* | 93% 20.60sec 4.81min"
#> [135] "|******************* | 97% 10.25sec 4.95min"
#> [136] "|******************* | 97% 10.27sec 4.96min"
#> [137] "|******************* | 97% 10.29sec 4.97min"
#> [138] "|******************* | 97% 10.29sec 4.97min"
#> [139] "|********************| 100% 0.00sec 5.12min"
#> [140] ""
#> [141] "|********************| 100% 0.00sec 5.13min"
#> [142] ""
#> [143] "|********************| 100% 0.00sec 5.15min"
#> [144] ""
#> [145] "|********************| 100% 0.00sec 5.15min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 9.994 0.181 10.878saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 2396.743 17.135 1332.385