- What is two-way ANOVA test?
- Two-way ANOVA test hypotheses
- Assumptions of two-way ANOVA test
- Compute two-way ANOVA test in R: balanced designs
- Compute two-way ANOVA test in R for unbalanced designs
- See also
- Read more
- Infos

The grouping variables are also known as factors. The different categories (groups) of a factor are called levels. The number of levels can vary between factors. The level combinations of factors are called cell.

When the sample sizes within cells are equal, we have the so-called

**balanced design**. In this case the standard two-way ANOVA test can be applied.- When the sample sizes within each level of the independent variables are not the same (case of
**unbalanced designs**), the ANOVA test should be handled differently.

This tutorial describes how to compute **two-way ANOVA test** in R software for balanced and unbalanced designs.

- There is no difference in the means of factor A
- There is no difference in means of factor B
- There is no interaction between factors A and B

The alternative hypothesis for cases 1 and 2 is: the means are not equal.

The alternative hypothesis for case 3 is: there is an interaction between A and B.

Two-way ANOVA, like all ANOVA tests, assumes that the observations within each cell are normally distributed and have equal variances. We’ll show you how to check these assumptions after fitting ANOVA.

Balanced designs correspond to the situation where we have equal sample sizes within levels of our independent grouping levels.

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use the built-in R data set named *ToothGrowth*. It contains data from a study evaluating the effect of vitamin C on tooth growth in Guinea pigs. The experiment has been performed on 60 pigs, where each animal received one of three dose levels of vitamin C (0.5, 1, and 2 mg/day) by one of two delivery methods, (orange juice or ascorbic acid (a form of vitamin C and coded as VC). Tooth length was measured and a sample of the data is shown below.

```
# Store the data in the variable my_data
my_data <- ToothGrowth
```

To get an idea of what the data look like, we display a random sample of the data using the function **sample_n**()[in **dplyr** package]. First, install dplyr if you don’t have it:

`install.packages("dplyr")`

```
# Show a random sample
set.seed(1234)
dplyr::sample_n(my_data, 10)
```

```
len supp dose
38 9.4 OJ 0.5
36 10.0 OJ 0.5
37 8.2 OJ 0.5
50 27.3 OJ 1.0
59 29.4 OJ 2.0
1 4.2 VC 0.5
13 15.2 VC 1.0
56 30.9 OJ 2.0
27 26.7 VC 2.0
53 22.4 OJ 2.0
```

```
# Check the structure
str(my_data)
```

```
'data.frame': 60 obs. of 3 variables:
$ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
$ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
$ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
```

From the output above, R considers “dose” as a numeric variable. We’ll convert it as a factor variable (i.e., grouping variable) as follow.

```
# Convert dose as a factor and recode the levels
# as "D0.5", "D1", "D2"
my_data$dose <- factor(my_data$dose,
levels = c(0.5, 1, 2),
labels = c("D0.5", "D1", "D2"))
head(my_data)
```

```
len supp dose
1 4.2 VC D0.5
2 11.5 VC D0.5
3 7.3 VC D0.5
4 5.8 VC D0.5
5 6.4 VC D0.5
6 10.0 VC D0.5
```

Question: We want to know if tooth length depends on supp and dose.

- Generate frequency tables:

`table(my_data$supp, my_data$dose)`

```
D0.5 D1 D2
OJ 10 10 10
VC 10 10 10
```

We have 2X3 design cells with the factors being supp and dose and 10 subjects in each cell. Here, we have a **balanced design**. In the next sections I’ll describe how to analyse data from balanced designs, since this is the simplest case.

Box plots and line plots can be used to visualize group differences:

**Box plot**to plot the data grouped by the combinations of the levels of the two factors.**Two-way interaction plot**, which plots the mean (or other summary) of the response for two-way combinations of factors, thereby illustrating possible interactions.To use R base graphs read this: R base graphs. Here, we’ll use the

**ggpubr**R package for an easy ggplot2-based data visualization.Install the latest version of ggpubr from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

- Visualize your data with ggpubr:

```
# Box plot with multiple groups
# +++++++++++++++++++++
# Plot tooth length ("len") by groups ("dose")
# Color box plot by a second group: "supp"
library("ggpubr")
ggboxplot(my_data, x = "dose", y = "len", color = "supp",
palette = c("#00AFBB", "#E7B800"))
```

```
# Line plots with multiple groups
# +++++++++++++++++++++++
# Plot tooth length ("len") by groups ("dose")
# Color box plot by a second group: "supp"
# Add error bars: mean_se
# (other values include: mean_sd, mean_ci, median_iqr, ....)
library("ggpubr")
ggline(my_data, x = "dose", y = "len", color = "supp",
add = c("mean_se", "dotplot"),
palette = c("#00AFBB", "#E7B800"))
```

If you still want to use R base graphs, type the following scripts:

```
# Box plot with two factor variables
boxplot(len ~ supp * dose, data=my_data, frame = FALSE,
col = c("#00AFBB", "#E7B800"), ylab="Tooth Length")
# Two-way interaction plot
interaction.plot(x.factor = my_data$dose, trace.factor = my_data$supp,
response = my_data$len, fun = mean,
type = "b", legend = TRUE,
xlab = "Dose", ylab="Tooth Length",
pch=c(1,19), col = c("#00AFBB", "#E7B800"))
```

Arguments used for the function **interaction.plot**():

**x.factor**: the factor to be plotted on x axis.**trace.factor**: the factor to be plotted as lines**response**: a numeric variable giving the response**type**: the type of plot. Allowed values include**p**(for point only),**l**(for line only) and**b**(for both point and line).

We want to know if tooth length depends on supp and dose.

The R function **aov**() can be used to answer this question. The function **summary.aov**() is used to summarize the analysis of variance model.

```
res.aov2 <- aov(len ~ supp + dose, data = my_data)
summary(res.aov2)
```

```
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 14.02 0.000429 ***
dose 2 2426.4 1213.2 82.81 < 2e-16 ***
Residuals 56 820.4 14.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

The output includes the columns *F value* and *Pr(>F)* corresponding to the p-value of the test.

From the ANOVA table we can conclude that both *supp* and *dose* are statistically significant. *dose* is the most significant factor variable. These results would lead us to believe that changing delivery methods (supp) or the dose of vitamin C, will impact significantly the mean tooth length.

Not the above fitted model is called **additive model**. It makes an assumption that the two factor variables are independent. If you think that these two variables might interact to create an synergistic effect, replace the plus symbol (+) by an asterisk (*), as follow.

```
# Two-way ANOVA with interaction effect
# These two calls are equivalent
res.aov3 <- aov(len ~ supp * dose, data = my_data)
res.aov3 <- aov(len ~ supp + dose + supp:dose, data = my_data)
summary(res.aov3)
```

```
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

It can be seen that the two main effects (supp and dose) are statistically significant, as well as their interaction.

Note that, in the situation where the interaction is not significant you should use the additive model.

From the ANOVA results, you can conclude the following, based on the p-values and a significance level of 0.05:

- the p-value of supp is 0.000429 (significant), which indicates that the levels of supp are associated with significant different tooth length.
- the p-value of dose is < 2e-16 (significant), which indicates that the levels of dose are associated with significant different tooth length.
- the p-value for the interaction between supp*dose is 0.02 (significant), which indicates that the relationships between dose and tooth length depends on the supp method.

- Compute mean and SD by groups using dplyr R package:

```
require("dplyr")
group_by(my_data, supp, dose) %>%
summarise(
count = n(),
mean = mean(len, na.rm = TRUE),
sd = sd(len, na.rm = TRUE)
)
```

```
Source: local data frame [6 x 5]
Groups: supp [?]
supp dose count mean sd
(fctr) (fctr) (int) (dbl) (dbl)
1 OJ D0.5 10 13.23 4.459709
2 OJ D1 10 22.70 3.910953
3 OJ D2 10 26.06 2.655058
4 VC D0.5 10 7.98 2.746634
5 VC D1 10 16.77 2.515309
6 VC D2 10 26.14 4.797731
```

- It’s also possible to use the function
**model.tables**() as follow:

`model.tables(res.aov3, type="means", se = TRUE)`

In ANOVA test, a significant p-value indicates that some of the group means are different, but we don’t know which pairs of groups are different.

It’s possible to perform multiple pairwise-comparison, to determine if the mean difference between specific pairs of group are statistically significant.

As the ANOVA test is significant, we can compute **Tukey HSD** (Tukey Honest Significant Differences, R function: **TukeyHSD**()) for performing multiple pairwise-comparison between the means of groups. The function **TukeyHD**() takes the fitted ANOVA as an argument.

We don’t need to perform the test for the “supp” variable because it has only two levels, which have been already proven to be significantly different by ANOVA test. Therefore, the Tukey HSD test will be done only for the factor variable “dose”.

`TukeyHSD(res.aov3, which = "dose")`

```
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = len ~ supp + dose + supp:dose, data = my_data)
$dose
diff lwr upr p adj
D1-D0.5 9.130 6.362488 11.897512 0.0e+00
D2-D0.5 15.495 12.727488 18.262512 0.0e+00
D2-D1 6.365 3.597488 9.132512 2.7e-06
```

**diff**: difference between means of the two groups**lwr**,**upr**: the lower and the upper end point of the confidence interval at 95% (default)**p adj**: p-value after adjustment for the multiple comparisons.

It can be seen from the output, that all pairwise comparisons are significant with an adjusted p-value < 0.05.

It’s possible to use the function **glht**() [in **multcomp** package] to perform multiple comparison procedures for an ANOVA. **glht** stands for general linear hypothesis tests. The simplified format is as follow:

`glht(model, lincft)`

**model**: a fitted model, for example an object returned by**aov**().**lincft**(): a specification of the linear hypotheses to be tested. Multiple comparisons in ANOVA models are specified by objects returned from the function**mcp**().

Use glht() to perform multiple pairwise-comparisons:

```
library(multcomp)
summary(glht(res.aov2, linfct = mcp(dose = "Tukey")))
```

```
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = len ~ supp + dose, data = my_data)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
D1 - D0.5 == 0 9.130 1.210 7.543 <1e-05 ***
D2 - D0.5 == 0 15.495 1.210 12.802 <1e-05 ***
D2 - D1 == 0 6.365 1.210 5.259 <1e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
```

The function **pairwise.t.test**() can be also used to calculate pairwise comparisons between group levels with corrections for multiple testing.

```
pairwise.t.test(my_data$len, my_data$dose,
p.adjust.method = "BH")
```

```
Pairwise comparisons using t tests with pooled SD
data: my_data$len and my_data$dose
D0.5 D1
D1 1.0e-08 -
D2 4.4e-16 1.4e-05
P value adjustment method: BH
```

ANOVA assumes that the data are normally distributed and the variance across groups are homogeneous. We can check that with some diagnostic plots.

The **residuals versus fits plot** is used to check the homogeneity of variances. In the plot below, there is no evident relationships between residuals and fitted values (the mean of each groups), which is good. So, we can assume the homogeneity of variances.

```
# 1. Homogeneity of variances
plot(res.aov3, 1)
```

Points 32 and 23 are detected as outliers, which can severely affect normality and homogeneity of variance. It can be useful to remove outliers to meet the test assumptions.

Use the **Levene’s test** to check the **homogeneity of variances**. The function **leveneTest**() [in **car** package] will be used:

```
library(car)
leveneTest(len ~ supp*dose, data = my_data)
```

```
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 1.7086 0.1484
54
```

From the output above we can see that the p-value is not less than the significance level of 0.05. This means that there is no evidence to suggest that the variance across groups is statistically significantly different. Therefore, we can assume the homogeneity of variances in the different treatment groups.

**Normality plot of the residuals**. In the plot below, the quantiles of the residuals are plotted against the quantiles of the normal distribution. A 45-degree reference line is also plotted.

The normal probability plot of residuals is used to verify the assumption that the residuals are normally distributed.

The normal probability plot of the residuals should approximately follow a straight line.

```
# 2. Normality
plot(res.aov3, 2)
```

As all the points fall approximately along this reference line, we can assume normality.

The conclusion above, is supported by the **Shapiro-Wilk test** on the ANOVA residuals (W = 0.98, p = 0.5) which finds no indication that normality is violated.

```
# Extract the residuals
aov_residuals <- residuals(object = res.aov3)
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals )
```

```
Shapiro-Wilk normality test
data: aov_residuals
W = 0.98499, p-value = 0.6694
```

An **unbalanced design** has unequal numbers of subjects in each group.

There are three fundamentally different ways to run an ANOVA in an unbalanced design. They are known as Type-I, Type-II and Type-III sums of squares. To keep things simple, note that The recommended method are the Type-III sums of squares.

The three methods give the same result when the design is balanced. However, when the design is unbalanced, they don’t give the same results.

The function **Anova**() [in **car** package] can be used to compute two-way ANOVA test for unbalanced designs.

First install the package on your computer. In R, type **install.packages**(“car”). Then:

```
library(car)
my_anova <- aov(len ~ supp * dose, data = my_data)
Anova(my_anova, type = "III")
```

```
Anova Table (Type III tests)
Response: len
Sum Sq Df F value Pr(>F)
(Intercept) 1750.33 1 132.730 3.603e-16 ***
supp 137.81 1 10.450 0.002092 **
dose 885.26 2 33.565 3.363e-10 ***
supp:dose 108.32 2 4.107 0.021860 *
Residuals 712.11 54
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

- Analysis of variance (ANOVA, parametric):
- Kruskal-Wallis Test in R (non parametric alternative to one-way ANOVA)

This analysis has been performed using **R software** (ver. 3.2.4).

- What is one-way ANOVA test?
- Assumptions of ANOVA test
- How one-way ANOVA test works?
- Visualize your data and compute one-way ANOVA in R
- Summary
- See also
- Read more
- Infos

The **one-way** **analysis of variance** (**ANOVA**), also known as *one-factor ANOVA*, is an extension of independent two-samples **t-test** for comparing means in a situation where there are more than two groups. In **one-way ANOVA**, the data is organized into several groups base on one single grouping variable (also called *factor* variable). This tutorial describes the basic principle of the **one-way** **ANOVA test** and provides practical **anova test examples** in R software.

**ANOVA test hypotheses**:

- Null hypothesis: the means of the different groups are the same
- Alternative hypothesis: At least one sample mean is not equal to the others.

Note that, if you have only two groups, you can use t-test. In this case the F-test and the t-test are equivalent.

Here we describe the requirement for **ANOVA test**. ANOVA test can be applied only when:

- The observations are obtained independently and randomly from the population defined by the factor levels
- The data of each factor level are normally distributed.
- These normal populations have a common variance. (
**Levene’s test**can be used to check this.)

Assume that we have 3 groups (A, B, C) to compare:

- Compute the common variance, which is called
**variance within samples**(\(S^2_{within}\)) or**residual variance**. - Compute the variance between sample means as follow:
- Compute the mean of each group
- Compute the
**variance between sample means**(\(S^2_{between}\))

- Produce F-statistic as the ratio of \(S^2_{between}/S^2_{within}\).

Note that, a lower ratio (ratio < 1) indicates that there are no significant difference between the means of the samples being compared. However, a higher ratio implies that the variation among group means are significant.

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use the built-in R data set named *PlantGrowth*. It contains the weight of plants obtained under a control and two different treatment conditions.

`my_data <- PlantGrowth`

To have an idea of what the data look like, we use the the function **sample_n**()[in **dplyr** package]. The **sample_n**() function randomly picks a few of the observations in the data frame to print out:

```
# Show a random sample
set.seed(1234)
dplyr::sample_n(my_data, 10)
```

```
weight group
19 4.32 trt1
18 4.89 trt1
29 5.80 trt2
24 5.50 trt2
17 6.03 trt1
1 4.17 ctrl
6 4.61 ctrl
16 3.83 trt1
12 4.17 trt1
15 5.87 trt1
```

In R terminology, the column “group” is called factor and the different categories (“ctr”, “trt1”, “trt2”) are named factor levels. **The levels are ordered alphabetically**.

```
# Show the levels
levels(my_data$group)
```

`[1] "ctrl" "trt1" "trt2"`

If the levels are not automatically in the correct order, re-order them as follow:

```
my_data$group <- ordered(my_data$group,
levels = c("ctrl", "trt1", "trt2"))
```

It’s possible to compute summary statistics (mean and sd) by groups using the **dplyr** package.

- Compute summary statistics by groups - count, mean, sd:

```
library(dplyr)
group_by(my_data, group) %>%
summarise(
count = n(),
mean = mean(weight, na.rm = TRUE),
sd = sd(weight, na.rm = TRUE)
)
```

```
Source: local data frame [3 x 4]
group count mean sd
(fctr) (int) (dbl) (dbl)
1 ctrl 10 5.032 0.5830914
2 trt1 10 4.661 0.7936757
3 trt2 10 5.526 0.4425733
```

To use R base graphs read this: R base graphs. Here, we’ll use the

**ggpubr**R package for an easy ggplot2-based data visualization.Install the latest version of ggpubr from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

- Visualize your data with ggpubr:

```
# Box plots
# ++++++++++++++++++++
# Plot weight by group and color by group
library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
order = c("ctrl", "trt1", "trt2"),
ylab = "Weight", xlab = "Treatment")
```

```
# Mean plots
# ++++++++++++++++++++
# Plot weight by group
# Add error bars: mean_se
# (other values include: mean_sd, mean_ci, median_iqr, ....)
library("ggpubr")
ggline(my_data, x = "group", y = "weight",
add = c("mean_se", "jitter"),
order = c("ctrl", "trt1", "trt2"),
ylab = "Weight", xlab = "Treatment")
```

If you still want to use R base graphs, type the following scripts:

```
# Box plot
boxplot(weight ~ group, data = my_data,
xlab = "Treatment", ylab = "Weight",
frame = FALSE, col = c("#00AFBB", "#E7B800", "#FC4E07"))
# plotmeans
library("gplots")
plotmeans(weight ~ group, data = my_data, frame = FALSE,
xlab = "Treatment", ylab = "Weight",
main="Mean Plot with 95% CI")
```

We want to know if there is any significant difference between the average weights of plants in the 3 experimental conditions.

The R function **aov**() can be used to answer to this question. The function **summary.aov**() is used to summarize the analysis of variance model.

```
# Compute the analysis of variance
res.aov <- aov(weight ~ group, data = my_data)
# Summary of the analysis
summary(res.aov)
```

```
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.766 1.8832 4.846 0.0159 *
Residuals 27 10.492 0.3886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

