gamlss (version 5.4-12)

wp: Worm plot

Description

Provides a single plot or multiple worm plots for a GAMLSS fitted or more general for any fitted models where the method resid() exist and the residuals are defined sensibly. The worm plot (a de-trended QQ-plot), van Buuren and Fredriks M. (2001), is a diagnostic tool for checking the residuals within different ranges (by default not overlapping) of the explanatory variable(s).

Usage

wp(object = NULL, xvar = NULL, resid = NULL, n.inter = 4,
xcut.points = NULL, overlap = 0, xlim.all = 4, 
xlim.worm = 3.5, show.given = TRUE, line = TRUE, 
ylim.all = 12 * sqrt(1/length(resid)), 
ylim.worm = 12 * sqrt(n.inter/length(resid)), 
cex = 1, cex.lab = 1, pch = 21, bg = "wheat", 
col = "red", bar.bg = c(num = "light blue"), ...)

Value

For multiple plots the xvar intervals and the coefficients of the fitted cubic polynomials to the residuals (within each xvar interval) are returned.

Arguments

object

a GAMLSS fitted object or any other fitted model where the resid() method works (preferably it should be standardised or quantile residuals)

xvar

the explanatory variable(s) against which the worm plots will be plotted. If only one variable is involved use xvar=x1 if two variables are involved use xvar=~x1*x2. See also note below for use of formula if the data argument is not found in the fitted model

resid

if object is missing this argument can be used to specify the residual vector (again it should a quantile residuals or it be assumed to come from a normal distribution)

n.inter

the number of intervals in which the explanatory variable xvar will be cut

xcut.points

the x-axis cut off points e.g. c(20,30). If xcut.points=NULL then the n.inter argument is activated

overlap

how much overlapping in the xvar intervals. Default value is overlap=0 for non overlapping intervals

xlim.all

for the single plot, this value is the x-variable limit, default is xlim.all=4

xlim.worm

for multiple plots, this value is the x-variable limit, default is xlim.worm=3.5

show.given

whether to show the x-variable intervals in the top of the graph, default is show.given=TRUE

line

whether to plot the polynomial line in the worm plot, default value is line=TRUE

ylim.all

for the single plot, this value is the y-variable limit, default value is ylim.all=12*sqrt(1/length(fitted(object)))

ylim.worm

for multiple plots, this values is the y-variable limit, default value is ylim.worm=12*sqrt(n.inter/length(fitted(object)))

cex

the cex plotting parameter for changing the side of worm with default cex=1

cex.lab

the cex plotting parameter for changing the size of the axis labels

pch

the pch plotting parameter with default pch=21

bg

The background colour of the worm plot points

col

the colour of the fitted (and horizontal and vertical) lines

bar.bg

the colour of the bars when xvar is used

...

for extra arguments

Author

Mikis Stasinopoulos and Bob Rigby

Details

If the xvar argument is not specified then a single worm plot is used. In this case a worm plot is a de-trended normal QQ-plot so departure from normality is highlighted.

If a single xvar is specified (with or without the use of a formula) i.e. xvar=x1 or xvar=~x1) then we have as many worm plot as n.iter. In this case the x-variable is cut into n.iter intervals with an equal number observations and de-trended normal QQ (i.e. worm) plots for each interval are plotted. This is a way of highlighting failures of the model within different ranges of the the single explanatory variable. The fitted coefficients from fitting cubic polynomials to the residuals (within each x-variable interval) can be obtain by e.g. coeffs<-wp(model1,xvar=x,n.iner=9). van Buuren and Fredriks M. (2001) used these residuals to identify regions (intervals) of the explanatory variable within which the model does not fit adequately the data (called "model violation")

Two variables can be displayed with the use of a formula, i.e. xvar=~x1*x2. In this case the n.inter can be a vector with two values.

References

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

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/).

van Buuren and Fredriks M. (2001) Worm plot: simple diagnostic device for modelling growth reference curves. Statistics in Medicine, 20, 1259--1277

See Also

gamlss, plot.gamlss

Examples

data(abdom)
# with data
a<-gamlss(y~pb(x),sigma.fo=~pb(x,1),family=LO,data=abdom)
wp(a)
coeff1<-wp(a,xvar=x)
coeff1
if (FALSE) {
# no data argument
b <- gamlss(abdom$y~pb(abdom$x),sigma.fo=~pb(abdom$x),family=LO)
wp(b) 
wp(b, xvar=abdom$x)# not wp(b, xvar=x)
# using  the argument resid
# this will work
wp(resid=resid(a),  xvar=abdom$x)
# not this
# wp(resid=resid(a),  xvar=x)
# this example uses the rent data
m1 <- gamlss(R~pb(Fl)+pb(A)+loc, sigma.fo=~pb(Fl)+pb(A), data=rent, family=GA)
# a single worm plot
wp(m1, ylim.all=0.5)
# a single continuous x variable 
wp(m1, xvar=Fl, ylim.worm=.8)
# a single x variable changing the default number of intervals
wp(m1, xvar=Fl, ylim.worm=1.5, n.inter=9)
# different x variable changing the default number of intervals
B1<-wp(m1, xvar=A, ylim.worm=1.2, n.inter=9) 
B1
# the number five plot has intervals 
# [5,] 1957.5 1957.5 
# rather disappoining 
# try formula for xvar
wp(m1, xvar=~A, ylim.worm=1.2, n.inter=9)
# better in this case using formula
# now using a factor included in the model
wp(m1, xvar=~loc, ylim.worm=1.2, n.inter=9)
# using a factor notin the model
wp(m1, xvar=~B, ylim.worm=1.5, n.inter=9)
# level 2 (with B=1) did not fit well
# trying two continuous variable 
wp(m1, xvar=~Fl*A, ylim.worm=1.5, n.inter=4)
# one continuous and one categorical 
wp(m1, xvar=~Fl*loc, ylim.worm=1.5, n.inter=4)
# two categorical
wp(m1, xvar=~B*loc, ylim.worm=1.5, n.inter=4)
}