function Barnardextest(a,b,c,d)
%BARNARDEXTEST Barnard's Exact Probability Test.
% BARNARDEXTEST, as the Fisher's exact test, performs the exact probability
% test for a table of frequency data cross-classified according to two categorical
% variables, each of which has two levels or subcategories (2x2). It is a
% non-parametric statistical test used to determine if there are nonrandom
% associations between the two categorical variables. Barnard's exact test is used
% to calculate an exact P-value with small number of expected frequencies, for
% which the Chi-square test is not appropriate (in case the total number of
% observations is less than 20 or the number of frequency cells are less than
% 5). The test was proposed by G. A. Barnard in two papers (1945 and 1947).
% While Barnard's test seems like a natural test to consider, it's not at all
% commonly used. This probably due that it is a little unknown.
%
% According to the next 2x2 table design,
%
% Var.1
% --------------
% a b r1=a+b
% Var.2
% c d r2=c+d
% --------------
% c1=a+c c2=b+d n=c1+c2
%
% The Barnard's exact test is a unconditioned test for it generates the exact
% distribution of the Wald statistic T(X),
%
% T(X) = abs((p(a) - p(b))/sqrt(p*(1-p)*((1/c1)+(1/c2)))),
% where,
% p(a) = a/c1, p(b) = b/c2 and p = (a+b)/n,
%
% by considering all tables X H and calculates P(PI) for all possible values of
% PI in (0, 1). The P-value is defined as the maximum value of p(PI),
%
% p(PI) = SUM_T(X)>=T(Xo) C(c1,a)*C(c2,b)*PI^(a+b)*(1-PI)^(n-a-b).
%
% Where T(Xo) is the actual (observed) Wald statistic.
%
% We calculate P(PI) for all possible values of PI in (0, 1) and choose the
% value that maximizes P(PI).
%
% So, the Barnard's exact P-value is evaluated as,
%
% P = sup{P(PI): PI (0, 1)}.
%
% Perhaps due to its computational difficulty it is not widely used until recently,
% where the computers make it feasible. It is considering that the Barnard's exact
% test is more powerful than the Fisher's one as a direct consequence of restricting
% the Fisher's test to 2x2 tables belonging to the conditional reference H(a+b),
% rather than the larger unconditional reference of Barnard's set for all H's.
% Consequently, the number of distinct P-values obtained with the Fisher's exact
% test is less than the corresponding number that one could get with the Barnard's exact
% test. So, for a restrict alpha-error, the Fisher's procedure will usually be more
% conservative, resulting in a loss of power.
%
% Syntax: function Barnardextest(a,b,c,d)
%
% Inputs:
% a,b,c,d - observed frequency cells
%
% Output:
% A table with:
% - Wald statistic, Nuisance parameter and P-value
% - Plot of the nuisance parameter PI against the corresponding P-value for
% all the PI in (0, 1). It shows the maximized PI where it attains the
% P-value. For a least conflict with figure, press the left mouse button on
% the legend and drag to the desired location.
%
% Example: From the example given on the Ina Parks S. Howell's internet homepage
% (http://www.fiu.edu/~howellip/Fisher.pdf). Suppose Crane and Egret are
% two very small collages, the results of the beginning physics course at
% each of the two schools are given in the follow table.
%
% Physics
% Pass Fail
% ---------------------
% Crane 8 14
% Collage
% Egret 1 3
% ---------------------
%
% Calling on Matlab the function:
% Barnardextest(8,14,1,3)
%
% Answer is:
%
% Table of the Barnard's exact test results for:
% a = 8, b = 14, c = 1, d = 3
% -------------------------------------------------
% Wald statistic Nuisance parameter P
% -------------------------------------------------
% 0.4394 0.9001 0.4159
% -------------------------------------------------
%
% [Note:For a least conflict with figure, press the left mouse button on the legend and drag to the desired location.]
%
% Created by A. Trujillo-Ortiz, R. Hernandez-Walls and A. Castro-Perez
% Facultad de Ciencias Marinas
% Universidad Autonoma de Baja California
% Apdo. Postal 453
% Ensenada, Baja California
% Mexico.
% atrujo@uabc.mx
% And the special collaboration of the post-graduate students of the 2004:2
% Multivariate Statistics Course: Laura Rodriguez-Cardozo, Norma Alicia Ramos-Delgado,
% and Rene Garcia-Sanchez.
% Copyright (C) October 25, 2004.
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A., R. Hernandez-Walls, A. Castro-Perez, L. Rodriguez-Cardozo
% N.A. Ramos-Delgado and R. Garcia-Sanchez. (2004). Barnardextest:Barnard's Exact
% Probability Test. A MATLAB file. [WWW document]. URL http://www.mathworks.com/
% matlabcentral/fileexchange/loadFile.do?objectId=6198
%
% References:
%
% Barnard, G.A. (1945), A new test for 2x2 tables. Nature, 156:177.
% Barnard, G.A. (1947), Significance tests for 2x2 tables. Biometrika, 34:123-138.
% Howell, I.P.S. (Internet homepage), http://www.fiu.edu/~howellip/Fisher.pdf
% Mehta, C.R. and Hilton, J.F. (1993), Exact power of conditional and unconditional
% tests: going beyond the 2x2 contingency table. American Statistician,
% 47:91-98.
%
if nargin < 4,
error('Too few arguments.');
end;
c1 = a+c;
c2 = b+d;
n = c1+c2;
pao = a/c1;
pbo = b/c2;
pxo = (a+b)/n;
warning off
TXO = abs(pao-pbo)/sqrt(pxo*(1-pxo)*((1/c1)+(1/c2)));
warning on
nans = isnan(TXO);
mo = find(nans);
TXO(mo) = zeros(size(mo));
n1 = prod(1:c1);
n2 = prod(1:c2);
P = [];
for p=.0001:.01:1
S = [];TX = [];
for i = 0:c1
for j = 0:c2
fac1 = prod(1:i);
fac2 = prod(1:j);
fac3 = prod(1:(c1-i));
fac4 = prod(1:(c2-j));
s = (n1*n2*(p.^(i+j))*(1-p).^(n-(i+j)))./(fac1*fac2*fac3*fac4);
S = [S s];
pa = i/c1;
pb = j/c2;
px = (i+j)/n;
warning off
tx = (pa-pb)/sqrt(px*(1-px)*((1/c1)+(1/c2)));
warning on
TX = [TX tx];
nans = isnan(TX);
m = find(nans);
TX(m) = zeros(size(m));
end;
end;
P = [P;sum(S(find(TX >= TXO)))];
end;
p = 0.0001:.01:1;
plot(p,P,'b-')
np = find(P==max(P));
[p(np) P(np)];
disp(' ')
disp('Table of the Barnard''s exact test results for:')
disp(['a = ' num2str(a) ', b = ' num2str(b) ', c = ' num2str(c) ', d = ' num2str(d) '']);
fprintf('---------------------------------------------------\n');
disp(' Wald statistic Nuisance parameter P ');
fprintf('---------------------------------------------------\n');
fprintf(' %10.4f %10.4f %10.4f\n',[TXO,p(np),P(np)].');
fprintf('---------------------------------------------------\n');
disp(' ')
disp('[Note:For a least conflict with figure, press the left mouse button on the legend and drag to the desired location.]')
hold on
plot(p(np),P(np),'ro')
title('Barnard''s exact P-value as a function of nuisance parameter PI.','FontSize',12);
xlabel('Nuisance parameter (PI)');
ylabel('P-value [P( TX >= TXO | PI )]');
legend(['P-value = ',num2str(P(np))],['Nuisance parameter (PI) = ',num2str(p(np))]);
return;
