Análisis de varianza: ANOVA. Efectos fijos y aleatorios en R.

En esta entrada evaluaremos el análisis de varianza con efectos fijos y efectos aleatorios, así como aquellas situaciones donde tenemos medidas con repeticiones. Como resumen teórico:
Resumen Anova Efectos Fijos y Aleatorios

Nuestros objetivos son:

• Deducir los valores esperados de las sumas de cuadrados medios (MS) entre tratamientos y dentro de cada tratamiento para el modelo de efectos fijos y efectos aleatorios.
• Evaluar cuál de los cuatro métodos estudiados para realizar comparaciones múltiples entre niveles medios de tratamientos por pares posee mejores propiedades.
• ¿Cómo afecta la violación de la hipótesis de normalidad a la prueba F en los modelos de efectos fijos y aleatorios?. ¿Influye de forma distinta en diseños con tamaños muestrales fijos y variables?.
• Elaborar un resumen sobre los métodos no paramétricos usuales en el análisis de varianza.

Aquí les dejo un resumen teórico del tema completo, y un ejemplo aplicado en R
anova_1_via

• Ejemplo: Análisis de datosen R (r-project)
  

########### ejemplo: ANOVA ##############

##Data
A_15=c(7,7,15,11,9)
B_20=c(12,17,12,18,18)
C_25=c(14,18,18,19,19)
D_30=c(19,25,22,19,23)
E_35=c(7,10,11,15,11)
##Examen gráfico
scores=data.frame(A_15,B_20,C_25,D_30,E_35)
boxplot(scores)
library(PASWR)
scores2=stack(scores) #preparación de los datos 
X<-scores2[,1]
INDEX<-scores2[,2]
oneway.plots(X,INDEX) #dotplot, boxplot y
design plot (means)

## Modelo con efectos fijos (FIXED MODEL)
# Las medias de los tratamientos son o no iguales: Ho: mu1=mu2=...=mua vs Ha: mui!=muj
# Cuando la Ho es cierta, se puede evaluar una afirmación equivalente en término de los efectos
de los tratamientos: Ho: tau1=tau2=...=taua vs Ha: taui!=0
##Estimaciones: E(MSerror)=sigma^2 y E(MStrat)=sigma^2 + sum((ni*taui^2)/(a-1))
#Prueba:
TreatMean<-tapply(X,INDEX,FUN=mean)
a<-length(TreatmentMean)
N<-length(X)
dft<-a-1
dfe<-N-a
GrandMean<-mean(X)
ni<-nrow(scores)
SStreat<-ni*sum((TreatMean-GrandMean)^2)
SStot<-sum((X-GrandMean)^2)
SSerror<-SStot-SStreat
MStreat<-SStreat/dft
MSerror<-SSerror/dfe
Fobs<-MStreat/MSerror
pvalue<-1-pf(Fobs,dft,dfe)
##equivalente
summary(aov(X~INDEX))
model.tables(aov(X~INDEX),type="means")
##Chequeo de supuestos; 3 supuestos en el componente de ERROR (se util los residuales como su estimador): independencia, distribución normal y varianza constante
mod.aov<-aov(X~INDEX)
library(MASS)
r<-stdres(mod.aov)
n<-length(X)
#Chequeo de independencia de errores
par(pty="s")
plot(1:n,r,ylab="Standardized
Residual",xlab="Ordered Value")
#Chequeo de normalidad: qqplot y shapiro.test()
par(pty="s")
qqnorm(r)
abline(a=0,b=1)
shapiro.test(r)
#Chequeo de varianza constante: plot de residuos (estandarizdos) vs valores ajustados; Levene
tm<-fitted(mod.aov)
plot(tm,r,xlab="Fitted Value",ylab="Standardized Residual")
med<-tapply(X,INDEX,median)
ZIJ<-abs(X-med[INDEX])
summary(aov(ZIJ~INDEX))
library(PASWR)
checking.plots(mod.aov) ##package PASWR
##comparación múltiple de medias
#Ho: Ho1 intersección Ho2 intersección ...HoK
library(multcomp)
library(multcompView)
CI<-TukeyHSD(aov(X~INDEX,which="INDEX"))
plot(CI,las=1)
INDEX.aov<-aov(X~INDEX)
MSE<-summary(aov(INDEX.aov))[[1]][2,3]
alpha.c<-0.05
ybari<-TreatmentMean
TcritLSD<-qt(1-alpha.c/2,dfe)
nn<-rep(ni,a)
LSD<-TcritLSD*sqrt(MSE)*sqrt(sum(1/nn))
TcritTUK<-qtukey(1-alpha.c/2,a,dfe)/sqrt(2)
HSD<-TcritTUK*sqrt(MSE)*sqrt(sum(1/nn)) #nn es un vector de ni y nj, con el length=número de tratamientos
library(gregmisc)
NS<-tapply(X,INDEX,length)
SE<-sqrt(MSE)/sqrt(NS)
t.v<-qt(.95,dfe)
ci.l<-ybari-t.v*SE
ci.u<-ybari+t.v*SE
barplot2(ybari,plot.ci=T,ci.l=ci.l,ci.u=ci.u,col="sky blue",ci.lwd=2)
title(main="Mean X por INDEX \n con CI individual 95%")
#multcompBoxplot(X~INDEX) Su gráfico no es fácilmente interpretable

## Modelo con efectos aleatorios (RANDOM MODEL)
##supuestos: eijNID(0,sigma), taui~NID(0,sigma), taui y eij son independientes
#Ho: sigma-subtau^2=0 vs Ha: sigmasubtau^2>0
#Estimaciones:
#cuando los a tratamientos tienen igual tam de muestreo: sig2=estim(sigma)^2=MSerror y sig2tau=estim(sigma)-subtau^2=(MStreat- MSerror)/n
#cuando los tamaños muestrales son desiguales, n se reemplaza por n`=1/(a-1)*sum(ni)-(sum(ni^2)/sum(ni))
summary(aov(X~INDEX))
MSC<-summary(aov(X~INDEX))[[1]][1,3]
MSE<-summary(aov(X~INDEX))[[1]][2,3]
#Estimación de los componentes de varianza
sig2tau<-(MSC-MSE)/n #nº de tratamientos