Growth model with t-distribution and no rainfall: 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.construct
library(gamlss.dist)
#> Loading required package: MASS
path_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")

source(file.path(path_proj, "src", "51-bamlss.R"))

bwdata_file <- file.path(path_processed, "bwdata_41_model.fst")
model_file <- file.path(path_modelled, "bw-25-growth-no-rain-t.rds")
form_file <- gsub("(\\.rds)$", "-form\\1", model_file)
model_file_burned <- gsub("(\\.rds)$", "-burned\\1", model_file)

bwdata_model <- fst::read_fst(bwdata_file)
form <- readRDS(form_file)
model <- readRDS(model_file_burned)

Summary

summary(model)
#> 
#> Call:
#> bamlss(formula = form, family = TF, data = bwdata_model, cores = 4, 
#>     combine = FALSE, light = TRUE, n.iter = 8000, burnin = 0)
#> ---
#> Family: TF 
#> Link function: mu = identity, sigma = log, nu = log
#> *---
#> Formula 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")
#> -
#> Parametric coefficients:
#>                                       Mean      2.5%       50%     97.5% parameters
#> (Intercept)                       3102.283  3041.208  3103.395  3157.566   3134.267
#> sexfemale                         -101.012  -104.264  -100.962   -97.606   -101.046
#> born_raceindigenous                -50.208   -55.586   -50.196   -44.773    -50.152
#> born_raceNA                          4.859   -14.839     4.964    23.863      4.666
#> birth_placeanother health center    11.923   -15.863    11.941    40.659     11.663
#> birth_placehome                    -62.748   -68.342   -62.687   -57.363    -62.693
#> birth_placeother                   -71.160   -90.462   -71.209   -52.012    -71.156
#> birth_placeNA                     -102.697  -123.639  -102.379   -81.951   -102.492
#> gestation_weeks32 - 36            -246.873  -253.327  -246.916  -240.382   -247.049
#> gestation_weeks28 - 31            -853.165  -879.892  -853.390  -827.124   -853.838
#> gestation_weeks22 - 27           -1772.918 -1818.968 -1772.764 -1723.943  -1772.845
#> marital_statusmarried               30.936    25.668    30.953    36.364     30.989
#> marital_statuswidow                -11.533   -63.998   -11.915    39.861    -12.393
#> marital_statusdivorced             -32.859   -79.202   -31.474    12.632    -32.821
#> marital_statusconsensual union      24.213    19.819    24.230    28.767     24.227
#> marital_statusNA                    22.372     8.155    22.190    38.418     22.043
#> study_years1 - 3                    60.665    51.959    60.824    68.938     60.675
#> study_years4 - 7                    81.690    74.028    81.658    89.794     81.721
#> study_years8 - 11                   75.396    67.095    75.211    84.003     75.487
#> study_years> 12                     72.494    62.341    72.384    83.240     72.660
#> study_yearsNA                       47.010    25.328    47.127    68.787     46.973
#> consult_num1 - 3                    66.813    58.554    66.920    74.139     66.866
#> consult_num4 - 6                   101.831    93.998   101.959   108.926    101.974
#> consult_num> 7                     146.946   138.967   147.016   154.609    147.100
#> consult_numNA                       21.169     2.092    21.141    41.225     21.102
#> remoteness                          -4.895   -66.232    -5.968    53.484    -10.860
#> prop_tap_toilet                    132.840    17.150   130.836   264.759    -70.476
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.8290 0.2997 0.9297     1
#> -
#> Smooth terms:
#>                                Mean      2.5%       50%     97.5% parameters
#> s(age).tau21              1.644e+05 4.203e+04 1.282e+05 5.043e+05  4.309e+04
#> s(age).alpha              9.809e-01 8.765e-01 9.995e-01 1.000e+00         NA
#> s(age).edf                6.516e+00 5.780e+00 6.519e+00 7.329e+00  5.784e+00
#> s(wk_ini).tau21           9.611e+04 2.649e+04 7.651e+04 2.897e+05  2.016e+05
#> s(wk_ini).alpha           9.966e-01 9.810e-01 9.999e-01 1.000e+00         NA
#> s(wk_ini).edf             8.461e+00 7.830e+00 8.501e+00 8.867e+00  8.787e+00
#> s(rivwk_conception).tau21 2.232e+00 1.887e-01 1.363e+00 1.039e+01  1.909e+00
#> s(rivwk_conception).alpha 9.991e-01 9.946e-01 1.000e+00 1.000e+00         NA
#> s(rivwk_conception).edf   3.481e+00 1.808e+00 3.437e+00 5.608e+00  3.749e+00
#> s(res_muni).tau21         1.908e+03 1.228e+03 1.841e+03 2.884e+03  4.152e+03
#> s(res_muni).alpha         9.857e-01 9.307e-01 9.999e-01 1.000e+00         NA
#> s(res_muni).edf           4.192e+01 4.142e+01 4.193e+01 4.232e+01  4.252e+01
#> ---
#> Formula sigma:
#> ---
#> sigma ~ sex + born_race + birth_place + study_years + consult_num + 
#>     s(age) + s(wk_ini) + s(res_muni, bs = "re")
#> -
#> Parametric coefficients:
#>                                        Mean       2.5%        50%      97.5% parameters
#> (Intercept)                       6.1077914  6.0893580  6.1077604  6.1272758      0.082
#> sexfemale                        -0.0410126 -0.0463576 -0.0410585 -0.0352277     -0.041
#> born_raceindigenous              -0.0566448 -0.0664529 -0.0566675 -0.0461373     -0.057
#> born_raceNA                      -0.0008566 -0.0385714 -0.0008478  0.0346444     -0.001
#> birth_placeanother health center  0.0630262  0.0114596  0.0632751  0.1117714      0.063
#> birth_placehome                   0.0522655  0.0421192  0.0523469  0.0623431      0.052
#> birth_placeother                  0.0610985  0.0267787  0.0612330  0.0954881      0.061
#> birth_placeNA                     0.0869311  0.0513823  0.0865297  0.1217677      0.087
#> study_years1 - 3                  0.0026571 -0.0131167  0.0026073  0.0183816      0.002
#> study_years4 - 7                  0.0044598 -0.0111246  0.0042728  0.0188575      0.004
#> study_years8 - 11                 0.0149874 -0.0011356  0.0150063  0.0308461      0.015
#> study_years> 12                   0.0239234  0.0034686  0.0237533  0.0446986      0.024
#> study_yearsNA                    -0.0222662 -0.0611596 -0.0229173  0.0162926     -0.023
#> consult_num1 - 3                 -0.0385886 -0.0520717 -0.0387949 -0.0246325     -0.038
#> consult_num4 - 6                 -0.0803982 -0.0930151 -0.0803824 -0.0670049     -0.080
#> consult_num> 7                   -0.1071433 -0.1215766 -0.1071288 -0.0933517     -0.107
#> consult_numNA                     0.0032929 -0.0318429  0.0034741  0.0375507      0.002
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.9628 0.8114 0.9932     1
#> -
#> Smooth terms:
#>                        Mean      2.5%       50%     97.5% parameters
#> s(age).tau21      9.550e-03 2.378e-04 4.198e-03 5.701e-02      0.005
#> s(age).alpha      9.922e-01 9.406e-01 9.997e-01 1.000e+00         NA
#> s(age).edf        3.054e+00 1.409e+00 2.994e+00 5.179e+00      3.141
#> s(wk_ini).tau21   1.318e-02 7.317e-04 8.639e-03 5.202e-02      0.018
#> s(wk_ini).alpha   9.959e-01 9.749e-01 9.999e-01 1.000e+00         NA
#> s(wk_ini).edf     5.290e+00 2.991e+00 5.344e+00 7.338e+00      6.173
#> s(res_muni).tau21 3.467e+01 2.243e+01 3.403e+01 5.190e+01     36.211
#> s(res_muni).alpha 9.359e-01 6.752e-01 9.982e-01 1.000e+00         NA
#> s(res_muni).edf   4.300e+01 4.300e+01 4.300e+01 4.300e+01     43.000
#> ---
#> Formula nu:
#> ---
#> nu ~ 1
#> -
#> Parametric coefficients:
#>              Mean  2.5%   50% 97.5% parameters
#> (Intercept) 2.117 2.087 2.117 2.144      2.121
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.8851 0.3133 0.9895     1
#> ---
#> Sampler summary:
#> -
#> runtime = 52206.79
#> ---
#> Optimizer summary:
#> -
#> AICc = 4408851 edf = 158.1527 logLik = -2204267
#> logPost = -2205082 nobs = 291479 runtime = 286.17

Parametric effects

par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot2d_bamlss(model, bwdata_model, model = "mu", term = "remoteness", grid = 50,
              FUN = c95)
plot2d_bamlss(model, bwdata_model, model = "mu", term = "prop_tap_toilet", grid = 50,
              FUN = c95)

Smooth effects

par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot(model, scale = 0, scheme = 2, spar = FALSE, image = TRUE, contour = TRUE)

Time to execute the task

Only useful when executed with Rscript.

proc.time()
#>    user  system elapsed 
#>   9.762   0.319  10.166