# NOT RUN {
# Example from Box(1991)
data(BoxM)
Gaptest(BoxM)
## The function is currently defined as
function (DesY)
{
ncheck <- dim(DesY)
ncheck <- ncheck[1]
tcnd = TRUE
if (ncheck == 8) {
tcnd = FALSE
}
if (ncheck == 16) {
tcnd = FALSE
}
if (ncheck == 32) {
tcnd = FALSE
}
if (tcnd) {
stop("This function only works for 8, 16, or 32 run designs",
"\n")
}
else {
if (ncheck == 8)
ncheck = 16
critg16 <- c(1.7884, 5.1009)
critg32 <- c(1.7297, 5.8758)
modf <- lm(y ~ (.)^4, x = TRUE, data = DesY)
nbeta <- dim(DesY)
nbeta <- nbeta[1]
he <- modf$coef
selcol <- which(!is.na(he))
he <- he[selcol]
he <- he[-1]
p <- length(he)
n <- p + 1
cn1 <- names(he)
ccn1 <- gsub("[^A-Z]", "", cn1)
names(he) <- ccn1
ahe <- abs(he)
s0 <- 1.5 * median(ahe)
selhe <- ahe < (2.5 * s0)
pse = 1.5 * median(ahe[selhe])
gap <- gapstat(he, pse)
if (ncheck == 16) {
test = (gap > critg16[1])
}
else {
test = (gap > critg32[1])
}
if (test) {
X <- modf$x
X <- X[, selcol]
X <- X[, -1]
se <- as.matrix(sign(he), nrow = 1)
sigef <- LGB(he, rpt = FALSE, plt = FALSE)
for (i in 1:length(he)) {
if (sigef[i] == "yes") {
se[i] = 0
}
}
sp <- X %*% se
asp <- abs(sp)
oo <- max.col(t(asp))
ae <- abs(he)
sae <- sort(ae)
nsmall <- round(length(he)/2)
bias <- 2 * sum(sae[1:nsmall])
y <- DesY$y
ycorr <- DesY$y
ycorr[oo] <- ycorr[oo] + (-1 * sign(sp[oo])) * bias
detect <- c(rep("no", n))
detect[oo] <- "yes"
cat("Initial Outlier Report", "\n")
cat("Standardized-Gap = ", gap, "Significant at 50th percentile",
"\n")
DesYc <- cbind(DesY[, 1:(dim(DesY)[2] - 1)], ycorr)
modf <- lm(ycorr ~ (.)^4, x = TRUE, data = DesYc)
che <- modf$coef
che <- che[!is.na(che)]
che <- che[-1]
p <- length(che)
n <- p + 1
cn <- names(che)
ccn <- gsub("[^A-Z]", "", cn)
names(che) <- ccn
ache <- abs(che)
s0 <- 1.5 * median(ache)
selche <- ache < (2.5 * s0)
psec = 1.5 * median(ache[selche])
gap <- gapstat(he, psec)
if (ncheck == 16)
test2 = (gap > critg16[2])
else test2 = (gap > critg32[2])
if (test2) {
cat("Final Outlier Report", "\n")
cat("Standardized-Gap = ", gap, "Significant at 99th percentile",
"\n")
cat(" ", "\n")
cat(" Corrrected Data Report ", "\n")
cat("Response Corrected Response Detect Outlier",
"\n")
cat(paste(format(DesY$y, width = 8), format(DesYc$ycorr,
width = 13), " ", format(detect,
width = 10), "\n"), sep = "")
tce <- LGB(che)
}
else {
cat("Final Outlier Report", "\n")
cat("No significant outlier detected in second pass",
"\n")
LGB(he)
cat(" ", "\n")
}
}
}
}
# }
Run the code above in your browser using DataLab