Introduction

Brief summary

Sparse canonical correlation analysis for neuroimaging (SCCAN) is a general purpose tool for “two-sided” multiple regression. This allows one to symmetrically compare one matrix of data to another and find linear relationships between them in a low-dimensional space. SCCAN derives from classic canonical correlation analysis and also relates to singular value decomposition. To handle data with \(p>>n\), SCCAN uses high-dimensional regularization methods common in \(\ell_1\) regression and spatial regularization to help ensure the biological plausibility of statistical maps in medical imaging. This problem is a difficult optimization (\(np\)-hard) and, to improve solution interpetability and stability, SCCAN allows one to to use prior knowledge to constrain the solution space. This tutorial is based on dimensionality reduction ideas outlined in the Eigenanatomy (Dhillon et al. (2014)) and SCCAN (Avants et al. (2014)) papers.

Examples

Perhaps the best way to understand how to use SCCAN is by running example data. The example data below consists of measurements from cortical gray matter and measurements from a diverse cognitive battery, the Philadelphia Brief Assessment of Cognition. The code and data below are at https://github.com/stnava/sccanTutorial and depend on ANTsR commit or more recent.

Read example data

We read in some neuroimaging and cognitive data below.

data(aal,package='ANTsR')
gfnl<-list.files(path=rootdir, pattern = glob2rx("pbac*mha"),
  full.names = T,recursive = T)
ptrainimg<-as.matrix(antsImageRead(gfnl[2],2))
ptestimg<-as.matrix(antsImageRead(gfnl[1],2))
gfnl<-list.files(path=rootdir, pattern = "gmask.nii.gz",
  full.names = T,recursive = T)
mask<-antsImageRead( gfnl[1], 3 )
afnl<-list.files(path=rootdir, pattern = "aal.nii.gz",
  full.names = T,recursive = T)
aalimg<-antsImageRead( afnl[1], 3 )
f1<-list.files(path =rootdir, pattern = "pbac_train_cog.csv",
  recursive=TRUE, full.names = TRUE, include.dirs=TRUE )
f2<-list.files(path = rootdir, pattern = "pbac_test_cog.csv",
  recursive=TRUE, full.names = TRUE )
ptraincog<-read.csv(f1)
ptestcog<-read.csv(f2)

We already divided the dataset into two different groups - one for testing and one for training.

Sparse regression

Use SCCAN to find brain regions relating to age. In this case, sparse CCA acts like a sparse regression. We impose a “cluster threshold” regularization to prevent isolated voxels from appearing in the solution. We will also compare the results in training with that in testing as a function of spareseness. This type of approach can be useful in parameter selection i.e. in choosing the optimization criterion based on training data.

agemat<-matrix( ptraincog$age, ncol=1)
paramsearch<-c(1:10)/(-100.0)
paramsearchcorrs<-rep(0,length(paramsearch))
paramsearchpreds<-rep(0,length(paramsearch))
ct<-1
for ( sp in paramsearch ) {
  ageresult<-sparseDecom2( inmatrix=list(ptrainimg,agemat), its=8, mycoption=1,
    sparseness=c(sp,0.9), inmask=c(mask,NA),nvecs=2, cthresh=c(50,0))
  # convert output images to matrix so we can validate in test data
  ccamat<-imageListToMatrix( ageresult$eig1, mask )
  agepred<-ptrainimg %*% t(ccamat)
  paramsearchcorrs[ct]<-cor( agepred[,1],  ptraincog$age )
  agepred<-ptestimg %*% t(ccamat)
  paramsearchpreds[ct]<-cor( agepred[,1],  ptestcog$age )
  ct<-ct+1
  }
mydf<-data.frame( sparseness=paramsearch, trainCorrs=paramsearchcorrs,
  testCorrs=paramsearchpreds )
mdl1<-lm( trainCorrs ~ stats::poly(sparseness,4), data=mydf )
mdl2<-lm( testCorrs ~ stats::poly(sparseness,4) , data=mydf )
visreg(mdl1)

plot of chunk sparreg

visreg(mdl2)

plot of chunk sparreg

SCCAN with prior initialization

Use SCCAN to find brain regions relating to a battery of tests that measure language-related cognitive function. We initialize SCCAN with left hemisphere regions. In this case, the initialization controls the sparseness parameters for each eigenvector.

langmat<-cbind(  ptraincog$speech_adj, ptraincog$writing_adj,
                 ptraincog$semantic_adj, ptraincog$reading_adj,
                 ptraincog$naming_adj )
colnames(langmat)<-c("speech","writing","semantic","reading","naming")
langmat2<-cbind( ptestcog$speech_adj, ptestcog$writing_adj,
                 ptestcog$semantic_adj, ptestcog$reading_adj,
                 ptestcog$naming_adj )