The output includes the columns *F value* and *Pr(>F)* corresponding to the p-value of the test.

As the p-value is less than the significance level 0.05, we can conclude that there are significant differences between the groups highlighted with “*" in the model summary.

In one-way ANOVA test, a significant p-value indicates that some of the group means are different, but we don’t know which pairs of groups are different.

It’s possible to perform multiple pairwise-comparison, to determine if the mean difference between specific pairs of group are statistically significant.

As the ANOVA test is significant, we can compute **Tukey HSD** (Tukey Honest Significant Differences, R function: **TukeyHSD**()) for performing multiple pairwise-comparison between the means of groups.

The function **TukeyHD**() takes the fitted ANOVA as an argument.

`TukeyHSD(res.aov)`

```
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = weight ~ group, data = my_data)
$group
diff lwr upr p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1 0.865 0.1737839 1.5562161 0.0120064
```

**diff**: difference between means of the two groups**lwr**,**upr**: the lower and the upper end point of the confidence interval at 95% (default)**p adj**: p-value after adjustment for the multiple comparisons.

It can be seen from the output, that only the difference between trt2 and trt1 is significant with an adjusted p-value of 0.012.

It’s possible to use the function **glht**() [in **multcomp** package] to perform multiple comparison procedures for an ANOVA. **glht** stands for general linear hypothesis tests. The simplified format is as follow:

