Code covered by the BSD License  

Highlights from
welchanova

welchanova

by

 

11 Jun 2012 (Updated )

Welch ANOVA Test for Unequal Variances.

welchanova.m
function welchanova(x,alpha)
%WELCHANOVA Welch ANOVA Test for Unequal Variances.
%The ANOVA F-test to compare the means of k normally distributed
%  populations is not applicable when the variances are unknown, and not
%  known to be equal. A spacial case, k=2, is the famous Behrens-Fisher 
%  problem (Behrens, 1929; Fisher, 1935). Welch (1951) test was proposed to
%  fill this void, a generalization to his 1947 previous paper (Welch, 
%  1947). 
%     
%  The Welch test for general k compares the statistic
%
%                  __
%                 \   
%                 /__ w_i*(m_i - M)^2/(k - 1)
%           FW = -----------------------------
%                      1 + 2/3*(k - 2)*L
%
%  to the F_[(k - 1),1/L] distribution. Where:
%                      __                 __             
%                     \                  \            
%  w_i = n_i/v_i; M = /__ w_i*m_i/W; W = /__ w_i; f_i = n_i - 1;
%          __
%         \
%      3* /__(1 - w_i/W)^2/f_i
%  L = ------------------------
%             (k^2 - 1)
%
%  [m_i = sample mean; v_i = sample variance; n_i = sample size]                                                     
%
%
%  Syntax: function welchanova(x,alpha)
%
%  Inputs:
%       x  data nx2 matrix (Col 1 = data; Col 2 = sample code)
%   alpha - significance level (default=0.05)
%
%  Outputs:
%       - Summary statistics from the samples
%       - Decision on the null-hypothesis tested
%
%  Taking the numerical example given at Research and Faculty Development,
%  University of Oregon, in his internet site [http://rfd.uoregon.edu/
%  files/rfd/StatisticalResources/glm10_homog_var.txt]. 
%  2. Testing Means with the Unequal Variance Model. For the testing for 
%  equality of variances for an ANOVA, first consider the following small
%  dataset which are measurements collected from 24 subjects randomly 
%  divided into one of three groups. Group 1 receives the standard 
%  treatment, group 2 receives a treatment similar to the standard, and 
%  group 3 gets a new and different treatment. The objective is to 
%  determine if the new treatment from group 3 is better than the 
%  treatments received by groups 1 and 2.
%
%  Data are:
%
%                  -----------------------------
%                              Group
%                  -----------------------------
%                     1          2          3 
%                  -----------------------------
%                   3.791      3.122      4.761
%                   5.174      2.514      3.884
%                   3.019      3.850      6.840
%                   3.218      2.564      9.150
%                   3.079      4.486      2.776
%                   4.054      3.199      5.398
%                   3.131      3.406      6.405
%                   2.822      3.986      2.115
%                  -----------------------------
%
%  Input data:
%
%  x=[3.791 1;5.174 1;3.019 1;3.218 1;3.079 1;4.054 1;3.131 1;2.822 1;
%  3.122 2;2.514 2;3.850 2;2.564 2;4.486 2;3.199 2;3.406 2;3.986 2;
%  4.761 3;3.884 3;6.840 3;9.150 3;2.776 3;5.398 3;6.405 3;2.115 3];
%
%  Calling on Matlab the function: 
%               welchanova(x,0.05)
%
%  Answer is:
%
%  Summary statistics from the samples.
%  --------------------------------------------------
%  Sample       Size        Mean           Variance 
%  --------------------------------------------------
%     1           8         3.5360           0.6096
%     2           8         3.3909           0.4752
%     3           8         5.1661           5.2988
%  --------------------------------------------------
% 
%  Welch's Analysis of Variance Table.
%  ----------------------------------------
%  SOV       df              F       P
%  ----------------------------------------
%  Treat.    2             2.075   0.1659
%  Error    12.7641
%  ----------------------------------------
%  The associated probability for the Welch's F test is equal or larger 
%  than 0.05. So, the assumption of sample means are equal was met.
%
%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%            Facultad de Ciencias Marinas
%            Universidad Autonoma de Baja California
%            Apdo. Postal 453
%            Ensenada, Baja California
%            Mexico.
%            atrujo@uabc.edu.mx
%
%  Copyright (C) May 15, 2012. 
%
%  --We thank Abdullah Chisti, University of Saskatchewan, Bangladesh, for
%   encourage us to produce this m-file.--
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2012). welchanova: Welch 
%     ANOVA Test for Unequal Variances. [WWW document]. URL http://
%     www.mathworks.com/matlabcentral/fileexchange/37121-welchanova
%
%  References:
%  Behrens, W. V. (1929), Ein beitrag zur Fehlerberechnung
%             bei wenigen Beobachtungen. (transl: A contribution to error
%             estimation with few observations). Landwirtschaftliche 
%             Jahrbcher, 68:80737.
%  Fisher, R. A. (1935), The fiducial argument in statistical inference.
%             Annals of Eugenics, 8:391398.
%  Welch, B. L. (1947), The generalization of Student's problem when 
%             several different population variances are involved. 
%             Biometrika, 34(12):2835  
%  Welch, B. L. (1951), On the comparision of several mean values: an
%             alternative approach. Biometrika, 38:330-336.
%