colnames(langmat2)<-colnames(langmat)
labels<-c(13,81,39,79)
print(aal$label_name[labels])
## [1] Frontal_Inf_Tri_L Temporal_Sup_L    ParaHippocampal_L Heschl_L         
## 116 Levels: Amygdala_L Amygdala_R Angular_L Angular_R ... Vermis_9
initmat<-matrix( rep(0,sum(mask==1)*length(labels)), nrow=length(labels) )
# fill the matrix with the aal region locations
for ( i in 1:length(labels) ) {
  vec<-( aalimg[ mask == 1 ] == labels[i] )
  initmat[i,]<-as.numeric( vec )
}
ccainit<-initializeEigenanatomy( initmat, mask )
pwsearch<-c(50,25,10)
langfn<-rep("",length(pwsearch))
langfn2<-rep("",length(pwsearch))
ct<-1
for ( pw in pwsearch ) {
langresult<-sparseDecom2( inmatrix=list(ptrainimg,langmat), its=15, mycoption=1,
  sparseness=c(sp,-0.5), inmask=c(mask,NA),nvecs=length(labels), cthresh=c(100,0),
  initializationList=ccainit$initlist, priorWeight=pw/100, smooth=0.0, ell1=-10 )
ccamat<-imageListToMatrix( langresult$eig1, mask )
langpred<-ptrainimg %*% t(ccamat)
colnames(langpred)<-paste("GM",c(1:ncol(langpred)),sep='')
cogpred<-langmat %*% data.matrix( langresult$eig2 )
bestpred<-which.max(abs(diag(cor(langpred,cogpred))))
mydf<-data.frame( cogpred, langpred )
myform<-as.formula( paste("Variate00",bestpred-1,"~GM1+GM2+GM3+GM4",sep='') )
mdltrain<-lm( myform, data=mydf )
langpred<-ptestimg %*% t(ccamat)
colnames(langpred)<-paste("GM",c(1:ncol(langpred)),sep='')
cogpred<-langmat2 %*% data.matrix( langresult$eig2 )
mydf<-data.frame( cogpred, langpred )
print(cor.test( mydf[,bestpred] ,predict(mdltrain,newdata=mydf)))
for ( i in 1:length(labels) )
  print( paste( "Dice: ",aal$label_name[labels[i]],
         sum( abs(ccamat[i,]) > 0 & initmat[i,] > 0 ) /
         sum( abs(ccamat[i,]) > 0 | initmat[i,] > 0 ) ) )
for ( x in langresult$eig1 ) {
  x[ mask == 1 ]<-abs( x[ mask == 1 ] )
  x[ mask == 1 ]<-x[ mask == 1 ]/max( x[ mask == 1 ] )
}
mycolors<-c("red","green","blue","yellow")
langfn[ct]<-paste(rootdir,'/figures/langSCCANRegression',pw,'.jpg',sep='')
langfn2[ct]<-paste(rootdir,'/figures/langSCCANRegression',pw,'.png',sep='')
plotANTsImage( mask, functional=(langresult$eig1), threshold='0.25x1',
  slices="12x50x1",color=mycolors,outname=langfn[ct] )
# cnt<-getCentroids( ntwkimage, clustparam = 100 )
brain<-renderSurfaceFunction( surfimg =list( aalimg ) , alphasurf=0.1 ,
  funcimg=langresult$eig1, smoothsval=1.5, smoothfval=0, mycol=mycolors )
id<-par3d("userMatrix")
rid<-rotate3d( id , -pi/2, 1, 0, 0 )
rid2<-rotate3d( id , pi/2, 0, 0, 1 )
rid3<-rotate3d( id , -pi/2, 0, 0, 1 )
par3d(userMatrix = id )
dd<-make3ViewPNG(  rid, id, rid2, paste(rootdir,'/figures/langSCCANRegression',pw,sep='') )
par3d(userMatrix = id )
ct<-ct+1
}
## 
##  Pearson's product-moment correlation
## 
## data:  mydf[, bestpred] and predict(mdltrain, newdata = mydf)
## t = 3.404, df = 81, p-value = 0.001036
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1494 0.5291
## sample estimates:
##    cor 
## 0.3537 
## 
## [1] "Dice:  Frontal_Inf_Tri_L 0.92090395480226"
## [1] "Dice:  Temporal_Sup_L 0.901474530831099"
## [1] "Dice:  ParaHippocampal_L 1"
## [1] "Dice:  Heschl_L 0"
## 
##  Pearson's product-moment correlation
## 
## data:  mydf[, bestpred] and predict(mdltrain, newdata = mydf)
## t = 6.223, df = 81, p-value = 2.032e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4024 0.6987
## sample estimates:
##    cor 
## 0.5687 
## 
## [1] "Dice:  Frontal_Inf_Tri_L 0.139437689969605"
## [1] "Dice:  Temporal_Sup_L 0.97645327446652"
## [1] "Dice:  ParaHippocampal_L 0.96010296010296"
## [1] "Dice:  Heschl_L 0"
## 
##  Pearson's product-moment correlation
## 
## data:  mydf[, bestpred] and predict(mdltrain, newdata = mydf)
## t = 6.748, df = 81, p-value = 2.059e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4413 0.7221
## sample estimates:
##    cor 
## 0.5999 
## 
## [1] "Dice:  Frontal_Inf_Tri_L 0.11939736346516"
## [1] "Dice:  Temporal_Sup_L 0.0382615156017831"
## [1] "Dice:  ParaHippocampal_L 0.51002004008016"
## [1] "Dice:  Heschl_L 0"

Strong use of prior

Strong prior