Smoothers are non-parametric techniques developed mostly in the 1980 and 90’s. The main advantage of a smoother is to suggest possible functional forms on how the explanatory variables effect the parameters of the response. That is, to let the data to detect non-linearities in the model. There are basically two types of smoothers;
the cubic splines smoothers and techniques associated with them like;
the local polynomial smoothing, Cleveland, Grosse, and Shyu (1993).
There are several ways to use smoothers within the original GAMLSS framework and they are described in Chapter 9 of Stasinopoulos et al. (2017). Local polynomials smoothers can be used with gamlss2 by connecting the R function loess() with gamlss2 using specials (extra additive terms). This connection is described in a different vignette (Mikis reference here). Here we give an explanation of to how use the smooth terms s() and pb() within gamlss2.
s() functionThe function s() is using penalized regression smoothers and it is identical to the function s() used in the R package mgcv. The first argument of the function specifies the explanatory term(s). For one dimensional smoother the function can be used as;
library(gamlss2)
library(gamlss)
gm11 <- gamlss2(rent~s(area)+s(yearc)+location+bath+kitchen|
s(area)+s(yearc)+location+bath+kitchen,
data=rent99, family=GA)GAMLSS-RS iteration 1: Global Deviance = 38525.1892 eps = 0.123444
GAMLSS-RS iteration 2: Global Deviance = 38471.5199 eps = 0.001393
GAMLSS-RS iteration 3: Global Deviance = 38471.1817 eps = 0.000008 We use smooth terms for the two continuous explanatory variables of the data, the size of the flats, area, and the year of construction, yearc. To visualise the fitted smooth functions use;
plot(gm11)
The first two plots are for the \(\mu\) model while the last two plots for the \(\sigma\) model. It seems that both area, and yearc, need smoothing functions for the \(\mu\) while only yearc need a smoother for the \(\sigma\).
Note that the first argument for s() is the continuous variables needed smoothing. Here we used single terms. Single terms represent main effects. More that one terms can be use to represent iteractions. Note that all argument of the s() function of the mgcv package apply here. Here we mention the most important ones;
k; the dimension of the basis;
bs; indicating the (penalized) smoothing basis to use, e.g. “tp” for thin plate regression spline, “cr” for cubic regression spline etc;
m; The order of the penalty e.g. 2 for cubic spline penalty;
by; a term for varying coefficient model (a specific form of iteraction);
to find information about all the argument please try;
?mgcv:::sTo model interactions within a additive smooth model we need two or more dimensional smoothers. The problem of course with more that two dimensional smoothers is that thery can not be visualised. Here we stick to two dimensional smoothers that is, two way interactions;
gm21 <- gamlss2(rent~s(area,yearc)+location+bath+kitchen|
s(area,yearc)+location+bath+kitchen,
data=rent99, family=GA)GAMLSS-RS iteration 1: Global Deviance = 38552.6223 eps = 0.122820
GAMLSS-RS iteration 2: Global Deviance = 38467.4636 eps = 0.002208
GAMLSS-RS iteration 3: Global Deviance = 38466.85 eps = 0.000015
GAMLSS-RS iteration 4: Global Deviance = 38466.8343 eps = 0.000000 The interaction plots can be visualised using;
plot(gm21)
One could check for the best model using AIC.
AIC(gm11, gm21) AIC df
gm11 38522.57 25.69202
gm21 38523.11 28.13821Note that the default value for m in s(), that is, the dimension of the basis for one dimensional smoother, is 10 which could be very small if a lot of degrees of freedom are needed for the smoother. This is important in the constructing of centile curves, see vignette Centile (Quantile) Estimation.
Mikis I think for two dimensional smoothers it uses \(m=5\) but I am not sure.
pb() functionThe pb() function uses P-splines, see Eilers and Marx (2021). It tries to imitate the equivalent function pb() in the older gamlss package but still uses the s() function with different options. The option are s(..., bs="ps", k=20) that is the number of basis is increased from 10 to 20. The way to estimate the hyper-parameter (smoothing-parameter) is also change, It uses criterion="ml", which amount to a REML estimation of the smoothing parameter see Rigby and Stasinopoulos (2013). The result should be similar but not necessarily identical to the pb() function of the gamlss package.
gm1 <- gamlss2(rent~pb(area)+pb(yearc)+location+bath+kitchen|
pb(area)+pb(yearc)+ location+bath+kitchen,
data=rent99, family=GA)GAMLSS-RS iteration 1: Global Deviance = 38439.2661 eps = 0.125399
GAMLSS-RS iteration 2: Global Deviance = 38429.5701 eps = 0.000252
GAMLSS-RS iteration 3: Global Deviance = 38429.4305 eps = 0.000003 plot(gm1)
The smooth functions for \(\mu\) using the pv() function are almost identical to the ones when s() was used. The smooth terms for the \(\sigma\) model especially for area are different. This could be to the fact that more degrees of freedom are allowed with pb() and possible this lead to overfit.
Note that there are no two dimensional smoothers with pb()
Mikis It looks that pb() in gamlss2 overfits \(\sigma\). Here is what we get in gamnlss
g1 <- gamlss(rent~pb(area)+pb(yearc)+location+bath+kitchen,~pb(area)+pb(yearc)+
location+bath+kitchen,data=rent99, family=GA)GAMLSS-RS iteration 1: Global Deviance = 38459.65
GAMLSS-RS iteration 2: Global Deviance = 38448.8
GAMLSS-RS iteration 3: Global Deviance = 38448.64
GAMLSS-RS iteration 4: Global Deviance = 38448.62
GAMLSS-RS iteration 5: Global Deviance = 38448.62 term.plot(g1, "sigma", pages=1, term=c(1,2), ask=FALSE)