[R] extracting p-values from Anova objects (from the car library)
David Winsemius
dwinsemius at comcast.net
Tue Aug 24 00:19:42 CEST 2010
On Aug 23, 2010, at 5:35 PM, Johan Steen wrote:
> Thanks for your replies,
>
> but unfortunately none of them seem to help.
> I do get p-values in the output, but can't seem to locate them
> anywhere in these objects via the str() function.
That is because in the case of Anova.mlm objects ... they aren't there.
> I also get very different output using str() than you obtained from
> the lm help page
>
> Here's my output:
>
> > A <- factor( rep(1:2,each=3) )
> > B <- factor( rep(1:3,times=2) )
> > idata <- data.frame(A,B)
> > idata
> A B
> 1 1 1
> 2 1 2
> 3 1 3
> 4 2 1
> 5 2 2
> 6 2 3
> >
> > fit <- lm( cbind(a1_b1,a1_b2,a1_b3,a2_b1,a2_b2,a2_b3) ~ sex,
> data=Data.wide)
> > result <- Anova(fit, type="III", test="Wilks", idata=idata,
> idesign=~A*B)
> > result
>
> Type III Repeated Measures MANOVA Tests: Wilks test statistic
> Df test stat approx F num Df den Df Pr(>F)
> (Intercept) 1 0.02863 610.81 1 18 2.425e-15
> sex 1 0.76040 5.67 1 18 0.02849
> A 1 0.91390 1.70 1 18 0.20925
> sex:A 1 0.99998 0.00 1 18 0.98536
> B 1 0.26946 23.05 2 17 1.443e-05
> sex:B 1 0.98394 0.14 2 17 0.87140
> A:B 1 0.27478 22.43 2 17 1.704e-05
> sex:A:B 1 0.98428 0.14 2 17 0.87397
Remember I suggested looking for print methods? methods(print)
produces a print.Anova.mlm citation. The above output is produced by
that function which is not visible unless you do:
> getAnywhere(print.Anova.mlm)
A single object matching ‘print.Anova.mlm’ was found
It was found in the following places
registered S3 method for print from namespace car
namespace:car
with value
function (x, ...)
{
test <- x$test
repeated <- x$repeated
ntests <- length(x$terms)
tests <- matrix(NA, ntests, 4)
if (!repeated)
SSPE.qr <- qr(x$SSPE)
for (term in 1:ntests) {
eigs <- Re(eigen(qr.coef(if (repeated) qr(x$SSPE[[term]])
else SSPE.qr,
x$SSP[[term]]), symmetric = FALSE)$values)
tests[term, 1:4] <- switch(test, Pillai = stats:::Pillai(eigs,
x$df[term], x$error.df), Wilks = stats:::Wilks(eigs,
x$df[term], x$error.df), `Hotelling-Lawley` =
stats:::HL(eigs,
x$df[term], x$error.df), Roy = stats:::Roy(eigs,
x$df[term], x$error.df))
}
ok <- tests[, 2] >= 0 & tests[, 3] > 0 & tests[, 4] > 0
ok <- !is.na(ok) & ok
tests <- cbind(x$df, tests, pf(tests[ok, 2], tests[ok, 3],
tests[ok, 4], lower.tail = FALSE))
rownames(tests) <- x$terms
colnames(tests) <- c("Df", "test stat", "approx F", "num Df",
"den Df", "Pr(>F)")
tests <- structure(as.data.frame(tests), heading = paste("\nType ",
x$type, if (repeated)
" Repeated Measures", " MANOVA Tests: ", test, " test
statistic",
sep = ""), class = c("anova", "data.frame"))
print(tests)
invisible(x)
}
<environment: namespace:car>
Notice that after printing "tests", those results are basically
thrown away, and the function invisibly returns its argument. Hacking
that function would be simple matter if you wanted to build a list of
results.
> > summary(result)
>
> snipped long and uninformative output.
David Winsemius, MD
West Hartford, CT
More information about the R-help
mailing list