Barnard’s unconditional exact test on 2×2 matrix

Barnard’s unconditional exact test on 2×2 matrix

There are two fundamentally different exact tests for comparing the equality of two binomial probabilities – Fisher’s exact test (Fisher, 1925), and Barnard’s exact test (Barnard, 1945). Fisher’s exact test (Fisher, 1925) is the more popular of the two. In fact, Fisher was bitterly critical of Barnard’s proposal for esoteric reasons that we will not go into here. For 2 × 2 tables, Barnard’s test is more powerful than Fisher’s, as Barnard noted in his 1945 paper, much to Fisher’s chagrin. Anyway, perhaps due to its computational difficulty the Barnard’s is not widely used.
Consider a general 2×2 matrix

Suppose that the common probability of infection under the null hypothesis is πA = πC = π. Then the probability of observing any generic table X is a product of two binomials:

The exact p-value is then the sum of such probabilities for all tables, X, that could have been observed which are at least as extreme as the observed table, X₀ , under the null hypothesis. Specifically, suppose T (X) is a “discrepancy measure” or test statistic that measures how discrepant any table X is relative to the type of table one would expect under the null hypothesis. Then, for any given π, the exact p-value of the observed table X₀ is:

Now there exist many different test statistics for determining if a table X is more or less extreme than the observed table X₀ under the null hypothesis: we will use the popular Wald statistic:

Although the Wald test statistic T (X) is well defined, equation (2) is not of much practical use for computing an exact p-value because it depends on knowing the value of π, the common probability under the null hypothesis. One could of course use the data to estimate π and then substitute this estimate into equation (2). But then the p-value would no longer be exact. The main difference between the Fisher’s and Barnard’s tests is the manner in which they eliminate this nuisance parameter from (2) without sacrificing the exactness.

Fisher’s exact test

Fisher’s exact test gets rid of π by restricting attention only to generic 2 × 2 tables X in which A + B = Rs₁, the sum of responses that were actually observed in the data. Accordingly define:

Fisher noted that, under the null hypothesis, provided we confine our attention only to tables X ∈ Γ(Rs₁), each table has a hypergeometric probability that no longer depends on the nuisance parameter π (see my post: Fisher’s conditional exact test on 2×2 matrix).

Barnard’s exact test

Barnard’s exact test is an unconditional test. It generates the exact distribution of T(X) by considering all the tables X ∈ Γ, and rather than just the tables X ∈ Γ(Rs₁). What shall we do with the unknown nuisance parameter π in the above p-value (equation 1)? Barnard suggested that we calculate p(π) for all possible values of π ∈ (0, 1) and choose the value, π^∗ , say, that maximizes p(π). And this is computationally very very difficult….

The algorithm

Now, let we see the algorithm in deep.

The first step is to compute the Wald’s statistic for each possible value of A and B (from now on I and J).

I∈[0:1:Cs₁] and J∈[0:1:Cs₂]. One could write:

1
2
3
4
5
6

Cs=sum(x); N=sum(Cs) %Columns sum and total observed elements
for I=0:1:Cs(1)
    for J=0:1:Cs(2)
         TX(I+1,J+1)=(I/Cs(1)-J/Cs(2))/sqrt(((I+J)/N)*(1-(I+J)/N)*sum(1./Cs))
    end
end

Notice that for (I=0 & J=0) or (I=Cs(1) & J=Cs(2)) T(I,J)=0/0=NaN. So there are two points where T(I,J) is not defined and NaN should be replaced with 0.
Two nested for…end loops are time consuming and are useless considering MatLab functions…

A more efficient code should be:

1
2
3
4
5
6

Cs=sum(x); N=sum(Cs) %Columns sum and total observed elements
[I,J]=ndgrid(0:1:Cs(1) 0:1:Cs(2));
warning('OFF','MATLAB:divideByZero') %disable the warning
TX=(I./Cs(1)-J./Cs(2))./realsqrt(((I+J)./N).*(1-((I+J)./N)).*sum(1./Cs)); %computes the Wald's statistics
warning('ON','MATLAB:divideByZero') %renable the warning
TX([1 end])=0; %resolve the singularities

