R апостериорных сравнений ezANOVA

Я выполняю следующие ezANOVA:

RMANOVAGHB1 <- ezANOVA(GHB1, dv=DIF.SCORE.STARTLE, wid=RAT.ID, within=TRIAL.TYPE, between=GROUP, detailed = TRUE, return_aov = TRUE)

Мой набор данных выглядит так:

RAT.ID DIF.SCORE.STARTLE GROUP TRIAL.TYPE
1       1            170.73   SAL       TONO
2       1             80.07   SAL       NOAL
3       2            456.40  PROP       TONO
4       2            290.40  PROP       NOAL
5       3            507.20   SAL       TONO
6       3            261.60   SAL       NOAL
7       4            208.67  PROP       TONO
8       4            137.60  PROP       NOAL
9       5            500.50   SAL       TONO
10      5            445.73   SAL       NOAL

до rat.id 16.

Мои руководители не работают с R, поэтому они не могут мне помочь. Мне нужен код, который даст мне все постфактум контрасты, но поиск его только сбивает меня с толку все больше и больше. Я уже пытался выполнить TukeyHSD на выводе aov ezANOVA, а затем попробовал pairwise.t.test (поскольку я обнаружил, что bonferroni является более подходящим исправлением в этом случае), но, похоже, ни один из них не работает. Я также нашел информацию об использовании линейной модели, а затем multcomp, но я не знаю, будет ли это хорошим решением в этом случае. Я чувствую, что проблема со всем, что я пробовал, заключается либо в том, что у меня есть переменные между и внутри, либо в том, что мой набор данных настроен неправильно. Еще одним усложняющим фактором является то, что я только начинаю работать с R, и мои статистические знания все еще довольно базовые, поскольку это один из моих первых практических опытов проведения анализа.

Если это важно, вот результат анова:

$ANOVA
            Effect DFn DFd       SSn       SSd         F           p p<.05         ges
1      (Intercept)   1  14 1233568.9 1076460.9 16.043280 0.001302172     * 0.508451750
2            GROUP   1  14  212967.9 1076460.9  2.769771 0.118272657       0.151521743
3       TRIAL.TYPE   1  14  137480.6  116097.9 16.578499 0.001143728     * 0.103365833
4 GROUP:TRIAL.TYPE   1  14   11007.2  116097.9  1.327335 0.268574391       0.009145489

$aov

Call:
aov(formula = formula(aov_formula), data = data)

Grand Mean: 196.3391

Stratum 1: RAT.ID

Terms:
                    GROUP Residuals
Sum of Squares   212967.9 1076460.9
Deg. of Freedom         1        14

Residual standard error: 277.2906
1 out of 2 effects not estimable
Estimated effects are balanced

Stratum 2: RAT.ID:TRIAL.TYPE

Terms:
                TRIAL.TYPE GROUP:TRIAL.TYPE Residuals
Sum of Squares    137480.6          11007.2  116097.9
Deg. of Freedom          1                1        14

Residual standard error: 91.0643
Estimated effects may be unbalanced

person Caeline    schedule 04.09.2017    source источник


Ответы (1)


Мое решение, учитывая ваш набор данных - первые 5 крыс:

<сильный>1. Давайте построим линейную модель:

model.lm = lm(DIF_SCORE_STARTLE ~ GROUP * TRIAL_TYPE, data = dat)

<сильный>2. Давайте проверим однородность дисперсии (leveneTest) и распределения нашей модели (Shapiro-Wilk). Мы ищем нормальное распределение, и наша дисперсия должна быть однородной. Два теста для этого:

>shapiro.test(resid(model.lm))

    Shapiro-Wilk normality test

data:  resid(model.lm)
W = 0.91783, p-value = 0.3392

> leveneTest(DIF_SCORE_STARTLE ~ GROUP * TRIAL_TYPE, data = dat)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3   0.066  0.976
       6 

Наши значения p выше 0,05 в обоих случаях, поэтому у нас нет доказательств того, что наша дисперсия различается между группами. В случае теста на нормальность мы также можем сделать вывод, что выборка не отклоняется от нормальности. Подводя итог, мы можем использовать параметрические тесты, такие как ANOVA или попарный t-тест.

3. Вы также можете бегать:

hist(resid(model.lm))

Чтобы проверить, как выглядит распределение наших данных. И проверьте модель:

plot(model.lm)

Здесь: https://stats.stackexchange.com/questions/58141/interpreting-plot-lm/65864 вы найдете интерпретацию графиков, созданных этой функцией. Как я видел, данные выглядят нормально.

4.Теперь, наконец, мы можем провести тест ANOVA и апостериорный тест HSD:

> anova(model.lm)
Analysis of Variance Table

