# Worm plot function
# Author: Stef van Buuren, TNO Prevention and Health, 2001
# For S-plus 5 and up
#
# See: van Buuren S, Fredriks M. Worm plot: a simple diagnostic device for modelling growth reference curves.
#      Statistics in Medicine 2001, 20, 1259-1277.

wp <- function(data, layout = c(4, 4), overlap = 0, worms = T, cubines = F, coefsave = F, labels = T, hor=T, vert=F, ci = T, sub = paste(deparse(substitute(
	data)), deparse(substitute(overlap))))
{
	panel <- function(x, y)
	{
		qq <- as.data.frame(qqnorm(y, plot = F))
		qq$y <- qq$y - qq$x
		plot(qq$x, qq$y, type = "n", ylim = c(-0.5, 0.5), xlim = c(-3, 3), lab = c(3, 5, 7), tck = -0.01)
		if (hor) abline(0, 0, lty = 2, col = 1)
		if (vert)	abline(0, 100000, lty = 2, col = 1)
		if(worms)
			points(qq$x, qq$y, col = 1, pch = 1, mkh = 0, cex = 0.25)
		if(cubines | coefsave)
			fit <- lm(y ~ x + x^2 + x^3, data = qq)
		if(cubines) {
			s <- spline(qq$x, fitted(fit))
			flags <- s$x > -2 & s$x < 2
			lines(list(x = s$x[flags], y = s$y[flags]))
		}
		if(coefsave) {
			est <- coef(fit)
			assign(".est", c(.est, est), frame = 0)
		}
		if (ci) ciplot(sum(!is.na(qq$y)))
	}
	agetext <- function(classes, layout = c(4,4), cex=0.6, dx = 0.06, dy=0.02)
	# function for adding age group text to the worm plot panels
	{
		txt <- apply(format(round(summary(classes)$intervals,1)),1,paste,collapse="-")
		x <- rep((0:(layout[1]-1))/layout[1]+dx,layout[2])
		y <- rep((1:layout[2])/layout[2]-dy,each=layout[1])
		text(x, y, txt, cex=cex)
	}
	assign("panel", panel, frame = 1)
	assign("worms", worms, frame = 1)
	assign("cubines", cubines, frame = 1)
	assign("coefsave", coefsave, frame = 1)
	assign("hor", hor, frame = 1)
	assign("vert", vert, frame = 1)
	assign("ci", ci, frame = 1)
	assign(".est", NULL, frame = 0)
	if(length(layout) == 1)
		layout <- rep(layout, 2)
	n <- prod(layout)
	classes <- equal.count(data$lft, n, overlap = overlap)
	if(n == 1)
		form <-  ~ data$z
	else form <-  ~ data$z | classes
	print.trellis(qqmath(form, layout = layout, aspect = 1, strip = F, sub = list(sub, cex = 0.5), xlab = list(
		"Unit normal quantile", cex = 0.75), ylab = list("Deviation", cex = 0.75), panel = panel))
	if (labels) agetext(classes, layout, cex=0.6, dx = 0.06, dy = 0.02)
	return(list(classes = classes, .est = get(".est", frame = 0)))
}

ciplot <- function(n, level=0.95, lz=-2.75, hz=2.75, dz=0.25){
	# 	adds confidence interval to Q-Q plot panel
		z <- seq(lz,hz,dz)
		p <- pnorm(z)	
		se <- (1/dnorm(z))*(sqrt(p*(1-p)/n))	
		low <- qnorm((1-level)/2)*se
		high <- -low
		lines(z, low, lty=2)
		lines(z, high, lty=2)
}

# ---- test it

# reproduce figure 
z <- wp(lnj0106r)

# reproduce analysis results in Table III for height for age (boys)
z <- wp(lnj0106r,coefsave=T,overlap=0.5)
matrix(.est,nrow=16,byrow=T,dimnames=list(NULL,c("b0","b1","b2","b3")))