# Help in Two-way ANOVA with Significant Interaction

I got the following results for an experiment of fermentation to produce Butanol:

Butanol =c(11.7462,11.7904,11.9162,11.8732,11.8583,11.8677,11.8697,
11.7289,11.9296,11.7722,11.9813,11.9873,11.8058,11.8711,11.937,
11.8628,11.7786,11.7649,11.8459,11.9139,12.0537,12.1359,11.9949,
12.0752,11.9993,12.2802,12.2227,12.1274,12.1408,11.9896,12.1362,
12.265,12.1353,12.0812,12.511,12.1871,11.7881,12.1962,12.2482,12.189)

Strain = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2)

Enzyme = c("AH","AH","AH","AH","AH","AA","AA","AA","AA","AA",
"AL","AL","AL","AL","AL","B","B","B","B","B","C","C","C","C","C",
"AA","AA","AA","AA","AA","B","B","B","B","B","C","C","C","C","C")

myData = data.frame(Butanol, Strain, Enzyme)

myData$$Strain = factor(myData$$Strain)

• The first column is the butanol produced (outcome) in the process;

• The second column is bacteria strain since I've conducted the experiment with two different strains;

• The third column is the enzyme to accelerate the butanol production, here I have five levels: AH (high dose of enzyme A), AA (average dose of enzyme A), AL (low dose of enzyme A), B (typical dose of enzyme B) and C (typical dose of enzyme C).

• 5 replicates per treatment

Then, I ran the following code:

anova(aov(Butanol ~ Strain*Enzyme,myData))

Only interaction and Strain were significant, so I'm a little bit confused how to proceed.

• Am I using the correct SS Type (in my case is 1)? Firstly, I was in doubt if my data was unbalanced because I did not perform every factorial combination, so I would need to use SS Type 3.

• What should I do next? I think I should filter treatments with strain 1 and treatments with strain 2 and then, run anova for these subsets, right? Or should I run post-hoc directly for interaction?

• I don't seem to get the interaction is signficant although the main effects are both signficant
– Emi
May 5 '19 at 9:56

Note: Revised with revised data

This may be somewhat of an R-specific answer.

By default, R uses Type I sums of squares. This is fine for balanced designs or when sequential sums of squares are desired, but people often want to use Type II or Type III sums of squares. As a matter of course I recommend using the lm function which fits a general linear model, and the Anova function in the car package for the anova table. This function defaults to Type II sums of squares, which is a good choice for a default. Another advantage of this approach is that it allows you to use the emmeans package for post-hoc comparisons.

In the case of your data, the effects are not completely crossed. That is, Enzymes AH and AL are measured only in Strain 1. This will create some problems in estimating effects.

The following code will run at: rdrr.io/snippets/, or in R.

Install packages and assemble data

if(!require(car)){install.packages("car")}
if(!require(FSA)){install.packages("FSA")}
if(!require(emmeans)){install.packages("emmeans")}

Butanol =c(11.7462,11.7904,11.9162,11.8732,11.8583,11.8677,11.8697,
11.7289,11.9296,11.7722,11.9813,11.9873,11.8058,11.8711,11.937,
11.8628,11.7786,11.7649,11.8459,11.9139,12.0537,12.1359,11.9949,
12.0752,11.9993,12.2802,12.2227,12.1274,12.1408,11.9896,12.1362,
12.265,12.1353,12.0812,12.511,12.1871,11.7881,12.1962,12.2482,12.189)

Strain = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2)

Enzyme = c("AH","AH","AH","AH","AH","AA","AA","AA","AA","AA",
"AL","AL","AL","AL","AL","B","B","B","B","B","C","C","C","C","C",
"AA","AA","AA","AA","AA","B","B","B","B","B","C","C","C","C","C")

myData = data.frame(Butanol, Strain, Enzyme)

myData$$Strain = factor(myData$$Strain)


Summarize data. Note that Enzyme AH and AL are measured only in Strain 1.

library(FSA)

Summarize(Butanol ~ Strain + Enzyme, myData)

#   Strain Enzyme n     mean        sd     min      Q1  median      Q3     max
# 1      1     AA 5 11.83362 0.0812620 11.7289 11.7722 11.8677 11.8697 11.9296
# 2      2     AA 5 12.15214 0.1101567 11.9896 12.1274 12.1408 12.2227 12.2802
# 3      1     AH 5 11.83686 0.0679207 11.7462 11.7904 11.8583 11.8732 11.9162
# 4      1     AL 5 11.91650 0.0773750 11.8058 11.8711 11.9370 11.9813 11.9873
# 5      1      B 5 11.83322 0.0616360 11.7649 11.7786 11.8459 11.8628 11.9139
# 6      2      B 5 12.22574 0.1732074 12.0812 12.1353 12.1362 12.2650 12.5110
# 7      1      C 5 12.05180 0.0583478 11.9949 11.9993 12.0537 12.0752 12.1359
# 8      2      C 5 12.12172 0.1881807 11.7881 12.1871 12.1890 12.1962 12.2482


Fit a general linear model. Despite the fact that R returns an anova table, caution should be used since the design is not fully crossed.

model = lm(Butanol ~ Strain*Enzyme, myData)

library(car)

Anova(model)

# Anova Table (Type II tests)
#
# Response: Butanol
#                Sum Sq Df F value    Pr(>F)
# Strain        0.50825  1 39.8805 4.373e-07 ***
# Enzyme        0.06587  4  1.2921  0.293837
# Strain:Enzyme 0.14279  2  5.6021  0.008203 **
# Residuals     0.40782 32


Compare estimated marginal means. Note some cannot be estimated.

