# Hypertension chez les souris
# La tension a été mesurée chez un ensemble de souris issues d’un rétrocroisement de deux
# souches de souris. On s’intéresse aux souris possédant le marqueur D4Mit214. On souhaite
# connaitre l’effet du marqueur sur la pression sanguine (systolique). Pour cela on note la
# tension des souris (BA = hétérozygote) et BB (homozygotes) pour ce marqueur.

# hétérozygotes (BA) 26 mesures
a = c(86, 88, 89, 89, 92, 93, 94, 94, 94, 95, 95, 96, 96, 97, 97, 98, 98, 99, 99, 101, 106, 107, 110, 113, 116, 118)


# homozygotes (BB) 24 mesures
 b = c(89, 90, 92, 93, 93, 96, 99, 99, 99, 102, 103, 104, 105, 106, 106, 107, 108, 108, 110, 110, 112, 114, 116, 116)

# On combine les deux vecteurs (a) et (b) en un seul (comb), car sous hypothèse nulle (H0)
# il n’existe pas de différences entre les deux groupes. C’est donc sans effet si H0 vraie.
comb = c(a,b)

# La difference observee entre les moyennes
diff.random = NULL       # on met à 0 une variable que l’on va utiliser
nombre_de_permutations = 1000

for (i in 1 : nombre_de_permutations)
{
	#  On brasse par tirage SANS remise l’ensemble comb, pour ensuite reconstruire a et b
	brasse = sample (comb, length(comb), replace = FALSE)

	a.random = brasse[1 : length(a)]
	b.random = brasse[(length(a) + 1) : length(comb)]

	# Difference entre les échantillons permutés
	diff.random[i] = mean(b.random) - mean(a.random)
}

#la moyenne des différences (en valeur absolue) entre les échantillons permutés
mean(abs(diff.random))
# à comparer avec :
diff.observe = mean(b) - mean(a)
print(diff.observe)
#L’ensemble des différences par permutations:
print(diff.random)


hist(diff.random, breaks=10, freq= FALSE, col="lightgreen",
      xlim=c(-10,10),  ylim=c(0, .20),main="permutations, 10 classes")
curve(dnorm(x, mean=mean(diff.random), sd=sd(diff.random)), add=TRUE, col="darkblue", lwd=2) 


# On calcule simplement le ratio du nombre de fois que la différence entre permutés est
# supérieure ou égale à celle observée sur le nombre total de permutations. C’est-à-dire la
# probabilité que l’on aie une différence plus grande que celle observée ALORS QUE l’on
# suppose H0 vraie. Si cette valeur est petite (< 0.05), on peut rejeter H0.

p = sum(abs(diff.random) >= abs(diff.observe)) / nombre_de_permutations
print(p)

# Alternative : trier diff.random et regarder si diff.observe est contenu (ou non) dans 
# l’intervalle [2.5, 97.5] de pourcentile de diff.random. 
# On rejette(accepte) H0 au seuil de 5%.

p <- c(0.025, 0.975)
quantile(sort(diff.random),p)
cat('rappel (diff observee) :', diff.observe,"\n")

#---------------------------------------------
# utilisation d'une librairie dediee
#---------------------------------------------

# Installer le package coin, puis 
# Charger le fichier « 2ech indep »

Indep <- read.table("2ech_indep.txt", header=TRUE, sep="\t", dec=",")
library(coin)
tapply(Indep$Mesures, Indep$Lots, median, na.rm=TRUE)
oneway_test(Mesures ~ Lots, alternative='two.sided', data=Indep)
 
#"""----------------EOF-----------------------------""""