Centile estimation includes methods for estimating the age related distribution of human growth.

The standard estimation of centile curves involves two continuous variables:

1. the response variable, that is, the variable we are interested in and for which we are trying to find the centile curves, e.g. weight, BMI, head circumference etc. and
2. the explanatory variable age.

The 100p centile of a random variable Y is the value y such that p(Y > y)=p, i.e. y= inv.cdf(p), where inv.cdf() the the inverse cumulative distribution function of  of Y applied to p. Here we consider the conditional centile of Y given explanatory variable x  (usually  age).  By varying x a 100p centile curve of y(x) against x is obtained. Centile curves can be obtained for different values of p. The World Health Organisation uses 100p=(3, 15, 50, 85, 97) in its charts and 100p=(1, 3, 5, 15, 25, 50, 75, 85, 95, 97, 99) in its tables.
This can be extended to more than one explanatory variable.

The methodology for creating growth centile references for individuals from a population comprises two different methods:

1. the non parametric method of quantile regression (Koenker, 2005; Koenker and Bassett, 1978, Koenker and Ng (2005), He and Ng (1999) and Ng and Maechler (2007))
2. the parametric LMS (i.e. Lambda, Mu and Sigma) method of Cole (1988), Cole and Green (1992)  and its extensions for example see Wright and Royston (1997),  van Buuren and Fredriks (2001),  and Rigby and Stasinopoulos (2004, 2006).

Here we are dealing with the LMS method and its extensions. The LMS method, within GAMLSS, is equivalent of assuming the Box- Cox Cole and Green distribution (BCCG) for the response variable  and fitting a smooth curves for μ, σ, and ν.  The BCCG distribution is derived by assuming that  Y, the response variable is a specific function  of a random variable Z which  has a (truncated) normal distribution. The BCCG distribution is suitable for positively or negatively skew data depending on the values of the parameter  ν.

Rigby and Stasinopoulos (2004, 2006) extended the LMS method (which allows for skewness and but not for kurtosis in the data), by introducing the Box-Cox power exponential (BCPE) and the Box-Cox t (BCT) distributions respectively and called the resulting methods LMSP and LMST respectively. The BCPE assumes that the transformed random variable Z has a (truncated) exponential power distribution while BCT assumes that Z has a (truncated) t distribution.

More recently the function lms() is introduced for fitting centile curves in gamlss package.

Two books on GAMLSS are in preparation:

1.  The Distribution Toolbox of GAMLSS
2.  GAMLSS Flexible Regression in R.

Draft versions of the two books will available in the web soon.

The gamlss packages reference card is available here

Each packages has its individual help files. Here is a list of what users may find useful.

• The second edition of the manual of the gamlss package in pdf form. This is a good starting point covering most of the topics of the original gamlss package but it is now dated since was created in 2008.
• The Journal of Statistical Software which has a brief introduction to GAMLSS and shows how the models can be used in practice.
• Course notes from the Utrecht 2008 short course on GAMLSS.
• Course notes from the Lancaster 2009 short course.
• The notes for the Athens 2010 short course can be downloaded here (12 MB) 

In the GAMLSS implementation in R, the function gamlss() allows modelling all the distribution parameters μ, σ, ν and τ as linear and/or non-linear and/or ‘non- parametric’ smoothing functions of the explanatory variables. This allow the explanatory variables to effect the predictors, (the η’s), of the specific parameters and therefore the parameters themselves. As a result the shape of the distribution of the response variable, (not only the mean), is effected by the explanatory variables.

All the standard linear terms as used in the lm() and glm() functions in R can be used here. In addition the following smoothing additive term functions can be used:

• pb(),  pvc() and cy(): based on P-splines,
• cs() and  scs(): based on cubic splines,
• fp(): fractional polynomials
• fk():  free knot smoothing (break points)
• lo():  local regression based on the loess() R function
• nn(): neural network based on the nnet() R function
• nl(), non-linear term fitting based on the nlm() R function
• random() : simple random effect
• ri(), ridge(): for ridge regression
• ga(): an iterface for the gam() function of Simon Wood in package mgcv

New  additive terms can be added relatively easy to the gamlss() function.

The R package gamlss.dist contains more than 70 distributions. We are refer to those distribution as “gamlss.family” distribution a name also given to the equivalent R objects. Each of those “gamlss.family” distributions has five related functions:

• the probability density function (d)
• the cumulative distribution function (p)
• the inverse of the cumulative distribution function and (q)
• a function which generates random numbers from the distribution (r) and
• the fitting gamlss function (what the gamlss() function uses for fitting)

1. continuous distributions
2. discrete distributions
3. mixed distributions

