1 ## Copyright (C) 1996-2012 Kurt Hornik
3 ## This file is part of Octave.
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.
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.
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/>.
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.
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}).
31 ## The LR test statistic (Wilks' Lambda) and approximate p-values are
32 ## computed and displayed.
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.)
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)
43 function manova (x, g)
50 error ("manova: Y must not be a vector");
55 if (!isvector (g) || (length (g) != n))
56 error ("manova: G must be a vector of length rows (Y)");
60 i = find (s (2:n) > s(1:(n-1)));
64 error ("manova: there should be at least 2 groups");
66 group_label = s ([1, (reshape (i, 1, k - 1) + 1)]);
69 x = x - ones (n, 1) * mean (x);
75 v = x (find (g == group_label (i)), :);
77 SSB = SSB + s' * s / rows (v);
84 l = real (eig (SSB / SSW));
86 if (isa (l, "single"))
87 l (l < eps ("single")) = 0;
95 Lambda = prod (1 ./ (1 + l));
97 delta = n_w + n_b - (p + n_b + 1) / 2;
99 W_pval_1 = 1 - chi2cdf (- delta * log (Lambda), df_num);
104 eta = sqrt ((p^2 * n_b^2 - 4) / (p^2 + n_b^2 - 5));
107 df_den = delta * eta - df_num / 2 + 1;
109 WT = exp (- log (Lambda) / eta) - 1;
110 W_pval_2 = 1 - fcdf (WT * df_den / df_num, df_num, df_den);
114 ## Hotelling-Lawley Test
115 ## =====================
119 theta = min (p, n_b);
120 u = (abs (p - n_b) - 1) / 2;
121 v = (n_w - p - 1) / 2;
123 df_num = theta * (2 * u + theta + 1);
124 df_den = 2 * (theta * v + 1);
126 HL_pval = 1 - fcdf (HL * df_den / df_num, df_num, df_den);
131 PB = sum (l ./ (1 + l));
133 df_den = theta * (2 * v + theta + 1);
134 PB_pval = 1 - fcdf (PB * df_den / df_num, df_num, df_den);
137 printf ("One-way MANOVA Table:\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);
150 printf ("MANOVA Results:\n");
152 printf ("# of groups: %d\n", k);
153 printf ("# of samples: %d\n", n);
154 printf ("# of variables: %d\n", p);
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);