find.hyper {gamlss}

A function to select values of hyper-parameters in a GAMLSS model

Description

This function selects the values of hyper parameters and/or non-linear parameters in a GAMLSS model. It uses the R function optim which then minimises the generalised Akaike information criterion (GAIC) with a user defined penalty.

Usage

find.hyper(model = NULL, parameters = NULL, other = NULL, k = 2,
        steps = c(0.1), lower = -Inf, upper = Inf, method = "L-BFGS-B",
        ...)

Arguments

model

this is a GAMLSS model in quote(). e.g.
quate(gamlss(y~cs(x,df=p[1]),sigma.fo=~cs(x,df=p[2]),data=abdom)) where p[1] and p[2] denote the parameters to be estimated

parameters

the starting values in the search of the optimum hyper-parameters and/or non-linear parameters e.g. parameters=c(3,3)

other

this is used to optimise other non-parameters, for example a transformation of the explanatory variable of the kind x^{p[3]}, others=quote(nx<-x^p[3]) where nx is now in the model formula

k

specifies the penalty in the GAIC, (the default is 2) e.g. k=3

steps

the steps taken in the optimisation procedure [see the ndeps option in optim()], by default is set to 0.1 for all hyper parameters and non-linear parameters

lower

the lower permissible level of the parameters i.e. lower=c(1,1) this does not apply if a method other than the default method "L-BFGS-B" is used

upper

the upper permissible level of the parameters i.e. upper=c(30,10), this is not apply if a method other than the default method "L-BFGS-B" is used

method

the method used in optim() to numerically minimise the GAIC over the hyper-parameters and/or non-linear parameters. By default this is "L-BFGS-B" to allow box-restriction on the parameters

...

for extra arguments to be passed to the R function optim() used in the optimisation

Details

This historically was an experimental function which worked well for the search of the optimum degrees of freedom and non-linear parameters (e.g. power parameter \lambda used to transform x to x^\lambda). With the introduction of the P-Spline smoothing function pb() the function find.hyper() became almost redundant. find.hyper() takes lot longer than pb() to find automatically the hyper parameters while both method produce similar results. See below the examples for a small demonstration.

Value

The function turns the same output as the function optim()

par

the optimum hyper-parameter values

value

the minimised value of the GAIC

counts

A two-element integer vector giving the number of calls to ‘fn’ and ‘gr’ respectively

convergence

An integer code. ‘0’ indicates successful convergence. see the function optim() for other errors

message

A character string giving any additional information returned by the optimiser, or ‘NULL’

Warning

It may be slow to find the optimum

Author(s)

Mikis Stasinopoulos

References

Rigby, R. A. and Stasinopoulos D. M. (2005). Generalized additive models for location, scale and shape,(with discussion), Appl. Statist., 54, part 3, pp 507-554.

Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC. An older version can be found in https://www.gamlss.com/.

Stasinopoulos D. M. Rigby R.A. (2007) Generalized additive models for location scale and shape (GAMLSS) in R. Journal of Statistical Software, Vol. 23, Issue 7, Dec 2007, https://www.jstatsoft.org/v23/i07/.

Stasinopoulos D. M., Rigby R.A., Heller G., Voudouris V., and De Bastiani F., (2017) Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC.

(see also https://www.gamlss.com/).

See Also

gamlss, plot.gamlss, optim

Examples

## Not run: 
data(abdom)
# Example estimating the smoothing parameters for mu and 
# the transformation parameters for x
# declare the model
mod1<-quote(gamlss(y~cs(nx,df=p[1]),family=BCT,data=abdom,
                        control=gamlss.control(trace=FALSE)))
# since we want also to find the transformation for x 
# we use the "other"" option
op <- find.hyper(model=mod1, other=quote(nx<-x^p[2]), parameters=c(3,0.5), 
            lower=c(1,0.001), steps=c(0.1,0.001))
op
# the optimum parameters found are 
# p = (p[1],p[2]) = (3.113218 0.001000) = (df for mu, lambda)
# so it needs df = 3 on top of the constant and linear 
# in  the cubic spline model for mu since p[1] is approximately  3
#  and log transformation for x since p[2] is approximately  0 
# here is an example with no data declaration in define the model
# we have to attach the data
attach(abdom)
mod2 <- quote(gamlss(y~cs(nx,df=p[1]),family=BCT,
                control=gamlss.control(trace=FALSE)))
op2<-find.hyper(model=mod2, other=quote(nx<-x^p[2]), parameters=c(3,0.5), 
                lower=c(1,0.001), steps=c(0.1,0.001))
op2
detach(abdom)
#--------------------------------------------------------------
# showing different ways of estimating the smoothing parameter
# get the df using local ML (PQL) 
m0 <- gamlss(y~pb(x), data=abdom)
# get the df using local GAIC 
m1<-gamlss(y~pb(x, method="GAIC", k=2), data=abdom)
# fiiting cubic splines with fixed df's at 3
m2<-gamlss(y~cs(x, df=3), data=abdom)
# fitting cubic splines using find hyper (global GAIC)
mod1 <- quote(gamlss(y~cs(x, df=p[1]),family=BCT,data=abdom,control=gamlss.control(trace=FALSE)))
op <- find.hyper(model=mod1, parameters=c(3), lower=c(1,0.001), steps=c(0.1,0.001))
# now fit final model
m3 <- gamlss(y~cs(x, df=op$par), data=abdom)
# effetive degrees of fredom for the 4 models
edf(m0);edf(m1); m2$mu.df; m3$mu.df
# deviances for the four models
deviance(m0); deviance(m1); deviance(m2); deviance(m3)
# their GAIC
GAIC(m0,m1,m2,m3)
# plotting  the models
plot(y~x, data=abdom, type="n")
lines(fitted(m3)~abdom$x, col="red")
lines(fitted(m1)~abdom$x, col="green")
lines(fitted(m0)~abdom$x, col="blue")
# almost identical

## End(Not run)