Penalized Likelihood Phenotyping: Unifying Voxelwise Analyses and Multi-Voxel Pattern Analyses in Neuroimaging  

Authors: Nagesh Adluru*, Bret M. Hanlon, Antoine Lutz, Janet E. Lainhart, Andrew L. Alexander, Richard J. Davidson

*Corresponding Author (adluru@wisc.edu)

 

Basic overview:

Neuroimage phenotyping for psychiatric and neurological disorders is performed using voxelwise analyses also known as voxel based analyses or morphometry (VBM). A typical voxelwise analysis treats measurements at each voxel (e.g. fractional anisotropy, gray matter probability) as outcome measures to study the effects of possible explanatory variables (e.g. age, group) in a linear regression setting. Furthermore, each voxel is treated independently until the stage of correction for multiple comparisons. Recently, multi-voxel pattern analyses (MVPA), such as classification, have arisen as an alternative to VBM. The main advantage of MVPA over VBM is that the former employ multivariate methods which can account for interactions among voxels in identifying significant patterns. They also provide ways for computer-aided diagnosis and prognosis at individual subject level. However, compared to VBM, the results of MVPA are often more difficult to interpret and prone to arbitrary conclusions. In this paper, first we use penalized likelihood modeling to provide a unified framework for understanding both VBM and MVPA. We then utilize statistical learning theory to provide practical methods for interpreting the results of MVPA beyond commonly used performance metrics, such as leave-one-out-cross validation accuracy and area under the receiver operating characteristic (ROC) curve. Additionally, we demonstrate that there are challenges in MVPA when trying to obtain image phenotyping information in the form of statistical parametric maps (SPMs), which are commonly obtained from VBM, and provide a bootstrap strategy as a potential solution for generating SPMs using MVPA. This technique also allows us to maximize the use of available training data. We illustrate the empirical performance of the proposed framework using two different neuroimaging studies that pose different levels of challenge for classification using MVPA.

 

 

 

Sample Results:

 

ROC optimal regionsz-score maps

 

For record-keeping sake, please email to the corresponding author for a copy of the MATLAB scripts used to produce the results published in the paper.

 

Snapshot of sample usage of the code:

MATLAB code for penalized likelihood phenotyping.

%Setting up parameters.
M.name='LIBLINEAR';
M.C=1;
M.no_of_folds=10;
M.no_of_iter=100;
M.gridfold=10;
M.gridsearch=false;
M.verbose=true;
M.grid=2.^(-4:1:6)';%Based on the Yuang et al. TR.
labels=diagnosis;
labels(labels==0)=-1;


%FA.
%--------------------------------------------------------------------------
feature_matrix=sparse([sfa siq]);
M.loss_regularization=1;%L2SVM-L2 100 iters, 10-fold.
[acc,spec,sens,conf,auroc,roc,loo_models,loo_acc,~,~,~,~,~,loo_scores,loo_labels]=get_k_fold_vals(feature_matrix,labels,M); save(sprintf('%s/l2svm_l2_fa_k_fold.mat',figures_root),'acc','spec','sens','conf','auroc','roc','loo_scores','loo_labels');
figure;hist(loo_acc(:));
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','w')
set(gcf,'color','w');set(gca,'fontweight','b','fontsize',12,'LineWidth',3);xlabel('ACCURACY OF PREDICTION');ylabel('FREQUENCY OF ATTEMPTS');
title(sprintf('SVM-FA, %d-fold, %d iter',M.no_of_folds,M.no_of_iter));
savefig(sprintf('%s/l2svm_l2_fa_k_fold.png',figures_root),0);

%ROC and optimal regions
N=0;P=0;
for i=1:length(loo_labels)
N=length(find(loo_labels{i}<0))+N;
P=length(find(loo_labels{i}>0))+P;
end

%Setting A
m1=tan(45/180*pi);
np=(0.5*(N/P))/m1;
cost=[0 np;0.5 0];
[X1,Y1,T1,auroc1,optrocpt1]=perfcurve(loo_labels(:),loo_scores(:),1,'Cost',cost);
c1=optrocpt1(2)-m1*optrocpt1(1);

%Setting B
m2=tan(60/180*pi);
np=(0.5*(N/P))/m2;
cost=[0 np;0.5 0];
[X2,Y2,T2,auroc2,optrocpt2]=perfcurve(loo_labels(:),loo_scores(:),1,'Cost',cost);
c2=optrocpt2(2)-m2*optrocpt2(1);

%Setting C
m3=tan(30/180*pi);
np=(0.8*(N/P))/m3;
cost=[0 np;0.8 0];
[X3,Y3,T3,auroc3,optrocpt3]=perfcurve(loo_labels(:),loo_scores(:),1,'Cost',cost);
c3=optrocpt3(2)-m3*optrocpt3(1);

%Visualization
figure;hold on;
plot(X1(:,1),Y1(:,1),'b','linewidth',3);
hold on;
plot([0;1],[0,1],'k-.','linewidth',3);
dx=0.1;
idx=find(X1(:,1)==optrocpt1(1)&Y1(:,1)==optrocpt1(2));
plot([X1(idx,1),X1(idx,1)],[0,Y1(idx,1)],'g-.','linewidth',3);
plot([0,X1(idx,1)],[Y1(idx,1),Y1(idx,1)],'g-.','linewidth',3);
x1=max(optrocpt1(1)-0.1,0);x2=optrocpt1(1)+0.1;
plot([x1,x2],[m1*x1+c1,m1*x2+c1],'m-.','linewidth',3);
errorbar(X1(idx,1),Y1(idx,1),Y1(idx,2)-Y1(idx,1),Y1(idx,3)-Y1(idx,1),'xr','linewidth',3);

