inst/doc/d_Common_Statistical_Approach.R

In mand: Multivariate Analysis for Neuroimaging Data

## ----echo=FALSE, error=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
options(warn=-1)
knitr::opts_chunk$set(eval = FALSE)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  n = 20; x = 1:n
#  set.seed(1); y = abs(rnorm(n)); names(y) = x

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  cuty = function(y1, cdt){
#  rbind(ifelse(y1-cdt<0, y1, cdt),
#  ifelse(y1-cdt<0, 0, y1-cdt))
#  }

## ----fig.width = 7, fig.height = 4------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  par(mfrow=c(1, 2), mar=c(4,3,4,3))
#  for(cdt in c(0.5, 1.5)){
#  y2 = cuty(y, cdt)
#  barplot(y2, col=c(8,2), main=paste("CDT =", cdt))
#  abline(h = cdt, lty=2, col=2)
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  cdt = 0.75

## ----fig.width = 7, fig.height = 4------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  par(mfrow=c(1,2), mar=c(4,3,4,3))
#  for(f1 in c(2,8)){
#  sy = gsmooth(x, y, f1); y2 = cuty(sy, cdt)
#  barplot(y2, col=c(8,2), main=paste("FWHM =", f1), ylim=c(0,2))
#  abline(h = cdt, lty=2, col=2)
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  exp(-0.53*(96/325*3.3*(3.3^2-1)^(2/3)))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  s=96; h=3.3; r=325;
#  Es=(4*log(2))^(3/2)*h*(h^2-1)/(2*pi)^(3/2);
#  d=(gamma(5/2)*Es)^(2/3);
#  exp(-d*(s/r)^(2/3))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (325*(-log(0.05)/0.53)^(3/2))/(3.3*(3.3^2-1))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  ss = c(1, 10, 50)
#  h = 4
#  b = (4*log(2)/(2*pi))*(gamma(5/2))^(2/3)
#  fs = 2:10

## ----fig.width = 5, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  par(mar=c(2,2,2,6))
#  for(sidx in 1:length(fs)){
#  p = exp(-b*(ss[sidx]/fs^3*h*((h)^2-1))^(2/3))
#  if(sidx == 1){
#  plot(fs, p, lty = sidx, type="l", ylim=c(0,1),
#  xlab="FWHM", ylab="p-value", main="")
#  }else{
#  points(fs, p, lty = sidx, type="l")
#  }}
#  abline(h=0.05, lty=5, lwd=2)
#  par(xpd=TRUE)
#  legend(par()$usr[2], par()$usr[4],
#  legend=paste("s =", ss[order(ss)]), lty=order(ss))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  TFCE0 = function(y1, cdt, E=0.5, H=2){
#  cy = 1*(y1 >= cdt)
#  
#  clsta0 = c(1, as.numeric(names(which(diff(cy)!=0))))
#  clend0 = c(clsta0[-1]-1, length(cy))
#  clsta = clsta0[cy[clsta0]==1]
#  clend = clend0[cy[clsta0]==1]
#  
#  clustidxs = lapply(1:length(clsta),
#  function(i) clsta[i]:clend[i])
#  clustsize = unlist(lapply(clustidxs, length))
#  
#  x = 1:length(y1)
#  clust = unlist(lapply(x, function(x1){
#  a = which(unlist(lapply(clustidxs,
#  function(cidx) x1 %in% cidx)))
#  ifelse(length(a)>0, a, NA)
#  }))
#  
#  a = clustsize[clust]^E * cdt^H
#  ifelse(is.na(a),0, a)
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  TFCE1 = function(f1){
#  for(cdt1 in cdts2){
#  sy = gsmooth(x, y, f1)
#  y2 = cuty(sy, cdt1)
#  barplot(y2, col=c(8,2), main=paste("CDT =", cdt1))
#  abline(h = cdt1, lty=2, col=2)
#  tfce = TFCE0(sy,cdt1)
#  barplot(tfce, main=paste("TFCE CDT =", cdt1), ylim=c(0,10))
#  }
#  }

## ----fig.width = 6, fig.height = 5------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  cdts2 = c(0.75, 0.5, 0.25)
#  par(mfrow=c(length(cdts2), 2), mar=c(3,4,2,4))
#  TFCE1(f1=2)

## ----fig.width = 7, fig.height = 6------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  f1s = c(0, 2, 4, 8)
#  par(mfrow=c(length(f1s),2), mar=c(3,4,2,4))
#  for(f1 in f1s){
#  if(f1 == 0){barplot(y, main="Original")}else{
#  sy = gsmooth(x, y, f1)
#  barplot(sy, main=paste("FWHM =", f1), ylim=c(0,2))
#  }
#  if(f1 > 0){ sy = gsmooth(x, y, f1) }else{ sy = y }
#  cdts = c(0, sort(unique(sy)))
#  tfce = colSums(do.call(rbind,lapply(cdts, function(cdt1){
#  TFCE0(sy,cdt1)
#  })))
#  barplot(tfce, main=paste("TFCE FWHM =", f1))
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  n = 5; x0 = rep(c(0,1), each = n)
#  set.seed(1); y = ifelse(x0==0, rnorm(n, 0), rnorm(n, 1))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  fitfunc = function(idx){
#  x = x0[idx]
#  fit = summary(lm(y~x))
#  coef(fit)[2,1:3]
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (fit0 = fitfunc(1:(2*n)))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1); idx = sample(2*n); round(fitfunc(idx), 3)
#  
#  set.seed(2); idx = sample(2*n); round(fitfunc(idx), 3)
#  
#  set.seed(1000); idx = sample(2*n); round(fitfunc(idx), 3)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  tststs = sapply(1:1000, function(s1){
#  set.seed(s1); idx = sample(2*n); fitfunc(idx)[3]
#  })

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (q1 = quantile(tststs, prob=0.975))
#  ifelse(abs(fit0[3])>q1, "reject", "not reject")

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (q2 = qt(0.975, 2*n-2) )
#  ifelse(abs(fit0[3])>q2, "reject", "not reject")

## ----fig.width = 4, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  hist(tststs, xlab="T stat", main="", freq=FALSE,
#  ylim=c(0, 0.4))
#  abline(v=fit0[3], col=2, lty=2)
#  abline(v=q1, col=4, lty=3); abline(v=q2, col=3, lty=4)
#  f1 = function(x)dt(x, 2*n-2)
#  curve(f1, -5, 5, add=TRUE, col=3)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (pp = 2*(mean(tststs > fit0[3])))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  (tp = 2*(1 - pt(fit0[3], 2*n-2)))

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  nreps = round(c(seq(10, 100, length=10),
#  seq(200, length(tststs)-100, length=10)))
#  q1s = sapply(nreps, function(x) {
#  sapply(1:1000, function(y) {
#  set.seed(y); quantile(tststs[sample(1:length(tststs), x)],
#  0.975)})})
#  colnames(q1s)=nreps

## ----fig.width = 4, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  boxplot(q1s, main="", xlab="Number of permutations",
#  ylab="Two-sided 95% point", type="b")
#  abline(h=q2, col=3, lty=2)

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  clustsize = function(y1, cdt){
#  cy = 1*(y1 >= cdt)
#  
#  if(all(cy==0)){
#  clustsize = 0
#  }else{
#  clsta0 = c(1, as.numeric(names(which(diff(cy)!=0))))
#  clend0 = c(clsta0[-1]-1, length(cy))
#  clsta = clsta0[cy[clsta0]==1]
#  clend = clend0[cy[clsta0]==1]
#  
#  clustidxs = lapply(1:length(clsta),
#  function(i) clsta[i]:clend[i])
#  clustsize = unlist(lapply(clustidxs, length))
#  }
#  clustsize
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  gsmooth = function(x, y, FWHM){
#  sigma = FWHM / 2*sqrt(2*log(2))
#  sy = sapply(x, function(x1)
#  weighted.mean(y,
#  dnorm(x, x1, sigma)/sum(dnorm(x, x1,sigma))) )
#  names(sy) = x
#  sy
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  cuty = function(y1, cdt){
#  rbind(ifelse(y1-cdt<0, y1, cdt),
#  ifelse(y1-cdt<0, 0, y1-cdt))
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  n = 20; x = 1:n; f1 = 2; cdt1 = 0.75

## ----fig.width = 7, fig.height = 6------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  sidxs = 1:4
#  par(mfrow=c(length(sidxs),1), mar=c(3,4,2,4))
#  for(sidx in sidxs){
#  sigma1 = ifelse(sidx==1, 1, 0.9)
#  set.seed(sidx); y = abs(rnorm(n,0,sigma1)); names(y) = x
#  sy = gsmooth(x, y, f1); y2 = cuty(sy, cdt1)
#  csize = clustsize(sy, cdt1)
#  barplot(y2, col=c(8,2),
#  main=paste("Max Cluster Size =", max(csize)), ylim=c(0,1.2))
#  abline(h = cdt1, lty=2, col=2)
#  }

## ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  sidxs2 = 1:1001
#  maxcsizes = sapply(sidxs2, function(sidx){
#  sigma1 = ifelse(sidx==1, 1, 0.9)
#  set.seed(sidx); y = abs(rnorm(n,0,sigma1)); names(y) = x
#  sy = gsmooth(x, y, f1)
#  csize = clustsize(sy, cdt1)
#  max(csize)
#  })
#  (obs = maxcsizes[1])
#  (q1 = quantile(maxcsizes[-1], 0.95))

## ----fig.width = 4, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  hist(maxcsizes, main="", xlab="Max Cluster Size")
#  abline(v = q1, col=2, lty=2); abline(v = obs, col=3, lty=2)

## ----fig.width = 4, fig.height = 3.5----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  nreps = seq(10, length(maxcsizes), length=100)
#  q1s = sapply(nreps, function(x)
#  quantile(maxcsizes[2:x], 0.95))
#  par(mfrow=c(1,1))
#  plot(nreps, q1s, main="", xlab="Number of permutations",
#  ylab="95% point", type="b")