This is my first time asking a question about R, so forgive me if I leave out details that are necessary. I'll try to be as thorough but concise as possible.
My data has 4 columns: score (values of -1, 0, or 1); species (at least 10 different ones); habitat (three different ones); and site (at least 8 different ones). So each row is a unique combination of score, species, habitat, and site. See below for top rows:
score species habitat site
0 NOCA NAG NAG1
0 BRTH NAG NAG1
0 BARS NAG NAG1
1 COYE NAG NAG1
0 HOWR NAG NAG1
0 SAVS NAG NAG1
0 CEDW NAG NAG1
1 CHSP NAG NAG1
0 EAKI NAG NAG1
1 MODO NAG NAG1
0 NOCA NAG NAG16
0 BRTH NAG NAG16
I ran an lme to evaluate the effect of species and habitat on score, with site being a random effect. Code was as follows:
anova(model<-lme(score~species*habitat, random=~1|site))
I then wanted to do pairwise comparisons to figure out where the significant differences are. TukeyHSD doesn't work with lme, so I used the multcomp package as follows:
summary(glht(model, linfct=mcp(habitat="Tukey")))
Which as far as I can tell would do pairwise comparisons of scores based on habitat type. I also tried substituting species in for habitat to do comparisons based on species. In both cases I got the following error:
Error in contrMat(table(mf[[nm]]), type = types[pm]) : less than two groups
I don't see how I can have less than two groups anywhere, since the minimum number of groups I have for any column is three. Maybe I am misinterpreting what "groups" R is referring to? Any suggestions on what might be going wrong and how to fix it?
EDIT As per recommendations (thank you), here is some more information about my data that might make it reproducible. First, the code I used in its entirety.
library(nlme)
library(multcomp)
null<-read.csv("baci_null_red.csv", head=TRUE)
attach(null)
anova(model<-lme(score~species*habitat, random=~1|site))
numDF denDF F-value p-value
(Intercept) 1 1392 1.929021 0.1651
species 29 1392 2.691207 <.0001
habitat 2 48 3.412485 0.0411
species:habitat 58 1392 1.239267 0.1099
summary(glht(model, linfct=mcp(habitat="Tukey")))
Error in contrMat(table(mf[[nm]]), type = types[pm]) :
less than two groups
Here are some summaries about the data that may help in reproducing it...it's just kind of big and I'm not sure how to make it minimal and yet still get the error.
str(null)
'data.frame': 1530 obs. of 4 variables:
$ score : int 0 0 -1 0 0 0 0 0 0 0 ...
$ species: Factor w/ 30 levels "AMGO","AMRO",..: 27 5 7 2 18 12 22 28 6 13 ...
$ habitat: Factor w/ 3 levels "CUM","IAG","NAG": 3 3 3 3 3 3 3 3 3 3 ...
$ site : Factor w/ 51 levels " CUM6","CUM1",..: 37 37 37 37 37 37 37 37 37 37 ..
summary(null)
score species habitat site
Min. :-1.00000 AMGO : 51 CUM:210 CUM6 : 30
1st Qu.: 0.00000 AMRO : 51 IAG:660 CUM1 : 30
Median : 0.00000 BARS : 51 NAG:660 CUM10 : 30
Mean :-0.01569 BCCH : 51 CUM12 : 30
3rd Qu.: 0.00000 BHCO : 51 CUM2 : 30
Max. : 1.00000 BLJA : 51 CUM3 : 30
(Other):1224 (Other):1350
dput(null[somerows,c("species","habitat","site")])
structure(list(species = structure(c(27L, 5L, 7L, 2L, 18L, 12L,
22L, 28L, 6L, 13L, 15L, 23L, 4L, 10L, 14L, 11L, 20L, 30L, 9L,
19L, 26L, 8L, 16L, 25L, 3L, 17L, 21L, 29L, 24L, 1L, 27L, 5L,
7L, 2L, 18L, 12L, 22L, 28L, 6L, 13L, 15L, 23L, 4L, 10L, 14L,
11L, 20L, 30L, 9L, 19L, 26L, 8L, 16L, 25L, 3L, 17L, 21L, 29L,
24L, 1L, 27L, 5L, 7L, 2L, 18L, 12L, 22L, 28L, 6L, 13L, 15L, 23L,
4L, 10L, 14L, 11L, 20L, 30L, 9L, 19L, 26L, 8L, 16L, 25L, 3L,
17L, 21L, 29L, 24L, 1L, 27L, 5L, 7L, 2L, 18L, 12L, 22L, 28L,
6L, 13L), .Label = c("AMGO", "AMRO", "BARS", "BCCH", "BHCO",
"BLJA", "BOBO", "BRTH", "CEDW", "CSWA", "EABL", "EAKI", "EAWP",
"FISP", "GCFL", "HOLA", "HOWR", "INBU", "KILL", "LEFL", "MODO",
"NOCA", "RBGR", "RWBL", "SAVS", "SOSP", "VESP", "WIFL", "YSFL",
"YWAR"), class = "factor"), habitat = structure(c(3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L), .Label = c("CUM", "IAG", "NAG"), class = "factor"), site = structure(c(37L,
37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L,
37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L,
37L, 37L, 37L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L,
46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L, 46L,
46L, 46L, 46L, 46L, 46L, 46L, 46L, 47L, 47L, 47L, 47L, 47L, 47L,
47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L,
47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 48L, 48L,
48L, 48L, 48L, 48L, 48L, 48L, 48L, 48L), .Label = c(" CUM6",
"CUM1", "CUM10", "CUM12", "CUM2", "CUM3", "CUM8", "IAG1", "IAG10",
"IAG13", "IAG14", "IAG15", "IAG16", "IAG18", "IAG19", "IAG21",
"IAG22", "IAG23", "IAG24", "IAG25", "IAG26", "IAG27", "IAG28",
"IAG3", "IAG4", "IAG5", "IAG6", "IAG8", "IAG9", "NAG10", "NAG11",
"NAG13", "NAG14", "NAG15", "NAG18", "NAG19", "NAG2", "NAG21",
"NAG22", "NAG23", "NAG24", "NAG25", "NAG26", "NAG27", "NAG28",
"NAG3", "NAG4", "NAG5", "NAG6", "NAG7", "NAG8"), class = "factor")), .Names =
c("species", "habitat", "site"), row.names = c(NA, 100L), class = "data.frame")
There you have it :s thanks for any guidance at all.