From what you describe, your response variable is ordinal. Since it is integer valued, it is probably not reasonable to assume normality and use GLM. You can use PROC LOGISTIC to model an ordinal multinomial response like yours. Using the data you provided, these statements fit the model allowing for interaction in the factors. It simultaneously models the cumulative probabilities - Pr(score=3), Pr(score=2 or 3), Pr(score=1,2, or 3). This ordering and grouping of level is the result of the DESC response variable option. The ILINK option in the LSMEANS statement provides the estimated cumulative probabilities for each A-B combination (in the Mean column). In the results you'll see that there is a significant interaction. As with any model, the interaction means that the effect of one factor changes with the level of the other factor. So, it is necessary to assess the effect of each factor at each level of the other factor. The effect of REP is not significant and could possibly be dropped from the model but is left in for the following.
proc logistic data=x;
class factora factorb / param=glm;
model score(desc)=factora|factorb rep;
lsmeans factora*factorb/ilink plots=none;
store mod2;
run;
You can visually plot the interaction effect using the fitted and stored model (in MOD2) and then using the EFFECTPLOT statement in PROC PLM to produce the plot. The following statements show the interaction two different but equivalent ways. Note that there is a plot for each of the three cumulative probabilities.
proc plm source=mod2;
effectplot interaction(x=factorb sliceby=factora);
effectplot interaction(x=factora sliceby=factorb);
run;
To assess the interaction, it is reasonable to use the approach you mentioned which is to include just the 4 level TRT variable in the model and then make suitable comparisons among the factor combinations. That is done by the following, but a better labeled variable, AB, is used instead of TRT.
data x; set x;
ab=cats(factora,factorb);
run;
proc logistic data=x;
class ab / param=glm;
model score(desc)=ab rep;
lsmeans ab / ilink plots=none;
store mod;
run;
Note that the model is equivalent to the first one - same likelihood value, same estimated cumulative probabilities. Using the order of the factor combinations shown in the LSMEANS table, you can use LSMESTIMATE statements to easily test the effects within the interaction as described above.
proc plm source=mod;
effectplot interaction(x=ab sliceby=score);
lsmestimate ab 'A1-A2 in Yes' 0 1 0 -1 / joint e;
lsmestimate ab 'A1-A2 in No' 1 0 -1 0 / joint e;
lsmestimate ab 'diff in diff' -1 1 1 -1 / joint e; *A1-A2 in Yes - A1-A2 in No;
run;
The EFFECTPLOT statement shows the cumulative probabilities for each of the 4 factor combinations. The first two LSMESTIMATE statements test the effect of A in each level of B. The last LSMESTIMATE statement tests the so-called "difference in difference" effect which is the difference of the two differences in the first two statements. This is, in fact, the interaction effect and you'll see that it has the same p-value as from the first model fit. The results form the first two statements show that the effect of A is significant in B=Yes but not in B=No. You can easily see that best in the plot, but also in the interaction plots earlier.
... View more