####   TEST DE PERMUTACIONES ######



## Introducimos los datos en R

Treat<-c(69, 24,63,87.5,77.5,40)
Cont<-c(9,12,36,77.5,-7.5,32.5)
Ratbp<-data.frame(Treat, Cont)
Blood<-c(Treat,Cont)

source("c:/......./Combinations.R") ##### Leemos el fichero fuente de combinaciones del sitio web 




pdT6<-t(Combinations(12,6))  ##### Matriz de combinaciones

# Alternativamente cargamos la librería BSDA de Alan


library("BSDA")
pdT6<-SRS(1:12,6) # o  pdT6<-t(Combinations(12,6))

## matriz de datos repetidos con 924 filas:

Comb<-matrix(rep(c(Treat,Cont),924),ncol=12,byrow=T) 

Comb[1:3,]           #mostramos las tres primeras filas de la matriz Comb


# Ayuda a entender lo que hacemos después

Comb[2,pdT6[2,]]

Comb[2, -pdT6[2,]]


########################### p-valor exacto ######################


Theta<-array(0,924)

for (i in 1:924) {
Theta[i]<-mean(Comb[i,pdT6[i,]])-mean(Comb[i, -pdT6[i,]])}

Theta.obs<-mean(Treat)-mean(Cont)
pval<-sum(Theta >=Theta.obs)/choose(12,6)   # Cálculo exacto del p-valor

> pval
[1] 0.03138528

########################## Utilizando la librería exactRankTests ########

######## Instalar la librería exactRankTests y cargarla en R

library(exactRankTests)

perm.test(Ratbp$Treat, Ratbp$Cont, alternative="greater")


  2-sample Permutation Test (scores mapped into 1:(m+n) using rounded
        scores)

data:  Ratbp$Treat and Ratbp$Cont
T = 52, p-value = 0.03355
alternative hypothesis: true mu is greater than 0


## Recordamos que esta librería ya NO conviene usarla.



########################## Utilizando la librería coin ########

######## Instalar la librería coin y cargarla en R


library(coin)
GR<-as.factor(c(rep("Treat",6), rep("Cont",6)))

oneway_test(Blood~GR, distribution="exact",alternative="less")





############################## Aproximación al p-valor exacto (muestras sin reemplazamiento) ##########

set.seed(13)
a0<-mean(Ratbp$Treat)-mean(Ratbp$Cont)

boot.blood<-sample(1:924,size=499, replace=F)

B<-499
a<-array(0,B)

for (i in 1:length(boot.blood)){
j<-boot.blood[i]
media1<-mean(Comb[j,pdT6[j,]])
media2<-mean(Comb[j,-pdT6[j,]])
a[i]<-media1-media2}
pval.aprox<-(sum(a >=a0)+1)/(B+1) ##pval.boot

#plot(boot.blood)

pval.aprox

 0.032


############################## Test permutaciones (aproximación con reemplazamiento)  utilizando boot ##########

library(boot)

set.seed(13)
B<-499

blood.fun<- function(datos,i){
d<-datos[i]
MD<-mean(d[1:6])-mean(d[7:12])
MD}

boot.blood.b<-boot(Blood,blood.fun, R=B, sim="permutation")

plot(boot.blood.b)

pval.boot.b<-(sum(boot.blood.b$t >= boot.blood.b$t0)+1)/(B+1)

pval.boot.b

 0.03









############################## Test bootstrap (aproximación con reemplazamiento) ##########
 


set.seed(13)
a0<-mean(Ratbp$Treat)-mean(Ratbp$Cont)

boot.blood<-sample(1:924,size=499, replace=T)

B<-499
a<-array(0,B)
for (i in 1:length(boot.blood)){
j<-boot.blood[i]
media1<-mean(Comb[j,pdT6[j,]])
media2<-mean(Comb[j,-pdT6[j,]])
a[i]<-media1-media2}
pval.boot<-(sum(a >=a0)+1)/(B+1) ##pval.boot

#plot(boot.blood)
pval.boot

 0.04



##################  Cálculo del p-valor bootstrap (Librería boot) ##################


library(boot)

set.seed(40)
B<-499

blood.fun<- function(datos,i){
d<-datos[i]
MD<-mean(d[1:6])-mean(d[7:12])
MD}

boot.blood.b<-boot(Blood,blood.fun, R=B, sim="ordinary")

plot(boot.blood.b)

pval.boot.b<-(sum(boot.blood.b$t >= boot.blood.b$t0)+1)/(B+1)

pval.boot.b

 0.032





===================================================




library(boot)

set.seed(13)
B<-499

blood.fun<- function(datos,i){
d<-datos[i]
MD<-mean(d[1:6])-mean(d[7:12])
MD}

boot.blood.b<-boot(Blood,blood.fun, R=B, sim="ordinary")

plot(boot.blood.b)

pval.boot.b<-(sum(boot.blood.b$t >= boot.blood.b$t0)+1)/(B+1)

pval.boot.b

 0.038