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



## Introducimos los datos en R

nuevom<-c(37,49,55,57)
mtradicional<-c(23,31,46)
test<-c(nuevom,mtradicional)



# Después hay que "instalar" la librería BSDA y  cargarla en R


library(BSDA)


pdT6<- SRS(1:7,4)
dim(pdT6)       # 35 combinaciones posibles de 7 elementos tomados de 4 en 4

pdT6[1:5,]       # mostramos las cinco primeras filas


# Alternativamente

pdT6=t(Combinations(7,4))

Comb<-matrix(rep(c(nuevom,mtradicional),35),ncol=7,byrow=T)
Comb[1:3,]           #mostramos las tres primeras filas

# Ayuda a entender lo que hacemos después

Comb[2,pdT6[2,]]

Comb[2, -pdT6[2,]]



Theta<-array(0,35)

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

Theta.obs<-mean(nuevom)-mean(mtradicional)
pval<-sum(Theta >=Theta.obs)/choose(7,4)   # Cálculo exacto del p-valor

pval

[1] 0.05714286





########## TEST DE PERMUTACIONES CON H1 GENERAL###############


### Si quisiéramos comprobar simplemente si ambos métodos difieren,
deberíamos calcular el p-valor

### para un contraste cuya alternativa es la más general posible F
distinto de G






Theta<-array(0,35)

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

Theta.obs<-mean(nuevom)-mean(mtradicional)
pval<-sum(abs(Theta) >=abs(Theta.obs))/choose(7,4)   # Cálculo exacto del p-valor

> pval
[1] 0.1142857

# Observemos que este p-valor es 0.05714286 * 2









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

## Sin utilizar la libreria boot



set.seed(40)
estadistico.obs<-mean(nuevom)-mean(mtradicional)

boot.test<-sample(1:35,size=500, replace=T)

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

qqnorm(estadistico)


> pval.boot
[1] 0.04391218