sábado, 12 de septiembre de 2009

Análisis de varianza de una sola vía

###   Un ejemplo ######

##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)

 ## 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
##Camino corto: equivalente al siguiente camino largo
summary(aov(X~INDEX))
model.tables(aov(X~INDEX),type="means")

##Camino largo
#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(TreatMean)
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)

##Chequeo de supuestos; 3 supuestos en el componente de ERROR (se utilizan 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: cuando la Ho del ANOVA es rechazada, analizamos cuáles son los tratamientos que producen  efectos diferentes en la variable respuesta
#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<-TreatMean
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="skyblue",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

 ##RANDOM MODEL
##supuestos: eijNID(0,sigma), taui~NID(0,sigma), taui y eij son independientes
#Ho: sigma-subtau^2=0 vs Ha: sigma-subtau^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