idx=find(X2(:,1)==optrocpt2(1)&Y2(:,1)==optrocpt2(2));
plot([X2(idx,1),X2(idx,1)],[0,Y2(idx,1)],'g-.','linewidth',3);
plot([0,X2(idx,1)],[Y2(idx,1),Y2(idx,1)],'g-.','linewidth',3);
x1=max(optrocpt2(1)-dx,0);x2=optrocpt2(1)+dx;
plot([x1,x2],[m2*x1+c2,m2*x2+c2],'m-.','linewidth',3);
errorbar(X2(idx,1),Y2(idx,1),Y2(idx,2)-Y2(idx,1),Y2(idx,3)-Y2(idx,1),'xr','linewidth',3);

idx=find(X3(:,1)==optrocpt3(1)&Y3(:,1)==optrocpt3(2));
plot([X3(idx,1),X3(idx,1)],[0,Y3(idx,1)],'g-.','linewidth',3);
plot([0,X3(idx,1)],[Y3(idx,1),Y3(idx,1)],'g-.','linewidth',3);
x1=max(optrocpt3(1)-dx,0);x2=optrocpt3(1)+dx;
plot([x1,x2],[m3*x1+c3,m3*x2+c3],'m-.','linewidth',3);
errorbar(X3(idx,1),Y3(idx,1),Y3(idx,2)-Y3(idx,1),Y3(idx,3)-Y3(idx,1),'xr','linewidth',3);
set(gcf,'color','w');set(gca,'fontweight','b','fontsize',12,'LineWidth',3);xlabel('FALSE POSITIVE RATE');ylabel('TRUE POSITIVE RATE');
title(sprintf('SVM-FA, Total area under ROC curve: %f',auroc1(1)),'fontweight','b','fontsize',12);
savefig(sprintf('%s/l2svm_l2_fa_k_fold_roc_new.png',figures_root),0);
save(sprintf('%s/l2svm_l2_fa_k_fold.mat',figures_root),'loo_acc','loo_scores','loo_labels','X1','Y1','T1','auroc1','optrocpt1','optrocpt2','optrocpt3');

%Manual interaction need to adjust the zoom factors and regions.
savefig(sprintf('%s/l2svm_l2_fa_k_fold_roc_new_zoom1.png',figures_root),0);
savefig(sprintf('%s/l2svm_l2_fa_k_fold_roc_new_zoom2.png',figures_root),0);
savefig(sprintf('%s/l2svm_l2_fa_k_fold_roc_new_zoom3.png',figures_root),0);

%Visualizing the weights.
weights=zeros(M.no_of_folds*M.no_of_iter,numel(loo_models{1}.w));
for i=1:length(loo_models)
weights(i,:)=loo_models{i}.w.^2;
end
mean_w=mean(weights,1);
std_w=std(weights,1);
mean_mean_w=mean(mean_w);
std_mean_w=std(mean_w);
z=(mean_w-mean_mean_w)./std_mean_w;
nii.vol=zeros(size(nii.vol));
nii.vol(wm_idx(fa_idx))=z(1:end-1);
save_nifti(nii,sprintf('%s/l2svm_l2_z_fa.nii',figures_root));
save(sprintf('%s/l2svm_l2_fa_k_fold.mat',figures_root),'acc','spec','sens','conf','auroc','roc','loo_acc','loo_scores','loo_labels','X','Y','T','auroc','optrocpt','weights','z');
rm_idx=find(srs==not_include);
srs(rm_idx)=[];
diagnosis(rm_idx)=[];
fiq(rm_idx)=[];
age(rm_idx)=[];
subjids(rm_idx)=[];
fa(rm_idx,:)=[];

asd_idx=find(diagnosis==1);
srs_asd=srs(asd_idx);
avgIQ_asd=fiq(asd_idx);
age_asd=age(asd_idx);

%Computing test-set bounds
%Shuffling and Splitting.
s=randperm(numel(labels));
labels=labels(s);
data=data(s,:);

train_percent=0.9;%0.8;%0.5;%0.6;%0.7;
train_idx=1:floor(train_percent*numel(labels));
test_idx=length(train_idx)+1:numel(labels);

%Training and testing.
model=train(labels(train_idx),data(train_idx,:),sprintf('-s %d -c %f',M.loss_regularization,M.C));
[train_preds,train_acc,out] = predict(labels(train_idx),data(train_idx,:),model);
if(labels(train_idx(1))==-1)
out=-out;
end
fprintf('Number of training errors: %d\n',sum(sign(labels(train_idx))~=sign(out)));
[test_preds,test_acc,out] = predict(labels(test_idx),data(test_idx,:),model);
if(labels(test_idx(1))==-1)
out=-out;
end
fprintf('Number of test errors: %d out of %d',sum(sign(labels(test_idx))~=sign(out)),length(out));

%Binomial test-set bounds.
a=0.05;
unix(sprintf('%s/upper_bound %f %d %d',bound_root,a,length(out),sum(sign(labels(test_idx))~=sign(out))));
unix(sprintf('%s/lower_bound %f %d %d',bound_root,a,length(out),sum(sign(labels(test_idx))~=sign(out))));

Acknowledgement: This work was supported by the NIMH R01 MH080826 (JEL) and R01 MH084795 (JEL) (University of Utah), the NIH Mental Retardation/ Developmental Disabilities Research Center (MRDDRC Waisman Center), NIMH 62015 (ALA), the Henry M. Jackson Foundation (ALA), the Autism Society of Southwestern Wisconsin, the NCCAM P01 AT004952-04 (RJD and AL) and the Waisman Core grant P30 HD003352-45 (RJD). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Mental Health, the National Institutes of Health or the Waisman Center.