Spatial model with random effects: fitting
Birthweight model including covariates at municipality level with linear effects.
Load packages, read data and source custom scripts
Paths are defined relative to the git repository location.
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-muni-16-sp-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
Now we define the same models as in the previous study.
form_sigma <- sigma ~ 1
form_mu <- born_weight ~ s(longitud, latitud) + 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 = 1000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4466655. logPost -2455255 logLik -2233322 edf 5.2735 eps 0.0328 iteration 1"
#> [2] "AICc 4457952. logPost -2231199 logLik -2228962 edf 14.041 eps 0.0105 iteration 2"
#> [3] "AICc 4456528. logPost -2228805 logLik -2228225 edf 39.119 eps 0.0034 iteration 3"
#> [4] "AICc 4456411. logPost -2228508 logLik -2228145 edf 60.255 eps 0.0011 iteration 4"
#> [5] "AICc 4456417. logPost -2228498 logLik -2228144 edf 64.079 eps 0.0000 iteration 5"
#> [6] "AICc 4456417. logPost -2228498 logLik -2228144 edf 64.079 eps 0.0000 iteration 5"
#> [7] "elapsed time: 28.67sec"
#> [8] "Starting the sampler..."
#> [9] "Starting the sampler..."
#> [10] "Starting the sampler..."
#> [11] "Starting the sampler..."
#> [12] ""
#> [13] "| | 0% 18.11min"
#> [14] "| | 0% 18.30min"
#> [15] "| | 0% 18.29min"
#> [16] "| | 0% 18.38min"
#> [17] ""
#> [18] "|* | 5% 16.67min 52.64sec|* | 5% 16.65min 52.59sec"
#> [19] "|* | 5% 16.64min 52.54sec"
#> [20] "|* | 5% 16.67min 52.65sec"
#> [21] "|** | 10% 15.58min 1.73min"
#> [22] "|** | 10% 15.58min 1.73min"
#> [23] "|** | 10% 15.59min 1.73min"
#> [24] "|** | 10% 15.61min 1.73min"
#> [25] "|*** | 15% 14.76min 2.60min"
#> [26] "|*** | 15% 14.78min 2.61min"
#> [27] "|*** | 15% 14.81min 2.61min"
#> [28] "|*** | 15% 14.82min 2.61min"
#> [29] "|**** | 20% 13.88min 3.47min"
#> [30] "|**** | 20% 13.89min 3.47min"
#> [31] "|**** | 20% 13.89min 3.47min"
#> [32] "|**** | 20% 13.91min 3.48min"
#> [33] ""
#> [34] "|***** | 25%|***** | 25% 12.99min12.98min 4.33min 4.33min"
#> [35] "|***** | 25% 13.03min 4.34min"
#> [36] "|***** | 25% 13.04min 4.35min"
#> [37] "|****** | 30% 12.16min 5.21min"
#> [38] "|****** | 30% 12.17min 5.21min"
#> [39] "|****** | 30% 12.21min 5.23min"
#> [40] "|****** | 30% 12.23min 5.24min"
#> [41] "|******* | 35% 11.38min 6.13min"
#> [42] "|******* | 35% 11.41min 6.15min"
#> [43] "|******* | 35% 11.45min 6.17min"
#> [44] "|******* | 35% 11.49min 6.19min"
#> [45] "|******** | 40% 10.62min 7.08min"
#> [46] "|******** | 40% 10.64min 7.10min"
#> [47] "|******** | 40% 10.65min 7.10min"
#> [48] "|******** | 40% 10.68min 7.12min"
#> [49] "|********* | 45% 9.77min 7.99min"
#> [50] "|********* | 45% 9.86min 8.07min"
#> [51] "|********* | 45% 9.87min 8.07min"
#> [52] "|********* | 45% 9.93min 8.12min"
#> [53] "|********** | 50% 8.95min 8.95min"
#> [54] "|********** | 50% 9.03min 9.03min"
#> [55] "|********** | 50% 9.05min 9.05min"
#> [56] "|********** | 50% 9.06min 9.06min"
#> [57] "|*********** | 55% 8.06min 9.85min"
#> [58] "|*********** | 55% 8.12min 9.92min"
#> [59] "|*********** | 55% 8.13min 9.94min"
#> [60] "|*********** | 55% 8.15min 9.96min"
#> [61] "|************ | 60% 7.18min 10.78min"
#> [62] "|************ | 60% 7.22min 10.83min"
#> [63] "|************ | 60% 7.24min 10.86min"
#> [64] "|************ | 60% 7.25min 10.88min"
#> [65] "|************* | 65% 6.31min 11.72min"
#> [66] "|************* | 65% 6.33min 11.76min"
#> [67] "|************* | 65% 6.35min 11.79min"
#> [68] "|************* | 65% 6.35min 11.80min"
#> [69] "|************** | 70% 5.41min 12.61min"
#> [70] "|************** | 70% 5.42min 12.65min"
#> [71] "|************** | 70% 5.46min 12.73min"
#> [72] "|************** | 70% 5.46min 12.74min"
#> [73] "|*************** | 75% 4.50min 13.51min"
#> [74] "|*************** | 75% 4.52min 13.55min"
#> [75] "|*************** | 75% 4.54min 13.63min"
#> [76] "|*************** | 75% 4.55min 13.64min"
#> [77] "|**************** | 80% 3.60min 14.40min"
#> [78] "|**************** | 80% 3.61min 14.43min"
#> [79] "|**************** | 80% 3.63min 14.52min"
#> [80] "|**************** | 80% 3.63min 14.52min"
#> [81] "|***************** | 85% 2.70min 15.28min"
#> [82] "|***************** | 85% 2.70min 15.31min"
#> [83] "|***************** | 85% 2.72min 15.40min"
#> [84] "|***************** | 85% 2.72min 15.40min"
#> [85] "|****************** | 90% 1.79min 16.15min"
#> [86] "|****************** | 90% 1.80min 16.19min"
#> [87] "|****************** | 90% 1.81min 16.27min"
#> [88] "|****************** | 90% 1.81min 16.27min"
#> [89] "|******************* | 95% 53.75sec 17.02min"
#> [90] "|******************* | 95% 53.89sec 17.06min"
#> [91] "|******************* | 95% 54.16sec 17.15min"
#> [92] "|******************* | 95% 54.17sec 17.15min"
#> [93] "|********************| 100% 0.00sec 17.91min"
#> [94] ""
#> [95] "|********************| 100% 0.00sec 17.94min"
#> [96] ""
#> [97] "|********************| 100% 0.00sec 18.00min"
#> [98] ""
#> [99] "|********************| 100% 0.00sec 18.00min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 2.244 0.004 2.258saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 4404.543 11.836 1142.992