Most of the experimental biologists measure a variable (for eg. cholesterol levels) under multiple conditions (for eg. measurements for three geographically different areas or under three different drug regimen). Now for the sake of simplicity, let us take in-build iris data, in R. Iris data has 4 different measurements (4 variables: Sepal length and width, Petal length and width) for 3 plants (Setosa, Versicolor and Viriginica). Now we want to execute anova between three species for each plant. For 4 variables, we would be running 4 different anova and corresponding code would use either loop or apply functions. While those approaches exist, we can do the same thing with purrr package in a different way. Let us write the code first and then explain (don't be distracted by size of code, much of it is for formatting and rearranging):
===========================================================
library(purrr)
library(dplyr)
library(tidyr)
library(broom)
library(knitr)
library(tibble)
library(kableExtra)
iris %>%
pivot_longer(-Species, names_to="measurement", values_to="value") %>%
nest(data=c(Species,value)) %>%
mutate(model=map(data, ~TukeyHSD(aov (value ~ Species, data=.))),
tidy=map(model,tidy)) %>%
unnest(c("tidy")) %>%
select(c(5,10)) %>%
add_column(Measurement="",.before="contrast") %>%
mutate(adj.p.value=-log10(adj.p.value)) %>%
kable(., "html", booktabs = T, escape = F, caption = "Flower measurement statistics",
col.names = c("","Comparison", "Ajdusted P value, (-log10p)")) %>%
kable_styling("striped",full_width = F) %>%
column_spec(2:3, width = "5cm") %>%
pack_rows("Sepal Length", 1, 3) %>%
pack_rows("Sepal Width", 4, 6) %>%
pack_rows("Petal Length", 7, 9) %>%
pack_rows("Petal Width", 10,12)
=============================================================
output:
Now this table has four measurements (Sepal Length, Width, Petal Length and width) for 3 different species (like 3 different conditions). Let us get to the code:============================================================
- iris %>% - Load the iris data
- pivot_longer(-Species, names_to="measurement", values_to="value") %>% - Converts iris data from wider format to long format. All measurement values are sent to column "value" and type of measurement is sent to column "measurement"
- nest(data=c(Species,value)) %>% -- Nest function creates list. But this list is named "data" and this list contains two columns: Species and value. For each measurement type, one list is created i.e for Sepal.Length one list is created. This list will have Plant type/Species and Sepal length measurement.
- mutate(model=map(data, ~TukeyHSD(aov (value ~ Species, data=.))),
tidy=map(model,tidy)) %>% -- Creates two columns: Model and tidy. Model column stores the output from TukeyHSD function for each list in data created in step 3. Tidy column will format output column "model" created in before step. - unnest(c("tidy")) %>% - Expand or flatten the tidy column
- select(c(5,10)) %>% - select the columns 5 and 10
- add_column(Measurement="",.before="contrast") %>% - Add a column (without values) for downstream beautification. Name the column "Measurement" and put it before column "contrast"
- mutate(adj.p.value=-log10(adj.p.value)) %>% -- Convert the p-values to -log10 p values
- col.names = c("","Comparison", "Ajdusted P value, (-log10p)")) %>%
kable_styling("striped",full_width = F) %>%column_spec(2:3, width = "5cm") %>%pack_rows("Sepal Length", 1, 3) %>%pack_rows("Sepal Width", 4, 6) %>%pack_rows("Petal Length", 7, 9) %>%pack_rows("Petal Width", 10,12)
Now step 9 generates the table for better viewing. Now how do we know that our code is correct. Let us try this on sepal width anova with species. Run following code:
sl_anova <- aov(Sepal.Width ~ Species, data=iris)
-log10(TukeyHSD(sl_anova)$Species[,4])
==============================================================
You should see 13.508986, 8.866226, 2.056495 and compare these values from table. They will be identical.