`glht(model, lincft)`

**model**: a fitted model, for example an object returned by**aov**().**lincft**(): a specification of the linear hypotheses to be tested. Multiple comparisons in ANOVA models are specified by objects returned from the function**mcp**().

Use glht() to perform multiple pairwise-comparisons for a one-way ANOVA:

```
library(multcomp)
summary(glht(res.aov, linfct = mcp(group = "Tukey")))
```

```
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = weight ~ group, data = my_data)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
trt1 - ctrl == 0 -0.3710 0.2788 -1.331 0.391
trt2 - ctrl == 0 0.4940 0.2788 1.772 0.198
trt2 - trt1 == 0 0.8650 0.2788 3.103 0.012 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
```

The function **pairewise.t.test**() can be also used to calculate pairwise comparisons between group levels with corrections for multiple testing.

```
pairwise.t.test(my_data$weight, my_data$group,
p.adjust.method = "BH")
```

```
Pairwise comparisons using t tests with pooled SD
data: my_data$weight and my_data$group
ctrl trt1
trt1 0.194 -
trt2 0.132 0.013
P value adjustment method: BH
```

The result is a table of p-values for the pairwise comparisons. Here, the p-values have been adjusted by the Benjamini-Hochberg method.

The ANOVA test assumes that, the data are normally distributed and the variance across groups are homogeneous. We can check that with some diagnostic plots.

The **residuals versus fits plot** can be used to check the homogeneity of variances.

In the plot below, there is no evident relationships between residuals and fitted values (the mean of each groups), which is good. So, we can assume the homogeneity of variances.

```
# 1. Homogeneity of variances
plot(res.aov, 1)
```

Points 17, 15, 4 are detected as outliers, which can severely affect normality and homogeneity of variance. It can be useful to remove outliers to meet the test assumptions.

It’s also possible to use **Bartlett’s test** or **Levene’s test** to check the **homogeneity of variances**.

We recommend **Levene’s test**, which is less sensitive to departures from normal distribution. The function **leveneTest**() [in **car** package] will be used:

```
library(car)
leveneTest(weight ~ group, data = my_data)
```

```
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 1.1192 0.3412
27
```

From the output above we can see that the p-value is not less than the significance level of 0.05. This means that there is no evidence to suggest that the variance across groups is statistically significantly different. Therefore, we can assume the homogeneity of variances in the different treatment groups.

The classical one-way ANOVA test requires an assumption of equal variances for all groups. In our example, the homogeneity of variance assumption turned out to be fine: the Levene test is not significant.

How do we save our ANOVA test, in a situation where the homogeneity of variance assumption is violated?

An alternative procedure (i.e.: **Welch one-way test**), that does not require that assumption have been implemented in the function **oneway.test**().

**ANOVA test with no assumption of equal variances**

`oneway.test(weight ~ group, data = my_data)`

**Pairwise t-tests with no assumption of equal variances**

```
pairwise.t.test(my_data$weight, my_data$group,
p.adjust.method = "BH", pool.sd = FALSE)
```

**Normality plot of residuals**. In the plot below, the quantiles of the residuals are plotted against the quantiles of the normal distribution. A 45-degree reference line is also plotted.

The normal probability plot of residuals is used to check the assumption that the residuals are normally distributed. It should approximately follow a straight line.

```
# 2. Normality
plot(res.aov, 2)
```

As all the points fall approximately along this reference line, we can assume normality.

The conclusion above, is supported by the **Shapiro-Wilk test** on the ANOVA residuals (W = 0.96, p = 0.6) which finds no indication that normality is violated.

```
# Extract the residuals
aov_residuals <- residuals(object = res.aov )
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals )
```

```
Shapiro-Wilk normality test
data: aov_residuals
W = 0.96607, p-value = 0.4379
```

Note that, a non-parametric alternative to one-way ANOVA is **Kruskal-Wallis** **rank sum test**, which can be used when ANNOVA assumptions are not met.

`kruskal.test(weight ~ group, data = my_data)`

```
Kruskal-Wallis rank sum test
data: weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842
```

- Import your data from a .txt tab file:
*my_data <- read.delim(file.choose())*. Here, we used*my_data <- PlantGrowth*. - Visualize your data:
*ggpubr::ggboxplot(my_data, x = “group”, y = “weight”, color = “group”)* - Compute one-way ANOVA test:
*summary(aov(weight ~ group, data = my_data))* - Tukey multiple pairwise-comparisons:
*TukeyHSD(res.aov)*

- Analysis of variance (ANOVA, parametric):
- Kruskal-Wallis Test in R (non parametric alternative to one-way ANOVA)

