Full model with t-distribution: 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-10-full-re-t.rds")
model_file_burned <- gsub("(\\.rds)$", "-burned\\1", model_file)

bwdata_model <- fst::read_fst(bwdata_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 = 1000, burnin = 0)
#> ---
#> Family: TF 
#> Link function: mu = identity, sigma = log, nu = log
#> *---
#> Formula mu:
#> ---
#> born_weight ~ sex + born_race + birth_place + 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") + 
#>     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)
#> -
#> Parametric coefficients:
#>                                      Mean     2.5%      50%    97.5% parameters
#> (Intercept)                      3098.171 3040.026 3097.444 3157.314   3127.974
#> sexfemale                        -100.385 -103.728 -100.400  -97.037   -100.296
#> born_raceindigenous               -57.535  -62.931  -57.499  -52.200    -57.616
#> born_raceNA                         4.253  -16.417    4.248   25.132      3.908
#> birth_placeanother health center    6.086  -24.108    6.249   36.214      6.186
#> birth_placehome                   -74.150  -79.353  -74.109  -68.995    -73.976
#> birth_placeother                  -87.903 -108.560  -87.588  -66.551    -88.372
#> birth_placeNA                    -114.229 -133.157 -114.416  -95.789   -113.776
#> marital_statusmarried              31.571   26.119   31.549   37.191     31.602
#> marital_statuswidow               -10.857  -59.360  -12.091   41.909    -10.646
#> marital_statusdivorced            -28.199  -73.846  -28.416   16.988    -28.633
#> marital_statusconsensual union     26.635   22.169   26.712   31.098     26.718
#> marital_statusNA                   23.514    7.630   23.654   39.569     23.824
#> study_years1 - 3                   63.025   54.600   63.152   71.531     63.155
#> study_years4 - 7                   84.305   76.040   84.312   92.507     84.357
#> study_years8 - 11                  77.240   68.753   77.242   85.401     77.242
#> study_years> 12                    75.968   65.154   75.837   87.558     76.017
#> study_yearsNA                      51.053   29.224   51.322   71.362     50.534
#> consult_num1 - 3                   49.722   42.390   49.681   57.625     49.774
#> consult_num4 - 6                   91.801   84.511   91.902   98.832     91.803
#> consult_num> 7                    150.700  142.876  150.730  158.690    150.667
#> consult_numNA                      16.570   -4.741   16.667   36.959     16.661
#> remoteness                        -12.001  -72.632  -12.886   49.346    -19.643
#> prop_tap_toilet                   106.267  -16.694  106.749  226.257    -75.855
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.9735 0.8658 0.9979     1
#> -
#> Smooth terms:
#>                                                           Mean      2.5%       50%     97.5%
#> s(age).tau21                                         1.812e+05 5.026e+04 1.389e+05 6.040e+05
#> s(age).alpha                                         9.762e-01 8.452e-01 9.992e-01 1.000e+00
#> s(age).edf                                           6.619e+00 5.917e+00 6.604e+00 7.400e+00
#> s(wk_ini).tau21                                      7.258e+03 1.651e+02 2.413e+03 4.212e+04
#> s(wk_ini).alpha                                      9.979e-01 9.873e-01 1.000e+00 1.000e+00
#> s(wk_ini).edf                                        5.400e+00 2.885e+00 5.207e+00 8.241e+00
#> s(rivwk_conception).tau21                            8.820e-01 4.163e-03 4.466e-01 4.319e+00
#> s(rivwk_conception).alpha                            9.993e-01 9.950e-01 1.000e+00 1.000e+00
#> s(rivwk_conception).edf                              2.445e+00 1.767e-01 2.414e+00 4.572e+00
#> s(res_muni).tau21                                    1.823e+03 1.172e+03 1.748e+03 2.845e+03
#> s(res_muni).alpha                                    9.822e-01 9.111e-01 9.978e-01 1.000e+00
#> s(res_muni).edf                                      4.186e+01 4.129e+01 4.188e+01 4.230e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21         5.197e+02 7.662e-01 3.501e+02 2.159e+03
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha         9.926e-01 9.594e-01 9.990e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf           1.083e+01 2.087e+00 1.063e+01 2.282e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21       1.543e+03 4.377e+02 1.355e+03 3.559e+03
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha       9.880e-01 9.382e-01 9.998e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf         2.207e+01 1.398e+01 2.191e+01 3.117e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 5.760e+04 2.945e+04 5.605e+04 9.661e+04
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 9.256e-01 6.866e-01 9.874e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf   5.229e+01 4.645e+01 5.244e+01 5.739e+01
#>                                                      parameters
#> s(age).tau21                                          3.874e+04
#> s(age).alpha                                                 NA
#> s(age).edf                                            5.714e+00
#> s(wk_ini).tau21                                       1.079e+05
#> s(wk_ini).alpha                                              NA
#> s(wk_ini).edf                                         8.619e+00
#> s(rivwk_conception).tau21                             3.250e-01
#> s(rivwk_conception).alpha                                    NA
#> s(rivwk_conception).edf                               2.160e+00
#> s(res_muni).tau21                                     3.514e+03
#> s(res_muni).alpha                                            NA
#> s(res_muni).edf                                       4.243e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21          7.230e+02
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha                 NA
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf            1.469e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21        2.620e+03
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha               NA
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf          2.800e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21  1.834e+05
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha         NA
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf    6.163e+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.1216358  6.1023620  6.1217264  6.1410523      0.112
#> sexfemale                        -0.0433255 -0.0492370 -0.0433804 -0.0372285     -0.043
#> born_raceindigenous              -0.0705560 -0.0814157 -0.0706281 -0.0597724     -0.071
#> born_raceNA                      -0.0064638 -0.0446406 -0.0060706  0.0340295     -0.006
#> birth_placeanother health center  0.0434365 -0.0047400  0.0443900  0.0977410      0.043
#> birth_placehome                   0.0255851  0.0154573  0.0255527  0.0352919      0.026
#> birth_placeother                  0.0715635  0.0380244  0.0710848  0.1079411      0.071
#> birth_placeNA                     0.0244067 -0.0098765  0.0245579  0.0618413      0.025
#> study_years1 - 3                  0.0002604 -0.0163886  0.0005253  0.0157620      0.000
#> study_years4 - 7                  0.0039563 -0.0115420  0.0041544  0.0186836      0.004
#> study_years8 - 11                 0.0218914  0.0053769  0.0218823  0.0374032      0.022
#> study_years> 12                   0.0386046  0.0193703  0.0382010  0.0592433      0.039
#> study_yearsNA                    -0.0099752 -0.0474218 -0.0103115  0.0281008     -0.010
#> consult_num1 - 3                 -0.0544704 -0.0686053 -0.0546685 -0.0398027     -0.054
#> consult_num4 - 6                 -0.1068352 -0.1205922 -0.1067169 -0.0920768     -0.107
#> consult_num> 7                   -0.1436410 -0.1582454 -0.1437191 -0.1278226     -0.144
#> consult_numNA                    -0.0151018 -0.0504646 -0.0155573  0.0189507     -0.015
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.9571 0.7741 0.9985     1
#> -
#> Smooth terms:
#>                        Mean      2.5%       50%     97.5% parameters
#> s(age).tau21       0.035074  0.003317  0.020955  0.140782      0.023
#> s(age).alpha       0.987775  0.903588  1.000000  1.000000         NA
#> s(age).edf         4.251902  2.802336  4.201711  5.931355      4.314
#> s(wk_ini).tau21    0.021786  0.002221  0.015113  0.079930      0.035
#> s(wk_ini).alpha    0.993307  0.960458  1.000000  1.000000         NA
#> s(wk_ini).edf      5.894105  3.965739  5.921776  7.738381      6.859
#> s(res_muni).tau21 33.630737 21.538703 32.615577 52.416922     36.391
#> s(res_muni).alpha  0.930053  0.678384  0.988716  1.000000         NA
#> s(res_muni).edf   42.999774 42.999655 42.999779 42.999860     43.000
#> ---
#> Formula nu:
#> ---
#> nu ~ 1
#> -
#> Parametric coefficients:
#>              Mean  2.5%   50% 97.5% parameters
#> (Intercept) 1.877 1.854 1.877 1.899      1.881
#> -
#> Acceptance probabilty:
#>         Mean   2.5%    50% 97.5%
#> alpha 0.8838 0.3908 0.9789     1
#> ---
#> Sampler summary:
#> -
#> runtime = 27704.27
#> ---
#> Optimizer summary:
#> -
#> AICc = 4420293 edf = 259.4192 logLik = -2209887
#> logPost = -2211647 nobs = 291479 runtime = 5561.42

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.987   0.338  10.359