### A graph for presenting individual growth with time in a population with:
### Histograms of size distribution with time,
### Quantile regressions for the median and 5% / 95% percentiles,
### A graph quantile residuals using boxplots around the median model.
###
### This graph is Fig. 2 of the paper: Grosjean, PH, Ch. Spirlet & M. Jangoux,
### 2003. A functional growth model with intraspecific competition applied to a
### sea urchin, Paracentrotus lividus. Can. J. Fish. Aquat. Sci., 60:237-246.

### Requires two libraries and a couple of custom functions + an example dataset
library(stats) 
library(quantreg) 

##### Custom functions required ####################################
# Fuzzy-remanent growth curve with two different kinetic parameters
# passing by {0, 0}
# See Grosjean et al, 2003, Can. J. Fish. Aquat. Sci. 60:237-246,
# for definition and discussion of this growth model
fuzremOrig2 <- function (input, Asym, lrc1, lrc2, c0) {
    .expr1 <- exp(lrc1) 
    .expr4 <- exp(-.expr1 * input) 
    .expr5 <- 1 - .expr4
    .expr6 <- Asym * .expr5
    .expr7 <- exp(lrc2) 
    .expr9 <- input - c0
    .expr11 <- exp(-.expr7 * .expr9) 
    .expr12 <- 1 + .expr11
    .expr22 <- .expr12^2
    .value <- .expr6/.expr12
    .actualArgs <- as.list(match.call()[c("Asym", "lrc1", "lrc2", "c0")]) 
    if (all(unlist(lapply(.actualArgs, is.name)))) { 
        .grad <- array(0, c(length(.value), 4), 
            list(NULL, c("Asym", "lrc1", "lrc2", "c0"))) 
        .grad[, "Asym"] <- .expr5 / .expr12
        .grad[, "lrc1"] <- Asym * (.expr4 * (.expr1 * input)) / .expr12
        .grad[, "lrc2"] <- .expr6 * (.expr11 * (.expr7 * .expr9)) / .expr22
        .grad[, "c0"] <- -.expr6 * (.expr11 * .expr7) / .expr22
        dimnames(.grad) <- list(NULL, .actualArgs) 
        attr(.value, "gradient") <- .grad 
    }
    .value
}
# Self-starting function for fuzremOrig2
fuzremOrig2.ival <- function (mCall, data, LHS) { 
    xy <- sortedXyData(mCall[["input"]], LHS, data) 
    if (nrow(xy) < 5) { 
        stop("Too few distinct input values to fit a fuzzy-remanent growth function with 2 kinetic parameter passing by {0, 0}") 
    }
    xydata <- c(as.list(xy), c0 = NLSstClosestX(xy, mean(range(xy[["y"]])))) 
    xydata <- as.list(xydata) 
    options(show.error.messages = FALSE) 
    # Sometimes, the following evaluation gives an error
    # (singular gradient, ...). Try another way to estimate initial parameters
    pars <- try(as.vector(coef(nls(y ~ (1 - exp(-exp(lrc) * x))/ 
        (1 + exp((c0 - x)*exp(lrc))), data = xydata, start = list(lrc = 0), 
        algorithm = "plinear"))))
    if (!is.null(class(pars)) && class(pars) == "try-error") { 
        cat("Using second alternative for initial parameters estimation\n") 
        options(show.error.messages = TRUE) 
        pars <- as.vector(coef(nls(y ~ SSlogis(x, Asym, xmid, scal), 
            data = xydata))) 
        xydata$c0 <- pars[2]
        pars[1] <- log(1/pars[3]) 
    } else { options(show.error.messages = TRUE)} 
    pars <- as.vector(coef(nls(y ~ (1 - exp(-exp(lrc1) * x))/ 
        (1 + exp((c0 - x)*exp(lrc2))), data = data.frame(xy), 
        start = list(c0 = xydata$c0, lrc1 = pars[1], lrc2 = pars[1]), 
        algorithm = "plinear")))
    value <- c(pars[4], pars[2], pars[3], pars[1])
    names(value) <- mCall[c("Asym", "lrc1", "lrc2", "c0")] 
    return(value) 
}
# Self-Starting function, combining fuzremOrig2 and fuzremOrig2.ival
SSfuzremOrig2 <- selfStart(fuzremOrig2, fuzremOrig2.ival, para = "", template = "") 
remove(fuzremOrig2, fuzremOrig2.ival) 
attr(SSfuzremOrig2,"pnames") <- c("Asym", "lrc1", "lrc2", "c0") 

