Summary and Analysis of Extension Program Evaluation in R

Salvatore S. Mangiafico

Factorial ANOVA: Main Effects, Interaction Effects, and Interaction Plots

 

Two-way or multi-way data often come from experiments with a factorial design.  A factorial design has at least two factor variables for its independent variables, and multiple observation for every combination of these factors.

 

The weight gain example below shows factorial data.  In this example, there are three observations for each combination of Diet and Country.

 

With this kind of data, we are usually interested in testing the effect of each factor variable (main effects) and then the effect of their combination (interaction effect).

 

For two-way data, an interaction plot shows the mean or median value for the response variable for each combination of the independent variables.  This type of plot, especially if it includes error bars to indicate the variability of data within each group, gives us some understanding of the effect of the main factors and their interaction. 

 

When main effects or interaction effects are statistically significant, post-hoc testing can be conducted to determine which groups differ significantly from other groups.  With a factorial experiment, there are a few guidelines for determining when to do post-hoc testing.  The following guidelines are presented for two-way data for simplicity.

 

•  When neither the main effects nor the interaction effect is statistically significant, no post-hoc mean-separation testing should be conducted.

 

•  When one or more of the main effects are statistically significant and the interaction effect is not, post-hoc mean-separation testing should be conducted on significant main effects only. 
(This is shown in the first weight gain example below)

 

•  When the interaction effect is statistically significant, post-hoc mean-separation testing should be conducted on the interaction effect only.  This is the case even when the main effects are also statistically significant.
(This is shown in the second weight gain example below)

 

The final guideline above is not always followed.  There are times when people will present the mean-separation tests for significant main effects even when the interaction effect is significant.  In general, though, if there is a significant interaction, the mean-separation tests for interaction will better explain the results of the analysis, and the mean-separation tests for the main effects will be of less interest.

 

Packages used in this chapter

 

The packages used in this chapter include:

•  psych

•  car

•  multcomp

•  emmeans

•  FSA

•  ggplot2

•  phia

 

The following commands will install these packages if they are not already installed:

if(!require(car)){install.packages("car")}
if(!require(psych)){install.packages("psych")}
if(!require(multcomp)){install.packages("multcomp")}
if(!require(emmeans)){install.packages("emmeans")}
if(!require(FSA)){install.packages("FSA")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(phia)){install.packages("phia")}

Two-way ANOVA example with interaction effect

 

Imagine for this example an experiment in which people were put on one of three diets to encourage weight gain.  The amount of weight gained will be the dependent variable and will be considered an interval/ratio variable.  For independent variables, there are three different diets and the country in which subjects live.

 

The model will be fit with the lm function, whose syntax is similar to that of the clm function.

 

This kind of analysis makes certain assumptions about the distribution of the data, but for simplicity, this example will ignore the need to determine that the data meet these assumptions.

 

