```{r setup, include=FALSE, size="huge"}
library(knitr)
## Default parameters for displaying the slides
knitr::opts_chunk$set(
  echo = TRUE, 
  eval=TRUE, 
  fig.width = 7, 
  fig.height = 5, 
  fig.align = "center", 
  fig.path = "figures/",
  size = "tiny", 
  warning = FALSE, 
  results = TRUE, 
  message = FALSE, 
  comment = "")

dir.main <- "~/stat1"
dir.current <- file.path(dir.main, "practicals", "02_yeast_annotations")
setwd(dir.current)

```


# But de ce TP

```{r echo=FALSE}
## Parameters for the exercise
R <- 10000     ## Number of replicates
alpha <- 0.05  ## Accepted risk for type I error
```

Au cours de ce TP nous effectuerons des tests de comparaison de moyennes sur deux types de données. 

1. **Données artificielles** générées selon des distributions normales, soit sous hypothèse nulle ($H_0$) soit sous hypothèse alternative ($H_1$).  Ceci nous permettra de réaliser des tests dans des situations où nous connaissons les  paramètres des populations à comparer ($\mu_1$, $\mu_2$, $\sigma_1$, $\sigma_2$), en connaissant donc la réponse correcte du test. Le but de cet exercice sera de 

    - nous familiariser avec les tests de comparaison de moyenne: choix d'un test en fonction des caractéristiques des données, choix des paramètres du test, interprétation des résultats; 
    - évaluer l'adéquation des tests en fonction des types de données;
    - mesurer empiriquement les taux d'erreurs de types I et II, et vérifier s'ils correspondent aux attentes théoriques.

