octave_packages/m/statistics/tests/manova.m

   1 ## Copyright (C) 1996-2012 Kurt Hornik
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} manova (@var{x}, @var{g})
  21 ## Perform a one-way multivariate analysis of variance (MANOVA).  The
  22 ## goal is to test whether the p-dimensional population means of data
  23 ## taken from @var{k} different groups are all equal.  All data are
  24 ## assumed drawn independently from p-dimensional normal distributions
  25 ## with the same covariance matrix.
  26 ##
  27 ## The data matrix is given by @var{x}.  As usual, rows are observations
  28 ## and columns are variables.  The vector @var{g} specifies the
  29 ## corresponding group labels (e.g., numbers from 1 to @var{k}).
  30 ##
  31 ## The LR test statistic (Wilks' Lambda) and approximate p-values are
  32 ## computed and displayed.
  33 ## @end deftypefn
  34
  35 ## Three test statistics (Wilks, Hotelling-Lawley, and Pillai-Bartlett)
  36 ## and corresponding approximate p-values are calculated and displayed.
  37 ## (Currently NOT because the fcdf respectively betai code is too bad.)
  38
  39 ## Author: TF <Thomas.Fuereder@ci.tuwien.ac.at>
  40 ## Adapted-By: KH <Kurt.Hornik@wu-wien.ac.at>
  41 ## Description: One-way multivariate analysis of variance (MANOVA)
  42
  43 function manova (x, g)
  44
  45   if (nargin != 2)
  46     print_usage ();
  47   endif
  48
  49   if (isvector (x))
  50     error ("manova: Y must not be a vector");
  51   endif
  52
  53   [n, p] = size (x);
  54
  55   if (!isvector (g) || (length (g) != n))
  56     error ("manova: G must be a vector of length rows (Y)");
  57   endif
  58
  59   s = sort (g);
  60   i = find (s (2:n) > s(1:(n-1)));
  61   k = length (i) + 1;
  62
  63   if (k == 1)
  64     error ("manova: there should be at least 2 groups");
  65   else
  66     group_label = s ([1, (reshape (i, 1, k - 1) + 1)]);
  67   endif
  68
  69   x = x - ones (n, 1) * mean (x);
  70   SST = x' * x;
  71
  72   s = zeros (1, p);
  73   SSB = zeros (p, p);
  74   for i = 1 : k;
  75     v = x (find (g == group_label (i)), :);
  76     s = sum (v);
  77     SSB = SSB + s' * s / rows (v);
  78   endfor
  79   n_b = k - 1;
  80
  81   SSW = SST - SSB;
  82   n_w = n - k;
  83
  84   l = real (eig (SSB / SSW));
  85
  86   if (isa (l, "single"))
  87     l (l < eps ("single")) = 0;
  88   else
  89     l (l < eps) = 0;
  90   endif
  91
  92   ## Wilks' Lambda
  93   ## =============
  94
  95   Lambda = prod (1 ./ (1 + l));
  96
  97   delta = n_w + n_b - (p + n_b + 1) / 2;
  98   df_num = p * n_b;
  99   W_pval_1 = 1 - chi2cdf (- delta * log (Lambda), df_num);
 100
 101   if (p < 3)
 102     eta = p;
 103   else
 104     eta = sqrt ((p^2 * n_b^2 - 4) / (p^2 + n_b^2 - 5));
 105   endif
 106
 107   df_den = delta * eta - df_num / 2 + 1;
 108
 109   WT = exp (- log (Lambda) / eta) - 1;
 110   W_pval_2 = 1 - fcdf (WT * df_den / df_num, df_num, df_den);
 111
 112   if (0)
 113
 114     ## Hotelling-Lawley Test
 115     ## =====================
 116
 117     HL = sum (l);
 118
 119     theta = min (p, n_b);
 120     u = (abs (p - n_b) - 1) / 2;
 121     v = (n_w - p - 1) / 2;
 122
 123     df_num = theta * (2 * u + theta + 1);
 124     df_den = 2 * (theta * v + 1);
 125
 126     HL_pval = 1 - fcdf (HL * df_den / df_num, df_num, df_den);
 127
 128     ## Pillai-Bartlett
 129     ## ===============
 130
 131     PB = sum (l ./ (1 + l));
 132
 133     df_den = theta * (2 * v + theta + 1);
 134     PB_pval = 1 - fcdf (PB * df_den / df_num, df_num, df_den);
 135
 136     printf ("\n");
 137     printf ("One-way MANOVA Table:\n");
 138     printf ("\n");
 139     printf ("Test             Test Statistic      Approximate p\n");
 140     printf ("**************************************************\n");
 141     printf ("Wilks            %10.4f           %10.9f \n", Lambda, W_pval_1);
 142     printf ("                                      %10.9f \n", W_pval_2);
 143     printf ("Hotelling-Lawley %10.4f           %10.9f \n", HL, HL_pval);
 144     printf ("Pillai-Bartlett  %10.4f           %10.9f \n", PB, PB_pval);
 145     printf ("\n");
 146
 147   endif
 148
 149   printf ("\n");
 150   printf ("MANOVA Results:\n");
 151   printf ("\n");
 152   printf ("# of groups:    %d\n", k);
 153   printf ("# of samples:   %d\n", n);
 154   printf ("# of variables: %d\n", p);
 155   printf ("\n");
 156   printf ("Wilks' Lambda:  %5.4f\n", Lambda);
 157   printf ("Approximate p:  %10.9f (chisquare approximation)\n", W_pval_1);
 158   printf ("                 %10.9f (F approximation)\n", W_pval_2);
 159   printf ("\n");
 160
 161 endfunction