function [stats] = mcnemart(matr, alpha)
%McNemar Test
% Use
% [stats] = mcnemart(matr, alpha)
% Input:    matr is 2x2 table matrix [A B; C D].
%                alpha - significance level, 1-alpha confidence level
% Output stats: 
% stats.chi2= chi^2 statistic;
% stats.pval = p-value; 
% stats.h=0/1 (null/alternative),
% the rest of stats  structure consists of estimators of p1, p2, RD, RR, and OR
% and their confidence intervals.

if  nargin==1, alpha=0.05 ; end   

A=matr(1,1); B= matr(1,2); C= matr(2,1); D=matr(2,2);
% If any entry is 0, add 0.5 to all cells
if( A*B*C*D==0 )
    matr = matr + 0.5;
end
N = A + B + C + D;   stats.N = N;
%TEST

stats.chi2 = (B-C)^2/(B+C);  %mc nemar chi2
stats.pval = 1- chi2cdf(stats.chi2, 1);   
stats.h = stats.chi2 > chi2inv(1-alpha, 1); %H0 rejected?

p1 = (A + B)/N; stats.p1=p1;% row prob   
p2 = (A + C)/N; stats.p2=p2;% column prob

p11 = A/N;
p12 = B/N; 
p21 = C/N; 
p22 = D/N;

% risk difference 
stats.rd = p1 - p2;
   delta = p11*p22 - p12*p21;
s = sqrt( p1*(1-p1)/N + p2*(1-p2)/N + 2*delta/N );
ssel = sqrt( ((B+C)*N  - (B-C)^2)/N^3 );
   lbd =  p1 - p2 - norminv(1-alpha/2)*s;
   ubd = p1 - p2 + norminv(1-alpha/2)*s;
stats.rdint = [lbd, ubd]; 
% risk ratio
stats.rr = (A+B)/(A+C);
varlrr = 1/(A+B) + 1/(A+C) - 2*A/((A+B)*(A+C));
    lblrr = log((A+B)/(A+C)) - norminv(1-alpha/2)*sqrt(varlrr) ;
   ublrr = log((A+B)/(A+C)) + norminv(1-alpha/2)*sqrt(varlrr) ;
%stats.rrlint = [lblrr, ublrr];
stats.rrint = [exp(lblrr), exp(ublrr)];
% odds ratio
stats.or = B/C;
% miettinen approx
s2lor1 = (log(B/C))^2 * (B+C)/((abs(B-C)-1)^2 + 0.0000001);
slor1 = sqrt(s2lor1);
    lblor1 = log(B/C) - norminv(1-alpha/2)*slor1 ;
    ublor1 = log(B/C) + norminv(1-alpha/2)*slor1 ;
%stats.miettlint = [lblor1, ublor1];
stats.ormintmiet  = [exp(lblor1), exp(ublor1)];

% RGB approx
slor2 = sqrt(1/B + 1/C);
lblor2 = log(B/C) - norminv(1-alpha/2)*slor2 ;
ublor2 = log(B/C) + norminv(1-alpha/2)*slor2 ;
rgblint = [lblor2, ublor2];
stats.orrgb  = [exp(lblor2), exp(ublor2)];
% --------- exact method(Liddell) ------------------
 df1l = 2*C +2; df2l = 2*B;
 F1 = finv(1-alpha/2, df1l, df2l);
 tl = B/(B+(C+1)*F1 );
 orexall = tl/(1-tl);
 
 df1u = 2*B +2; df2u = 2*C;
 F2 = finv(1-alpha/2, df1u, df2u);
 tu = (B+1)*F2/(C+(B+1)*F2 );
 orexalu = tu/(1-tu);
 stats.orliddell  = [orexall, orexalu];

% %--------- approximate prob ---------
% phat = B/(B+C);
% pl = phat - norminv(1-alpha/2)*sqrt(phat * (1-phat)/(B+C)) - 1/(2*B+2*C);
% pu = phat + norminv(1-alpha/2)*sqrt(phat * (1-phat)/(B+C)) + 1/(2*B+2*C);
% stats.intor = [pl/(1-pl), pu/(1-pu)];

