In statlab/permuter: permuter: Permutation Tests

Vignette Info

library(permuter)
data(survival)

Unstratified Analysis

g <- three_sample[, 1]
str <- 1
data <- three_sample[, -1]
stratum <- rep(1, dim(data)[1])

B <- 100
p <- dim(data)[2]
T <- array(0, dim = c((B + 1), p, str))

for (s in 1:str) {
    g.str <- g[stratum == s]
    nu <- unlist(apply(data[stratum == s, ], 2, function(x) {
        table(g[stratum == s], is.na(x))
    }))
    n1 <- (6 * p - length(nu))/3
    NU <- array(0, dim = c(3, 2, p))
    NU[, 1, 1:n1] <- nu[1:(3 * n1)]
    NU[, , -c(1:n1)] <- nu[-c(1:(3 * n1))]

    ni1 <- array(0, dim = c(3, p))
    S1 <- array(0, dim = c(3, p))
    S2 <- array(0, dim = c(3, p))
    D <- apply(data[stratum == s, ], 2, function(x) {
        ifelse(is.na(x), 0, x)
    })
    for (nn in 1:3) {
        ni1[nn, ] <- apply(NU, 3, function(x) {
            sqrt((sum(x[, 1]) - x[nn, 1])/x[nn, 1])
        })
        S1[nn, ] <- apply(D, 2, function(x) {
            sum(x[g.str == (nn - 1)])
        })
        S2[nn, ] <- apply(D, 2, function(x) {
            sum(x) - sum(x[g.str == (nn - 1)])
        })
    }
    ni2 <- 1/ni1
    T[1, , s] <- apply(((ni1 * S1 - ni2 * S2)^2), 2, sum)
}



for (bb in 2:(B + 1)) {
    for (s in 1:str) {
        data.star <- data[stratum == s, ]
        n <- dim(data.star)[1]
        data.star <- data.star[sample(1:n), ]

        g.str <- g[stratum == s]
        nu <- unlist(apply(data.star, 2, function(x) {
            table(g[stratum == s], is.na(x))
        }))
        n1 <- (6 * p - length(nu))/3
        NU <- array(0, dim = c(3, 2, p))
        NU[, 1, 1:n1] <- nu[1:(3 * n1)]
        NU[, , -c(1:n1)] <- nu[-c(1:(3 * n1))]
        ni1 <- array(0, dim = c(3, p))
        S1 <- array(0, dim = c(3, p))
        S2 <- array(0, dim = c(3, p))
        D <- apply(data.star, 2, function(x) {
            ifelse(is.na(x), 0, x)
        })

        for (nn in 1:3) {
            ni1[nn, ] <- apply(NU, 3, function(x) {
                sqrt((sum(x[, 1]) - x[nn, 1])/x[nn, 1])
            })
            S1[nn, ] <- apply(D, 2, function(x) {
                sum(x[g.str == (nn - 1)])
            })
            S2[nn, ] <- apply(D, 2, function(x) {
                sum(x) - sum(x[g.str == (nn - 1)])
            })
        }
        ni2 <- 1/ni1
        T[bb, , s] <- apply(((ni1 * S1 - ni2 * S2)^2), 2, sum)
    }
}


P <- t2p_old(T)
T1 <- apply(P, 1, function(x) {
    -2 * log(prod(x))
})
t2p_old(T1)[1, ]

Plots

library(survival)

census <- function(x) {
    p <- length(x)
    time <- seq(1:p)
    t <- min(time[x == 0 | is.na(x) == TRUE])
    if (t == Inf) {
        t <- p
    }
    status <- ifelse(sum(is.na(x)) > 0, 0, 1)
    return(c(t, status))
}

par(mfrow = c(1, str))
options(warn = -1)

for (s in 1:str) {
    data.str <- data[stratum == s, ]
    g.str <- g[stratum == s]
    D <- t(apply(data.str, 1, census))
    D <- data.frame(g.str, D)
    colnames(D)[2:3] <- c("time", "status")
    plot(survfit(Surv(D$time, D$status) ~ g.str), col = c(2, 3, 4), xlab = "Time", 
        ylab = "S(t)", conf.int = FALSE)
    title(main = paste("Stratum: ", s, sep = ""))
}

Stratified Analysis

data(three_sample_strata)
g <- three_sample_strata[, 1]
stratum <- three_sample_strata[, 2]
str <- length(unique(stratum))
data <- three_sample_strata[, -c(1, 2)]

B <- 100

p <- dim(data)[2]
T <- array(0, dim = c((B + 1), p, str))

for (s in 1:str) {
    g.str <- g[stratum == s]
    nu <- unlist(apply(data[stratum == s, ], 2, function(x) {
        table(g[stratum == s], is.na(x))
    }))
    n1 <- (6 * p - length(nu))/3
    NU <- array(0, dim = c(3, 2, p))
    NU[, 1, 1:n1] <- nu[1:(3 * n1)]
    NU[, , -c(1:n1)] <- nu[-c(1:(3 * n1))]
    ni1 <- array(0, dim = c(3, p))
    S1 <- array(0, dim = c(3, p))
    S2 <- array(0, dim = c(3, p))
    D <- apply(data[stratum == s, ], 2, function(x) {
        ifelse(is.na(x), 0, x)
    })
    for (nn in 1:3) {
        ni1[nn, ] <- apply(NU, 3, function(x) {
            sqrt((sum(x[, 1]) - x[nn, 1])/x[nn, 1])
        })
        S1[nn, ] <- apply(D, 2, function(x) {
            sum(x[g.str == (nn - 1)])
        })
        S2[nn, ] <- apply(D, 2, function(x) {
            sum(x) - sum(x[g.str == (nn - 1)])
        })
    }
    ni2 <- 1/ni1
    T[1, , s] <- apply(((ni1 * S1 - ni2 * S2)^2), 2, sum)
}


for (bb in 2:(B + 1)) {
    for (s in 1:str) {
        data.star <- data[stratum == s, ]
        n <- dim(data.star)[1]
        data.star <- data.star[sample(1:n), ]
        g.str <- g[stratum == s]
        nu <- unlist(apply(data.star, 2, function(x) {
            table(g[stratum == s], is.na(x))
        }))
        n1 <- (6 * p - length(nu))/3
        NU <- array(0, dim = c(3, 2, p))
        NU[, 1, 1:n1] <- nu[1:(3 * n1)]
        NU[, , -c(1:n1)] <- nu[-c(1:(3 * n1))]
        ni1 <- array(0, dim = c(3, p))
        S1 <- array(0, dim = c(3, p))
        S2 <- array(0, dim = c(3, p))
        D <- apply(data.star, 2, function(x) {
            ifelse(is.na(x), 0, x)
        })
        for (nn in 1:3) {
            ni1[nn, ] <- apply(NU, 3, function(x) {
                sqrt((sum(x[, 1]) - x[nn, 1])/x[nn, 1])
            })
            S1[nn, ] <- apply(D, 2, function(x) {
                sum(x[g.str == (nn - 1)])
            })
            S2[nn, ] <- apply(D, 2, function(x) {
                sum(x) - sum(x[g.str == (nn - 1)])
            })
        }
        ni2 <- 1/ni1
        T[bb, , s] <- apply(((ni1 * S1 - ni2 * S2)^2), 2, sum)
    }
}

P <- t2p_old(T)
T1 <- apply(P, c(1, 3), function(x) {
    -2 * log(prod(x))
})
t2p_old(T1)[1, ]