In statlab/permuter: permuter: Permutation Tests

Vignette Info

library(permuter)
data(survival)

Unstratified Analysis

g <- two_sample_surv[, 1]
str <- 1
data <- two_sample_surv[, -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 <- (4 * p - length(nu))/2
    NU <- array(0, dim = c(2, 2, p))
    NU[, 1, 1:n1] <- nu[1:(2 * n1)]
    NU[, , -c(1:n1)] <- nu[-c(1:(2 * n1))]

    ni1 <- apply(NU, 3, function(x) {
        sqrt(x[2, 1]/x[1, 1])
    })
    ni2 <- 1/ni1

    D <- apply(data[stratum == s, ], 2, function(x) {
        ifelse(is.na(x), 0, x)
    })
    D1 <- apply(D, 2, function(x) {
        sum(x[g.str == 0])
    })
    D2 <- apply(D, 2, function(x) {
        sum(x[g.str == 1])
    })
    T[1, , s] <- ni1 * D1 - ni2 * D2
}



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 <- (4 * p - length(nu))/2
        NU <- array(0, dim = c(2, 2, p))
        NU[, 1, 1:n1] <- nu[1:(2 * n1)]
        NU[, , -c(1:n1)] <- nu[-c(1:(2 * n1))]

        ni1 <- apply(NU, 3, function(x) {
            sqrt(x[2, 1]/x[1, 1])
        })
        ni2 <- 1/ni1

        D <- apply(data.star, 2, function(x) {
            ifelse(is.na(x), 0, x)
        })
        D1 <- apply(D, 2, function(x) {
            sum(x[g.str == 0])
        })
        D2 <- apply(D, 2, function(x) {
            sum(x[g.str == 1])
        })

        T[bb, , s] <- ni1 * D1 - ni2 * D2
    }
}

P <- t2p_old(abs(T))  ## two.sided
# P=t2p_old(-T) ## one.sided: H1 g1 > g0
T1 <- apply(P, 1, function(x) {
    -2 * log(prod(x))
})  # no strata 
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, 4), xlab = "Time", 
        ylab = "S(t)", conf.int = FALSE)
    title(main = paste("Stratum: ", s, sep = ""))
}

Stratified Analysis

g <- two_sample_surv_strata[, 1]
stratum <- two_sample_surv_strata[, 2]
str <- length(unique(stratum))
data <- two_sample_surv_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 <- (4 * p - length(nu))/2
    NU <- array(0, dim = c(2, 2, p))
    NU[, 1, 1:n1] <- nu[1:(2 * n1)]
    NU[, , -c(1:n1)] <- nu[-c(1:(2 * n1))]

    ni1 <- apply(NU, 3, function(x) {
        sqrt(x[2, 1]/x[1, 1])
    })
    ni2 <- 1/ni1

    D <- apply(data[stratum == s, ], 2, function(x) {
        ifelse(is.na(x), 0, x)
    })
    D1 <- apply(D, 2, function(x) {
        sum(x[g.str == 0])
    })
    D2 <- apply(D, 2, function(x) {
        sum(x[g.str == 1])
    })
    T[1, , s] <- ni1 * D1 - ni2 * D2
}



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 <- (4 * p - length(nu))/2
        NU <- array(0, dim = c(2, 2, p))
        NU[, 1, 1:n1] <- nu[1:(2 * n1)]
        NU[, , -c(1:n1)] <- nu[-c(1:(2 * n1))]

        ni1 <- apply(NU, 3, function(x) {
            sqrt(x[2, 1]/x[1, 1])
        })
        ni2 <- 1/ni1

        D <- apply(data.star, 2, function(x) {
            ifelse(is.na(x), 0, x)
        })
        D1 <- apply(D, 2, function(x) {
            sum(x[g.str == 0])
        })
        D2 <- apply(D, 2, function(x) {
            sum(x[g.str == 1])
        })

        T[bb, , s] <- ni1 * D1 - ni2 * D2
    }
}

P <- t2p_old(abs(T))  ## two.sided
# P=t2p_old(-T) ## one.sided: H1 g1 > g0
T1 <- apply(P, c(1, 3), function(x) {
    -2 * log(prod(x))
})
t2p_old(T1)[1, ]