WRS-Ex.R

In WRS: A package of R.R. Wilcox' robust statistics functions

pkgname <- "WRS"
source(file.path(R.home("share"), "R", "examples-header.R"))
options(warn = 1)
library('WRS')

base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
cleanEx()
nameEx("corb")
### * corb

flush(stderr()); flush(stdout())

### Name: corb
### Title: Compute a .95 confidence interval for a correlation
### Aliases: corb

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, y, corfun = pbcor, nboot = 599, SEED = T, ...) 
{
    est <- corfun(x, y, ...)$cor
    if (SEED) 
        set.seed(2)
    data <- matrix(sample(length(y), size = length(y) * nboot, 
        replace = T), nrow = nboot)
    bvec <- apply(data, 1, corbsub, x, y, corfun, ...)
    ihi <- floor(0.975 * nboot + 0.5)
    ilow <- floor(0.025 * nboot + 0.5)
    bsort <- sort(bvec)
    corci <- 1
    corci[1] <- bsort[ilow]
    corci[2] <- bsort[ihi]
    phat <- sum(bvec < 0)/nboot
    sig <- 2 * min(phat, 1 - phat)
    list(cor.ci = corci, p.value = sig, cor.est = est)
  }



cleanEx()
nameEx("outpro")
### * outpro

flush(stderr()); flush(stdout())

### Name: outpro
### Title: Detect outliers using a modification of the Stahel-Donoho
###   projection method
### Aliases: outpro

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (m, gval = NA, center = NA, plotit = T, op = T, MM = F, 
    cop = 3, xlab = "VAR 1", ylab = "VAR 2") 
{
    library(MASS)
    m <- as.matrix(m)
    if (ncol(m) == 1) {
        dis <- (m - median(m))^2/mad(m)^2
        dis <- sqrt(dis)
        crit <- sqrt(qchisq(0.975, 1))
        chk <- ifelse(dis > crit, 1, 0)
        vec <- c(1:nrow(m))
        outid <- vec[chk == 1]
        keep <- vec[chk == 0]
    }
    if (ncol(m) > 1) {
        if (is.na(gval) && cop == 1) 
            gval <- sqrt(qchisq(0.95, ncol(m)))
        if (is.na(gval) && cop != 1) 
            gval <- sqrt(qchisq(0.975, ncol(m)))
        m <- elimna(m)
        if (cop == 1 && is.na(center[1])) {
            if (ncol(m) > 2) 
                center <- dmean(m, tr = 0.5, cop = 1)
            if (ncol(m) == 2) {
                tempd <- NA
                for (i in 1:nrow(m)) tempd[i] <- depth(m[i, 1], 
                  m[i, 2], m)
                mdep <- max(tempd)
                flag <- (tempd == mdep)
                if (sum(flag) == 1) 
                  center <- m[flag, ]
                if (sum(flag) > 1) 
                  center <- apply(m[flag, ], 2, mean)
            }
        }
        if (cop == 2 && is.na(center[1])) {
            center <- cov.mcd(m)$center
        }
        if (cop == 4 && is.na(center[1])) {
            center <- cov.mve(m)$center
        }
        if (cop == 3 && is.na(center[1])) {
            center <- apply(m, 2, median)
        }
        if (cop == 5 && is.na(center[1])) {
            center <- tbs(m)$center
        }
        if (cop == 6 && is.na(center[1])) {
            center <- rmba(m)$center
        }
        if (cop == 7 && is.na(center[1])) {
            center <- spat(m)
        }
        flag <- rep(0, nrow(m))
        outid <- NA
        vec <- c(1:nrow(m))
        for (i in 1:nrow(m)) {
            B <- m[i, ] - center
            dis <- NA
            BB <- B^2
            bot <- sum(BB)
            if (bot != 0) {
                for (j in 1:nrow(m)) {
                  A <- m[j, ] - center
                  temp <- sum(A * B) * B/bot
                  dis[j] <- sqrt(sum(temp^2))
                }
                temp <- idealf(dis)
                if (!MM) 
                  cu <- median(dis) + gval * (temp$qu - temp$ql)
                if (MM) 
                  cu <- median(dis) + gval * mad(dis)
                outid <- NA
                temp2 <- (dis > cu)
                flag[temp2] <- 1
            }
        }
        if (sum(flag) == 0) 
            outid <- NA
        if (sum(flag) > 0) 
            flag <- (flag == 1)
        outid <- vec[flag]
        idv <- c(1:nrow(m))
        keep <- idv[!flag]
        if (ncol(m) == 2) {
            if (plotit) {
                plot(m[, 1], m[, 2], type = "n", xlab = xlab, 
                  ylab = ylab)
                points(m[keep, 1], m[keep, 2], pch = "*")
                if (length(outid) > 0) 
                  points(m[outid, 1], m[outid, 2], pch = "o")
                if (op) {
                  tempd <- NA
                  keep <- keep[!is.na(keep)]
                  mm <- m[keep, ]
                  for (i in 1:nrow(mm)) tempd[i] <- depth(mm[i, 
                    1], mm[i, 2], mm)
                  mdep <- max(tempd)
                  flag <- (tempd == mdep)
                  if (sum(flag) == 1) 
                    center <- mm[flag, ]
                  if (sum(flag) > 1) 
                    center <- apply(mm[flag, ], 2, mean)
                  m <- mm
                }
                points(center[1], center[2], pch = "+")
                x <- m
                temp <- fdepth(m, plotit = F)
                flag <- (temp >= median(temp))
                xx <- x[flag, ]
                xord <- order(xx[, 1])
                xx <- xx[xord, ]
                temp <- chull(xx)
                xord <- order(xx[, 1])
                xx <- xx[xord, ]
                temp <- chull(xx)
                lines(xx[temp, ])
                lines(xx[c(temp[1], temp[length(temp)]), ])
            }
        }
    }
    list(out.id = outid, keep = keep)
  }



