context("Testing statistics")
library(amp)
test_that("summary statistics work without grouping", {
smry <- summary_statistics(mark_nas(example_set, 0),
grouping_cols = NA)
ex <- exprs(mark_nas(example_set, 0))
for (fun in c("finite_mean", "finite_sd", "finite_median", "finite_mad")) {
expect_equal(unname(apply(ex, 1, fun)), smry[, gsub("finite_", "", fun)])
}
expect_equal(unname(apply(ex, 1, finite_quantile, probs = 0.25)), smry$Q25)
expect_equal(unname(apply(ex, 1, finite_quantile, probs = 0.75)), smry$Q75)
})
test_that("summary statistics work with grouping", {
smry <- summary_statistics(mark_nas(example_set, 0))
exa <- exprs(mark_nas(example_set[, example_set$Group == "A"], 0))
for (fun in c("finite_mean", "finite_sd", "finite_median", "finite_mad")) {
expect_equal(unname(apply(exa, 1, fun)), smry[, gsub("finite", "Group_A", fun)])
}
expect_equal(unname(apply(exa, 1, finite_quantile, probs = 0.25)), smry$Group_A_Q25)
expect_equal(unname(apply(exa, 1, finite_quantile, probs = 0.75)), smry$Group_A_Q75)
exb <- exprs(mark_nas(example_set[, example_set$Group == "B"], 0))
for (fun in c("finite_mean", "finite_sd", "finite_median", "finite_mad")) {
expect_equal(unname(apply(exb, 1, fun)), smry[, gsub("finite", "Group_B", fun)])
}
expect_equal(unname(apply(exb, 1, finite_quantile, probs = 0.25)), smry$Group_B_Q25)
expect_equal(unname(apply(exb, 1, finite_quantile, probs = 0.75)), smry$Group_B_Q75)
})
test_that("summary statistics work with all NA features", {
ex_set_na <- mark_nas(example_set, 0)
exprs(ex_set_na)[1, ] <- NA
smry <- summary_statistics(ex_set_na)
expect_equal(nrow(smry), nrow(exprs(ex_set_na)))
expect_equal(smry$Feature_ID, featureNames(ex_set_na))
expect_true(all(is.na(smry[1, 2:ncol(smry)])))
})
test_that("Cohen's d works", {
ex <- example_set %>%
drop_qcs() %>%
mark_nas(0)
cd <- ex %>%
combined_data()
cd1 <- cd[cd$Time == 1, ]
cd2 <- cd[cd$Time == 2, ]
d <- c()
for (feature in featureNames(ex)) {
tdiff <- data.frame(feature = cd2[, feature] - cd1[, feature],
group = cd1$Group)
mean_a <- finite_mean(tdiff$feature[tdiff$group == "A"])
sd_a <- finite_sd(tdiff$feature[tdiff$group == "A"])
mean_b <- finite_mean(tdiff$feature[tdiff$group == "B"])
sd_b <- finite_sd(tdiff$feature[tdiff$group == "B"])
d <- c(d, (mean_b - mean_a)/mean(c(sd_a, sd_b)))
}
cohd <- cohens_d(ex)
df <- data.frame(Feature_ID = featureNames(ex),
Cohen_d = d,
stringsAsFactors = FALSE)
rownames(df) <- df$Feature_ID
expect_equal(cohd, df)
})
test_that("Cohen's d work with all NA features", {
ex_set_na <- drop_qcs(mark_nas(example_set, 0))
exprs(ex_set_na)[1:2, ] <- NA
cohd <- cohens_d(ex_set_na)
expect_equal(nrow(cohd), nrow(exprs(ex_set_na)))
expect_equal(cohd$Feature_ID, featureNames(ex_set_na))
expect_true(all(is.na(cohd[1:2, 2:ncol(cohd)])))
})
test_that("Cohen's d is not run with multiple time or group levels", {
expect_error(cohens_d(example_set), "Column Group")
ex <- example_set
ex$Group <- c(1,2)
expect_error(cohens_d(ex), "Column Time")
ex$Time <- c(1,2)
ex$Group <- 1
expect_error(cohens_d(ex), "Column Group")
})
test_that("Fold change works", {
ex <- example_set %>%
mark_nas(0)
cd <- ex %>%
combined_data()
cd1 <- cd[cd$Time == 1, ]
cd2 <- cd[cd$Time == 2, ]
fc <- data.frame(Feature_ID = featureNames(ex),
FC_B_vs_A = 1,
FC_QC_vs_A = 1,
FC_QC_vs_B = 1,
stringsAsFactors = FALSE)
rownames(fc) <- fc$Feature_ID
for (i in seq_len(nrow(fc))) {
feature <- fc$Feature_ID[i]
mean_a <- finite_mean(cd[cd$Group == "A", feature])
mean_b <- finite_mean(cd[cd$Group == "B", feature])
mean_qc <- finite_mean(cd[cd$Group == "QC", feature])
fc$FC_B_vs_A[i] <- mean_b/mean_a
fc$FC_QC_vs_A[i] <- mean_qc/mean_a
fc$FC_QC_vs_B[i] <- mean_qc/mean_b
}
foldc <- fold_change(ex)
expect_equal(foldc, fc)
})
test_that("Fold change works with all NA features", {
ex_set_na <- drop_qcs(mark_nas(example_set, 0))
exprs(ex_set_na)[1:2, ] <- NA
foldc <- fold_change(ex_set_na)
expect_equal(nrow(foldc), nrow(exprs(ex_set_na)))
expect_equal(foldc$Feature_ID, featureNames(ex_set_na))
expect_true(all(is.na(foldc[1:2, 2:ncol(foldc)])))
})
test_that("P-value correction works", {
ps <- data.frame(x = letters,
x_P = runif(26),
P_x = runif(26),
y_P = runif(26))
adj <- adjust_p_values(ps, flags = rep(NA, 26))
expect_equal(adj$x_P_FDR, p.adjust(ps$x_P, method = "BH"))
expect_equal(adj[colnames(ps)], ps)
expect_equal(ncol(adj), ncol(ps) + 2)
adj2 <- adjust_p_values(ps, flags = c("a", "a", "a",
rep(NA, 23)))
expect_equal(adj2$x_P_FDR,
c(rep(NA, 3), p.adjust(ps$x_P[4:26], method = "BH")))
})
test_that("Linear model works", {
cd <- combined_data(drop_qcs(example_set))
lm_fit <- lm(HILIC_neg_259_9623a4_4322 ~ Time,
data = cd)
smry <- summary(lm_fit)
# Works for a simple example
lm_res <- perform_lm(drop_qcs(example_set),
formula_char = "Feature ~ Time")
expect_equal(lm_res$Time2_Estimate[1], smry$coefficients[2,1])
expect_equal(lm_res$Time2_Std_Error[1], smry$coefficients[2,2])
expect_equal(lm_res$Time2_t_value[1], smry$coefficients[2,3])
expect_equal(lm_res$Time2_P[1], smry$coefficients[2,4])
expect_equal(lm_res$R2[1], smry$r.squared)
expect_equal(lm_res$Adj_R2[1], smry$adj.r.squared)
# Works with column with only NA values
ex_set_na <- drop_qcs(mark_nas(example_set, 0))
exprs(ex_set_na)[1:2, ] <- NA
lm_res <- perform_lm(ex_set_na,
formula_char = "Feature ~ Time")
expect_equal(nrow(lm_res), nrow(exprs(example_set)))
expect_equal(lm_res$Feature_ID, featureNames(example_set))
expect_true(all(is.na(lm_res[1:2, 2:ncol(lm_res)])))
# FDR correction ignored for flagged compounds
ex_set_na <- flag_quality(ex_set_na)
lm_res2 <- perform_lm(ex_set_na, formula_char = "Feature ~ Group")
flag_idx <- !is.na(flag(ex_set_na))
expect_true(all(is.na(lm_res2$Group_P_FDR[flag_idx])))
})
test_that("Logistic regression works", {
cd <- combined_data(drop_qcs(example_set))
glm_fit <- glm(Group ~ HILIC_neg_259_9623a4_4322,
data = cd,
family = binomial())
smry <- summary(glm_fit)
glm_res <- perform_logistic(drop_qcs(example_set),
formula_char = "Group ~ Feature")
expect_equal(glm_res$Feature_Estimate[1], smry$coefficients[2,1])
expect_equal(glm_res$Feature_Std_Error[1], smry$coefficients[2,2])
expect_equal(glm_res$Feature_z_value[1], smry$coefficients[2,3])
expect_equal(glm_res$Feature_P[1], smry$coefficients[2,4])
expect_equal(glm_res$Feature_LCI95[1], confint(glm_fit)[2,1])
expect_equal(glm_res$Feature_UCI95[1], confint(glm_fit)[2,2])
ex_set_na <- drop_qcs(mark_nas(example_set, 0))
exprs(ex_set_na)[1:2, ] <- NA
glm_res <- perform_logistic(ex_set_na,
formula_char = "Group ~ Feature")
expect_equal(nrow(glm_res), nrow(exprs(example_set)))
expect_equal(glm_res$Feature_ID, featureNames(example_set))
expect_true(all(is.na(glm_res[1:2, 2:ncol(glm_res)])))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.