Rについて質問するのはこれが初めてなので、必要な詳細を省略していることをお許しください。私はできるだけ徹底的に、しかしできるだけ簡潔にしようとします。
私のデータには 4 つの列があります。スコア (-1、0、または 1 の値)。種（少なくとも10種）; 生息地（3つの異なるもの）; およびサイト (少なくとも 8 つの異なるもの)。したがって、各行は、スコア、種、生息地、およびサイトの一意の組み合わせです。上の行については、以下を参照してください。
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")))
私が知る限り、生息地のタイプに基づいてスコアのペアワイズ比較を行います。また、種に基づく比較を行うために、生息地を種に置き換えてみました。どちらの場合も、次のエラーが発生しました。
contrMat(table(mf[[nm]])、type = types[pm]) のエラー: 2 つ未満のグループ
任意の列のグループの最小数は 3 であるため、どこでも 2 つ未満のグループを持つ方法がわかりません。Rが参照している「グループ」を誤解しているのかもしれません。何がうまくいかないのか、それを修正する方法について何か提案はありますか?
EDIT推奨事項 (ありがとう) に従って、再現可能になる可能性のある私のデータに関する詳細情報を次に示します。まず、全体で使用したコード。
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(生息地="Tukey")))
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")
これで終わりです。ご指導いただきありがとうございます。