#get the mean of our root length
mean(growth_data$Rootlength)[1] NA7 Statistics in R
WarningWork in progress
This chapter was recently updated to new example datasets. The code runs, but the explanations of the results are drafts that are still being checked. Please read the interpretations with care.
A lot of helpful information that that was used here was inspired by a more detailed site that can be found here
And please take note, that the writer of this tutorial is no statistician and if you spot issues feel free to let me know!
7.1 Summary stats
To do some basic stats we use the growth data set.
mean()= print the mean of a certain columnsummary()= print the mean, median, max, min, … of a certain column.
The result is NA. Our root length column contains missing values, and the mean of a list of numbers that contains a missing value is unknown. With the argument na.rm = TRUE (“NA remove”), mean() ignores the missing values. Many other functions, such as sd(), median() or sum(), have the same argument.
#get the mean of our root length, ignoring missing values
mean(growth_data$Rootlength, na.rm = TRUE)[1] 75.19195#return summary stats for the root length
summary(growth_data$Rootlength) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
28.73 56.60 75.64 75.19 92.25 135.49 15 We can also summarize data of several columns
#return summary stats for the root length
summary(growth_data[,c("Rootlength", "FW_shoot_mg")]) Rootlength FW_shoot_mg
Min. : 28.73 Min. : 1.56
1st Qu.: 56.60 1st Qu.:12.10
Median : 75.64 Median :15.20
Mean : 75.19 Mean :16.58
3rd Qu.: 92.25 3rd Qu.:19.63
Max. :135.49 Max. :47.60
NA's :15 NA's :8 We can also run summary() on our complete dataframe.
#run summary on all our data
summary(growth_data) Nutrient Condition Experiment Rootlength
Length:594 Length:594 Length:594 Min. : 28.73
Class :character Class :character Class :character 1st Qu.: 56.60
Mode :character Mode :character Mode :character Median : 75.64
Mean : 75.19
3rd Qu.: 92.25
Max. :135.49
NA's :15
FW_shoot_mg
Min. : 1.56
1st Qu.:12.10
Median :15.20
Mean :16.58
3rd Qu.:19.63
Max. :47.60
NA's :8 For number columns, summary() also tells us how many values are missing (NA's). For text columns, it only tells us how many values there are.
7.2 The table command
The function table builds a contingency table of the counts at each factor level, or combinations of factor levels. This is quite useful to count the number of data points with certain metadata associated. So for example for our growth data we can ask how many measurements we made for each condition.
#summarize how many measurements we have for each condition
table(growth_data$Condition)
129E 13C mock
210 175 209 If we give table() two columns, we get the counts for each combination:
#summarize how many measurements we have for each combination of nutrient and condition
table(growth_data$Nutrient, growth_data$Condition)
129E 13C mock
high_P 103 87 104
low_P 107 88 105Here, we see that we have a slightly different number of measurements for each combination of nutrient and condition. For this tutorial we can ignore them, but this might be relevant if we have huge differences for some statistical analyses. This approach allows you to for larger datasets to easily check for samples that might be outliers in terms of measurements per sample.
7.3 The ddply command
The table command can get rather slow and there are some useful packages to speed up things and run more complicated mathematical operations. One example is the ddply() function of the plyr package. A useful feature of ddply is that this tool stores the output as a dataframe instead of a table.
Below we see examples were we summarize the data for the root length across different nutrient conditions. I.e. we do not want to have the mean for all root lengths but we want to see if roots have different lengths under low and high P conditions.
#load our package
library(plyr)
#calculate the mean root length across our Sample IDs
growth_data_summarized <- ddply(growth_data, .(Nutrient), summarize, Mean_RootLength = mean(Rootlength, na.rm = TRUE))
#view data
head(growth_data_summarized)The structure of the command is as follows:
ddply(Input_Data, .(Colum_we_want_to_sum_arcoss), summarize, New_Column_name = mean(Column_we_want_to_calc_the_mean_for))
We can also summarize the data combining different conditions (i.e. nutrient and condition).
#Summarize our data across Nutrients and Conditions
growth_data_summarized <- ddply(growth_data, .(Nutrient, Condition), summarize, Mean_RootLength = mean(Rootlength, na.rm = TRUE))
#view data
head(growth_data_summarized)We can also calculate the mean, sd and se in one line of code
#and we can use even fancier functions (i.e. get the se and sd), check the plyr package for more details
#for the se we need the number of measurements; sum(!is.na(Rootlength)) counts only the values that are not missing
growth_data_summarized <- ddply(growth_data, .(Nutrient, Condition), summarize, RootLength = mean(Rootlength, na.rm = TRUE), sd = sd(Rootlength, na.rm = TRUE), se = sd(Rootlength, na.rm = TRUE) / sqrt(sum(!is.na(Rootlength))))
#view data
head(growth_data_summarized)7.4 tidyr
Tydir is another package to summarize data but also to transform data. For now we will just discuss the basics to summarize data but we will try to add an extended section on this package later.
You notice here that the syntax is a bit unusual and we use the %>% symbol (a so-called forward pipe operator). This symbol is commonly used in the dplyr and tidyr packages. This function passes the left hand side of the operator to the first argument of the right hand side of the operator.
In the following example, the a subset of the growth_data data (only 3 columns, not the whole dataframe) gets passed to count()
New functions:
count()= A function of the dplyr package. Here, we count the number of rows for each condition (or each combination of the columns we give it)summarize()creates new data frame into one (or more) rows for each combination of grouping variables; if there are no grouping variables, the output will have a single row summarising all observations in the input. It will contain one column for each grouping variable and one column for each of the summary statistics that you have specified.mean()= calculate the arithmetic mean.
#read in package
library(tidyverse)#count how many measurements we have per conditions
growth_data %>% count(Condition, sort = FALSE)#count how many measurements we have per nutrient and conditions
growth_data %>% count(Nutrient, Condition, sort = FALSE)#get more comprehensive data stats and summarize for the whole dataset
growth_data %>% summarise(
count = n(),
mean_root = mean(Rootlength, na.rm = TRUE))Good things about tydr are:
- it is extremely fast, so is important for dealing with larger datasets
- we can combine several commands by using the pipe
I.e. below we can see that we first group the input data and then summarized the group data.
#get more comprehensive data stats and summarize and grouping for the different conditions
growth_data %>%
group_by(Nutrient, Condition) %>%
summarise(
count = n(),
mean_root = mean(Rootlength, na.rm = TRUE) ,.groups = 'drop')We can also do other stats than just calculating the mean:
#get more comprehensive data stats and summarize and grouping for the different conditions
growth_data %>%
group_by(Nutrient, Condition) %>%
summarise(
count = n(),
mean_root = mean(Rootlength, na.rm = TRUE),
sd_root = sd(Rootlength, na.rm = TRUE),.groups = 'drop')Other useful summary functions:
mean(x): sum of x divided by the lengthmedian(x): 50% of x is above and 50% is belowsd(x): standard deviationIQR(x): interquartile range (robust equivalent of sd when outliers are present in the data)mad(x): median absolute deviation (robust equivalent of sd when outliers are present in the data)min(x): minimum value of xmax(x): maximum value of xquantile(x, 0.25): 25% of x is below this valuefirst(x): equivalent to x[1]nth(x, 2): equivalent to n<-2; x[n]last(x): equivalent to x[length(x)]n(x): the number of elements in xsum(!is.na(x)): count non-missing valuesn_distinct(x): count the number of unique valuesum(x > 10): count the number of elements where x > 10mean(y == 0): proportion of elements where y = 0
7.5 More stats: tapply
Tapply is a base R function and allows to apply a function to each cell of a ragged array, that is to each (non-empty) group of values given by a unique combination of the levels of certain factors. It can be used as an alternative to ddply.
Suppose now we want to estimate the mean root length for each growth condition. Notice, how we can vary the index?
#index = factors we want to select
tapply(X = growth_data$Rootlength, INDEX = growth_data$Nutrient, FUN = mean,na.rm = TRUE) high_P low_P
90.48900 59.94766 #same but for two-way tables (this is not so useful here, but might be handy when you have different conditions for the same organism or a timecourse)
tapply(growth_data$Rootlength, INDEX = list(growth_data$Nutrient, growth_data$Condition),FUN = mean, na.rm = TRUE) 129E 13C mock
high_P 98.35854 97.68835 76.40475
low_P 69.91808 64.19465 47.068487.6 Linear models
7.6.1 Basics
It is very simple to investigate linear relationships among variables in R. We want to estimate how a quantitative dependent variable changes according to the levels of one or more categorical independent variables. In the command below, the linear relationship between Rootlength (the dependent variable, i.e. the one we’re trying to predict) and FW_shoot_mg (the independent variable or the predictor) is calculated.
We need to use the function summary to see the results of that command; coef extracts the best-fit coefficients, anova performs an analysis of variance; there are many other extractor functions.
In this case lets work with the data were we do not have to deal with the different time points to simplify things.
#read in data were we have several measurements
growth_data <- read.table("data/Growth_Data.txt", sep = "\t", header = TRUE, quote = "")
#is there a correlation between freshweight and root length?
linearmodel <- lm(Rootlength ~ FW_shoot_mg, data = growth_data)
linearmodel
Call:
lm(formula = Rootlength ~ FW_shoot_mg, data = growth_data)
Coefficients:
(Intercept) FW_shoot_mg
45.471 1.782 #let's extract the entire t-table
summary(linearmodel)
Call:
lm(formula = Rootlength ~ FW_shoot_mg, data = growth_data)
Residuals:
Min 1Q Median 3Q Max
-59.467 -16.619 1.126 13.202 50.139
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 45.4713 2.2970 19.80 <2e-16 ***
FW_shoot_mg 1.7818 0.1279 13.94 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.64 on 569 degrees of freedom
(23 observations deleted due to missingness)
Multiple R-squared: 0.2544, Adjusted R-squared: 0.2531
F-statistic: 194.2 on 1 and 569 DF, p-value: < 2.2e-16Here one of the values is the model p-Value (bottom last line) and the p-Value of individual predictor variables (extreme right column under ‘Coefficients’). The p-Values are very important because, we can consider a linear model to be statistically significant only when both these p-Values are less that the pre-determined statistical significance level, which is ideally 0.05. This is visually interpreted by the significance stars at the end of the row. The more the stars beside the variable’s p-Value, the more significant the variable.
- Residuals: The section summarizes the residuals, the error between the prediction of the model and the actual results. Smaller residuals are better.
- Coefficients: For each variable and the intercept, a weight is produced and that weight has other attributes like the standard error, a t-test value and significance.
- Estimate: This is the weight given to the variable. In the simple regression case (one variable plus the intercept), it tells you by how much the root length changes, according to the model, for every increase of the shoot fresh weight by 1 mg. In our data, the estimate is 1.78: for every additional mg of shoot fresh weight, the model predicts a root that is 1.78 mm longer.
lm() automatically leaves out all plants with a missing root length or a missing shoot weight. You can see how many it removed in the line “observations deleted due to missingness” of the summary. * Std. Error: Tells you how precisely was the estimate measured. It’s really only useful for calculating the t-value. * t-value and Pr(>[t]): The t-value is calculated by taking the coefficient divided by the Std. Error. It is then used to test whether or not the coefficient is significantly different from zero. If it isn’t significant, then the coefficient really isn’t adding anything to the model and could be dropped or investigated further. Pr(>|t|) is the significance level.
Performance Measures:
Three sets of measurements are provided.
- Residual Standard Error: This is the standard deviation of the residuals. Smaller is better.
- Multiple / Adjusted R-Square: For one variable, the distinction doesn’t really matter. R-squared shows the amount of variance explained by the model. Adjusted R-Square takes into account the number of variables and is most useful for multiple-regression.
- F-Statistic: The F-test checks if at least one variable’s weight is significantly different than zero. This is a global test to help asses a model. If the p-value is not significant (e.g. greater than 0.05) than your model is essentially not doing anything.
We also can just print parts of the data:
#print only the coefficients
coef(summary(linearmodel)) Estimate Std. Error t value Pr(>|t|)
(Intercept) 45.471292 2.2969633 19.79626 9.470283e-67
FW_shoot_mg 1.781836 0.1278697 13.93478 3.480146e-38#print only the anova stats
anova(linearmodel)#plot add the best line to a plot
with(growth_data, plot(Rootlength ~ FW_shoot_mg, col = 2))
abline(linearmodel)
If we look at the stats and the p value we see a decent correlation. However, since we tested different nutrient conditions it also can be useful to separate data to get a clearer picture. Something else to consider is to ask whether the data is normally distributed and based on that what statistical test to choose.
#color the points by nutrient condition to check this
with(growth_data, plot(Rootlength ~ FW_shoot_mg, col = factor(Nutrient)))
legend("topleft", legend = levels(factor(growth_data$Nutrient)), col = 1:2, pch = 1)
7.6.2 Analysing residuals
Anyone can fit a linear model in R. The real test is analyzing the residuals (the error or the difference between actual and predicted results).
There are four things we are looking for when analyzing residuals.
- The mean of the errors is zero (and the sum of the errors is zero)
- The distribution of the errors are normal.
- All of the errors are independent.
- Variance of errors is constant (Homoscedastic)
In R, you pull out the residuals with resid() and the values the model predicts (the fitted values) with fitted(). Using the simple linear regression model (simple.fit) we’ll plot a few graphs to help illustrate any problems with the model.
Below some examples:
simple.fit <- linearmodel
layout(matrix(c(1,1,2,3),2,2,byrow=TRUE))
#Fitted values x Residuals Plot
plot(resid(simple.fit) ~ fitted(simple.fit),
main="Fitted values x Residuals\nfor Simple Regression",
xlab="Fitted root length", ylab="Residuals")
abline(h=0,lty=2)
#Histogram of Residuals
hist(simple.fit$resid, main="Histogram of Residuals",
ylab="Residuals")
#Q-Q Plot
qqnorm(simple.fit$resid)
qqline(simple.fit$resid)
The histogram and QQ-plot are the ways to visually evaluate if the residual fit a normal distribution.
- If the histogram looks like a bell-curve it might be normally distributed.
- If the QQ-plot has the vast majority of points on or very near the line, the residuals may be normally distributed.
7.7 A new dataset for statistical tests
For the statistical tests in the rest of this chapter, we use a smaller and cleaner dataset from the same study (Garrido-Oter et al. 2018). Plants of three Arabidopsis genotypes were grown with or without the strain 129E, and the authors counted the number of cells in the meristem of the root tip, the zone where new root cells are made. The genotypes are:
Col-0: the wild type, i.e. the standard Arabidopsis line most labs useahk3-3andarr1_arr12: two mutant lines
We use this dataset because the groups are small, similar in size and well-behaved, which makes it easier to see how the tests work. A second file with the root length of the same genotypes is used in the exercise at the end of this chapter.
Download the two files in the same way as in the first steps chapter:
# download the meristem data, but only if we don't have it yet
if (!file.exists("data/Genotype_Meristem.txt")) {
download.file(url = "https://raw.githubusercontent.com/ndombrowski/R_cookbook/main/data/Genotype_Meristem.txt",
destfile = "data/Genotype_Meristem.txt",
mode = "wb")
}
# download the root length data of the same genotypes, but only if we don't have it yet
if (!file.exists("data/Genotype_Rootlength.txt")) {
download.file(url = "https://raw.githubusercontent.com/ndombrowski/R_cookbook/main/data/Genotype_Rootlength.txt",
destfile = "data/Genotype_Rootlength.txt",
mode = "wb")
}# read in the meristem data
meristem <- read.table("data/Genotype_Meristem.txt", sep = "\t", header = TRUE, quote = "")
# check the data
head(meristem)dim(meristem)[1] 80 3# how many plants do we have per genotype and condition?
table(meristem$Genotype, meristem$Condition)
129E mock
ahk3-3 14 15
arr1_arr12 13 15
Col-0 11 127.8 Comparing two groups: the t-test
Let’s start with the simplest question: does 129E change the number of meristem cells in the wild type? For this, we only need the Col-0 plants:
# keep only the wild type
col0 <- meristem[meristem$Genotype == "Col-0", ]
# look at the data
boxplot(Meristem_cells ~ Condition, data = col0, ylab = "Number of meristem cells")
A t-test compares the means of two groups. It asks: if there was really no difference between the groups, how likely is a difference between the means as large as the one we see? This likelihood is the p-value. If it is smaller than 0.05, we usually call the difference statistically significant.
The t-test assumes that the values in each group come from a normal distribution. We will look at this in more detail below, but for now we check it with the Shapiro-Wilk test, once per group:
# test for normality in each group
shapiro.test(col0$Meristem_cells[col0$Condition == "mock"])
Shapiro-Wilk normality test
data: col0$Meristem_cells[col0$Condition == "mock"]
W = 0.96225, p-value = 0.8154shapiro.test(col0$Meristem_cells[col0$Condition == "129E"])
Shapiro-Wilk normality test
data: col0$Meristem_cells[col0$Condition == "129E"]
W = 0.92528, p-value = 0.3651The Shapiro-Wilk test works the other way round: here, a p-value above 0.05 is what we want, because it means that there is no evidence against a normal distribution.
Now the t-test itself. The ~ reads as “depends on”: we test whether the number of meristem cells depends on the condition.
# compare mock and 129E in the wild type
t.test(Meristem_cells ~ Condition, data = col0)
Welch Two Sample t-test
data: Meristem_cells by Condition
t = 3.4714, df = 19.311, p-value = 0.002508
alternative hypothesis: true difference in means between group 129E and group mock is not equal to 0
95 percent confidence interval:
2.988992 12.041311
sample estimates:
mean in group 129E mean in group mock
38.18182 30.66667 The most important lines of the output are:
tandp-value: the test statistic and the p-value.95 percent confidence interval: the range in which the true difference between the two means probably lies.mean in group ...: the mean of each group.
By default, R runs the Welch version of the t-test, which does not assume that both groups have the same variance. That is why the title says “Welch Two Sample t-test”.
In our data, wild-type plants treated with 129E have on average 38.2 meristem cells, and mock plants 30.7. The difference is significant (p = 0.0025), and the confidence interval tells us that 129E adds roughly 3 to 12 cells. Both Shapiro-Wilk tests above gave p-values well above 0.05 (0.82 and 0.37), so we can use a t-test here.
If your data is clearly not normally distributed, you can use the Wilcoxon rank sum test (also called Mann-Whitney U test) instead. It does not compare the means, but uses the ranks of the values, so it does not need a normal distribution:
# a test that does not assume a normal distribution
wilcox.test(Meristem_cells ~ Condition, data = col0)Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
compute exact p-value with ties
Wilcoxon rank sum test with continuity correction
data: Meristem_cells by Condition
W = 114, p-value = 0.003382
alternative hypothesis: true location shift is not equal to 0You may see a warning that R “cannot compute exact p-value with ties”. Ties are values that occur more than once, which is common for counts such as our cell numbers. R then calculates an approximate p-value instead, which is fine here. The Wilcoxon test comes to the same conclusion as the t-test (p = 0.003).
7.9 Checking for a normal distribution
Many statistical tests, such as the t-test and the ANOVA below, assume that the values in each group come from a normal distribution. It is important to check this per group: if you put all plants together, you mix groups with different means, and the combined data will not look normal even if every group is.
7.9.1 Density plots
A density plot is a smoothed version of a histogram. The ggdensity() function of the ggpubr package draws one curve per group. With facet.by we get one panel per genotype:
library(ggpubr)
ggdensity(meristem, x = "Meristem_cells",
add = "mean", rug = TRUE,
color = "Condition", palette = c("#2A6FDB", "#8C8C8C"),
facet.by = "Genotype")
In all three genotypes, the curve of the 129E plants lies to the right of the mock curve, so treated plants tend to have more meristem cells. The shift is largest in arr1_arr12. Some curves have small bumps, for example Col-0 with 129E, but with only 11 to 15 plants per group, a density curve is never perfectly smooth.
7.9.2 Q-Q plots
A Q-Q plot (quantile-quantile plot) compares the values in a group with the values you would expect from a normal distribution. If the data is normally distributed, the points lie close to the diagonal line. The grey band shows the range in which we would expect the points.
ggqqplot(meristem, x = "Meristem_cells",
color = "Condition", palette = c("#2A6FDB", "#8C8C8C"),
facet.by = "Genotype")
In all six groups, the points lie close to the line and inside the grey band, so there is no sign of a strong deviation from a normal distribution.
7.9.3 Shapiro-Wilk test for each group
Instead of running shapiro.test() six times, we can use the normality() function of the dlookr package together with group_by(). It runs one test for each combination of genotype and condition:
library(dplyr)
library(dlookr)
meristem %>%
group_by(Genotype, Condition) %>%
normality()The Shapiro-Wilk test checks the null hypothesis that the values come from a normal distribution. If the p-value is smaller than 0.05, we reject this hypothesis and should not assume a normal distribution.
In our data, all six p-values are between 0.37 and 0.87. None of the groups deviates significantly from a normal distribution, which fits what we saw in the plots.
Keep in mind that the result of the test depends on the number of values. With small groups like ours (11 to 15 plants), the test can easily miss a deviation from the normal distribution. With very large groups, it flags even tiny deviations that don’t matter in practice. So always look at the plots as well, not only at the p-value.
7.10 ANOVA
With a t-test we can compare two groups. Our data has six groups: three genotypes, each with and without 129E. To compare more than two groups, we use an analysis of variance (ANOVA).
Let’s first look at the data:
ggboxplot(meristem, x = "Genotype", y = "Meristem_cells", color = "Condition",
palette = c("#2A6FDB", "#8C8C8C"))
In every genotype, plants treated with 129E have more meristem cells than mock plants. The arr1_arr12 plants with 129E stand out, with clearly more cells than all other groups.
We start with an ANOVA that tests whether the number of meristem cells depends on the genotype and on the condition:
# two-way ANOVA, additive model
aov_add <- aov(Meristem_cells ~ Genotype + Condition, data = meristem)
summary(aov_add) Df Sum Sq Mean Sq F value Pr(>F)
Genotype 2 2958 1479.1 38.21 3.28e-12 ***
Condition 1 2646 2646.0 68.36 3.38e-12 ***
Residuals 76 2942 38.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The output has one row per factor. Pr(>F) is the p-value, and the stars at the end of the row show how small it is (see the Signif. codes line).
Both the genotype and the condition have a significant effect on the number of meristem cells (both p < 0.001).
The model above is called an additive model. It assumes that the effect of 129E is the same in every genotype. But maybe 129E has a larger effect in one genotype than in another. This is called an interaction between the two factors. To include it, we replace the + by a *:
# Two-way ANOVA with interaction effect
# These two calls are equivalent
aov_int <- aov(Meristem_cells ~ Genotype * Condition, data = meristem)
aov_int <- aov(Meristem_cells ~ Genotype + Condition + Genotype:Condition, data = meristem)
# summarize the ANOVA
summary(aov_int) Df Sum Sq Mean Sq F value Pr(>F)
Genotype 2 2958.1 1479.1 49.83 1.96e-14 ***
Condition 1 2646.0 2646.0 89.14 2.45e-14 ***
Genotype:Condition 2 745.3 372.7 12.55 2.02e-05 ***
Residuals 74 2196.6 29.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The new row Genotype:Condition tests the interaction. It is significant (p = 2e-05): the effect of 129E on the number of meristem cells depends on the genotype.
If the interaction is not significant, use the additive model. If it is significant, as here, keep the interaction in the model, and be careful with statements about the main effects alone: “129E increases the number of meristem cells” is then only part of the story, because the size of the effect depends on the genotype.
An interaction plot shows this well. Each line is one condition, and the points are the group means. If there is no interaction, the lines are roughly parallel:
interaction.plot(x.factor = meristem$Genotype, trace.factor = meristem$Condition,
response = meristem$Meristem_cells,
xlab = "Genotype", ylab = "Mean number of meristem cells", trace.label = "Condition")
Here the lines are clearly not parallel. In arr1_arr12, 129E adds about 20 cells on average, in ahk3-3 and Col-0 only about 7.
7.11 Which groups differ? Tukey’s test
A significant ANOVA tells us that some group means differ, but not which ones. To find out, we compare all pairs of groups. Because we make many comparisons at the same time, the chance that at least one of them is significant just by chance becomes large. Tukey’s HSD (Honest Significant Differences) test corrects for this. It is a post-hoc test, which means that it is used after an ANOVA.
# compare all pairs of groups
TukeyHSD(aov_int) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Meristem_cells ~ Genotype + Condition + Genotype:Condition, data = meristem)
$Genotype
diff lwr upr p adj
arr1_arr12-ahk3-3 12.060345 8.607844 15.512845 0.0000000
Col-0-ahk3-3 -1.428786 -5.067208 2.209637 0.6173667
Col-0-arr1_arr12 -13.489130 -17.156178 -9.822083 0.0000000
$Condition
diff lwr upr p adj
mock-129E -11.51505 -13.94553 -9.08457 0
$`Genotype:Condition`
diff lwr upr p adj
arr1_arr12:129E-ahk3-3:129E 19.2417582 13.1022369 25.3812796 0.0000000
Col-0:129E-ahk3-3:129E -0.9610390 -7.3834481 5.4613701 0.9978852
ahk3-3:mock-ahk3-3:129E -6.6761905 -12.5996807 -0.7527003 0.0180753
arr1_arr12:mock-ahk3-3:129E -0.6095238 -6.5330140 5.3139664 0.9996550
Col-0:mock-ahk3-3:129E -8.4761905 -14.7469530 -2.2054279 0.0023220
Col-0:129E-arr1_arr12:129E -20.2027972 -26.7329888 -13.6726056 0.0000000
ahk3-3:mock-arr1_arr12:129E -25.9179487 -31.9581308 -19.8777666 0.0000000
arr1_arr12:mock-arr1_arr12:129E -19.8512821 -25.8914642 -13.8110999 0.0000000
Col-0:mock-arr1_arr12:129E -27.7179487 -34.0990557 -21.3368417 0.0000000
ahk3-3:mock-Col-0:129E -5.7151515 -12.0426642 0.6123612 0.0999201
arr1_arr12:mock-Col-0:129E 0.3515152 -5.9759975 6.6790278 0.9999836
Col-0:mock-Col-0:129E -7.5151515 -14.1688831 -0.8614200 0.0177130
arr1_arr12:mock-ahk3-3:mock 6.0666667 0.2462015 11.8871318 0.0360994
Col-0:mock-ahk3-3:mock -1.8000000 -7.9735355 4.3735355 0.9562146
Col-0:mock-arr1_arr12:mock -7.8666667 -14.0402022 -1.6931311 0.0048678The output has one block per factor. The last block, $Genotype:Condition, compares all six groups with each other. For each pair, diff is the difference between the means, lwr and upr give the 95% confidence interval of this difference, and p adj is the p-value, adjusted for the number of comparisons.
The comparisons we are most interested in are mock against 129E within the same genotype. All three are significant: ahk3-3 (p adj = 0.018), arr1_arr12 (p adj < 0.001) and Col-0 (p adj = 0.018). Without 129E, arr1_arr12 has slightly more cells than the other two genotypes (p adj = 0.036 against ahk3-3 and 0.005 against Col-0), while ahk3-3 and Col-0 do not differ (p adj = 0.96).
With six groups, this table is long. In figures, the result is often shown with letters instead: groups that share a letter are not significantly different. The HSD.test() function of the agricolae package does this. It needs one factor, so we first make a column that combines genotype and condition:
library(agricolae)
# one column with the six groups
meristem$Group <- paste(meristem$Genotype, meristem$Condition, sep = "_")
# ANOVA with the six groups as one factor
aov_group <- aov(Meristem_cells ~ Group, data = meristem)
# Tukey test with letters
HSD.test(aov_group, "Group", group = TRUE, console = TRUE)
Study: aov_group ~ "Group"
HSD Test for Meristem_cells
Mean Square Error: 29.68326
Group, means
Meristem_cells std r se Min Max Q25 Q50 Q75
ahk3-3_129E 39.14286 5.375974 14 1.456102 30 49 36.50 40.0 41.75
ahk3-3_mock 32.46667 3.997618 15 1.406728 25 39 30.00 32.0 35.00
arr1_arr12_129E 58.38462 7.228860 13 1.511068 48 71 53.00 58.0 62.00
arr1_arr12_mock 38.53333 5.449334 15 1.406728 29 48 35.00 40.0 42.50
Col-0_129E 38.18182 5.671299 11 1.642704 31 47 33.50 37.0 43.00
Col-0_mock 30.66667 4.599078 12 1.572770 22 40 28.75 29.5 33.00
Alpha: 0.05 ; DF Error: 74
Critical Value of Studentized Range: 4.137591
Groups according to probability of means differences and alpha level( 0.05 )
Treatments with the same letter are not significantly different.
Meristem_cells groups
arr1_arr12_129E 58.38462 a
ahk3-3_129E 39.14286 b
arr1_arr12_mock 38.53333 b
Col-0_129E 38.18182 bc
ahk3-3_mock 32.46667 cd
Col-0_mock 30.66667 darr1_arr12 with 129E gets its own letter (a), so it differs from all other groups. Within each genotype, the mock and the 129E group share no letter, which matches the Tukey table above: 129E changes the number of meristem cells in all three genotypes.
7.12 Checking the assumptions of the ANOVA
Besides a normal distribution, the ANOVA assumes that all groups have a similar variance (homogeneity of variances). We can check both assumptions on the residuals of the model, i.e. the differences between each value and the mean of its group.
7.12.1 Homogeneity of variances
The residuals versus fitted values plot shows the residuals of each group (the fitted values are the group means). The spread of the points should be similar in all groups:
plot(aov_int, 1)
R labels the three most extreme values with their row numbers. In our data, the spread is similar in all six groups, and the red line stays close to 0. Plants 39, 48 and 51 are marked as the most extreme values, but they are not far away from the rest.
We can also test this with Levene’s test from the car package. Its null hypothesis is that all groups have the same variance, so here again a p-value above 0.05 is what we want:
library(car)
leveneTest(Meristem_cells ~ Genotype * Condition, data = meristem)With p = 0.30 there is no evidence that the variances differ between the groups.
7.12.2 Normality of the residuals
The Q-Q plot of the residuals should roughly follow the straight line:
plot(aov_int, 2)
The points follow the line closely, with a small deviation at the upper end (plants 39 and 48). Together with Levene’s test, this means that the assumptions of the ANOVA are met for this dataset.
7.13 Exercise
The file Genotype_Rootlength.txt contains the root length (in mm) of plants from the same genotypes, with and without 129E.
- Read the file and check how many plants you have per genotype and condition.
- Does 129E change the root length of the wild type (
Col-0)? Use a t-test. - Run a two-way ANOVA with an interaction. Is the interaction significant?
- Check the assumptions with Levene’s test. What do you notice?
TipSolution
# 1. read the data and count the plants
genotype_root <- read.table("data/Genotype_Rootlength.txt", sep = "\t", header = TRUE, quote = "")
table(genotype_root$Genotype, genotype_root$Condition)
129E mock
ahk3-3 22 26
arr1_arr12 27 25
Col-0 18 25# 2. t-test for the wild type
col0_root <- genotype_root[genotype_root$Genotype == "Col-0", ]
t.test(Rootlength ~ Condition, data = col0_root)
Welch Two Sample t-test
data: Rootlength by Condition
t = 2.2912, df = 34.99, p-value = 0.02809
alternative hypothesis: true difference in means between group 129E and group mock is not equal to 0
95 percent confidence interval:
0.5924654 9.8066013
sample estimates:
mean in group 129E mean in group mock
42.85833 37.65880 # 3. two-way ANOVA with interaction
aov_root <- aov(Rootlength ~ Genotype * Condition, data = genotype_root)
summary(aov_root) Df Sum Sq Mean Sq F value Pr(>F)
Genotype 2 12129 6064 68.67 < 2e-16 ***
Condition 1 4865 4865 55.08 1.10e-11 ***
Genotype:Condition 2 1969 984 11.15 3.27e-05 ***
Residuals 137 12099 88
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 4. Levene's test
leveneTest(Rootlength ~ Genotype * Condition, data = genotype_root)As for the meristem cells, 129E has a significant effect on the root length of the wild type (p = 0.028), and the interaction between genotype and condition is significant (p = 3e-05). But this time, Levene’s test is significant too (p = 0.008): the variances are not the same in all groups, so one assumption of the ANOVA is not met. With groups of similar size, the ANOVA is usually still fairly reliable, but you should mention this when you report the result. An alternative is a test that does not assume equal variances, for example Welch’s ANOVA with oneway.test().
Garrido-Oter, Ruben, Ryohei Thomas Nakano, Nina Dombrowski, Ka-Wai Ma, Alice C. McHardy, and Paul Schulze-Lefert. 2018. “Modular Traits of the Rhizobiales Root Microbiota and Their Evolutionary Relationship with Symbiotic Rhizobia.” Cell Host & Microbe 24 (1): 155–167.e5. https://doi.org/10.1016/j.chom.2018.06.006.