# Expand a "class"-"frequencies at time" data frame into
# a "time"-"mean size in class" data frame
expandFreq <- function(ages, freqs) {
	for (i in 2:ncol(freqs)){ 
		size <- rep(freqs[, 1], freqs[, i]) 
		age <- rep(ages[i-1], length(size)) 
		if (i == 2){
			res <- cbind(age, size) 
		} else {
			res <- rbind(res, cbind(age, size)) 
		}
	}
	res <- as.data.frame(res) 
	return(res) 
}

# Plot of size distributions, quantiles regressions and residuals diagnosis
# (boxplot) for curves tau = 0.05, 0.5 & 0.95 in the case of individual size
# distributions are known at each sampled time
rqFreqPlot <- function(ages, freqs, agecurves, curves, xlim = c(min(ages), 
    max(ages)), ylim = c(min(curves), max(curves)), barscale = 100, barcol = 8, 
    boxwex = 50, ylab1 = "", ylab2 = "", lty = c(2, 1, 2), ...) {
 	# Verify parameters ages and agecurves
	if (!all(ages %in% agecurves)) 
		stop("ages must be a subset of agecurves!") 

	# Create a layout for the histograms and the boxplots
	nf <- layout(matrix(c(2,1),2,1,byrow = TRUE), widths = 7, 
        heights = c(2, 5), respect = TRUE)
	layout.show(nf) 
	par(mar = c(5, 4, 0, 2), lab = c(5, 5, 7)) 
	# Draw the curves
	par(new = FALSE) 
	X <- matrix(rep(agecurves, ncol(curves)), nrow = length(agecurves), 
        ncol = ncol(curves)) 
	Y <- curves
	matplot(X, Y, type = "l", xaxs = "i", lty = lty, col = 1, lwd = 2, 
        bty = "l", xlim = xlim, ylim = ylim, ylab = ylab1, ...)
	# Draw the histograms
	for (i in 1:length(ages)){ 
		par(new = TRUE) 
		xmin <- -ages[i] + xlim[1]
		xmax <- xlim[2] - ages[i]
		ser <- freqs[, i+1]
		ser <- ser/max(ser) * barscale 
		barplot(ser, horiz = TRUE, axes = FALSE, xlim = c(xmin, xmax), 
            ylim = ylim, col = barcol, space = 0) 
	}
	# Add a residual analysis above the graph
	spreadCalc <- function(x, na.rm = FALSE){ 
		# Substract median from expanded frequencies data
		# And narrow range to 0.05 - 0.95 quantiles
		res <- x
		q <- quantile(res, c(0.05, 0.95), na.rm = na.rm) 
		res[res < q[1]] <- q[1]
		res[res > q[2]] <- q[2]
		return(res) 
	}
	## Expand one or several frequency columns into individual measurements
	expandFreq2 <- function(freqs) {
		# Consider mean size at each class for expanding the observations
		# with classes being first column and frequencies in each class
        # in the following columns
		# A list of vectors is returned
		res <- list(NULL) 
		nCols <- ncol(freqs) 
		for (i in 2:nCols)
			res[[i-1]] <- rep(freqs[, 1], freqs[, i]) 
		return(res) 
	}

	# A modified boxplot.default() function with different parameters
    # (no box around the graph and plain lines for hinges)
	boxplot2 <- function (x, ..., range = 1.5, width = NULL, varwidth = FALSE, 
        notch = FALSE, outline = TRUE, names, boxwex = 0.8, plot = TRUE, 
        border = par("fg"), col = NULL, log = "", pars = NULL, 
        horizontal = FALSE, add = FALSE, at = NULL) {
	    args <- list(x, ...) 
	    namedargs <- if (!is.null(attributes(args)$names)) 
	        attributes(args)$names != "" 
	    else rep(FALSE, length = length(args)) 
	    pars <- c(args[namedargs], pars) 
	    groups <- if (is.list(x)) 
	        x
	    else args[!namedargs] 
	    if (0 == (n <- length(groups))) 
	        stop("invalid first argument") 
	    if (length(class(groups))) 
	        groups <- unclass(groups) 
	    if (!missing(names)) 
	        attr(groups, "names") <- names 
	    else {
	        if (is.null(attr(groups, "names"))) 
	            attr(groups, "names") <- 1:n 
	        names <- attr(groups, "names") 
	    }
	    for (i in 1:n) groups[i] <- list(boxplot.stats(groups[[i]], range)) 
	    stats <- matrix(0, nr = 5, nc = n) 
	    conf <- matrix(0, nr = 2, nc = n) 
	    ng <- out <- group <- numeric(0) 
	    ct <- 1
	    for (i in groups) {
	        stats[, ct] <- i$stats
	        conf[, ct] <- i$conf
	        ng <- c(ng, i$n)
	        if ((lo <- length(i$out))) { 
	            out <- c(out, i$out)
	            group <- c(group, rep(ct, lo)) 
	        }
	        ct <- ct + 1
	    }
	    z <- list(stats = stats, n = ng, conf = conf, out = out, group = group, names = names) 
	    if (plot) { 
	        bxp2(z, width, varwidth = varwidth, notch = notch, boxwex = boxwex, border = border, col = col, log = log, pars = pars, outline = outline, horizontal = horizontal, add = add, at = at) 
	        invisible(z) 
	    }
	    else z
	}

	# A modified bxp() function for plotting the boxplot generated by boxplot2()
	bxp2 <- function (z, notch = FALSE, width = NULL, varwidth = FALSE,
        outline = TRUE, notch.frac = 0.5, boxwex = 0.8, border = par("fg"), 
        col = NULL, log = "", pars = NULL, frame.plot = axes, 
        horizontal = FALSE, add = FALSE, at = NULL, show.names = NULL, ...) { 
	    pars <- c(pars, list(...)) 
	    bplt <- function(x, wid, stats, out, conf, notch, border,
	        col, horizontal) { 
	        if (!any(is.na(stats))) { 
	            wid <- wid/2
	            if (notch) {
	                xx <- x + wid * c(-1, 1, 1, notch.frac, 1, 1,
	                  -1, -1, -notch.frac, -1)
	                yy <- c(stats[c(2, 2)], conf[1], stats[3], conf[2],
	                  stats[c(4, 4)], conf[2], stats[3], conf[1])
	            }
	            else {
	                xx <- x + wid * c(-1, 1, 1, -1)
	                yy <- stats[c(2, 2, 4, 4)]
	            }
	            if (!notch)
	                notch.frac <- 1
	            wntch <- notch.frac * wid
	            if (horizontal) {
	                polygon(yy, xx, col = col, border = border) 
	                segments(stats[3], x - wntch, stats[3], x + wntch, 
	                  col = border) 
	                segments(stats[c(1, 5)], rep(x, 2), stats[c(2, 
	                  4)], rep(x, 2), lty = 1, col = border) 
	                segments(stats[c(1, 5)], rep(x - wid/2, 2), stats[c(1, 
	                  5)], rep(x + wid/2, 2), col = border) 
	                do.call("points", c(list(out, rep(x, length(out))), 
	                  pt.pars)) 
	            }
	            else {
	                polygon(xx, yy, col = col, border = border) 
	                segments(x - wntch, stats[3], x + wntch, stats[3], 
	                  col = border) 
	                segments(rep(x, 2), stats[c(1, 5)], rep(x, 2), 
	                  stats[c(2, 4)], lty = 1, col = border) 
	                segments(rep(x - wid/2, 2), stats[c(1, 5)], rep(x + 
	                  wid/2, 2), stats[c(1, 5)], col = border) 
	                do.call("points", c(list(rep(x, length(out)), 
	                  out), pt.pars)) 
	            }
	            if (any(inf <- !is.finite(out))) { 
	                warning(paste("Outlier (", paste(unique(out[inf]), 
	                  collapse = ", "), ") in ", paste(x, c("st", 
	                  "nd", "rd", "th")[pmin(4, x)], sep = ""), " boxplot are NOT drawn", 
	                  sep = ""))
	            }
	        }
	    }
	    if (!is.list(z) || 0 == (n <- length(z$n))) 
	        stop("invalid first argument") 
	    if (is.null(at)) 
	        at <- 1:n
	    else if (length(at) != n) 
	        stop(paste("`at' must have same length as `z $ n', i.e.", 
	            n))
	    if (is.null(z$out)) 
	        z$out <- numeric() 
	    if (is.null(z$group) || !outline) 
	        z$group <- integer() 
	    if (is.null(pars$ylim)) 
	        ylim <- range(z$stats[is.finite(z$stats)], z$out[is.finite(z$out)], 
	            if (notch)
	                z$conf[is.finite(z$conf)]) 
	    else {
	        ylim <- pars$ylim
	        pars$ylim <- NULL
	    }
	    width <- if (!is.null(width)) { 
	        if (length(width) != n | any(is.na(width)) | any(width <= 0)) 
	            stop("invalid boxplot widths") 
	        boxwex * width/max(width) 
	    }
	    else if (varwidth)
	        boxwex * sqrt(z$n/max(z$n)) 
	    else if (n == 1)
	        0.5 * boxwex
	    else rep(boxwex, n) 
	    if (missing(border) || length(border) == 0) 
	        border <- par("fg") 
	    pt.pars <- c(pars[names(pars) %in% c("pch", "cex", "bg")], 
	        col = border) 
	    for (i in 1:n) bplt(at[i], wid = width[i], stats = z$stats[,
	        i], out = z$out[z$group == i], conf = z$conf[, i], notch = notch,
	        border = border[(i - 1)%%length(border) + 1], col = if (is.null(col)) 
	            col 
	        else col[(i - 1)%%length(col) + 1], horizontal = horizontal) 
	    invisible(at) 
	}

	# Calculate and plot residuals
	baselevels <- curves[agecurves %in% ages, 2]
	expfreqs <- expandFreq2(freqs)
	for (i in 1:length(baselevels)) 
		expfreqs[[i]] <- expfreqs[[i]] - baselevels[i]
	# Draw it
	par(mar = c(1, 4, 1, 2), lab = c(5, 3, 7)) 
	matplot(X, Y - Y[, 2], type = "l", xaxs = "i", xaxt = "n", lty = lty, 
        col = 1, lwd = 2, bty = "l", xlim = xlim, ylim = c(-12, 20), 
        ylab = ylab2,...)
	# Some options cannot be changed in the original boxplot function!
    # So, we use our own implementation here
	boxplot2(lapply(expfreqs, spreadCalc), col = "gray90", boxwex = boxwex, 
        range = 0, add = TRUE, at = ages) 
	# Draw tick marks for the X axis
	axis(1, labels = FALSE) 
}


