1

I perform following ezANOVA:

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

My dataset looks like this: 

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

up until rat.id 16.

My supervisors don't work with R, so they can't help me. I need code that will give me all post hoc contrasts, but looking it up only confuses me more and more. I already tried to do TukeyHSD on the aov output of ezANOVA and tried pairwise.t.test next (as I found out bonferroni is a more appropriate correction in this case), but none seem to work. I've also found things about using a linear model and then multcomp, but I don't know if that would be a good solution in this case. I feel like the problem with everything I tried is either that I have between and within variables or that my dataset is not set up right. Another complicating factor is that I'm just a beginner with R and my statistical knowledge is still pretty basic as this is one of my first practical experiences with doing analyses.

If it's important, this is the output of the anova: 

$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
Caeline
  • 119
  • 3
  • 13

1 Answers1

2

My solution, considering your dataset - first 5 rats:

1. Let's build the linear model:

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

2. Let's chceck the homogeneity of variance (leveneTest) and distribution of our model (Shapiro-Wilk). We are looking for normal distribution and our variance should be homogenic. Two tests for this:

>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

Our p-values are higher than 0.05 in both cases so we don't have proof that our variance differs between groups. In case of normality test we can also conclude that the sample doesn't deviate from normality. Summarizing we can use parametrical tests such as ANOVA or pairwise t-test.

3.Yo can also run:

hist(resid(model.lm))

To check how does distribution of our data look like. And check the model:

plot(model.lm)

Here: https://stats.stackexchange.com/questions/58141/interpreting-plot-lm/65864 you'll find interpretation of plots produced by this function. As I saw, data looks fine.

4.Now finally we can do ANOVA test and post hoc HSD test:

> 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

As you see, our 'pairs' have been grouped in one big group a that means that there are not significant difference between them. However there's some difference between NOAL and TONO no matter of SAL and PROP.

Adamm
  • 2,150
  • 22
  • 30
  • Thank you for such a clear explanation (thank you for taking the time)! I was wondering though, it gives me different p values, and the p values that I posted are the correct ones according to my supervisors... Is what you posted keeping into account that I have repeated measurements? – Caeline Sep 05 '17 at 06:58
  • I used what you gave (data) that's why our results are different. In case of p-values in ANOVA (analysis of variance) we have F value and it should be below 0.05 - that means one of groups significantly differ from other based of variance. In my calculations above we can conclude that there's not significant difference between rats (SAL vs PROP) or (TONO vs NOAL) but if we chceck some combinations ex. (SAL and TONO vs SAL and NOAL) there is some difference, hence I performed post hoc test to find out which combination differs from other. – Adamm Sep 05 '17 at 07:12
  • I used your code on the whole dataset, which gives me different results than the ezANOVA (for example Group has a significant effect according to the anova on the lm model, but according to ezANOVA it doesn't). According to my supervisor, the results that ezANOVA gives me, are the right results. I think that the anova on the lm model doesn't take into account that every rat/subject is measured twice and that I thus need a repeated measurements anova and not a normal one. Or maybe it doesn't recognize that group is a between-subjects variable, while trial.type is a within-subjects variable. – Caeline Sep 05 '17 at 09:04
  • Yes, as I see ezANOVA is designed for data with replications. You put flag: `return_aov = TRUE` so you are able to do some post-hoc test. Maybe something like `TukeyHSD(returned_aov)` will work? – Adamm Sep 05 '17 at 13:57
  • Sadly, it doesn't. Tried to do it on the ANOVA and the aov of the output giving me these errors: `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')"`. – Caeline Sep 07 '17 at 14:08
  • Hmm it might be a little difficult to do some post hoc with anova with replicates TukeyHSD is designed for independent data. However I found and answer here: [link](https://stackoverflow.com/questions/14707725/post-hoc-tests-with-ezanova-output) – Adamm Sep 08 '17 at 05:59
  • `LMGHB1 <- lme(DIF.SCORE.STARTLE ~ TRIAL.TYPE * GROUP, data = GHB1, random = ~1|RAT.ID/TRIAL.TYPE)`. This gives me the same, correct F and p values as ezANOVA. Only now I'm stuck with using multcomp and glht. I'm using `summary(glht(LMGHB1, linfct=mcp(TRIAL.TYPE = "Tukey")), test = adjusted(type = "bonferroni"))` but this gives only gives me TONO - NOAL or if I put in GROUP, it gives me SAL - PROP, but I want to see comparisons like: `TONO SAL – TONO PROP; TONO SAL – NOAL SAL; TONO SAL – NOAL PROP; TONO PROP – NOAL SAL; ...` Thanks for all your help, couldn't have done without it! – Caeline Sep 13 '17 at 17:31