CODEWITHSUNDEEP
ranova
Caeline
Caeline

Reputation: 119

R post hoc comparisons of ezANOVA

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

Upvotes: 1

Views: 4901

Answers (1)

Adamm
Adamm

Reputation: 2306

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.

Upvotes: 2

Related Questions