data step1;

** This program is used to calculate one-sided lower confidence limit;

*** Warning ***;

** The approximate lower limit should be used cautiously if any of the cell frequency is less than 5;

** If the OR of drug B versus drug A is our interest;

** n1 represents the total number patients with discordant response in group I receiving A-B sequence;

** x1 represents the number of patients among n1 with a positive in period 2 and a negative in period 1;

** n2 represents the total number of patients with discordant response in group II receiving B-A sequence;

** x2 represents the number of patients among n2 with a positive in period 2 and a negative in period 1;

** If the OR of drug A versus drug B is our interest;

** n1 represents the total number patients with discordant response in group II receiving B-A sequence;

** x1 represents the number of patients among n1 with a positive in period 2 and a negative in period 1;

** n2 represents the total number of patients with discordant response in group I receiving A-B sequence;

** x2 represents the number of patients among n2 with a positive in period 2 and a negative in period 1;

alpha=0.05; ** one-sided;

za=1.645; ** one-sided;

n1=48;

x1=32;

n2=56;

x2=15;

xt=x1+x2;

a=max(xt-n2, 0);

b=min(xt, n1);

if (n1-x1) gt 0 and x1 gt 0 and (n2-x2) gt 0 and x2 gt 0 then do;

mleor=sqrt((x1*(n2-x2))/((n1-x1)*x2));

estdor=sqrt((1/x1+1/(n1-x1)+1/x2+1/(n2-x2))/4);

lowmle=mleor*exp(-1*za*estdor);

put 'the mle of the OR' mleor 60-70 3;

put 'the approx lower limit of the OR without adding 0.5' lowmle 60-70 3;

end;

mleq=(x1+0.5)*(n2-x2+0.5)/((x2+0.5)*(n1-x1+0.5));

adjmle=sqrt(mleq);

estd=sqrt((1/(x1+0.5)+1/(n1-x1+0.5)+1/(x2+0.5)+1/(n2-x2+0.5))/4);

apxlow=adjmle*exp(-1*za*estd);

put 'the MLE with adjustment of adding 0.5' adjmle 60-70 3;

put 'the approx lower limit of the OR with adding 0.5' apxlow 60-70 3;

p=mleq;

snum=0;

sden=0;

do x=a to b;

prob=exp(lgamma(n1+1)-lgamma(x+1)-lgamma(n1-x+1))*

exp(lgamma(n2+1)-lgamma(xt-x+1)-lgamma(n2-xt+x+1))*p**x;

if x ge x1 then snum+prob;

sden+prob;

end;

cx=snum/sden;

check1=cx-alpha;

llimu=mleq;

lliml=0;

diff=llimu-lliml;

do while (diff gt 0.001);

p=(llimu+lliml)/2;

snum=0;

sden=0;

do x=a to b;

prob=exp(lgamma(n1+1)-lgamma(x+1)-lgamma(n1-x+1))*

exp(lgamma(n2+1)-lgamma(xt-x+1)-lgamma(n2-xt+x+1))*p**x;

if x ge x1 then snum+prob;

sden+prob;

end;

cx=snum/sden;

check5=cx-alpha;

if check5*check1 gt 0 then llimu=p;

if check5*check1 le 0 then lliml=p;

diff=llimu-lliml;

end;

lowlimit=sqrt(p);

put 'the exact conditional lower limit of the OR' lowlimit 60-70 3;