Let us return to the greenhouse example with plant species as a predictive factor in addition to fertilizer type. The study then becomes a \(2 \times 4\) factorial as 2 types of plant species and 4 types of fertilizers are investigated.  The total number of experimental units (plants) that are needed now is 48, as \(r=6\) and there are 8 plant species and fertilizer type combinations.

The data might look like this:

The ANOVA table would now be constructed as follows:

The data presented in the table above are in unstacked format. One needs to convert this into a stacked format when attempting to use statistical software. The SAS code is as follows.

data greenhouse_2way;
                input fert $ species $ height;
                datalines;
                control SppA      21.0
                control SppA      19.5
                control SppA      22.5
                control SppA      21.5
                control SppA      20.5
                control SppA      21.0
                control SppB      23.7
                control SppB      23.8
                control SppB      23.8
                control SppB      23.7
                control SppB      22.8
                control SppB      24.4
                f1      SppA      32.0
                f1      SppA      30.5
                f1      SppA      25.0
                f1      SppA      27.5
                f1      SppA      28.0
                f1      SppA      28.6
                f1      SppB      30.1
                f1      SppB      28.9
                f1      SppB      30.9
                f1      SppB      34.4
                f1      SppB      32.7
                f1      SppB      32.7
                f2      SppA      22.5
                f2      SppA      26.0
                f2      SppA      28.0
                f2      SppA      27.0
                f2      SppA      26.5
                f2      SppA      25.2
                f2      SppB      30.6
                f2      SppB      31.1
                f2      SppB      28.1
                f2      SppB      34.9
                f2      SppB      30.1
                f2      SppB      25.5
                f3      SppA      28.0
                f3      SppA      27.5
                f3      SppA      31.0
                f3      SppA      29.5
                f3      SppA      30.0
                f3      SppA      29.2
                f3      SppB      36.1
                f3      SppB      36.6
                f3      SppB      38.7
                f3      SppB      37.1
                f3      SppB      36.8
                f3      SppB      37.1
                ;
                run;
                /*The code to generate the boxplot 
                for distribution of height by species organized by fertilizer 
                in Figure 5.1*/
                proc sort data=greenhouse_2way; by fert species;
                proc boxplot data=greenhouse_2way;
                plot height*species (fert); 
                run;
                

As a preliminary step in Exploratory Data Analysis (EDA), a side-by-side boxplot display of height vs. species organized by fertilizer type would be an ideal graphic. The plot below illustrates, the height differences between species are variable among fertilizer types; see for example the difference in height between SppA and SppB for Control is much less than that for F3. This indicates that fert*species could be a significant interaction prompting a factorial model with interaction.

To run the two-factor factorial model with interaction in SAS proc mixed, we can use:

/*Runs the two-factor factorial model with interaction*/
                proc mixed data=greenhouse_2way method=type3;
                class fert species;
                model height = fert species fert*species;
                store out2way;
                run;
                

Similar to when running the single factor ANOVA, the name of the dataset is specified in the proc mixed statement as well as the method=type3 option that specifies the way the F-test is calculated. The categorical factors fert and species are included in the class statement. The terms (or effects) in the model statement are consistent with the source effects in the layout of the "theoretical" ANOVA table illustrated in 5.1. Finally, the store command stores the elements necessary for the generation of the LS-means interval plot.

Recall the two ANOVA rules, applicable to any model: (1) the df values add up to total df and (2) the sums of squares add up to total sums of squares. As seen by the output below, the df values and sums of squares follow these rules. (It is easy to confirm that the total sum of squares = 1168.7325, by the 2^nd ANOVA rule.)

Applying the above rule for this example, the small p-value of 0.0051 displayed in the table above indicates that the interaction effect is significant suggesting the main effects of either fert and species should not be considered individually. It is the average response differences among the fert and species combinations that matter. In order to determine the statistically significant treatment combinations a suitable multiple comparison procedure (such as the Tukey procedure) can be performed on the LS-means of the interaction effect. 

The necessary follow-up SAS code to perform this procedure is given below. 

ods graphics on;
                proc plm restore=out2way;
                lsmeans fert*species / adjust=tukey plot=(diffplot(center) 
                meanplot(cl ascending)) cl lines;
                /* Because the 2-factor interaction is significant, we work with 
                the means for treatment combination*/
                run; 
                

SAS Output for the LS-means:

Note that the p-values here (Pr > t) are testing the hypotheses that the fert and species combination means = 0. This may be of very little interest. However, a comparison of mean response values for different species and fertilizer combinations may be more beneficial and can be derived from the diffogram shown in Figure 5.2. Again recall that if the confidence interval does not contain zero, then the difference between the two associated means is statistically significant.

Notice also that we see a single value for the standard error based on the MSE from the ANOVA rather than a separate standard error for each mean as we would get from proc summary for the sample means. In this example, with equal sample sizes and no covariates, the LS-means will be identical to the ordinary means displayed in the summary procedure.

There are a total of 8 fert*species combinations resulting in a total of \({8 \choose 2} = 28\) pairwise comparisons. From the diffogram for differences in fert*species combinations, we see that 10 of them are not significant and 18 of them are significant at a 5% level after Tukey adjustment (more about diffograms). The information used to generate the diffogram is presented in the table for differences of fert*species least squares means in the SAS output (this table is not displayed here).  

We can save the differences estimated in SAS proc mixed and utilize proc sgplot to create the plot of differences in mean response for the fert*species combinations as shown in Figure 5.3 (SAS code not provided in these notes). The CIs shown are the Tukey adjusted CIs and the interpretations are similar to what we observed from the diffogram in Figure 5.2.

In addition to comparing differences in mean responses for the fert*species combinations, the SAS code above will also produce the line plot for multiple comparisons of means for fert*species combinations (Figure 5.4) and the plot of means responses organized in the ascending order with 95% CIs for fert*species combinations (Figure 5.5).

The line plot Figure 5.4 connects groups in which the LS-means are not statistically different and displays a summary of which groups have similar means. The plot of means with 95% CIs in Figure 5.5 illustrates the same result, although it uses unadjusted CIs. The organization of the estimated mean responses in ascending order facilitates interpretation of the results. 

Using LS-means subsequent to performing an ANOVA will help to identify the significantly different treatment level combinations. In other words, the analysis of a factorial model does not end with a p-value for an F-test. Rather, a small p-value signals the need for a mean comparison procedure.