- (Quick-R: ANOVA/MANOVA)[http://www.statmethods.net/stats/anova.html]
- (Quick-R: (M)ANOVA Assumptions)[http://www.statmethods.net/stats/anovaAssumptions.html]
- (R and Analysis of Variance)[http://personality-project.org/r/r.guide/r.anova.html

This analysis has been performed using **R software** (ver. 3.2.4).

Previously, we described the essentials of R programming and provided quick start guides for importing data into **R**. Additionally, we described how to compute descriptive or summary statistics and correlation analysis using R software.

This chapter contains articles describing **statistical tests** to use for **comparing means**. These tests include:

**T-test****Wilcoxon test****ANOVA test**and**Kruskal-Wallis test**

- Comparing one-sample mean to a standard known mean:
- Comparing the means of two independent groups:
- Comparing the means of paired samples:
- Comparing the means of more than two groups
- Analysis of variance (ANOVA, parametric):
- Kruskal-Wallis Test in R (non parametric alternative to one-way ANOVA)

- What is one-sample t-test?
- Research questions and statistical hypotheses
- Formula of one-sample t-test
- Visualize your data and compute one-sample t-test in R
- R function to compute one-sample t-test
- Visualize your data using box plots
- Preliminary test to check one-sample t-test assumptions
- Compute one-sample t-test
- Interpretation of the result

Read more: —> One-Sample T-test.

- What’s one-sample Wilcoxon signed rank test?
- Research questions and statistical hypotheses
- Visualize your data and compute one-sample Wilcoxon test in R
- R function to compute one-sample Wilcoxon test
- Visualize your data using box plots
- Compute one-sample Wilcoxon test

Read more: —> One-Sample Wilcoxon Test (non-parametric).

- What is unpaired two-samples t-test?
- Research questions and statistical hypotheses
- Formula of unpaired two-samples t-test
- Visualize your data and compute unpaired two-samples t-test in R
- R function to compute unpaired two-samples t-test
- Visualize your data using box plots
- Preliminary test to check independent t-test assumptions
- Compute unpaired two-samples t-test

- Interpretation of the result

Read more: —> Unpaired Two Samples T-test (parametric).

- R function to compute Wilcoxon test
- Visualize your data using box plots
- Compute unpaired two-samples Wilcoxon test

Read more: —> Unpaired Two-Samples Wilcoxon Test (non-parametric).

Read more: —> Paired Samples Wilcoxon Test (non-parametric).

An extension of independent two-samples **t-test** for comparing means in a situation where there are more than two groups.

- What is one-way ANOVA test?
- Assumptions of ANOVA test
- How one-way ANOVA test works?
- Visualize your data and compute one-way ANOVA in R
- Visualize your data
- Compute one-way ANOVA test
- Interpret the result of one-way ANOVA tests
- Multiple pairwise-comparison between the means of groups
- Tukey multiple pairewise-comparisons
- Multiple comparisons using multcomp package
- Pairwise t-test

- Check ANOVA assumptions: test validity?
- Check the homogeneity of variance assumption
- Relaxing the homogeneity of variance assumption
- Check the normality assumption

- Non-parametric alternative to one-way ANOVA test

Read more: —> One-Way ANOVA Test in R.

- What is two-way ANOVA test?
- Two-way ANOVA test hypotheses
- Assumptions of two-way ANOVA test
- Compute two-way ANOVA test in R: balanced designs
- Visualize your data
- Compute two-way ANOVA test
- Interpret the results
- Compute some summary statistics
- Multiple pairwise-comparison between the means of groups
- Tukey multiple pairewise-comparisons
- Multiple comparisons using multcomp package
- Pairwise t-test

- Check ANOVA assumptions: test validity?
- Check the homogeneity of variance assumption

- Check the normality assumption

- Compute two-way ANOVA test in R for unbalanced designs

Read more: —> Two-Way ANOVA Test in R.

- What is MANOVA test?
- Assumptions of MANOVA
- Interpretation of MANOVA
- Compute MANOVA in R

Read more: —> MANOVA Test in R: Multivariate Analysis of Variance.

- What is Kruskal-Wallis test?
- Visualize your data and compute Kruskal-Wallis test in R
- Visualize the data using box plots
- Compute Kruskal-Wallis test
- Multiple pairwise-comparison between groups

Read more: —> Kruskal-Wallis Test in R (non parametric alternative to one-way ANOVA).

This analysis has been performed using **R statistical software** (ver. 3.2.4).

In the situation where there multiple response variables you can test them simultaneously using a **multivariate analysis of variance** (**MANOVA**). This article describes how to compute **manova** in R.

For example, we may conduct an experiment where we give two treatments (A and B) to two groups of mice, and we are interested in the weight and height of mice. In that case, the weight and height of mice are two dependent variables, and our hypothesis is that both together are affected by the difference in treatment. A multivariate analysis of variance could be used to test this hypothesis.

MANOVA can be used in certain conditions:

The dependent variables should be normally distribute within groups. The R function

**mshapiro.test**( )[in the**mvnormtest**package] can be used to perform the Shapiro-Wilk test for multivariate normality. This is useful in the case of MANOVA, which assumes**multivariate normality**.Homogeneity of variances across the range of predictors.

Linearity between all pairs of dependent variables, all pairs of covariates, and all dependent variable-covariate pairs in each cell

If the global multivariate test is significant, we conclude that the corresponding effect (treatment) is significant. In that case, the next question is to determine if the treatment affects only the weight, only the height or both. In other words, we want to identify the specific dependent variables that contributed to the significant global effect.

To answer this question, we can use one-way ANOVA (or univariate ANOVA) to examine separately each dependent variable.

**Prepare your data**as specified here: [url=/wiki/best-practices-for-preparing-your-data-set-for-r]Best practices for preparing your data set for R[/url]**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use iris data set:

```
# Store the data in the variable my_data
my_data <- iris
```

The R code below display a random sample of our data using the function **sample_n**()[in **dplyr** package]. First, install dplyr if you don’t have it:

`install.packages("dplyr")`

```
# Show a random sample
set.seed(1234)
dplyr::sample_n(my_data, 10)
```

```
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
94 5.0 2.3 3.3 1.0 versicolor
91 5.5 2.6 4.4 1.2 versicolor
93 5.8 2.6 4.0 1.2 versicolor
127 6.2 2.8 4.8 1.8 virginica
150 5.9 3.0 5.1 1.8 virginica
2 4.9 3.0 1.4 0.2 setosa
34 5.5 4.2 1.4 0.2 setosa
96 5.7 3.0 4.2 1.2 versicolor
74 6.1 2.8 4.7 1.2 versicolor
98 6.2 2.9 4.3 1.3 versicolor
```

Question: We want to know if there is any significant difference, in *sepal* and *petal* length, between the different species.

The function **manova**() can be used as follow:

```
sepl <- iris$Sepal.Length
petl <- iris$Petal.Length
# MANOVA test
res.man <- manova(cbind(Sepal.Length, Petal.Length) ~ Species, data = iris)
summary(res.man)
```

```
Df Pillai approx F num Df den Df Pr(>F)
Species 2 0.9885 71.829 4 294 < 2.2e-16 ***
Residuals 147
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

```
# Look to see which differ
summary.aov(res.man)
```

```
Response Sepal.Length :
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 63.212 31.606 119.26 < 2.2e-16 ***
Residuals 147 38.956 0.265
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Response Petal.Length :
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 437.10 218.551 1180.2 < 2.2e-16 ***
Residuals 147 27.22 0.185
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

From the output above, it can be seen that the two variables are highly significantly different among Species.

- Analysis of variance (ANOVA, parametric):
- [url=/wiki/one-way-anova-test-in-r]One-Way ANOVA Test in R[/url]
- [url=/wiki/two-way-anova-test-in-r]Two-Way ANOVA Test in R[/url]

- [url=/wiki/kruskal-wallis-test-in-r]Kruskal-Wallis Test in R (non parametric alternative to one-way ANOVA)[/url]

This analysis has been performed using **R software** (ver. 3.2.4).

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use the built-in R data set named *PlantGrowth*. It contains the weight of plants obtained under a control and two different treatment conditions.

`my_data <- PlantGrowth`

```
# print the head of the file
head(my_data)
```

```
weight group
1 4.17 ctrl
2 5.58 ctrl
3 5.18 ctrl
4 6.11 ctrl
5 4.50 ctrl
6 4.61 ctrl
```

In R terminology, the column “group” is called factor and the different categories (“ctr”, “trt1”, “trt2”) are named factor levels. **The levels are ordered alphabetically**.

```
# Show the group levels
levels(my_data$group)
```

`[1] "ctrl" "trt1" "trt2"`

If the levels are not automatically in the correct order, re-order them as follow:

```
my_data$group <- ordered(my_data$group,
levels = c("ctrl", "trt1", "trt2"))
```

It’s possible to compute summary statistics by groups. The dplyr package can be used.

- To install
**dplyr**package, type this:

`install.packages("dplyr")`

- Compute summary statistics by groups:

```
library(dplyr)
group_by(my_data, group) %>%
summarise(
count = n(),
mean = mean(weight, na.rm = TRUE),
sd = sd(weight, na.rm = TRUE),
median = median(weight, na.rm = TRUE),
IQR = IQR(weight, na.rm = TRUE)
)
```

```
Source: local data frame [3 x 6]
group count mean sd median IQR
(fctr) (int) (dbl) (dbl) (dbl) (dbl)
1 ctrl 10 5.032 0.5830914 5.155 0.7425
2 trt1 10 4.661 0.7936757 4.550 0.6625
3 trt2 10 5.526 0.4425733 5.435 0.4675
```

To use R base graphs read this: R base graphs. Here, we’ll use the

**ggpubr**R package for an easy ggplot2-based data visualization.Install the latest version of ggpubr from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

- Visualize your data with ggpubr:

```
# Box plots
# ++++++++++++++++++++
# Plot weight by group and color by group
library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
order = c("ctrl", "trt1", "trt2"),
ylab = "Weight", xlab = "Treatment")
```

```
# Mean plots
# ++++++++++++++++++++
# Plot weight by group
# Add error bars: mean_se
# (other values include: mean_sd, mean_ci, median_iqr, ....)
library("ggpubr")
ggline(my_data, x = "group", y = "weight",
add = c("mean_se", "jitter"),
order = c("ctrl", "trt1", "trt2"),
ylab = "Weight", xlab = "Treatment")
```

We want to know if there is any significant difference between the average weights of plants in the 3 experimental conditions.

The test can be performed using the function **kruskal.test**() as follow:

`kruskal.test(weight ~ group, data = my_data)`

```
Kruskal-Wallis rank sum test
data: weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842
```

As the p-value is less than the significance level 0.05, we can conclude that there are significant differences between the treatment groups.

From the output of the Kruskal-Wallis test, we know that there is a significant difference between groups, but we don’t know which pairs of groups are different.

It’s possible to use the function **pairwise.wilcox.test**() to calculate pairwise comparisons between group levels with corrections for multiple testing.

```
pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group,
p.adjust.method = "BH")
```

```
Pairwise comparisons using Wilcoxon rank sum test
data: PlantGrowth$weight and PlantGrowth$group
ctrl trt1
trt1 0.199 -
trt2 0.095 0.027
P value adjustment method: BH
```

The pairwise comparison shows that, only trt1 and trt2 are significantly different (p < 0.05).

- Analysis of variance (ANOVA, parametric):

This analysis has been performed using **R software** (ver. 3.2.4).

The **paired samples Wilcoxon test** (also known as **Wilcoxon signed-rank test**) is a **non-parametric** alternative to paired t-test used to compare paired data. It’s used when your data are not normally distributed. This tutorial describes how to compute paired samples Wilcoxon test in **R**.

Differences between paired samples should be distributed symmetrically around the median.

The R function **wilcox.test**() can be used as follow:

`wilcox.test(x, y, paired = TRUE, alternative = "two.sided")`

**x,y**: numeric vectors**paired**: a logical value specifying that we want to compute a paired Wilcoxon test**alternative**: the alternative hypothesis. Allowed value is one of “two.sided” (default), “greater” or “less”.

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use an example data set, which contains the weight of 10 mice before and after the treatment.

```
# Data in two numeric vectors
# ++++++++++++++++++++++++++
# Weight of the mice before treatment
before <-c(200.1, 190.9, 192.7, 213, 241.4, 196.9, 172.2, 185.5, 205.2, 193.7)
# Weight of the mice after treatment
after <-c(392.9, 393.2, 345.1, 393, 434, 427.9, 422, 383.9, 392.3, 352.2)
# Create a data frame
my_data <- data.frame(
group = rep(c("before", "after"), each = 10),
weight = c(before, after)
)
```

We want to know, if there is any significant difference in the median weights before and after treatment?

```
# Print all data
print(my_data)
```

```
group weight
1 before 200.1
2 before 190.9
3 before 192.7
4 before 213.0
5 before 241.4
6 before 196.9
7 before 172.2
8 before 185.5
9 before 205.2
10 before 193.7
11 after 392.9
12 after 393.2
13 after 345.1
14 after 393.0
15 after 434.0
16 after 427.9
17 after 422.0
18 after 383.9
19 after 392.3
20 after 352.2
```

Compute summary statistics (median and inter-quartile range (IQR)) by groups using the dplyr package can be used.

- Install
**dplyr**package:

`install.packages("dplyr")`

- Compute summary statistics by groups:

```
library("dplyr")
group_by(my_data, group) %>%
summarise(
count = n(),
median = median(weight, na.rm = TRUE),
IQR = IQR(weight, na.rm = TRUE)
)
```

```
Source: local data frame [2 x 4]
group count median IQR
(fctr) (int) (dbl) (dbl)
1 after 10 392.95 28.800
2 before 10 195.30 12.575
```

**ggpubr**R package for an easy ggplot2-based data visualization.Install the latest version of ggpubr from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

- Visualize your data:

```
# Plot weight by group and color by group
library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("#00AFBB", "#E7B800"),
order = c("before", "after"),
ylab = "Weight", xlab = "Groups")
```

Box plots show you the increase, but lose the paired information. You can use the function **plot.paired**() [in **pairedData** package] to plot paired data (“before - after” plot).

- Install pairedData package:

`install.packages("PairedData")`

- Plot paired data:

```
# Subset weight data before treatment
before <- subset(my_data, group == "before", weight,
drop = TRUE)
# subset weight data after treatment
after <- subset(my_data, group == "after", weight,
drop = TRUE)
# Plot paired data
library(PairedData)
pd <- paired(before, after)
plot(pd, type = "profile") + theme_bw()
```

Question : Is there any significant changes in the weights of mice before after treatment?

**1) Compute paired Wilcoxon test - Method 1**: The data are saved in two different numeric vectors.

```
res <- wilcox.test(before, after, paired = TRUE)
res
```

```
Wilcoxon signed rank test
data: before and after
V = 0, p-value = 0.001953
alternative hypothesis: true location shift is not equal to 0
```

**2) Compute paired Wilcoxon-test - Method 2**: The data are saved in a data frame.

```
# Compute t-test
res <- wilcox.test(weight ~ group, data = my_data, paired = TRUE)
res
```

```
Wilcoxon signed rank test
data: weight by group
V = 55, p-value = 0.001953
alternative hypothesis: true location shift is not equal to 0
```

```
# print only the p-value
res$p.value
```

`[1] 0.001953125`

As you can see, the two methods give the same results.

The **p-value** of the test is 0.001953, which is less than the significance level alpha = 0.05. We can conclude that the median weight of the mice before treatment is significantly different from the median weight after treatment with a **p-value** = 0.001953.

Note that:

- if you want to test whether the median weight before treatment is less than the median weight after treatment, type this:

```
wilcox.test(weight ~ group, data = my_data, paired = TRUE,
alternative = "less")
```

- Or, if you want to test whether the median weight before treatment is greater than the median weight after treatment, type this

```
wilcox.test(weight ~ group, data = my_data, paired = TRUE,
alternative = "greater")
```

You can perform **paired-sample Wilcoxon test**, **online**, without any installation by clicking the following link:

This analysis has been performed using **R software** (ver. 3.2.4).

The **paired samples t-test** is used to compare the means between two related groups of samples. In this case, you have two values (i.e., pair of values) for the same samples. This article describes how to compute **paired samples t-test** using **R software**.

As an example of data, 20 mice received a treatment X during 3 months. We want to know whether the treatment X has an impact on the weight of the mice.

To answer to this question, the weight of the 20 mice has been measured before and after the treatment. This gives us 20 sets of values before treatment and 20 sets of values after treatment from measuring twice the weight of the **same mice**.

In such situations, **paired t-test** can be used to compare the mean weights before and after treatment.

Paired t-test analysis is performed as follow:

- Calculate the difference (\(d\)) between each pair of value
- Compute the mean (\(m\)) and the standard deviation (\(s\)) of \(d\)
- Compare the average difference to 0. If there is any significant difference between the two pairs of samples, then the mean of d (\(m\)) is expected to be far from 0.

Paired t-test can be used only when the difference \(d\) is normally distributed. This can be checked using **Shapiro-Wilk test**.

Typical research questions are:

- whether the mean difference (\(m\))
*is equal*to 0? - whether the mean difference (\(m\))
*is less than*0? - whether the mean difference (\(m\))
*is greather than*0?

In statistics, we can define the corresponding *null hypothesis* (\(H_0\)) as follow:

- \(H_0: m = 0\)
- \(H_0: m \leq 0\)
- \(H_0: m \geq 0\)

The corresponding *alternative hypotheses* (\(H_a\)) are as follow:

- \(H_a: m \ne 0\) (different)
- \(H_a: m > 0\) (greater)
- \(H_a: m < 0\) (less)

Note that:

- Hypotheses 1) are called
**two-tailed tests** - Hypotheses 2) and 3) are called
**one-tailed tests**

