#r #anova #posthoc
#r #anova #posthoc
Вопрос:
Я пытаюсь понять, как работать с ANOVA и post-hoc тестами в R. До сих пор я использовал aov () и TukeyHSD() для анализа своих данных. Пример:
uni2.anova <- aov(Sum_Uni ~ Micro, data= uni2)
uni2.anova
Call:
aov(formula = Sum_Uni ~ Micro, data = uni2)
Terms:
Micro Residuals
Sum of Squares 0.04917262 0.00602925
Deg. of Freedom 15 48
Residual standard error: 0.01120756
Estimated effects may be unbalanced
Моя проблема в том, что теперь у меня есть огромный список попарных сравнений, но я ничего не могу с этим поделать:
TukeyHSD(uni2.anova)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Sum_Uni ~ Micro, data = uni2)
$Micro
diff lwr upr p adj
Act_Glu2-Act_Ala2 -0.0180017863 -0.046632157 0.0106285840 0.6448524
Ana_Ala2-Act_Ala2 -0.0250134285 -0.053643799 0.0036169417 0.1493629
NegI_Ala2-Act_Ala2 0.0702274527 0.041597082 0.0988578230 0.0000000
Этот набор данных содержит 40 строк…
В идеале я хотел бы получить набор данных, который выглядит примерно так:
- Act_Glu2: a
- Act_Ala2: a
- NegI_Ala2: b…
Надеюсь, вы поняли суть. До сих пор я не нашел ничего сопоставимого в Интернете… Я также попытался выбрать только значимые пары в файле, полученном из TukeyHSD, но файл не «подтверждает», что он состоит из строк и столбцов, что делает выбор невозможным…
Может быть, в моем подходе есть что-то принципиально неправильное?
Комментарии:
1. Что означает «Act_Glu2: a»? Чем это отличается от «Act_Glu2-Act_Ala2»
2. @John Ооо, возможно, мы ошибаемся. OP упоминает «классифицировать» в названии, но нигде в сообщении. Если она действительно хочет классифицировать (кластер?) тогда она может писать это, чтобы показать, что ей нужен список аминокислот и кластер, к которому они были отнесены (т.Е. Act_Glu2 и Act_Ala2 находятся в кластере «a»). Я не знаю, хотя я могу быть совершенно неправ. Во всяком случае, Кэролин, ты можешь прояснить некоторые моменты?
3. @ John Colby: Да, я думаю, вы понимаете, что я имею в виду. Act_Glu2 и Act_Ala2 не показывают существенной разницы в тесте Тьюки, следовательно, они будут классифицированы (или сгруппированы, если это правильный термин) в одну группу. NegI_Ala существенно отличается по крайней мере от одного из них, поэтому, если я построю график данных, я бы показал это значение, добавив «a» к первым двум и «b» к третьей точке данных. Но, поскольку существует так много наборов данных, я бы предпочел не делать это вручную…
Ответ №1:
Я думаю, что OP хочет, чтобы буквы отображали сравнения.
library(multcompView)
multcompLetters(extract_p(TukeyHSD(uni2.anova)))
Это даст вам буквы.
Вы также можете использовать пакет multcomp
library(multcomp)
cld(glht(uni2.anova, linct = mcp(Micro = "Tukey")))
Надеюсь, это то, что вам нужно.
Комментарии:
1. С разъяснением Кэролин, я думаю, что это правильный путь.
2. Отлично! Это именно то, что я искал! Большое вам спасибо.
3.С небольшой поправкой 🙂
hsd <- TukeyHSD(uni2.anova)multcompLetters(extract_p(hsd$Micro))Потому что TukeyHSD(uni2.anova) приводит к большему, чем просто список попарных сравнений amp; в этом случае hsd $ Micro — это просто список.
Ответ №2:
Результаты из TukeyHSD представляют собой список. Используйте str , чтобы посмотреть на структуру. В вашем случае вы увидите, что это список из одного элемента, и этот элемент в основном является матрицей. Итак, чтобы извлечь первый столбец, вам нужно сохранить результат TukeyHSD
hsd <- TukeyHSD(uni2.anova)
Если вы посмотрите на str(hsd) то, что вы можете, то сможете получить биты…
hsd$Micro[,1]
Это даст вам столбец ваших различий. Вы должны быть в состоянии извлечь то, что хотите сейчас.
Комментарии:
1. О, здорово! Я пробовал что-то вроде:
TukeyHSD(uni2.anova)[4,]который вернул «Неправильное количество измерений»… Спасибо!2. Хм, если я хочу выбрать строки с определенными атрибутами, подобными этому:
hsd$Micro[hsd$Micro[,4] < 0.05]я не получаю все столбцы hsd $ Micro, только 4-й.3. Исправлено!
hsd$Micro[hsd$Micro[,4] < 0.05,]
Ответ №3:
Трудно сказать без примеров данных, но предполагается Micro , что это всего лишь фактор с 4 уровнями и uni2 выглядит примерно так
n = 40
Micro = c('Act_Glu2', 'Act_Ala2', 'Ana_Ala2', 'NegI_Ala2')[sample(4, 40, rep=T)]
Sum_Uni = rnorm(n, 5, 0.5)
Sum_Uni[Micro=='Act_Glu2'] = Sum_Uni[Micro=='Act_Glu2'] 0.5
uni2 = data.frame(Sum_Uni, Micro)
> uni2
Sum_Uni Micro
1 4.964061 Ana_Ala2
2 4.807680 Ana_Ala2
3 4.643279 NegI_Ala2
4 4.793383 Act_Ala2
5 5.307951 NegI_Ala2
6 5.171687 Act_Glu2
...
тогда я думаю, что вы на самом деле пытаетесь получить базовый результат множественной регрессии:
fit = lm(Sum_Uni ~ Micro, data = uni2)
summary(fit)
anova(fit)
> summary(fit)
Call:
lm(formula = Sum_Uni ~ Micro, data = uni2)
Residuals:
Min 1Q Median 3Q Max
-1.26301 -0.35337 -0.04991 0.29544 1.07887
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.8364 0.1659 29.157 < 2e-16 ***
MicroAct_Glu2 0.9542 0.2623 3.638 0.000854 ***
MicroAna_Ala2 0.1844 0.2194 0.841 0.406143
MicroNegI_Ala2 0.1937 0.2158 0.898 0.375239
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4976 on 36 degrees of freedom
Multiple R-squared: 0.2891, Adjusted R-squared: 0.2299
F-statistic: 4.88 on 3 and 36 DF, p-value: 0.005996
> anova(fit)
Analysis of Variance Table
Response: Sum_Uni
Df Sum Sq Mean Sq F value Pr(>F)
Micro 3 3.6254 1.20847 4.8801 0.005996 **
Residuals 36 8.9148 0.24763
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Вы можете получить доступ к числам в любой из этих таблиц, например,
> summary(fit)$coef[2,4]
[1] 0.0008536287
Чтобы просмотреть список того, что хранится в каждом объекте, используйте names() :
> names(summary(fit))
[1] "call" "terms" "residuals" "coefficients"
[5] "aliased" "sigma" "df" "r.squared"
[9] "adj.r.squared" "fstatistic" "cov.unscaled"
В дополнение к TukeyHSD() функции, которую вы нашли, есть много других вариантов для дальнейшего просмотра парных тестов и исправления p-значений, если это необходимо. К ним относятся pairwise.table() , estimable() в gmodels , resampling boot пакеты и и другие…