data step1;
***** This program is for simulation (Table 1)in testing equality;
***** under an incomplete block crossover design with ordinal responses;
options ls=80;
alpha=0.05;
nsimul=10000; ** the number of repeated samples;
do sigma=0.5,1.0; ** this is the standard deviation of normal random effects;
do dda=0.0,0.50; ** this is the effect for treatment A versus placebo;
do delta=0.0,1.0; ** this is the effect for treatment B versus treatment A;
ddb=dda+delta; ** this is the effect for treatment B versus placebo;
do n =10,15,25; ** this is the number of subjects per group;
do period=0.10; ** this is the effect for period 2 versus period 1;
sumrej1=0; sumrej2=0; sumrej3=0;
sumrej4=0; sumrej5=0; sumrej6=0;
sumreja=0; sumrejb=0;
sumf1=0; sumf2=0; sumf3=0;
sumf4=0; sumf5=0; sumf6=0;
do num=1 to nsimul;
** yg10 and yg01 (where g = 1, 2, 3, 4, 5, 6) are the counters;
** of patients with response at period 1 larger than response at period 2;
** and of patients with response at period 1 less than response at period 2;
** in group g respectively;
y110=0; y101=0;
y210=0; y201=0;
y310=0; y301=0;
y410=0; y401=0;
y510=0; y501=0;
y610=0; y601=0;
do i=1 to n;
mu=rannor(1453161)*sigma;
p10=(exp(mu)/(1+exp(mu)))*(1/(1+exp(mu+dda+period)));
p01=(1/(1+exp(mu)))*(exp(mu+dda+period)/(1+exp(mu+dda+period)));
response=rantbl(1352121,p10,p01,1-p10-p01);
if response eq 1 then y110+1;
if response eq 2 then y101+1;
mu=rannor(1453161)*sigma;
p10=(exp(mu+dda)/(1+exp(mu+dda)))*(1/(1+exp(mu+period)));
p01=(1/(1+exp(mu+dda)))*(exp(mu+period)/(1+exp(mu+period)));
response=rantbl(1352121,p10,p01,1-p10-p01);
if response eq 1 then y210+1;
if response eq 2 then y201+1;
mu=rannor(1453161)*sigma;
p10=(exp(mu)/(1+exp(mu)))*(1/(1+exp(mu+ddb+period)));
p01=(1/(1+exp(mu)))*(exp(mu+ddb+period)/(1+exp(mu+ddb+period)));
response=rantbl(1352121,p10,p01,1-p10-p01);
if response eq 1 then y310+1;
if response eq 2 then y301+1;
mu=rannor(1453161)*sigma;
p10=(exp(mu+ddb)/(1+exp(mu+ddb)))*(1/(1+exp(mu+period)));
p01=(1/(1+exp(mu+ddb)))*(exp(mu+period)/(1+exp(mu+period)));
response=rantbl(1352121,p10,p01,1-p10-p01);
if response eq 1 then y410+1;
if response eq 2 then y401+1;
mu=rannor(1453161)*sigma;
p10=(exp(mu+dda)/(1+exp(mu+dda)))*(1/(1+exp(mu+ddb+period)));
p01=(1/(1+exp(mu+dda)))*(exp(mu+ddb+period)/(1+exp(mu+ddb+period)));
response=rantbl(1352121,p10,p01,1-p10-p01);
if response eq 1 then y510+1;
if response eq 2 then y501+1;
mu=rannor(1453161)*sigma;
p10=(exp(mu+ddb)/(1+exp(mu+ddb)))*(1/(1+exp(mu+dda+period)));
p01=(1/(1+exp(mu+ddb)))*(exp(mu+dda+period)/(1+exp(mu+dda+period)));
response=rantbl(1352121,p10,p01,1-p10-p01);
if response eq 1 then y610+1;
if response eq 2 then y601+1;
end; ** end of do i;
*** Weighted least squares test for treatment A versus placebo;
snap=0;
sdap=0;
if y110 gt 0 and y101 gt 0 and y210 gt 0 and y201 gt 0 then do;
wt=4/(1/y110+1/y101+1/y210+1/y201);
snap+wt*log((y101*y210)/(y110*y201))/2;
sdap+wt; end;
else do;
wt=4/(1/(y110+0.5)+1/(y101+0.5)+1/(y210+0.5)+1/(y201+0.5));
snap+wt*log(((y101+0.5)*(y210+0.5))/((y110+0.5)*(y201+0.5)))/2;
sdap+wt; end;
if y310 gt 0 and y301 gt 0 and y510 gt 0 and y501 gt 0 then do;
wt=1/(1/y310+1/y301+1/y510+1/y501);
snap+wt*log((y301*y510)/(y310*y501));
sdap+wt; end;
else do;
wt=1/(1/(y310+0.5)+1/(y301+0.5)+1/(y510+0.5)+1/(y501+0.5));
snap+wt*log(((y301+0.5)*(y510+0.5))/((y310+0.5)*(y501+0.5)));
sdap+wt; end;
if y610 gt 0 and y601 gt 0 and y410 gt 0 and y401 gt 0 then do;
wt=1/(1/y610+1/y601+1/y410+1/y401);
snap+wt*log((y601*y410)/(y610*y401));
sdap+wt; end;
else do;
wt=1/(1/(y610+0.5)+1/(y601+0.5)+1/(y410+0.5)+1/(y401+0.5));
snap+wt*log(((y601+0.5)*(y410+0.5))/((y610+0.5)*(y401+0.5)));
sdap+wt; end;
twls=snap**2/sdap;
pap=1-probchi(twls,1);
if pap le alpha then sumreja+1;
** Mantel-Hanszel Test for treatment A versus placebo;
a1=y101; t1=y110+y101;
a2=y201; t2=y210+y201;
a3=y301; t3=y310+y301;
a4=y401; t4=y410+y401;
a5=y501; t5=y510+y501;
a6=y601; t6=y610+y601;
a12=a1+a2;
a35=a3+a5;
a46=a4+a6;
t12=t1+t2;
t35=t3+t5;
t46=t4+t6;
soap=0;
seap=0;
svap=0;
if t1 gt 0 and t2 gt 0 then do;
soap+a1;
seap+a12*t1/t12;
vap1=a12*(t12-a12)*t1*t2/(t12**2*(t12-1));
svap+vap1;
end;
if t3 gt 0 and t5 gt 0 then do;
soap+a3;
seap+a35*t3/t35;
vap2=a35*(t35-a35)*t3*t5/(t35**2*(t35-1));
svap+vap2;
end;
if t4 gt 0 and t6 gt 0 then do;
soap+a6;
seap+a46*t6/t46;
vap3=a46*(t46-a46)*t4*t6/(t46**2*(t46-1));
svap+vap3;
end;
if svap gt 0 then do;
tap=(soap-seap)**2/svap;
pap=1-probchi(tap,1);
if pap le alpha then sumrej2+1;
end; ** end of if svap gt 0;
if svap le 0 then sumf2+1;
*** exact test for treatment A versus placebo;
num1=lgamma(t1+1)-lgamma(a1+1)-lgamma(t1-a1+1)+lgamma(t2+1)-lgamma(a12-a1+1)-lgamma(t2-a12+a1+1);
den1=lgamma(t1+t2+1)-lgamma(a12+1)-lgamma(t1+t2-a12+1);
prob1=exp(num1-den1);
num2=lgamma(t5+1)-lgamma(a35-a3+1)-lgamma(t5-a35+a3+1)+lgamma(t3+1)-lgamma(a3+1)-lgamma(t3-a3+1);
den2=lgamma(t3+t5+1)-lgamma(a35+1)-lgamma(t3+t5-a35+1);
prob2=exp(num2-den2);
num3=lgamma(t6+1)-lgamma(a6+1)-lgamma(t6-a6+1)+lgamma(t4+1)-lgamma(a46-a6+1)-lgamma(t4-a46+a6+1);
den3=lgamma(t4+t6+1)-lgamma(a46+1)-lgamma(t4+t6-a46+1);
prob3=exp(num3-den3);
prob0=prob1*prob2*prob3;
cumpro=0;
do i=max(0,a12-t2) to min(a12,t1);
do j=max(0,a35-t5) to min(a35,t3);
do k=max(0,a46-t4) to min(a46,t6);
num1=lgamma(t1+1)-lgamma(i+1)-lgamma(t1-i+1)+lgamma(t2+1)-lgamma(a12-i+1)-lgamma(t2-a12+i+1);
den1=lgamma(t1+t2+1)-lgamma(a12+1)-lgamma(t1+t2-a12+1);
prob1=exp(num1-den1);
num2=lgamma(t5+1)-lgamma(a35-j+1)-lgamma(t5-a35+j+1)+lgamma(t3+1)-lgamma(j+1)-lgamma(t3-j+1);
den2=lgamma(t3+t5+1)-lgamma(a35+1)-lgamma(t3+t5-a35+1);
prob2=exp(num2-den2);
num3=lgamma(t6+1)-lgamma(k+1)-lgamma(t6-k+1)+lgamma(t4+1)-lgamma(a46-k+1)-lgamma(t4-a46+k+1);
den3=lgamma(t4+t6+1)-lgamma(a46+1)-lgamma(t4+t6-a46+1);
prob3=exp(num3-den3);
prob=prob1*prob2*prob3;
if prob le prob0 then cumpro+prob;
end;
end;
end;
if cumpro le alpha then sumrej1+1;
*** Weighted least squares test for treatment B versus placebo;
snap=0;
sdap=0;
if y310 gt 0 and y301 gt 0 and y410 gt 0 and y401 gt 0 then do;
wt=4/(1/y310+1/y301+1/y410+1/y401);
snap+wt*log((y301*y410)/(y310*y401))/2;
sdap+wt; end;
else do;
wt=4/(1/(y310+0.5)+1/(y301+0.5)+1/(y410+0.5)+1/(y401+0.5));
snap+wt*log(((y301+0.5)*(y410+0.5))/((y310+0.5)*(y401+0.5)))/2;
sdap+wt; end;
if y510 gt 0 and y501 gt 0 and y210 gt 0 and y201 gt 0 then do;
wt=1/(1/y510+1/y501+1/y210+1/y201);
snap+wt*log((y501*y210)/(y510*y201));
sdap+wt; end;
else do;
wt=1/(1/(y510+0.5)+1/(y501+0.5)+1/(y210+0.5)+1/(y201+0.5));
snap+wt*log(((y501+0.5)*(y210+0.5))/((y510+0.5)*(y201+0.5)));
sdap+wt; end;
if y110 gt 0 and y101 gt 0 and y610 gt 0 and y601 gt 0 then do;
wt=1/(1/y110+1/y101+1/y610+1/y601);
snap+wt*log((y101*y610)/(y110*y601));
sdap+wt; end;
else do;
wt=1/(1/(y110+0.5)+1/(y101+0.5)+1/(y610+0.5)+1/(y601+0.5));
snap+wt*log(((y101+0.5)*(y610+0.5))/((y110+0.5)*(y601+0.5)));
sdap+wt; end;
twls=snap**2/sdap;
pap=1-probchi(twls,1);
if pap le alpha then sumrejb+1;
** Mantel-Hanszel test for treatment B versus placebo;
a34=a3+a4;
a25=a2+a5;
a16=a1+a6;
t34=t3+t4;
t25=t2+t5;
t16=t1+t6;
sobp=0;
sebp=0;
svbp=0;
if t1 gt 0 and t6 gt 0 then do;
sobp+a1;
sebp+a16*t1/t16;
vbp1=a16*(t16-a16)*t1*t6/(t16**2*(t16-1));
svbp+vbp1;
end;
if t3 gt 0 and t4 gt 0 then do;
sobp+a3;
sebp+a34*t3/t34;
vbp2=a34*(t34-a34)*t3*t4/(t34**2*(t34-1));
svbp+vbp2;
end;
if t2 gt 0 and t5 gt 0 then do;
sobp+a5;
sebp+a25*t5/t25;
vbp3=a25*(t25-a25)*t2*t5/(t25**2*(t25-1));
svbp+vbp3;
end;
if svbp gt 0 then do;
tbp=(sobp-sebp)**2/svbp;
pbp=1-probchi(tbp,1);
if pbp le alpha then sumrej4+1;
end; ** end of if svbp gt 0;
if svbp le 0 then sumf4+1;
** exact test for treatment B verus placebo;
num1=lgamma(t3+1)-lgamma(a3+1)-lgamma(t3-a3+1)+lgamma(t4+1)-lgamma(a34-a3+1)-lgamma(t4-a34+a3+1);
den1=lgamma(t3+t4+1)-lgamma(a34+1)-lgamma(t3+t4-a34+1);
prob1=exp(num1-den1);
num2=lgamma(t5+1)-lgamma(a5+1)-lgamma(t5-a5+1)+lgamma(t2+1)-lgamma(a25-a5+1)-lgamma(t2-a25+a5+1);
den2=lgamma(t2+t5+1)-lgamma(a25+1)-lgamma(t2+t5-a25+1);
prob2=exp(num2-den2);
num3=lgamma(t6+1)-lgamma(a16-a1+1)-lgamma(t6-a16+a1+1)+lgamma(t1+1)-lgamma(a1+1)-lgamma(t1-a1+1);
den3=lgamma(t1+t6+1)-lgamma(a16+1)-lgamma(t1+t6-a16+1);
prob3=exp(num3-den3);
prob0=prob1*prob2*prob3;
cumpro=0;
do i=max(0,a34-t4) to min(a34,t3);
do j=max(0,a25-t2) to min(a25,t5);
do k=max(0,a16-t6) to min(a16,t1);
num1=lgamma(t3+1)-lgamma(i+1)-lgamma(t3-i+1)+lgamma(t4+1)-lgamma(a34-i+1)-lgamma(t4-a34+i+1);
den1=lgamma(t3+t4+1)-lgamma(a34+1)-lgamma(t3+t4-a34+1);
prob1=exp(num1-den1);
num2=lgamma(t5+1)-lgamma(j+1)-lgamma(t5-j+1)+lgamma(t2+1)-lgamma(a25-j+1)-lgamma(t2-a25+j+1);
den2=lgamma(t2+t5+1)-lgamma(a25+1)-lgamma(t2+t5-a25+1);
prob2=exp(num2-den2);
num3=lgamma(t6+1)-lgamma(a16-k+1)-lgamma(t6-a16+k+1)+lgamma(t1+1)-lgamma(k+1)-lgamma(t1-k+1);
den3=lgamma(t1+t6+1)-lgamma(a16+1)-lgamma(t1+t6-a16+1);
prob3=exp(num3-den3);
prob=prob1*prob2*prob3;
if prob le prob0 then cumpro+prob;
end;
end;
end;
if cumpro le alpha then sumrej3+1;
end; ** end of do nsimul;
err1=sumrej1/nsimul;
err2=sumrej2/nsimul;
erra=sumreja/nsimul;
err3=sumrej3/nsimul;
err4=sumrej4/nsimul;
errb=sumrejb/nsimul;
put sigma 1-5 1 dda 6-10 2 ddb 11-16 2 n 17-21 erra 22-30 3 err2 31-39 3 err1 40-48 3 errb 49-57 3 err4 58-66 3 err3 67-75 3;
end;
end;
end;
end;
end;