if nargin < 2 || isempty(alpha)
    alpha = 0.05; %default
elseif numel(alpha) ~= 1 || alpha <= 0 || alpha >= 1
    error('welchanova:BadAlpha','ALPHA must be a scalar between 0 and 1.');
end

X = x;
c = size(X,2);
if c ~= 2
    error('stats:welchanova:BadData','X must have two colums.');
end

%Remove NaN values, if any
X = X(~any(isnan(X),2),:);

k = max(X(:,2));

indice = X(:,2);
for i = 1:k
    Xe = indice == i;
    d(i).X = X(Xe,1);
    d(i).m = mean(d(i).X);
    d(i).v = var(d(i).X);
    d(i).n = length(d(i).X);
    d(i).f = d(i).n - 1;
    d(i).W = d(i).n/d(i).v;
    d(i).N = d(i).W*d(i).m;
end
m=cat(1,d.m);v=cat(1,d.v);n=cat(1,d.n);f=cat(1,d.f);W=cat(1,d.W);
N=cat(1,d.N);

M = sum(N)/sum(W);

for i = 1:k
    Xe = indice == i;
    d(i).A = ((1 - d(i).W/sum(W))^2)/d(i).f;
    d(i).B = d(i).W*(d(i).m - M)^2;
end
A=cat(1,d.A);B=cat(1,d.B);

L = 3*sum(A)/(k^2 - 1);

FW = (sum(B)/(k - 1))/(1 + 2/3*(k - 2)*L);  %Welch's F-statistic

v1 = k-1;  %numerator degrees of freedom
v2 = 1/L;  %denominator degrees of freedom

disp(' ')
disp('Summary statistics from the samples.')
disp('--------------------------------------------------')
disp(' Sample       Size        Mean           Variance ')
disp('--------------------------------------------------')
for i = 1:k
   fprintf('   %d           %i        %7.4f          %7.4f\n',i,n(i),m(i),v(i))
end
disp('--------------------------------------------------')
disp(' ')

P = 1-fcdf(FW,v1,v2);  %P-value

disp(' ')
disp('Welch''s Analysis of Variance Table.')
fprintf('----------------------------------------\n');
disp('SOV       df              F       P')
fprintf('----------------------------------------\n');
fprintf('Treat.  %3i%18.3f%9.4f\n\n',v1,FW,P);
fprintf('Error%11.4f\n',v2);
fprintf('----------------------------------------\n');

if P >= alpha;
    fprintf('The associated probability for the Welch''s F test is equal or larger than% 3.2f\n', alpha);
    fprintf('So, the assumption of sample means are equal was met.\n');
else
    fprintf('The associated probability for the Welch''s F test is smaller than% 3.2f\n', alpha);
    fprintf('So, the assumption of sample means are equal was not met.\n');
end

return,

Contact us