search

gamlss-dev / gamlss

gamlss: Generalized Additive Models for Location Scale and Shape
https://CRAN.R-project.org/package=gamlss
10 stars 4 forks source link

pbz smoother doesn't work with other optimization criteria than maximum likelihood #11

Open wdkrnls opened 1 month ago

wdkrnls commented 1 month ago

I noticed that pbz.control lets the user choose a different objective criterion for the smooth than the default method = "ML". However, in my experiments I couldn't get anything other than that to work. When I try I get an error which seems to be buried deep in the code:

Error in regpen(y = y, X = X, w = w, lambda = lambda) : 
  unused argument (lambda = lambda)

Here is an example dset:

dset = structure(list(x = c(0, 7, 14, 21, 28, 35, 42, 49, 56, 63, 
70, 77, 84, 91, 98, 105, 112, 119, 126, 133, 140, 147, 154, 161, 
168, 175, 182, 189, 196, 203, 210, 217, 224, 231, 238, 245, 252, 
259, 266, 273, 280, 287, 294, 301, 308, 315, 322, 329, 336, 343, 
350, 357, 364, 371, 378, 385, 392, 399, 406, 413, 420, 427, 434, 
441, 448, 455, 462, 469, 476, 483, 490, 497, 504, 511, 518, 525, 
532, 539, 546, 553, 560, 567, 574, 581, 588, 595, 602, 609, 616, 
623, 630, 637, 644, 651, 658, 665, 672, 679, 686, 693, 700, 707, 
714, 721, 728, 735, 742, 749, 756, 763, 770, 777, 784, 791, 798, 
805, 812, 819, 826, 833, 840, 847, 854, 861, 868, 875, 882, 889, 
896, 903, 910, 917, 924, 931, 938, 945, 952, 959, 966, 973, 980, 
987, 994, 1001, 1008, 1015, 1022, 1029, 1036, 1043, 1050, 1057, 
1064, 1071, 1078, 1085), y = c(1.01444077846316, 1.05456445827494, 
0.97287921791233, 1.00980800370113, 1.02615781424756, 1.03140374983378, 
0.983131497160004, 1.04500297376479, 1.01854433320661, 0.991399176954732, 
0.970930716127422, 1.03258528600965, 1.02914821513961, 1.02360052667278, 
1.00281367640781, 1.02611617858857, 1.02594954783436, 0.922286760174066, 
0.978756402669224, 1.03806069701517, 1.01494516090153, 1.03775267538644, 
0.94727619047619, 1.01443738438028, 0.993187309383293, 0.94292314739211, 
0.964825077202927, 1.02083105162859, 1.01461900038812, 1.05809749941175, 
1.04300316525877, 0.803824787565569, 0.958407036163522, 0.958861172783074, 
1.00605678579259, 1.11262752333406, 1.02003119212399, 0.957224357291194, 
1.05332801714471, 0.999162559219027, 0.96284019120367, 1.10368058651607, 
0.991145809623526, 0.991963917761802, 0.995061104474869, 0.971623788607683, 
0.968167403421103, 0.963117976815564, 0.988181986804102, 0.997126047230431, 
1.00526850149549, 1.04205305154195, 0.950200193651375, 0.90526572139842, 
1.00994041922638, 1.06532177694561, 1.00272556924649, 1.05618843683084, 
0.942649876813578, 1.05000722647781, 1.06416240850711, 1.00366905155017, 
0.990322412353923, 0.933763157894737, 0.967994513345145, 0.988439306358382, 
1.05282652255753, 0.888038802436701, 0.949315068493151, 0.975207746720692, 
1.03482030485247, 0.953257169817679, 0.991668708649841, 1.03873734533183, 
1.02807504741263, 1.15178542043772, 0.949932067181222, 1.03739862280031, 
0.250636282923942, 0.970545804279875, 0.99794064203513, 1.01112197506722, 
1.00181312703977, 0.984020987359886, 0.969773299748111, 1.00493461633358, 
0.995219416523658, 0.980188211986132, 0.973210868733257, 0.991774383078731, 
0.969243113131854, 1.06958304853042, 0.921628889742605, 1.00780487804878, 
1.00041373603641, 1.01449473540271, 1.02618014107434, 1.12915617962644, 
0.97591488366655, 0.956417624521073, 1.03860569715142, 0.978103946102021, 
1.03173431734317, 0.996063461767864, 0.990453460620525, 0.949860724233983, 
1.00358010484593, 0.996229358991029, 1.06775670372793, 0.995686999383857, 
1.03655384615385, 1.06401713469776, 0.939244479318462, 0.982542149946192, 
1.03166259168704, 0.934878719082328, 1.020366598778, 0.996962025316456, 
0.915370889625287, 1.01782752902156, 0.948051948051948, 0.986521739130435, 
1.00486008836524, 1.00474284867348, 0.94331566548881, 0.934547166849788, 
0.980791715383331, 0.979053133514986, 1.070301542777, 1.04745540828015, 
1.01286274509804, 1.06617761823006, 1.10756449497158, 1.01359832635983, 
0.993561679114087, 0.930741031683687, 1.01376597836775, 0.988334504567814, 
1.02382294897096, 0.913977005125364, 1.03095129722349, 1.05154185022026, 
1.018052594172, 1.03636616970879, 1.03629520096787, 0.944516801250326, 
0.986479028697572, 0.981318834518333, 0.946855123674912, 0.992120130835563, 
0.923580136782634, 1.0243044451551, 1.20034166796086, 1.04426948890036, 
1.02027027027027, 0.970683647723177)), row.names = c(NA, -156L
), class = c("tbl_df", "tbl", "data.frame"))

Here is a minimal example. This works:

gamlss(formula = y ~ pbz(x), sigma.formula = ~ pbz(x), nu.formula = ~ pbz(x, control = pbz.control(method = "ML")), family = BCPEo(), control = gamlss.control(nu.step = 0.01), method = mixed(10, 20), data = dset)

However, this triggers the error:

gamlss(formula = y ~ pbz(x), sigma.formula = ~ pbz(x), nu.formula = ~ pbz(x, control = pbz.control(method = "GAIC")), family = BCPEo(), control = gamlss.control(nu.step = 0.01), method = mixed(10, 20), data = dset)

I was interested in "GAIC" for parsimoniously describing the quantiles of the data set. Note that there is a stark drop in the response, which has trouble with the default optimization settings.