daewr (version 1.0-1)

LGB: This function uses the LGB Method to detect significant effects in unreplicated fractional factorials.

Description

This function uses the LGB Method to detect significant effects in unreplicated fractional factorials.

Usage

LGB(Beta, alpha = 0.05, rpt = TRUE, plt = TRUE, pltl = TRUE)

Arguments

Beta

input - this is the numeric vector of effects or coefficients to be tested

alpha

input - This is the significance level of the test

rpt

input - this is a logical variable that controls whether the report is written (default is TRUE)

plt

input - this is a logical variable that controls whether a half-normal plot is made (default is TRUE)

pltl

input - this is a logical variable that controls whether the significance limit line is drawn on the half-normal plot (default is TRUE)

References

Lawson, J., Grimshaw, S., Burt, J. (1998) "A quantitative method for identifying active contrasts in unreplicated factorial experiments based on the half-normal plot", Computational Statistics and Data Analysis, 26, 425-436.

Examples

data(chem)
modf<-lm(y~A*B*C*D,data=chem)
LGB(coef(modf)[-1],rpt=FALSE)


## The function is currently defined as
function (Beta, alpha = 0.05, rpt = TRUE, plt = TRUE, pltl = TRUE) 
{
    siglev <- c(0.1, 0.05, 0.025, 0.01)
    df <- c(7, 8, 11, 15, 16, 17, 26, 31, 32, 35, 63, 127)
    crittab <- matrix(c(1.265, 1.196, 1.161, 1.122, 1.11, 1.106, 
        1.072, 1.063, 1.06, 1.059, 1.037, 1.023, 1.534, 1.385, 
        1.291, 1.201, 1.186, 1.178, 1.115, 1.099, 1.093, 1.091, 
        1.056, 1.034, 1.889, 1.606, 1.449, 1.297, 1.274, 1.26, 
        1.165, 1.14, 1.13, 1.127, 1.074, 1.043, 2.506, 2.026, 
        1.74, 1.447, 1.421, 1.377, 1.232, 1.197, 1.185, 1.178, 
        1.096, 1.058), ncol = 4, byrow = FALSE)
    colind <- which(siglev == alpha, arr.ind = TRUE)
    if (length(colind) == 0) {
        stop("this function works only when alpha= .1, .05, .025 or .01")
    }
    rowind <- which(df == length(Beta), arr.ind = TRUE)
    if (length(rowind) == 0) {
        stop("this function works only for coefficent vectors of 
		length 7,8,11,15,16,26,31,32,35,63,or 127")
    }
    critL <- crittab[rowind, colind]
    acj <- abs(Beta)
    ranks <- rank(acj, ties.method = "first")
    s0 <- 1.5 * median(acj)
    p <- (ranks - 0.5)/length(Beta)
    z <- qnorm((p + 1)/2)
    moda <- lm(acj ~ -1 + z)
    beta1 <- moda$coef
    sel <- acj < 2.5 * s0
    modi <- lm(acj[sel] ~ -1 + z[sel])
    beta2 <- modi$coef
    Rn <- beta1/beta2
    pred <- beta2 * z
    n <- length(acj[sel])
    df <- n - 1
    sig <- sqrt(sum(modi$residuals^2)/df)
    se.pred <- sig * (1 + 1/n + (z^2)/sum(z[sel]^2))^0.5
    pred.lim <- pred + qt(0.975, df) * se.pred
    sigi <- c(rep("no", length(Beta)))
    sel2 <- acj > pred.lim
    sigi[sel2] <- "yes"
    if (plt) {
        plot(z, acj, xlab = "Half Normal Scores", ylab = "Absoulute Effects")
        lines(sort(z), sort(pred), lty = 1)
        for (i in 1:length(Beta)) {
            if (sigi[i] == "yes") 
                text(z[i], acj[i], names(Beta)[i], pos = 1)
        }
        if (pltl) {
            lines(sort(z), sort(pred.lim), lty = 3)
        }
    }
    if (rpt) {
        cat("Effect Report", "")
        cat("", "")
        cat("Label     Half Effect    Sig(.05)", "")
        cat(paste(format(names(Beta), width = 8), format(Beta, 
            width = 8), "", format(sigi, width = 10), ""), 
            sep = "")
        cat("", "")
        cat("Lawson, Grimshaw & Burt Rn Statistic = ", Rn, "")
        cat("95th percentile of Rn = ", critL, "")
    }
    return(sigi)
  }