library(permuter) data(survival)
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, ]
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 = "")) }
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, ]
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.