# Anova unifactorial y multifactorial

# Librerías "lsr" ("Learning Statistics with R"), "effects", "car" (para el test de Levene) y Deducer
# install.packages("lsr")
# install.packages("effects")
# install.packages("Deducer")
library(lsr)
library(effects)
library(car)
library(Deducer)

# Matriz de datos AQU.rda
load(url("http://pagines.uab.cat/plopez/sites/pagines.uab.cat.plopez/files/AQU.rda"))

# Asignación de nombres a las variables
Y=AQU$ICO
A=AQU$Area
B=AQU$Sexo

############ Anova unifactorial
############ ICO - Area

# Tabla de medias
summary(Y ~ A)
descriptive.table(vars = d(Y), strata = d(A),data= Valora, 
                  func.names =c("Mean","St. Deviation","Valid N"))
library(psych)
describe.by(Y, A)

# Gráficos: de medias y diagramas de caja
ggplot() + geom_line(aes(x=Area, y=ICO, group=1), data=na.omit(AQU),fun.data=mean_sdl,
                     fun.args=list(mult=1), stat='summary')
plot(Y ~ A)

# Análisis de varianza
AV=aov(Y ~ A)
summary(AV)

# Condicionbes de aplicación
hist(AV$residuals) # Histograma de los residuos
plot(AV) # Gráficos de ajuste del modelo
leveneTest(Y ~ A, center=mean) # Prueba de igualdad de varianzas: Levene Test

# Comparación de medias: Tukey's Honest Significance Test
summary(effect(term="A",mod=AV))
plot(effect(term="A",mod=AV,default.levels=20),rescale.axis=FALSE)
TukeyHSD(AV)
plot(TukeyHSD(AV),las=3)

# Capacidad explicativa y coeficientes
etaSquared(AV)
summary.lm(AV)

############ Anova unifactorial
############ ICO - Sexo
summary(Y ~ B)
ggplot() + geom_line(aes(x=Sexo, y=ICO, group=1), data=na.omit(AQU),fun.data=mean_sdl,
                     fun.args=list(mult=1), stat='summary')
plot(Y ~ B)
AV=aov(Y ~ B)
summary(AV)
hist(AV$residuals)
plot(AV)
leveneTest(Y ~ A, center=mean)
summary(effect(term="B",mod=AV))
plot(effect(term="B",mod=AV,default.levels=20),rescale.axis=FALSE)
TukeyHSD(AV)
plot(TukeyHSD(AV),las=3)
etaSquared(AV)
summary.lm(AV)

############ Anova multifactorial
############ ICO - Area + Sexo 

# Tabla de medias
aggregate(Y ~ A + B, FUN="mean")

# Gráficos
ggplot() + geom_line(aes(x=Area,y=ICO,group=Sexo),data=na.omit(AQU),fun.data=mean_sdl,fun.args=list(mult = 1),stat = 'summary')
         + geom_line(aes(x=Area,y=ICO,group=Sexo),data=na.omit(AQU),fun.data=mean_sdl,fun.args=list(mult = 1),stat = 'summary')
plot(Y ~ A + B, col=rainbow(6))
plot(A:B,Y, col=rainbow(2))
boxplot(Y ~ A + B, col=rainbow(6))
interaction.plot(A,B,Y)

# Anova
AV=aov(Y ~ A + B + A:B, contrasts=list(A=contr.sum, B=contr.sum))
plot(AV) # Gráficos de ajuste del modelo
leveneTest(Y ~ A*B, center=mean) # Test de Levene
summary(AV) # SC Tipo I
Anova(AV, type="III") # SC Tipo III

# Comparación de medias: Tukey's Honest Significance Test
summary(effect(term="A:B",mod=AV))
plot(effect(term="A:B",mod=AV,default.levels=20),rescale.axis=FALSE)
TukeyHSD(AV)
plot(TukeyHSD(AV),las=1)

# Capacidad explicativa y coeficientes
etaSquared(AV, type=3)
summary.lm(AV)