### The example dataset ########################################################
# A contingency table of observed frequencies by size classes and ages
freqs <- structure(list(MeanD = c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 
        8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,
        19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
        30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5,
        41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5,
        52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5,
        63.5, 64.5, 65.5, 66.5), F0.5y = c(65, 142, 114, 115, 65, 62,
        31, 34, 31, 20, 17, 18, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
        ), F1.0y = c(0, 11, 22, 19, 27, 32, 42, 46, 44, 38, 28, 22, 25,
        16, 16, 15, 23, 15, 13, 11, 5, 6, 10, 5, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), F1.5y = c(0, 0,
        0, 0, 0, 0, 0, 0, 2, 4, 6, 10, 22, 31, 30, 35, 29, 35, 25, 31,
        25, 14, 21, 27, 19, 17, 7, 8, 12, 3, 8, 7, 6, 8, 5, 7, 4, 3,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0), F2.0y = c(0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 2, 0, 3, 10, 6, 12, 20, 25, 26, 36, 36, 35, 10, 23,
        27, 18, 18, 16, 7, 4, 7, 7, 10, 9, 8, 7, 6, 11, 1, 2, 0, 1, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
        ), F2.5y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 1, 0, 0, 0, 2, 0, 5, 6, 15, 18, 15, 12, 23, 29, 22, 33, 23,
        17, 16, 19, 11, 15, 12, 13, 3, 16, 11, 12, 13, 10, 7, 3, 2, 1,
        2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), F3.0y = c(0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        1, 0, 0, 1, 0, 2, 1, 4, 9, 16, 9, 16, 9, 13, 17, 15, 13, 28,
        21, 22, 19, 14, 17, 12, 13, 10, 8, 11, 14, 8, 0, 0, 0, 1, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0), F3.5y = c(0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
        0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 5, 5, 9, 8, 19, 16, 17, 31, 28,
        24, 30, 22, 16, 15, 8, 13, 11, 12, 8, 4, 1, 2, 0, 0, 0, 0, 0,
        0, 0), F4.0y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 1, 0, 0, 1, 8, 3, 1, 15, 13, 15, 16, 9, 34, 22, 22, 16, 6,
        14, 13, 8, 7, 2, 4, 2, 0, 1, 0, 0, 0, 0), F4.5y = c(0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 6,
        3, 7, 8, 9, 12, 21, 25, 28, 13, 21, 17, 9, 9, 5, 11, 7, 4, 2,
        0, 2, 0, 0), F5.0y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 2, 3, 2, 2, 6, 10, 8, 10, 14, 16,
        19, 12, 9, 6, 5, 4, 0, 5, 1, 0, 0, 1, 0), F5.5y = c(0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
        2, 1, 4, 3, 1, 5, 6, 8, 12, 14, 10, 4, 6, 6, 3, 1, 3, 1, 1, 0,
        0, 0), F6.0y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 2, 3, 2, 2, 3, 4, 7, 1, 13, 12,
        7, 8, 6, 5, 3, 1, 1, 0, 1, 2, 0), F6.5y = c(0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 3,
        0, 2, 1, 1, 5, 4, 8, 13, 7, 13, 7, 5, 4, 2, 1, 0, 2, 1, 0), F7.0y = c(0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
        0, 0, 0, 1, 1, 1, 0, 2, 2, 2, 3, 12, 5, 11, 8, 5, 5, 2, 3, 1,
        0, 1, 0, 1)), .Names = c("MeanD", "F0.5y", "F1.0y", "F1.5y", 
        "F2.0y", "F2.5y", "F3.0y", "F3.5y", "F4.0y", "F4.5y", "F5.0y",
        "F5.5y", "F6.0y", "F6.5y", "F7.0y"), class = "data.frame", row.names = c("1", 
        "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
        "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",
        "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35",
        "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46",
        "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57",
        "58", "59", "60", "61", "62", "63", "64", "65", "66", "67"))
