Gaussian growth model: 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.constructpath_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-22-growth-re.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, data = bwdata_model,
n.iter = 3000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4428095. logPost -2442899 logLik -2213933 edf 114.68 eps 0.0376 iteration 1"
#> [2] "AICc 4420736. logPost -2214323 logLik -2210241 edf 126.70 eps 0.0079 iteration 2"
#> [3] "AICc 4417630. logPost -2210556 logLik -2208655 edf 159.28 eps 0.0044 iteration 3"
#> [4] "AICc 4416793. logPost -2209856 logLik -2208189 edf 207.34 eps 0.0026 iteration 4"
#> [5] "AICc 4416672. logPost -2209829 logLik -2208088 edf 247.65 eps 0.0009 iteration 5"
#> [6] "AICc 4416646. logPost -2209834 logLik -2208065 edf 257.43 eps 0.0004 iteration 6"
#> [7] "AICc 4416642. logPost -2209834 logLik -2208063 edf 257.70 eps 0.0001 iteration 7"
#> [8] "AICc 4416642. logPost -2209834 logLik -2208063 edf 257.70 eps 0.0001 iteration 7"
#> [9] "elapsed time: 18.10min"
#> [10] "Starting the sampler..."
#> [11] "Starting the sampler..."
#> [12] "Starting the sampler..."
#> [13] "Starting the sampler..."
#> [14] ""
#> [15] "| | 0% 426.64min"
#> [16] "| | 0% 428.33min"
#> [17] "| | 0% 426.28min"
#> [18] "| | 0% 427.44min"
#> [19] "|* | 5% 401.91min 21.15min"
#> [20] "|* | 5% 402.07min 21.16min"
#> [21] "|* | 5% 402.89min 21.20min"
#> [22] "|* | 5% 403.13min 21.22min"
#> [23] "|** | 10% 381.38min 42.38min"
#> [24] "|** | 10% 381.85min 42.43min"
#> [25] "|** | 10% 381.82min 42.42min"
#> [26] "|** | 10% 382.09min 42.45min"
#> [27] "|*** | 15% 360.09min 63.54min"
#> [28] "|*** | 15% 360.05min 63.54min"
#> [29] "|*** | 15% 360.47min 63.61min"
#> [30] "|*** | 15% 360.52min 63.62min"
#> [31] "|**** | 20% 338.78min 84.70min"
#> [32] "|**** | 20% 338.87min 84.72min"
#> [33] "|**** | 20% 338.96min 84.74min"
#> [34] "|**** | 20% 339.20min 84.80min"
#> [35] "|***** | 25% 317.58min 105.86min"
#> [36] "|***** | 25% 317.72min 105.91min"
#> [37] "|***** | 25% 317.75min 105.92min"
#> [38] "|***** | 25% 317.73min 105.91min"
#> [39] "|****** | 30% 296.37min 127.01min"
#> [40] "|****** | 30% 296.38min 127.02min"
#> [41] "|****** | 30% 296.49min 127.07min"
#> [42] "|****** | 30% 296.62min 127.12min"
#> [43] "|******* | 35% 275.16min 148.17min"
#> [44] "|******* | 35% 275.25min 148.21min"
#> [45] "|******* | 35% 275.26min 148.22min"
#> [46] "|******* | 35% 275.38min 148.28min"
#> [47] "|******** | 40% 253.98min 169.32min"
#> [48] "|******** | 40% 254.05min 169.36min"
#> [49] "|******** | 40% 254.13min 169.42min"
#> [50] "|******** | 40% 254.18min 169.45min"
#> [51] "|********* | 45% 232.78min 190.46min"
#> [52] "|********* | 45% 232.86min 190.52min"
#> [53] "|********* | 45% 232.96min 190.60min"
#> [54] "|********* | 45% 232.96min 190.60min"
#> [55] "|********** | 50% 211.60min 211.60min"
#> [56] "|********** | 50% 211.69min 211.69min"
#> [57] "|********** | 50% 211.72min 211.72min"
#> [58] "|********** | 50% 211.78min 211.78min"
#> [59] "|*********** | 55% 190.48min 232.81min"
#> [60] "|*********** | 55% 190.51min 232.84min"
#> [61] "|*********** | 55% 190.51min 232.84min"
#> [62] "|*********** | 55% 190.59min 232.95min"
#> [63] "|************ | 60% 169.27min 253.90min"
#> [64] "|************ | 60% 169.34min 254.01min"
#> [65] "|************ | 60% 169.38min 254.07min"
#> [66] "|************ | 60% 169.39min 254.08min"
#> [67] "|************* | 65% 148.13min 275.09min"
#> [68] "|************* | 65% 148.13min 275.09min"
#> [69] "|************* | 65% 148.20min 275.23min"
#> [70] "|************* | 65% 148.25min 275.31min"
#> [71] "|************** | 70% 126.95min 296.22min"
#> [72] "|************** | 70% 126.98min 296.28min"
#> [73] "|************** | 70% 127.00min 296.34min"
#> [74] "|************** | 70% 127.08min 296.52min"
#> [75] "|*************** | 75% 105.77min 317.31min"
#> [76] "|*************** | 75% 105.82min 317.46min"
#> [77] "|*************** | 75% 105.85min 317.54min"
#> [78] "|*************** | 75% 105.88min 317.65min"
#> [79] "|**************** | 80% 84.62min 338.47min"
#> [80] "|**************** | 80% 84.66min 338.66min"
#> [81] "|**************** | 80% 84.67min 338.67min"
#> [82] "|**************** | 80% 84.71min 338.82min"
#> [83] "|***************** | 85% 63.46min 359.59min"
#> [84] "|***************** | 85% 63.49min 359.78min"
#> [85] "|***************** | 85% 63.51min 359.87min"
#> [86] "|***************** | 85% 63.54min 360.04min"
#> [87] "|****************** | 90% 42.30min 380.71min"
#> [88] "|****************** | 90% 42.32min 380.89min"
#> [89] "|****************** | 90% 42.34min 381.07min"
#> [90] "|****************** | 90% 42.36min 381.27min"
#> [91] "|******************* | 95% 21.15min 401.87min"
#> [92] "|******************* | 95% 21.16min 402.09min"
#> [93] "|******************* | 95% 21.17min 402.24min"
#> [94] "|******************* | 95% 21.18min 402.36min"
#> [95] "|********************| 100% 0.00sec 423.02min"
#> [96] ""
#> [97] "|********************| 100% 0.00sec 423.28min"
#> [98] ""
#> [99] "|********************| 100% 0.00sec 423.44min"
#> [100] ""
#> [101] "|********************| 100% 0.00sec 423.48min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 46.546 0.152 46.703saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 103496.989 72.666 27025.316