Как изменить порядок образцов в тесте Тьюки в R?
Проблема: я хотел бы узнать, как я могу изменить порядок выборок, для которых тест Тьюки в R вычисляет средние значения и назначает соответствующие буквы. Очень простой пример ниже.
Я играл с данными радужной оболочки и обнаружил, что есть различия в Sepal.Length среди разных видов. Вот коробочный сюжет:
Я провел тест ANOVA и обнаружил, что различия статистически значимы.
> fit <- lm(Sepal.Length ~ Species, data = iris)
> summary(aov(fit))
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 63.21 31.606 119.3 <2e-16 ***
Residuals 147 38.96 0.265
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Затем я провел тест Тьюки и получил следующее:
> library(agricolae)
> HSD.test(fit, "Species", group=T, console=T)
Study: fit ~ "Species"
HSD Test for Sepal.Length
Mean Square Error: 0.2650082
Species, means
Sepal.Length std r Min Max
setosa 5.006 0.3524897 50 4.3 5.8
versicolor 5.936 0.5161711 50 4.9 7.0
virginica 6.588 0.6358796 50 4.9 7.9
alpha: 0.05 ; Df Error: 147
Critical Value of Studentized Range: 3.348424
Honestly Significant Difference: 0.2437727
Means with the same letter are not significantly different.
Groups, Treatments and means
a virginica 6.588
b versicolor 5.936
c setosa 5.006
В соответствии с таблицей групп функция HSD.test сортирует средства в порядке убывания, а затем присваивает буквы. Таким образом, "virginica" имеют наибольшее среднее значение, поэтому оно является первым в таблице.
Вопросы: Есть ли способ изменить сортировку и назначение букв по умолчанию? Могу ли я отсортировать образцы в порядке возрастания средств, а затем назначить буквы. Ожидаемый результат следующий:
a setosa 5.006
b versicolor 5.936
c virginica 6.588
Возможное решение: в пакете multcomp есть две функции, которые могут сделать это, работая вместе:
1 - glht сделать тест Тьюки
> an <- aov(fit)
> library(multcomp)
> glht(an, linfct = mcp(Species = "Tukey"))
General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Linear Hypotheses:
Estimate
versicolor - setosa == 0 0.930
virginica - setosa == 0 1.582
virginica - versicolor == 0 0.652
2 - cld можете предоставить мне письма, назначенные Species в соответствии с уровнями фактора iris$Species
> cld(glht(an, linfct = mcp(Species = "Tukey")))
setosa versicolor virginica
"a" "b" "c"
К несчастью, glht Функция не отображает другие данные, которые могут быть полезны и необходимы для создания барплота (средние значения, std, p-значения). Конечно, я могу сделать это отдельно с другими специальными функциями, или просто использовать оба HSD.test а также cld, Но я бы предпочел решить проблему с сортировкой средств в HSD.test функционировать и использовать только этот.
3 ответа
Я заметил, что уже поздно отвечать на этот вопрос. Однако я столкнулся с точно такой же проблемой и хотел бы поделиться своим решением в будущем. Надеюсь, это когда-нибудь кому-нибудь поможет.
первый вариант
Можно использовать multcompLetters() например, с результатами из TukeyHSD(), Однако это не позволяет произвольно упорядочить результат и не так просто использовать.
второй вариант
Так как мне нужен был произвольный порядок, я написал свою собственную функцию, которая принимает вектор букв, возвращаемых из HSD.test и меняет буквы таким образом, чтобы результат был хорошим. Значение букв в алфавите появляется первым.
library(agricolae)
reorder<-function(inV){
collapsed <- paste(inV,sep="",collapse = "")
u <- unique(strsplit(collapsed,"")[[1]])
if(length(u)<2){
return(inV)
}
u <- u[order(u)]
m <- matrix(nrow=NROW(inV),ncol=length(u))
m[]<-F
for(i in 1:length(inV)){
s <- strsplit(inV[i],"")[[1]]
index <- match(s,u)
m[i,index] <- T
}
for(i in 1:(length(u)-1)){
firstColT <- match(T,m[,i])[1] #first row with true in current column
firstT <- match(T,rowSums(m[,i:length(u)] > 0))[1] #first row with true in rest
if(firstT < firstColT){
colT <- match(T,m[firstT,i:length(u)])[1]
colT <- colT + i - 1 #correct index for leftout columns in match
tmp <- m[,colT]
m[,colT] <- m[,i]
m[,i] <- tmp
}
}
res <- vector(mode = "character", length=length(trt))
for(i in 1:length(inV)){
l <- u[m[i,]]
res[i] <- paste(l,sep="",collapse = "")
}
return(res)
}
fit <- lm(Sepal.Length ~ Species, data = iris)
a <- HSD.test(fit, "Species", group=T, console=F)$groups
a <- a[rev(rownames(a)),] #order the result the way you want
a$M <- reorder(as.character(a$M))
Для примера это немного излишне, но оно должно работать и для более сложных случаев.
Также возможно решить с помощью multcompLetters() и TukeyHSD(). Вы должны изменить параметр «обратный»
library(multcompView)
fit <- aov(Sepal.Length ~ Species, data = iris)
tukey<-TukeyHSD(fit, ordered = T)
tukey_1<-multcompLetters2(Sepal.Length ~ Species,
tukey$Species[,"p adj"],
iris,reversed = T)
tukey_2<-multcompLetters2(Sepal.Length ~ Species,
tukey$Species[,"p adj"],
iris,reversed = F)
tukey_1
tukey_2
tapply(iris$Sepal.Length, iris$Species, mean)
Прежде всего, спасибо за функцию. Это было то, что я искал. Но я думаю, что есть ошибка в
res <- vector(mode = "character", length=length(trt)),
так должно быть
res <- vector(mode = "character", length=length("trt"))