Inverse smoothing parameters
In this script we evaluate the sensitivity to hyper-parameters defined for the prior
distribution of the inverse smoothing parameters in bamlss models. We focus on the
analysis using only inverse gamma distributions as priors.
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(invgamma)Explore three types of gamma priors
igamma_mode <- function(a, b) b / (a + 1)
a <- 1e-04; b <- 1e-04
a1 <- 5; b1 <- 1000
a2 <- 10; b2 <- 0.0001
igamma_mode(a, b)#> [1] 9.999e-05igamma_mode(a1, b1)#> [1] 166.6667igamma_mode(a2, b2)#> [1] 9.090909e-06from <- 0
to <- 0.0001
n <- 20000
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
plot(function(x) dinvgamma(x, a, b), from, 0.005, n = n, main = "a) Default")
plot(function(x) dinvgamma(x, a1, b1), from, 1000, n = n, main = "b) Flexible")
plot(function(x) dinvgamma(x, a2, b2), from, 0.0001, n = n, main = "c) Smoother")
Evaluate estimation under three types of priors
d <- GAMart()
model <- bamlss(num ~ s(x2, xt = list("a" = a, "b" = a)), data = d)#> AICc 394.1316 logPost -225.594 logLik -187.075 edf 9.7746 eps 0.9573 iteration 1
#> AICc 393.8751 logPost -225.943 logLik -186.651 edf 10.058 eps 0.0316 iteration 2
#> AICc 393.8742 logPost -225.994 logLik -186.618 edf 10.089 eps 0.0031 iteration 3
#> AICc 393.8742 logPost -225.995 logLik -186.618 edf 10.090 eps 0.0000 iteration 4
#> AICc 393.8742 logPost -225.995 logLik -186.618 edf 10.090 eps 0.0000 iteration 4
#> elapsed time: 0.21sec
#> Starting the sampler...
#>
#> | | 0% 9.76sec
#> |* | 5% 9.56sec 0.50sec
#> |** | 10% 9.05sec 1.01sec
#> |*** | 15% 8.55sec 1.51sec
#> |**** | 20% 8.17sec 2.04sec
#> |***** | 25% 7.87sec 2.62sec
#> |****** | 30% 7.44sec 3.19sec
#> |******* | 35% 6.99sec 3.76sec
#> |******** | 40% 6.48sec 4.32sec
#> |********* | 45% 5.97sec 4.88sec
#> |********** | 50% 5.44sec 5.44sec
#> |*********** | 55% 4.90sec 5.99sec
#> |************ | 60% 4.44sec 6.65sec
#> |************* | 65% 3.88sec 7.22sec
#> |************** | 70% 3.33sec 7.78sec
#> |*************** | 75% 2.78sec 8.33sec
#> |**************** | 80% 2.22sec 8.90sec
#> |***************** | 85% 1.67sec 9.45sec
#> |****************** | 90% 1.11sec 10.01sec
#> |******************* | 95% 0.56sec 10.57sec
#> |********************| 100% 0.00sec 11.15secmodel1 <- bamlss(num ~ s(x2, xt = list("a" = a1, "b" = b1)), data = d)#> AICc 394.1316 logPost -292.515 logLik -187.075 edf 9.7746 eps 0.9573 iteration 1
#> AICc 393.8751 logPost -279.884 logLik -186.651 edf 10.058 eps 0.0316 iteration 2
#> AICc 393.8742 logPost -278.650 logLik -186.618 edf 10.089 eps 0.0031 iteration 3
#> AICc 393.8742 logPost -278.628 logLik -186.618 edf 10.090 eps 0.0000 iteration 4
#> AICc 393.8742 logPost -278.628 logLik -186.618 edf 10.090 eps 0.0000 iteration 4
#> elapsed time: 0.18sec
#> Starting the sampler...
#>
#> | | 0% 9.76sec
#> |* | 5% 9.52sec 0.50sec
#> |** | 10% 9.04sec 1.00sec
#> |*** | 15% 8.52sec 1.50sec
#> |**** | 20% 8.19sec 2.05sec
#> |***** | 25% 7.82sec 2.61sec
#> |****** | 30% 7.42sec 3.18sec
#> |******* | 35% 6.97sec 3.75sec
#> |******** | 40% 6.48sec 4.32sec
#> |********* | 45% 6.00sec 4.91sec
#> |********** | 50% 5.47sec 5.47sec
#> |*********** | 55% 4.93sec 6.02sec
#> |************ | 60% 4.39sec 6.58sec
#> |************* | 65% 3.86sec 7.16sec
#> |************** | 70% 3.31sec 7.73sec
#> |*************** | 75% 2.76sec 8.29sec
#> |**************** | 80% 2.21sec 8.84sec
#> |***************** | 85% 1.66sec 9.41sec
#> |****************** | 90% 1.11sec 9.97sec
#> |******************* | 95% 0.55sec 10.52sec
#> |********************| 100% 0.00sec 11.08secmodel2 <- bamlss(num ~ s(x2, xt = list("a" = a2, "b" = b2)), data = d)#> AICc 394.1316 logPost -344.746 logLik -187.075 edf 9.7746 eps 0.9573 iteration 1
#> AICc 393.8751 logPost -346.646 logLik -186.651 edf 10.058 eps 0.0316 iteration 2
#> AICc 393.8742 logPost -346.866 logLik -186.618 edf 10.089 eps 0.0031 iteration 3
#> AICc 393.8742 logPost -346.870 logLik -186.618 edf 10.090 eps 0.0000 iteration 4
#> AICc 393.8742 logPost -346.870 logLik -186.618 edf 10.090 eps 0.0000 iteration 4
#> elapsed time: 0.18sec
#> Starting the sampler...
#>
#> | | 0% 9.64sec
#> |* | 5% 9.63sec 0.51sec
#> |** | 10% 9.03sec 1.00sec
#> |*** | 15% 8.57sec 1.51sec
#> |**** | 20% 8.22sec 2.06sec
#> |***** | 25% 7.93sec 2.64sec
#> |****** | 30% 7.49sec 3.21sec
#> |******* | 35% 7.01sec 3.78sec
#> |******** | 40% 6.50sec 4.33sec
#> |********* | 45% 5.99sec 4.90sec
#> |********** | 50% 5.46sec 5.46sec
#> |*********** | 55% 4.93sec 6.03sec
#> |************ | 60% 4.40sec 6.60sec
#> |************* | 65% 3.86sec 7.17sec
#> |************** | 70% 3.32sec 7.74sec
#> |*************** | 75% 2.77sec 8.31sec
#> |**************** | 80% 2.22sec 8.88sec
#> |***************** | 85% 1.67sec 9.44sec
#> |****************** | 90% 1.11sec 10.02sec
#> |******************* | 95% 0.56sec 10.58sec
#> |********************| 100% 0.00sec 11.16secpar(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
plot(model, main = "a) Default", spar = FALSE)
plot(model1, main = "b) Flexible", spar = FALSE)
plot(model2, main = "c) Smoother", spar = FALSE)
Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 42.609 0.612 43.188