############################################################################################# ### Modelos lineales: analisis de covarianza y de dos vias ### ### Clase 3 2011/10/12 ### ############################################################################################# #Iniciando R y cargando paquete R library(ISwR) library(gplots) library(effects) library(bbmle) ############################################################################################## #1. Analisis de covarianza #Accesando base de datos, tomados de Crawley "The R Book", pag. 505 bd <-read.table("contrastes.txt", header=T, sep= "\t") str(bd) levels(bd$sexo) <- c("hembra","macho") #Inspeccion visual par(mfrow=c(1,2)) boxplot(peso~sexo, data=bd) dotchart(bd$edad, groups=bd$sexo) #Ajuste del modelo m1 <-lm(peso~sexo*edad, data=bd) #Validacion del modelo par(mfrow=c(2,2)) plot(m1) #Analisis anova(m1) summary(m1) # que representan los coeficientes, intercepto y pendiente para cada grupo int_f <- coef(m1)[1] pend_f <- coef(m1)[3] int_m <- coef(m1)[1] + coef(m1)[2] pend_m <-coef(m1)[3] + coef(m1)[4] #Disenho matricial model.matrix(m1) #Graficando par(las=1) plot(bd$edad,bd$peso,type="n", xlab= "Edad", ylab="Peso") points(bd$edad[bd$sexo=="hembra"], bd$peso[bd$sexo=="hembra"], pch=19) #hembras circulos negros points(bd$edad[bd$sexo=="macho"], bd$peso[bd$sexo=="macho"], pch=21) #machos circulos blancos #Intercepto y pendiente hembras abline(int_f,pend_f, col="blue") abline(int_m,pend_m, col="black",lty=2) # 2. Analisis de varianza de dos vias: modelo aditivo e interactivo data(ToothGrowth) help(ToothGrowth) tg <- ToothGrowth str(tg) summary(tg) names(tg) <- c("longi", "supl", "dosis") tg$dosis <- as.factor(tg$dosis) #Inspeccion visual par(las=1) boxplot(longi ~ supl + dosis, data = tg, ylab= "Longitud de los odontoblastos") ##Modelo aditivo m2 <- lm(longi ~ supl+dosis,data=tg) par(mfrow=c(2,2)) plot(m2) anova(m2) summary(m2) #Obteniendo los valores estimados para cada combinacion de los dos tratamientos levels(tg$supl) levels(tg$dosis) #Como es el disenho de la matriz model.matrix(m2) coef(m2) jugd05 <- coef(m2)[1] jugd1 <- coef(m2)[1] + coef(m2)[3] jugd2 <- coef(m2)[1] + coef(m2)[4] aad05 <- coef(m2)[1] + coef(m2)[2] aad1 <- coef(m2)[1] + coef(m2)[2] + coef(m2)[3] aad2 <- coef(m2)[1] + coef(m2)[2] + coef(m2)[4] #Valores predichos para cada combinacion valest <- c(jugd05,jugd1,jugd2,aad05,aad1,aad2) names(valest) <- c("jugd05","jugd1","jugd2","aad05","aad1","aad2") valest #Generando una nueva base de datos para los valores predichos a<-expand.grid(dosis=levels(tg$dosis),supl=levels(tg$supl)) #Comparando los valores estimados "a mano" con los obtenidos con la funcion predict predict(m2,newdata=a);valest #Si queremos incluir los errores estandar en la nueva base de datos de los valores predichos a2 <-data.frame(a,predict(m2,newdata=a, se=T)) a2 ##Modelo interactivo, analis de varianza con interaccion #Obtencion de las medias y errores estandar para cada combinacion mpg<-tapply(tg$longi,list(tg$supl,tg$dosis),mean) mpg #No existe una funcion default para estimar el error estandar es <- function(x) sqrt(var(x)/length(x)) ees <- tapply(tg$longi, list(tg$supl,tg$dosis),es) #Promedios +- errores estandar bi <- mpg-ees bs <- mpg+ees #Graficando los valores promedios de las posibles combinaciones de los dos factores barplot2(mpg, beside=TRUE) colores <- rep(c("white","darkgray"),3) nombres <- c("Aa-0.5","Jugo-0.5","Aa-1","Jugo-1","Aa-2","Jugo-2") par(las=1) barplot2(mpg, plot.ci=TRUE, ci.l = bi, ci.u = bs, beside=TRUE, col=colores, ylim=c(0,30), names=nombres, cex.names=0.9, ylab = "Longitud de los dientes") #Anova de dos vias m3 <- lm(longi ~ supl*dosis,data=tg) #Validacion par(mfrow=c(2,2)) plot(m3) #Analisis de varianza anova(m3) #Soporte relativo con relacion al modelo aditivo AICtab(m2,m3, base=TRUE,weights=TRUE) #Coeficientes del modelo coef(m3) #Disenho matricial model.matrix(m3) #Valores predichos para cada grupo coef(m3) jugd05i <- coef(m3)[1] jugd1i <- coef(m3)[1] + coef(m3)[3] jugd2i <- coef(m3)[1] + coef(m3)[4] aad05i <- coef(m3)[1] + coef(m3)[2] aad1i <- coef(m3)[1] + coef(m3)[2] + coef(m3)[3]+ coef(m3)[5] aad2i <- coef(m3)[1] + coef(m3)[2] + coef(m3)[4] + coef(m3)[6] #Valores predichos para cada combinacion valesti <- c(jugd05i,jugd1i,jugd2i,aad05i,aad1i,aad2i) names(valesti) <- c("jugd05","jugd1","jugd2","aad05","aad1","aad2") valesti #Valores predichos con la función allEffects efdien<- allEffects(m3, xlevels=list(supl=levels(tg$supl), dosis=levels(tg$dosis))) efdien #Valores predichos con los intervalos de confianza al 95 % summary(efdien) #Ejercicio #chickwts {datasets} #Realiza un analisis de varianza con la base de datos chickwts del paquete datasets #Efectua un grafico de barras con los valores promedios de cada grupo y sus errores estandar #Obten los valores predichos para cada nivel a partir de los coeficientes y usando el paquete allEffects. data(chickwts) str(chickwts) m4 <- lm(weight~feed, data=chickwts) par(mfrow=c(2,2)) plot(m4)