.packageName <- "BHH2"
"anovaPlot" <-
function (obj, stacked = TRUE, base = TRUE, axes = TRUE, faclab = TRUE, 
    labels = FALSE, cex = par("cex"), cex.lab = par("cex.lab"), 
    ...) 
{
    if (!any(class(obj) == "lm")) 
        stop(paste("Object", deparse(substitute(obj)), "should be of class 'lm'"))
    if (!any(class(obj) == "aov")) 
        obj <- do.call("aov", as.list(obj$call[-1]))
    tables <- model.tables(obj, type = "effects")[[1]]
    Residuals <- resid(obj)
    factors <- names(tables)
    k <- length(factors)
    Df <- anova(obj)[, "Df"]
    Scale.Factor <- sqrt(Df[length(Df)]/Df[-length(Df)])
    lst <- lst.dev <- lst.names <- list()
    for (i in 1:k) {
        nf <- length(unlist(strsplit(names(tables[i]), ":")))
        if (nf == 1) {
            label <- dimnames(tables[i][[1]])[[1]]
            eff <- as.numeric(tables[i][[1]])
            names(eff) <- label
        }
        else {
            label <- dimnames(tables[i][[1]])[[1]]
            eff <- as.numeric(tables[i][[1]])
            for (j in 2:nf) {
                lab <- dimnames(tables[i][[1]])[[j]]
                label <- paste(rep(label, length(lab)), rep(lab, 
                  each = length(label)), sep = ":")
            }
            names(eff) <- label
        }
        lst.dev[[factors[i]]] <- Scale.Factor[i] * eff
    }
    xmax <- max(abs(range(c(unlist(lst.dev), Residuals))))
    xlim <- c(-xmax, +xmax)
    xpd <- par("xpd")
    on.exit(par(xpd = xpd))
    par(xpd = TRUE)
    plot(c(0, 1), c(0, 1), xlim = xlim, type = "n", xlab = "", 
        ylab = "", frame = FALSE, axes = FALSE, ...)
    if (stacked) 
        h <- 1/(k + 2)
    else h <- 1/(k + 1)
    hinc <- h/(k + 1)
    for (i in 1:k) {
        y <- (k - i + 2) * h
        if (labels) 
            lab <- names(lst.dev[[i]])
        else lab <- NULL
        lst[[factors[i]]] <- dots(lst.dev[[i]], y = y + hinc, 
            xlim = xlim, hmax = y + hinc + h, stacked = stacked, 
            base = base, axes = FALSE, labels = lab, cex = cex)
        if (faclab) 
            text(0, y, labels = paste("scaled <", factors[i], 
                "> deviations"), adj = 0.5, cex = cex.lab)
    }
    lst[["Rediduals"]] <- dots(Residuals, y = 0, xlim = xlim, 
        hmax = 2 * h, stacked = stacked, base = FALSE, axes = FALSE, 
        labels = NULL)
    if (axes) 
        axis(1, at = pretty(xlim), mgp = c(1.5, 0.5, 0), line = -1 + 
            2 * hinc)
    if (faclab & axes) 
        mtext("Residuals", side = 1, cex = cex.lab, line = 1 + 
            hinc)
    if (faclab & !axes) 
        mtext("Residuals", side = 1, cex = cex.lab, line = 0)
    invisible(lst)
}
"dotPlot" <-
function (x, y = 0, xlim = range(x, na.rm = TRUE), xlab = NULL, 
    scatter = FALSE, hmax = 1, base = TRUE, axes = TRUE, frame = FALSE, 
    pch = 21, pch.size = "x", labels = NULL, hcex = 1, cex = par("cex"), 
    cex.axis = par("cex.axis"), ...) 
{
    if (is.null(xlab)) 
        xlab <- deparse(substitute(x))
    x <- x[!is.na(x)]
    xpd <- par("xpd")
    par(xpd = TRUE)
    on.exit(par(xpd = xpd))
    plot(c(0, 1), c(0, 1), xlim = xlim, type = "n", axes = FALSE, 
        cex = cex, cex.axis = cex.axis, frame = frame, xlab = xlab, 
        ylab = "", ...)
    if (axes) 
        axis(1, cex.axis = cex.axis)
    if (scatter) {
        dots(x, y = y, xlim = xlim, stacked = FALSE, hmax = hmax, 
            base = base, axes = FALSE, pch = pch, pch.size = pch.size, 
            labels = labels, hcex = hcex, cex = cex, cex.axis = cex.axis)
        y <- y + 2 * strheight(pch.size, units = "user")
        xlab <- ""
        axes <- FALSE
        base = FALSE
    }
    coord <- dots(x, y, xlim = xlim, stacked = TRUE, hmax = hmax, 
        base = base, axes = FALSE, pch = pch, pch.size = pch.size, 
        labels = labels, hcex = hcex, cex = cex, cex.axis = cex.axis)
    invisible(coord)
}
"dots" <-
function (x, y = 0.1, xlim = range(x, na.rm = TRUE), stacked = FALSE, 
    hmax = 0.5, base = TRUE, axes = FALSE, pch = 21, pch.size = "x", 
    labels = NULL, hcex = 1, cex = par("cex"), cex.axis = par("cex.axis")) 
{
    x <- x[!is.na(x)]
    hdots <- y
    xmin <- xlim[1]
    xmax <- xlim[2]
    x <- x[(x >= xmin) & (x <= xmax)]
    b <- strwidth(pch.size, units = "user", cex = cex)
    h <- strheight(pch.size, units = "user", cex = hcex * cex)
    if (stacked) {
        if (xmax - xmin < b) {
            error("x-dimension resolution problem")
        }
        else {
            xu <- seq(xmin, xmax, by = b)
        }
        m <- length(xu)
        tab <- data.frame(j = 1:m, k = rep(0, m), xu = xu)
        n <- length(x)
        y <- rep(0, n)
        for (i in 1:n) {
            l <- max(tab$j[tab$xu <= x[i]])
            x[i] <- xu[l] + b/2
            tab$k[l] <- 1 + tab$k[l]
            y[i] <- tab$k[l]
        }
        y <- y * h
        u <- hdots + max(y)
        if (hmax <= hdots) 
            warning(paste("dot base <hdots=", hdots, "> higher than maximum column height <hmax=", 
                hmax, ">...", sep = ""))
        if (u > hmax) 
            y <- (hmax - hdots) * y/u
        y <- hdots + y
    }
    else {
        y <- rep(hdots, length(x))
    }
    if (!is.null(labels)) 
        text(x, y, labels = labels, cex = cex)
    else points(x, y, pch = pch, cex = cex)
    points.coord <- data.frame(x, y)
    if (axes) {
        segments(xmin - b, hdots - h/4, xmax + b, hdots - h/4)
        x <- pretty(x, n = 3, h = 0.5)
        x <- x[(x > xmin - b) & (x < xmax + b)]
        for (i in seq(x)) segments(x[i], hdots - h/4, x[i], hdots - 
            h/2)
        y <- rep(hdots - h, length(x))
        text(x, y, labels = x, cex = cex.axis)
    }
    if (base && !axes) 
        segments(xmin - b, hdots - h/4, xmax + b, hdots - h/4)
    invisible(points.coord)
}
"ffDesMatrix" <-
function (k, gen = NULL) 
{
    N <- 2^k
    X <- matrix(NA, nrow = N, ncol = k)
    for (j in 1:k) X[, j] <- rep(sort(rep(c(-1, 1), N/2^j)), 
        2^(j - 1))
    X <- X[, ncol(X):1]
    if (is.null(gen)) 
        return(X)
    for (i in 1:length(gen)) {
        ind <- trunc(gen[[i]])
        if (any(abs(ind) > k)) 
            stop(paste("generator:", paste(ind[1], "=", paste(ind[-1], 
                collapse = "*")), "includes undefined columns"))
        x <- rep(sign(ind[1]), N)
        for (j in ind[-1]) x <- x * X[, j]
        X[, abs(ind[1])] <- x
    }
    X <- unique(X)
    X
}
"ffFullMatrix" <-
function (X, x, maxInt, blk = NULL) 
{
    if (!is.data.frame(X)) 
        X <- as.data.frame(X)
    Z <- data.frame(one = rep(1, nrow(X)))
    k <- 0
    if (!is.null(blk)) {
        blk <- as.matrix(blk)
        if (nrow(blk) != nrow(X)) 
            stop("Matrix and block should have the same number of rows.")
        k <- ncol(blk)
        Z <- cbind(Z, blk)
        names(Z)[-1] <- paste("bk", seq(k), sep = "")
    }
    ord <- min(length(x), maxInt)
    nT <- rep(0, ord)
    for (i in seq(ord)) {
        tt <- subsets(length(x), i, x)
        if (is.null(dim(tt))) 
            tt <- matrix(tt, nrow = 1)
        for (j in 1:nrow(tt)) {
            nT[i] <- nT[i] + 1
            Z <- cbind(Z, eval(parse(text = (paste("X[", tt[j, 
                ], "]", collapse = "*", sep = "")))))
            names(Z)[ncol(Z)] <- paste("x", tt[j, ], collapse = "*", 
                sep = "")
        }
    }
    nT <- c(k, nT)
    if (ord > 1) 
        names(nT) <- c("blk", "main", paste("int", 2:ord, sep = "."))
    else names(nT) <- c("blk", "main")
    list(Xa = as.matrix(Z), x = x, maxInt = ord, nTerms = nT)
}
"lambdaPlot" <-
function (mod, lambda = seq(-1, 1, by = 0.1), stat = "F", global = TRUE, 
    cex = par("cex"), ...) 
{
    if (stat == "F") {
        org.fit <- mod
        y <- org.fit$model[, 1]
        resp <- names(org.fit$model)[1]
        print(resp)
        dat <- org.fit$model
        form <- as.formula(org.fit$call$formula)
        tt <- lm(form, data = dat)
        sav <- anova(tt)
        n <- nrow(sav)
        numdf <- sum(sav[-n, "Df"])
        dendf <- sav[n, "Df"]
        k <- length(lambda)
        if (global) {
            f.lambda <- matrix(NA, nrow = 1, ncol = k)
            dimnames(f.lambda) <- list("Model", paste("l", round(lambda, 
                2), sep = ""))
            for (j in seq(lambda)) {
                l <- lambda[j]
                if (l == 0) 
                  dat[, resp] <- log(y)
                else dat[, resp] <- (y^l - 1)/l
                tt <- lm(form, data = dat)
                f.lambda[1, j] <- (sum(anova(tt)[-n, "Sum Sq"])/numdf)/(anova(tt)[n, 
                  "Sum Sq"]/dendf)
            }
        }
        else {
            f.lambda <- matrix(NA, nrow = n, k)
            dimnames(f.lambda) <- list(dimnames(sav)[[1]], paste("l", 
                round(lambda, 2), sep = ""))
            for (j in seq(lambda)) {
                l <- lambda[j]
                if (l == 0) 
                  dat[, resp] <- log(y)
                else dat[, resp] <- (y^l - 1)/l
                tt <- lm(form, data = dat)
                f.lambda[, j] <- anova(tt)[, "F value"]
            }
            f.lambda <- f.lambda[-n, ]
        }
        Labels <- data.frame(term = dimnames(f.lambda)[[1]], 
            label = LETTERS[seq(nrow(f.lambda))])
        plot(lambda, lambda, xlim = range(lambda), ylim = range(f.lambda), 
            type = "n", xlab = "lambda", ylab = "F", ...)
        for (i in 1:nrow(f.lambda)) lines(lambda, f.lambda[i, 
            ])
        xlab <- lambda[k]
        ylab <- f.lambda[, k]
        lab <- paste(" ", as.character(Labels[, "label"]), sep = "")
        for (i in 1:nrow(f.lambda)) text(lambda[k], f.lambda[, 
            k], labels = lab, adj = 0)
        lab <- paste(as.character(Labels[, "label"]), " ", sep = "")
        for (i in 1:nrow(f.lambda)) text(lambda[1], f.lambda[, 
            1], labels = lab, adj = 1)
        print(Labels)
        invisible(list(lambda = lambda, f.lambda = f.lambda))
    }
    else if (stat == "t") {
        y <- mod$model[, 1]
        org.fit <- lm(y ~ ., qr = TRUE, data = mod$model[, -1])
        QR <- org.fit$qr
        n <- length(y)
        p <- length(coef(org.fit))
        idx <- 1:p
        rdf <- n - p
        coef.lambda <- matrix(NA, nrow = p, ncol = length(lambda))
        dimnames(coef.lambda) <- list(names(coef(org.fit)), paste("l", 
            round(lambda, 2), sep = ""))
        t.lambda <- se.lambda <- coef.lambda
        for (j in seq(lambda)) {
            l <- lambda[j]
            if (l == 0) 
                y.lambda <- log(y)
            else y.lambda <- (y^l - 1)/l
            resvar <- sum(qr.resid(QR, y.lambda)^2)/rdf
            coef.lambda[, j] <- qr.coef(QR, y.lambda)
            R <- chol2inv(QR$qr[idx, idx, drop = FALSE])
            se.lambda[, j] <- sqrt(diag(R) * resvar)
        }
        t.lambda <- coef.lambda/se.lambda
        Labels <- data.frame(term = names(coef(org.fit)), label = c(" ", 
            LETTERS[seq(nrow(t.lambda) - 1)]))
        plot(lambda, lambda, xlim = range(lambda), ylim = range(t.lambda[-1, 
            ]), type = "n", xlab = "lambda", ylab = "t", ...)
        for (i in 2:nrow(t.lambda)) lines(lambda, t.lambda[i, 
            ])
        xlab <- lambda[length(lambda)]
        ylab <- t.lambda[, ncol(t.lambda)]
        lab <- paste(" ", as.character(Labels[, "label"]), sep = "")
        for (i in 2:nrow(t.lambda)) text(lambda[length(lambda)], 
            t.lambda[, ncol(t.lambda)], labels = lab, adj = 0)
        lab <- paste(as.character(Labels[, "label"]), " ", sep = "")
        for (i in 2:nrow(t.lambda)) text(lambda[1], t.lambda[, 
            1], labels = lab, adj = 1)
        print(Labels)
        invisible(list(lambda = lambda, coef = coef.lambda, se = se.lambda))
    }
    else {
        warning("argument stat should be either \"F\" or \"t\"")
        invisible(NULL)
    }
}
"permtest" <-
function (x, y = NULL) 
{
    if (is.null(y)) {
        mx <- mean(x)
        n <- length(x)
        t.obs <- mx/sqrt(var(x)/n)
        N <- 2^n
        mat <- matrix(NA, nrow = N, ncol = n)
        for (j in 1:n) {
            m <- 2^j
            mat[, j] <- rep(c(rep(-1, N/m), rep(+1, N/m)), m/2)
        }
        d <- as.numeric(mat %*% x/n)
        k <- length(which(d > mx))
        return(c(N = N, t.obs = t.obs, "t-Dist-P(>t)" = 1 - pt(t.obs, 
            n - 1), "PermDist-P(>t)" = k/N))
    }
    else {
        x <- x[!is.na(x)]
        y <- y[!is.na(y)]
        nx <- length(x)
        ny <- length(y)
        S2x <- sum(x^2) - sum(x)^2/nx
        S2y <- sum(y^2) - sum(y)^2/ny
        t.stat <- function(x, y) {
            (mean(x) - mean(y))/sqrt((S2x + S2y)/(nx + ny - 2) * 
                (1/nx + 1/ny))
        }
        f.stat <- function(x, y) {
            (S2x/(nx - 1))/(S2y/(ny - 1))
        }
        t.obs <- t.stat(x, y)
        f.obs <- f.stat(x, y)
        z <- c(x, y)
        n <- nx + ny
        mat <- subsets(n, nx)
        N <- nrow(mat)
        kt <- kf <- 0
        for (i in 1:nrow(mat)) {
            x <- z[mat[i, ]]
            y <- z[-mat[i, ]]
            S2x <- sum(x^2) - sum(x)^2/nx
            S2y <- sum(y^2) - sum(y)^2/ny
            if (t.obs < t.stat(x, y)) 
                kt <- kt + 1
            if (f.obs < f.stat(x, y)) 
                kf <- kf + 1
        }
        return(c(N = N, t.obs = t.obs, "t-Dist:P(>t)" = 1 - pt(t.obs, 
            nx + ny - 2), "PermDist:P(>t)" = kt/N, F.obs = f.obs, 
            "F-Dist:P(>F)" = 1 - pf(f.obs, nx - 1, ny - 1), "PermDist:P(>F)" = kf/N))
    }
}
"subsets" <-
function (n, r, v = 1:n) 
if (r <= 0) vector(mode(v), 0) else if (r >= n) v[1:n] else {
    rbind(cbind(v[1], Recall(n - 1, r - 1, v[-1])), Recall(n - 
        1, r, v[-1]))
}
