data fuelcell;
input vapor $ voltage current @@;
datalines;
Trt 0.50 32.3 Ctl 0.50 28.6
Trt 0.50 32.8 Ctl 0.50 28.5
Trt 0.60 20.3 Ctl 0.60 19.5
Trt 0.60 21.2 Ctl 0.60 19.7
Trt 0.70 10.9 Ctl 0.70 10.6
Trt 0.70 11.1 Ctl 0.70 10.6
Trt 0.80 2.53 Ctl 0.80 2.23
Trt 0.80 2.43 Ctl 0.80 2.20
;
run;
* initial fit of model V;
proc glm data=fuelcell plots=all;
class vapor voltage;
model current = vapor voltage vapor*voltage/ solution;
run;

*Runs the effects model;
proc genmod data=fuelcell plots=all;
*class vapor  voltage (param=effect); *Effect coding for both factors;
*Can also control reference level, e.g;
class vapor (param=ref ref='Ctl') voltage (param=effect);
model current=vapor voltage vapor*voltage/dist=normal link=identity type3;
*Type 3 test is replaced by a LR test;
run;

proc glm data=fuelcell;
class vapor voltage;
model current=vapor|voltage / solution;
lsmeans voltage vapor / pdiff adjust=tukey alpha=0.05 cl;
run;

proc glimmix data=fuelcell;
class vapor voltage;
model current=vapor|voltage / solution;
*Cell means for both levels of vapor at voltage=0.5 and their contrast;
estimate "Trt 0.50V" intercept 1 vapor 0 1 voltage 1 0 0 0 vapor*voltage  0 0 0 0 1 0 0 0,
"Ctl 0.50V" intercept 1 vapor 1 0 voltage 1 0 0 0 vapor*voltage  1 0 0 0 0 0 0 0 ,
"Trt vs Ctl 0.50V" vapor 1 -1 voltage 0 0 0 0 vapor*voltage 1 0 0 0 -1 0 0 0/ adjust=bon alpha=0.1 cl;
run;

data predict;
input days exper years @@;
datalines;
240.0 1 1 206.0 1 1 217.0 1 1 225.0 1 1 110.0 1 2 118.0 1 2 103.0 1 2 95.0 1 2
56.0 1 3 60.0 1 3 68.0 1 3 58.0 1 3 71.0 2 1 53.0 2 1 68.0 2 1 57.0 2 1
47.0 2 2 52.0 2 2 31.0 2 2 49.0 2 2 37.0 2 3 33.0 2 3 40.0 2 3 45.0 2 3
;
run;

* gives all a*b choose 2 pairwise comparisons via Tukey;
* slice subcommand only gives F-test within each slice;
proc glm data=predict plots=all; class exper years;
model days=exper|years;
lsmeans exper*years / adjust=tukey slice=exper; run;
* slice command in glimmix better;
* gives PW comparisons within each slice, i.e.;
* b choose 2 comparisons within each of a slices;
proc glimmix; class exper years;
model days=exper|years;
slice exper*years / sliceby=exper adjust=tukey cl;
slice exper*years / sliceby=years adjust=tukey cl; run;

proc glimmix; class exper years;
model days=exper|years; * adjust=tukey not needed;
slice exper*years / sliceby=years adjust=tukey cl; run;