The ndgrid function creates two square matrix:

So TX will be a Cs(1)+1 x Cs(2)+1 matrix. Using .* and ./ all elements of the matrix will be correctly multiplied or divided.

To set up the observed matrix (TX_o) and the index of all TX≥TX_o

1
2

TX0=abs(TX(x(1)+1,x(3)+1)); %catch the observed Table (TXo), taking the 0 in account
idx=TX>=TX0; %set the index of all T(X)>=TXo

Now, we must compute P(X|π) for each I, J and π∈(0, 1) . I and J were just defined by ndgrid function. Expanding the binomials we obtain:

As I explained in my previous post Fisher’s conditional exact test on 2×2 matrix, it is useful to use the logarithms of P(X|π): the factorials can be obtained using the MatLab Gammaln function; the multiplications became sums; the divisions became subtractions and the power became multiplications.

CF is a matrix with the same dimensions of TX because it is a linear combination of I and J.

The S matrix is a 3D matrix: two dimensions are Cs₁+1 and Cs₂+1; the third dimension is set by the user. In fact in the arguments of this function, the user can choose to set the length (the Tbx variable) of the vector π (by default it is 100). Let we see the code:

1
2
3
4
5
6
7
8
9
10

%Set the basic parameters...
A=[1 1 Tbx];
B=Cs+1;
npa = linspace(0.0001,0.9999,Tbx);
LP=log(npa);
ALP=log(1-npa);
E=repmat(I+J,A); F=N-E;
%Generate a 3D matrix (Cs1+1 x Cs2+1 x Tbx): this is a "box" where each of Tbx 2D
%slice is the Cs1+1 x Cs2+1 matrix CF
CF=repmat(sum(gammaln(B))-(gammaln(I+1)+gammaln(J+1)+gammaln(B(1)-I)+gammaln(B(2)-J)),A);

Now we have two 1xTbx arrays (LP and ALP) and three Cs₁+1 x Cs₂+1 x Tbx “boxes” (E, F and CF).
As Peter J. Acklam wrote in his “Matlab array manipulation tips and tricks” , there is a no for…loop solution to multiply each 2D slice with the corresponding
element of a vector. He wrote:
Assume X is an m-by-n-by-p array and v is a row vector with length p.
How does one write:

Y = zeros(m, n, p);
for i = 1:p
Y(:,:,i) = X(:,:,i) * v(i);
end

with no for-loop? One way is to use:
Y = X .* repmat(reshape(v, [1 1 p]), [m n]);

So we can construct the two others 3D boxes:
Box1=CF
Box2=E.*repmat(reshape(LP,[1 1 Tbx]),[Cs1+1 Cs2+1])=E.*repmat(reshape(LP,A),B)
Box3=F.*repmat(reshape(ALP,[1 1 Tbx]),[Cs1+1 Cs2+1])=F.*repmat(reshape(ALP,A),B)
Finally, we can obtain the S box that is the sum of these three boxes (remember that we are using logarithms, so we must convert them using the exp function)

1

S=exp(CF+E.*repmat(reshape(LP,A),B)+F.*repmat(reshape(ALP,A),B));

Now we are at the last point: for each 2D slice of the S box we must sum the values corresponding to TX>=TXo, obtaining a new vector P.
To use the logical indexing tecnique, we must replicate the idx matrix:
S(repmat(idx,A)) is a column vector that has L*Tbx row;
L=sum(idx(idx==1)) is the number of 1 in the idx matrix.
Using the reshape function we can obtain a new LxTbx matrix: in each column we’ll have the P(X|π & TX>=TXo) and using the function sum we’ll, finally,
obtain the P vector.

1

P=sum(reshape(S(repmat(idx, A)),sum(idx(idx==1)),Tbx));

Now we can find the p-value and the nuisance parameter:

1
2

PV=max(P); %The p-value is tha max value of the P vector;
np=npa(P==PV); %The nuisance parameter is the value of the npa array coinciding with PMax

Display the results…

An example

Supposing to have this matrix

