library(testthat)
library(dplyr)
context("mincLm - two group test")
if(!exists("dataPath"))
dataPath <- tempdir()
getRMINCTestData(dataPath)
dataPath <- file.path(dataPath, "rminctestdata/")
gf <- read.csv(file.path(dataPath, "test_data_set.csv"))
voxel_left <- mincGetVoxel(gf$jacobians_fixed_2[1:10], 0,0,0)
voxel_right <- mincGetVoxel(gf$jacobians_fixed_2[11:20], 0,0,0)
Sex <- gf$Sex[1:10]
Scale <- gf$scale[1:10]
Coil <- as.factor(gf$coil[1:10])
gf$coil <- as.factor(gf$coil)
gftest <- gf[1:10,]
gftest$voxel_right <- (gf$jacobians_fixed_2[11:20])
gftest$voxel_left_file <- gf$jacobians_fixed_2[1:10]
rmod <- lm(voxel_left ~ Sex)
rLm = summary(rmod)
# silence the output of mincLm, in order to make the test output information is more clear to read
rmincLm = verboseRun("mincLm(voxel_left_file ~ Sex, gftest)",getOption("verbose"))
test_that("mincLm Two Factors",{
expect_equal(rmincLm[1,1:6]
, with(rLm
,c(fstatistic[1]
, r.squared[1]
, coefficients[,"Estimate"]
, coefficients[,"t value"])), check.attributes = FALSE)
expect_that(attr(rmincLm,"df")[[2]],is_equivalent_to(rLm$df[2]))
})
test_that("Likelihood and information criteria are computed correctly", {
expect_equal(as.numeric(rmincLm[1,"logLik"]), as.numeric(logLik(rmod)))
expect_equal(as.numeric(AIC(rmincLm)[1]), as.numeric(AIC(rmod)))
expect_equal(as.numeric(BIC(rmincLm)[1]), as.numeric(BIC(rmod)))
})
# now test that findPeaks works
context("mincFindPeaks - ensure we can find correct peaks")
verboseRun({
# get rid of NAs
rmincLm[is.na(rmincLm)] <- 0
# find peaks
peaks <- mincFindPeaks(rmincLm, "tvalue-SexM", minDistance = 1)
# find the min and max by hand
rmincLmArray <- mincArray(rmincLm, "tvalue-SexM")
maxPeak <- arrayInd(which.max(rmincLmArray), .dim=c(10,10,10))
minPeak <- arrayInd(which.min(rmincLmArray), .dim=c(10,10,10))
minPeakFromPeaks <- as.integer(arrange(peaks, value)[1,1:3])
maxPeakFromPeaks <- as.integer(peaks[1,1:3])
})
test_that("mincFindPeaks min and max", {
expect_that(maxPeak[1],is_equivalent_to(maxPeakFromPeaks[1]))
expect_that(maxPeak[2],is_equivalent_to(maxPeakFromPeaks[2]))
expect_that(maxPeak[3],is_equivalent_to(maxPeakFromPeaks[3]))
expect_that(minPeak[1],is_equivalent_to(minPeakFromPeaks[1]))
expect_that(minPeak[2],is_equivalent_to(minPeakFromPeaks[2]))
expect_that(minPeak[3],is_equivalent_to(minPeakFromPeaks[3]))
})
context("vertexFindPeaks - ensure matches mincFindPeaks")
test_that("vertexFindPeaks matches mincFindPeaks", {
adj <- RMINC:::neighbour_list(10,10,10, 6)
g <- igraph::graph.adjlist(lapply(adj, function(nebs) nebs + 1))
pos_peaks <-
vertexFindPeaks(rmincLm[,"tvalue-SexM"], g, direction = "positive", mindist = 1, output = "indices")
pos_peak <- arrayInd(pos_peaks[which.max(rmincLm[pos_peaks,"tvalue-SexM"])], .dim = rep(10,3))
neg_peaks <-
vertexFindPeaks(rmincLm[,"tvalue-SexM"], g, direction = "negative", mindist = 1, output = "indices")
neg_peak <- arrayInd(neg_peaks[which.min(rmincLm[neg_peaks,"tvalue-SexM"])], .dim = rep(10,3))
expect_equal(pos_peak[1:3], maxPeak[1:3])
expect_equal(neg_peak[1:3], minPeak[1:3])
})
context("mincLm - two group test with interaction")
# silence the output of mincLm, in order to make the test output information is more clear to read
rmincLm2 = verboseRun("mincLm(voxel_left_file~Sex*Scale,gftest)",getOption("verbose"))
gftest$voxel_left = voxel_left
rLm2 = summary(lm(voxel_left~Sex*Scale,gftest))
test_that("mincLm interaction",{
expect_equal(rmincLm2[1,1:10]
, with(rLm2
,c(fstatistic[1]
, r.squared[1]
, coefficients[,"Estimate"]
, coefficients[,"t value"])), check.attributes = FALSE)
expect_that(attr(rmincLm2,"df")[[2]],is_equivalent_to(rLm2$df[2]))
})
context("mincLm - three group test")
# silence the output of mincLm, in order to make the test output information is more clear to read
rmincLm3 <- verboseRun("mincLm(voxel_left_file~coil,gftest)",getOption("verbose"))
gftest$voxel_left = voxel_left
rLm3 = summary(lm(voxel_left~Coil,gftest))
test_that("mincLm Three Factors",{
expect_equal(rmincLm3[1,1:8]
, with(rLm3
,c(fstatistic[1]
, r.squared[1]
, coefficients[,"Estimate"]
, coefficients[,"t value"])), check.attributes = FALSE)
expect_that(attr(rmincLm3,"df")[[2]],is_equivalent_to(rLm3$df[2]))
})
context("mincLm - three group test with interaction")
# silence the output of mincLm, in order to make the test output information is more clear to read
rmincLm4 = verboseRun("mincLm(voxel_left_file~Scale*Coil,gftest)",getOption("verbose"))
gftest$voxel_left = voxel_left
rLm4 = summary(lm(voxel_left~Scale*Coil,gftest))
test_that("mincLm Three Factors Interaction",{
expect_equal(rmincLm4[1,1:14]
, with(rLm4
,c(fstatistic[1]
, r.squared[1]
, coefficients[,"Estimate"]
, coefficients[,"t value"])), check.attributes = FALSE)
expect_that(attr(rmincLm4,"df")[[2]],is_equivalent_to(rLm4$df[2]))
})
context("mincLm - Model Selection")
test_that("Model Selection Works", {
comp1 <- compare_models(rmincLm, rmincLm2, rmincLm3, rmincLm4, metric = AIC)
comp2 <- compare_models(rmincLm, rmincLm2, rmincLm3, rmincLm4)
params <- sapply(list(rmincLm, rmincLm2, rmincLm3, rmincLm4)
, function(mod) ncol(attr(mod, "model")) + 1)
n <- nrow(gftest)
aicc_corr <- 2*(params+1)*(params + 2) / (n - params - 2)
aicc_corr_mat <- t(matrix(rep(aicc_corr, nrow(comp1)), ncol = nrow(comp1)))
expect_equal(unclass(comp2), unclass(comp1 + aicc_corr_mat), check.attributes = FALSE)
expect_equal(as.numeric(summary(comp2)$wins), c(100, 100, 119, 981))
})
context("mincLm - Parallel")
test_that("mincLm local multicore works", {
skip_on_cran()
skip_on_travis()
if(Sys.getenv("TEST_Q_MINC") != "yes")
skip("qMinc tests disabled")
verboseRun(prlm <- mincLm(voxel_left_file~Scale*Coil,gftest, parallel = c("local", 3)))
expect_equal(rmincLm4, prlm, check.attributes = FALSE)
})
# test_that("mincLm queue parallel works", {
# verboseRun(qrlm <- mincLm(voxel_left_file~Scale*Coil,gftest, parallel = c("sge", 2)))
#
# expect_equal(rmincLm, qrlm, check.attributes = FALSE)
# })
## Need a much better test here.
context("mincLm - Randomize")
test_that("mincLm randomize works", {
verboseRun(rlm <- mincLm(voxel_left_file~Sex, gftest))
verboseRun(rand <- RMINC:::mincRandomize.mincLm(rlm, R = 10))
thresh <- thresholds(rand)
expect_equal(dim(thresholds(rand)), c(4,2))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.