Home > demo2 > demo22.m

demo22

PURPOSE ^

% Nine triplets identified by Baranzini et al. to be discriminative

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

% Nine triplets identified by Baranzini et al. to be discriminative

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% Nine triplets identified by Baranzini et al. to be discriminative
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 % Clean up the data
0027 ix = [153, 211];
0028 Xtmp(ix,:)=[]; id_sbj(ix)=[]; Ytmp(ix)=[]; Ttmp(ix)=[];
0029 
0030 Xtmp(isnan(Xtmp))=0;
0031 
0032 
0033 % Collect subjects into cells
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 %% Reject the last example because it lacks sample for t=0
0052 X(end)=[];
0053 Y(end)=[];
0054 
0055 %% Perform cross-validation
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 %% Plot accuracy
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 %% Plot kernel weights
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')

Generated on Sat 22-Aug-2009 22:15:36 by m2html © 2003