# Capítulo 7. GLM para variables respuesta de tipo conteo # 7.1. GLM con distribuciones de errores Poisson library(ADER) data(ara) str(ara) # 7.1.1. Análisis exploratorio y ajuste del modelo library(lattice) xyplot(seedlings ~ dungs|property, ara, xlab="N. de bostas", ylab="N. de plántulas") library(lattice) xyplot(log(seedlings+1) ~ dungs|property, ara, xlab="N. de bostas", ylab="log(N. de plántulas)") glm.ara <- glm(seedlings ~ dungs*property, ara, family=poisson) # 7.1.2. Inferencia en GLM: resolviendo hipótesis sobre causali- # dad library(car) Anova(glm.ara, type="II") # 7.1.3. Interpretación del modelo # A) La función summary() summary(glm.ara) # B) Devianza explicada por el modelo 1- (deviance(glm.ara)/deviance(update(glm.ara, .~1))) library(modEvA) Dsquared(glm.ara) # C) Gráficos de predicciones library(visreg) visreg(glm.ara, "dungs", by = "property", overlay = T, scale = "response", ylab = "N. de plántulas", xlab="N. de bostas") xv <- 0:2500 yv.camp <- predict(glm.ara, list(dungs=xv, property = rep("Campesino", length(xv))), type="response", se.fit=T) yv.Ef <- predict(glm.ara, list(dungs=xv, property = rep("Empresa Forestal", length(xv))), type="response", se.fit=TRUE) library(scales) plot(seedlings~dungs, ara, pch=16, col=c("steelblue", "firebrick"), ylab="N. de plántulas", xlab="N. de bostas", ylim=c(0,25)) polygon(x=c(xv, rev(xv)), y=c(yv.camp$fit + 1.96*yv.camp$se.fit, rev(yv.camp$fit - 1.96*yv.camp$se.fit)), col=alpha("firebrick", .6), border=NA) lines(xv, yv.camp$fit, lwd=2, col="firebrick") polygon(x=c(xv, rev(xv)), y=c(yv.Ef$fit + 1.96*yv.Ef$se.fit, rev(yv.Ef$fit - 1.96*yv.Ef$se.fit)), col=alpha("steelblue", .6), border=NA) lines(xv, yv.Ef$fit, lwd=2, col="steelblue") legend("topright", lwd=c(2,2), col=c("firebrick", "steelblue"), legend=c("Pequeño propietario", "Gran empresa forestal"), bty="n") # 7.1.4. Revisión de los supuestos del modelo e idoneidad de la # familia de distribución de errores glm.ara.Q <- glm(seedlings ~ dungs*property, ara, family=quasipoisson) summary(glm.ara.Q)$dispersion library(car) Anova(glm.ara.Q, type="II") # B) ¿Cómo podemos reconocer la sobredispersión y la subdispersión? library(aods3) gof(glm.ara) # Cuadro 7.1. (continuación) obs <- ara$seedlings fit <- glm.ara$fitted obs - fit residuals(glm.ara, type="response") (obs - fit)/fit residuals(glm.ara, type="working") (obs - fit)/sqrt(fit) residuals(glm.ara, type="pearson") poisson.dev <- function (y, mu) { sqrt(2 * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))) * ifelse(y > mu, 1, -1) } poisson.dev(y = obs, mu = fit) residuals(glm.ara) glm.ara$residuals # 7.2. GLM con distribución de errores binomial nega- # tiva. library(MASS) glm.ara.nb <- glm.nb(seedlings ~ dungs*property, ara) summary(glm.ara.nb)$theta library(car) Anova(glm.ara.nb, type="II") glm.ara.nb2 <- glm.nb(seedlings ~ dungs+property, ara) Anova(glm.ara.nb2, type="II") par(mfrow=c(2,2)) plot(glm.ara.nb2) library(aods3) gof(glm.ara.nb2)