Yes thank you! The experience was very good, its very useful! I think it just need to learn how check, verify "goodness of fit".
Then, you can directly access the Markov chain Monte Carlo simulation results. Utilize the functions in MBA.R in your AFNI binary directory to construct uncertainty intervals with your preferred quantile. For instance, information about region pairs can be found under the heading:
Ah, very cool! so by running the following one can access all the subthreshold details?
rp_summary <- function(ps, ns, nR) {
mm <- apply(ps, c(2,3), mean)
for(ii in 1:nR) for(jj in 1:nR) ps[,ii,jj] <- sqrt(2)*(ps[,ii,jj] - mm[ii,jj]) + mm[ii,jj]
RP <- array(NA, dim=c(nR, nR, 8))
RP[,,1] <- apply(ps, c(2,3), mean)
RP[,,2] <- apply(ps, c(2,3), sd)
RP[,,3] <- apply(ps, c(2,3), cnt, ns)
RP[,,4:8] <- aperm(apply(ps, c(2,3), quantile, probs=c(0.025, 0.05, 0.5, 0.95, 0.975)), dim=c(2,3,1))
dimnames(RP)[[1]] <- dimnames(ps)[[2]]
dimnames(RP)[[2]] <- dimnames(ps)[[3]]
dimnames(RP)[[3]] <- c('mean', 'SD', 'P+', '2.5%', '5%', '50%', '95%', '97.5%')
return(RP)
}
# extract region-pair posterior samples for an effect 'tm'
region_pair <- function(pe, ge, tm, nR) {
ps0 <- array(apply(ge[['mmROI1ROI2']][,,tm], 2, "+", ge[['mmROI1ROI2']][,,tm]), c(ns, nR, nR)) # obtain (ROIi, ROIj) matrix for region-specific main effects
ps <- apply(ps0, c(2,3), '+', pe[,tm]) # add the common effect to the matrix
dimnames(ps) <- list(1:ns, rL, rL)
tmp <- ps
sel1 <- match(dimnames(ge$`ROI1:ROI2`)[[2]], outer(dimnames(ps)[[2]],dimnames(ps)[[3]], function(x,y) paste(x,y,sep="_")))
sel2 <- match(dimnames(ge$`ROI1:ROI2`)[[2]], outer(dimnames(ps)[[2]],dimnames(ps)[[3]], function(x,y) paste(y,x,sep="_")))
ad <- function(tt,ge,s1,s2) {tt[s1] <- tt[s1] + ge; tt[s2] <- tt[s2] + ge; return(tt)}
for(ii in 1:ns) tmp[ii,,] <- ad(tmp[ii,,], ge$`ROI1:ROI2`[ii,,tm], sel1, sel2) # add the interaction
ps <- tmp
return(ps)
}
ns <- lop$iterations*lop$chains/2
#nR <- nlevels(lop$dataTable$ROI1)
pe <- fixef(fm, summary = FALSE) # Population-Level Estimates
ge <- ranef(fm, summary = FALSE) # Extract Group-Level (or random-effect) Estimates
nR <- length(union(levels(lop$dataTable$ROI1), levels(lop$dataTable$ROI2)))
if(!is.na(lop$ROIlist)) {
lop$ROI$order <- as.character(lop$ROI$order)
if(!setequal(lop$ROI$order, union(levels(lop$dataTable$ROI1), levels(lop$dataTable$ROI2))))
errex.AFNI(paste("The ROIs listed under option -ROIlist do not fully match those specified in the data table!"))
# change the factor levels to characters
}
rp_summary(region_pair(pe, ge, lop$EOIq[ii], nR), ns, nR)
You can find an R plot function here. Simply provide the posterior samples as input in the format of the data.frame dat .
Hm... so that would the output of the region_summary function right?
glimpse(region_pair(pe, ge, lop$EOIq[ii], nR))
num [1:2000, 1:8, 1:8] 0.05109 0.00483 -0.0305 -0.05837 -0.00657 ...
- attr(*, "dimnames")=List of 3
..$ : chr [1:2000] "1" "2" "3" "4" ...
..$ : chr [1:8] "CA1" "CA2" "CA3" "CA4" ...
..$ : chr [1:8] "CA1" "CA2" "CA3" "CA4" ...
but it will have to be altered so as to only include unique pairwise combination of N rois right? I basically need to select the "upper triangle" of the NxNx2000 matrix?
Could you illustrate how one ought to do the posterior predictive checks?
