How to extract factors names from anova function

#1

Dear R experts
I am quite confortable with R for Anova processes, but I found an unexpected difficulty to label "barplot2" for having a convenient sensitivity diagram of terms influences.
When typing :

y_TWI=aov(Donnees[,k] ~ (H1+H2+H3+H4+H5)*(H1+H2+H3+H4+H5))
anov_TWI=anova(y_TWI)

I get the right info:

Analysis of Variance Table
Response: Donnees[, k]
          Df  Sum Sq Mean Sq   F value    Pr(>F)    
H1         1    1.27    1.27    1.8089  0.181942    
H2         1   22.50   22.50   31.9387 1.768e-07 ***
H3         1 2229.56 2229.56 3164.7206 < 2.2e-16 ***
H4         1  626.99  626.99  889.9756 < 2.2e-16 ***
H5         1    7.50    7.50   10.6513  0.001545 ** 
H1:H2      1    0.01    0.01    0.0091  0.924039    
H1:H3      1    0.08    0.08    0.1145  0.735859    
H1:H4      1    0.06    0.06    0.0793  0.778848 

etc...
BUT, in order to automate what ever the model, I would like to extract the names of the terms "H1", "H2",..."H1.H2".... to feed the function argument "names.arg" in :

barplot2(anov_TWI[,4],xlab="Co-Variables",ylab="Signifiance ",main=Title,col="green",names.arg= ???

I did not succeed even after a lot of structure conversion attempts....
Thanks a lot for any help

0 Likes

#2

I would recommend broom.

a <- anova(aov(Sepal.Length ~ ., iris))

a <- broom::tidy(a)

a
#> # A tibble: 5 x 6
#>   term            df  sumsq  meansq statistic   p.value
#>   <chr>        <int>  <dbl>   <dbl>     <dbl>     <dbl>
#> 1 Sepal.Width      1  1.41   1.41       15.0   1.63e- 4
#> 2 Petal.Length     1 84.4   84.4       897.    1.01e-63
#> 3 Petal.Width      1  1.88   1.88       20.0   1.56e- 5
#> 4 Species          2  0.889  0.444       4.72  1.03e- 2
#> 5 Residuals      144 13.6    0.0941     NA    NA

a$term
#> [1] "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"     
#> [5] "Residuals"

Created on 2019-02-14 by the reprex package (v0.2.1)

0 Likes

#3

Hi Gilles, it's hard to figure out why you have problems to get your names of interest. A reprex would make it may more easy for us! Please, try to use 'reprex' (FAQ) in future. It's a cool tool to post your questions to the R-communtiy. It helps us to help you! :slight_smile:

As Kamil mentioned the broom package gives you a nice and easy way to summaries or handel different model objects.

Below you will find a solution similar to your code (without the use of the broom package, but I would try broom if you have to handle different model types)...

library(gplots)

#calculate anova example (see aov() help...)
npk.aov <- aov(yield ~ block + N*P*K, npk)
#get anova model summary
npk.summary <- anova(npk.aov)
npk.summary
#> Analysis of Variance Table
#> 
#> Response: yield
#>           Df Sum Sq Mean Sq F value   Pr(>F)   
#> block      5 343.29  68.659  4.4467 0.015939 * 
#> N          1 189.28 189.282 12.2587 0.004372 **
#> P          1   8.40   8.402  0.5441 0.474904   
#> K          1  95.20  95.202  6.1657 0.028795 * 
#> N:P        1  21.28  21.282  1.3783 0.263165   
#> N:K        1  33.13  33.135  2.1460 0.168648   
#> P:K        1   0.48   0.482  0.0312 0.862752   
#> Residuals 12 185.29  15.441                    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#drop residual information for barplot2 visualization
# as hint: the anova object can be handled like a data frame, because it is one! :)
fig.data <- subset(npk.summary, !rownames(npk.summary) %in% 'Residuals')
fig.data
#>       Df Sum Sq Mean Sq F value   Pr(>F)   
#> block  5 343.29  68.659  4.4467 0.015939 * 
#> N      1 189.28 189.282 12.2587 0.004372 **
#> P      1   8.40   8.402  0.5441 0.474904   
#> K      1  95.20  95.202  6.1657 0.028795 * 
#> N:P    1  21.28  21.282  1.3783 0.263165   
#> N:K    1  33.13  33.135  2.1460 0.168648   
#> P:K    1   0.48   0.482  0.0312 0.862752   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#show rownames or colnames
rownames(fig.data)
#> [1] "block" "N"     "P"     "K"     "N:P"   "N:K"   "P:K"
colnames(fig.data)
#> [1] "Df"      "Sum Sq"  "Mean Sq" "F value" "Pr(>F)"

#plot your stats of interest
Title <- 'Barplot'
barplot2(height = fig.data$'Pr(>F)', xlab = "Co-Variables", ylab = "p-value",
         main = Title, col = "green", names.arg = rownames(fig.data))
#add your level of significance
abline(h=0.05, lty=2)
#add legend
legend('topleft', 'Significance level',lty=2)

Created on 2019-02-14 by the reprex package (v0.2.1)

0 Likes

#4

Thanks Adam
I suceeded in using Broom package but I will take advantage from some tips you mentionned.
Concerning Reprex, I am not that much familiar with with it but I included it as the package to be installed in my script for further purpose..
Thanks again..

0 Likes

#5

Hi Slowkow
I had replied this by mail two days ago, but the message failed.. I copy it here !

Absolutely great !
Thanks a lot...
I would not have the idea to find such a function...
Have a nice day !

0 Likes

#6

If your question's been answered (even by you!), would you mind choosing a solution? It helps other people see which questions still need help, or find solutions if they have similar problems. Here’s how to do it:

0 Likes

#7

Again, thanks for the "broom" function allowing to extract the names of terms of the Anova model .
Just a remark concerning the limitation of this in case you use "SO" , "PQ" or "TWI" keyword to specify the type of model to be used. The naming is so heavy ( for instance ... SO(A,B,C,D)B:D to name the interaction term between B and D....) But this is another problem to solve on which I am testing some cases...
Have a nice WE

1 Like

closed #8

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.

0 Likes