All the distributions in gamlss.family are shown in this distribution table. In addition any gamlss.family distribution can be:

• truncated using the function gen.trun() of the package gamlss.tr
• censored using the function gen.cens() of the package gamlss.cens
• finite mixed using the package gamlss.mx

Also any continuous distribution defined on real line can be log or logit transformed.

The ACEGES-based scenarios, such as oil and gas scenarios, are analysed within the GAMLSS framework.  GAMLSS is also used for the statistical modelling of the agents of the ACEGES model.The statistical modelling of the agents is in addition to the heuristic modelling of the agents. For a theoretical discussion of the integration of ACEGES and GAMLSS, see here,

The diagnostics for GAMLSS models are based on the residuals of the fitted model.The GAMLSS models use the  normalised quantile residuals for continuous response variables and randomised normalised quantile residuals for discrete response variables.

The main advantage of the normalised (randomised) quantile residuals is that, whatever the distribution of the response variable their true values r always have a standard normal distribution given that the assumptions the model is correct. Since within the statistical literature checking the normality assumption is well established the normalised (randomised) quantile residuals provide us an easy way to check the adequacy of a GAMLSS fitted model.

There are several functions in the R implementation of GAMLSS using the residuals.

• The plot() function displaying the residuals 
• The worm plot function wp()
• The detrended Own’s plot function dtop() 
• The Q-statistics function Q.stats()
• The randomised residual plot function rqres.plot()

package gamlss:

• The functions prof.dev() and prof.term() are improved. The argument step is not anymore compulsory and if not set the argument length is used instead. For most cases there is no need to have a fine grid since the function is approximated using splinefun(). The output is saved as an ”profDeviance.gamlss” object. Note that more testing is needed for the reliability of the function.
• The function fitDist() allows extra arguments “…” to be passed to gamlssML() and gamlss().
• The prof.dev() can be used now in conjunction with gamlssML()
• The function Q.stats() now allows plotting the resulting matrix for easy identification of the parts of the model which do not fit well.
• The cs() and scs() are calling now the R function smooth.spline() rather than the FORTRAN code to comply with R regulations
• The function vc() is disfuction the user is advised to use the equivalent function pvc()
• The function lo() is rewritten to comply with R regulation. Now it takes a formula as its first argument rather than a list of explanatory variables. Also no standard errors for the smooth function are provided since the R function loess()  do not provide this information at the moment.
• The logic link function in the package gamlss.dist is amended so it does not call the R .C function.
• summary.gamlss() now have an argument “save” for saving the output, thanks to  Wilmar Igl
• gamlssML(): a bug with vcov.gamlssML() function is fixed also “nlminb” is now the default maximisation procedure rather than “optim”
• The default values of the argument cent on lms() and calibrartion() is change to 100*pnorm((-4:4)*2/3) as suggested by Tim Cole. Also a bug which did not allow term.plot() to work with lms() is fixed.

• All the FORTRAN routines have their REAL change to “double precession”
• The distribution function SICHEL has been amended to so in the limited case where sigma is bigger than 10000 and nu > zero it switches to the negative binomial.  (the old version has sigma>1000 which creates problem with the prof.dev() function when was used for nu in the lice data see the Stasinopoulos and Rigby (2007) paper in JSS. Thanks to Ivailo Stoyanov who point out the problem.)
• The skew Normal 1, SN1, skew Normal 2, SN2,  distributions are introduced
• The link function “logshiftto2″ is added in make.link.gamlss() in package gamlss.dist. The reason for this is to prevent  the degrees of freedom parameter nu in TF2 (see below) to be less than 2.
• The distributions SST and TF2 are introduced in package gamlss.dist. The distributions are reparametrisations of the ST3 and TF respectively. The sigma parameter in SST and TF2 is the standard deviation of the distribution. Not that the standard deviation is not defined for degrees of freedom less than 2. The “logshiftto2″ link function (see above) prevents this.
• The logit normal, LOGITNO and the log normal 2 (with mu as the median), LOGNO2, are introduced in gamlss.dist.
• The functions gen.Family(), Family(), Family.d(), Family.p(), Family.q() and  Family.p() are introduced in gamlss.dist for generating “log” and “logic” versions of continuous  gamlss.family distributions in the real line.
• The generalised Pareto (GP) distribution (a re-parameterisation of PARETO2 and PARETO2o) is introduced in gamlss.dist
• A bug in the GB2 and  EGB2 q functions is fixed.

• The functions fitFixBP() and fitFreeKnots() are moved to package gamlss.add
• garmaFit() for fitting generalised ARMA models is introduced
• centile.ts() for giving centiles for time series data is introduced
• lagPlot() for scatter plotting of lags is introduced.