0001
0002 ixtrip = {[36, 44, 49],...
0003 [36 37 54],...
0004 [36 41 54],...
0005 [4 36 54],...
0006 [44 54 62],...
0007 [10 44 62],...
0008 [28 36 52],...
0009 [4 36 37],...
0010 [36 54 56]};
0011
0012
0013 addpath ../SpicyMKL/
0014
0015 S=importdata('journal.pbio.0030002.sd001.xls');
0016
0017 Ytmp=cell2mat(foreach(@isequal,S.textdata.Sheet1(2:end,1),[],'good'));
0018
0019 Gene=S.textdata.Sheet1(1,7:end-6);
0020
0021 id_sbj=S.data.Sheet1(:,1); [tmp,I]=unique(id_sbj); code=id_sbj(sort(I));
0022
0023 Ttmp=S.textdata.Sheet1(2:end,3);
0024 Xtmp=S.data.Sheet1(:,6:end-6);
0025
0026
0027 ix = [153, 211];
0028 Xtmp(ix,:)=[]; id_sbj(ix)=[]; Ytmp(ix)=[]; Ttmp(ix)=[];
0029
0030 Xtmp(isnan(Xtmp))=0;
0031
0032
0033
0034 Y=zeros(1,length(code));
0035 X=cell(1,length(code));
0036 T=cell(1,length(code));
0037 L=zeros(1,length(code));
0038 for ii=1:length(code)
0039 I=find(id_sbj==code(ii));
0040 Y(ii)=mean(Ytmp(I));
0041 X{ii}=Xtmp(I,:);
0042 T{ii}=Ttmp(I);
0043 L(ii)=length(I);
0044 if ~isequal(T{ii}{1},'t0')
0045 display(I);
0046 display(T{ii});
0047 warning('Index seems to be wrong!');
0048 end
0049 end
0050
0051
0052 X(end)=[];
0053 Y(end)=[];
0054
0055
0056 xTrials = [10 4];
0057 lambda=exp(linspace(log(0.8),log(0.008),10));
0058 [Am, As,divTr,divTe,output,memo]=xvalidation(X,Y,'SpicyMKL',lambda,xTrials,'preproc','best_triplets','ixtr',ixtrip);
0059
0060
0061
0062
0063 figure, h=errorbar(log(lambda), Am, As);
0064 set(gca,'fontsize',16);
0065 set(h,'linewidth',2);
0066 grid on;
0067 logticks('x');
0068 set(gca,'position',[0.13, 0.2 0.775, 0.75])
0069 xlabel('Regularization constant');
0070 ylabel('Accuracy');
0071 format_ticks(gca,get(gca,'xticklabel'))
0072 xlim(log([0.007 1]))
0073
0074 hold on;
0075 plot(xlim,0.878*ones(1,2),'--','color',[0 .5 0],'linewidth',2);
0076
0077 legend('SpicyMKl','Single best (Baranzini et al.)')
0078
0079
0080 ds=zeros(prod(xTrials), length(ixtrip));
0081 for kk=1:prod(xTrials)
0082 [ii,jj]=ind2sub(xTrials,kk);
0083 ds(kk,:)=memo(ii,jj,5).C.d;
0084 end
0085
0086 for ii=1:length(ixtrip), lab{ii}=[Gene{ixtrip{ii}(1)} '/' Gene{ixtrip{ii}(2)} '/' Gene{ixtrip{ii}(3)} ' ']; end
0087 figure, plot(mean(ds),'linewidth',2);
0088 set(gca,'fontsize',12)
0089 set(gca,'position',[0.25 0.4 0.7 0.5])
0090 set(gca,'xtick',1:length(ixtrip))
0091 format_ticks(gca,lab,[],[],[],45);
0092 grid on;
0093 title('Kernel weights')
0094 set(gcf,'paperpositionmode','auto')