############################################################################################# ### Modelos lineales: Heterogeneidad de varianza, minimos cuadrados generalizados ### ### Rems 2011/10/27 y 2011/11/03 ### ############################################################################################# R #Cargando paquetes library(MASS) library(nlme) library(bbmle) library(car) #de ser necesario instalar los paquetes sandwich y lmtest library(sandwich) library(lmtest) #Heterocedasticidad, cuando la variable predictora es continua, #Incremento en el numero de enfermedades registradas por municipio al anho # x numero de habitantes por 100 set.seed(1010) n <- 80 a <- 20 b <- 0.2 x <- runif(n, 1,100) y <- a + b*x + 0.3*x*rnorm(n, sd=0.5) #Ajuste del modelo plot(y~x) m1 <- lm(y~x) abline(m1, col='blue') summary(m1) coef(m1) #Validacion par(mfrow=c(2,2)) plot(m1) ###Posibles alternativas para lidiar con la heterocedasticidad, sugerencia de transformacion spreadLevelPlot(m1) ## a) Transformacion de la variable de respuesta yt <- y^(-1.8) m1a <- lm(yt ~ x) par(mfrow=c(2,2)) plot(m1a) #Heterogeneidad de varianza, aunque menor, grafica de los quantiles de los residuos asimetricos #evitemos esta sugerencia por la perdida de interpretabilidad de los coefs. ## b) Minimos cuadrados ponderados m2 <-update(m1, weights=1/x) AICtab(m1,m2, base=TRUE, weights=TRUE) #El modelo con mayor soporte es el de minimos cuadrados ponderados par(mfrow=c(1,2)) plot(fitted(m1), rstandard(m1)) plot(fitted(m2), rstandard(m2)) #Mejoro aunque aun se percibe heterogeneidad ## c) Correccion a la heterocedasticidad con el metodo de White, ver libro de John Fox (2002). coeftest(m1, vcov = vcovHC(m1, type = "HC3")) #compararlo con el modelo lineal convencional summary(m1) ## d) Minimos cuadrados generalizados, primero ajustando el equivalente al modelo lineales m1gls <- gls(y~x) # 1. Aumento de la varianza proporcional a la variable predictora m2gls <- update(m1gls, weights=varFixed(~ x)) # Varianza elevada a la potencia, no usar cuando x puede tener un valor de 0 m3gls <- update(m1gls, weights=varPower(form=~ x)) # incremento de la varianza exponencial m4gls <- update(m1gls, weights=varExp(form=~ x)) AICtab(m1gls,m2gls,m3gls,m4gls, base=T, weights=T) #Coeficientes del modelo con mayor soporte, comparandolo con el modelo de varianza homogenea summary(m3gls) summary(m1gls) #Los coeficientes estimados del modelo de minimos cuadrados presentan menor sesgo con relacion a los parametros poblacionales. #Validacion del modelo con mayor soporte rst <- resid(m3gls, type="normalized") ajus <- fitted(m3gls) plot(x=ajus, y=rst, xlab = "Valores predichos", ylab = "Residuos Estandarizados") cgls <- coef(m3gls) clm <- coef(m1) par(las=1) plot(x,y) abline(clm[1],clm[2], col="blue", lty=2) abline(cgls[1],cgls[2], col="red") legend("topleft", c("gls", "lm"), col= c("red","blue"), lty=c(1,2)) ## EJERCICIO, base de datos (Mandible) en paquete MASS. ## tamanho de la mandibula en fetos de diferente edad #1. Ajusta y el valida modelo lineal #2. Comparalo con un modelo de cuadrados minimos ponderados #3. Utiliza los modelos de minimos cuadrados generalizados con diferentes estructuras de la varianza residual y encuentra la que presenta mayor soporte #4. Valida el modelo con mejor soporte, obten los coeficientes, las pruebas estadisticas y los respectivos valores de probabilidad #5. Grafica los valores predichos tanto los obtenidos del modelo lineal como del modelo de minimos cuadrados data(Mandible) #help(Mandible) bd <- Mandible str(bd) names(bd) <- c("edad","longi") plot(bd$edad, bd$longi) m5 <- lm(longi~edad, data=bd) par(mfrow=c(2,2)) plot(m5) m6 <- lm(longi~edad, weights=1/edad, data=bd) AICtab(m5,m6) summary(m5);summary(m6) par(mfrow=c(1,2)) plot(fitted(m5), rstandard(m5)) plot(fitted(m6), rstandard(m6)) #Mejoro aunque aun se percibe heterogeneidad m5gls <- gls(longi~edad, data=bd) ev1 <- varFixed(~ edad) ev2 <- varPower(form =~ edad) ev3 <- varExp(form =~ edad) m6gls <- gls(longi~edad, weights=ev1, data=bd) m7gls <- update(m6gls, weights=ev2) m8gls <- update(m7gls, weights=ev3) AICtab(m5gls,m6gls,m7gls,m8gls, base=T, weights=T) summary(m8gls) #Validacion rst <- resid(m8gls, type="normalized") ajus <- fitted(m8gls) plot(x=ajus, y=rst, xlab = "Valores predichos", ylab = "Residuos Estandarizados") cgls <- coef(m8gls) clm <- coef(m5) #Valores predichos par(las=1) plot(bd$edad,bd$longi) abline(clm[1],clm[2], col="blue", lty=2) abline(cgls[1],cgls[2], col="red") legend("topleft", c("gls", "lm"), col= c("red","blue"), lty=c(1,2)) ##Ejercicio, sobre longitud del lobulo frontal cangrejos en relacion al sexo y a a la longitud del carapacho ## data (crabs) en el paquete MASS, por fines didacticos omitimos que los datos preceden 2 especies diferentes data(crabs) str(crabs) cb <-crabs #Graficando las observaciones por el sexo al que pertenecen plot(cb$CL, cb$FL, type="n") points(cb$CL[cb$sex=="F"], cb$FL[cb$sex=="F"], pch=19) points(cb$CL[cb$sex=="M"], cb$FL[cb$sex=="M"], pch=21) #Inspeccion visual boxplot(CL~sex, data=cb) #Ajuste del modelo m9 <- lm(FL ~ sex*CL, data=cb) #Validacion del modelo con una grafica condicional de los residuos por sexo e <- resid(m9) coplot(e~CL|sex, data=cb) #Grafico de cajas de los residuos por sexo boxplot(e~sex, data=cb) #Diferentes estructuras para modelar la heterogeneidad de varianza ev4 <- varFixed(~ CL) ev5 <- varPower(form =~ CL) ev6 <- varExp(form =~ CL) ev7 <- varIdent(form=~ 1|sex) #Ajuste de los modelos m9 <- gls(FL ~ sex*CL, data=cb) m10 <- update(m9, weights=ev4) m11 <- update(m9, weights=ev5) m12 <- update(m9, weights=ev6) m13 <- update(m9, weights=ev7) #Seleccion de modelos AICtab(m9, m10,m11,m12,m13, base=TRUE, weights=TRUE) ##Solo con fines demonstrativos utilizar una funcion que combina un par de estructura de la varianza #Aumento de la varianza de forma exponencial diferenciada por sexo ev8 <- varComb(varIdent(form =~1|sex), varExp(form=~CL)) m14 <- update(m9, weights=ev8) AICtab(m9, m10,m11,m12,m13,m14, base=TRUE, weights=TRUE) m11 #Comparando el modelo global con interaccion y el modelo aditivo, utilizar maxima verosimilitud m11a <- update(m11, method="ML") m11a m12 <- gls(FL ~ sex+CL, weights= ev5, method="ML", data=cb) #Comparando los modelos anova(m11a,m12) #Validacion del modelo e2 <-resid(m11, type="normalized") coplot(e2~CL|sex, data=cb)