eeg_regress

 eeg_regress - plot multiple linear regression results
 
 [H] = eeg_regress(data,[robust],[GP],labels)
 
 data is NxP matrix, with N observations, P-1 predictors
 and the last column is the observed values to predict.
 So, if only 2 columns are given it does a simple linear
 regression.
 
 robust = 1, use robust least-squares regression
 robust = 0, use least-squares regression
 note, current robust is forced to 1.
 
 GP is Nx1 grouping variable

 labels = cell strings with variable names, keep to
 5 characters or less for best format of figure text
 
 H is an array of handles to figures
 
 All the results are plotted to new figures, which can
 be saved/exported using print(H(i)...)

0001 function [H] = eeg_regress(data,robust,GP,labels)
0002 
0003 % eeg_regress - plot multiple linear regression results
0004 %
0005 % [H] = eeg_regress(data,[robust],[GP],labels)
0006 %
0007 % data is NxP matrix, with N observations, P-1 predictors
0008 % and the last column is the observed values to predict.
0009 % So, if only 2 columns are given it does a simple linear
0010 % regression.
0011 %
0012 % robust = 1, use robust least-squares regression
0013 % robust = 0, use least-squares regression
0014 % note, current robust is forced to 1.
0015 %
0016 % GP is Nx1 grouping variable
0017 %
0018 % labels = cell strings with variable names, keep to
0019 % 5 characters or less for best format of figure text
0020 %
0021 % H is an array of handles to figures
0022 %
0023 % All the results are plotted to new figures, which can
0024 % be saved/exported using print(H(i)...)
0025 %
0026 
0027 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:33 $
0028 
0029 % Licence:  GNU GPL, no express or implied warranties
0030 % History:  11/2002, Darren.Weber_at_radiology.ucsf.edu
0031 %
0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0033 
0034 
0035 if ~exist('data','var'),
0036   fprintf('...no input data.\n\n');
0037   return
0038 end
0039 
0040 if ~exist('robust','var'),
0041   robust = 1;
0042 else,
0043   warning('must use robust method, at present');
0044   robust = 1; % cheat for now, as some plotting stuff below requires it
0045 end
0046 
0047 if ~exist('labels','var'),
0048   for i = 1:size(data,2),
0049     if i == size(data,2),
0050       labels{i} = 'Y';
0051     else
0052       labels{i} = sprintf('X%d',i);
0053     end
0054   end
0055 end
0056 
0057 if ~exist('GP','var'),
0058   error('GP input variable is undefined\n');
0059 end
0060 
0061 
0062 plotmatrix(data); % to explore relationships
0063 H(1) = gcf;
0064 
0065 
0066 y = data(:,end);           % values to fit from last column
0067 e = ones(length(data),1);  % for constant term in model
0068 X = [e data(:,1:end-1)];   % predictor values from first columns
0069 
0070 
0071 if ~robust,
0072   %beta = X\y;              % 3 x 1 matrix, const, beta1 & beta2
0073   [B,BCI,R,RCI,STATS] = regress(y,X,.05);
0074   P = X*B;      % fitted values
0075   
0076   % use B slope, with BCI intercepts, not strictly correct
0077   CI = [ [BCI(1,1); B(2:end)]  [BCI(1,2); B(2:end)] ];
0078   YCI = [X*CI(:,1),  X*CI(:,2) ];
0079   
0080   % Note if BCI(x,:) < or > 0, B(x) is significant
0081   
0082   % define figure rows & columns
0083   cols = 6;
0084   rows = size(X,2);
0085   
0086 else
0087   [B,BCI,R,RCI,STATS] = regress(y,X,.05);
0088   clear B BCI R RCI;
0089   
0090   R2 = STATS(1);
0091   F  = STATS(2);
0092   Fp = STATS(3); % for total regression
0093   clear STATS;
0094   
0095   [B,STATS] = robustfit(X(:,2:end),y);
0096   
0097   P = X*B;           % predicted values
0098   R = STATS.resid;   % residuals
0099   
0100   BCI = [B - 2*STATS.se, B + 2*STATS.se]; % beta CI
0101   YCI = [P - 2*STATS.s,  P + 2*STATS.s ]; % predicted CI
0102   
0103   % define figure rows & columns
0104   cols = 5;
0105   rows = size(X,2);
0106   
0107 end
0108 
0109 
0110 
0111 for fig = 1:2,
0112   
0113   % Setup figure
0114   H(end+1) = figure('color',[0 0 0]);
0115   pos = get(gcf,'Position');
0116   set(gcf,'Position',[pos(1)-(pos(3)/2) pos(2)-pos(4) pos(3)*2 pos(4)*2]);  % double width/height
0117   
0118   colormap(prism);
0119   fontsize = 8;
0120   fontcolor = [1 1 1]; % white
0121   
0122   % plot observed, predicted and 95% CI
0123   if fig < 2,
0124     subplot(2,3,1);
0125   else
0126     subplot(rows,cols,1);
0127   end
0128   if max(GP),
0129     scatter(P,y,5,GP, 'filled'); hold on;
0130   else,
0131     scatter(P,y,5,'m','filled'); hold on;
0132   end
0133   %scatter(P,P,5,'w','filled');
0134   plot(P,P,'w-',...
0135     P,YCI(:,1),'w:',...
0136     P,YCI(:,2),'w:');
0137   %legend('predicted','95% CI',0);
0138   set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
0139   
0140   title{1} = sprintf('y'' = %8.2f',B(1));
0141   for j = 2:length(B),
0142     title{1} = sprintf('%s + %6.2f*X%d',title{1},B(j),j-1);
0143   end
0144   Htitle = get(gca,'title');
0145   set(Htitle,'string',title{1}, 'fontsize',fontsize,'color',fontcolor);
0146   
0147   xlab = get(gca,'xlabel');
0148   set(xlab,'string','Y predicted','fontsize',fontsize,'color',fontcolor);
0149   ylab = get(gca,'ylabel');
0150   set(ylab,'string','Y observed', 'fontsize',fontsize,'color',fontcolor);
0151   
0152   
0153   
0154   % plot observed, predicted and 95% CI against Xmean
0155   if fig < 2,
0156     subplot(2,3,2);
0157   else
0158     subplot(rows,cols,2);
0159   end
0160   % calculate mean X values
0161   Xmean = mean(X(:,2:end),2);
0162   % plot the observed & predicted values
0163   if max(GP),
0164     scatter(Xmean,y,5,GP, 'filled'); hold on;
0165   else,
0166     scatter(Xmean,y,5,'m','filled'); hold on;
0167   end
0168   scatter(Xmean,P,5,'w','filled');
0169   set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
0170   xlab = get(gca,'xlabel');
0171   set(xlab,'string','mean Xi', 'fontsize',fontsize,'color',fontcolor);
0172   ylab = get(gca,'ylabel');
0173   set(ylab,'string','Y observed','fontsize',fontsize,'color',fontcolor);
0174   
0175   
0176   % curve fitting, cubic polynomial
0177   [ Xsort, k ] = sort(Xmean);
0178   [p,S] = polyfit(Xsort,P(k),3);
0179   polfit = polyval(p,Xsort,S);
0180   % % plot observed, polyfit and CI
0181   plot(Xsort,polfit,'w-');
0182   % legend('data','polyfit',0);
0183   
0184   
0185   % plot residuals against predicted
0186   if fig < 2,
0187     subplot(2,3,4);
0188   else
0189     subplot(rows,cols,3);
0190   end
0191   if max(GP),
0192     scatter(P,R,5,GP,'filled'); hold on
0193   else,
0194     scatter(P,R,5,'m','filled'); hold on
0195   end
0196   line(P,zeros(size(R)),'linestyle','-','color','w');
0197   set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
0198   xlab = get(gca,'xlabel');
0199   set(xlab,'string','Y predicted','fontsize',fontsize,'color',fontcolor);
0200   ylab = get(gca,'ylabel');
0201   set(ylab,'string','residual', 'fontsize',fontsize,'color',fontcolor);
0202   
0203   % plot normal QQ plot of residuals
0204   if fig < 2,
0205     subplot(2,3,5);
0206   else
0207     subplot(rows,cols,4);
0208   end
0209   qqplot(R);
0210   QQtitle = get(gca,'title');
0211   set(QQtitle,'string','QQ plot','fontsize',fontsize,'color',fontcolor);
0212   xlab = get(gca,'xlabel');
0213   set(xlab,'string','','fontsize',fontsize,'color',fontcolor);
0214   ylab = get(gca,'ylabel');
0215   set(ylab,'string','','fontsize',fontsize,'color',fontcolor);
0216   set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
0217   
0218   
0219   
0220   
0221   
0222   % print model statistics
0223   if robust,
0224     
0225     %    [B,STATS] = ROBUSTFIT(...) also returns a STATS structure
0226     %     containing the following fields:
0227     %         stats.ols_s     sigma estimate (rmse) from least squares fit
0228     %         stats.robust_s  robust estimate of sigma
0229     %         stats.mad_s     MAD estimate of sigma; used for scaling
0230     %                         residuals during the iterative fitting
0231     %         stats.s         final estimate of sigma, the larger of robust_s
0232     %                         and a weighted average of ols_s and robust_s
0233     %         stats.se        standard error of coefficient estimates
0234     %         stats.t         ratio of b to stats.se
0235     %         stats.p         p-values for stats.t
0236     %         stats.coeffcorr estimated correlation of coefficient estimates
0237     %         stats.w         vector of weights for robust fit
0238     %         stats.h         vector of leverage values for least squares fit
0239     %         stats.dfe       degrees of freedom for error
0240     %         stats.R         R factor in QR decomposition of X matrix
0241     
0242     if fig < 2,
0243       subplot(2,3,3);
0244     else
0245       subplot(rows,cols,5);
0246     end
0247     axis off;
0248     txt{1} = sprintf('%6s %8s %8s','','Beta','SE');
0249     for i = 1:length(B),
0250       if i == 1,
0251         predictor = 'CONST';
0252       else
0253         if i-1 <= length(labels),
0254           predictor = labels{i-1};
0255         end
0256       end
0257       txt{i+1} = sprintf('%6s %8.4f %8.4f',predictor,B(i),STATS.se(i));
0258     end
0259     text(-.2,.5,txt,'Fontname','courier','fontsize',fontsize,'color',fontcolor)
0260     
0261     txt = '';
0262     txt{1} = sprintf('%6s %8s %6s %6s %5s','','F','DF','sig','r^2');
0263     txt{2} = sprintf('%6s %8.4f %6.2f %6.4f %4.2f','ALL',F,STATS.dfe,Fp,R2);
0264     for i = 1:length(B),
0265       if i == 1,
0266         predictor = 'CONST';
0267       else
0268         if i-1 <= length(labels),
0269           predictor = labels{i-1};
0270         end
0271       end
0272       txt{i+2} = sprintf('%6s %8.4f %6.2f %6.4f',predictor,STATS.t(i),STATS.dfe,STATS.p(i));
0273     end
0274     
0275     % output the correlation matrix
0276     datacorr = corrcoef(data);
0277     rtxt = num2str(datacorr,' %+5.2f');
0278     txt2{1} = '    Correlation Matrix';
0279     txt2{2} = '   ';
0280     for r = 1:length(datacorr),
0281       if r < length(datacorr),
0282         txt2{2}   = sprintf('%5s %5s',txt2{2},sprintf('X%d',r));
0283         txt2{r+2} = sprintf('%5s %s',labels{r},rtxt(r,:));
0284       else
0285         txt2{2}   = sprintf('%5s %5s',txt2{2},'Y');
0286         txt2{r+2} = sprintf('%5s %s',labels{r},rtxt(r,:));
0287       end
0288     end
0289     
0290     if fig < 2,
0291       text(-.2,.1,txt ,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
0292       text(-.1,1 ,txt2,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
0293     else,
0294       text(-.2,.1,txt ,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
0295       text(-.1,1 ,txt2,'Fontname','courier','fontsize',fontsize,'color',fontcolor);
0296     end
0297     txt = '';
0298     txt2 = '';
0299   else
0300     
0301     % stats(1) = R^2, stats(2) = F, stats(3) = p; % for total regression
0302     if fig < 2,
0303       subplot(2,3,3);
0304     else
0305       subplot(rows,cols,5);
0306     end
0307     axis off;
0308     txt{1} = sprintf('%12s %10s  %10s  %10s\n','','R^2','F','sig');
0309     txt{2} = sprintf('%12s%10.4f  %10.4f  %10.4f','MODEL',STATS(1),STATS(2),STATS(3));
0310     text(-.2,.9,txt,'Fontname','courier','fontsize',fontsize,'color',fontcolor)
0311     
0312     if fig > 1,
0313       subplot(rows,cols,6); rcoplot(R, RCI);
0314     end
0315   end
0316   
0317 end
0318 
0319 
0320 
0321 
0322 
0323 % plot observed, predicted and 95% CI against individual predictors
0324 for i = 2:size(X,2),
0325   subplot(rows,cols,cols*i-(cols-1));
0326   % plot the observed values
0327   if max(GP),
0328     scatter(X(:,i),y,5,GP, 'filled'); hold on;
0329   else
0330     scatter(X(:,i),y,5,'w','filled'); hold on;
0331   end
0332   set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
0333   
0334   % Calculate predicted values for just this predictor
0335   beta = [B(1);B(i)];
0336   P = [X(:,1),X(:,i)]*beta;
0337   R = y - P;
0338   YCI = [P - 2*STATS.s,  P + 2*STATS.s ]; % predicted CI
0339   
0340   %scatter(X(:,i),P,5,'w','filled');
0341   plot(X(:,i),P,'w-',...
0342     X(:,i),YCI(:,1),'w:',...
0343     X(:,i),YCI(:,2),'w:');
0344   %legend('predicted','95% CI',0);
0345   
0346   title{1} = sprintf('y'' = %8.2f',B(1));
0347   title{1} = sprintf('%s + %6.2f*{\\bfX%d}',title{1},B(i),i-1);
0348   %Htitle = get(gca,'title');
0349   %set(Htitle,'string',title{1},'interpreter','tex');
0350   ylab = get(gca,'ylabel');
0351   set(ylab,'string',title{1},'fontsize',fontsize,'color',fontcolor);
0352   
0353   
0354   
0355   
0356   % plot residuals against predicted
0357   subplot(rows,cols,cols*i-(cols-3));
0358   if max(GP),
0359     scatter(X(:,i),R,5,GP, 'filled'); hold on;
0360   else
0361     scatter(X(:,i),R,5,'w','filled'); hold on;
0362   end
0363   line(X(:,i),zeros(size(R)),'linestyle','-','color','w');
0364   %xlab = get(gca,'xlabel');
0365   %set(xlab,'string','Xi','fontsize',fontsize,'color',fontcolor);
0366   ylab = get(gca,'ylabel');
0367   set(ylab,'string','residual','fontsize',fontsize,'color',fontcolor);
0368   set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
0369   
0370   
0371   
0372   % plot normal QQ plot of residuals
0373   subplot(rows,cols,cols*i-(cols-4));
0374   qqplot(R);
0375   QQtitle = get(gca,'title');
0376   set(QQtitle,'string','QQ plot','fontsize',fontsize,'color',fontcolor);
0377   xlab = get(gca,'xlabel');
0378   set(xlab,'string','','fontsize',fontsize,'color',fontcolor);
0379   ylab = get(gca,'ylabel');
0380   set(ylab,'string','','fontsize',fontsize,'color',fontcolor);
0381   set(gca,'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1]);
0382 end
0383 
0384 
0385 return