Data = read.table(header=TRUE, stringsAsFactors=TRUE, text="

 Diet    Country  Weight_change
 A       USA      0.120
 A       USA      0.125
 A       USA      0.112
 A       UK       0.052
 A       UK       0.055
 A       UK       0.044
 B       USA      0.096
 B       USA      0.100
 B       USA      0.089
 B       UK       0.025
 B       UK       0.029
 B       UK       0.019
 C       USA      0.149
 C       USA      0.150
 C       USA      0.142
 C       UK       0.077
 C       UK       0.080
 C       UK       0.066
")


### Order levels of the factor; otherwise R will alphabetize them

Data$Country = factor(Data$Country,
                      levels=unique(Data$Country))


###  Check the data frame

library(psych)

headTail(Data)

str(Data)

summary(Data)

Simple interaction plot

The interaction.plot function in the native stats package creates a simple interaction plot for two-way data.  The options shown indicate which variables will used for the x-axis, trace variable, and response variable.  The fun=mean option indicates that the mean for each group will be plotted.  For the meaning of other options, see ?interaction.plot.

 

This style of interaction plot does not show the variability of each group mean, so it is difficult to use this style of plot to determine if there are significant differences among groups.

 

The plot shows that mean weight gain for each diet was lower for the UK compared with USA.  And that this difference was relatively constant for each diet, as is evidenced by the lines on the plot being parallel.  This suggests that there is no large or significant interaction effect.  That is, the difference among diets is consistent across countries.  And vice-versa, the difference in countries is consistent across diets.

 

A couple of other styles of interaction plot are shown at the end of this chapter.

 

interaction.plot(x.factor     = Data$Country,
                 trace.factor = Data$Diet,
                 response     = Data$Weight_change,
                 fun = mean,
                 type="b",
                 col=c("black","red","green"),  ### Colors for levels of trace var.
                 pch=c(19, 17, 15),             ### Symbols for levels of trace var.
                 fixed=TRUE,                    ### Order by factor order in data
                 leg.bty = "o")

Specify the linear model and conduct an analysis of variance

A linear model is specified with the lm function.  Weight_change is the dependent variable.  Country and Diet are the independent variables, and including Country:Diet in the formula adds the interaction term for Country and Diet to the model.

 

The ANOVA table indicates that the main effects are significant, but that the interaction effect is not.

 

model = lm(Weight_change ~ Country + Diet + Country:Diet,
           data = Data)

library(car)

Anova(model,
      type = "II")

Anova Table (Type II tests)

               Sum Sq Df  F value    Pr(>F)   
Country      0.022472  1 657.7171 7.523e-12 ***
Diet         0.007804  2 114.2049 1.547e-08 ***
Country:Diet 0.000012  2   0.1756    0.8411   
Residuals    0.000410 12                      

Post-hoc testing with emmeans

Because the main effects were significant, we will want to perform post-hoc mean separation tests for each main effect factor variable.

 

For this, we will use the emmeans package.  The linear model under consideration is called model, created the lm function above.  The formula in the emmeans function indicates that comparisons should be conducted for the variable Country in the first call, and for the variable Diet in the second call.

library(emmeans)

marginal = emmeans(model,
                   ~ Country)

pairs(marginal,
      adjust="tukey")

 contrast  estimate      SE  df  t.ratio  p.value
 USA - UK    0.0707 0.00276  12   25.646   <.0001

Results are averaged over the levels of: Diet

library(emmeans)

marginal = emmeans(model,
                   ~ Diet)

pairs(marginal,
      adjust="tukey")

 
contrast estimate      SE df t.ratio p.value
 A - B       0.025 0.00337 12   7.408  <.0001
 A - C      -0.026 0.00337 12  -7.704  <.0001
 B - C      -0.051 0.00337 12 -15.112  <.0001

Results are averaged over the levels of: Country
P value adjustment: tukey method for comparing a family of 3 estimates

Extended example with additional country

For the following example, the hypothetical data have been amended to include a third country, New Zealand.

Data = read.table(header=TRUE, stringsAsFactors=TRUE, text="

 Diet    Country  Weight_change
 A       USA      0.120
 A       USA      0.125
 A       USA      0.112
 A       UK       0.052
 A       UK       0.055
 A       UK       0.044
 A       NZ       0.080
 A       NZ       0.090
 A       NZ       0.075
 B       USA      0.096
 B       USA      0.100
 B       USA      0.089
 B       UK       0.025
 B       UK       0.029
 B       UK       0.019
 B       NZ       0.055
 B       NZ       0.065
 B       NZ       0.050
 C       USA      0.149
 C       USA      0.150
 C       USA      0.142
 C       UK       0.077
 C       UK       0.080
 C       UK       0.066
 C       NZ       0.055
 C       NZ       0.065
 C       NZ       0.050
 C       NZ       0.054
")


### Order levels of the factor; otherwise R will alphabetize them

Data$Country = factor(Data$Country,
                      levels=unique(Data$Country))


###  Check the data frame

library(psych)

headTail(Data)

str(Data)

summary(Data)

Simple interaction plot

The plot suggests that the effect of diet is not consistent across all three countries.  While Diet C showed the greatest mean weight gain for USA and UK, for NZ it has a lower mean than Diet A.  This suggests there may be a meaningful or significant interaction effect, but we will need to do a statistical test to confirm this hypothesis.

interaction.plot(x.factor     = Data$Country,
                 trace.factor = Data$Diet,
                 response     = Data$Weight_change,
                 fun = mean,
                 type="b",
                 col=c("black","red","green"),  ### Colors for levels of trace var.
                 pch=c(19, 17, 15),             ### Symbols for levels of trace var.
                 fixed=TRUE,                    ### Order by factor order in data
                 leg.bty = "o")

Specify the linear model and conduct an analysis of variance

The ANOVA table indicates that the interaction effect is significant, as are both main effects.

model = lm(Weight_change ~ Country + Diet + Country:Diet,
            data = Data)

library(car)

Anova(model,
      type = "II")

Anova Table (Type II tests)

                Sum Sq Df F value    Pr(>F)   
Country      0.0256761  2 318.715 2.426e-15 ***
Diet         0.0051534  2  63.969 3.634e-09 ***
Country:Diet 0.0040162  4  24.926 2.477e-07 ***
Residuals    0.0007653 19                     

Post-hoc testing with emmeans

Because the interaction effect was significant, we would like to compare all group means from the interaction.  Even though the main effects were significant, the typical advice is to not conduct pairwise comparisons for main effects when their interaction is significant.  This is because focusing on the groups in the interaction better describe the results of the analysis.

 

We will use the emmeans package and ask for a compact letter display with the cld function.  First we create an object, named marginal, with the results of the call to emmeans.  Notice here that the formula indicates that comparisons should be conducted for the interaction of Country and Diet, indicated with Country:Diet.

 

Be sure to read the Estimated Marginal Means for Multiple Comparisons chapter for correct interpretation of estimated marginal means.  For lm model objects, the values for emmean, SE, LCL, and UCL values are meaningful.

library(emmeans)

marginal = emmeans(model,
                   ~ Country:Diet)

pairs(marginal,
      adjust="tukey")

library(multcomp)

cld(marginal,
    alpha=0.05,
    Letters=letters,      ### Use lower-case letters for .group
    adjust="tukey")       ### Tukey-adjusted comparisons 

 Country Diet emmean      SE df lower.CL upper.CL .group 
 UK      B    0.0243 0.00366 19   0.0129   0.0358  a     
 UK      A    0.0503 0.00366 19   0.0389   0.0618   b    
 NZ      C    0.0560 0.00317 19   0.0461   0.0659   b    
 NZ      B    0.0567 0.00366 19   0.0452   0.0681   bc   
 UK      C    0.0743 0.00366 19   0.0629   0.0858    cd  
 NZ      A    0.0817 0.00366 19   0.0702   0.0931     de 
 USA     B    0.0950 0.00366 19   0.0836   0.1064      e 
 USA     A    0.1190 0.00366 19   0.1076   0.1304       f
 USA     C    0.1470 0.00366 19   0.1356   0.1584        g

Confidence level used: 0.95
Conf-level adjustment: sidak method for 9 estimates
P value adjustment: tukey method for comparing a family of 9 estimates
significance level used: alpha = 0.05

### Groups sharing a letter are not significantly different
###   at the alpha = 0.05 level.

Interaction plot with error bars using ggplot2

The interaction plots created with the interaction.plot function above are handy to investigate trends in the data, but have the disadvantage of not showing the variability in data within each group.

 

The package ggplot2 can be used to create attractive interaction plots with error bars.  Here, we will use standard error of each mean for the error bars. 

### Create a data frame called Sum with means and standard deviations

library(FSA)

Sum = Summarize(Weight_change ~ Country + Diet,
                data=Data,
                digits=3)

### Add standard error of the mean to the Sum data frame

Sum$se = Sum$sd / sqrt(Sum$n)

Sum$se = signif(Sum$se, digits=3)

Sum

  Country Diet n nvalid  mean    sd   min    Q1 median    Q3   max percZero      se
1     USA    A 3      3 0.119 0.007 0.112 0.116  0.120 0.122 0.125        0 0.00404
2      UK    A 3      3 0.050 0.006 0.044 0.048  0.052 0.054 0.055        0 0.00346
3      NZ    A 3      3 0.082 0.008 0.075 0.078  0.080 0.085 0.090        0 0.00462
4     USA    B 3      3 0.095 0.006 0.089 0.092  0.096 0.098 0.100        0 0.00346
5      UK    B 3      3 0.024 0.005 0.019 0.022  0.025 0.027 0.029        0 0.00289
6      NZ    B 3      3 0.057 0.008 0.050 0.052  0.055 0.060 0.065        0 0.00462
7     USA    C 3      3 0.147 0.004 0.142 0.146  0.149 0.150 0.150        0 0.00231
8      UK    C 3      3 0.074 0.007 0.066 0.072  0.077 0.078 0.080        0 0.00404
9      NZ    C 4      4 0.056 0.006 0.050 0.053  0.054 0.058 0.065        0 0.00300

### Order levels of the factor; otherwise R will alphabetize them

Sum$Country = factor(Sum$Country,
                     levels=unique(Sum$Country))


### Produce interaction plot

library(ggplot2)

pd = position_dodge(.2)

ggplot(Sum, aes(x = Country,
                y = mean,
                color = Diet)) +
    geom_errorbar(aes(ymin = mean - se,
                      ymax = mean + se),
                   width=.2, size=0.7, position=pd) +
    geom_point(shape=15, size=4, position=pd) +
    theme_bw() +
    theme(axis.title = element_text(face = "bold")) +

    scale_colour_manual(values= c("black","red","green")) +
    ylab("Mean weight change")

 

Interaction plot with mean separation letters manually added

It is common to add mean separation letters from post-hoc analyses to interaction plots.  One option is to add letters manually in either image manipulation software like Photoshop or GIMP, or in a word processor or other software that can handle graphic manipulation. 

 

Letters can also be added to the plot by ggplot2 with the annotate or geom_text options.  See “Optional:  Interaction plot of estimated marginal means with mean separation letters” in the Estimated Marginal Means for Multiple Comparisons chapter for examples.

 

In some cases it is desirable for means to be lettered so that the greatest mean is indicated with a.  However, emmeans by default labels the least mean with a.  The order of letters can be reversed manually, or the reversed=TRUE option can be used.

 

Plot of mean weight change for three diets in three countries.  Means sharing a letter are not significantly different according to pairwise comparisons of estimated marginal means with Tukey adjustment for multiple comparisons.

Interaction plot with the phia package

The phia package can be used to create interaction plots quickly.  By default the error bars indicate standard error of the means.  For additional options, see ?interactionMeans.

 

Note that model is the linear model specified above.

library(phia)

IM = interactionMeans(model)

IM

  Country Diet adjusted mean  std. error
1     USA    A    0.11900000 0.003664274
2      UK    A    0.05033333 0.003664274
3      NZ    A    0.08166667 0.003664274
4     USA    B    0.09500000 0.003664274
5      UK    B    0.02433333 0.003664274
6      NZ    B    0.05666667 0.003664274
7     USA    C    0.14700000 0.003664274
8      UK    C    0.07433333 0.003664274
9      NZ    C    0.05600000 0.003173354

### Produce interaction plot

plot(IM)


### Return the graphics device to its default 1-plot-per-window state

par(mfrow=c(1,1))

©2016 by Salvatore S. Mangiafico.
Rutgers Cooperative Extension, New Brunswick, NJ.

Non-commercial reproduction of this content, with attribution, is permitted.
For-profit reproduction without permission is prohibited.

If you use the code or information in this site in a published work, please cite it as a source.  Also, if you are an instructor and use this book in your course, please let me know.   My contact information is on the About the Author of this Book page.

Citation

Mangiafico, S.S. 2016. Summary and Analysis of Extension Program Evaluation in R, version 1.20.05, revised 2023. rcompanion.org/handbook/. (Pdf version: rcompanion.org/documents/RHandbookProgramEvaluation.pdf.)