**t-test statistisc value** can be calculated using the following formula:

\[ t = \frac{m}{s/\sqrt{n}} \]

where,

**m**is the mean differences**n**is the sample size (i.e., size of d).**s**is the standard deviation of d

We can compute the p-value corresponding to the absolute value of the **t-test statistics** (|t|) for the **degrees of freedom** (df): \(df = n - 1\).

If the p-value is inferior or equal to 0.05, we can conclude that the difference between the two paired samples are significantly different.

To perform paired samples t-test comparing the means of two paired samples (x & y), the R function **t.test**() can be used as follow:

`t.test(x, y, paired = TRUE, alternative = "two.sided")`

**x,y**: numeric vectors**paired**: a logical value specifying that we want to compute a paired t-test**alternative**: the alternative hypothesis. Allowed value is one of “two.sided” (default), “greater” or “less”.

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use an example data set, which contains the weight of 10 mice before and after the treatment.

```
# Data in two numeric vectors
# ++++++++++++++++++++++++++
# Weight of the mice before treatment
before <-c(200.1, 190.9, 192.7, 213, 241.4, 196.9, 172.2, 185.5, 205.2, 193.7)
# Weight of the mice after treatment
after <-c(392.9, 393.2, 345.1, 393, 434, 427.9, 422, 383.9, 392.3, 352.2)
# Create a data frame
my_data <- data.frame(
group = rep(c("before", "after"), each = 10),
weight = c(before, after)
)
```

We want to know, if there is any significant difference in the mean weights after treatment?

```
# Print all data
print(my_data)
```

```
group weight
1 before 200.1
2 before 190.9
3 before 192.7
4 before 213.0
5 before 241.4
6 before 196.9
7 before 172.2
8 before 185.5
9 before 205.2
10 before 193.7
11 after 392.9
12 after 393.2
13 after 345.1
14 after 393.0
15 after 434.0
16 after 427.9
17 after 422.0
18 after 383.9
19 after 392.3
20 after 352.2
```

Compute summary statistics (mean and sd) by groups using the dplyr package.

- To install
**dplyr**package, type this:

`install.packages("dplyr")`

- Compute summary statistics by groups:

```
library("dplyr")
group_by(my_data, group) %>%
summarise(
count = n(),
mean = mean(weight, na.rm = TRUE),
sd = sd(weight, na.rm = TRUE)
)
```

```
Source: local data frame [2 x 4]
group count mean sd
(fctr) (int) (dbl) (dbl)
1 after 10 393.65 29.39801
2 before 10 199.16 18.47354
```

**ggpubr** R package for an easy ggplot2-based data visualization.

- Install the latest version of ggpubr from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

- Visualize your data:

```
# Plot weight by group and color by group
library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("#00AFBB", "#E7B800"),
order = c("before", "after"),
ylab = "Weight", xlab = "Groups")
```

Box plots show you the increase, but lose the paired information. You can use the function **plot.paired**() [in **pairedData** package] to plot paired data (“before - after” plot).

- Install pairedData package:

`install.packages("PairedData")`

- Plot paired data:

```
# Subset weight data before treatment
before <- subset(my_data, group == "before", weight,
drop = TRUE)
# subset weight data after treatment
after <- subset(my_data, group == "after", weight,
drop = TRUE)
# Plot paired data
library(PairedData)
pd <- paired(before, after)
plot(pd, type = "profile") + theme_bw()
```

Assumption 1: Are the two samples paired?

Yes, since the data have been collected from measuring twice the weight of the same mice.

Assumption 2: Is this a large sample?

No, because n < 30. Since the sample size is not large enough (less than 30), we need to check whether the differences of the pairs follow a normal distribution.

How to check the normality?

Use Shapiro-Wilk normality test as described at: Normality Test in R.

- Null hypothesis: the data are normally distributed
- Alternative hypothesis: the data are not normally distributed

```
# compute the difference
d <- with(my_data,
weight[group == "before"] - weight[group == "after"])
# Shapiro-Wilk normality test for the differences
shapiro.test(d) # => p-value = 0.6141
```

From the output, the p-value is greater than the significance level 0.05 implying that the distribution of the differences (d) are not significantly different from normal distribution. In other words, we can assume the normality.

Note that, if the data are not normally distributed, it’s recommended to use the non parametric **paired two-samples Wilcoxon test**.

Question : Is there any significant changes in the weights of mice after treatment?

**1) Compute paired t-test - Method 1**: The data are saved in two different numeric vectors.

```
# Compute t-test
res <- t.test(before, after, paired = TRUE)
res
```

```
Paired t-test
data: before and after
t = -20.883, df = 9, p-value = 6.2e-09
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-215.5581 -173.4219
sample estimates:
mean of the differences
-194.49
```

**2) Compute paired t-test - Method 2**: The data are saved in a data frame.

```
# Compute t-test
res <- t.test(weight ~ group, data = my_data, paired = TRUE)
res
```

```
Paired t-test
data: weight by group
t = 20.883, df = 9, p-value = 6.2e-09
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
173.4219 215.5581
sample estimates:
mean of the differences
194.49
```

As you can see, the two methods give the same results.

In the result above :

**t**is the**t-test statistic**value (t = 20.88),**df**is the degrees of freedom (df= 9),**p-value**is the significance level of the**t-test**(p-value = 6.210^{-9}).**conf.int**is the**confidence interval**(conf.int) of the mean differences at 95% is also shown (conf.int= [173.42, 215.56])**sample estimates**is the mean differences between pairs (mean = 194.49).

Note that:

- if you want to test whether the average weight before treatment is less than the average weight after treatment, type this:

```
t.test(weight ~ group, data = my_data, paired = TRUE,
alternative = "less")
```

- Or, if you want to test whether the average weight before treatment is greater than the average weight after treatment, type this

```
t.test(weight ~ group, data = my_data, paired = TRUE,
alternative = "greater")
```

The **p-value** of the test is 6.210^{-9}, which is less than the significance level alpha = 0.05. We can then reject null hypothesis and conclude that the average weight of the mice before treatment is significantly different from the average weight after treatment with a **p-value** = 6.210^{-9}.

The result of **t.test()** function is a list containing the following components:

**statistic**: the value of the**t test statistics****parameter**: the**degrees of freedom**for the**t test statistics****p.value**: the**p-value**for the test**conf.int**: a**confidence interval**for the mean appropriate to the specified**alternative hypothesis**.**estimate**: the means of the two groups being compared (in the case of**independent t test**) or difference in means (in the case of**paired t test**).

The format of the **R** code to use for getting these values is as follow:

```
# printing the p-value
res$p.value
```

`[1] 6.200298e-09`

```
# printing the mean
res$estimate
```

```
mean of the differences
194.49
```

```
# printing the confidence interval
res$conf.int
```

```
[1] 173.4219 215.5581
attr(,"conf.level")
[1] 0.95
```

You can perform **paired-samples t-test**, **online**, without any installation by clicking the following link:

This analysis has been performed using **R software** (ver. 3.2.4).

The **unpaired two-samples Wilcoxon test** (also known as **Wilcoxon rank sum test** or **Mann-Whitney** test) is a non-parametric alternative to the unpaired two-samples t-test, which can be used to compare two independent groups of samples. It’s used when your data are not normally distributed.

This article describes how to compute two samples Wilcoxon test in R.

To perform two-samples Wilcoxon test comparing the means of two independent samples (x & y), the R function **wilcox.test**() can be used as follow:

`wilcox.test(x, y, alternative = "two.sided")`

**x,y**: numeric vectors**alternative**: the alternative hypothesis. Allowed value is one of “two.sided” (default), “greater” or “less”.

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use an example data set, which contains the weight of 18 individuals (9 women and 9 men):

```
# Data in two numeric vectors
women_weight <- c(38.9, 61.2, 73.3, 21.8, 63.4, 64.6, 48.4, 48.8, 48.5)
men_weight <- c(67.8, 60, 63.4, 76, 89.4, 73.3, 67.3, 61.3, 62.4)
# Create a data frame
my_data <- data.frame(
group = rep(c("Woman", "Man"), each = 9),
weight = c(women_weight, men_weight)
)
```

We want to know, if the median women’s weight differs from the median men’s weight?

`print(my_data)`

```
group weight
1 Woman 38.9
2 Woman 61.2
3 Woman 73.3
4 Woman 21.8
5 Woman 63.4
6 Woman 64.6
7 Woman 48.4
8 Woman 48.8
9 Woman 48.5
10 Man 67.8
11 Man 60.0
12 Man 63.4
13 Man 76.0
14 Man 89.4
15 Man 73.3
16 Man 67.3
17 Man 61.3
18 Man 62.4
```

It’s possible to compute summary statistics (median and interquartile range (IQR)) by groups. The dplyr package can be used.

- To install
**dplyr**package, type this:

`install.packages("dplyr")`

- Compute summary statistics by groups:

```
library(dplyr)
group_by(my_data, group) %>%
summarise(
count = n(),
median = median(weight, na.rm = TRUE),
IQR = IQR(weight, na.rm = TRUE)
)
```

```
Source: local data frame [2 x 4]
group count median IQR
(fctr) (int) (dbl) (dbl)
1 Man 9 67.3 10.9
2 Woman 9 48.8 15.0
```

You can draw R base graphs as described at this link: R base graphs. Here, we’ll use the **ggpubr** R package for an easy ggplot2-based data visualization

- Install the latest version of ggpubr from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

- Visualize your data:

```
# Plot weight by group and color by group
library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("#00AFBB", "#E7B800"),
ylab = "Weight", xlab = "Groups")
```

Question : Is there any significant difference between women and men weights?

**1) Compute two-samples Wilcoxon test - Method 1**: The data are saved in two different numeric vectors.

```
res <- wilcox.test(women_weight, men_weight)
res
```

```
Wilcoxon rank sum test with continuity correction
data: women_weight and men_weight
W = 15, p-value = 0.02712
alternative hypothesis: true location shift is not equal to 0
```

It will give a warning message, saying that “cannot compute exact p-value with tie”. It comes from the assumption of a Wilcoxon test that the responses are continuous. You can suppress this message by adding another argument **exact = FALSE**, but the result will be the same.

**2) Compute two-samples Wilcoxon test - Method 2**: The data are saved in a data frame.

```
res <- wilcox.test(weight ~ group, data = my_data,
exact = FALSE)
res
```

```
Wilcoxon rank sum test with continuity correction
data: weight by group
W = 66, p-value = 0.02712
alternative hypothesis: true location shift is not equal to 0
```

```
# Print the p-value only
res$p.value
```

`[1] 0.02711657`

As you can see, the two methods give the same results.

The **p-value** of the test is 0.02712, which is less than the significance level alpha = 0.05. We can conclude that men’s median weight is significantly different from women’s median weight with a **p-value** = 0.02712.

Note that:

- if you want to test whether the median men’s weight is less than the median women’s weight, type this:

```
wilcox.test(weight ~ group, data = my_data,
exact = FALSE, alternative = "less")
```

- Or, if you want to test whether the median men’s weight is greater than the median women’s weight, type this

```
wilcox.test(weight ~ group, data = my_data,
exact = FALSE, alternative = "greater")
```

You can perform unpaired **two-samples Wilcoxon test**, **online**, without any installation by clicking the following link:

- Compare one-sample mean to a standard known mean
- Compare the means of two independent groups

This analysis has been performed using **R software** (ver. 3.2.4).

- What is unpaired two-samples t-test?
- Research questions and statistical hypotheses
- Formula of unpaired two-samples t-test
- Visualize your data and compute unpaired two-samples t-test in R
- Install ggpubr R package for data visualization
- R function to compute unpaired two-samples t-test
- Import your data into R
- Check your data
- Visualize your data using box plots
- Preleminary test to check independent t-test assumptions
- Compute unpaired two-samples t-test
- Interpretation of the result
- Access to the values returned by t.test() function

- Online unpaired two-samples t-test calculator
- See also
- Infos

The **unpaired two-samples t-test** is used to compare the **mean** of two independent groups.

For example, suppose that we have measured the weight of 100 individuals: 50 women (group A) and 50 men (group B). We want to know if the mean weight of women (\(m_A\)) is significantly different from that of men (\(m_B\)).

