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 >
Wrap
Text File
|
1997-02-26
|
5KB
|
155 lines
## Copyright (C) 1996, 1997 Kurt Hornik
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this file. If not, write to the Free Software Foundation,
## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
## usage: manova (Y, g)
##
## Performs a one-way multivariate analysis of variance (MANOVA). The
## goal is to test whether the p-dimensional population means of data
## taken from k different groups are all equal. All data are assumed
## drawn independently from p-dimensional normal distributions with the
## same covariance matrix.
##
## Y is the data matrix. As usual, rows are observations and columns
## are variables. g is the vector of corresponding group labels (e.g.,
## numbers from 1 to k), so that necessarily, length (g) must be the
## same as rows (Y).
##
## The LR test statistic (Wilks' Lambda) and approximate p-values are
## computed and displayed.
## Three test statistics (Wilks, Hotelling-Lawley, and Pillai-Bartlett)
## and corresponding approximate p-values are calculated and displayed.
## (Currently NOT because the f_cdf respectively betai code is too bad.)
## Author: TF <Thomas.Fuereder@ci.tuwien.ac.at>
## Adapted-By: KH <Kurt.Hornik@ci.tuwien.ac.at>
## Description: One-way multivariate analysis of variance (MANOVA)
function manova (Y, g)
if (nargin != 2)
usage ("manova (Y, g)");
endif
if (is_vector (Y))
error ("manova: Y must not be a vector");
endif
[n, p] = size (Y);
if (!is_vector (g) || (length (g) != n))
error ("manova: g must be a vector of length rows (Y)");
endif
s = sort (g);
i = find (s (2:n) > s(1:(n-1)));
k = length (i) + 1;
if (k == 1)
error ("manova: there should be at least 2 groups");
else
group_label = s ([1, reshape (i, 1, k - 1) + 1]);
endif
Y = Y - ones (n, 1) * mean (Y);
SST = Y' * Y;
s = zeros (1, p);
SSB = zeros (p, p);
for i = 1 : k;
v = Y (find (g == group_label (i)), :);
s = sum (v);
SSB = SSB + s' * s / rows (v);
endfor
n_b = k - 1;
SSW = SST - SSB;
n_w = n - k;
l = real (eig (SSB / SSW));
l (l < eps) = 0;
## Wilks' Lambda
## =============
Lambda = prod (1 ./ (1 + l));
delta = n_w + n_b - (p + n_b + 1) / 2
df_num = p * n_b
W_pval_1 = 1 - chisquare_cdf (- delta * log (Lambda), df_num);
if (p < 3)
eta = p;
else
eta = sqrt ((p^2 * n_b^2 - 4) / (p^2 + n_b^2 - 5))
endif
df_den = delta * eta - df_num / 2 + 1
WT = exp (- log (Lambda) / eta) - 1
W_pval_2 = 1 - f_cdf (WT * df_den / df_num, df_num, df_den);
if (0)
## Hotelling-Lawley Test
## =====================
HL = sum (l);
theta = min (p, n_b);
u = (abs (p - n_b) - 1) / 2;
v = (n_w - p - 1) / 2;
df_num = theta * (2 * u + theta + 1);
df_den = 2 * (theta * v + 1);
HL_pval = 1 - f_cdf (HL * df_den / df_num, df_num, df_den);
## Pillai-Bartlett
## ===============
PB = sum (l ./ (1 + l));
df_den = theta * (2 * v + theta + 1);
PB_pval = 1 - f_cdf (PB * df_den / df_num, df_num, df_den);
printf ("\n");
printf ("One-way MANOVA Table:\n");
printf ("\n");
printf ("Test Test Statistic Approximate p\n");
printf ("**************************************************\n");
printf ("Wilks %10.4f %10.9f \n", Lambda, W_pval_1);
printf (" %10.9f \n", W_pval_2);
printf ("Hotelling-Lawley %10.4f %10.9f \n", HL, HL_pval);
printf ("Pillai-Bartlett %10.4f %10.9f \n", PB, PB_pval);
printf ("\n");
endif
printf ("\n");
printf ("MANOVA Results:\n");
printf ("\n");
printf ("# of groups: %d\n", k);
printf ("# of samples: %d\n", n);
printf ("# of variables: %d\n", p);
printf ("\n");
printf ("Wilks' Lambda: %5.4f\n", Lambda);
printf ("Approximate p: %10.9f (chisquare approximation)\n", W_pval_1);
printf (" %10.9f (F approximation)\n", W_pval_2);
printf ("\n");
endfunction