octave_packages/nan-2.5.5/bland_altman.m

   1 function RES = bland_altman(data,group,arg3)
   2 % BLAND_ALTMANN shows the Bland-Altman plot of two columns of measurements
   3 %   and computes several summary results.
   4 %
   5 %   bland_altman(m1, m2 [,group])
   6 %   bland_altman(data [, group])
   7 %   R = bland_altman(...)
   8 %
   9 %   m1,m2 are two colums with the same number of elements
  10 %       containing the measurements. m1,m2 can be also combined
  11 %       in a single two column data matrix.
  12 %   group [optional] indicates which measurements belong to the same group
  13 %       This is useful to account for repeated measurements.
  14 %
  15 %
  16 % References:
  17 % [1] JM Bland and DG Altman, Measuring agreement in method comparison studies.
  18 %       Statistical Methods in Medical Research, 1999; 8; 135.
  19 %       doi:10.1177/09622802990080204
  20 % [2] P.S. Myles, Using the Bland– Altman method to measure agreement with repeated measures
  21 %       British Journal of Anaesthesia 99(3):309–11 (2007)
  22 %       doi:10.1093/bja/aem214
  23
  24 %       $Id$
  25 %       Copyright (C) 2010,2011 by Alois Schloegl <alois.schloegl@gmail.com>
  26 %       This function is part of the NaN-toolbox
  27 %       http://pub.ist.ac.at/~schloegl/matlab/NaN/
  28
  29 % This program is free software; you can redistribute it and/or
  30 % modify it under the terms of the GNU General Public License
  31 % as published by the Free Software Foundation; either version 3
  32 % of the  License, or (at your option) any later version.
  33 %
  34 % This program is distributed in the hope that it will be useful,
  35 % but WITHOUT ANY WARRANTY; without even the implied warranty of
  36 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  37 % GNU General Public License for more details.
  38 %
  39 % You should have received a copy of the GNU General Public License
  40 % along with this program; if not, write to the Free Software
  41 % Foundation, Inc., 51 Franklin Street - Fifth Floor, Boston, MA 02110-1301, USA.
  42
  43 if nargin<2, group = []; end;
  44 if nargin<3, arg3 = []; end;
  45
  46 if (size(data,2)==1)
  47         data  = [data, group];
  48         group = arg3;
  49 end;
  50
  51
  52 D = data * [1;-1];
  53 M = data * [1;1]/2;
  54
  55 RES.corrcoef = corrcoef(data(:,1),data(:,2),'spearman');
  56 [REs.cc,RES.p] = corrcoef(M,D,'spearman');
  57 if (RES.p<0.05)
  58         warning('A regression model according to section 3.2 [1] should be used');
  59         %% TODO: implement support for this type of data.
  60         RES.a = [ones(size(data,1),1),D]\M;
  61         RES.b = [ones(size(data,1),1),M]\D;
  62 end;
  63
  64 if isempty(group)
  65         G = [1:size(data,1)]';
  66         m = ones(size(data,1),1);
  67         d = D;
  68         RES.Bias = mean(d,1);
  69         RES.Var  = var(d);
  70
  71 elseif ~isempty(group)
  72         %% TODO: this is not finished
  73         warning('analysis of data with repetitions is experimental - it might yield incorrect results - you are warned.!')
  74         [G,I,J] = unique (group);
  75         R       = zeros(size(data));
  76         m       = repmat(NaN,length(G),1);
  77         n       = repmat(NaN,length(G),1);
  78         d       = repmat(NaN,length(G),1);
  79         d2      = repmat(NaN,length(G),1);
  80         data2   = repmat(NaN,length(G),size(data,2));
  81         SW2     = repmat(NaN,length(G),size(data,2));
  82         for i   = 1:length(G),
  83                 ix         = find(group==G(i));
  84                 n(i)       = length(ix);
  85 %               IX((i-1)*N+1:i*N) = ix(ceil(rand(N,1)*n(i)));
  86
  87                 [R(ix,:), data2(i,:)] = center(data(ix,:),1);
  88                 d(i)       = mean(D(ix,:),1);
  89                 m(i)       = mean(M(ix,:),1);
  90                 d2(i)      = mean(D(ix,:).^2,1);
  91                 RES.SW2(i,:)   = var(data(ix,:),[],1);
  92                 RES.avg(i,:)   = mean(data(ix,:),1);
  93         end;
  94
  95         W = 1./n(J);
  96         RES.SSW = sumskipnan(R.^2,1,W);
  97         RES.SSB = var(data,[],1,W)*sum(W)*(sum(W)-1);
  98         RES.sigma2_w= RES.SSW/(sum(W)*(length(G)-1));
  99         RES.sigma2_u= RES.SSB/(sum(W)*(length(G)-1)) - RES.sigma2_w/(length(G));
 100         RES.group = bland_altman(data2);        % FIXME: this plot shows incorrect interval, it does not account for the group/repeated samples.
 101         RES.repeatability_coefficient1 = 2.77*sqrt(var(R,1,1));         % variance with factor group removed
 102         RES.repeatability_coefficient = 2.77*sqrt(mean(SW2,1));         % variance with factor group removed
 103         RES.std_d_ = std(d);
 104         RES.std_D_ = std(D);
 105         RES.std_m_ = std(m);
 106
 107         RES.n = n;
 108         return;
 109
 110         D = d;
 111         M = m;
 112 %       RES.sigma2_dw =
 113
 114         RES.Bias = mean(d,1,[],n);
 115 end;
 116
 117
 118 plot(M,D,'o', [min(M),max(M)]', [0,0]','k--', [min(M),max(M)]', [1,1,1; 0,1.96,-1.96]'*[RES.Bias;std(D)]*[1,1], 'k-');
 119 xlabel('mean');
 120 ylabel('difference');
 121