#
#    A continuación, un ejemplo de análisis MANOVA como el
#    solicitado en la Tarea 4.
#
#
aritm <- read.table("aritm.dat",header=TRUE)
#
#  Veamos que pinta tienen los datos
#
aritm[1:3,]
summary(aritm)
#
#   Hmmm, "Metodo" ha sido leida como una variable numérica,
#   y, siendo como es una variable cualitativa, nos
#   interesará convertirla en un factor.
#
aritm[,"Metodo"] <- as.factor(aritm[,"Metodo"])
aritm[1:3,]
#
#   Aparentemente nada ha cambiado...
#
summary(aritm)
#
#   ...pero "summary" muestra a las claras la
#   diferencia.
#
boxplot(NotaSumas ~ Metodo, data=aritm)
boxplot(NotaMult  ~ Metodo, data=aritm)
#
#   Representaciòn de los datos sin desagregar
#
pairs(NotaSumas ~ NotaMult, data=aritm)
#
#   Representación de los datos desagregados
#
coplot(NotaSumas ~ NotaMult | Metodo, data=aritm)
#
#   Una alternativa más laboriosa:
#   Vamos a representar las observaciones, cosa
#   fácil dado que son bivariantes. Emplearemos
#   un carácter diferente para cada método.
#
plot(aritm[,"NotaSumas"],aritm[,"NotaMult"],
     type="n",xlab="NotaSumas",ylab="NotaMult")
title(main="Resultado de cuatro métodos de enseñanza \n de operaciones elementales",
      sub="Cada + representa el par de notas obtenidas por uno de los 80 niños en la muestra.",
      xlab="NotaSumas",ylab="NotaMult")
desglose <- split(aritm[,2:3],aritm[,"Metodo"])
for (i in 1:length(desglose)) {
    points(desglose[[i]][,1],desglose[[i]][,2],
           pch="+",col=i)
    Xbar <- apply(desglose[[i]],2,mean)
    text(Xbar[1],Xbar[2],
         eval(substitute(expression(bold(bar(X)[grupo])),
                         list(grupo=i))),
         cex=1.2,col=i)
}
#   Claramente, parece habe diferencias entre los
#   métodos, que un MANOVA difícilmente dejará de
#   poner de manifiesto.
#
manova.niños <- manova(cbind(NotaSumas,NotaMult) ~ Metodo,
                       data=aritm)
summary(manova.niños,test="Pillai")
summary(manova.niños,test="Wilks")
summary(manova.niños,test="Hot")
summary(manova.niños,test="Roy")
#
#   No obstante, los resultados están sujetos a crítica:
#   la hipotesis de normalidad no es crítica, y no parece
#   que sea clamorosamente incumplida. La hipótesis de
#   homogeneidad de matrices de covarianzas, en cambio,
#   si parece completamente fuera de lugar, como la simple
#   contemplación del gr-afico pone de manifiesto.
#
#   Estimemos las matrices de covarianzas...
#
mat.cov <- lapply(desglose,FUN=cov.wt)
mat.cov <- lapply(mat.cov,FUN=function(x) { x$cov })
#
#   ...y de correlaciones
#
lapply(mat.cov,FUN=cov2cor)
#
#   Podemos hacer un contraste formal de la hipótesis
#   de igualdad de las matrices de covarianzas: el determinante
#   de la matriz de covarianzas estimada sobre toda la
#   muestra elevado a -N/2 es...
#
S <- det( cov.wt(aritm[,2:3])$cov ) ^ (- 80 / 2)
#
#   Análogamente podemos calcular lo mismo para cada una de las
#   cuatro submuestras:
#
Si <- lapply(mat.cov,FUN=function(x) { det(x) ^ (- 20 / 2) })
#
#   El estadístico será,
#
    lambda <- -2 * log( S / prod(unlist(Si)))
#
#   Debemos compararlo con una chi^2 cuyos grados de libertad
#   sean la diferencia de parámetros entre numerador y denominador:
#
#   Grados de libertad: 5*3 - 3 = 12
#
lambda
1 - pchisq(lambda,12)
#
#   La probabilidad en la cola es prácticamente cero, y la
#   hipótesis es contundentemente rechazada.
#