In this case, we have two unrelated (i.e., independent or unpaired) groups of samples. Therefore, it’s possible to use an **independent t-test** to evaluate whether the means are different.

Note that, unpaired two-samples t-test can be used only under certain conditions:

- when the two groups of samples (A and B), being compared, are
**normally distributed**. This can be checked using**Shapiro-Wilk test**. - and when the
**variances**of the two groups are equal. This can be checked using**F-test**.

This article describes the formula of the independent t-test and provides pratical examples in R.

Typical research questions are:

- whether the mean of group A (\(m_A\))
*is equal*to the mean of group B (\(m_B\))? - whether the mean of group A (\(m_A\))
*is less than*the mean of group B (\(m_B\))? - whether the mean of group A (\(m_A\))
*is greather than*the mean of group B (\(m_B\))?

In statistics, we can define the corresponding *null hypothesis* (\(H_0\)) as follow:

- \(H_0: m_A = m_B\)
- \(H_0: m_A \leq m_B\)
- \(H_0: m_A \geq m_B\)

The corresponding *alternative hypotheses* (\(H_a\)) are as follow:

- \(H_a: m_A \ne m_B\) (different)
- \(H_a: m_A > m_B\) (greater)
- \(H_a: m_A < m_B\) (less)

Note that:

- Hypotheses 1) are called
**two-tailed tests** - Hypotheses 2) and 3) are called
**one-tailed tests**

**Classical t-test**:

If the variance of the two groups are equivalent (**homoscedasticity**), the **t-test value**, comparing the two samples (\(A\) and \(B\)), can be calculated as follow.

\[ t = \frac{m_A - m_B}{\sqrt{ \frac{S^2}{n_A} + \frac{S^2}{n_B} }} \]

where,

- \(m_A\) and \(m_B\) represent the mean value of the group A and B, respectively.
- \(n_A\) and \(n_B\) represent the sizes of the group A and B, respectively.
- \(S^2\) is an estimator of the pooled
**variance**of the two groups. It can be calculated as follow :

\[ S^2 = \frac{\sum{(x-m_A)^2}+\sum{(x-m_B)^2}}{n_A+n_B-2} \]

with degrees of freedom (df): \(df = n_A + n_B - 2\).

2.**Welch t-statistic**:

If the variances of the two groups being compared are different (**heteroscedasticity**), it’s possible to use the **Welch t test**, an adaptation of Student t-test.

**Welch t-statistic** is calculated as follow :

\[ t = \frac{m_A - m_B}{\sqrt{ \frac{S_A^2}{n_A} + \frac{S_B^2}{n_B} }} \]

where, \(S_A\) and \(S_B\) are the standard deviation of the the two groups A and B, respectively.

Unlike the classic Student’s t-test, **Welch t-test formula** involves the variance of each of the two groups (\(S_A^2\) and \(S_B^2\)) being compared. In other words, it does not use the pooled variance\(S\).

The **degrees of freedom** of **Welch t-test** is estimated as follow :

\[ df = (\frac{S_A^2}{n_A}+ \frac{S_B^2}{n_B^2}) / (\frac{S_A^4}{n_A^2(n_B-1)} + \frac{S_B^4}{n_B^2(n_B-1)} ) \]

A p-value can be computed for the corresponding absolute value of t-statistic (|t|).

Note that, the Welch t-test is considered as the safer one. Usually, the results of the **classical t-test** and the **Welch t-test** are very similar unless both the group sizes and the standard deviations are very different.

How to interpret the results?

If the p-value is inferior or equal to the significance level 0.05, we can reject the null hypothesis and accept the alternative hypothesis. In other words, we can conclude that the mean values of group A and B are significantly different.

You can draw R base graphs as described at this link: R base graphs. Here, we’ll use the **ggpubr** R package for an easy ggplot2-based data visualization

- Install the latest version from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

To perform two-samples t-test comparing the means of two independent samples (x & y), the R function **t.test**() can be used as follow:

`t.test(x, y, alternative = "two.sided", var.equal = FALSE)`

**x,y**: numeric vectors**alternative**: the alternative hypothesis. Allowed value is one of “two.sided” (default), “greater” or “less”.**var.equal**: a logical variable indicating whether to treat the two variances as being equal. If TRUE then the pooled variance is used to estimate the variance otherwise the Welch test is used.

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use an example data set, which contains the weight of 18 individuals (9 women and 9 men):

```
# Data in two numeric vectors
women_weight <- c(38.9, 61.2, 73.3, 21.8, 63.4, 64.6, 48.4, 48.8, 48.5)
men_weight <- c(67.8, 60, 63.4, 76, 89.4, 73.3, 67.3, 61.3, 62.4)
# Create a data frame
my_data <- data.frame(
group = rep(c("Woman", "Man"), each = 9),
weight = c(women_weight, men_weight)
)
```

We want to know, if the average women’s weight differs from the average men’s weight?

```
# Print all data
print(my_data)
```

```
group weight
1 Woman 38.9
2 Woman 61.2
3 Woman 73.3
4 Woman 21.8
5 Woman 63.4
6 Woman 64.6
7 Woman 48.4
8 Woman 48.8
9 Woman 48.5
10 Man 67.8
11 Man 60.0
12 Man 63.4
13 Man 76.0
14 Man 89.4
15 Man 73.3
16 Man 67.3
17 Man 61.3
18 Man 62.4
```

It’s possible to compute summary statistics (mean and sd) by groups. The dplyr package can be used.

- To install
**dplyr**package, type this:

`install.packages("dplyr")`

- Compute summary statistics by groups:

```
library(dplyr)
group_by(my_data, group) %>%
summarise(
count = n(),
mean = mean(weight, na.rm = TRUE),
sd = sd(weight, na.rm = TRUE)
)
```

```
Source: local data frame [2 x 4]
group count mean sd
(fctr) (int) (dbl) (dbl)
1 Man 9 68.98889 9.375426
2 Woman 9 52.10000 15.596714
```

```
# Plot weight by group and color by group
library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("#00AFBB", "#E7B800"),
ylab = "Weight", xlab = "Groups")
```

**Assumption 1**: Are the two samples independents?

Yes, since the samples from men and women are not related.

**Assumtion 2**: Are the data from each of the 2 groups follow a normal distribution?

Use Shapiro-Wilk normality test as described at: Normality Test in R. - Null hypothesis: the data are normally distributed - Alternative hypothesis: the data are not normally distributed

We’ll use the functions **with**() and **shapiro.test**() to compute Shapiro-Wilk test for each group of samples.

```
# Shapiro-Wilk normality test for Men's weights
with(my_data, shapiro.test(weight[group == "Man"]))# p = 0.1
# Shapiro-Wilk normality test for Women's weights
with(my_data, shapiro.test(weight[group == "Woman"])) # p = 0.6
```

From the output, the two p-values are greater than the significance level 0.05 implying that the distribution of the data are not significantly different from the normal distribution. In other words, we can assume the normality.

Note that, if the data are not normally distributed, it’s recommended to use the non parametric two-samples Wilcoxon rank test.

**Assumption 3**. Do the two populations have the same variances?

We’ll use **F-test** to test for homogeneity in variances. This can be performed with the function **var.test**() as follow:

```
res.ftest <- var.test(weight ~ group, data = my_data)
res.ftest
```

```
F test to compare two variances
data: weight by group
F = 0.36134, num df = 8, denom df = 8, p-value = 0.1714
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.08150656 1.60191315
sample estimates:
ratio of variances
0.3613398
```

**F-test** is p = 0.1713596. It’s greater than the significance level alpha = 0.05. In conclusion, there is no significant difference between the variances of the two sets of data. Therefore, we can use the classic **t-test** witch assume equality of **the two variances**.

Question : Is there any significant difference between women and men weights?

**1) Compute independent t-test - Method 1**: The data are saved in two different numeric vectors.

```
# Compute t-test
res <- t.test(women_weight, men_weight, var.equal = TRUE)
res
```

```
Two Sample t-test
data: women_weight and men_weight
t = -2.7842, df = 16, p-value = 0.01327
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-29.748019 -4.029759
sample estimates:
mean of x mean of y
52.10000 68.98889
```

**2) Compute independent t-test - Method 2**: The data are saved in a data frame.

```
# Compute t-test
res <- t.test(weight ~ group, data = my_data, var.equal = TRUE)
res
```

```
Two Sample t-test
data: weight by group
t = 2.7842, df = 16, p-value = 0.01327
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
4.029759 29.748019
sample estimates:
mean in group Man mean in group Woman
68.98889 52.10000
```

As you can see, the two methods give the same results.

In the result above :

**t**is the**t-test statistic**value (t = 2.784),**df**is the degrees of freedom (df= 16),**p-value**is the significance level of the**t-test**(p-value = 0.01327).**conf.int**is the**confidence interval**of the mean at 95% (conf.int = [4.0298, 29.748]);**sample estimates**is he mean value of the sample (mean = 68.9888889, 52.1).

Note that:

- if you want to test whether the average men’s weight is less than the average women’s weight, type this:

```
t.test(weight ~ group, data = my_data,
var.equal = TRUE, alternative = "less")
```

- Or, if you want to test whether the average men’s weight is greater than the average women’s weight, type this