2. **Données de transcriptome** obtenues au moyen de puces à ADN (DNA microarrays) chez 190 patients souffrant de leucémies lymphoblastiques aigues (**LLA**), classifiés en différents groupes selon les perturbations chromosomiques (hyperploidie) ou génétiques (mutations d'un gène particulier) supposées être à l'origine du cancer. 

Le but de cette analyse sera de détecter les gènes exprimés différentiellement (***differentially expressed genes***, **DEG**) entre deux sous-types particuliers de cancer. 

****************************************
# Jeux de données artificiels

Pour les jeux de données artificiels, nous poserons arbitrairement un seuil $\alpha = `r alpha`$ sur le risque d'erreur de première espèce. 

## Sous hypothèse nulle

Dans un premier temps, nous allons délibérément générer des données sous hypothèse nulle ($H_0$) c'est-à-dire en tirant des échantillons dans deux populations de taille égale. Les données seront générées dans les **conditions d'applicabilité du test de Student** (***hypothèses de travail***, c'est-à-dore hypothèses préalables à la réalisation du test): les 2 populations dont les échantillons sont extraits sont supposées:

1. suivre des distributions normales (hypothèse de normalité);
2. avoir des variances identques (hypothèse d'homoscédasticité). 

Le but de l'exercice sera de mesurer le **taux de faux posiitfs**, c'est-à-dire la proportion des tests déclarés (à tort) positifs, alors qu'il n'existe pas de différence entre les moyennes des deux populations.

### 1. Génération des données

Au moyen de la fonction `rnorm()`, tirez deux échantillons de taille $n_1 = n_2 = 10$ à partir de populations normales de moyennes $\mu_1 = \mu_2 = 7$ et d'écart-types $\sigma_1 = \sigma_2 = 2$. 

```{r random_numbers}
n1 <- 10 ## Taille du premier échantillon
n2 <- 10 ## Taille du second échantillojn
mu1 <- 7 ## Moyenne de la 1ère population
mu2 <- 7 
sigma1 <- 2 ## Ecart-type de la popumation 1
sigma2 <- 2

## On tire des échantillons des deux populations respectives
x1 <- rnorm(n = n1, mean = mu1, sd = sigma1)
x2 <- rnorm(n = n2, mean = mu2, sd = sigma2)
```

<!--
Calculez 

- les moyennes des deux échantillons ($\bar{x_1}$, $\bar{x_2}$),
- les écart-types ($s_1$, $s_2$) des échantillons (en utilisant les fonctions `sum()`et `sqrt()`), 
- les estimateurs de maximum de vraisemblance des moyennes et écarts-types des populations ($\hat{\mu}_1$, $\hat{\mu}_2$, $\hat{\sigma}_1$, $\hat{\sigma}_1$), en utilisant soit `sum()` soit `sd()`.

```{r mean_and_sd}
## Moyennes d'échantillons
m1 <- mean(x1)
m2 <- mean(x2)

## Ecarts-types d'échantillons
s1 <- sqrt(sum((x1 - m1)^2)/n1)
s2 <- sqrt(sum((x2 - m2)^2)/n2)


## Estimateurs de moyennes des populations
mu1.est <- m1
mu2.est <- m2

## Estimateurs des écarts-types d'échantillons
sigma1.est <- sqrt(sum((x1 - m1)^2)/(n1-1))
sigma2.est <- sqrt(sum((x2 - m2)^2)/(n2-1))


sigma1.est == sd(x1)
sigma2.est == sd(x2)

```

-->
## 2. Application du test avec la fonction `t.test()`


Effectuez un test de comparaison de moyenne en utilisant la fonction `t.test()`. Lisez attentivement l'aide de cette fonction, et choisissez les paramètres en tenant compte des caractéristiques de vos données. Nous effectuerons ici un test bilatéral (*two-tailed*). 

Interprétez le résultat de la fonction `t.test()`  (paramètres, décision, interprétation de la p-valeur). 

```{r}

## Seuil critique sur le risque d'erreur de type I
alpha <- 0.05

## Utiliser la fonction t.test() pour réaliser le test de Student
t.test(x1, x2, 
       var.equal = TRUE,  ## On postule les variances égales
       alternative = "two.sided", ## Test bilatéral
       conf.level = 1 - alpha ## Niveau de confiance
       )

```


### 3. Calcul manuel des statistiques de Student

a. Au moyen de la fonction R `sum()`, calculez les paramètres de vos échantillons nécessaires au test de Student ($\bar{x}_1, \bar{x}_2, s_1, s_2$). 

```{r sample_parameters}

n1 <- length(x1) ## Taille d'échantillon
m1 <- sum(x1) / n1 ## Moenne du premier échantillon
s1 <- sqrt(sum((x1 - m1)^2)/n1) ## Ecart-type du premier échantillon

## Idem pour le second échantillon
n2 <- length(x2) 
m2 <- sum(x2) / n2
s2 <- sqrt(sum((x2 - m2)^2)/n2)

## Afficher les paramètres des échantillons
cat("n1 =", n1, "; m1 =", round(digits=2, m1), "; s1 = ", round(digits=2, s1))
cat("n2 =", n2, "; m2 =", round(digits=2, m2), "; s2 = ", round(digits=2, s2))

```


b. Au moyen de cette même fonction `sum()`, estimez les estimateurs de ces paramètres pour les populations  ($\hat{\mu}_1, \hat{\mu}_2, \hat{\sigma}_1, \hat{\sigma}_2$).

```{r}
## Estimation des moyennes de population
mu1.est <- m1
mu2.est <- m2

## Estimation des écarts-types
sigma1.est <- s1 * sqrt(n1/(n1-1))
sigma2.est <- s2 * sqrt(n2/(n2-1))

## Afficher les paramètres estimés pour la population
cat("mu1.est = ", round(digits=2, mu1.est), "sigma1.est", round(digits=3, sigma1.est))
cat("mu2.est = ", round(digits=2, mu2.est), "sigma2.est", round(digits=3, sigma2.est))

```


c. Recalculez ces paramètres au moyen des fonctions R `mean()` et `sd()`.

```{r}
mean1 <- mean(x1)
mean2 <- mean(x2)
sd1 <- sd(x1)
sd2 <- sd(x2)

cat("mean1 = ", round(digits=2, mean1), "sd1", round(digits=2, sd1))
cat("mean2 = ", round(digits=2, mean2), "sd2", round(digits=2, sd2))

## Print the result for the knitr report
result <- t(data.frame(
  "moyennes de populations" = c(mu1, mu2),
  "moyennes d'échantillons" = c(m1, m2),
  "moyennes estimées" = c(mu1.est, mu2.est),
  "écarts-types de population" = c(sigma1, sigma2),
  "écarts-types d'échantillons" = c(s1, s2),
  "écarts-types estimés" = c(sigma1.est, sigma2.est)))
colnames(result) <- c("Population 1", "Population 2")

kable(as.data.frame(result), digits = 2, caption = "Estimation des paramètres de centralité (moyenne) et dispersion (écart-type) des populations à partir des échantillons. ")


```


**Attention**, les fonctions R `var()` et `sd()` ne calculent pas les paramètres d'échantillons ($s^2$, $s$) mais les estimateurs des paramètres de population  correspondants ($\hat{\sigma}^2$, $\hat{\sigma}$). 

d. Sur base de ces paramètres, calculez la statistique $t_{obs}$ de Student. Justifiez le choix des paramètres que vous choisirez pour cette fonction. 

$$t_{S} = \frac{\hat{\delta}}{\hat{\sigma}_\delta} =  \frac{\bar{x}_{2} - \bar{x}_{1}}{\sqrt{\frac{n_1 s_{1}^2 + n_2 s_{2}^2}{n_1+n_2-2} \left(\frac{1}{n_1}+ \frac{1}{n_2}\right)}}$$

```{r}

## Difference between sample means
diff <- m1 - m2

## Estimation of the standard error on the difference between sample means
Student.diff.err <- sqrt((1/n1 + 1/n2) * (n1 * s1^2 + n2 * s2^2) / (n1 + n2 - 2))

## Student statistics
Student.t <- diff / Student.diff.err


## Print Student statistics with 5 significant digits
cat("Student statistics  t =", signif(digits=5, Student.t))

```


e. Avec la fonction $pt()$, calculez la p-valeur de cette valeur $t$. 

```{r}

```


f. Interprétez les résultats (décision, interprétation de la p-valeur). 


### 4. Mesure empirique du taux de faux-positifs


Dans cet exercice, nous allons réaliser un grand nombre de tests de Student en nous plaçant sous hypothèse nulle, et compter le nombre de tests retournant une réponse positive. 

a. Avant de commencer l'expérience, indiquez le nombre de faux positifs attendus *a prior* si l'one ffectue $R = `r R`$ tests sous hypothèse nulle, avec un seuil critique de $\alpha= `r  alpha`$.

b. Répétez `r R` fois le test de Student au moyen de la fonction `t.test()` et récupérez dans deux vecteurs séparés les valeurs rapportées pour la statistique $t$ et pour la $p$ valeur. 


c. Dessinez l'histogramme des valeurs $t$ obtenues empiriquement.

d. Calculez la proportion de faux-positifs. Cette proportion correspond-elle à vos attentes ?


## 5. Test sous hypothèse alternative

Effectuez $R = `r R`$ tests de comparaison de moyenne sur des échantillons aléatoires tirés dans des populations de moyennes respectives $\mu_1 = 6$ et $\mu_2  8$, ayant toutes deux un écart-type $\sigma=1$. 

a. Choisissez le test et ses paramètres en fonction des caractéristiqus de vos données. 

b. Comptez le nombre de résultats déclarés positifs et négatifs avec un seuil $\alpha = `r alpha`$.

c. Interprétez ce résultat.

## 6. Test avec variances inégales

Effectuez les mêmes tests (`r R` répliques, d'abord sous hypothèse nulle puis sous hypothèse alternative) avec des données tirées de population de variances inégales: $\sigma_1^2 = 4$, $\sigma_2^2 = 25$. 

a. Justiifez le choix du test et des paramètres. 

b. Comptez le nombre de résultats déclarés positifs ou négatifs.

c. Interprétez le résultat. 

## 7. Test avec données non-normales

Au moyen de la fonction R `rpois()`, générez des données selon une loi de Poisson dont l'espérance vaut $\lambda_1 = 4$ pour la première population et $\lambda_2 = 6$ pour la seconde. 

a. Choisissez le test approprié.

b. Mesurez les proportions de tests respectivement déclarés positifs et négatifs. 

c. Sur ces mêmes données, effectuez un test paramétrique de comparaison de moyennes (Student ou Welch, à vous de choisi en le justifiant). 

d. Interprétez le résultat. 

****************************************
# Données transcriptomiques obtenues par puces à ADN

Si vous arrivez ici, suivez ce tutoriel:

<http://pedagogix-tagc.univ-mrs.fr/ASG/practicals/microarrays_student_test/DenBoer_Student_test.html>

****************************************
# Formules mathématiques

- Les syboles grecs ($\mu$, $\sigma$) correspondent aux statistiques de population, les symboles romains ($\bar{x}$, $s$) aux statistiques d'échantillon.
- L'accent circonflexe ($\hat{ }$) indique les estimateurs de paramètres de population calculés à partir de paramètres d'échantillons. 

| Symbol       | Description   |
|--------------|----------------------------------------------------|
| $\mu_{1},  \mu_{2}$ | Moyennes respectives des populations 1 et 2. |
| $\sigma_{1}, \sigma_{2}$ | Ecarts-types respectifs des populations 1 and 2. |
| $N_1$, $N_2$ | Tailles (nombre d'individus) des populations 1 et 2. |
| $n_1$, $n_2$ | Effectifs (nombre d'individus) des échantillons prélevés sur les populations 1 et 2. |
| $\bar{x}_{1}, \bar{x}_{2}$ | Moyennes d'échantillons. |
| $\delta = \mu_{2} - \mu{1} | Différence entre les moyennes des populations. |
| $d = \hat{\delta} = \hat{\mu}_2 - \hat{\mu}_1  = \bar{x}_2 - \bar{x}_1$ | $d$ = **Taille d'effet**: dans un test de comparaison de moyennes, il s'agit de la différence entre les moyennes d'échantillons, utilisée comme estimateur de $\delta$. |
| $s^2_{1}, s^2_{2}$ | Variances mesurées sur les échantillons. |
| $\hat{\sigma}_p = \sqrt{\frac{n_1 s_1^2 + n_2 s_2^2}{n_1+n_2-2}}$ | Ecart-type groupé (*pooled standard deviation*), utilisé comme estimateur de l'écart-type commun des deux populations, en supposant leurs variances égales (hypothèse de travail d'homoscédasticité). |
| $\hat{\sigma}_\delta = \hat{\sigma}_p \sqrt{\left(\frac{1}{n_1}+ \frac{1}{n_2}\right)}$  | Erreur standard sur la différence entre moyennes, en supposant que les populations ont la même variance (test de Student). |
| $t_{S} = \frac{\hat{\delta}}{\hat{\sigma}_\delta} =  \frac{\bar{x}_{2} - \bar{x}_{1}}{\sqrt{\frac{n_1 s_{1}^2 + n_2 s_{2}^2}{n_1+n_2-2} \left(\frac{1}{n_1}+ \frac{1}{n_2}\right)}}$ | statistique $t$ de Student |
| $t_{W}=\frac{\bar{x}_{1} - \bar{x}_{2}}{\sqrt{\frac {s^2_{1}}{n_1} + \frac{s^2_{2}}{n_2}}}$ | statistique $t$ de Welch |
| | |

