#############################################################################################
### Alternativas para cuando los modelos no presentan una distribucion normal en los      ###
### residuos: Transformaciones a logaritmos, Box-Cox y tecnicas de muestreo               ###
### Roberto Munguia-Steyer                                                                ###
### Practica 5 Modelos lineales generales, generalizados y mixtos               20111024  ###
#############################################################################################

R


#Instalar el paquete (lmPerm)
library(lmPerm)
library(faraway)
library(ISwR)
library(lattice)

#Base de datos malaria
data(malaria)

#help(malaria)

str(malaria)

#Cambiando los nombres de las columnas
names(malaria) <- c("individuo","edad", "ac", "mal")
head(malaria)

#Podemos dejar malaria como una variacle indicadora o como un factor ya que solo tiene 2 niveles
malaria$mal2 <- factor(malaria$mal)
str(malaria)

#resumen de los datos por columnas
summary(malaria)

#Inspeccion visual de la distribucion de los datos
par(las=1)
boxplot(ac ~ mal, data=malaria, ylac= "Conteo de anticuerpos") #Heterogeneidad de varianza


par(las=1)
dotplot(ac~edad|mal, data=malaria) #Distribucion de los datos


#Ajuste del modelo
m1 <- lm(ac~mal*edad, data=malaria)

#Cuantos parametros
coef(m1) 

#Validacion
par(mfrow=c(2,2))
plot(m1) #Heterogeneidad de varianza y distribucion poco parecida a la normal

##Transformaciones mas usadas, logaritmos (ya vista), raiz cuadrada y box-cox, en caso de proporciones raiz cuadrada de arco coseno.

#Transformacion de concentracion de ac a raiz cuadrada

malaria$sqac <- sqrt(malaria$ac)

m2 <- lm(sqac~mal*edad, data=malaria)
par(mfrow=c(2,2))
plot(m2) #Mejor pero todavia la distribucion de los residuos no se parece a la d. normal

#Transformacion de box-cox a la concentracion de ac 
boxcox(m1)

#lambda=0, equivale a una transformacion a logaritmos

#Transformacion de concentracion de ac a logaritmos
malaria$lac <- log(malaria$ac)

m4 <- lm(lac~mal*edad, data=malaria)

par(mfrow=c(2,2))
plot(m4)

##Metodos de remuestreo, permutaciones, barajando las observaciones
#Crecimiento de plantulas en condiciones de luminosidad contrastantes, una semana, medicion de la longitud

#fijar la generacion de valores pseudoaleatorios para fines de contar con los mismos estadisticos y valores de p. en clase
set.seed(123)

#observaciones por categoria
n<-25
l<-rnorm(n,mean=8.5,sd=0.8)
s<-rnorm(n,mean=8.0,sd=0.8)

lumi<-c(l,s)
cat <-c(rep("l",n),rep("s",n))
plantula <- data.frame(lumi,cat)

tapply(plantula$lumi,plantula$cat,mean)


boxplot(lumi~cat, data=plantula,ylab="Longitud de la plantula (mm)")

pla <-t.test(lumi~cat, data=plantula, var.equal=TRUE)
pla

#viendo que podemos extraer del objeto pla
names(pla)

#valor de t
t0 <-pla$statistic



#Barajando las observaciones, simulacion y obtencion de los valores de t
ns <- 5000
tstats <- numeric(ns)
 
for (i in 1:ns){
 m <- t.test(lumi~sample(cat), data=plantula)
 tstats [i] <- m$statistic
 }

#Distribucion de los valores de t, barajando los datos, H0
hist(tstats, main="Distribucion de los valores de t")

#Agregando el valor de t obtenido en nuestro analisis
abline(v=t0, col="red")

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, e.g. las plantulas con luz crecen mas que las sombreadas
pla2 <-t.test(lumi~cat, data=plantula, var.equal=TRUE, alternative="greater")
pla2

#Valor de probabilidad para prueba de una cola obtenida a partir de las permutacionees
length(tstats[tstats > t0])/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 > t0])/ns

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

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


## Retornando al primer modelo
## Permutaciones Utilizando el paquete lpermTest 
m5 <- lmp(ac~mal*edad, data=malaria, center=FALSE)
summary(m1);summary(m5)

#Ojo las permutaciones pueden evadir el supuesto de normalidad, pero no al de la heterogeneidad de varianza.


#Ejercicio
#1. Usando la base de la energia calorica consumida por mujeres obesas y delgadas, que se encuentra en el la base de datos energy realiza:
#2. Inpeccion visual de las variables
#3. Ajuste del modelo utilizando la funcion t.test 
#4. Realizando  10000 permutaciones con la funcion sample, obtenlos valores de probabilidad considerando :
#4a que las mujeres obesas consumen mas energia, prueba de una cola
#4b que las mujeres de los dos niveles tienen gastos energeticos diferentes,prueba de 2 colas
#5  Coteja tus resultados con la funcion lmPerm, modula la notacion para obtener los valores de probabilidad con 4 digitos
#6  OPCIONAL, realiza las permutaciones con la funcion lm si te sientes mas ambicioso (+1 en la calificacion final al primero que lo entregue o mande
#   por correo)