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