################################################################################################################################ ## Modelos lineales generales y disenho matricial. Casos particulares: analisis de varianza, covarianza. ## ## Funciones tapply, anova, model.matrix, expand.grid, glht y confint ## ## Clase 3 2011/10/11 Roberto Munguia-Steyer ## ################################################################################################################################ # Iniciando R R # Cargando paquetes library(effects) library(gplots) library(multcomp) # Prueba de t #generando valores, dimorfismo sexual en un aracnido set.seed(1234) tm <- 8 th <- 6 de <- 2 n <- 60 psex <- 0.5 sexo <- rbinom(n,1,0.5) table(sexo) nh<- table(sexo)[1] nm<- table(sexo)[2] sexo <-sort(sexo) sexo.cat <- ifelse(sexo==1,"m", "h") sexo.cat <- as.factor(sexo.cat) tam.h <- rnorm(nh,th,de) tam.m <- rnorm(nm,tm,de) tam <- ifelse(sexo.cat=="m",rnorm(nm,tm, de), rnorm(nh,th,de)) par(las=1) boxplot(tam~sexo.cat, names=c("Hembras", "Machos"), ylab= "Longitud del escudo dorsal (mm)") #Obteniendo las medias para cada grupo tapply(tam,sexo,mean) ## Prueba de t pdt <- t.test(tam~sexo, var.equal=T) names(pdt) pdt # Modelo lineal m1 <- lm(tam~sexo) coef(m1) m1b <- lm(tam~sexo.cat) coef(m1b) #Validacion del modelo par(mfrow=c(2,2)) plot(m1b) # Parametrizacion de efectos, coeficientes de la regresion summary(m1) #Analisis de varianza anova(m1) #Valores esperados para el tamanho de machos y hembras # Hembras y_hat_h = b0 tpred.h <- coef(m1)[1] # Machos y_hat_m = b0 + b1 tpred.m <- coef(m1)[1] + coef(m1)[2] #Cotejando los varlores predichos de la prueba de t y del modelo lineal pdt$estimate tpred.h;tpred.m #Funcion model.matrix model.matrix(m1) #Parametrizacion alterna: de los efectos m1c <- lm(tam~sexo.cat-1) coefi<- coef(m1c) coefi model.matrix(m1c) #Intervalos de confianza para los valores predichos para la longitud de cada grupo ic <- confint(m1c) ic coefiv <- as.vector(coefi) ici <- as.vector(ic[,1]) ics <-as.vector(ic[,2]) #Graficando los valores predichos y los intervalos de confianza barplot2(coefi, ci.l = ici, ci.u =ics, plot.ci=TRUE, ylim= c(0,10),names=c("Hembras", "Machos"),ylab= "Longitud del escudo dorsal (mm)", col="white") #Usando el paquete effects efi<- allEffects(m1b) efi #Viendo el objeto efi como lista str(efi) #Obteniendo los valores predichos y los intervalos de confianza summary(efi) ## Analisis de varianza data(cholesterol) #help(cholesterol) str(cholesterol) head(cholesterol) names(cholesterol) <- c("ndosis","colesterol") levels(cholesterol$ndosis)<- c("1d", "2d", "3d", "c1", "c2") str(cholesterol) #Colocando el nivel control 1 como referencia cholesterol$ndosis <-relevel(cholesterol$ndosis, ref="c1") str(cholesterol) # Grafico de cajas par(las=1) boxplot(colesterol~ndosis, data= cholesterol) #Analisis m2 <- lm(colesterol~ndosis,data=cholesterol) #Validacion par(mfrow=c(2,2)) plot(m2) resi <-resid(m2) boxplot(resi~ndosis,data=cholesterol) #Menor dispersion en los residuos del segundo grupo, verificar en el topico de gls #Analisis de varianza anova(m2) #Coeficientes del modelo #Los grupos difieren del primer control? summary(m2) #Comparaciones multiples, prueba de Tukey cht <- glht(m2, linfct = mcp(ndosis = "Tukey")) cht #Intervalos de confianza confint(cht) #Grafica plot(confint(cht)) summary(cht, test = univariate()) #Correccion de bonferroni summary(cht, test = adjusted("bonferroni")) #Valores predichos efi2 <- allEffects(m2) summary(efi2) #Otra posibilidad, utilizar las funciones expand.grid para generar las combinaciones de interes y la funcion predict a<-expand.grid(group=levels(cholesterol$ndosis)) ei<-data.frame(a,predict(m2,se=T)) #Retirando el intercepto, valores predichos medios para cada nivel del factor m2b <- lm(colesterol~ndosis-1,data=cholesterol) coef(m2b) confint(m2b) #Disenho de la matriz model.matrix(m2) # Ejercicio niveles de coagulacion en la sangre (base de datos coagulation en el paquete faraway) #1. Carga y revisa la base de datos coagulation #2. Inspecciona visualmente la variable predictora #3. Ajusta el modelo y realiza un analisis de varianza #4. Valida y en su caso especifica cuales son los supuestos que potencialmente estan siendo violados, discute si el tomarlos en cuenta pudieran cambiar las conclusiones #5. Cuales son los valores predichos y los intervalos de confianza al 95 % para cada nivel del la variable predictora categorica, usa el paquete effects y la funcion predict #6. Realiza comparaciones multiples entre los niveles del factor usando, graficando las diferencias entre las medias predichas para cada nivel y sus intervalos de confianza al 95% # asi como ajustando al valor de probabilidad de los contrastes con el ajuste de Holm library(faraway) data(coagulation) help(coagulation) boxplot(coag~diet,data=coagulation) m3 <-lm(coag~diet,data=coagulation) par(mfrow=c(2,2)) plot(m3) anova(m3) summary(m3) options(scipen=10) summary(m3) #a m3b <- lm(coag~diet-1, data=coagulation) coef(m3b) #b b <- expand.grid(group=levels(coagulation$diet)) ei2 <- data.frame(b,predict(m3,se=T)) #c efi3 <-allEffects(m3) summary(efi3) #Comparaciones multiples, prueba de Tukey cht2 <- glht(m3, linfct = mcp(diet = "Tukey")) cht2 #Grafica de los intervalos de confianza plot(confint(cht2)) summary(cht2, test = adjusted("holm"))