cleanEx()
nameEx("pb2gen")
### * pb2gen

flush(stderr()); flush(stdout())

### Name: pb2gen
### Title: Compute a bootstrap confidence interval for the difference
###   between any two parameters corresponding to independent groups
### Aliases: pb2gen

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, y, alpha = 0.05, nboot = 2000, est = onestep, SEED = TRUE, 
    pr = TRUE, ...) 
{
    x <- x[!is.na(x)]
    y <- y[!is.na(y)]
    if (SEED) 
        set.seed(2)
    if (pr) 
        print("Taking bootstrap samples. Please wait.")
    datax <- matrix(sample(x, size = length(x) * nboot, replace = T), 
        nrow = nboot)
    datay <- matrix(sample(y, size = length(y) * nboot, replace = T), 
        nrow = nboot)
    bvecx <- apply(datax, 1, est, ...)
    bvecy <- apply(datay, 1, est, ...)
    bvec <- sort(bvecx - bvecy)
    low <- round((alpha/2) * nboot) + 1
    up <- nboot - low
    temp <- sum(bvec < 0)/nboot + sum(bvec == 0)/(2 * nboot)
    sig.level <- 2 * (min(temp, 1 - temp))
    se <- var(bvec)
    list(ci = c(bvec[low], bvec[up]), p.value = sig.level, sq.se = se)
  }



cleanEx()
nameEx("pball")
### * pball

flush(stderr()); flush(stdout())

### Name: pball
### Title: Percentage Bend Correlation (matrix)
### Aliases: pball

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (m, beta = 0.2) 
{
    if (!is.matrix(m)) 
        stop("Data must be stored in an n by p matrix")
    pbcorm <- matrix(0, ncol(m), ncol(m))
    temp <- matrix(1, ncol(m), ncol(m))
    siglevel <- matrix(NA, ncol(m), ncol(m))
    cmat <- matrix(0, ncol(m), ncol(m))
    for (i in 1:ncol(m)) {
        ip1 <- i
        for (j in ip1:ncol(m)) {
            if (i < j) {
                pbc <- pbcor(m[, i], m[, j], beta)
                pbcorm[i, j] <- pbc$cor
                temp[i, j] <- pbcorm[i, j]
                temp[j, i] <- pbcorm[i, j]
                siglevel[i, j] <- pbc$siglevel
                siglevel[j, i] <- siglevel[i, j]
            }
        }
    }
    tstat <- pbcorm * sqrt((nrow(m) - 2)/(1 - pbcorm^2))
    cmat <- sqrt((nrow(m) - 2.5) * log(1 + tstat^2/(nrow(m) - 
        2)))
    bv <- 48 * (nrow(m) - 2.5)^2
    cmat <- cmat + (cmat^3 + 3 * cmat)/bv - (4 * cmat^7 + 33 * 
        cmat^5 + 240^cmat^3 + 855 * cmat)/(10 * bv^2 + 8 * bv * 
        cmat^4 + 1000 * bv)
    H <- sum(cmat^2)
    df <- ncol(m) * (ncol(m) - 1)/2
    h.siglevel <- 1 - pchisq(H, df)
    list(pbcorm = temp, siglevel = siglevel, H = H, H.siglevel = h.siglevel)
  }