```
t.test(weight ~ group, data = my_data,
var.equal = TRUE, alternative = "greater")
```

The **p-value** of the test is 0.01327, which is less than the significance level alpha = 0.05. We can conclude that men’s average weight is significantly different from women’s average weight with a **p-value** = 0.01327.

The result of **t.test()** function is a list containing the following components:

**statistic**: the value of the**t test statistics****parameter**: the**degrees of freedom**for the**t test statistics****p.value**: the**p-value**for the test**conf.int**: a**confidence interval**for the mean appropriate to the specified**alternative hypothesis**.**estimate**: the means of the two groups being compared (in the case of**independent t test**) or difference in means (in the case of**paired t test**).

The format of the **R** code to use for getting these values is as follow:

```
# printing the p-value
res$p.value
```

`[1] 0.0132656`

```
# printing the mean
res$estimate
```

```
mean in group Man mean in group Woman
68.98889 52.10000
```

```
# printing the confidence interval
res$conf.int
```

```
[1] 4.029759 29.748019
attr(,"conf.level")
[1] 0.95
```

You can perform unpaired **two-samples t-test**, **online**, without any installation by clicking the following link:

- Compare one-sample mean to a standard known mean
- Compare the means of two independent groups

This analysis has been performed using **R software** (ver. 3.2.4).

- What is one-sample t-test?
- Research questions and statistical hypotheses
- Formula of one-sample t-test
- Visualize your data and compute one-sample t-test in R
- Install ggpubr R package for data visualization
- R function to compute one-sample t-test
- Import your data into R
- Check your data
- Visualize your data using box plots
- Preleminary test to check one-sample t-test assumptions
- Compute one-sample t-test
- Interpretation of the result
- Access to the values returned by t.test() function

- Online one-sample t-test calculator
- See also
- Infos

Generally, the theoretical mean comes from:

- a previous experiment. For example, compare whether the mean weight of mice differs from 200 mg, a value determined in a previous study.
- or from an experiment where you have control and treatment conditions. If you express your data as “percent of control”, you can test whether the average value of treatment condition differs significantly from 100.

Note that, one-sample t-test can be used only, when the data are normally distributed . This can be checked using **Shapiro-Wilk test** .

Typical research questions are:

- whether the mean (\(m\)) of the sample
*is equal*to the theoretical mean (\(\mu\))? - whether the mean (\(m\)) of the sample
*is less than*the theoretical mean (\(\mu\))? - whether the mean (\(m\)) of the sample
*is greater than*the theoretical mean (\(\mu\))?

In statistics, we can define the corresponding *null hypothesis* (\(H_0\)) as follow:

- \(H_0: m = \mu\)
- \(H_0: m \leq \mu\)
- \(H_0: m \geq \mu\)

The corresponding *alternative hypotheses* (\(H_a\)) are as follow:

- \(H_a: m \ne \mu\) (different)
- \(H_a: m > \mu\) (greater)
- \(H_a: m < \mu\) (less)

Note that:

- Hypotheses 1) are called
**two-tailed tests** - Hypotheses 2) and 3) are called
**one-tailed tests**

The t-statistic can be calculated as follow:

\[ t = \frac{m-\mu}{s/\sqrt{n}} \]

where,

**m**is the sample**mean****n**is the sample**size****s**is the sample**standard deviation**with \(n-1\)**degrees of freedom**- \(\mu\) is the
**theoretical value**

We can compute the p-value corresponding to the absolute value of the **t-test statistics** (|t|) for the **degrees of freedom** (df): \(df = n - 1\).

How to interpret the results?

If the p-value is inferior or equal to the significance level 0.05, we can reject the null hypothesis and accept the alternative hypothesis. In other words, we conclude that the sample mean is significantly different from the theoretical mean.

You can draw R base graps as described at this link: R base graphs. Here, we’ll use the **ggpubr** R package for an easy ggplot2-based data visualization

- Install the latest version from GitHub as follow (recommended):

```
# Install
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
```

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

To perform one-sample t-test, the R function **t.test**() can be used as follow:

`t.test(x, mu = 0, alternative = "two.sided")`

**x**: a numeric vector containing your data values**mu**: the theoretical mean. Default is 0 but you can change it.**alternative**: the alternative hypothesis. Allowed value is one of “two.sided” (default), “greater” or “less”.

**Prepare your data**as specified here: Best practices for preparing your data set for R**Save your data**in an external .txt tab or .csv files**Import your data into R**as follow:

```
# If .txt tab file, use this
my_data <- read.delim(file.choose())
# Or, if .csv file, use this
my_data <- read.csv(file.choose())
```

Here, we’ll use an example data set containing the weight of 10 mice.

We want to know, if the average weight of the mice differs from 25g?

```
set.seed(1234)
my_data <- data.frame(
name = paste0(rep("M_", 10), 1:10),
weight = round(rnorm(10, 20, 2), 1)
)
```

```
# Print the first 10 rows of the data
head(my_data, 10)
```

```
name weight
1 M_1 17.6
2 M_2 20.6
3 M_3 22.2
4 M_4 15.3
5 M_5 20.9
6 M_6 21.0
7 M_7 18.9
8 M_8 18.9
9 M_9 18.9
10 M_10 18.2
```

```
# Statistical summaries of weight
summary(my_data$weight)
```

```
Min. 1st Qu. Median Mean 3rd Qu. Max.
15.30 18.38 18.90 19.25 20.82 22.20
```

**Min.**: the minimum value**1st Qu.**: The first quartile. 25% of values are lower than this.**Median**: the median value. Half the values are lower; half are higher.**3rd Qu.**: the third quartile. 75% of values are higher than this.**Max.**: the maximum value

```
library(ggpubr)
ggboxplot(my_data$weight,
ylab = "Weight (g)", xlab = FALSE,
ggtheme = theme_minimal())
```

**Is this a large sample**? - No, because n < 30.- Since the sample size is not large enough (less than 30, central limit theorem), we need to
**check whether the data follow a normal distribution**.

How to check the normality?

Read this article: Normality Test in R.

Briefly, it’s possible to use the **Shapiro-Wilk normality test** and to look at the **normality plot**.

**Shapiro-Wilk test**:- Null hypothesis: the data are normally distributed
- Alternative hypothesis: the data are not normally distributed

`shapiro.test(my_data$weight) # => p-value = 0.6993`

From the output, the p-value is greater than the significance level 0.05 implying that the distribution of the data are not significantly different from normal distribtion. In other words, we can assume the normality.

**Visual inspection**of the data normality using**Q-Q plots**(quantile-quantile plots). Q-Q plot draws the correlation between a given sample and the normal distribution.

```
library("ggpubr")
ggqqplot(my_data$weight, ylab = "Men's weight",
ggtheme = theme_minimal())
```

From the normality plots, we conclude that the data may come from normal distributions.

Note that, if the data are not normally distributed, it’s recommended to use the non parametric one-sample Wilcoxon rank test.

We want to know, if the average weight of the mice differs from 25g (two-tailed test)?

```
# One-sample t-test
res <- t.test(my_data$weight, mu = 25)
# Printing the results
res
```

```
One Sample t-test
data: my_data$weight
t = -9.0783, df = 9, p-value = 7.953e-06
alternative hypothesis: true mean is not equal to 25
95 percent confidence interval:
17.8172 20.6828
sample estimates:
mean of x
19.25
```

In the result above :

**t**is the**t-test statistic**value (t = -9.078),**df**is the degrees of freedom (df= 9),**p-value**is the significance level of the**t-test**(p-value = 7.95310^{-6}).**conf.int**is the**confidence interval**of the mean at 95% (conf.int = [17.8172, 20.6828]);**sample estimates**is he mean value of the sample (mean = 19.25).

Note that:

- if you want to test whether the mean weight of mice is less than 25g (one-tailed test), type this:

```
t.test(my_data$weight, mu = 25,
alternative = "less")
```

- Or, if you want to test whether the mean weight of mice is greater than 25g (one-tailed test), type this:

```
t.test(my_data$weight, mu = 25,
alternative = "greater")
```

The **p-value** of the test is 7.95310^{-6}, which is less than the significance level alpha = 0.05. We can conclude that the mean weight of the mice is significantly different from 25g with a **p-value** = 7.95310^{-6}.

The result of **t.test()** function is a list containing the following components:

**statistic**: the value of the**t test statistics****parameter**: the**degrees of freedom**for the**t test statistics****p.value**: the**p-value**for the test**conf.int**: a**confidence interval**for the mean appropriate to the specified**alternative hypothesis**.**estimate**: the means of the two groups being compared (in the case of**independent t test**) or difference in means (in the case of**paired t test**).

The format of the **R** code to use for getting these values is as follow:

```
# printing the p-value
res$p.value
```

`[1] 7.953383e-06`

```
# printing the mean
res$estimate
```

```
mean of x
19.25
```

```
# printing the confidence interval
res$conf.int
```

```
[1] 17.8172 20.6828
attr(,"conf.level")
[1] 0.95
```

You can perform **one-sample t-test**, **online**, without any installation by clicking the following link:

This analysis has been performed using **R software** (ver. 3.2.4).