ancbbmed:

Usage Arguments Examples

Usage

1
ancbbmed(x1, y1, x2, y2, fr1 = 1, fr2 = 1, nboot = 100, pts = NA, plotit = TRUE, SEED = TRUE, alpha = 0.05)

Arguments

x1
y1
x2
y2
fr1
fr2
nboot
pts
plotit
SEED
alpha

Examples

  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
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
##---- 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, fr1 = 1, fr2 = 1, nboot = 100, pts = NA, 
    plotit = TRUE, SEED = TRUE, alpha = 0.05) 
{
    if (SEED) 
        set.seed(2)
    if (is.na(pts[1])) {
        isub <- c(1:5)
        test <- c(1:5)
        xorder <- order(x1)
        y1 <- y1[xorder]
        x1 <- x1[xorder]
        xorder <- order(x2)
        y2 <- y2[xorder]
        x2 <- x2[xorder]
        n1 <- 1
        n2 <- 1
        vecn <- 1
        for (i in 1:length(x1)) n1[i] <- length(y1[near(x1, x1[i], 
            fr1)])
        for (i in 1:length(x1)) n2[i] <- length(y2[near(x2, x1[i], 
            fr2)])
        for (i in 1:length(x1)) vecn[i] <- min(n1[i], n2[i])
        sub <- c(1:length(x1))
        isub[1] <- min(sub[vecn >= 12])
        isub[5] <- max(sub[vecn >= 12])
        isub[3] <- floor((isub[1] + isub[5])/2)
        isub[2] <- floor((isub[1] + isub[3])/2)
        isub[4] <- floor((isub[3] + isub[5])/2)
        mat <- matrix(NA, 5, 7)
        dimnames(mat) <- list(NULL, c("X", "n1", "n2", "DIF", 
            "ci.low", "ci.hi", "p.value"))
        gv1 <- vector("list")
        for (i in 1:5) {
            j <- i + 5
            temp1 <- y1[near(x1, x1[isub[i]], fr1)]
            temp2 <- y2[near(x2, x1[isub[i]], fr2)]
            temp1 <- temp1[!is.na(temp1)]
            temp2 <- temp2[!is.na(temp2)]
            mat[i, 1] <- x1[isub[i]]
            mat[i, 2] <- length(temp1)
            mat[i, 3] <- length(temp2)
            mat[, 4] <- runmbo(x1, y1, pts = x1[isub], pyhat = TRUE, 
                plotit = FALSE, SEED = FALSE, est = tmean) - 
                runmbo(x2, y2, pts = x1[isub], pyhat = TRUE, 
                  plotit = FALSE, SEED = FALSE, est = median)
            gv1[[i]] <- temp1
            gv1[[j]] <- temp2
        }
        I1 <- diag(5)
        I2 <- 0 - I1
        con <- rbind(I1, I2)
        estmat1 <- matrix(nrow = nboot, ncol = length(isub))
        estmat2 <- matrix(nrow = nboot, ncol = length(isub))
        data1 <- matrix(sample(length(y1), size = length(y1) * 
            nboot, replace = TRUE), nrow = nboot)
        data2 <- matrix(sample(length(y2), size = length(y2) * 
            nboot, replace = TRUE), nrow = nboot)
        for (ib in 1:nboot) {
            estmat1[ib, ] = runmbo(x1[data1[ib, ]], y1[data1[ib, 
                ]], pts = x1[isub], pyhat = TRUE, plotit = FALSE, 
                SEED = FALSE, est = median)
            estmat2[ib, ] = runmbo(x2[data2[ib, ]], y2[data2[ib, 
                ]], pts = x1[isub], pyhat = TRUE, plotit = FALSE, 
                SEED = FALSE, est = median)
        }
        dif <- (estmat1 < estmat2)
        dif0 <- (estmat1 == estmat2)
        pvals = apply(dif, 2, mean, na.rm = TRUE) + 0.5 * apply(dif0, 
            2, mean, na.rm = TRUE)
        tmat <- rbind(pvals, 1 - pvals)
        pvals = 2 * apply(tmat, 2, min)
        mat[, 7] <- pvals
        for (ij in 1:length(isub)) {
            dif <- estmat1[, ij] - estmat2[, ij]
            dif <- elimna(dif)
            nbad <- length(dif)
            lo <- round(nbad * alpha/2)
            hi <- nbad - lo
            dif <- sort(dif)
            mat[ij, 5] <- dif[lo]
            mat[ij, 6] <- dif[hi]
        }
    }
    if (!is.na(pts[1])) {
        n1 <- 1
        n2 <- 1
        vecn <- 1
        for (i in 1:length(pts)) {
            n1[i] <- length(y1[near(x1, pts[i], fr1)])
            n2[i] <- length(y2[near(x2, pts[i], fr2)])
            if (n1[i] <= 5) 
                print(paste("Warning, there are", n1[i], " points corresponding to the design point X=", 
                  pts[i]))
            if (n2[i] <= 5) 
                print(paste("Warning, there are", n2[i], " points corresponding to the design point X=", 
                  pts[i]))
        }
        mat <- matrix(NA, length(pts), 7)
        dimnames(mat) <- list(NULL, c("X", "n1", "n2", "DIF", 
            "ci.low", "ci.hi", "p.value"))
        gv <- vector("list", 2 * length(pts))
        for (i in 1:length(pts)) {
            g1 <- y1[near(x1, pts[i], fr1)]
            g2 <- y2[near(x2, pts[i], fr2)]
            g1 <- g1[!is.na(g1)]
            g2 <- g2[!is.na(g2)]
            j <- i + length(pts)
            gv[[i]] <- g1
            gv[[j]] <- g2
        }
        I1 <- diag(length(pts))
        I2 <- 0 - I1
        con <- rbind(I1, I2)
        isub = c(1:length(pts))
        estmat1 <- matrix(nrow = nboot, ncol = length(isub))
        estmat2 <- matrix(nrow = nboot, ncol = length(isub))
        data1 <- matrix(sample(length(y1), size = length(y1) * 
            nboot, replace = TRUE), nrow = nboot)
        data2 <- matrix(sample(length(y2), size = length(y2) * 
            nboot, replace = TRUE), nrow = nboot)
        est1 = runmbo(x1, y1, pts = pts, pyhat = TRUE, plotit = FALSE, 
            SEED = FALSE, est = median)
        est2 = runmbo(x2, y2, pts = pts, pyhat = TRUE, plotit = FALSE, 
            SEED = FALSE, est = median)
        mat[, 4] <- est1 - est2
        for (ib in 1:nboot) {
            estmat1[ib, ] = runmbo(x1[data1[ib, ]], y1[data1[ib, 
                ]], pts = pts, pyhat = TRUE, plotit = FALSE, 
                SEED = FALSE, est = median)
            estmat2[ib, ] = runmbo(x2[data2[ib, ]], y2[data2[ib, 
                ]], pts = pts, pyhat = TRUE, plotit = FALSE, 
                SEED = FALSE, est = median)
        }
        dif <- (estmat1 < estmat2)
        dif0 <- (estmat1 == estmat2)
        pvals = apply(dif, 2, mean, na.rm = TRUE) + 0.5 * apply(dif0, 
            2, mean, na.rm = TRUE)
        tmat <- rbind(pvals, 1 - pvals)
        pvals = 2 * apply(tmat, 2, min)
        mat[, 1] <- pts
        mat[, 2] <- n1
        mat[, 3] <- n2
        mat[, 7] <- pvals
        for (ij in 1:length(pts)) {
            dif <- sort(estmat1[, ij] - estmat2[, ij])
            dif <- elimna(dif)
            nbad <- length(dif)
            lo <- round(nbad * alpha/2)
            hi <- nbad - lo
            mat[ij, 5] <- dif[lo]
            mat[ij, 6] <- dif[hi]
        }
    }
    if (plotit) 
        runmean2g(x1, y1, x2, y2, fr = fr1, est = median, sm = T)
    list(output = mat)
  }

musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.