# ages are the time at which these freqs are measured (in years)
ages <- c(0.49, 0.91, 1.41, 1.91, 2.41, 2.92, 3.42, 3.92, 4.41, 4.91,
         5.42, 5.91, 6.42, 6.91)
# To draw the curves, we use more data points on the time axis
agecurves <- c(0.00, 0.49, 0.76, 0.91, 1.17, 1.41, 1.66, 1.91, 2.16, 2.41, 2.67,
         2.92, 3.17, 3.42, 3.67, 3.92, 4.16, 4.41, 4.91, 5.42, 5.91, 6.42, 6.91)
# Since quantile() and nlrq() fonctions have no weight argument,
# we have to expand the data before using these functions
agesize <- expandFreq(ages, freqs)


### Here is the actual calculations and graph plotting ########################
# Quantile regressions for tau=0.05, 0.50 & 0.95 with the choosen growth model
rq.05 <- nlrq(size ~ SSfuzremOrig2(age, Asym, lrc1, lrc2, c0),
    data = agesize, tau = 0.05, trace = FALSE) 
rq.50 <- nlrq(size ~ SSfuzremOrig2(age, Asym, lrc1, lrc2, c0),
    data = agesize, tau = 0.50, trace = FALSE) 
rq.95 <- nlrq(size ~ SSfuzremOrig2(age, Asym, lrc1, lrc2, c0),
    data = agesize, tau = 0.95, trace = FALSE) 
# Predict sizes with these three growth models
c.05 = predict(rq.05, newdata = list(age = agecurves)) 
c.50 = predict(rq.50, newdata = list(age = agecurves)) 
c.95 = predict(rq.95, newdata = list(age = agecurves)) 
curves = data.frame(c.05, c.50, c.95) 
# Plot the graph (histograms with time, curves for the three quantile
# regressions and residuals as boxplots)
rqFreqPlot(ages, freqs, agecurves, curves, barscale = .35, barcol = "gray90",
    boxwex = .15, xlim = c(0, 7), ylim = c(0,67), main = "",
    xlab = expression(paste(italic("t"), " (years)")), 
    ylab1 = expression(paste(italic("D"), " (mm)")), 
    ylab2 = expression(paste(Delta, italic("D"), " (mm)")), las = 1, lty = 1) 
title("Sea urchin growth modeled using quantile regression")