Nothing
tests/testthat/test-test_only_hierarchy.R
In hierinf: Hierarchical Inference
require("testthat")
## random number generator
##### Change RNGversion in the future. Change Description such that #####
##### Depends: R (>= 3.6.0) #####
suppressWarnings(RNGversion("3.5.0"))
RNGkind("L'Ecuyer-CMRG")
test_that("test_only_hierarchy: check input", {
expect_error(test_only_hierarchy(x = NULL, y = NULL, dendr = NULL,
res.multisplit = NULL,
family = NULL),
"The response y is required to be a vector or a list of vectors if multiple data sets are present.")
expect_error(test_only_hierarchy(x = NULL, y = matrix(1:4, ncol = 2),
dendr = NULL, res.multisplit = NULL,
family = NULL),
"The elements of the list y are required to be numeric vectors or matrices with only one column. In the case of only one data set, it is enough that y is a numeric vector or matrix with only one column but it can as well be a list with one element.")
expect_error(test_only_hierarchy(x = NULL, y = list(matrix(1:4, ncol = 2),
matrix(1:4, ncol = 2)),
dendr = NULL, res.multisplit = NULL,
family = NULL),
"The elements of the list y are required to be numeric vectors or matrices with only one column.")
expect_error(test_only_hierarchy(x = NULL, y = 1:2, dendr = NULL,
res.multisplit = NULL, family = NULL),
"The input x is required to be a matrix or a list of matrices if multiple data sets are present.")
expect_error(test_only_hierarchy(x = matrix(1:4, ncol = 2), y = 1:2,
dendr = NULL, res.multisplit = NULL,
family = NULL),
"The matrix x is required to have column names. If there is no natural naming convention, then one can set them to some integer, say, 1 to p.",
fixed = TRUE)
tt <- matrix(1:4, ncol = 2)
colnames(tt) <- c("a", "b")
expect_error(test_only_hierarchy(x = tt, y = 1:2, dendr = NULL,
res.multisplit = NULL, family = NULL),
"The input res.multisplit is required to be a list.")
set.seed(88)
x <- matrix(rnorm(100), ncol = 2)
colnames(x) <- c("col1", "col2")
y <- 1:50
res.multisplit <- multisplit(x = x, y = y)
expect_error(test_only_hierarchy(x = x, y = y, dendr = NULL,
res.multisplit = res.multisplit,
family = NULL),
"The input dendr is required to be a list of dendrograms.")
# Please note that family = NULL results in taking the default value.
dendr <- cluster_var(x = x)
# expect_error(test_only_hierarchy(x = x, y = y, dendr = dendr,
# res.multisplit = res.multisplit,
# family = "alsdkk"),
# "'arg' should be one of \"gaussian\", \"binomial\"")
expected_result_1 <- data.frame(block = NA, p.value = NA,
significant.cluster = NA)
expected_result_1$significant.cluster <- list(NA)
attr(expected_result_1, "class") <- c("data.frame")
expected_result <- list(res.multisplit = res.multisplit,
res.hierarchy = expected_result_1)
attr(expected_result, "class") <- c("hierT", "list")
expect_equal(test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
expected_result)
# The column names of x or each element of x (list containing data sets)
# are required to have unique column names.
# (The tree and the output of the function multisplit are only fitted on x
# and y. This is irrelevant for this test.)
expect_error(test_only_hierarchy(x = cbind(x, x), y = y,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
"The matrix x is required to have unique column names.",
fixed = TRUE)
expect_error(test_only_hierarchy(x = list(cbind(x, x), x),
y = list(y, y),
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
"Each of the matrices which are stored in x are required to have unique column names.",
fixed = TRUE)
# If the trees (argument block is supplied) were fit on less variables,
# then we try to fit the model.
require("MASS")
p <- 20
n <- 80
B <- 50
sim.geno1 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno1) <- paste0("rsid", 1:p)
sim.geno2 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno2) <- paste0("rsid", 1:p)
dendr <- cluster_var(x = list(sim.geno1[, paste0("rsid", 1:10)],
sim.geno2[, paste0("rsid", 1:10)]),
block = data.frame(paste0("rsid", 1:10),
rep(1:2, each = 5),
stringsAsFactors = FALSE))
# plot(dendr$res.tree[[1]])
# plot(dendr$res.tree[[2]])
set.seed(144)
data.dim <- c(80, 20)
ind.active <- sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y1 <- sim.geno1 %*% beta + rnorm(data.dim[1])
y2 <- sim.geno2 %*% beta + rnorm(data.dim[1])
res.multisplit <- multisplit(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
family = "gaussian", B = B)
expect_error(test_only_hierarchy(x = list(sim.geno1, sim.geno2),
y = list(y1, y2), dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
"There are column name of x which have no corresponding values in the first column of block (column names of x).",
fixed = TRUE)
})
#### Check output with one data set ####
# This function is used to calculate the p-value of a given split for the
# binomial case.
MEL2 <- function(x, y, maxit, delta = 0.01, epsilon = 1e-6) {
# mean of y but we bound it awy from 0 and 1. See Equation (10) on page 8.
pi.hat <- max(delta, min(1 - delta, mean(y)))
# The two parameters delta.0 and delta.1 are constrained such that the average
# of the pseudo-observation is equal to pi.hat. See Equation (9) on page 8.
delta.0 <- (pi.hat * delta) / (1 + delta)
delta.1 <- (1 + pi.hat * delta) / (1 + delta)
# Pseudo-observations. See Equation (3) on page 4.
y.tilde <- delta.0 * (1 - y) + delta.1 * y
pseudo.y <- cbind(y.tilde, 1 - y.tilde)
# Suppress warning that "non-integer counts in a binomial glm!" because the
# function glm expects in the first column to be the number of successes and
# the second column to be the number of failures
suppressWarnings(res <- glm(pseudo.y ~ x, family = "binomial",
control = glm.control(epsilon = epsilon,
maxit = maxit),
model = FALSE, y = FALSE))
return(res)
}
# This function calculates the p-value for each of the notes in the tree.
check_test_hierarchy <- function(x, y, clvar, res.multisplit, B, cluster_test,
binomial = FALSE){
CT_colnames <- lapply(X = cluster_test, function(x) paste(x, collapse = "_"))
res <- matrix(NA, ncol = length(cluster_test), nrow = B)
colnames(res) <- CT_colnames
RES <- rep(NA, length(cluster_test))
names(RES) <- CT_colnames
for (i in cluster_test) { # for each cluster
for (b in 1:B) { # for each split (multi-sample splitting)
# selected coefficients
sel.coef <- res.multisplit[[1]]$sel.coef[b, ][!is.na(res.multisplit[[1]]$sel.coef[b, ])]
# other half of the samples
ind <- res.multisplit[[1]]$out.sample[b, ]
# combined data set (clvar plus x)
clvar_x <- cbind(clvar, x[, sel.coef, drop = FALSE])[ind, , drop = FALSE]
# intersection and set difference of selected coefficients and the given cluster
intersect_i <- intersect(sel.coef, i)
setdiff_i <- setdiff(sel.coef, i)
# columnnames for the model \hat{S}^{(b)} \setminus C
sel_i_clvar <- c(colnames(clvar), setdiff_i)
# browser()
# design matrix of the reduced model: variables of \hat{S}^{(b)} \setminus C & clvar
clvar_x_reduced <- clvar_x[, sel_i_clvar, drop = FALSE]
if (ncol(clvar_x_reduced) == 0) {
clvar_x_reduced <- rep(1, length(y[ind]))
}
res[b, paste(i, collapse = "_")] <-
if (length(intersect_i) == 0) { # Equation (2) on page 333 of Mandozzi and Buehlmann (2016)
1
} else {
if (binomial) {
# min(1, ... * |\hat{S}^{(b)}| / |\hat{S}^{(b)} \setminus C|) => Equation (2) & (3)
# on page 333 of Mandozzi and Buehlmann (2016)
min(1, anova(
# full model: variables of \hat{S}^{(b)} & clvar
MEL2(y = y[ind], x = clvar_x, maxit = 100),
# reduced model: variables of \hat{S}^{(b)} \setminus C & clvar
MEL2(y = y[ind], x = clvar_x_reduced, maxit = 100),
test = "Chisq")$"Pr(>Chi)"[2] *
length(sel.coef) / length(intersect_i)
)
} else {
# min(1, ... * |\hat{S}^{(b)}| / |\hat{S}^{(b)} \setminus C|) => Equation (2) & (3)
# on page 333 of Mandozzi and Buehlmann (2016)
min(1, anova(
# full model: variables of \hat{S}^{(b)} & clvar
lm(y[ind] ~ clvar_x),
# reduced model: variables of \hat{S}^{(b)} \setminus C & clvar
lm(y[ind] ~ clvar_x_reduced),
test = "F")$P[2] *
length(sel.coef) / length(intersect_i)
)
}
}
}
# Equation (4) on page 333 of Mandozzi and Buehlmann (2016)
RES[paste(i, collapse = "_")] <- adj_pval(res[, paste(i, collapse = "_")], B = B)
}
# hierarchical adjustment has to be done by hand (Equation below Equation (4))
return(RES)
}
adj_pval <- function(pvals, B) {
# define the sequence of gamma values
gamma_min <- 0.05
gamma_step <- 0.01
gamma_seq <- seq(gamma_min, 1, gamma_step)
# compute the empirical quantile vector
gamma_step <- vector("numeric", length = length(gamma_seq))
for (g in 1:length(gamma_seq)) {
gamma_step[g] <- min(1, quantile(pvals / gamma_seq[g], gamma_seq[g],
na.rm = TRUE))
}
# compute the adjusted p value
# Equation 4 on page 333 in Mandozzi and Buehlmann (2016)
return(min(1, (1 - log(gamma_min)) * min(gamma_step)))
}
### Example I ###
test_that("test_only_hierarchy: check output (Example I)", {
## simulate index
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y <- sim.geno %*% beta + rnorm(data.dim[1])
# cluster the data
dendr <- cluster_var(x = sim.geno)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = sim.geno, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = NULL,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-120)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-88)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-120)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-88)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example II ###
test_that("test_only_hierarchy: check output (Example II)", {
## simulate index
n <- 800
p <- 5
B <- 5
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
sim.clvar <- matrix(rnorm(n * 3), ncol = 3)
colnames(sim.clvar) <- paste0("clvar", 1:3)
set.seed(144)
data.dim <- dim(sim.geno) # first entry corresponds to rows and second to columns
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y <- sim.geno %*% beta + sim.clvar %*% c(0.25, 0.5, 1) + rnorm(data.dim[1])
# cluster the data
dendr <- cluster_var(x = sim.geno)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(555)
res.multisplit <- multisplit(x = sim.geno, y = y, clvar = sim.clvar,
family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y, clvar = sim.clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y, clvar = sim.clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", agg.method = "Stouffer")
## test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = sim.clvar,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-115)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-85)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-115)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-85)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Check output with multiple data sets ####
# Function for generating the data
require(MASS)
gen_one <- function(n, p, seed1, ind.a, seed3, num_clvar = NULL,
coef_clvar = NULL) {
set.seed(seed1)
x <- mvrnorm(n = n, mu = rep(0, p), Sigma = toeplitz(0.8^(seq(0, p - 1))) )
colnames(x) <- paste0("rsid", 1:p)
if (!is.null(num_clvar)) {
clvar <- matrix(rnorm(n * 3), ncol = 3)
colnames(clvar) <- paste0("clvar", 1:3)
} else {
clvar <- NULL
}
data.dim <- dim(x) # first entry corresponds to rows and second to columns
ind.active <- ind.a # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
set.seed(seed3)
if (!is.null(num_clvar)) {
y <- x %*% beta + rnorm(data.dim[1]) + clvar %*% coef_clvar
} else {
y <- x %*% beta + rnorm(data.dim[1])
}
return(list(x = x, y = y, clvar = clvar))
}
### Example III ###
test_that("test_only_hierarchy: check output (Example III multiple data sets)", {
skip_on_bioc()
## simulate index
n <- 800
p <- 5
B <- 50
## simulate data
r1 <- gen_one(n = n, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8)
r2 <- gen_one(n = n, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99)
r3 <- gen_one(n = n, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100)
r4 <- gen_one(n = n, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111)
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
# clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar)
# cluster the data
dendr <- cluster_var(x = x)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(744)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid4", "rsid1")])
expected_result$significant.cluster <- list(c("rsid4"), c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-145)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
# c(0.5, 0.5, 0.5, 0.5)
stouffer_weights <- sqrt(c(800, 800, 800, 800) / sum(c(800, 800, 800, 800)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid4", "rsid1")])
expected_result$significant.cluster <- list(c("rsid4"), c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-200)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example IV unbalanced data sets ###
test_that("test_only_hierarchy: check output (Example IV multiple data sets)", {
skip_on_bioc()
## simulate index
p <- 5
B <- 50
## simulate data
require(MASS)
r1 <- gen_one(n = 800, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8)
r2 <- gen_one(n = 200, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99)
r3 <- gen_one(n = 350, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100)
r4 <- gen_one(n = 50, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111)
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
# clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar)
# cluster the data
dendr <- cluster_var(x = x)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(6)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-132)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-135)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 200, 350, 50) / sum(c(800, 200, 350, 50)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-220)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-190)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example V with co-variables ###
###### The result of this test will change under RNGversion("3.6.0"). #####
###### TODO Adjust example calcualted by hand when switching to #####
###### RNGversion("3.6.0"). #####
test_that("test_only_hierarchy: check output (Example V multiple data sets)", {
skip_on_bioc()
## simulate index
p <- 5
B <- 50
## simulate data
require(MASS)
r1 <- gen_one(n = 800, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r2 <- gen_one(n = 200, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r3 <- gen_one(n = 350, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r4 <- gen_one(n = 50, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar) # with co-variables
# cluster the data
dendr <- cluster_var(x = x)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(3)
res.multisplit <- multisplit(x = x, y = y, clvar = clvar, family = "gaussian",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, clvar = clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, clvar = clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid4", "rsid5"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = clvar[[1]],
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = clvar[[2]],
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = clvar[[3]],
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = clvar[[4]],
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-134)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-138)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 200, 350, 50) / sum(c(800, 200, 350, 50)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-220)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-180)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
### Example VI, same data as in Example V but without clvar ###
# test hierarchy: no clvar but the response was created incl. clvar
# test_only_hierarchy: check output (Example VI multiple data sets)
# multisplit
set.seed(4)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = TRUE)
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = TRUE,
agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid4", "rsid5"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-70)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-80)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid3_rsid4_rsid5", "rsid1")])
expected_result$significant.cluster <- list(c("rsid3", "rsid4", "rsid5"), c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-120)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-100)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example VII: similar as Example V but with a block ###
test_that("test_only_hierarchy: check output (Example V multiple data sets)", {
skip_on_bioc()
## simulate index
p <- 10
B <- 50
## simulate data
require(MASS)
r1 <- gen_one(n = 800, p = p, seed1 = 9229, ind.a = c(7, 3), seed3 = 8,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r2 <- gen_one(n = 200, p = p, seed1 = 929, ind.a = c(7, 3), seed3 = 99,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r3 <- gen_one(n = 350, p = p, seed1 = 99, ind.a = c(7, 3), seed3 = 100,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r4 <- gen_one(n = 50, p = p, seed1 = 9, ind.a = c(7, 3), seed3 = 1111,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar) # with co-variables
# Block
block <- data.frame("var.names" = paste0("rsid", 1:10),
"blocks" = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2),
stringsAsFactors = FALSE)
# cluster the data
dendr <- cluster_var(x = x, block = block)
# plot(dendr$res.tree[[1]])
# plot(dendr$res.tree[[2]])
# multisplit
set.seed(3)
res.multisplit <- multisplit(x = x, y = y, clvar = clvar, family = "gaussian",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit, family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit, family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5", "rsid6",
"rsid7", "rsid8", "rsid9", "rsid10"),
c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid6", "rsid7", "rsid8", "rsid9", "rsid10"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid2", "rsid3"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5",
c("rsid6", "rsid7", "rsid8"),
c("rsid9", "rsid10"),
c("rsid7", "rsid8"),
"rsid6",
"rsid7",
"rsid8",
"rsid9",
"rsid10")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = clvar[[1]],
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = clvar[[2]],
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = clvar[[3]],
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = clvar[[4]],
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c("1", "2"),
p.value = compare_with[c("rsid3", "rsid7")],
stringsAsFactors = FALSE)
expected_result$significant.cluster <- list(c("rsid3"), c("rsid7"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_identical(res.T$res.hierarchy$block, expected_result$block)
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-90)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-100)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Tippett with global = FALSE
res.T <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = FALSE)
expect_identical(res.T$res.hierarchy$block, expected_result$block)
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-90)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-100)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 200, 350, 50) / sum(c(800, 200, 350, 50)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c("1", "2"),
p.value = compare_with[c("rsid3", "rsid7")],
stringsAsFactors = FALSE)
expected_result$significant.cluster <- list(c("rsid3"), c("rsid7"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_identical(res.S$res.hierarchy$block, expected_result$block)
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-125)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-140)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer with global = FALSE
res.S <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit, family = "gaussian",
agg.method = "Stouffer", global.test = FALSE)
expect_identical(res.S$res.hierarchy$block, expected_result$block)
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-125)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-140)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing on a tree which is build on less variables ####
test_that("test_only_hierarchy: check output (Example with smaller tree)", {
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y <- sim.geno %*% beta + rnorm(data.dim[1])
# cluster the data
dendr <- cluster_var(x = sim.geno[, c(1, 2, 5)])
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = sim.geno, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = NULL,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA),
p.value = compare_with[c("rsid1")])
expected_result$significant.cluster <- list(c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# test_only_hierarchy: check output (Example with smaller tree; no global)
# no global test
res.T <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = FALSE)
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA),
p.value = compare_with[c("rsid1")])
expected_result$significant.cluster <- list(c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# test_only_hierarchy: check output (Example with smaller tree; no global)
# no global test
res.S <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = FALSE,
agg.method = "Stouffer")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing with three data sets not measuring all the same variables ####
test_that("test_only_hierarchy: check return object for multiple data sets not measuring the same variables", {
B <- 50
set.seed(938)
tt1 <- matrix(rnorm(40), ncol = 2)
tt2 <- matrix(rnorm(40), ncol = 2)
tt3 <- matrix(rnorm(40), ncol = 2)
colnames(tt1) <- c("c1", "c2")
colnames(tt2) <- c("c1", "c5")
colnames(tt3) <- c("c2", "c5")
set.seed(144)
ind.active <- 2 # sample(1:3, 1)
beta <- rep(0, 3)
beta[ind.active] <- 2
y1 <- tt1 %*% beta[c(1, 2)] + rnorm(20)
y2 <- tt2 %*% beta[c(1, 3)] + rnorm(20)
y3 <- tt3 %*% beta[c(2, 3)] + rnorm(20)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = list(tt1, tt2, tt3), d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = list(tt1, tt2, tt3), y = list(y1, y2, y3),
family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = list(tt1, tt2, tt3), y = list(y1, y2, y3),
dendr = res_d, res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = list(tt1, tt2, tt3), y = list(y1, y2, y3),
dendr = res_d, res.multisplit = res.multisplit,
family = "gaussian", agg.method = "Stouffer")
cluster_test <- list(c("c1", "c2", "c5"),
c("c1", "c2"),
"c1",
"c2",
"c5")
res1 <- check_test_hierarchy(x = tt1, y = y1, clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = tt2, y = y2, clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = tt3, y = y3, clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 3)
# expected_result <- data.frame(block = c(NA),
# p.value = NA)
# expected_result$significant.cluster <- list(NA)
# rownames(expected_result) <- NULL
# For RNGversion("3.6.0"), we need to apply the rule that the p-value can
# only increase by going from top to bottom through the tree.
expected_result <- data.frame(block = c(NA),
p.value = max(compare_with[c("c2")],
compare_with[c("c1_c2")],
compare_with[c("c1_c2_c5")]))
expected_result$significant.cluster <- list(c("c2"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-15)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(20, 20, 20) / sum(c(20, 20, 20)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 3, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA),
p.value = NA)
expected_result$significant.cluster <- list(NA)
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### TODO We would have to apply the hierarchical constraint to the result
#### calculated by hand in order to receive the exact same results.
#### Perform testing with a data set that contains colinear variables ####
test_that("test_only_hierarchy: check return object for data set containing colinear variables", {
skip_on_bioc()
B <- 50
n <- 200
p <- 500
set.seed(3)
x <- mvrnorm(n, mu = rep(0, p), Sigma = diag(p))
colnames(x) <- paste0("Var", 1:p)
beta <- rep(0, p)
beta[c(5, 20, 46)] <- 1
y <- x %*% beta + rnorm(n)
x <- cbind(x, x) # duplicate all the variables (so now 1000 variables)
colnames(x)[501:1000] <- paste0("D2_", colnames(x)[501:1000])
dendr1 <- cluster_var(x = x)
set.seed(68)
res.multisplit1 <- multisplit(x = x, y = y, B = B)
# # both variables are included
# for (i in seq_len(nrow(res.multisplit1[[1]]$sel.coef))) {
# print(i)
# print(c("Var5", "D2_Var5") %in% res.multisplit1[[1]]$sel.coef[i, ])
# print(c("Var20", "D2_Var20") %in% res.multisplit1[[1]]$sel.coef[i, ])
# print(c("Var46", "D2_Var46") %in% res.multisplit1[[1]]$sel.coef[i, ])
# }
# library(glmnet)
# fit.lasso <- glmnet(x = x, y = y)
# # unequal zero
# for (i in seq(0.01225, 1.11600, by = 0.005)) {
# print(i)
# coef.lasso <- coef(fit.lasso, s = i)
# print(c("Var5", "D2_Var5") %in% rownames(coef.lasso)[as.vector(coef.lasso) != 0])
# print(c("Var20", "D2_Var20") %in% rownames(coef.lasso)[as.vector(coef.lasso) != 0])
# print(c("Var46", "D2_Var46") %in% rownames(coef.lasso)[as.vector(coef.lasso) != 0])
# }
# # larer than 0.01 (some fixed constant)
# for (i in seq(0.01225, 1.11600, by = 0.005)) {
# print(i)
# coef.lasso <- coef(fit.lasso, s = i)
# print(c("Var5", "D2_Var5") %in% rownames(coef.lasso)[as.vector(coef.lasso) > 0.01])
# print(c("Var20", "D2_Var20") %in% rownames(coef.lasso)[as.vector(coef.lasso) > 0.01])
# print(c("Var46", "D2_Var46") %in% rownames(coef.lasso)[as.vector(coef.lasso) > 0.01])
# }
# plot(fit.lasso)
# sum(coef(fit.lasso, s = 0.01225) > 0.5)
# rownames(coef.lasso)[as.vector(coef(fit.lasso, s = 0.01225) > 0.01)]
# We need to apply the rule that the p-value can only increase by going
# from top to bottom through the tree. This is NOT done here!!!!
suppressWarnings(sign.clusters1 <- test_only_hierarchy(x = x, y = y, dendr = dendr1,
res.multisplit = res.multisplit1,
family = "gaussian"))
# attr(sign.clusters1$res.hierarchy, "warningMsgs")
cluster_test <- list(c("Var46", "D2_Var46"),
c("Var5", "D2_Var5"),
c("Var20", "D2_Var20"))
# # and the single variables
# c("Var46"), c("D2_Var46"),
# c("Var5"), c("D2_Var5"),
# c("Var20"), c("D2_Var20"))
expected_result <- check_test_hierarchy(x = x, y = y, clvar = NULL,
res.multisplit = res.multisplit1[1],
B = B, cluster_test = cluster_test)
expect_true(all(sign.clusters1$res.hierarchy$p.value >= expected_result))
expect_equal(sign.clusters1$res.hierarchy$significant.cluster,
cluster_test)
})
#### Perform testing with binary response ####
test_that("test_only_hierarchy: check return object for a data set with binary response", {
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
eta <- sim.geno %*% beta
pr <- 1 / (1 + exp(-eta))
y <- rbinom(n, 1, prob = pr)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = sim.geno, d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = sim.geno, y = y,
family = "binomial", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y,
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y,
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial", agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = NULL,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test,
binomial = TRUE)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-25)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-20)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-25)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-20)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing two data sets with binary response ####
test_that("test_only_hierarchy: check return object for two data sets with binary response", {
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno1 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno1) <- paste0("rsid", 1:p)
sim.geno2 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno2) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno1)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
eta1 <- sim.geno1 %*% beta
pr1 <- 1 / (1 + exp(-eta1))
y1 <- rbinom(n, 1, prob = pr1)
eta2 <- sim.geno2 %*% beta
pr2 <- 1 / (1 + exp(-eta2))
y2 <- rbinom(n, 1, prob = pr2)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = list(sim.geno1, sim.geno2), d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
family = "binomial", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial", agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = sim.geno1, y = y1, clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test,
binomial = TRUE)
res2 <- check_test_hierarchy(x = sim.geno2, y = y2, clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test,
binomial = TRUE)
pvals_to_be <- rbind(res1, res2)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 2)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 800) / sum(c(800, 800)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 2, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing three data sets with binary response and unequal sample size ####
test_that("test_only_hierarchy: check return object for three data sets with binary response", {
skip_on_bioc()
n1 <- 800
n2 <- 200
n3 <- 150
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno1 <- mvrnorm(n = n1, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno1) <- paste0("rsid", 1:p)
sim.geno2 <- mvrnorm(n = n2, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno2) <- paste0("rsid", 1:p)
sim.geno3 <- mvrnorm(n = n3, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno3) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno1)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
eta1 <- sim.geno1 %*% beta
pr1 <- 1 / (1 + exp(-eta1))
y1 <- rbinom(n1, 1, prob = pr1)
eta2 <- sim.geno2 %*% beta
pr2 <- 1 / (1 + exp(-eta2))
y2 <- rbinom(n2, 1, prob = pr2)
eta3 <- sim.geno3 %*% beta
pr3 <- 1 / (1 + exp(-eta3))
y3 <- rbinom(n3, 1, prob = pr3)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = list(sim.geno1, sim.geno2, sim.geno3), d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = list(sim.geno1, sim.geno2, sim.geno3),
y = list(y1, y2, y3), family = "binomial",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = list(sim.geno1, sim.geno2, sim.geno3),
y = list(y1, y2, y3), dendr = res_d,
res.multisplit = res.multisplit,
family = "binomial")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = list(sim.geno1, sim.geno2, sim.geno3),
y = list(y1, y2, y3), dendr = res_d,
res.multisplit = res.multisplit,
family = "binomial", agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = sim.geno1, y = y1, clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test,
binomial = TRUE)
res2 <- check_test_hierarchy(x = sim.geno2, y = y2, clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test,
binomial = TRUE)
res3 <- check_test_hierarchy(x = sim.geno3, y = y3, clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test,
binomial = TRUE)
pvals_to_be <- rbind(res1, res2, res3)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 3)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(n1, n2, n3) / sum(c(n1, n2, n3)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 2, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Test parallel computations for multicore and sown ####
test_that("Parallel computations (in test_only_hierarchy): Test parallel computations for multicore and sown", {
skip_on_bioc()
## simulate index
n <- 800
p <- 10
B <- 50
## simulate data
r1 <- gen_one(n = n, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8)
r2 <- gen_one(n = n, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99)
r3 <- gen_one(n = n, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100)
r4 <- gen_one(n = n, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111)
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
# clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar)
# Block
block <- data.frame("var.names" = paste0("rsid", 1:10),
"blocks" = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2),
stringsAsFactors = FALSE)
### parallel = "multicore"
# cluster the data
dendr <- cluster_var(x = x, block = block, parallel = "multicore", ncpus = 2)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(744)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B,
parallel = "multicore", ncpus = 2)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", parallel = "multicore",
ncpus = 2)
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(1, 1, 2),
p.value = c(compare_with[c("rsid1", "rsid4")], NA))
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"), c(NA))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-40)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
### parallel = "snow"
# cluster the data
dendr <- cluster_var(x = x, block = block, parallel = "snow", ncpus = 2)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(744)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B,
parallel = "snow", ncpus = 2)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", parallel = "snow",
ncpus = 2)
## Test
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-40)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
Try the hierinf package in your browser
Any scripts or data that you put into this service are public.
hierinf documentation built on Nov. 8, 2020, 7:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
require("testthat")
## random number generator
##### Change RNGversion in the future. Change Description such that #####
##### Depends: R (>= 3.6.0) #####
suppressWarnings(RNGversion("3.5.0"))
RNGkind("L'Ecuyer-CMRG")
test_that("test_only_hierarchy: check input", {
expect_error(test_only_hierarchy(x = NULL, y = NULL, dendr = NULL,
res.multisplit = NULL,
family = NULL),
"The response y is required to be a vector or a list of vectors if multiple data sets are present.")
expect_error(test_only_hierarchy(x = NULL, y = matrix(1:4, ncol = 2),
dendr = NULL, res.multisplit = NULL,
family = NULL),
"The elements of the list y are required to be numeric vectors or matrices with only one column. In the case of only one data set, it is enough that y is a numeric vector or matrix with only one column but it can as well be a list with one element.")
expect_error(test_only_hierarchy(x = NULL, y = list(matrix(1:4, ncol = 2),
matrix(1:4, ncol = 2)),
dendr = NULL, res.multisplit = NULL,
family = NULL),
"The elements of the list y are required to be numeric vectors or matrices with only one column.")
expect_error(test_only_hierarchy(x = NULL, y = 1:2, dendr = NULL,
res.multisplit = NULL, family = NULL),
"The input x is required to be a matrix or a list of matrices if multiple data sets are present.")
expect_error(test_only_hierarchy(x = matrix(1:4, ncol = 2), y = 1:2,
dendr = NULL, res.multisplit = NULL,
family = NULL),
"The matrix x is required to have column names. If there is no natural naming convention, then one can set them to some integer, say, 1 to p.",
fixed = TRUE)
tt <- matrix(1:4, ncol = 2)
colnames(tt) <- c("a", "b")
expect_error(test_only_hierarchy(x = tt, y = 1:2, dendr = NULL,
res.multisplit = NULL, family = NULL),
"The input res.multisplit is required to be a list.")
set.seed(88)
x <- matrix(rnorm(100), ncol = 2)
colnames(x) <- c("col1", "col2")
y <- 1:50
res.multisplit <- multisplit(x = x, y = y)
expect_error(test_only_hierarchy(x = x, y = y, dendr = NULL,
res.multisplit = res.multisplit,
family = NULL),
"The input dendr is required to be a list of dendrograms.")
# Please note that family = NULL results in taking the default value.
dendr <- cluster_var(x = x)
# expect_error(test_only_hierarchy(x = x, y = y, dendr = dendr,
# res.multisplit = res.multisplit,
# family = "alsdkk"),
# "'arg' should be one of \"gaussian\", \"binomial\"")
expected_result_1 <- data.frame(block = NA, p.value = NA,
significant.cluster = NA)
expected_result_1$significant.cluster <- list(NA)
attr(expected_result_1, "class") <- c("data.frame")
expected_result <- list(res.multisplit = res.multisplit,
res.hierarchy = expected_result_1)
attr(expected_result, "class") <- c("hierT", "list")
expect_equal(test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
expected_result)
# The column names of x or each element of x (list containing data sets)
# are required to have unique column names.
# (The tree and the output of the function multisplit are only fitted on x
# and y. This is irrelevant for this test.)
expect_error(test_only_hierarchy(x = cbind(x, x), y = y,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
"The matrix x is required to have unique column names.",
fixed = TRUE)
expect_error(test_only_hierarchy(x = list(cbind(x, x), x),
y = list(y, y),
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
"Each of the matrices which are stored in x are required to have unique column names.",
fixed = TRUE)
# If the trees (argument block is supplied) were fit on less variables,
# then we try to fit the model.
require("MASS")
p <- 20
n <- 80
B <- 50
sim.geno1 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno1) <- paste0("rsid", 1:p)
sim.geno2 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno2) <- paste0("rsid", 1:p)
dendr <- cluster_var(x = list(sim.geno1[, paste0("rsid", 1:10)],
sim.geno2[, paste0("rsid", 1:10)]),
block = data.frame(paste0("rsid", 1:10),
rep(1:2, each = 5),
stringsAsFactors = FALSE))
# plot(dendr$res.tree[[1]])
# plot(dendr$res.tree[[2]])
set.seed(144)
data.dim <- c(80, 20)
ind.active <- sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y1 <- sim.geno1 %*% beta + rnorm(data.dim[1])
y2 <- sim.geno2 %*% beta + rnorm(data.dim[1])
res.multisplit <- multisplit(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
family = "gaussian", B = B)
expect_error(test_only_hierarchy(x = list(sim.geno1, sim.geno2),
y = list(y1, y2), dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian"),
"There are column name of x which have no corresponding values in the first column of block (column names of x).",
fixed = TRUE)
})
#### Check output with one data set ####
# This function is used to calculate the p-value of a given split for the
# binomial case.
MEL2 <- function(x, y, maxit, delta = 0.01, epsilon = 1e-6) {
# mean of y but we bound it awy from 0 and 1. See Equation (10) on page 8.
pi.hat <- max(delta, min(1 - delta, mean(y)))
# The two parameters delta.0 and delta.1 are constrained such that the average
# of the pseudo-observation is equal to pi.hat. See Equation (9) on page 8.
delta.0 <- (pi.hat * delta) / (1 + delta)
delta.1 <- (1 + pi.hat * delta) / (1 + delta)
# Pseudo-observations. See Equation (3) on page 4.
y.tilde <- delta.0 * (1 - y) + delta.1 * y
pseudo.y <- cbind(y.tilde, 1 - y.tilde)
# Suppress warning that "non-integer counts in a binomial glm!" because the
# function glm expects in the first column to be the number of successes and
# the second column to be the number of failures
suppressWarnings(res <- glm(pseudo.y ~ x, family = "binomial",
control = glm.control(epsilon = epsilon,
maxit = maxit),
model = FALSE, y = FALSE))
return(res)
}
# This function calculates the p-value for each of the notes in the tree.
check_test_hierarchy <- function(x, y, clvar, res.multisplit, B, cluster_test,
binomial = FALSE){
CT_colnames <- lapply(X = cluster_test, function(x) paste(x, collapse = "_"))
res <- matrix(NA, ncol = length(cluster_test), nrow = B)
colnames(res) <- CT_colnames
RES <- rep(NA, length(cluster_test))
names(RES) <- CT_colnames
for (i in cluster_test) { # for each cluster
for (b in 1:B) { # for each split (multi-sample splitting)
# selected coefficients
sel.coef <- res.multisplit[[1]]$sel.coef[b, ][!is.na(res.multisplit[[1]]$sel.coef[b, ])]
# other half of the samples
ind <- res.multisplit[[1]]$out.sample[b, ]
# combined data set (clvar plus x)
clvar_x <- cbind(clvar, x[, sel.coef, drop = FALSE])[ind, , drop = FALSE]
# intersection and set difference of selected coefficients and the given cluster
intersect_i <- intersect(sel.coef, i)
setdiff_i <- setdiff(sel.coef, i)
# columnnames for the model \hat{S}^{(b)} \setminus C
sel_i_clvar <- c(colnames(clvar), setdiff_i)
# browser()
# design matrix of the reduced model: variables of \hat{S}^{(b)} \setminus C & clvar
clvar_x_reduced <- clvar_x[, sel_i_clvar, drop = FALSE]
if (ncol(clvar_x_reduced) == 0) {
clvar_x_reduced <- rep(1, length(y[ind]))
}
res[b, paste(i, collapse = "_")] <-
if (length(intersect_i) == 0) { # Equation (2) on page 333 of Mandozzi and Buehlmann (2016)
1
} else {
if (binomial) {
# min(1, ... * |\hat{S}^{(b)}| / |\hat{S}^{(b)} \setminus C|) => Equation (2) & (3)
# on page 333 of Mandozzi and Buehlmann (2016)
min(1, anova(
# full model: variables of \hat{S}^{(b)} & clvar
MEL2(y = y[ind], x = clvar_x, maxit = 100),
# reduced model: variables of \hat{S}^{(b)} \setminus C & clvar
MEL2(y = y[ind], x = clvar_x_reduced, maxit = 100),
test = "Chisq")$"Pr(>Chi)"[2] *
length(sel.coef) / length(intersect_i)
)
} else {
# min(1, ... * |\hat{S}^{(b)}| / |\hat{S}^{(b)} \setminus C|) => Equation (2) & (3)
# on page 333 of Mandozzi and Buehlmann (2016)
min(1, anova(
# full model: variables of \hat{S}^{(b)} & clvar
lm(y[ind] ~ clvar_x),
# reduced model: variables of \hat{S}^{(b)} \setminus C & clvar
lm(y[ind] ~ clvar_x_reduced),
test = "F")$P[2] *
length(sel.coef) / length(intersect_i)
)
}
}
}
# Equation (4) on page 333 of Mandozzi and Buehlmann (2016)
RES[paste(i, collapse = "_")] <- adj_pval(res[, paste(i, collapse = "_")], B = B)
}
# hierarchical adjustment has to be done by hand (Equation below Equation (4))
return(RES)
}
adj_pval <- function(pvals, B) {
# define the sequence of gamma values
gamma_min <- 0.05
gamma_step <- 0.01
gamma_seq <- seq(gamma_min, 1, gamma_step)
# compute the empirical quantile vector
gamma_step <- vector("numeric", length = length(gamma_seq))
for (g in 1:length(gamma_seq)) {
gamma_step[g] <- min(1, quantile(pvals / gamma_seq[g], gamma_seq[g],
na.rm = TRUE))
}
# compute the adjusted p value
# Equation 4 on page 333 in Mandozzi and Buehlmann (2016)
return(min(1, (1 - log(gamma_min)) * min(gamma_step)))
}
### Example I ###
test_that("test_only_hierarchy: check output (Example I)", {
## simulate index
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y <- sim.geno %*% beta + rnorm(data.dim[1])
# cluster the data
dendr <- cluster_var(x = sim.geno)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = sim.geno, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = NULL,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-120)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-88)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-120)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-88)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example II ###
test_that("test_only_hierarchy: check output (Example II)", {
## simulate index
n <- 800
p <- 5
B <- 5
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
sim.clvar <- matrix(rnorm(n * 3), ncol = 3)
colnames(sim.clvar) <- paste0("clvar", 1:3)
set.seed(144)
data.dim <- dim(sim.geno) # first entry corresponds to rows and second to columns
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y <- sim.geno %*% beta + sim.clvar %*% c(0.25, 0.5, 1) + rnorm(data.dim[1])
# cluster the data
dendr <- cluster_var(x = sim.geno)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(555)
res.multisplit <- multisplit(x = sim.geno, y = y, clvar = sim.clvar,
family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y, clvar = sim.clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y, clvar = sim.clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", agg.method = "Stouffer")
## test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = sim.clvar,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-115)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-85)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-115)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-85)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Check output with multiple data sets ####
# Function for generating the data
require(MASS)
gen_one <- function(n, p, seed1, ind.a, seed3, num_clvar = NULL,
coef_clvar = NULL) {
set.seed(seed1)
x <- mvrnorm(n = n, mu = rep(0, p), Sigma = toeplitz(0.8^(seq(0, p - 1))) )
colnames(x) <- paste0("rsid", 1:p)
if (!is.null(num_clvar)) {
clvar <- matrix(rnorm(n * 3), ncol = 3)
colnames(clvar) <- paste0("clvar", 1:3)
} else {
clvar <- NULL
}
data.dim <- dim(x) # first entry corresponds to rows and second to columns
ind.active <- ind.a # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
set.seed(seed3)
if (!is.null(num_clvar)) {
y <- x %*% beta + rnorm(data.dim[1]) + clvar %*% coef_clvar
} else {
y <- x %*% beta + rnorm(data.dim[1])
}
return(list(x = x, y = y, clvar = clvar))
}
### Example III ###
test_that("test_only_hierarchy: check output (Example III multiple data sets)", {
skip_on_bioc()
## simulate index
n <- 800
p <- 5
B <- 50
## simulate data
r1 <- gen_one(n = n, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8)
r2 <- gen_one(n = n, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99)
r3 <- gen_one(n = n, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100)
r4 <- gen_one(n = n, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111)
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
# clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar)
# cluster the data
dendr <- cluster_var(x = x)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(744)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid4", "rsid1")])
expected_result$significant.cluster <- list(c("rsid4"), c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-145)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
# c(0.5, 0.5, 0.5, 0.5)
stouffer_weights <- sqrt(c(800, 800, 800, 800) / sum(c(800, 800, 800, 800)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid4", "rsid1")])
expected_result$significant.cluster <- list(c("rsid4"), c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-200)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example IV unbalanced data sets ###
test_that("test_only_hierarchy: check output (Example IV multiple data sets)", {
skip_on_bioc()
## simulate index
p <- 5
B <- 50
## simulate data
require(MASS)
r1 <- gen_one(n = 800, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8)
r2 <- gen_one(n = 200, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99)
r3 <- gen_one(n = 350, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100)
r4 <- gen_one(n = 50, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111)
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
# clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar)
# cluster the data
dendr <- cluster_var(x = x)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(6)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-132)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-135)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 200, 350, 50) / sum(c(800, 200, 350, 50)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-220)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-190)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example V with co-variables ###
###### The result of this test will change under RNGversion("3.6.0"). #####
###### TODO Adjust example calcualted by hand when switching to #####
###### RNGversion("3.6.0"). #####
test_that("test_only_hierarchy: check output (Example V multiple data sets)", {
skip_on_bioc()
## simulate index
p <- 5
B <- 50
## simulate data
require(MASS)
r1 <- gen_one(n = 800, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r2 <- gen_one(n = 200, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r3 <- gen_one(n = 350, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r4 <- gen_one(n = 50, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar) # with co-variables
# cluster the data
dendr <- cluster_var(x = x)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(3)
res.multisplit <- multisplit(x = x, y = y, clvar = clvar, family = "gaussian",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, clvar = clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, clvar = clvar,
dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid4", "rsid5"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = clvar[[1]],
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = clvar[[2]],
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = clvar[[3]],
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = clvar[[4]],
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-134)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-138)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 200, 350, 50) / sum(c(800, 200, 350, 50)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-220)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-180)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
### Example VI, same data as in Example V but without clvar ###
# test hierarchy: no clvar but the response was created incl. clvar
# test_only_hierarchy: check output (Example VI multiple data sets)
# multisplit
set.seed(4)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = TRUE)
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = TRUE,
agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid4", "rsid5"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-70)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-80)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid3_rsid4_rsid5", "rsid1")])
expected_result$significant.cluster <- list(c("rsid3", "rsid4", "rsid5"), c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-120)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-100)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
### Example VII: similar as Example V but with a block ###
test_that("test_only_hierarchy: check output (Example V multiple data sets)", {
skip_on_bioc()
## simulate index
p <- 10
B <- 50
## simulate data
require(MASS)
r1 <- gen_one(n = 800, p = p, seed1 = 9229, ind.a = c(7, 3), seed3 = 8,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r2 <- gen_one(n = 200, p = p, seed1 = 929, ind.a = c(7, 3), seed3 = 99,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r3 <- gen_one(n = 350, p = p, seed1 = 99, ind.a = c(7, 3), seed3 = 100,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
r4 <- gen_one(n = 50, p = p, seed1 = 9, ind.a = c(7, 3), seed3 = 1111,
num_clvar = 3, coef_clvar = c(0.5, 0.25, 1.25))
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar) # with co-variables
# Block
block <- data.frame("var.names" = paste0("rsid", 1:10),
"blocks" = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2),
stringsAsFactors = FALSE)
# cluster the data
dendr <- cluster_var(x = x, block = block)
# plot(dendr$res.tree[[1]])
# plot(dendr$res.tree[[2]])
# multisplit
set.seed(3)
res.multisplit <- multisplit(x = x, y = y, clvar = clvar, family = "gaussian",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit, family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit, family = "gaussian",
agg.method = "Stouffer")
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5", "rsid6",
"rsid7", "rsid8", "rsid9", "rsid10"),
c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid6", "rsid7", "rsid8", "rsid9", "rsid10"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid2", "rsid3"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5",
c("rsid6", "rsid7", "rsid8"),
c("rsid9", "rsid10"),
c("rsid7", "rsid8"),
"rsid6",
"rsid7",
"rsid8",
"rsid9",
"rsid10")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = clvar[[1]],
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = clvar[[2]],
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = clvar[[3]],
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = clvar[[4]],
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c("1", "2"),
p.value = compare_with[c("rsid3", "rsid7")],
stringsAsFactors = FALSE)
expected_result$significant.cluster <- list(c("rsid3"), c("rsid7"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_identical(res.T$res.hierarchy$block, expected_result$block)
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-90)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-100)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Tippett with global = FALSE
res.T <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = FALSE)
expect_identical(res.T$res.hierarchy$block, expected_result$block)
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-90)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-100)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 200, 350, 50) / sum(c(800, 200, 350, 50)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 4, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c("1", "2"),
p.value = compare_with[c("rsid3", "rsid7")],
stringsAsFactors = FALSE)
expected_result$significant.cluster <- list(c("rsid3"), c("rsid7"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_identical(res.S$res.hierarchy$block, expected_result$block)
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-125)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-140)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer with global = FALSE
res.S <- test_only_hierarchy(x = x, y = y, clvar = clvar, dendr = dendr,
res.multisplit = res.multisplit, family = "gaussian",
agg.method = "Stouffer", global.test = FALSE)
expect_identical(res.S$res.hierarchy$block, expected_result$block)
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-125)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-140)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing on a tree which is build on less variables ####
test_that("test_only_hierarchy: check output (Example with smaller tree)", {
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
y <- sim.geno %*% beta + rnorm(data.dim[1])
# cluster the data
dendr <- cluster_var(x = sim.geno[, c(1, 2, 5)])
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = sim.geno, y = y, family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian",
agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = NULL,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA),
p.value = compare_with[c("rsid1")])
expected_result$significant.cluster <- list(c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# test_only_hierarchy: check output (Example with smaller tree; no global)
# no global test
res.T <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = FALSE)
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA),
p.value = compare_with[c("rsid1")])
expected_result$significant.cluster <- list(c("rsid1"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# test_only_hierarchy: check output (Example with smaller tree; no global)
# no global test
res.S <- test_only_hierarchy(x = sim.geno, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", global.test = FALSE,
agg.method = "Stouffer")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing with three data sets not measuring all the same variables ####
test_that("test_only_hierarchy: check return object for multiple data sets not measuring the same variables", {
B <- 50
set.seed(938)
tt1 <- matrix(rnorm(40), ncol = 2)
tt2 <- matrix(rnorm(40), ncol = 2)
tt3 <- matrix(rnorm(40), ncol = 2)
colnames(tt1) <- c("c1", "c2")
colnames(tt2) <- c("c1", "c5")
colnames(tt3) <- c("c2", "c5")
set.seed(144)
ind.active <- 2 # sample(1:3, 1)
beta <- rep(0, 3)
beta[ind.active] <- 2
y1 <- tt1 %*% beta[c(1, 2)] + rnorm(20)
y2 <- tt2 %*% beta[c(1, 3)] + rnorm(20)
y3 <- tt3 %*% beta[c(2, 3)] + rnorm(20)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = list(tt1, tt2, tt3), d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = list(tt1, tt2, tt3), y = list(y1, y2, y3),
family = "gaussian", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = list(tt1, tt2, tt3), y = list(y1, y2, y3),
dendr = res_d, res.multisplit = res.multisplit,
family = "gaussian")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = list(tt1, tt2, tt3), y = list(y1, y2, y3),
dendr = res_d, res.multisplit = res.multisplit,
family = "gaussian", agg.method = "Stouffer")
cluster_test <- list(c("c1", "c2", "c5"),
c("c1", "c2"),
"c1",
"c2",
"c5")
res1 <- check_test_hierarchy(x = tt1, y = y1, clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = tt2, y = y2, clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = tt3, y = y3, clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 3)
# expected_result <- data.frame(block = c(NA),
# p.value = NA)
# expected_result$significant.cluster <- list(NA)
# rownames(expected_result) <- NULL
# For RNGversion("3.6.0"), we need to apply the rule that the p-value can
# only increase by going from top to bottom through the tree.
expected_result <- data.frame(block = c(NA),
p.value = max(compare_with[c("c2")],
compare_with[c("c1_c2")],
compare_with[c("c1_c2_c5")]))
expected_result$significant.cluster <- list(c("c2"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-15)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(20, 20, 20) / sum(c(20, 20, 20)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 3, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA),
p.value = NA)
expected_result$significant.cluster <- list(NA)
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-120)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### TODO We would have to apply the hierarchical constraint to the result
#### calculated by hand in order to receive the exact same results.
#### Perform testing with a data set that contains colinear variables ####
test_that("test_only_hierarchy: check return object for data set containing colinear variables", {
skip_on_bioc()
B <- 50
n <- 200
p <- 500
set.seed(3)
x <- mvrnorm(n, mu = rep(0, p), Sigma = diag(p))
colnames(x) <- paste0("Var", 1:p)
beta <- rep(0, p)
beta[c(5, 20, 46)] <- 1
y <- x %*% beta + rnorm(n)
x <- cbind(x, x) # duplicate all the variables (so now 1000 variables)
colnames(x)[501:1000] <- paste0("D2_", colnames(x)[501:1000])
dendr1 <- cluster_var(x = x)
set.seed(68)
res.multisplit1 <- multisplit(x = x, y = y, B = B)
# # both variables are included
# for (i in seq_len(nrow(res.multisplit1[[1]]$sel.coef))) {
# print(i)
# print(c("Var5", "D2_Var5") %in% res.multisplit1[[1]]$sel.coef[i, ])
# print(c("Var20", "D2_Var20") %in% res.multisplit1[[1]]$sel.coef[i, ])
# print(c("Var46", "D2_Var46") %in% res.multisplit1[[1]]$sel.coef[i, ])
# }
# library(glmnet)
# fit.lasso <- glmnet(x = x, y = y)
# # unequal zero
# for (i in seq(0.01225, 1.11600, by = 0.005)) {
# print(i)
# coef.lasso <- coef(fit.lasso, s = i)
# print(c("Var5", "D2_Var5") %in% rownames(coef.lasso)[as.vector(coef.lasso) != 0])
# print(c("Var20", "D2_Var20") %in% rownames(coef.lasso)[as.vector(coef.lasso) != 0])
# print(c("Var46", "D2_Var46") %in% rownames(coef.lasso)[as.vector(coef.lasso) != 0])
# }
# # larer than 0.01 (some fixed constant)
# for (i in seq(0.01225, 1.11600, by = 0.005)) {
# print(i)
# coef.lasso <- coef(fit.lasso, s = i)
# print(c("Var5", "D2_Var5") %in% rownames(coef.lasso)[as.vector(coef.lasso) > 0.01])
# print(c("Var20", "D2_Var20") %in% rownames(coef.lasso)[as.vector(coef.lasso) > 0.01])
# print(c("Var46", "D2_Var46") %in% rownames(coef.lasso)[as.vector(coef.lasso) > 0.01])
# }
# plot(fit.lasso)
# sum(coef(fit.lasso, s = 0.01225) > 0.5)
# rownames(coef.lasso)[as.vector(coef(fit.lasso, s = 0.01225) > 0.01)]
# We need to apply the rule that the p-value can only increase by going
# from top to bottom through the tree. This is NOT done here!!!!
suppressWarnings(sign.clusters1 <- test_only_hierarchy(x = x, y = y, dendr = dendr1,
res.multisplit = res.multisplit1,
family = "gaussian"))
# attr(sign.clusters1$res.hierarchy, "warningMsgs")
cluster_test <- list(c("Var46", "D2_Var46"),
c("Var5", "D2_Var5"),
c("Var20", "D2_Var20"))
# # and the single variables
# c("Var46"), c("D2_Var46"),
# c("Var5"), c("D2_Var5"),
# c("Var20"), c("D2_Var20"))
expected_result <- check_test_hierarchy(x = x, y = y, clvar = NULL,
res.multisplit = res.multisplit1[1],
B = B, cluster_test = cluster_test)
expect_true(all(sign.clusters1$res.hierarchy$p.value >= expected_result))
expect_equal(sign.clusters1$res.hierarchy$significant.cluster,
cluster_test)
})
#### Perform testing with binary response ####
test_that("test_only_hierarchy: check return object for a data set with binary response", {
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
eta <- sim.geno %*% beta
pr <- 1 / (1 + exp(-eta))
y <- rbinom(n, 1, prob = pr)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = sim.geno, d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = sim.geno, y = y,
family = "binomial", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = sim.geno, y = y,
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = sim.geno, y = y,
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial", agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
pvals_to_be <- check_test_hierarchy(x = sim.geno, y = y, clvar = NULL,
res.multisplit = res.multisplit,
B = B, cluster_test = cluster_test,
binomial = TRUE)
# Tippett
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-25)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-20)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
compare_with <- pvals_to_be
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-25)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-20)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing two data sets with binary response ####
test_that("test_only_hierarchy: check return object for two data sets with binary response", {
n <- 800
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno1 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno1) <- paste0("rsid", 1:p)
sim.geno2 <- mvrnorm(n = n, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno2) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno1)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
eta1 <- sim.geno1 %*% beta
pr1 <- 1 / (1 + exp(-eta1))
y1 <- rbinom(n, 1, prob = pr1)
eta2 <- sim.geno2 %*% beta
pr2 <- 1 / (1 + exp(-eta2))
y2 <- rbinom(n, 1, prob = pr2)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = list(sim.geno1, sim.geno2), d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
family = "binomial", B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = list(sim.geno1, sim.geno2), y = list(y1, y2),
dendr = res_d, res.multisplit = res.multisplit,
family = "binomial", agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = sim.geno1, y = y1, clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test,
binomial = TRUE)
res2 <- check_test_hierarchy(x = sim.geno2, y = y2, clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test,
binomial = TRUE)
pvals_to_be <- rbind(res1, res2)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 2)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(800, 800) / sum(c(800, 800)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 2, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Perform testing three data sets with binary response and unequal sample size ####
test_that("test_only_hierarchy: check return object for three data sets with binary response", {
skip_on_bioc()
n1 <- 800
n2 <- 200
n3 <- 150
p <- 5
B <- 50
## simulate data
require(MASS)
set.seed(9229)
sim.geno1 <- mvrnorm(n = n1, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno1) <- paste0("rsid", 1:p)
sim.geno2 <- mvrnorm(n = n2, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno2) <- paste0("rsid", 1:p)
sim.geno3 <- mvrnorm(n = n3, mu = rep(0, p),
Sigma = toeplitz(0.8^(seq(0, p - 1))))
colnames(sim.geno3) <- paste0("rsid", 1:p)
set.seed(144)
data.dim <- dim(sim.geno1)
ind.active <- c(4, 1) # sample(1:data.dim[2], 2)
beta <- rep(0, data.dim[2])
beta[ind.active] <- 2
eta1 <- sim.geno1 %*% beta
pr1 <- 1 / (1 + exp(-eta1))
y1 <- rbinom(n1, 1, prob = pr1)
eta2 <- sim.geno2 %*% beta
pr2 <- 1 / (1 + exp(-eta2))
y2 <- rbinom(n2, 1, prob = pr2)
eta3 <- sim.geno3 %*% beta
pr3 <- 1 / (1 + exp(-eta3))
y3 <- rbinom(n3, 1, prob = pr3)
# use = "pairwise.complete.obs" (default)
res_d <- cluster_var(x = list(sim.geno1, sim.geno2, sim.geno3), d = NULL,
method = "average",
block = NULL)
# plot(res_d$res.tree[[1]])
# multisplit
set.seed(2)
res.multisplit <- multisplit(x = list(sim.geno1, sim.geno2, sim.geno3),
y = list(y1, y2, y3), family = "binomial",
B = B)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = list(sim.geno1, sim.geno2, sim.geno3),
y = list(y1, y2, y3), dendr = res_d,
res.multisplit = res.multisplit,
family = "binomial")
# test hierarchy: Stouffer
res.S <- test_only_hierarchy(x = list(sim.geno1, sim.geno2, sim.geno3),
y = list(y1, y2, y3), dendr = res_d,
res.multisplit = res.multisplit,
family = "binomial", agg.method = "Stouffer")
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2"),
c("rsid3", "rsid4", "rsid5"),
c("rsid3", "rsid4"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = sim.geno1, y = y1, clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test,
binomial = TRUE)
res2 <- check_test_hierarchy(x = sim.geno2, y = y2, clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test,
binomial = TRUE)
res3 <- check_test_hierarchy(x = sim.geno3, y = y3, clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test,
binomial = TRUE)
pvals_to_be <- rbind(res1, res2, res3)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 3)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.T$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
# Stouffer
stouffer_weights <- sqrt(c(n1, n2, n3) / sum(c(n1, n2, n3)))
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y, stouffer_weights) {
pnorm(sum(stouffer_weights * qnorm(x)))
},
len_y = 2, stouffer_weights = stouffer_weights)
expected_result <- data.frame(block = c(NA, NA),
p.value = compare_with[c("rsid1", "rsid4")])
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.S$res.hierarchy$p.value[1], expected_result$p.value[1],
tol = 1e-40)
expect_equal(res.S$res.hierarchy$p.value[2], expected_result$p.value[2],
tol = 1e-30)
expect_equal(res.S$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
#### Test parallel computations for multicore and sown ####
test_that("Parallel computations (in test_only_hierarchy): Test parallel computations for multicore and sown", {
skip_on_bioc()
## simulate index
n <- 800
p <- 10
B <- 50
## simulate data
r1 <- gen_one(n = n, p = p, seed1 = 9229, ind.a = c(4, 1), seed3 = 8)
r2 <- gen_one(n = n, p = p, seed1 = 929, ind.a = c(4, 1), seed3 = 99)
r3 <- gen_one(n = n, p = p, seed1 = 99, ind.a = c(4, 1), seed3 = 100)
r4 <- gen_one(n = n, p = p, seed1 = 9, ind.a = c(4, 1), seed3 = 1111)
x <- list(r1$x, r2$x, r3$x, r4$x)
y <- list(r1$y, r2$y, r3$y, r4$y)
# clvar <- list(r1$clvar, r2$clvar, r3$clvar, r4$clvar)
# Block
block <- data.frame("var.names" = paste0("rsid", 1:10),
"blocks" = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2),
stringsAsFactors = FALSE)
### parallel = "multicore"
# cluster the data
dendr <- cluster_var(x = x, block = block, parallel = "multicore", ncpus = 2)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(744)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B,
parallel = "multicore", ncpus = 2)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", parallel = "multicore",
ncpus = 2)
## Test
# This list encodes the tree structure
cluster_test <- list(c("rsid1", "rsid2", "rsid3", "rsid4", "rsid5"),
c("rsid1", "rsid2", "rsid3"),
c("rsid4", "rsid5"),
c("rsid1", "rsid2"),
"rsid1",
"rsid2",
"rsid3",
"rsid4",
"rsid5")
res1 <- check_test_hierarchy(x = x[[1]], y = y[[1]], clvar = NULL,
res.multisplit = res.multisplit[1],
B = B, cluster_test = cluster_test)
res2 <- check_test_hierarchy(x = x[[2]], y = y[[2]], clvar = NULL,
res.multisplit = res.multisplit[2],
B = B, cluster_test = cluster_test)
res3 <- check_test_hierarchy(x = x[[3]], y = y[[3]], clvar = NULL,
res.multisplit = res.multisplit[3],
B = B, cluster_test = cluster_test)
res4 <- check_test_hierarchy(x = x[[4]], y = y[[4]], clvar = NULL,
res.multisplit = res.multisplit[4],
B = B, cluster_test = cluster_test)
pvals_to_be <- rbind(res1, res2, res3, res4)
# Tippett
compare_with <- apply(X = pvals_to_be, MARGIN = 2,
FUN = function(x, len_y) {
max(1 - (1 - min(x))^(len_y), .Machine$double.neg.eps)
},
len_y = 4)
expected_result <- data.frame(block = c(1, 1, 2),
p.value = c(compare_with[c("rsid1", "rsid4")], NA))
expected_result$significant.cluster <- list(c("rsid1"), c("rsid4"), c(NA))
rownames(expected_result) <- NULL
attr(expected_result, "class") <- c("data.frame")
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-40)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
### parallel = "snow"
# cluster the data
dendr <- cluster_var(x = x, block = block, parallel = "snow", ncpus = 2)
# plot(dendr$res.tree[[1]])
# multisplit
set.seed(744)
res.multisplit <- multisplit(x = x, y = y, family = "gaussian", B = B,
parallel = "snow", ncpus = 2)
# test hierarchy: Tippett
res.T <- test_only_hierarchy(x = x, y = y, dendr = dendr,
res.multisplit = res.multisplit,
family = "gaussian", parallel = "snow",
ncpus = 2)
## Test
expect_equal(res.T$res.hierarchy$p.value, expected_result$p.value,
tol = 1e-40)
expect_equal(res.T$res.hierarchy$significant.cluster,
expected_result$significant.cluster)
})
Try the hierinf package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.