I was recently running a series of ANOVA analyses, and I used the aov() function because it had a few options that I preferred. Much like lm(), the function returns an object that you typically pass to summary() to view and interpret the output.

It took me a bit of playing to figure out how to extract the information I needed. My aov object is called "fullmod", and here is the summary() output:

> summary(fullmod)

Error: sample

Df Sum Sq Mean Sq

Type 1 1.4535 1.4535

Error: sample:treatment

Df Sum Sq Mean Sq

treatment 1 0.0086021 0.0086021

Error: Within

Df Sum Sq Mean Sq F value Pr(>F)

treatment 1 0.0077 0.00769 0.0325 0.857505

Type 3 0.3159 0.10529 0.4443 0.722043

age 1 0.0351 0.03512 0.1482 0.701367

sex 1 2.2166 2.21661 9.3547 0.003122 **

hybridization 1 0.1125 0.11255 0.4750 0.492923

Residuals 72 17.0605 0.23695

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I assigned the summary output to an object s so I could examine it a bit further.

> s <- summary(fullmod)

Then I peeked under the hood using the str() function.

> str(s)

List of 3

$ Error: sample :List of 1

..$ :Classes ‘anova’ and 'data.frame': 1 obs. of 3 variables:

.. ..$ Df : num 1

.. ..$ Sum Sq : num 1.45

.. ..$ Mean Sq: num 1.45

..- attr(*, "class")= chr [1:2] "summary.aov" "listof"

$ Error: sample:treatment:List of 1

..$ :Classes ‘anova’ and 'data.frame': 1 obs. of 3 variables:

.. ..$ Df : num 1

.. ..$ Sum Sq : num 0.0086

.. ..$ Mean Sq: num 0.0086

..- attr(*, "class")= chr [1:2] "summary.aov" "listof"

$ Error: Within :List of 1

..$ :Classes ‘anova’ and 'data.frame': 6 obs. of 5 variables:

.. ..$ Df : num [1:6] 1 3 1 1 1 72

.. ..$ Sum Sq : num [1:6] 0.0077 0.3159 0.0351 2.2166 0.1125 ...

.. ..$ Mean Sq: num [1:6] 0.0077 0.1053 0.0351 2.2166 0.1125 ...

.. ..$ F value: num [1:6] 0.0325 0.4443 0.1482 9.3547 0.475 ...

.. ..$ Pr(>F) : num [1:6] 0.8575 0.72204 0.70137 0.00312 0.49292 ...

..- attr(*, "class")= chr [1:2] "summary.aov" "listof"

- attr(*, "class")= chr "summary.aovlist"

This tells me that s is a list of 3 objects ("Error: sample", "Error: sample:treatment:", and "Error: within"). Each of these objects is a list of 1 element for (some reason!?) with no name or identifier. For "Error: sample" and "Error: sample:treatment:", there are three single numbers. For "Error: Within", each element contains a list of 6 numbers (one for each coefficient in my model).

So, to access the degrees of freedom ("Df") for the "Error: sample", I do the following:

> s[[1]][[1]][[1]]

[1] 1

Or if I prefer to use labels (for the elements that have them)

> s[["Error: sample"]][[1]][["Df"]]

[1] 1

Just to clarify:

s [["Error: sample"]] [[1]] [["Df"]]

^ ^ ^

1st layer 2nd unnamed layer 3rd layer

It is VERY important to use [[ instead of [ here. In the R syntax, [[ retrieves a SINGLE element of the object/array and [ retrieves the full object as a list.

> s["Error: sample"][1]["Df"]

$

NULL

See: that didn't work.

To retrieve the p-values I'm interested in, I'll have to dig one more layer deeper to pull the p-value for the coefficient I want (in this case, the 4th element).

> s[[3]][[1]][[5]][[4]]

[1] 0.003122076

OR

> s[["Error: Within"]][[1]][["Pr(>F)"]][[4]]

[1] 0.003122076

Using this basic strategy, you should be able to pull any value from a summary object, but remember that using offsets as I have can be dangerous. For this to work consistently across many models, they must all have the same number of coefficients in the model.

Fantastic! Thank you very much!

ReplyDeleteHello

ReplyDeleteI was able to extract the proportion of variance from a PCA, so thank you very much.

I would like to use these instead of "PC1" and "PC2" on a plot, but I couldn't figure out how to do it. Could you help me?

Great! Thank U from Schiedam!

ReplyDeleteExactly what I needed, and clear instructions. Thanks!

ReplyDelete