Using the Fisher’s conditional exact test the p-value (1-tailed) is 0.6407: you will conclude that the vaccine is not effective.
Using the Barnard’s unconditional exact test the results will be:

2×2 matrix Barnard’s exact test: 100 16×16 tables were evaluated
——————————————————————————–
Wald statistic = 1.8943 – Nuisance parameter = 0.6666
p-values 1-tailed = 0.034074 2-tailed = 0.068148
——————————————————————————–

So with the Barnard’s test you will conclude that the vaccine is effective.
Anyway, using the Fisher’s conditional exact test with the Lancaster’s correction (mid p) you will obtain a 2-tailed p-value=0.04774 ….

If any problems occurs in execution, or if you found a bug, have a suggestion or question just contact me at: giuseppe dot cardillo-edta at poste dot it

You can visit my homepage http://home.tele2.it/cardillo

My profile on XING http://www.xing.com/go/invita/13675097

My profile on LinkedIN http://it.linkedin.com/in/giuseppecardillo

Download Now

Popularity: 1% [?]

Fisher’s conditional exact test on 2×2 matrix

Fisher’s conditional exact test on 2×2 matrix

Fisher’s exact test is a statistical significance test used in the analysis of contingency tables where sample sizes are small. It is named after its inventor, R. A. Fisher, and is one of a class of exact tests, so called because the significance of the deviation from a null hypothesis can be calculated exactly, rather than relying on an approximation that becomes exact in the limit as the sample size grows to infinity, as with many statistical tests. Fisher is said to have devised the test following a comment from Muriel Bristol, who claimed to be able to detect whether the tea or the milk was added first to her cup.

The test is useful for categorical data that result from classifying objects in two different ways; it is used to examine the significance of the association (contingency) between the two kinds of classification. So in Fisher’s original example, one criterion of classification could be whether milk or tea was put in the cup first; the other could be whether Ms Bristol thinks that the milk or tea was put in first. We want to know whether these two classifications are associated – that is, whether Ms Bristol really can tell whether milk or tea was poured in first. Most uses of the Fisher test involve, like this example, a 2 x 2 contingency table. The p-value from the test is computed as if the margins of the table are fixed, i.e. as if, in the tea-tasting example, Ms. Bristol knows the number of cups with each treatment (milk or tea first) and will therefore provide guesses with the correct number in each category. As pointed out by Fisher, this leads under a null hypothesis of independence to a hypergeometric distribution of the numbers in the cells of the table.

Consider a general 2×2 matrix

What does the Fisher’s algorithm plan?

1. Enumerate all possible tables, given RS1 RS2 CS1 CS2;
2. For each table compute p-value using hypergeometric distribution;
3. Sum all p-values≤p-value of the observed table

Now we’ll analyze the algorithm step-by-step, observing the optimizations that can be implemented using MatLab.

A 2×2 matrix has only one degree of freedom. What does this mean? This means that if you change the cell A all the other cells are determined: B=RS1-A; C=CS1-A; D=RS2-CS1+A.

Question: Which is the range of variation for A?
Answer: 0≤A≤min(RS1,CS1)

First optimization: rearrange the matrix (if needed) so as to obtain: RS1≤RS2 and CS1≤CS2

1
2
3
4
5
6
7
8
9
10
11
12
13

Rs=sum(x,2); %rows sum
Cs=sum(x); %columns sum
N=sum(Rs); %sum of all elements
 
%Rearrange the matrix if necessary.
if ~issorted(Rs)
    x=flipud(x);
    Rs=sort(Rs);
end
if ~issorted(Cs)
    x=fliplr(x);
    Cs=sort(Cs);
end

Now you can enumerate all the possible tables:

1
2

A=0:1:min(Rs(1),Cs(1)); %all possible values for the first cell
z=[A;Rs(2)-Cs(1)+A;Rs(1)-A;Cs(1)-A;]; %all possible tables

Note that in the z matrix the cells are disposed in this order [A D B C] (in the next step you will understand the reason…).

Step 2: For each table compute the p-value using the hypergeometric distribution

