gamlss (version 4.2-4)

prof.term: Plotting the Profile: deviance or information criterion for one of the terms (or hyper-parameters) in a GAMLSS model

Description

This functions plots the profile deviance for a chosen parameter included in the linear predictor of any of the mu,sigma, nu or tau models so profile confidence intervals can be obtained. In can also be used to plot the profile of a specified information criterion for any hyper-parameter when smooth additive terms are used.

Usage

prof.term(model = NULL, criterion = c("GD", "GAIC"), penalty = 2.5, 
          other = NULL, min = NULL, max = NULL, step = NULL, 
          length = 7, xlabel = NULL, plot = TRUE, perc = 95, 
          start.prev = TRUE, ...)

Arguments

model

this is a GAMLSS model, e.g. model=gamlss(y~cs(x,df=this),sigma.fo=~cs(x,df=3),data=abdom), where this indicates the (hyper)parameter to be profiled

criterion

whether global deviance ("GD") or information criterion ("GAIC") is profiled. The default is global deviance criterion="GD"

penalty

The penalty value if information criterion is used in criterion, default penalty=2.5

other

this can be used to evaluate an expression before the actual fitting of the model (Make sure that those expressions are well define in the global environment)

min

the minimum value for the parameter e.g. min=1

max

the maximum value for the parameter e.g. max=20

step

how often to evaluate the global deviance (defines the step length of the grid for the parameter) e.g. step=1

length

if the step is left NULL then length is considered for evaluating the grid for the parameter. It has a default value of 11

xlabel

if a label for the axis is required

plot

whether to plot, plot=TRUE the resulting profile deviance (or GAIC)

perc

what % confidence interval is required

start.prev

whether to start from the previous fitted model parameters values or not (default is TRUE)

...

for extra arguments

Value

• Return a profile plot (if the argument plot=TRUE) and an ProfLikelihood.gamlss object if saved. The object contains:
• values
  the values at the grid where the parameter was evaluated
• fun
  the function which approximates the points using splines
• min
  the minimum values in the grid
• max
  the maximum values in the grid
• max.value
  the value of the parameter maximising the Profile deviance (or GAIC)
• CI
  the profile confidence interval (if global deviance is used)
• criterion
  which criterion was used

Warning

A dense grid (i.e. small step) evaluation of the global deviance can take a long time, so start with a sparse grid (i.e. large step) and decrease gradually the step length for more accuracy.

Details

This function can be use to provide likelihood based confidence intervals for a parameter involved in terms in the linear predictor(s). These confidence intervals are more accurate than the ones obtained from the parameters' standard errors. The function can also be used to plot a profile information criterion (with a given penalty) against a hyper-parameter. This can be used to check the uniqueness in hyper-parameter determination using for example find.df.

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.

Stasinopoulos D. M., Rigby R.A. and Akantziliotou C. (2006) Instructions on how to use the GAMLSS package in R. Accompanying documentation in the current GAMLSS help files, (see also http://www.gamlss.org/).

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, http://www.jstatsoft.org/v23/i07.

See Also

gamlss, prof.dev

Examples

data(aids)
# fitting a linear model
gamlss(y~x+qrt,family=NBI,data=aids)
# testing the linear beta parameter
mod<-quote(gamlss(y ~ offset(this * x) + qrt, data = aids, family = NBI))
prof.term(mod, min=0.06, max=0.11)
# find the hyper parameter using cubic splines smoothing
mod1<-quote(gamlss(y ~ cs(x,df=this) + qrt, data = aids, family = NBI))
prof.term(mod1, min=1, max=15, step=1, criterion="GAIC", penalty=log(45))
# find a break point in x
mod2 <- quote(gamlss(y ~ x+I((x>this)*(x-this))+qrt,family=NBI,data=aids))
prof.term(mod2, min=1, max=45, step=1, criterion="GD")
rm(mod,mod1,mod2)