How to estimate the standard error of the percentiles with gamlss::centiles.pred

Question

I have fitted a model with

m0 <- lms(weight,age, data=df)

After, I use the model to predict percentiles for same ages:

newages <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)
centiles.pred(m0, xname="age", xvalues= newages, cent= c(5,25,50,75,95))

Is it possible to estimate the error of the percentiles predicted? How?

Below a reproducible examples with data include in the gamlss package:

library("gamlss")
library("gamlss.data")
data("dbbmi")
m0 <- lms(bmi,age, data=dbbmi, trans.x=TRUE, k=2)
newages <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)
percentiles <- centiles.pred(m0, xname="age", xvalues= newages, cent= c(5,25,50,75,95))

Answer 1

In principle, you could try bootstrapping them. I would parallelize the code since the package is rather slow. On Linux we can conveniently use parallel::mclapply, which is the parallel version of lapply. We might silence the function a little using capture.output. I demonstrate this with example from ?centiles.pred.

library('gamlss')
a <- gamlss(y ~ pb(x), sigma.fo=~pb(x), data=abdom, family=BCT)
newx <- seq.int(12, 40, 2)
cent <- centiles.pred(a, xname="x", xvalues=newx)

library(parallel)
set.seed(42)
system.time(
  b <- mclapply(seq_len(299L), \(i) {
    X <- abdom[sample.int(nrow(abdom), nrow(abdom), replace=TRUE), ]
    nul <- capture.output(a <- gamlss(y ~ pb(x), sigma.fo=~pb(x), data=X, family=BCT))
    newx <- seq.int(12, 40, 2)
    nul <- capture.output(r <- as.matrix(centiles.pred(a, xname="x", xvalues=newx)[-1]))
    `rownames<-`(r, newx)
  }, mc.cores=7L) |> simplify2array()
)
#    user  system elapsed 
# 743.861   7.524 147.565 

Using simplify2array gives us an "array" on which we can apply sd to get bootstrap standard errors,

se <- apply(b, c(1, 2), sd)

which we can put into a neat list.

(res_lst <- list(centiles=`rownames<-`(cent[-1], cent$x), error=se))
# $centiles
#          0.4         2        10        25        50        75        90        98      99.6
# 12  42.18638  45.94808  50.55402  54.10164  58.10869  62.34625  66.52695  72.67661  78.43236
# 14  63.42588  68.02193  73.57773  77.80866  82.54190  87.49937  92.34611  99.40311 105.93746
# 16  85.55128  90.63694  96.72240 101.31493 106.41336 111.71203 116.85444 124.28015 131.09580
# 18 107.59336 113.00638 119.43694 124.25863 129.58205 135.08375 140.39497 148.01855 154.97127
# 20 128.88564 134.65890 141.48729 146.58710 152.19862 157.97822 163.53955 171.49265 178.71729
# 22 149.23815 155.45232 162.78365 168.24664 174.24610 180.41304 186.33591 194.78788 202.44813
# 24 168.66036 175.43397 183.41582 189.35719 195.87604 202.57063 208.99454 218.15225 226.44314
# 26 187.33960 194.83997 203.67733 210.25491 217.47125 224.88154 231.99167 242.12678 251.30172
# 28 205.35844 213.75354 223.65182 231.02355 239.11539 247.42914 255.41020 266.79334 277.10445
# 30 222.60262 231.96987 243.02502 251.26548 260.31763 269.62507 278.56649 291.32982 302.90129
# 32 239.39290 249.66070 261.78656 270.83041 280.77006 290.99522 300.82306 314.85943 327.59258
# 34 255.89293 266.91077 279.92414 289.63108 300.30058 311.27771 321.82933 336.90114 350.57522
# 36 271.73865 283.47606 297.34085 307.68389 319.05355 330.75203 341.99796 358.06301 372.63969
# 38 286.19116 298.83107 313.77356 324.92838 337.19773 349.82961 361.97990 379.34828 395.11866
# 40 299.29953 313.05326 329.33575 341.50653 354.90804 368.72097 382.02141 401.05688 418.36327
# 
# $error
#         0.4        2       10        25        50        75        90        98      99.6
# 12 8.754073 9.058913 9.772938 10.445019 11.292740 12.268047 13.285140 14.890101 16.699377
# 14 8.302389 8.584919 9.173224  9.690382 10.314486 11.007560 11.707848 12.786696 14.049017
# 16 8.225677 8.403785 8.823087  9.191237  9.630658 10.112131 10.589999 11.319983 12.232993
# 18 8.104683 8.211071 8.524018  8.802718  9.131756  9.485537  9.829230 10.351928 11.051896
# 20 7.913490 7.973818 8.212644  8.426547  8.673893  8.932846  9.178813  9.556925 10.111166
# 22 7.804297 7.807309 7.971309  8.128343  8.311364  8.503114  8.687357  8.989410  9.481017
# 24 7.724357 7.680904 7.792952  7.915895  8.065365  8.227356  8.390109  8.678840  9.171750
# 26 7.649060 7.594180 7.693274  7.808808  7.952233  8.109245  8.266915  8.543402  9.011338
# 28 7.640950 7.599286 7.706000  7.824723  7.970978  8.128637  8.279786  8.520554  8.915502
# 30 7.692231 7.660621 7.777694  7.906453  8.069103  8.248826  8.420955  8.675671  9.053343
# 32 7.811806 7.792124 7.930340  8.081476  8.278466  8.505167  8.731095  9.068910  9.523329
# 34 7.863307 7.907566 8.128554  8.343744  8.613907  8.919423  9.224228  9.679546 10.248512
# 36 7.954630 8.077143 8.398259  8.690299  9.045676  9.439037  9.828830 10.409909 11.104851
# 38 8.265592 8.432389 8.813371  9.153125  9.560848 10.005845 10.443078 11.091570 11.853296
# 40 8.834803 9.000554 9.402092  9.768263 10.208644 10.686466 11.154407 11.851206 12.664367