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

Strong prior 3D

Medium use of prior

Medium prior

Medium prior 3D

Weak use of prior

Weak prior

Weak prior 3D

Identifying the anatomical network

The best results are initialized by the prior but, in the end, drift away from that initialization. Where in the brain do the solution vectors end up? We write a quick function to answer this question.

reportAnatomy<-function( eigIn, maskIn, wt=0.3 )
  {
  data('aal',package='ANTsR')
  ccaanat<-list()
  for ( img in eigIn ) {
    nzind<-abs(img[ maskIn == 1 ]) > 0
    aalvals<-aalimg[ maskIn == 1 ][ nzind ]
    ccaanat<-lappend( ccaanat, aalvals )
  }
  ccaanat<-unlist( ccaanat )
  anatcount<-hist(ccaanat,breaks=0:100, plot = F)$count
  anatcount[ anatcount < wt*max(anatcount) ]<-0
  anatcount<-which( anatcount > 0 )
  return( toString(aal$label_name[anatcount] ) )
  }
ccaaal<-reportAnatomy( langresult$eig1 , mask )

The SCCAN predictors include: Frontal_Inf_Tri_L, Frontal_Inf_Tri_R, ParaHippocampal_L, ParaHippocampal_R, Fusiform_R, Precuneus_R, Temporal_Mid_R, Temporal_Inf_R.

How good were our original hypothetical regions as predictors?

Associating classes to SCCAN predictors

Recalling: CCA maximizes \(PearsonCorrelation( XW^T, ZY^T )\), where \(X\) and \(Z\) are data matrices, we can study matrix \(Y\) (or \(W\)) which contrasts or combines columns of the associated data matrix. In this example, \(Y\) operates on the cognition/design matrix.

rownames(langresult$eig2)<-colnames(langmat)
temp<-(langresult$eig2)
temp[ abs(langresult$eig2) < 0.03 ]<-0

pheatmap(temp)

Sparse regression with nuisance variables

Often, we want to control for the presence of nuisance variables. As usual, there are several options: (1) control after you do dimensionality reduction; (2) orthogonalize the predictors. (3) Use alternative SCCAN formulations (e.g. set mycoption to 0 or 2). Let’s try the first 2 choices as they are more traditional.

