reg1wayMC:

In musto101/wilcox_R: What the package does (short line)

Usage

reg1wayMC(x, y, regfun = tsreg, nboot = 100, SEED = TRUE, xout = FALSE, outfun = outpro, STAND = TRUE, alpha = 0.05, pr = TRUE, AD = FALSE, ...)

Arguments

x
y
regfun
nboot
SEED
xout
outfun
STAND
alpha
pr
AD
...

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, regfun = tsreg, nboot = 100, SEED = TRUE, xout = FALSE, 
    outfun = outpro, STAND = TRUE, alpha = 0.05, pr = TRUE, AD = FALSE, 
    ...) 
{
    library(parallel)
    if (pr) {
        if (!xout) 
            print("Might want to consider xout=T to  remove leverage points")
    }
    if (SEED) 
        set.seed(2)
    if (!is.list(x)) 
        stop("Argument x should have list mode")
    if (pr) {
        if (xout) 
            print("xout=T, so an adjusted critical is not computed and apparently not needed")
    }
    J = length(x)
    x = lapply(x, as.matrix)
    pchk = lapply(x, ncol)
    temp = matl(pchk)
    if (var(as.vector(temp)) != 0) 
        stop("Something is wrong. Number of covariates differs among the groups being compared")
    nv = NULL
    nv.keep = NULL
    nv.all = NULL
    p = ncol(x[[1]])
    p1 = p + 1
    for (j in 1:J) {
        xy = elimna(cbind(x[[j]], y[[j]]))
        x[[j]] = xy[, 1:p]
        y[[j]] = xy[, p1]
        x[[j]] = as.matrix(x[[j]])
        nv.all[j] = c(nv, nrow(x[[j]]))
    }
    if (xout) {
        temp = lapply(x, outfun, plotit = FALSE, STAND = STAND, 
            ...)
        for (j in 1:J) {
            x[[j]] = x[[j]][temp[[j]]$keep, ]
            y[[j]] = y[[j]][temp[[j]]$keep]
        }
    }
    x = lapply(x, as.matrix)
    p1 = ncol(x[[1]]) + 1
    K = p1
    est = matrix(NA, nrow = J, ncol = p1)
    hlabs = NULL
    vlabs = "Intercept"
    for (j in 1:J) hlabs[j] = paste("Group", j)
    for (j in 2:p1) vlabs[j] = paste("Slope", j - 1)
    dimnames(est) <- list(hlabs, vlabs)
    nv = NA
    ecov = list()
    ecovinv = list()
    W = rep(0, p1)
    gmean = rep(0, p1)
    for (j in 1:J) {
        est[j, ] = regfun(x[[j]], y[[j]])$coef
        nv.keep[j] = nrow(x[[j]])
        nv[j] = nv.keep[j]
        vals = matrix(NA, nrow = nboot, ncol = p1)
        data <- matrix(sample(nv[j], size = nv[j] * nboot, replace = TRUE), 
            ncol = nboot)
        data = listm(data)
        bvec <- mclapply(data, regbootMC, x[[j]], y[[j]], regfun, 
            ...)
        vals = t(matl(bvec))
        ecov[[j]] = var(vals)
        ecovinv[[j]] = solve(ecov[[j]])
        gmean = gmean + ecovinv[[j]] %*% est[j, ]
        W = W + ecovinv[[j]]
    }
    estall = solve(W) %*% gmean
    F = 0
    for (k in 1:K) {
        for (m in 1:K) {
            for (j in 1:J) {
                F = F + ecovinv[[j]][k, m] * (est[j, k] - estall[k]) * 
                  (est[j, m] - estall[m])
            }
        }
    }
    df = K * (J - 1)
    pvalad = NULL
    df = K * (J - 1)
    if (!xout || AD) {
        iden = diag(p1)
        Aw = 0
        for (j in 1:J) {
            temp = iden - solve(W) %*% ecovinv[[j]]
            tempsq = temp %*% temp
            Aw = Aw + (sum(diag(tempsq)) + (sum(diag(temp)))^2)/(nv[j] - 
                1)
        }
        Aw = Aw/2
        alval <- c(1:999)/1000
        for (i in 1:999) {
            irem <- i
            crit = qchisq(alval[i], df)
            critad = crit + (crit/(2 * df)) * (Aw + 3 * Aw * 
                crit/(df + 2))
            if (F < critad) 
                break
        }
        pavida = 1 - irem/1000
    }
    pval = 1 - pchisq(F, df)
    crit = qchisq(1 - alpha, df)
    critad = NULL
    if (AD) 
        critad = crit + (crit/(2 * df)) * (Aw + 3 * Aw * crit/(df + 
            2))
    est = data.frame(est)
    list(n = nv.all, n.keep = nv.keep, test.stat = F, crit.value = crit, 
        adjusted.crit = critad, p.value = pval, adjusted.p.value = pvalad, 
        est = est)
  }

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