* conditional logistic regression;
* separate fixed effect for each clinic;
proc logistic data=cream; 
 class center; 
 model y/n = treat center;
run;

* conditional logistic regression;
* stratified logistic regression;
proc logistic data=cream; 
 class center; 
 model y/n = treat;
 strata center;
run;

* conditional analysis using glimmix;
proc glimmix data=cream3; class center treat;
 model y(event='1') = treat / dist=binary link=logit solution;
 random intercept / subject=center;
run;

* marginal analysis using GEE;
proc genmod data=cream3 descending;
 class center treat; 
 model y = treat / dist=bin link=logit type3;
 repeated subject=center / type=exch;
 estimate "OR" treat 1 -1 / exp;
run;

* marginal analysis using glimmix;
proc glimmix data=cream3; class center treat;
 model y(event='1') = treat / dist=binary link=logit solution;
 random _residual_ / subject=center type=cs;
run;

proc print data=cream2; run;

data cream; * as binomial counts y out of n;
input center$ treat y n; treat=abs(treat-2);
datalines;
 a 1 11 36
 a 2 10 37
 b 1 16 20
 b 2 22 32
 c 1 14 19
 c 2  7 19
 d 1  2 16
 d 2  1 17
 e 1  6 17
 e 2  0 12
 f 1  1 11
 f 2  0 10
 g 1  1  5
 g 2  1  9
 h 1  4  6
 h 2  6  7
 ;

data cream2; * as counts in 2 x 2 x 8 table;
input center$ treat response count;
response=2-response;
datalines; 
 a 1 1 11
 a 1 2 25
 a 2 1 10 
 a 2 2 27
 b 1 1 16
 b 1 2 4 
 b 2 1 22
 b 2 2 10
 c 1 1 14
 c 1 2 5 
 c 2 1 7 
 c 2 2 12
 d 1 1 2 
 d 1 2 14
 d 2 1 1 
 d 2 2 16
 e 1 1 6 
 e 1 2 11
 e 2 1 0 
 e 2 2 12
 f 1 1 1 
 f 1 2 10
 f 2 1 0 
 f 2 2 10
 g 1 1 1 
 g 1 2 4 
 g 2 1 1 
 g 2 2 8 
 h 1 1 4 
 h 1 2 2 
 h 2 1 6 
 h 2 2 1 
;

data cream3; set cream2; * suitable for genmod;
drop response count total i;
do i=1 to count;
 y=response; treat=treat; center=center; output;
end;
proc print;
run;