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)
Sample Results:
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))));