# 1. control for age and mmse after the dimensionality reduction
langresult<-sparseDecom2( inmatrix=list(ptrainimg,langmat), its=15, mycoption=1,
  sparseness=c(sp,-0.9), inmask=c(mask,NA),nvecs=length(labels), cthresh=c(50,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, mmse=ptraincog$mmse,age=ptraincog$age)
myform<-as.formula( paste("Variate00",bestpred-1,
  "~GM1+GM2+GM3+GM4+mmse+age",sep='') )
mdltrain<-lm( myform, data=mydf )
print(summary(mdltrain))

## 
## Call:
## lm(formula = myform, data = mydf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1455 -0.0569 -0.0101  0.0426  0.2427 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.77952    0.24416    3.19    0.002 ** 
## GM1           3.35687    1.86484    1.80    0.076 .  
## GM2         -14.66344    2.27809   -6.44  7.7e-09 ***
## GM3          -1.96884    1.73766   -1.13    0.260    
## GM4           6.12872    2.97341    2.06    0.042 *  
## mmse         -0.00784    0.00149   -5.27  1.1e-06 ***
## age           0.00124    0.00115    1.08    0.284    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0848 on 82 degrees of freedom
## Multiple R-squared:  0.641,  Adjusted R-squared:  0.615 
## F-statistic: 24.4 on 6 and 82 DF,  p-value: <2e-16

langpred2<-ptestimg %*% t(ccamat)
colnames(langpred2)<-paste("GM",c(1:ncol(langpred)),sep='')
cogpred2<-langmat2 %*% data.matrix( langresult$eig2 )
mydf<-data.frame( cogpred2, langpred2,mmse=ptestcog$mmse,age=ptestcog$age )
print(cor.test( mydf[,bestpred] ,predict(mdltrain,newdata=mydf)))

## 
##  Pearson's product-moment correlation
## 
## data:  mydf[, bestpred] and predict(mdltrain, newdata = mydf)
## t = 7.932, df = 81, p-value = 1.023e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5196 0.7674
## sample estimates:
##    cor 
## 0.6612

Now, the second option.

# 2. orthogonalize the matrices against mmse and education
rlangmat<-residuals(lm(langmat~ptraincog$mmse+ptraincog$age))
rptrainimg<-residuals(lm(ptrainimg~ptraincog$mmse+ptraincog$age))
langresult<-sparseDecom2( inmatrix=list(rptrainimg,rlangmat), its=15, mycoption=1,
  sparseness=c(sp,-0.9), inmask=c(mask,NA),nvecs=length(labels), cthresh=c(50,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, mmse=ptraincog$mmse,age=ptraincog$age)
myform<-as.formula( paste("Variate00",bestpred-1,
  "~GM1+GM2+GM3+GM4+mmse+age",sep='') )
mdltrain<-lm( myform, data=mydf )
print(summary(mdltrain))

## 
## Call:
## lm(formula = myform, data = mydf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1793 -0.0569  0.0033  0.0547  0.2156 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.06350    0.17514   -6.07  3.8e-08 ***
## GM1          4.38088    1.44269    3.04   0.0032 ** 
## GM2         22.92672    2.07253   11.06  < 2e-16 ***
## GM3         -4.30491    1.61736   -2.66   0.0094 ** 
## GM4         -0.15167    2.10902   -0.07   0.9428    
## mmse         0.01071    0.00125    8.58  4.9e-13 ***
## age          0.00218    0.00114    1.91   0.0600 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0806 on 82 degrees of freedom
## Multiple R-squared:  0.779,  Adjusted R-squared:  0.763 
## F-statistic: 48.2 on 6 and 82 DF,  p-value: <2e-16

langpred2<-ptestimg %*% t(ccamat)
colnames(langpred2)<-paste("GM",c(1:ncol(langpred)),sep='')
cogpred2<-langmat2 %*% data.matrix( langresult$eig2 )
mydf<-data.frame( cogpred2, langpred2,mmse=ptestcog$mmse,age=ptestcog$age )
print(cor.test( mydf[,bestpred] ,predict(mdltrain,newdata=mydf)))

## 
##  Pearson's product-moment correlation
## 
## data:  mydf[, bestpred] and predict(mdltrain, newdata = mydf)
## t = 6.923, df = 81, p-value = 9.505e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4537 0.7295
## sample estimates:
##    cor 
## 0.6097

Predicting the full cognitive battery from the neuroimaging data

Try to predict all the demographic variability from the imaging data. We use mycoption 0 to try to reduce correlation in low-dimensional space. This enforces a new SCCAN constraint (not previously reported).

nv<-11
nfn<-rep("",nv)
cognames<-rep("",nv)
cogmat<-data.matrix(ptraincog)
rcogmat<-residuals( lm( data.matrix(ptraincog) ~ ptraincog$mmse + ptraincog$age ) )
rptrainimg<-residuals( lm( ptrainimg ~ ptraincog$mmse ) )
batt<-sparseDecom2( inmatrix=list(rptrainimg,rcogmat), its=15,
  sparseness=c(0.015, -0.05), inmask=c(mask,NA), nvecs=nv, cthresh=c(50,0),
  smooth=0.25, ell1=-10, mycoption=1 )
ccamat<-imageListToMatrix( batt$eig1, mask )
gvars<-paste("GM",c(1:nrow(ccamat)),sep='',collapse='+')

Now let’s use our previously developed reporting capabilities.

render<-TRUE
for ( bestpred in 1:nrow(ccamat)) {
  battpred<-ptrainimg %*% t(ccamat)
  colnames(battpred)<-paste("GM",c(1:ncol(battpred)),sep='')
  cogpred<-( rcogmat %*% data.matrix( batt$eig2 ) )[,bestpred]
  mydf<-data.frame( cogpred, battpred )
  myform<-as.formula( paste("cogpred~",gvars,sep='') )
  mdltrain<-lm( myform, data=mydf )
  mdlinterp<-bigLMStats( mdltrain )
  battpred<-ptestimg %*% t(ccamat)
  colnames(battpred)<-paste("GM",c(1:ncol(battpred)),sep='')
  cogpred<-(data.matrix(ptestcog) %*% data.matrix( batt$eig2 ))[,bestpred]
  mydf<-data.frame( cogpred, battpred )
  cat(paste("Eig",bestpred,"is related to:\n"))
  mycog<-colnames(ptraincog)[ abs(batt$eig2[,bestpred]) > 0 ]
  cat( mycog )
  cat("\nwith weights\n")
  cat( abs(batt$eig2[,bestpred])[ abs(batt$eig2[,bestpred]) > 0 ])
  cat(paste("\nwith predictive correlation:",
    cor( cogpred,predict(mdltrain,newdata=mydf))))
  cat("\nAnatomy:")
  for ( x in which.min(p.adjust(mdlinterp$beta.pval)) )  {
    myanat<-reportAnatomy( list( batt$eig1[[x]]) , mask , 0.5 )
    cat(myanat)
    if ( render ) {
    vizimg<-antsImageClone( batt$eig1[[x]] )
    ImageMath(3,vizimg,'abs',vizimg)
    brain<-renderSurfaceFunction( surfimg =list( aalimg ) , alphasurf=0.1 ,
      funcimg=list(vizimg), smoothsval = 1.5 )
    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 )
    ofn<-paste(rootdir,'/figures/battery',bestpred,sep='')
    nfn[ bestpred ]<-paste(ofn,'.png',sep='')
    cognames[ bestpred ]<-paste(mycog,collapse='+')
    dd<-make3ViewPNG(  rid, id, rid2, ofn )
    par3d(userMatrix = id )
    }
    cat("\n")
  }
  cat("\n")
}

## Eig 1 is related to:
## recog_adj
## with weights
## 0.1066
## with predictive correlation: 0.0251593498713502
## Anatomy:ParaHippocampal_L, ParaHippocampal_R, Temporal_Inf_L, Temporal_Inf_R
## 
## Eig 2 is related to:
## socialcomportment
## with weights
## 0.1066
## with predictive correlation: 0.14257952784209
## Anatomy:Frontal_Inf_Oper_R, Rectus_L, Insula_L, Insula_R, Caudate_R
## 
## Eig 3 is related to:
## rey_copy_adj
## with weights
## 0.1066
## with predictive correlation: 0.553698776973147
## Anatomy:Occipital_Mid_L
## 
## Eig 4 is related to:
## fluency_adj
## with weights
## 0.1066
## with predictive correlation: 0.0476808144149163
## Anatomy:Cuneus_R, Precuneus_R
## 
## Eig 5 is related to:
## disinhibition
## with weights
## 0.1066
## with predictive correlation: 0.554939219976453
## Anatomy:Frontal_Med_Orb_L, Rectus_L, Cingulum_Ant_L
## 
## Eig 6 is related to:
## delay_free_adj
## with weights
## 0.1066
## with predictive correlation: 0.134661661896197
## Anatomy:Cingulum_Mid_R, Fusiform_R
## 
## Eig 7 is related to:
## apathy
## with weights
## 0.1066
## with predictive correlation: 0.492559320532837
## Anatomy:Temporal_Sup_R
## 
## Eig 8 is related to:
## rey_recall_adj
## with weights
## 0.1066
## with predictive correlation: 0.326893909184624
## Anatomy:Frontal_Inf_Tri_R, Hippocampus_L, Fusiform_R, Temporal_Inf_L
## 
## Eig 9 is related to:
## JOLO_adj
## with weights
## 0.1066
## with predictive correlation: 0.288852475356899
## Anatomy:Occipital_Mid_L
## 
## Eig 10 is related to:
## rey_recall_adj
## with weights
## 0.1066
## with predictive correlation: 0.326893909184624
## Anatomy:Frontal_Inf_Tri_R, Hippocampus_L, Fusiform_R, Temporal_Inf_L
## 
## Eig 11 is related to:
## empathy
## with weights
## 0.1066
## with predictive correlation: 0.272080219991265
## Anatomy:Temporal_Sup_R

Anatomy related with recog_adj

Select results

Anatomy related with socialcomportment

Select results

Anatomy related with rey_copy_adj

Select results

Anatomy related with fluency_adj

Anatomy related with disinhibition

Select results

Anatomy related with delay_free_adj

Select results

Anatomy related with apathy

Select results

Anatomy related with rey_recall_adj

Select results

Anatomy related with JOLO_adj

Select results

Anatomy related with rey_recall_adj

Select results

Anatomy related with empathy

Select results

Can the neuroimaging data predict the full cognitive battery?

# use cca to transform cortical signal to the cognitive battery
batt2<-sparseDecom2( inmatrix=list(rptrainimg,rcogmat), its=15,
  sparseness=c(0.02, -0.9), inmask=c(mask,NA), nvecs=nv, cthresh=c(100,0),
  smooth=0.0, ell1=-10, mycoption=1 )
ccamat<-imageListToMatrix( batt2$eig1, mask )
predictedBattery<-data.frame( vox=ptrainimg %*% t(ccamat) %*% t(batt2$eig2) )
print(diag(cor(cogmat,predictedBattery)))

##  [1] 0.04707 0.18585 0.12448 0.08079 0.05763 0.01225 0.51996 0.44400
##  [9] 0.46830 0.38188 0.39129 0.48106 0.37810 0.52096 0.50097 0.10328
## [17] 0.39308 0.42130 0.43802 0.19642 0.32807

predictedBattery<-data.frame( vox=ptestimg %*% t(ccamat) %*% t(batt2$eig2) )
print(diag(cor(data.matrix(ptestcog),predictedBattery)))

##  [1] -0.22095  0.07252 -0.04368 -0.10955 -0.01810 -0.04819  0.39080
##  [8]  0.08382  0.14221  0.08864  0.36601  0.50798  0.31882  0.41619
## [15]  0.34914  0.27937  0.27407  0.44326  0.10207  0.12704  0.37066

qv<-rep(NA,ncol(ptestcog) )
for ( i in 1:ncol(ptestcog) ) {
 qv[i]<-cor.test(data.matrix(ptestcog)[,i],predictedBattery[,i])$p.value
 ttl<-paste(  colnames(ptestcog)[i],
      cor(data.matrix(ptestcog)[,i],predictedBattery[,i]) )
 mdl<-data.frame( realCog=data.matrix(ptestcog)[,i],
                  predCog=predictedBattery[,i] )
 mylm<-lm( predCog ~ realCog , data=mdl )
 visreg( mylm , main=ttl)
 Sys.sleep(1)
}

plot of chunk mvarpred

qv[ is.na(qv) ]<-1
qv<-p.adjust(qv,method='BH')

The following univariate columns may be predicted using SCCAN multivariate mapping: naming_adj, semantic_adj, delay_free_adj, recog_adj, rey_recall_adj, JOLO_adj, rey_copy_adj, apathy, disinhibition, empathy.

Avants, Brian B., David J. Libon, Katya Rascovsky, Ashley Boller, Corey T. McMillan, Lauren Massimo, H Branch Coslett, Anjan Chatterjee, Rachel G. Gross, and Murray Grossman. 2014. “Sparse Canonical Correlation Analysis Relates Network-Level Atrophy to Multivariate Cognitive Measures in a Neurodegenerative Population.” Neuroimage 84 (Jan). Department of Radiology, University of Pennsylvania School of Medicine, Philadelphia, PA, USA.: 698–711. doi:10.1016/j.neuroimage.2013.09.048.

Dhillon, Paramveer S., David A. Wolk, Sandhitsu R. Das, Lyle H. Ungar, James C. Gee, and Brian B. Avants. 2014. “Subject-Specific Functional Parcellation via Prior Based Eigenanatomy.” Neuroimage, May. Penn Image Computing; Science Laboratory (PICSL), Department of Radiology, University of Pennsylvania, Philadelphia, PA, USA. doi:10.1016/j.neuroimage.2014.05.026.

SCCAN tutorial and tests

Introduction

Brief summary

Examples

Read example data

Sparse regression

SCCAN with prior initialization

Identifying the anatomical network

Associating classes to SCCAN predictors

Sparse regression with nuisance variables

Predicting the full cognitive battery from the neuroimaging data