##Permutaciones de una prueba de t, prueba de 1 y 2 colas con la funcion lm
##Ejercicio resuelto por Laura Paulina Osorio, modificaciones rems

R

#Cargando paquetes 
library(ISwR)
library(lmPerm)
data(energy)


# Ajustando modelo
mas1<-lm(expend~stature, data=energy)
#H0 no existe diferencia entre las medias de los grupos

# Obteniendo el valor de t a partir del modelo mas1
g <- summary(mas1)
f0 <- g[[4]][[6]]
##alternativamente 
#names(g)
#f0 <- g$coefficients[2,3]

# Señalando el número de permutaciones
ns <- 10000
Tstats <- array(0,ns)
# Realizando las permutaciones
for (i in 1:ns){
    m <- lm(expend ~ sample(stature), data= energy)
    g1 <- summary (m)
    Tstats [i] <- g1[[4]][[6]]
}
# Realizando el histograma de las permutaciones 
hist(Tstats, main="Distribución de los valores de t", xlim=c(-4,8))
abline (v=f0, col="red")#Agregando el valor de t obtenido en nuestro analisis

#Agregando el valor de t obtenido en nuestro analisis
q1<-quantile(Tstats,probs=c(0.95)) 
#Region de rechazo prueba 1 cola
q1

abline(v=q1, col="blue")

q2<-quantile(Tstats, probs=c(0.025,0.975)) #Region de rechazo prueba 2 colas
q2

abline(v=q2, col="black", lty=2)

colores <- c("red","blue","black")
legend("topright", c("t estimada","1 cola","2 colas"),lty=c(1,1,2), col= colores)


#Prueba de una cola,
## Ojo aqui empleamos la opcion  "menor" en vez de mayor porque R usa al nivel "lean" como referencia, observen los valores negativos del valor de t y de la diferencia entre las medias de los grupos
## Es la misma hipotesis: delgadas < obesas que obesas > delgadas, la otra opcion seria cambiar el nivel de referencia con la funcion relevel

mt <-t.test(expend~stature, data=energy, var.equal=TRUE, alternative="less") 
#Valor de probabilidad para prueba de una cola obtenida a partir de las permutacionees
options(scipen=10)
length(Tstats[Tstats > f0])/ns


#Valor absoluto de t, prueba de 2 colas, las plantulas de los grupos luz y sombra difieren en su longitud
bTstats<-abs(Tstats)

#Valor de probabilidad
length(bTstats[bTstats > f0])/ns

#comparandolo con el valor de t
names(mt)
mt$p.value

#Utilizando del paquete lmPerm, la funcion lmp que te proporciona los valores de probabilidad obtenidas por permutaciones
tlpem<-lmp(expend~stature, data=energy, perm="prob")
summary(tlpem)