cleanEx()
nameEx("pbcor")
### * pbcor

flush(stderr()); flush(stdout())

### Name: pbcor
### Title: Percentage Bend Correlation (bivariate)
### Aliases: pbcor

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, y, beta = 0.2) 
{
    if (length(x) != length(y)) 
        stop("The vectors do not have equal lengths")
    if (sum(is.na(c(x, y))) != 0) {
        m1 <- matrix(c(x, y), length(x), 2)
        m1 <- elimna(m1)
        x <- m1[, 1]
        y <- m1[, 2]
    }
    temp <- sort(abs(x - median(x)))
    omhatx <- temp[floor((1 - beta) * length(x))]
    temp <- sort(abs(y - median(y)))
    omhaty <- temp[floor((1 - beta) * length(y))]
    a <- (x - pbos(x, beta))/omhatx
    b <- (y - pbos(y, beta))/omhaty
    a <- ifelse(a <= -1, -1, a)
    a <- ifelse(a >= 1, 1, a)
    b <- ifelse(b <= -1, -1, b)
    b <- ifelse(b >= 1, 1, b)
    pbcor <- sum(a * b)/sqrt(sum(a^2) * sum(b^2))
    test <- pbcor * sqrt((length(x) - 2)/(1 - pbcor^2))
    sig <- 2 * (1 - pt(abs(test), length(x) - 2))
    list(cor = pbcor, test = test, siglevel = sig)
  }



cleanEx()
nameEx("trimpb")
### * trimpb

flush(stderr()); flush(stdout())

### Name: trimpb
### Title: Compute a 1-alpha confidence interval for a trimmed mean.
### Aliases: trimpb

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, tr = 0.2, alpha = 0.05, nboot = 2000, WIN = F, win = 0.1, 
    plotit = F, pop = 1, null.value = 0, pr = T, xlab = "X", 
    fr = NA) 
{
    if (pr) {
        print("The p-value returned by the this function is based on the")
        print("null value specified by the argument null.value, which defaults to 0")
    }
    x <- x[!is.na(x)]
    if (WIN) {
        if (win > tr) 
            stop("The amount of Winsorizing must be <= to the amount of trimming")
        x <- winval(x, win)
    }
    crit <- alpha/2
    icl <- round(crit * nboot) + 1
    icu <- nboot - icl
    bvec <- NA
    set.seed(2)
    print("Taking bootstrap samples. Please wait.")
    data <- matrix(sample(x, size = length(x) * nboot, replace = T), 
        nrow = nboot)
    bvec <- apply(data, 1, mean, tr)
    bvec <- sort(bvec)
    p.value <- mean(bvec < null.value) + 0.5 * mean(bvec == null.value)
    p.value <- 2 * min(p.value, 1 - p.value)
    ci <- NA
    ci[1] <- bvec[icl]
    ci[2] <- bvec[icu]
    if (plotit) {
        if (pop == 1) 
            rdplot(as.vector(bvec), fr = fr, xlab = xlab)
        if (pop == 2) 
            kdplot(as.vector(bvec), rval = rval)
        if (pop == 3) 
            boxplot(as.vector(bvec))
        if (pop == 4) 
            stem(as.vector(bvec))
        if (pop == 5) 
            hist(as.vector(bvec))
        if (pop == 6) 
            akerd(as.vector(bvec), xlab = xlab)
    }
    list(ci = ci, p.value = p.value)
  }



cleanEx()
nameEx("twocor")
### * twocor

flush(stderr()); flush(stdout())

