options(mc.cores=21)
facedes <- read.table(system.file("extdata", "face_design.txt", package = "fmrireg"), header=TRUE)
facedes$repnum <- factor(facedes$rep_num)
gen_mask_file <- function(d, perc) {
arr = array(0,d)
vals <- ifelse(runif(prod(d)) > .5, 1, 0)
vol <- NeuroVol(vals, NeuroSpace(d))
fname <- paste0(tempfile(), ".nii")
write_vol(vol, fname)
fname
}
gen_fake_dataset <- function(d, nscans) {
onames <- vector(length=nscans, mode="list")
for (i in 1:nscans) {
arr <- array(rnorm(prod(d)), d)
bspace <- neuroim2::NeuroSpace(dim=d)
vec <- neuroim2::NeuroVec(arr, bspace)
fname <- paste0(tempfile(), ".nii")
write_vec(vec, fname)
onames[i] <- fname
}
onames
}
## test that latent and fmri_mem_dataset of same underlying latent dataset produce the same betas
test_that("can construct and run a simple fmri glm from in memory dataset", {
scans <- lapply(1:length(unique(facedes$run)), function(i) {
arr <- array(rnorm(10*10*10*244), c(10,10,10, 244))
bspace <- neuroim2::NeuroSpace(dim=c(10,10,10,244))
neuroim2::NeuroVec(arr, bspace)
})
mask <- neuroim2::LogicalNeuroVol(array(rnorm(10*10*10), c(10,10,10)) > 0, neuroim2::NeuroSpace(dim=c(10,10,10)))
#scans <- list.files("test_data/images_study/epi/", "rscan0.*nii", full.names=TRUE)
dset <- fmri_mem_dataset(scans=scans,
mask=mask,
TR=1.5,
event_table=facedes)
mod <- fmri_lm(onset ~ hrf(repnum), block = ~ run, dataset=dset, durations=0, strategy="chunkwise", nchunks=4)
expect_true(!is.null(mod))
})
test_that("can construct and run a simple fmri glm from in memory dataset and one contrast", {
scans <- lapply(1:length(unique(facedes$run)), function(i) {
arr <- array(rnorm(10*10*10*244), c(10,10,10, 244))
bspace <- neuroim2::NeuroSpace(dim=c(10,10,10,244))
neuroim2::NeuroVec(arr, bspace)
})
mask <- neuroim2::LogicalNeuroVol(array(rnorm(10*10*10), c(10,10,10)) > 0, neuroim2::NeuroSpace(dim=c(10,10,10)))
#scans <- list.files("test_data/images_study/epi/", "rscan0.*nii", full.names=TRUE)
dset <- fmri_mem_dataset(scans=scans,
mask=mask,
TR=1.5,
event_table=facedes)
con <<- contrast_set(pair_contrast( ~ repnum == 1, ~ repnum == 2, name="rep2_rep1"))
mod1 <- fmri_lm(onset ~ hrf(repnum, contrasts=con), block = ~ run, dataset=dset, durations=0)
mod1a <- fmri_lm(onset ~ hrf(repnum, contrasts=con), block = ~ run, dataset=dset, durations=0,
meta_weighting="equal")
mod2 <- fmri_lm(onset ~ hrf(repnum, contrasts=con), block = ~ run, dataset=dset, durations=0,
strategy="chunkwise", nchunks=10)
expect_true(!is.null(mod1))
expect_true(!is.null(mod1a))
expect_true(!is.null(mod2))
expect_equal(ncol(mod1$result$contrasts$estimate), 1)
expect_equal(ncol(mod1a$result$contrasts$estimate), 1)
expect_equal(ncol(mod2$result$contrasts$estimate), 1)
expect_equal(nrow(mod1$result$contrasts$estimate), nrow(mod2$result$contrasts$estimate))
c1 <- cor(mod1$result$contrasts$estimate[,1], mod2$result$contrasts$estimate[,1])
expect_true(c1> .97)
})
test_that("can construct and run a simple fmri glm from a matrix_dataset with 1 column", {
vals <- rep(rnorm(244),6)
dset <- matrix_dataset(as.matrix(vals),TR=1.5, run_length=rep(244,6), event_table=facedes)
c1 <- pair_contrast( ~ repnum == 1, ~ repnum == 2, name="rep2_rep1")
c2 <- pair_contrast( ~ repnum == 3, ~ repnum == 4, name="rep3_rep4")
con <<- contrast_set(c1,c2)
mod1 <- fmri_lm(onset ~ hrf(repnum, contrasts=con), block = ~ run, dataset=dset, durations=0)
mod2 <- fmri_lm(onset ~ hrf(repnum, contrasts=con), block = ~ run, dataset=dset, durations=0,
strategy="chunkwise", nchunks=1)
expect_true(!is.null(mod1))
expect_equal(ncol(mod1$result$contrasts$estimate), 2)
expect_equal(ncol(mod2$result$contrasts$estimate), 2)
})
test_that("fmri glm for multivariate matrix and complex contrast ", {
vals <- do.call(cbind, lapply(1:100, function(i) rnorm(244*6)))
fd <- subset(facedes, null == 0 & rt < 2)
fd$letter <- sample(factor(rep(letters[1:4], length.out=nrow(fd))))
dset <- matrix_dataset(vals,TR=1.5, run_length=rep(244,6), event_table=fd)
cset <<- contrast_set(pair_contrast( ~ letter %in% c("a", "b"),
~ letter %in% c("c", "d"),
name="abcd_efgh"),
pair_contrast( ~ letter %in% c("a", "c"),
~ letter %in% c("b", "d"),
name="ijkl_mnop"),
unit_contrast(~ letter, "letter"))
#c3 <- unit_contrast(~ letter, "letter")
# bmod <- baseline_model(basis="constant", degree=1, intercept="none", sframe=dset$sampling_frame)
mod1 <- fmri_lm(onset ~ hrf(letter, contrasts=cset),
#baseline_model=bmod,
block = ~ run, dataset=dset, durations=0, nchunks=1,strategy="chunkwise")
zz <- stats(mod1, "contrasts")
betas <- mod1$result$betas$estimate
expect_true(!is.null(mod1))
})
test_that("can construct and run a simple fmri glm from a matrix_dataset with 2 columns", {
vals <- cbind(rep(rnorm(244),6), rep(rnorm(244),6))
dset <- matrix_dataset(as.matrix(vals),TR=1.5, run_length=rep(244,6), event_table=facedes)
c1 <- pair_contrast( ~ repnum == 1, ~ repnum == 2, name="rep2_rep1")
c2 <- pair_contrast( ~ repnum == 2, ~ repnum == 3, name="rep3_rep2")
con <<- contrast_set(c1,c2)
mod1 <- fmri_lm(onset ~ hrf(repnum, contrasts=con), block = ~ run, dataset=dset, durations=0)
mod2 <- fmri_lm(onset ~ hrf(repnum, contrasts=con), block = ~ run, dataset=dset, durations=0,
strategy="chunkwise", nchunks=1)
expect_true(!is.null(mod1))
expect_true(!is.null(mod2))
expect_equal(ncol(mod1$result$contrasts$estimate), 2)
expect_equal(ncol(mod2$result$contrasts$estimate), 2)
})
test_that("can construct and run a simple fmri glm two terms and prefix args", {
vals <- cbind(rep(rnorm(244),6), rep(rnorm(244),6))
dset <- matrix_dataset(as.matrix(vals),TR=1.5, run_length=rep(244,6), event_table=facedes)
mod1 <- fmri_lm(onset ~ hrf(repnum, subset=repnum %in% c(1,2), prefix="r12")+
hrf(repnum, subset=repnum %in% c(3,4), prefix="r34"),
block = ~ run, dataset=dset, durations=0)
expect_true(!is.null(mod1))
#expect_true(!is.null(mod2))
expect_equal(ncol(mod1$result$betas$estimate), 4)
#expect_equal(ncol(mod2$result$contrasts$estimate()), 2)
})
test_that("can run video fmri design with matrix_dataset", {
des <- read.table(system.file("extdata", "video_design.txt", package = "fmrireg"), header=TRUE)
events <- rep(320,7)
sframe <- sampling_frame(rep(320, length(events)), TR=1.5)
evmod <- event_model(Onset ~ hrf(Video, Condition, basis="spmg1"),
block = ~ run, sampling_frame=sframe, data=des)
bmod <- baseline_model(basis="bs", degree=4, sframe=sframe)
fmod <- fmri_model(evmod, bmod)
dset <- matrix_dataset(matrix(rnorm(320*7*100), 320*7, 100),TR=1.5, run_length=rep(320,7), event_table=des)
#conset <- fmrireg::one_against_all_contrast(levels(des$Video), "Video")
conset <<- do.call("contrast_set", lapply(levels(factor(des$Video)), function(v) {
f1 <- as.formula(paste("~ Video == ", paste0('"', v, '"')))
f2 <- as.formula(paste("~ Video != ", paste0('"', v, '"')))
pair_contrast(f1, f2, name=paste0(v, "_vsall"))
}))
res1 <- fmrireg:::fmri_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec"), block= ~ run, dataset=dset,
strategy="runwise")
res2 <- fmrireg:::fmri_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec"), block= ~ run, dataset=dset,
strategy="chunkwise", nchunks=12)
res3 <- fmrireg:::fmri_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec"), block= ~ run, dataset=dset,
strategy="chunkwise", nchunks=1)
expect_true(!is.null(coef(res1)))
expect_true(!is.null(coef(res2)))
expect_true(!is.null(coef(res3)))
expect_true(!is.null(coef(res1, "contrasts")))
expect_true(!is.null(coef(res2, "contrasts")))
expect_true(!is.null(coef(res3, "contrasts")))
})
test_that("can run video fmri design with fmri_file_dataset", {
library(neuroim2)
des <- read.table(system.file("extdata", "video_design.txt", package = "fmrireg"), header=TRUE)
events <- rep(320,7)
sframe <- sampling_frame(rep(320, length(events)), TR=1.5)
scans <- gen_fake_dataset(c(10,10,10,320), 7)
maskfile <- gen_mask_file(c(10,10,10))
dset <- fmri_dataset(scans, maskfile,TR=1.5, rep(320,7), base_path="/", mode="bigvec", event_table=as_tibble(des))
evmod <- event_model(Onset ~ hrf(Video, Condition, basis="spmg1"),
block = ~ run, sampling_frame=sframe, data=des)
bmod <- baseline_model(basis="bs", degree=4, sframe=sframe)
fmod <- fmri_model(evmod, bmod)
conset <<- NULL
conset <<- do.call("contrast_set", lapply(levels(factor(des$Video)), function(v) {
f1 <- as.formula(paste("~ Video == ", paste0('"', v, '"')))
f2 <- as.formula(paste("~ Video != ", paste0('"', v, '"')))
pair_contrast(f1, f2, name=paste0(v, "_vsall"))
}))
res2 <- fmrireg:::fmri_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec"), block= ~ run,
dataset=dset,
strategy="chunkwise", nchunks=22)
res3 <- fmrireg:::fmri_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec"), block= ~ run, dataset=dset,
strategy="chunkwise", nchunks=1)
expect_true(!is.null(coef(res2)))
expect_true(!is.null(coef(res3)))
expect_true(!is.null(coef(res2, "contrasts")))
expect_true(!is.null(coef(res3, "contrasts")))
})
test_that("can run video fmri design with latent_dataset", {
#library(multivarious)
des <- read.table(system.file("extdata", "video_design.txt", package = "fmrireg"), header=TRUE)
events <- rep(320,7)
sframe <- sampling_frame(rep(320, length(events)), TR=1.5)
scans <- gen_fake_dataset(c(10,10,10,320), 7)
vecs <- lapply(scans, read_vec)
maskfile <- gen_mask_file(c(10,10,10))
mask <- read_vol(maskfile)
mats <- lapply(vecs, function(v) series(v, mask!=0))
mat <- do.call(rbind, mats)
pres <- multivarious::pca(mat, ncomp=488, preproc=multivarious::pass())
lvec <- neuroim2::LatentNeuroVec(pres$s, pres$v, add_dim(space(mask), nrow(mat)),
mask=mask)
ldset <- latent_dataset(lvec, 1.5, run_length=rep(320,7), des)
evmod <- event_model(Onset ~ hrf(Video, Condition, basis="spmg1"),
block = ~ run, sampling_frame=sframe, data=des)
conset <<- NULL
conset <<- do.call("contrast_set", lapply(levels(factor(des$Video)), function(v) {
f1 <- as.formula(paste("~ Video == ", paste0('"', v, '"')))
f2 <- as.formula(paste("~ Video != ", paste0('"', v, '"')))
pair_contrast(f1, f2, name=paste0(v, "_vsall"))
}))
conset2 <<- do.call("contrast_set", lapply(levels(factor(des$Video)), function(v) {
f1 <- as.formula(paste("~ rec_Video == ", paste0('"', v, '"')))
f2 <- as.formula(paste("~ rec_Video != ", paste0('"', v, '"')))
pair_contrast(f1, f2, name=paste0("rec_", v, "_vsall"))
}))
res2 <- fmrireg:::fmri_latent_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec", contrasts=conset2),
block= ~ run,
autocor="none", dataset=ldset)
se1 <- standard_error.fmri_latent_lm(res2, "contrasts", recon=TRUE)
con1 <- stats.fmri_latent_lm(res2, "contrasts", recon=TRUE)
dset <- fmri_dataset(scans, maskfile,TR=1.5, rep(320,7), base_path="/", event_table=des)
res2 <- fmrireg:::fmri_latent_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec", contrasts=conset2),
block= ~ run,
autocor="none", bootstrap=TRUE, nboot=50, dataset=ldset)
res2a <- fmrireg:::fmri_latent_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec", contrasts=conset2),
block= ~ run,
autocor="ar1", dataset=ldset)
res3 <- fmrireg:::fmri_lm(Onset ~ hrf(Video, subset=Condition=="Encod", contrasts=conset) +
hrf(Video, subset=Condition=="Recall", prefix="rec"), block= ~ run,
strategy="chunkwise", nchunks=1, dataset=dset)
se2 <- standard_error(res3, "contrasts")
con2 <- stats(res3, "contrasts")
expect_true(!is.null(se2))
expect_true(!is.null(con2))
})
# test_that("a one-run, one-contrast linear model analysis", {
# df1 <- subset(imagedes,run==1)
# df1 <- subset(df1, !is.na(onsetTime))
#
# df1$sdur <- scale(df1$duration)[,1]
#
# dmat <- matrix(rnorm(400*100), 400, 100)
# md <- matrix_dataset(dmat, TR=1.5, run_length=400, event_table=df1)
# con <- contrast_set(contrast( ~ Thorns - Massage, name="Thorns_Massage"))
# mod <- fmri_lm(onsetTime ~ hrf(imageName, subset = !is.na(onsetTime), contrasts=con), ~ run, dataset=md, durations=sdur)
#
# })
# test_that("a two-run, one contrast linear model analysis", {
# df1 <- subset(imagedes,run %in% c(1,2))
# df1 <- subset(df1, !is.na(onsetTime))
#
# df1$sdur <- scale(df1$duration)[,1]
#
# dmat <- matrix(rnorm(800*100), 800, 100)
# md <- matrix_dataset(dmat, TR=1.5, run_length=c(400,400), event_table=df1)
# con <- contrast_set(contrast( ~ Thorns - Massage, name="Thorns_Massage"))
# mod <- fmri_lm(onsetTime ~ hrf(imageName, contrasts=con), ~ run, dataset=md, durations=sdur)
#
#
# })
# test_that("can load and run a simple config file", {
# config <- read_fmri_config("test_data/images_study/config.R")
# dset <- fmri_dataset(config$scans, config$mask, config$TR,
# config$run_length,
# config$event_table,
# config$aux_table,
# base_path=config$base_path)
#
# frame <- sampling_frame(dset$run_length, config$TR)
# mod <- fmri_model(config$event_model, config$baseline_model, config$design, dset$aux_table,
# basis=HRF_SPMG1, dset$runids, dset$run_length, config$TR, drop_empty=TRUE)
#
#
# mod <- fmri_glm(config$event_model,
# dataset=dset, durations=0)
# })
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.