First of all, now we know how many tables we wrote down and so we can use the preallocation (see http://www.advancedmcode.org/writing-fast-matlab-code-3-preallocation.html).

1

np=zeros(1,length(A)); %p-values vector preallocation

In this step we must compute the p-value for each table using the hypergeometric distribution. In general, the p-value is computed as:

Considering that we are using a 2×2 matrix:

Remember that: and x is a natural number.

As you can see, the factorial (!) growth very quickly. Moreover the computation of a factorial is very time expensive. Luckly, there is a very important relation between the factorial and another function: the Gamma function.

Matlab can compute the Gamma function using a matrix as input. Anyway, we didn’t yet solve the factorial growing (20!=2.4329e+18). The management of so big number is very troubling. The solution is to use the logarithm function: in fact, log(20!)=42.3356 a more tractable number.

Matlab, of course, can compute the gammaln function. Now, we are computing log(p-value), so when we’ll finish we must convert into p-value using exp function.

Remember two log properties (they will be useful later): LOG(A*B)=LOG(A)+LOG(B)   and   LOG(A/B)=LOG(A)-LOG(B).

Now, it is the time to observe another useful property of 2×2 matrix.
When you go from the Table T_i to the Table T_i+1:

Substituting you will find a recursive relation:

And, more in general:

At this point, the only problem is to compute the p-value of the first table. Consider the first Table:

And so:

Finally, the code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

np=zeros(1,length(A)); %p-values vector preallocation
lz=log(z);
%LOG(X!)=GAMMALN(x+1)
np(1)=sum(gammaln([Rs(2)+1 Cs(2)+1])-gammaln([N+1 z(2)+1]));
%remember that
%np(i+1)=np(i)*[B(i)*C(i)]/(A(i+1)*D(i+1)]=np(i)*f(i)
%This formula is vectorizable and is recursive!
%using Logarithm to improve computation
%log(np(i+1))=log(np(i))+[log(B(i))+log(C(i))]-[log(A(i+1))+log(D(i+1))]=log(np(i))+f(i)
%J>1
%log(np(J))=log(np(1))+sum(f(1:J-1))=log(np(1))+cumsum(f)
f=sum(lz(3:4,1:end-1))-sum(lz(1:2,2:end));
np(2:end)=np(1)+cumsum(f);
np=exp(np);

Now, to obtain the 2-tail p-value we must sum all p-values≤p-value of the observed table. In the z matrix, the observed table will be in the A+1 row (because you must also consider the 0) and, of course, in the np vector the p-value of the observed table will be in A+1 row.
On the contrary, if you desire the 1-tail p-value you must sum all A+1 p-values.
The code:

1
2
3
4
5
6

W=x(1)+1;
if x(1)<=round(A(end)/2) %choose direction
     P=[sum(np(1:W)) sum(np(W:end)) sum(np(np<=np(W)))];
else
     P=[sum(np(W:end)) sum(np(1:W)) sum(np(np<=np(W)))];
end

Finally, there is another “little” question. The Fisher’s exact test is a “conditional” method and is too conservative in comparison to an “unconditional” method like the Barnard’s test; in other words is less powerful than Barnard’s test. On the contrary, the Barnard’s test is so difficult to implement that none uses it (but me… see mi post on Barnard’s unconditional exact test).
The Lancaster’s correction (mid-p) to Fisher’s exact test makes it less conservative and more powerful and the results are similar to Barnard’s Test.

The Lancaster’s correction envisages that the p-values of all tables more extreme than the observed are summed to one half of the observed p-value. In code

1

P-mid=0.5*np(W)+sum(np(np<np(W)));

Don’t forget to display the results!

The routine computes the Power and, if necessary, the sample sizes needed to achieve a power=0.80 using a modified asymptotic normal method with continuity correction as described by Hardeo Sahai and Anwer Khurshid in Statistics in Medicine, 1996, Vol. 15, Issue 1: 1-21. Because there aren’t interesting tricks, if you are interested you can directly read the code.

If any problems occurs in execution, or if you found a bug, have a suggestion or question just contact me at: giuseppe dot cardillo-edta at poste dot it

You can visit my homepage http://home.tele2.it/cardillo

My profile on XING http://www.xing.com/go/invita/13675097

My profile on LinkedIN http://it.linkedin.com/in/giuseppecardillo

Download Now

Popularity: 1% [?]

Coins

Coins

How do you choose the coins in your pocket to pay something? There is an algorithm to use the max numbers of coins? Of course there is, as shown in the book “L’algoritmo del parcheggio” (The parking algorithm) by Prof. Furio Honsell. This GUI help you to choose the exact combination of coins that you have in your pocket to pay the bill.

For example: Suppose you are in Naples (Italy, my city) and have in your pocket this situation. Now you are in “Piazza del Plebiscito” and you want to drink a coffee in the famous “Gambrinus” Bar. You go to the cash, order the coffee and the beautiful cashier give you the sales cash asking 0.80 euro (I wish…). All these coins in your pocket are very heavy to carry, and so you want to give as many as possible coins to the cashier (that will be happiest to have all these coins in the cash to give them as change). The situation can be resumed into a mathematical model:

The bill to pay is stored into the T variable. As you can see, it is possible to us all coins less than 0.20€ to pay the bill. But what is the most efficient methods?

1
2
3
4
5
6
7
8

while T>0 %saddle the cashier with coins
    V0=find(V>T,1,'first'); %find the first value in the V array greater than T
    coinsout(V0)=coinsout(V0)+1; %give the coin to the cashier
    coinsin(V0)=coinsin(V0)-1;%this coin is no more in your pocket
    coinval(V0)=coinval(V0)-values(V0); %update the array
    V=cumsum(coinval); %update the array
    T=T-values(V0); %update the bill to pay
end

But is it truly so simple? No, of course. The floating point representation of cents sometimes creates troubles: for example: 1.0000-1 could not be 0 but 1.45782e-17 (that is >0…. ). To fix this problem, multiply (and after divide before printing results) by 100: values, coinval and T.

Finally, this is the result:

If any problems occurs in execution, or if you found a bug, have a suggestion or question just contact me at: giuseppe dot cardillo-edta at poste dot it

You can visit my homepage http://home.tele2.it/cardillo

My profile on XING http://www.xing.com/go/invita/13675097

My profile on LinkedIN http://it.linkedin.com/in/giuseppecardillo

Download Now

Popularity: 1% [?]

Next Page »

• Recent Posts

  • Gabor filters for Face Recognition
  • Free Open Source FEM Softwares
  • Facial Gender Recognition Using GA
  • Debugging, Profiling, & Style Guidelines in Matlab
  • Text-Independent Speaker Recognition Based on Neural Networks
  • Unstructured 2D Triangular Mesh Generation in Matlab
  • Neural Network Technology
  • Strings, Cells, Structs, and Sets in Matlab
  • High Definition Image Compression Technology
  • RandSphere: Generating Random Points On The Surface Of An nD Sphere

• Most rated

  • 7.1: MATLAB for Digita ... (13 votes)
  • 8.3: ROC curve (8 votes)
  • 9.9: Plotting in Matla ... (8 votes)
  • 6.6: MATPOWER: A MATLA ... (8 votes)
  • 10.0: Poker Predictor: ... (7 votes)
  • 10.0: Object Oriented P ... (6 votes)
  • 8.2: The Future of MAT ... (5 votes)
  • 8.2: Poker Predictor: ... (5 votes)
  • 8.8: ASK, OOK, FSK, BP ... (5 votes)
  • 8.2: Flow of Control a ... (5 votes)
  • 9.0: Quick Delaunay Tr ... (5 votes)
  • 10.0: GLTree:Fast K-Nea ... (5 votes)
  • 8.8: Surface Reconstru ... (5 votes)
  • 8.8: Writing Fast Matl ... (5 votes)
  • 10.0: Fast Surface Reco ... (4 votes)

• Popular Posts

  • Plotting in Mat... 97 view(s) per day
  • Object Oriented... 52 view(s) per day
  • MATLAB for Digi... 50 view(s) per day
  • Free Open Sourc... 47 view(s) per day
  • The Matlab Matr... 37 view(s) per day
  • Gabor filters f... 33 view(s) per day
  • How to create a... 33 view(s) per day
  • Surface Reconst... 30 view(s) per day
  • ASK, OOK, FSK, ... 28 view(s) per day
  • Interfacing MAT... 27 view(s) per day
  • Writing Fast Ma... 25 view(s) per day
  • Optical Charact... 25 view(s) per day
  • ROC curve 25 view(s) per day
  • Strings, Cells,... 24 view(s) per day
  • Matlab and Micr... 23 view(s) per day
  • Poker Predictor... 23 view(s) per day
  • Surface Reconst... 20 view(s) per day
  • Unstructured 2D... 20 view(s) per day
  • MATPOWER: A MAT... 19 view(s) per day
  • How to plot a c... 19 view(s) per day
  • Writing and Usi... 19 view(s) per day
  • Calling Java, R... 15 view(s) per day
  • Fast Convex Hul... 15 view(s) per day
  • Quick Delaunay ... 15 view(s) per day
  • Making, reading... 15 view(s) per day

• Tags

  2d 3D advanced matlab code ann cards Convex Hull copyright. dct delaunay Delaunay2.5D Delaunay2_5D editbox face matching Face recognition fast matlab code feature extraction FEM filter function GLTree Graphics GUI Example identification image Image processing iris manifold matching Matlab source code mesh MRI nearest neighbor nearest neighbour neural networks Poker poker odds popup pushbutton. recognition scanner 3D strain stress String surface surface reconstruction System Texas Holdem text verification watermark

• Category Cloud

  Media Logistics Java Fun Genetic MEX science Whatever Web Medicine Neural Network Linear Algebra Electronic Electric Physics Mechanics Sound technology C Economy Biotecnology Games Numerical Analysis FEM Geometry Probability CAD Optimization Mathematics Statistics Signal Processing NewsGroup Graphics Utility General Code Optimization GUI Algorithms Computational Geometry Image processing Tutorials

• Recent Comments

  • Li-aung Yip ("Lewis") on Debugging, Profiling, & Style Guidelines in Matlab
  • priyanka pawar on CONV2 Overlap-add Method
  • TB on Surface Reconstruction from Scattered Points Cloud: MyCrust Robust
  • Luigi Giaccari on Surface Reconstruction from Scattered Points Cloud: MyCrust Robust
  • TB on Surface Reconstruction from Scattered Points Cloud: MyCrust Robust
  • Admin on Free Open Source FEM Softwares
  • Thomas on Free Open Source FEM Softwares
  • Luigi Giaccari on Volume Enclosed by a Triangulated Surface
  • Christopher Clark on Volume Enclosed by a Triangulated Surface
  • Cam Straatsma on In-Polyhedron Test

• Search

• Translator

• Make a contribution

  

• Log-In/Register

  Start posting on Advanced M-code!

  
  

  Register

  
  

  Log In

• Follow Us

  Facebook

  Twitter

  By Email

  Entries RSS

  Comment RSS

• Blogroll

  • advancedsourcecode
  • Online Foreign exchange
  • Statistical Consulting
  • Undocumented Matlab
  • walkingrandomly

• Video Posts

  

• Archives

  • August 2010 (2)
  • July 2010 (7)
  • June 2010 (11)
  • May 2010 (17)
  • April 2010 (11)
  • March 2010 (13)
  • February 2010 (29)
  • January 2010 (21)
  • December 2009 (16)
  • November 2009 (18)
  • October 2009 (21)
  • September 2009 (12)
  • August 2009 (28)

• Pages

  • Contents
    • Books
    • Commercial Softwares
    • Papers
    • Posts Directory
  • Infos
    • FAQ
    • Writing Post (Tutorial)
  • NewsGroup
  • Submit

• StatPress

  Visits Today: 128
  Visits Last Month: 12981
  Visitors on line: 4
  

• Meta

  • Register
  • Log in
  • Entries RSS
  • Comments RSS
  • WordPress.org
• Free SEO Tools