### Name: twocor
### Title: Compute a .95 confidence interval for the difference between two
###   correlation coefficients corresponding to two independent groups
### Aliases: twocor

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x1, y1, x2, y2, corfun = pbcor, nboot = 599, alpha = 0.05, 
    SEED = T, ...) 
{
    if (SEED) 
        set.seed(2)
    print("Taking bootstrap samples. Please wait.")
    data1 <- matrix(sample(length(y1), size = length(y1) * nboot, 
        replace = T), nrow = nboot)
    bvec1 <- apply(data1, 1, corbsub, x1, y1, corfun, ...)
    data2 <- matrix(sample(length(y2), size = length(y2) * nboot, 
        replace = T), nrow = nboot)
    bvec2 <- apply(data2, 1, corbsub, x2, y2, corfun, ...)
    bvec <- bvec1 - bvec2
    bsort <- sort(bvec)
    term <- alpha/2
    ihi <- floor((1 - term) * nboot + 0.5)
    ilow <- floor(term * nboot + 0.5)
    corci <- 1
    corci[1] <- bsort[ilow]
    corci[2] <- bsort[ihi]
    r1 <- corfun(x1, y1)$cor
    r2 <- corfun(x2, y2)$cor
    reject <- "NO"
    if (corci[1] > 0 || corci[2] < 0) 
        reject = "YES"
    list(r1 = r1, r2 = r2, ci.dif = corci, reject = reject)
  }



cleanEx()
nameEx("twopcor")
### * twopcor

flush(stderr()); flush(stdout())

### Name: twopcor
### Title: Compute a .95 confidence interval for the difference between two
###   Pearson correlations corresponding to two independent groups
### Aliases: twopcor

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x1, y1, x2, y2, SEED = T) 
{
    nboot <- 599
    if (SEED) 
        set.seed(2)
    X <- elimna(cbind(x1, y1))
    x1 <- X[, 1]
    y1 <- X[, 2]
    X <- elimna(cbind(x2, y2))
    x2 <- X[, 1]
    y2 <- X[, 2]
    print("Taking bootstrap samples; please wait")
    data1 <- matrix(sample(length(y1), size = length(y1) * nboot, 
        replace = T), nrow = nboot)
    bvec1 <- apply(data1, 1, pcorbsub, x1, y1)
    data2 <- matrix(sample(length(y2), size = length(y2) * nboot, 
        replace = T), nrow = nboot)
    bvec2 <- apply(data2, 1, pcorbsub, x2, y2)
    bvec <- bvec1 - bvec2
    ilow <- 15
    ihi <- 584
    if (length(y1) + length(y2) < 250) {
        ilow <- 14
        ihi <- 585
    }
    if (length(y1) + length(y2) < 180) {
        ilow <- 11
        ihi <- 588
    }
    if (length(y1) + length(y2) < 80) {
        ilow <- 8
        ihi <- 592
    }
    if (length(y1) + length(y2) < 40) {
        ilow <- 7
        ihi <- 593
    }
    bsort <- sort(bvec)
    r1 <- cor(x1, y1)
    r2 <- cor(x2, y2)
    ci <- c(bsort[ilow], bsort[ihi])
    list(r1 = r1, r2 = r2, ci = ci)
  }



cleanEx()
nameEx("yuen")
### * yuen

flush(stderr()); flush(stdout())

### Name: yuen
### Title: Yuen's test for trimmed means
### Aliases: yuen

### ** Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, y, tr = 0.2, alpha = 0.05) 
{
    if (tr == 0.5) 
        stop("Using tr=.5 is not allowed; use a method designed for medians")
    if (tr > 0.25) 
        print("Warning: with tr>.25 type I error control might be poor")
    x <- x[!is.na(x)]
    y <- y[!is.na(y)]
    h1 <- length(x) - 2 * floor(tr * length(x))
    h2 <- length(y) - 2 * floor(tr * length(y))
    q1 <- (length(x) - 1) * winvar(x, tr)/(h1 * (h1 - 1))
    q2 <- (length(y) - 1) * winvar(y, tr)/(h2 * (h2 - 1))
    df <- (q1 + q2)^2/((q1^2/(h1 - 1)) + (q2^2/(h2 - 1)))
    crit <- qt(1 - alpha/2, df)
    dif <- mean(x, tr) - mean(y, tr)
    low <- dif - crit * sqrt(q1 + q2)
    up <- dif + crit * sqrt(q1 + q2)
    test <- abs(dif/sqrt(q1 + q2))
    yuen <- 2 * (1 - pt(test, df))
    list(ci = c(low, up), p.value = yuen, dif = dif, se = sqrt(q1 + 
        q2), teststat = test, crit = crit, df = df)
  }



### * <FOOTER>
###
options(digits = 7L)
base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
grDevices::dev.off()
###
### Local variables: ***
### mode: outline-minor ***
### outline-regexp: "\\(> \\)?### [*]+" ***
### End: ***
quit('no')