2013-02-12 69 views
1

这是我第一次问一个关于R的问题,所以如果我遗漏了必要的细节,请原谅我。我会尽量做到尽可能彻底而简洁。我的数据有4列:分数(值-1,0或1);物种(至少10种);栖息地(三个不同的);和网站(至少8个不同的网站)。所以每一行都是得分,物种,栖息地和场地的独特组合。请参阅下面的顶行:错误,当试图进行两两比较lme

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 

我跑了LME评估的物种和栖息地得分的影响,与网站是一个随机效应。代码如下:

anova(model<-lme(score~species*habitat, random=~1|site)) 

然后,我想做两两比较来找出重大差异在哪里。 TukeyHSD不与LME的工作,所以我用了multcomp包如下:

summary(glht(model, linfct=mcp(habitat="Tukey"))) 

其中,据我可以告诉会做基于栖息地类型分数的成对比较。我也尝试用物种替代栖息地来进行基于物种的比较。在这两种情况下,我得到了以下错误:

Error in contrMat(table(mf[[nm]]), type = types[pm]) : less than two groups

我不知道怎样才能有不到两个组的任何地方,因为群体的最小数目我有任何列是三。也许我误解了R所指的是什么“组”?有什么可能出错的建议以及如何解决它?

编辑根据建议(谢谢),这里是一些关于我的数据的更多信息,可能会使它具有可重现性。首先,我整个使用的代码。

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 

摘要(glht(模型,linfct = MCP(境= “杜克”)))

Error in contrMat(table(mf[[nm]]), type = types[pm]) : 
    less than two groups 

以下是有关可能再现它有助于数据的一些汇总......它只是有点大,我不知道如何使它最小化,但仍然得到错误。

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") 

你有它:感谢任何指导。

+1

这不是[再现的](http://stackoverflow.com/questions/5963269/how-to-make-a-great- r-reproducible-example),因为在你的例子中'habitat'或'site'没有变化。 – 2013-02-12 21:10:42

+0

您确实需要提供有关您的数据的更多信息。看看这里的其他问题,了解如何使您的产品具有可重复性。至少要给你的数据str()和summary()。我制作了一套与您的描述相匹配的虚拟集合,并且运行良好,所以我期望有关于您的数据的信息会导致此问题。 – alexwhan 2013-02-13 01:47:58

+0

这不是我。这就是说:我也无法重现它('score'不包含在你的'dput()'中)。制作一个可重复使用的例子通常非常困难,但这是唯一的出路(您是否阅读了@ user1317221_G的评论中的链接?)您可以在某处发布'baci_null_red.csv'数据集或其匿名/随机版本的数据集可公开访问?你可以检查'桌子(物种,栖息地)',并确保你有每个栖息地的物种? – 2013-03-11 04:04:45

回答

0

我刚刚通过避免附加数据框来解决此问题。我在网上阅读,当试图将一个合适的模型传递给一个函数时,附加数据框会产生问题。

我建议使用这样的代码,而不是使用attach()

summary(model<-lme(fix~treat*week, random=~1|rep/week, data=dataframename)) 
summary(glht(model, linfct=mcp(treat="Tukey")))