**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).

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).

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
```

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:

```
# 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 **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).

The **one-sample Wilcoxon signed rank test** is a non-parametric alternative to **one-sample t-test** when the data cannot be assumed to be normally distributed. It’s used to determine whether the median of the sample is equal to a known standard value (i.e. theoretical value).

Note that, the data should be distributed symmetrically around the median. In other words, there should be roughly the same number of values above and below the median.

Typical research questions are:

- whether the median (\(m\)) of the sample
*is equal*to the theoretical value (\(m_0\))? - whether the median (\(m\)) of the sample
*is less than*to the theoretical value (\(m_0\))? - whether the median (\(m\)) of the sample
*is greater than*to the theoretical value(\(m_0\))?

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

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

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

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

Note that:

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

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 one-sample Wilcoxon-test, the R function **wilcox.test**() can be used as follow:

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

**x**: a numeric vector containing your data values**mu**: the theoretical mean/median value. 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 median 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())
```

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

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

```
Wilcoxon signed rank test with continuity correction
data: my_data$weight
V = 0, p-value = 0.005793
alternative hypothesis: true location is not equal to 25
```

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

`[1] 0.005793045`

The **p-value** of the test is 0.005793, which is less than the significance level alpha = 0.05. We can reject the null hypothesis and conclude that the average weight of the mice is significantly different from 25g with a **p-value** = 0.005793.

Note that:

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

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

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

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

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).

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
```

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:

```
# 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).

- 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.**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).

Many of statistical tests including correlation, regression, t-test, and analysis of variance (ANOVA) assume some certain characteristics about the data. They require the data to follow a **normal distribution** or **Gaussian distribution**. These tests are called **parametric tests**, because their validity depends on the distribution of the data.

Normality and the other assumptions made by these tests should be taken seriously to draw reliable interpretation and conclusions of the research.

Before using a parametric test, we should perform some **preleminary tests** to make sure that the test assumptions are met. In the situations where the assumptions are violated, **non-paramatric** tests are recommended.

Here, we’ll describe how to check the normality of the data by visual inspection and by significance tests.

**dplyr**for data manipulation

`install.packages("dplyr")`

**ggpubr**for an easy ggplot2-based data visualization

- Install the latest version from GitHub as follow:

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

- Or, install from CRAN as follow:

`install.packages("ggpubr")`

```
library("dplyr")
library("ggpubr")
```

**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.

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

We start by displaying a random sample of 10 rows using the function **sample_n**()[in **dplyr** package].

Show 10 random rows:

```
set.seed(1234)
dplyr::sample_n(my_data, 10)
```

```
len supp dose
7 11.2 VC 0.5
37 8.2 OJ 0.5
36 10.0 OJ 0.5
58 27.3 OJ 2.0
49 14.5 OJ 1.0
57 26.4 OJ 2.0
1 4.2 VC 0.5
13 15.2 VC 1.0
35 14.5 OJ 0.5
27 26.7 VC 2.0
```

We want to test if the variable *len* (tooth length) is normally distributed.

If the sample size is large enough (n > 30), we can ignore the distribution of the data and use parametric tests.

**The central limit theorem** tells us that no matter what distribution things have, the sampling distribution tends to be normal if the sample is large enough (n > 30).

However, to be consistent, normality can be checked by visual inspection [**normal plots (histogram)**, **Q-Q plot** (quantile-quantile plot)] or by **significance tests**].

**Density plot** and **Q-Q plot** can be used to check normality visually.

**Density plot**: the**density**plot provides a visual judgment about whether the distribution is bell shaped.

```
library("ggpubr")
ggdensity(my_data$len,
main = "Density plot of tooth length",
xlab = "Tooth length")
```

**Q-Q plot**:**Q-Q plot**(or quantile-quantile plot) draws the correlation between a given sample and the normal distribution. A 45-degree reference line is also plotted.

```
library(ggpubr)
ggqqplot(my_data$len)
```

It’s also possible to use the function **qqPlot**() [in **car** package]:

```
library("car")
qqPlot(my_data$len)
```

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

Visual inspection, described in the previous section, is usually unreliable. It’s possible to use a **significance test** comparing the sample distribution to a normal one in order to ascertain whether data show or not a serious deviation from normality.

There are several methods for **normality test** such as **Kolmogorov-Smirnov (K-S) normality test** and **Shapiro-Wilk’s test**.

The null hypothesis of these tests is that “sample distribution is normal”. If the test is **significant**, the distribution is non-normal.

**Shapiro-Wilk’s method** is widely recommended for normality test and it provides better power than K-S. It is based on the correlation between the data and the corresponding normal scores.

Note that, normality test is sensitive to sample size. Small samples most often pass normality tests. Therefore, it’s important to combine visual inspection and significance test in order to take the right decision.

The R function **shapiro.test**() can be used to perform the Shapiro-Wilk test of normality for one variable (univariate):

`shapiro.test(my_data$len)`

```
Shapiro-Wilk normality test
data: my_data$len
W = 0.96743, p-value = 0.1091
```

From the output, the p-value > 0.05 implying that the distribution of the data are not significantly different from normal distribution. In other words, we can assume the normality.

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