I am comparing a fertilizer experiment where I have a response variable (growth rate) with two independent variables (genotype and rate). I ran a 2-way ANOVA using the lsmeans, car, and multcompView packages in R. My interaction effects are not significant, but my main effect variables of genotype and rate are significant. I wish to view letters indicating differences for a mean separation. How can I do that?
Here is what I have tried:
Dataset
demo <- data.frame(genotype=c(1,
1,
1,
1,
2,
2,
2,
2,
3,
3,
3,
3,
2,
2,
1,
3,
3,
3,
1,
3,
2,
2,
1,
1,
1,
2,
2,
1,
3,
3,
3,
3,
2,
1,
2,
1),
rate = c(0,
1,
2,
3,
0,
1,
2,
3,
0,
1,
2,
3,
1,
3,
0,
0,
2,
1,
2,
3,
0,
2,
3,
1,
1,
0,
2,
2,
3,
0,
2,
1,
3,
0,
1,
3),
response = c(0.636008,
0.520401,
0.511821,
0.200255,
0.535433,
0.272521,
0.192943,
0.238736,
0.568422,
0.497654,
0.433995,
0.316064,
0.336663,
0.187805,
0.696257,
0.517592,
0.244133,
0.469655,
0.277319,
0.188577,
0.534307,
0.204349,
0.263632,
0.651846,
0.46279,
0.499629,
0.150601,
0.168777,
0.227221,
0.518858,
0.35837,
0.516537,
0.221283,
0.753765,
0.446882,
0.301356))
demo$genotype <- as.factor(demo$genotype)
demo$rate <- as.factor(demo$rate)
#Load packages
library(car)
library(multcompView)
library(lsmeans)
I previously checked the interaction terms, and they are not significant. Here's the main effect model:
mod1 <- lm(response ~ genotype + rate, data = demo)
Anova(mod1, type = "III")
#Results
Response: response
Sum Sq Df F value Pr(>F)
(Intercept) 2.50289 1 389.5155 < 2.2e-16 ***
genotype 0.11256 2 8.7589 0.001009 **
rate 0.70042 3 36.3345 4.106e-10 ***
Residuals 0.19277 30
---
Looking at differences between means for genotype:
genotype <- lsmeans(mod1, pairwise ~ genotype)
genotype
#Results:
Results are averaged over the levels of: rate
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
1 - 2 0.1353 0.0327 30 4.133 0.0008
1 - 3 0.0489 0.0327 30 1.495 0.3075
2 - 3 -0.0863 0.0327 30 -2.638 0.0340
Results are averaged over the levels of: rate
P value adjustment: tukey method for comparing a family of 3 estimates
Same process for rate:
rate <- lsmeans(mod1, pairwise ~ rate)
rate
#Results:
Results are averaged over the levels of: genotype
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
0 - 1 0.1206 0.0378 30 3.191 0.0165
0 - 2 0.3020 0.0378 30 7.992 <.0001
0 - 3 0.3461 0.0378 30 9.160 <.0001
1 - 2 0.1814 0.0378 30 4.801 0.0002
1 - 3 0.2256 0.0378 30 5.969 <.0001
2 - 3 0.0442 0.0378 30 1.168 0.6509
Results are averaged over the levels of: genotype
P value adjustment: tukey method for comparing a family of 4 estimates
Now, I attempt to get letters by the genotype output, but I get an error message:
cld(genotype,
alpha = 0.05,
Letters = letters)
#Error message thrown here:
Error in cld(genotype, alpha = 0.05, Lettes = letters) :
could not find function "cld"
Is there a better way to be getting mean separation letters, or am I just making a simple mistake?
The solution was to use the agricolae package:
install.packages("agricolae")
library(agricolae)
mod1 <- lm(response ~ genotype + rate, data = demo)
Anova(mod1, type = "III")
genotype <- HSD.test(mod1, "genotype")
genotype
#Part of the output:
response groups
1 0.4536856 a
3 0.4047565 a
2 0.3184293 b
rate <- HSD.test(mod1, "rate")
rate
#Part of the output:
response groups
0 0.5844746 a
1 0.4638832 b
2 0.2824787 c
3 0.2383254 c