octave_packages/statistics-1.1.3/anovan.m

   1 ## Copyright (C) 2003-2005 Andy Adler <adler@ncf.ca>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{pval}, @var{f}, @var{df_b}, @var{df_e}] =} anovan (@var{data}, @var{grps})
  18 ## @deftypefnx {Function File} {[@var{pval}, @var{f}, @var{df_b}, @var{df_e}] =} anovan (@var{data}, @var{grps}, 'param1', @var{value1})
  19 ## Perform a multi-way analysis of variance (ANOVA).  The goal is to test
  20 ## whether the population means of data taken from @var{k} different
  21 ## groups are all equal.
  22 ##
  23 ## Data is a single vector @var{data} with groups specified by
  24 ## a corresponding matrix of group labels @var{grps}, where @var{grps}
  25 ## has the same number of rows as @var{data}. For example, if
  26 ## @var{data} = [1.1;1.2]; @var{grps}= [1,2,1; 1,5,2];
  27 ## then data point 1.1 was measured under conditions 1,2,1 and
  28 ## data point 1.2 was measured under conditions 1,5,2.
  29 ## Note that groups do not need to be sequentially numbered.
  30 ##
  31 ## By default, a 'linear' model is used, computing the N main effects
  32 ## with no interactions. this may be modified by param 'model'
  33 ##
  34 ## p= anovan(data,groups, 'model', modeltype)
  35 ## - modeltype = 'linear': compute N main effects
  36 ## - modeltype = 'interaction': compute N effects and
  37 ##                               N*(N-1) two-factor interactions
  38 ## - modeltype = 'full': compute interactions at all levels
  39 ##
  40 ## Under the null of constant means, the statistic @var{f} follows an F
  41 ## distribution with @var{df_b} and @var{df_e} degrees of freedom.
  42 ##
  43 ## The p-value (1 minus the CDF of this distribution at @var{f}) is
  44 ## returned in @var{pval}.
  45 ##
  46 ## If no output argument is given, the standard one-way ANOVA table is
  47 ## printed.
  48 ##
  49 ## BUG: DFE is incorrect for modeltypes != full
  50 ## @end deftypefn
  51
  52 ## Author: Andy Adler <adler@ncf.ca>
  53 ## Based on code by: KH <Kurt.Hornik@ci.tuwien.ac.at>
  54 ## $Id: anovan.m 10203 2012-04-12 13:47:32Z carandraug $
  55 ##
  56 ## TESTING RESULTS:
  57 ##  1. ANOVA ACCURACY: www.itl.nist.gov/div898/strd/anova/anova.html
  58 ##     Passes 'easy' test. Comes close on 'Average'. Fails 'Higher'.
  59 ##     This could be fixed with higher precision arithmetic
  60 ##  2. Matlab anova2 test
  61 ##      www.mathworks.com/access/helpdesk/help/toolbox/stats/anova2.html
  62 ##     % From web site:
  63 ##      popcorn= [  5.5 4.5 3.5; 5.5 4.5 4.0; 6.0 4.0 3.0;
  64 ##                  6.5 5.0 4.0; 7.0 5.5 5.0; 7.0 5.0 4.5];
  65 ##     % Define groups so reps = 3
  66 ##      groups = [  1 1;1 2;1 3;1 1;1 2;1 3;1 1;1 2;1 3;
  67 ##                  2 1;2 2;2 3;2 1;2 2;2 3;2 1;2 2;2 3 ];
  68 ##      anovan( vec(popcorn'), groups, 'model', 'full')
  69 ##     % Results same as Matlab output
  70 ##  3. Matlab anovan test
  71 ##      www.mathworks.com/access/helpdesk/help/toolbox/stats/anovan.html
  72 ##    % From web site
  73 ##      y = [52.7 57.5 45.9 44.5 53.0 57.0 45.9 44.0]';
  74 ##      g1 = [1 2 1 2 1 2 1 2];
  75 ##      g2 = {'hi';'hi';'lo';'lo';'hi';'hi';'lo';'lo'};
  76 ##      g3 = {'may'; 'may'; 'may'; 'may'; 'june'; 'june'; 'june'; 'june'};
  77 ##      anovan( y', [g1',g2',g3'])
  78 ##    % Fails because we always do interactions
  79
  80 function [PVAL, FSTAT, DF_B, DFE] = anovan (data, grps, varargin)
  81
  82     if nargin <= 1
  83         usage ("anovan (data, grps)");
  84     end
  85
  86     # test supplied parameters
  87     modeltype= 'linear';
  88     for idx= 3:2:nargin
  89        param= varargin{idx-2};
  90        value= varargin{idx-1};
  91
  92        if strcmp(param, 'model')
  93           modeltype= value;
  94 #      elseif strcmp(param    # add other parameters here
  95        else
  96           error(sprintf('parameter %s is not supported', param));
  97        end
  98     end
  99
 100     if ~isvector (data)
 101           error ("anova: for `anova (data, grps)', data must be a vector");
 102     endif
 103
 104     nd = size (grps,1); # number of data points
 105     nw = size (grps,2); # number of anova "ways"
 106     if (~ isvector (data) || (length(data) ~= nd))
 107       error ("anova: grps must be a matrix of the same number of rows as data");
 108     endif
 109
 110     [g,grp_map]   = relabel_groups (grps);
 111     if strcmp(modeltype, 'linear')
 112        max_interact  = 1;
 113     elseif strcmp(modeltype,'interaction')
 114        max_interact  = 2;
 115     elseif strcmp(modeltype,'full')
 116        max_interact  = rows(grps);
 117     else
 118        error(sprintf('modeltype %s is not supported', modeltype));
 119     end
 120     ng = length(grp_map);
 121     int_tbl       = interact_tbl (nw, ng, max_interact );
 122     [gn, gs, gss] = raw_sums(data, g, ng, int_tbl);
 123
 124     stats_tbl = int_tbl(2:size(int_tbl,1),:)>0;
 125     nstats= size(stats_tbl,1);
 126     stats= zeros( nstats+1, 5); # SS, DF, MS, F, p
 127     for i= 1:nstats
 128         [SS, DF, MS]= factor_sums( gn, gs, gss, stats_tbl(i,:), ng, nw);
 129         stats(i,1:3)= [SS, DF, MS];
 130     end
 131
 132     # The Mean squared error is the data - avg for each possible measurement
 133     # This calculation doesn't work unless there is replication for all grps
 134 #   SSE= sum( gss(sel) ) - sum( gs(sel).^2 ./ gn(sel) );
 135     SST= gss(1) - gs(1)^2/gn(1);
 136     SSE= SST - sum(stats(:,1));
 137     sel = select_pat( ones(1,nw), ng, nw); %incorrect for modeltypes != full
 138     DFE= sum( (gn(sel)-1).*(gn(sel)>0) );
 139     MSE= SSE/DFE;
 140     stats(nstats+1,1:3)= [SSE, DFE, MSE];
 141
 142     for i= 1:nstats
 143         MS= stats(i,3);
 144         DF= stats(i,2);
 145         F= MS/MSE;
 146         pval = 1 - fcdf (F, DF, DFE);
 147         stats(i,4:5)= [F, pval];
 148     end
 149
 150     if nargout==0;
 151         printout( stats, stats_tbl );
 152     else
 153         PVAL= stats(1:nstats,5);
 154         FSTAT=stats(1:nstats,4);
 155         DF_B= stats(1:nstats,2);
 156         DF_E= DFE;
 157     end
 158 endfunction
 159
 160
 161 # relabel groups to a mapping from 1 to ng
 162 # Input
 163 #   grps    input grouping
 164 # Output
 165 #   g       relabelled grouping
 166 #   grp_map map from output to input grouping
 167 function [g,grp_map] = relabel_groups(grps)
 168     grp_vec= vec(grps);
 169     s= sort (grp_vec);
 170     uniq = 1+[0;find(diff(s))];
 171     # mapping from new grps to old groups
 172     grp_map = s(uniq);
 173     # create new group g
 174     ngroups= length(uniq);
 175     g= zeros(size(grp_vec));
 176     for i = 1:ngroups
 177         g( find( grp_vec== grp_map(i) ) ) = i;
 178     end
 179     g= reshape(g, size(grps));
 180 endfunction
 181
 182 # Create interaction table
 183 #
 184 # Input:
 185 #    nw            number of "ways"
 186 #    ng            number of ANOVA groups
 187 #    max_interact  maximum number of interactions to consider
 188 #                  default is nw
 189 function int_tbl =interact_tbl(nw, ng, max_interact)
 190     combin= 2^nw;
 191     inter_tbl= zeros( combin, nw);
 192     idx= (0:combin-1)';
 193     for i=1:nw;
 194        inter_tbl(:,i) = ( rem(idx,2^i) >= 2^(i-1) );
 195     end
 196
 197     # find elements with more than max_interact 1's
 198     idx = ( sum(inter_tbl',1) > max_interact );
 199     inter_tbl(idx,:) =[];
 200     combin= size(inter_tbl,1); # update value
 201
 202     #scale inter_tbl
 203     # use ng+1 to map combinations of groups to integers
 204     # this would be lots easier with a hash data structure
 205     int_tbl = inter_tbl .* (ones(combin,1) * (ng+1).^(0:nw-1) );
 206 endfunction
 207
 208 # Calculate sums for each combination
 209 #
 210 # Input:
 211 #    g             relabelled grouping matrix
 212 #    ng            number of ANOVA groups
 213 #    max_interact
 214 #
 215 # Output (virtual (ng+1)x(nw) matrices):
 216 #    gn            number of data sums in each group
 217 #    gs            sum of data in each group
 218 #    gss           sumsqr of data in each group
 219 function    [gn, gs, gss] = raw_sums(data, g, ng, int_tbl);
 220     nw=    size(g,2);
 221     ndata= size(g,1);
 222     gn= gs= gss=  zeros((ng+1)^nw, 1);
 223     for i=1:ndata
 224         # need offset by one for indexing
 225         datapt= data(i);
 226         idx = 1+ int_tbl*g(i,:)';
 227         gn(idx)  +=1;
 228         gs(idx)  +=datapt;
 229         gss(idx) +=datapt^2;
 230     end
 231 endfunction
 232
 233 # Calcualte the various factor sums
 234 # Input:
 235 #    gn            number of data sums in each group
 236 #    gs            sum of data in each group
 237 #    gss           sumsqr of data in each group
 238 #    select        binary vector of factor for this "way"?
 239 #    ng            number of ANOVA groups
 240 #    nw            number of ways
 241
 242 function [SS,DF]= raw_factor_sums( gn, gs, gss, select, ng, nw);
 243    sel= select_pat( select, ng, nw);
 244    ss_raw=   gs(sel).^2 ./ gn(sel);
 245    SS= sum( ss_raw( ~isnan(ss_raw) ));
 246    if length(find(select>0))==1
 247        DF= sum(gn(sel)>0)-1;
 248    else
 249        DF= 1; #this isn't the real DF, but needed to multiply
 250    end
 251 endfunction
 252
 253 function [SS, DF, MS]= factor_sums( gn, gs, gss, select, ng, nw);
 254    SS=0;
 255    DF=1;
 256
 257    ff = find(select);
 258    lff= length(ff);
 259    # zero terms added, one term subtracted, two added, etc
 260    for i= 0:2^lff-1
 261        remove= find( rem( floor( i * 2.^(-lff+1:0) ), 2) );
 262        sel1= select;
 263        if ~isempty(remove)
 264            sel1( ff( remove ) )=0;
 265        end
 266        [raw_sum,raw_df]= raw_factor_sums(gn,gs,gss,sel1,ng,nw);
 267
 268        add_sub= (-1)^length(remove);
 269        SS+= add_sub*raw_sum;
 270        DF*= raw_df;
 271    end
 272
 273    MS=  SS/DF;
 274 endfunction
 275
 276 # Calcualte the various factor sums
 277 # Input:
 278 #    select        binary vector of factor for this "way"?
 279 #    ng            number of ANOVA groups
 280 #    nw            number of ways
 281 function sel= select_pat( select, ng, nw);
 282    # if select(i) is zero, remove nonzeros
 283    # if select(i) is zero, remove zero terms for i
 284    field=[];
 285
 286    if length(select) ~= nw;
 287        error("length of select must be = nw");
 288    end
 289    ng1= ng+1;
 290
 291    if isempty(field)
 292        # expand 0:(ng+1)^nw in base ng+1
 293        field= (0:(ng1)^nw-1)'* ng1.^(-nw+1:0);
 294        field= rem( floor( field), ng1);
 295        # select zero or non-zero elements
 296        field= field>0;
 297    end
 298    sel= find( all( field == ones(ng1^nw,1)*select(:)', 2) );
 299 endfunction
 300
 301
 302 function printout( stats, stats_tbl );
 303   nw= size( stats_tbl,2);
 304   [jnk,order]= sort( sum(stats_tbl,2) );
 305
 306   printf('\n%d-way ANOVA Table (Factors A%s):\n\n', nw, ...
 307          sprintf(',%c',toascii('A')+(1:nw-1)) );
 308   printf('Source of Variation        Sum Sqr   df      MeanSS    Fval   p-value\n');
 309   printf('*********************************************************************\n');
 310   printf('Error                  %10.2f  %4d %10.2f\n', stats( size(stats,1),1:3));
 311
 312   for i= order(:)'
 313       str=  sprintf(' %c x',toascii('A')+find(stats_tbl(i,:)>0)-1 );
 314       str= str(1:length(str)-2); # remove x
 315       printf('Factor %15s %10.2f  %4d %10.2f  %7.3f  %7.6f\n', ...
 316          str, stats(i,:) );
 317   end
 318   printf('\n');
 319 endfunction
 320
 321 #{
 322 # Test Data from http://maths.sci.shu.ac.uk/distance/stats/14.shtml
 323 data=[7  9  9  8 12 10 ...
 324       9  8 10 11 13 13 ...
 325       9 10 10 12 10 12]';
 326 grp = [1,1; 1,1; 1,2; 1,2; 1,3; 1,3;
 327        2,1; 2,1; 2,2; 2,2; 2,3; 2,3;
 328        3,1; 3,1; 3,2; 3,2; 3,3; 3,3];
 329 data=[7  9  9  8 12 10  9  8 ...
 330       9  8 10 11 13 13 10 11 ...
 331       9 10 10 12 10 12 10 12]';
 332 grp = [1,4; 1,4; 1,5; 1,5; 1,6; 1,6; 1,7; 1,7;
 333        2,4; 2,4; 2,5; 2,5; 2,6; 2,6; 2,7; 2,7;
 334        3,4; 3,4; 3,5; 3,5; 3,6; 3,6; 3,7; 3,7];
 335 # Test Data from http://maths.sci.shu.ac.uk/distance/stats/9.shtml
 336 data=[9.5 11.1 11.9 12.8 ...
 337      10.9 10.0 11.0 11.9 ...
 338      11.2 10.4 10.8 13.4]';
 339 grp= [1:4,1:4,1:4]';
 340 # Test Data from http://maths.sci.shu.ac.uk/distance/stats/13.shtml
 341 data=[7.56  9.68 11.65  ...
 342       9.98  9.69 10.69  ...
 343       7.23 10.49 11.77  ...
 344       8.22  8.55 10.72  ...
 345       7.59  8.30 12.36]';
 346 grp = [1,1;1,2;1,3;
 347        2,1;2,2;2,3;
 348        3,1;3,2;3,3;
 349        4,1;4,2;4,3;
 350        5,1;5,2;5,3];
 351 # Test Data from www.mathworks.com/
 352 #                access/helpdesk/help/toolbox/stats/linear10.shtml
 353 data=[23  27  43  41  15  17   3   9  20  63  55  90];
 354 grp= [ 1    1   1   1   2   2   2   2   3   3   3   3;
 355        1    1   2   2   1   1   2   2   1   1   2   2]';
 356 #}
 357
 358
 359