roc_post_regress.m 1.7 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182
  1. function [roc_p,roc_v,D]=roc_post_regress(Y, X, G, varargin)
  2. % [roc_p,roc_v,D]=roc_post_regress(Y, X, G, varargin)
  3. % Performs a sliding ROC analysis on the residuals of the regression of Y
  4. % with X after excluding non-significant regressors [at p<0.01 by default].
  5. % Y an M x N data matrix
  6. % X a cell array of M x N data matrices.
  7. % G an M x 1 binary array, describing which rows of Y belong to which
  8. % category
  9. %
  10. % Note: the evaluation of the significance of regressors assumes that the
  11. % bins of Y are independent. If the bins of Y are not independent because
  12. % of smoothing, the regressors are more likely to be significant.
  13. pairs={'alph' 0.01;
  14. 'quick', 0;
  15. };
  16. parseargs(varargin, pairs);
  17. [m,n]=size(Y);
  18. Y=Y(:);
  19. if iscell(X)
  20. Xp=ones(numel(X{1}), numel(X)+1);
  21. for cx=1:numel(X)
  22. Xp(:,cx)=X{cx}(:);
  23. end
  24. else
  25. Xp=X;
  26. end
  27. clear X
  28. % do the regression with p<0.01 for each coefficient
  29. [B,BINT,R,RINT,STATS] = regress(Y,Xp,alph);
  30. % remove parameters that are not significant
  31. keeps=prod(BINT,2)>0;
  32. keeps(end)=true; % This is probably not necessary, but just in case, we don't want to break the regression.
  33. if all(keeps==false)
  34. fprintf('No regressors are significant. Performing ROC on Y\n');
  35. R=Y;
  36. elseif any(keeps==false)
  37. Xpr=Xp(:,keeps);
  38. [B,BINT,R,RINT,STATS] = regress(Y,Xpr,alph);
  39. end
  40. Rm=reshape(R,m,n);
  41. aM=Rm(G==0,:);
  42. bM=Rm(G==1,:);
  43. if quick
  44. % This is MUCH faster, but less accurate.
  45. roc_v=auc(aM,bM);
  46. [~,roc_p]=ttest2(aM,bM);
  47. else
  48. [roc_v,roc_p]=slidingROC(aM,bM);
  49. end
  50. D.keeps=keeps;
  51. D.B=B;
  52. D.BINT=BINT;
  53. D.Rm=Rm;
  54. D.RINT=RINT;
  55. D.STATS=STATS;
  56. D.roc_v=roc_v;
  57. D.roc_p=roc_p;
  58. [roc_p,pi]=min(roc_p);
  59. roc_v=roc_v(pi);