CODEWITHSUNDEEP
rstatistics
Docconcoct
Docconcoct

Reputation: 2050

ezANOVA and pairwise.t.test in R: output

I've been running ezANOVA using R followed by pairwise comparisons (pairwise.t.test...)

There seems to be output not included(?) once run so I have the following questions in relation to this:

Is it possible to generate Mean Square Error in ezANOVA for repeated measures and if so how?

And also how does one generate output for df, t and SE in pairwise.t.test?

This may be very obvious to some but I've consulted books and online resources and come up empty.

Thanks as always for your time. It is very much appreciated.

Update:

Calculating the MSE in ezANOVA output:

: > Exp2DModel
$ANOVA
        Effect DFn DFd        SSn       SSd          F            p p<.05       ges
1  (Intercept)   1  24 159.231492 3.1947477 1196.19955 5.453251e-22     * 0.9468257
2       Visual   2  48   1.603228 1.9257323   19.98069 4.861610e-07     * 0.1520259
3        Audio   4  96   1.831761 0.9966656   44.10934 5.788074e-21     * 0.1700122
4 Visual:Audio   8 192  14.230991 2.8253824  120.88409 1.010389e-70     * 0.6141057

: > Exp2DModel$ANOVA$SSd[2]/Exp2DModel$ANOVA$DFd[2]

MSE = [1] 0.04011942

Upvotes: 2

Views: 5656

Answers (2)

Marcus Morrisey
Marcus Morrisey

Reputation: 166

In response to your second question, how to extract t-values from pairwise.t.test. It occurred to me that having the degrees of freedom from your data, and the precise p-values would allow you to use the qt() function to determine the t-values. Example:

myTTests<-pairwise.t.test(vector1,vector2
                ,p.adjust.method = "none"
                )
qt(myTTests[[3]],degreesOfFreedom)

Upvotes: 2

Jota
Jota

Reputation: 17621

Maybe this information from the help information for ezANOVA will be useful.

library(ez)
data(ANT)
rt_anova = ezANOVA(
    data = ANT[ANT$error==0,]
    , dv = rt
    , wid = subnum
    , within = .(cue,flank)
    , between = group
)

#Show the ANOVA and assumption tests.
print(rt_anova)

See this link for more examples.

The help page indicates the following Values are part of the output object:

DFn          Degrees of Freedom in the numerator (a.k.a. DFeffect).
DFd          Degrees of Freedom in the denominator (a.k.a. DFerror).
SSn          Sum of Squares in the numerator (a.k.a. SSeffect).
SSd          Sum of Squares in the denominator (a.k.a. SSerror).
F            F-value.
p            p-value (probability of the data given the null hypothesis).
p<.05        Highlights p-values less than the traditional alpha level of .05.
ges          Generalized Eta-Squared measure of effect size (see in references below: Bakeman, 2005).
GGe          Greenhouse-Geisser epsilon.
p[GGe]       p-value after correction using Greenhouse-Geisser epsilon.
p[GGe]<.05   Highlights p-values (after correction using Greenhouse-Geisser epsilon) less than the traditional alpha level of .05.
HFe          Huynh-Feldt epsilon.
p[HFe]       p-value after correction using Huynh-Feldt epsilon.
p[HFe]<.05   Highlights p-values (after correction using Huynh-Feldt epsilon) less than the traditional alpha level of .05.
W            Mauchly's W statistic

These columns should contain most of the raw material needed for calculating MSE... Though, I don't actually know how to calculate MSE in a repeated measures ANOVA...

As far as I know pairwise.t.test provides output that contains only p-values. From the help page:

attach(airquality)
Month <- factor(Month, labels = month.abb[5:9])
pairwise.t.test(Ozone, Month)

I did find this post on extracting t-values from pairwise.t.test, which may help. Essentially, in the post, the code for pairwise.t.test is modified. The line that calculates p-values is commented out (i.e. # 2 * pt(-abs(t.val), total.degf) and they add a line to return the t.val instead of p-values. You can see in that code that degrees of freedom are calculated and an object called se.diff, thus, you should be able to return them, too, by modifying the code.

Using the modified paired.t.test from the post I linked, I added a little bit to return total.degf. I'll leave it to you to include se.dif, if you want. Note that this code is for example purposes and needs some cleaning (e.g. the code recognizes the t.val values as being p-values. This wouldn't be hard to change.

my.pairded.t.test <- function (x, g, p.adjust.method = 
p.adjust.methods, pool.sd = TRUE,
...)
{
DNAME <- paste(deparse(substitute(x)), "and", 
deparse(substitute(g)))
g <- factor(g)
p.adjust.method <- match.arg(p.adjust.method)
if (pool.sd) {
    METHOD <- "t tests with pooled SD"
    xbar <- tapply(x, g, mean, na.rm = TRUE)
    s <- tapply(x, g, sd, na.rm = TRUE)
    n <- tapply(!is.na(x), g, sum)
    degf <- n - 1
    total.degf <- sum(degf)
    pooled.sd <- sqrt(sum(s^2 * degf)/total.degf)
    compare.levels <- function(i, j) {
        dif <- xbar[i] - xbar[j]
        se.dif <- pooled.sd * sqrt(1/n[i] + 1/n[j])
        t.val <- dif/se.dif
        #  2 * pt(-abs(t.val), total.degf)  this is commented out
        t.val    # this is added
    }
}
else {
    METHOD <- "t tests with non-pooled SD"
    compare.levels <- function(i, j) {
        xi <- x[as.integer(g) == i]
        xj <- x[as.integer(g) == j]
        t.test(xi, xj, ...)$statistic       # this is changed in case 
                            # pool.sd=F
    }
}
PVAL <- pairwise.table(compare.levels, levels(g), 
p.adjust.method)
ans <- list(method = METHOD, data.name = DNAME, p.value 
= PVAL,
    p.adjust.method = p.adjust.method, total.degf=total.degf) # I added the total.degf part
class(ans) <- "pairwise.htest"
ans
}

Put it to use.

a<-my.pairded.t.test(Ozone, Month)
a$total.degf

Upvotes: 1

Related Questions