Response: DIF_SCORE_STARTLE
                 Df Sum Sq Mean Sq F value Pr(>F)
GROUP             1   7095    7095  0.2323 0.6469
TRIAL_TYPE        1  39451   39451  1.2920 0.2990
GROUP:TRIAL_TYPE  1     84      84  0.0027 0.9600
Residuals         6 183215   30536
> (result.hsd = HSD.test(model.lm, list('GROUP', 'TRIAL_TYPE')))
$statistics
    Mean       CV  MSerror      HSD r.harmonic
  305.89 57.12684 30535.91 552.2118        2.4

$parameters
  Df ntr StudentizedRange alpha  test           name.t
   6   4         4.895599  0.05 Tukey GROUP:TRIAL_TYPE

$means
          DIF_SCORE_STARTLE      std r    Min    Max
PROP:NOAL          214.0000 108.0459 2 137.60 290.40
PROP:TONO          332.5350 175.1716 2 208.67 456.40
SAL:NOAL           262.4667 182.8315 3  80.07 445.73
SAL:TONO           392.8100 192.3561 3 170.73 507.20

$comparison
NULL

$groups
        trt    means M
1 SAL:TONO  392.8100 a
2 PROP:TONO 332.5350 a
3 SAL:NOAL  262.4667 a
4 PROP:NOAL 214.0000 a

Как видите, наши «пары» сгруппированы в одну большую группу a, что означает, что между ними нет существенной разницы. Однако есть некоторая разница между NOAL и TONO независимо от SAL и PROP.

person Adamm    schedule 04.09.2017
comment
Спасибо за такое ясное объяснение (спасибо, что нашли время)! Однако мне было интересно, это дает мне разные значения p, и значения p, которые я опубликовал, являются правильными, согласно моим руководителям ... То, что вы опубликовали, учитывает, что я повторил измерения? - person Caeline; 05.09.2017
comment
Я использовал то, что вы дали (данные), поэтому наши результаты разные. В случае p-значений в ANOVA (дисперсионном анализе) у нас есть значение F, и оно должно быть ниже 0,05 - это означает, что одна из групп значительно отличается от другой на основе дисперсии. В моих расчетах выше мы можем заключить, что нет существенной разницы между крысами (SAL против PROP) или (TONO против NOAL), но если мы проверим некоторые комбинации, например. (SAL и TONO против SAL и NOAL) есть некоторая разница, поэтому я провел апостериорный тест, чтобы выяснить, какая комбинация отличается от другой. - person Adamm; 05.09.2017
comment
Я использовал ваш код для всего набора данных, что дает мне результаты, отличные от ezANOVA (например, группа оказывает значительное влияние в соответствии с дисперсионным анализом на модель lm, но согласно ezANOVA это не так). По словам моего руководителя, результаты, которые дает мне ezANOVA, являются правильными. Я думаю, что анова в модели lm не учитывает, что каждая крыса/субъект измеряется дважды, и поэтому мне нужен повторный анова, а не обычный. Или, может быть, он не распознает, что группа является переменной между субъектами, в то время как Trial.type является переменной внутри субъектов. - person Caeline; 05.09.2017
comment
Да, насколько я понимаю, ezANOVA предназначен для данных с репликациями. Вы ставите флаг: return_aov = TRUE, чтобы вы могли выполнить постфактум тест. Может что-то вроде TukeyHSD(returned_aov) подойдет? - person Adamm; 05.09.2017
comment
К сожалению, это не так. Пытался сделать это на ANOVA, и aov вывода дал мне эти ошибки: no applicable method for 'TukeyHSD' applied to an object of class "data.frame" и no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')". - person Caeline; 07.09.2017
comment
Хм, может быть немного сложно сделать некоторые постфактум с анова с репликами. TukeyHSD предназначен для независимых данных. Однако я нашел и ответил здесь: ссылка - person Adamm; 08.09.2017
comment
LMGHB1 <- lme(DIF.SCORE.STARTLE ~ TRIAL.TYPE * GROUP, data = GHB1, random = ~1|RAT.ID/TRIAL.TYPE). Это дает мне те же правильные значения F и p, что и ezANOVA. Только теперь я застрял с использованием multcomp и glht. Я использую summary(glht(LMGHB1, linfct=mcp(TRIAL.TYPE = "Tukey")), test = adjusted(type = "bonferroni")), но это дает мне только TONO - NOAL или, если я добавляю GROUP, это дает мне SAL - PROP, но я хочу видеть сравнения, такие как: TONO SAL – TONO PROP; TONO SAL – NOAL SAL; TONO SAL – NOAL PROP; TONO PROP – NOAL SAL; ... Спасибо за всю вашу помощь, без этого не обошлось. ! - person Caeline; 13.09.2017