ols1way2g:

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

Usage

ols1way2g(x, y, grp = c(1, 2), iv = 1, xout = FALSE, outfun = outpro, STAND = TRUE, alpha = 0.05, pr = TRUE, BLO = FALSE, ...)

Arguments

x
y
grp
iv
xout
outfun
STAND
alpha
pr
BLO
...

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, grp = c(1, 2), iv = 1, xout = FALSE, outfun = outpro, 
    STAND = TRUE, alpha = 0.05, pr = TRUE, BLO = FALSE, ...) 
{
    if (!is.list(x)) 
        stop("Argument x should have list mode")
    iv1 = iv + 1
    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
    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 = c(nv, nrow(x[[j]]))
    }
    nv.keep = nv
    critrad = NULL
    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]
        }
        if (BLO) {
            for (j in 1:J) {
                temp = reglev(x[[j]], y[[j]], plotit = FALSE)
                ad1 = c(temp1[[j]]$out.id, temp$regout)
                flag1 = duplicated(ad1)
                if (sum(flag1) > 0) {
                  flag1 = ad1[flag1]
                  x[[j]] = as.matrix(x[[j]])
                  x[[j]] = x[[j]][-flag1, ]
                  y[[j]] = y[[j]][-flag1]
                }
            }
        }
    }
    x = lapply(x, as.matrix)
    K = p1
    est = matrix(NA, nrow = J, ncol = p1)
    grpnum = NULL
    for (j in 1:J) grpnum[j] = paste("Group", j)
    vlabs = "Intercept"
    for (j in 2:p1) vlabs[j] = paste("Slope", j - 1)
    dimnames(est) = list(grpnum, vlabs)
    ecov = list()
    ecovinv = list()
    W = rep(0, p1)
    gmean = rep(0, K)
    for (j in 1:J) {
        est[j, ] = ols(x[[j]], y[[j]], xout = FALSE, plotit = FALSE, 
            ...)$coef
        nv.keep[j] = nrow(x[[j]])
        ecov[[j]] = olshc4(x[[j]], y[[j]])$cov
        ecovinv[[j]] = solve(ecov[[j]])
        gmean = gmean + ecovinv[[j]] %*% est[j, ]
        W = W + ecovinv[[j]]
    }
    estall = solve(W) %*% gmean
    F = 0
    for (j in 1:2) {
        F = F + ecovinv[[grp[j]]][iv1, iv1] * (est[grp[j], iv1] - 
            estall[iv1]) * (est[grp[j], iv1] - estall[iv1])
    }
    pvalad = NULL
    df = 1
    iden = 1
    Aw = 0
    for (j in 1:J) {
        temp = iden - solve(W[iv1, iv1]) %*% ecovinv[[grp[j]]][iv1, 
            iv1]
        tempsq = temp %*% temp
        Aw = Aw + (sum(diag(tempsq)) + (sum(diag(temp)))^2)/(nv[grp[j]] - 
            1)
    }
    Aw = Aw/2
    crit = qchisq(alpha, df)
    crit = crit + (crit/(2 * df)) * (Aw + 3 * Aw * crit/(df + 
        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
        pvalad = 1 - irem/1000
    }
    pval = 1 - pchisq(F, df)
    if (is.null(pvalad)) 
        pvalad = pval
    crit = qchisq(1 - alpha, df)
    critad = NULL
    critad = crit + (crit/(2 * df)) * (Aw + 3 * Aw * crit/(df + 
        2))
    est = data.frame(est)
    list(n = nv, 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.