library(emmeans)

marginal = emmeans(model, ~ Strain:Enzyme)

pairs(marginal)

CLD(marginal, Letters=letters)

#  Strain Enzyme emmean     SE df lower.CL upper.CL .group
#  1      B        11.8 0.0505 32     11.7     11.9  a
#  1      AA       11.8 0.0505 32     11.7     11.9  a
#  1      AH       11.8 0.0505 32     11.7     11.9  a
#  1      AL       11.9 0.0505 32     11.8     12.0  ab
#  1      C        12.1 0.0505 32     11.9     12.2  abc
#  2      C        12.1 0.0505 32     12.0     12.2   bc
#  2      AA       12.2 0.0505 32     12.0     12.3   bc
#  2      B        12.2 0.0505 32     12.1     12.3    c
#  2      AH     nonEst     NA NA       NA       NA
#  2      AL     nonEst     NA NA       NA       NA


Also note that the joint_tests function warns of problems.

library(emmeans)

joint_tests(model)

#  model term    df1 df2 F.ratio p.value note
#  Enzyme          2  32   1.757 0.1888     e
#  Strain:Enzyme   2  32   5.602 0.0082     e
#
#  e: df1 reduced due to non-estimability


Perhaps one idea is to look at Enzymes only within Strains.

library(emmeans)

marginal = emmeans(model, ~ Enzyme|Strain)

pairs(marginal)

CLD(marginal, Letters=letters)

#  Strain = 1:
#   Enzyme emmean     SE df lower.CL upper.CL .group
#   B        11.8 0.0505 32     11.7     11.9  a
#   AA       11.8 0.0505 32     11.7     11.9  a
#   AH       11.8 0.0505 32     11.7     11.9  a
#   AL       11.9 0.0505 32     11.8     12.0  ab
#   C        12.1 0.0505 32     11.9     12.2   b
#
#  Strain = 2:
#   Enzyme emmean     SE df lower.CL upper.CL .group
#   C        12.1 0.0505 32     12.0     12.2  a
#   AA       12.2 0.0505 32     12.0     12.3  a
#   B        12.2 0.0505 32     12.1     12.3  a
#   AH     nonEst     NA NA       NA       NA
#   AL     nonEst     NA NA       NA       NA
#
#  Confidence level used: 0.95
#  P value adjustment: tukey method for comparing a family of 5 estimates
#  significance level used: alpha = 0.05


Check some model assumptions

hist(residuals(model), col="darkgray")

plot(predict(model), residuals(model))


Type III anova appears to not produce output in this case

model3 = lm(Butanol ~ Strain*Enzyme, myData,
contrasts=list(Strain=contr.sum, Enzyme=contr.sum))

library(car)

Anova(model3, type=3)

# Error in Anova.III.lm(mod, error, singular.ok = singular.ok, ...) :
#   there are aliased coefficients in the model


Perhaps one approach would be to fit a model with only the interaction as an independent variable. In general, this isn't an ideal approach, but in this case it might be a reasonable approach.

myData$$Int = interaction(myData$$Strain, myData$Enzyme) model2 = lm(Butanol ~ Int, myData) Anova(model2) # Anova Table (Type II tests) # # Response: Butanol # Sum Sq Df F value Pr(>F) # Int 0.90294 7 10.121 1.294e-06 *** # Residuals 0.40782 32 marginal2 = emmeans(model2, ~ Int) CLD(marginal2, Letters=letters) # Int emmean SE df lower.CL upper.CL .group # 1.B 11.8 0.0505 32 11.7 11.9 a # 1.AA 11.8 0.0505 32 11.7 11.9 a # 1.AH 11.8 0.0505 32 11.7 11.9 a # 1.AL 11.9 0.0505 32 11.8 12.0 ab # 1.C 12.1 0.0505 32 11.9 12.2 abc # 2.C 12.1 0.0505 32 12.0 12.2 bc # 2.AA 12.2 0.0505 32 12.0 12.3 c # 2.B 12.2 0.0505 32 12.1 12.3 c # # Confidence level used: 0.95 # P value adjustment: tukey method for comparing a family of 8 estimates # significance level used: alpha = 0.05  Check some model assumptions hist(residuals(model2), col="darkgray") plot(predict(model2), residuals(model2))  • Sal, thank you for your reply. I typed some results wrong, that's why everybody couldn't see the significant interaction. May I edit things? Did you consider Anova SS type 2 because of the not fully crossed? May 5 '19 at 15:03 • @user2501348 , Yes, if you change the data I can update the code and results... The reason I chose type 2 is because that's usually what people want. I think it's unfortunate that R defaults to type 1. May 5 '19 at 15:26 • Thank you! I've changed the data. Could I compare if there were differences in doses of Enzyme A first and then, consider it as one group? If the answer is yes, should I include the Strain 2 with AA together with Strain 1 with AA, AL, AH for one-way Anova? I have problems when thinking on subsets because I am afraid of change the results. May 5 '19 at 17:43 • I updated the results. With the revised data, the model behaves a little better for whatever reason.... Within Strain 1, the A enzymes are statistically similar. You don't know if the different levels of A have an effect in Strain 2. I would resist the urge to combine the treatments. Perhaps report the result that the different levels of A were similar in Strain 1, and then refit the model without AH and AL. I think you will find some difference in opinions in how to proceed in cases like this. But this gives you a fully crossed model. May 5 '19 at 20:39 • You can create the reduced data frame with myData2 = myData[myData$Enzyme=="B"|myData$Enzyme=="C"|myData$Enzyme=="AA",] May 5 '19 at 20:39