diff --git a/DESCRIPTION b/DESCRIPTION
index 818371eff285bd5ae03c6fc312fdd84999d8bfaa..221860da0b664ae9825aaac88dafe38a56f9b7f2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -45,6 +45,7 @@ Suggests:
reshape2,
knitr,
rmarkdown,
- shiny
+ shiny,
+ rpanel
VignetteBuilder: knitr
RoxygenNote: 5.0.1
diff --git a/vignettes/v04_poisson_generelizada.Rmd b/vignettes/v04_poisson_generelizada.Rmd
index 2570a2ed1a1d4c06e1034b26e1659b642307d265..181eecacb95b30f4cbff0aaa5750631f4cba10e8 100644
--- a/vignettes/v04_poisson_generelizada.Rmd
+++ b/vignettes/v04_poisson_generelizada.Rmd
@@ -44,10 +44,11 @@ f(y) =
\end{cases}
$$
- - $m$ é maior inteiro positivo para o qual $\theta + m\gamma >
- 0$ quando $\gamma$ é negativo. Usa-se $m \geq 4$ para se ter pelo
- menos 4 pontos no suporte com probabilidade não nula (verificar).
- $\theta > 0$;
+ - $m$ é maior inteiro positivo para o qual $\theta + m\gamma > 0$
+ quando $\gamma$ é negativo. A literatura recomenda usar $m \geq 4$
+ para se ter pelo menos 4 pontos no suporte com probabilidade não
+ nula (verificar).
- $\max\{-1, -\theta/m\} < \gamma < 1$;
- Resolvendo a iniqualdade, tem-se que $m = \left\lfloor
\frac{-\theta}{\gamma} \right\rfloor$ quando $\gamma < 0$;
@@ -231,6 +232,7 @@ levelplot(sum ~ theta + gamma,
layer(panel.abline(h = 0, lty = 2))
#-----------------------------------------------------------------------
+# Gráfico do espaço paramétrico de lambda x alpha.
fun <- Vectorize(vectorize.args = c("lambda", "alpha"),
FUN = function(lambda, alpha) {
@@ -262,7 +264,7 @@ levelplot(sum ~ lambda + alpha,
MRDCr::llpgnz
#-----------------------------------------------------------------------
-# Gerando uma amostra aleatória da Poisson.
+# Gerando uma amostra aleatória da Poisson, para teste.
# Offset = 2, lambda = 3.
y <- rpois(100, lambda = 2 * 3)
@@ -292,19 +294,25 @@ plot(profile(n1))
# Covariância.
V <- cov2cor(vcov(n1))
-corrplot.mixed(V, upper = "ellipse", col = "gray50")
+corrplot.mixed(V, lower = "number", upper = "ellipse",
+ diag = "l", tl.pos = "lt", tl.col = "black",
+ tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()
```
## Número de Vagens Produzidas em Soja ##
-Experimento fatorial (3 x 5), em delineamento de blocos casualizados,
-que estudou a influência de níveis de potássio na adubação em combinação
-com irrigação na produção de soja. As variáveis de contagem registradas
-nesse experimento foram o número de vagens viáveis (e não viáveis) e o
-número total de sementes por parcela.
+Resultados de um experimento fatorial (3 x 5), em delineamento de blocos
+casualizados, que estudou a influência de níveis de potássio na adubação
+de soja em combinação com irrigação em casa de vegetação. As variáveis
+de contagem registradas nesse experimento foram o número de vagens
+viáveis (e não viáveis) e o número total de sementes por parcela com
+duas plantas de soja.
```{r}
+#-----------------------------------------------------------------------
+# Carregando e explorando os dados.
+
data(soja, package = "MRDCr")
str(soja)
@@ -313,6 +321,8 @@ soja <- soja[-74, ]
xyplot(nvag ~ K | umid, data = soja, layout = c(NA, 1),
type = c("p", "smooth"),
+ ylab = "Número de vagens por parcela",
+ xlab = expression("Dose de potássio aplicada"~(mg ~ dm^{3})),
strip = strip.custom(strip.names = TRUE, var.name = "Umidade"))
soja <- transform(soja, K = factor(K))
@@ -321,21 +331,24 @@ soja <- transform(soja, K = factor(K))
# Modelo Poisson.
m0 <- glm(nvag ~ bloc + umid * K, data = soja, family = poisson)
+m1 <- update(m0, family = quasipoisson)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Medidas de decisão.
-anova(m0, test = "Chisq")
-summary(m0)
+# anova(m0, test = "Chisq")
+anova(m1, test = "F")
+summary(m1)
#-----------------------------------------------------------------------
# Modelo Poisson Generalizado.
-L <- with(soja, list(y = nvag, offset = 1, X = model.matrix(m0)))
+L <- with(soja,
+ list(y = nvag, offset = 1, X = model.matrix(m0)))
-# Usa as estimativas do Poisson como valore iniciais.
+# Usa as estimativas do Poisson como valores iniciais para a PGen.
start <- c(alpha = 0, coef(m0))
parnames(llpgnz) <- names(start)
@@ -353,6 +366,7 @@ m3 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE)
# Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0).
anova(m3, m2)
+# Estimativas dos coeficientes.
cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"PGeneraliz" = coef(m3))
@@ -362,12 +376,17 @@ plot(profile(m3, which = "alpha"))
abline(v = 0, lty = 2)
V <- cov2cor(vcov(m3))
-corrplot.mixed(V, upper = "ellipse", col = "gray50")
+corrplot.mixed(V, lower = "number", upper = "ellipse",
+ diag = "l", tl.pos = "lt", tl.col = "black",
+ tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()
# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))
+#-----------------------------------------------------------------------
+# Testes de hipótese.
+
# Teste de Wald para a interação.
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a == max(a)
@@ -383,7 +402,7 @@ crossprod(L %*% coef(m3),
L %*% coef(m3)))
# Teste de Wald para interação (poderia ser LRT, claro).
-# É necessário um objeto glm.
+# É necessário passar um objeto glm mesmo fornecendo o restante a parte.
linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
@@ -416,7 +435,6 @@ str(pred$pois)
# Matrix de covariância completa e sem o alpha (marginal).
V <- vcov(m3)
V <- V[-1, -1]
-
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
@@ -439,7 +457,7 @@ xyplot(fit ~ K | umid, data = pred,
layout = c(NA, 1), as.table = TRUE,
xlim = extendrange(range(pred$K), f = 0.1),
key = key, pch = pred$modelo,
- xlab = expression("Dose de potássion"~(mg~dm^{-3})),
+ xlab = expression("Dose de potássio"~(mg~dm^{-3})),
ylab = "Número de vagens por parcela",
ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
prepanel = prepanel.cbH,
@@ -456,27 +474,32 @@ Análise do número de grãos por pacela do experimento com soja.
xyplot(ngra ~ K | umid, data = soja, layout = c(NA, 1),
type = c("p", "smooth"),
+ ylab = "Número de grãos por parcela",
+ xlab = expression("Dose de potássio aplicada"~(mg ~ dm^{3})),
strip = strip.custom(strip.names = TRUE, var.name = "Umidade"))
#-----------------------------------------------------------------------
# Modelo Poisson.
m0 <- glm(ngra ~ bloc + umid * K, data = soja, family = poisson)
+m1 <- update(m0, family = quasipoisson)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Medidas de decisão.
-anova(m0, test = "Chisq")
-summary(m0)
+# anova(m0, test = "Chisq")
+anova(m1, test = "Chisq")
+summary(m1)
#-----------------------------------------------------------------------
# Modelo Poisson Generalizado.
-L <- with(soja, list(y = ngra, offset = 1, X = model.matrix(m0)))
+L <- with(soja,
+ list(y = ngra, offset = 1, X = model.matrix(m0)))
-# Usa as estimativas do Poisson como valore iniciais.
+# Usa as estimativas do Poisson como valores iniciais.
start <- c(alpha = 0, coef(m0))
parnames(llpgnz) <- names(start)
@@ -493,6 +516,7 @@ m3 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE)
# Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0).
anova(m3, m2)
+# Estimaitvas dos parâmetros.
cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"PGeneraliz" = coef(m3))
@@ -502,7 +526,9 @@ plot(profile(m3, which = "alpha"))
abline(v = 0, lty = 2)
V <- cov2cor(vcov(m3))
-corrplot.mixed(V, upper = "ellipse", col = "gray50")
+corrplot.mixed(V, lower = "number", upper = "ellipse",
+ diag = "l", tl.pos = "lt", tl.col = "black",
+ tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()
# Tamanho das covariâncias com \alpha.
@@ -533,7 +559,7 @@ pred <- attr(X, "grid")
pred <- transform(pred,
K = as.integer(K),
umid = factor(umid))
-pred <- list(pois = pred, pgen = pred)
+pred <- list(pois = pred, quasi = pred, pgen = pred)
# Quantil normal.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
@@ -543,12 +569,16 @@ aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
-str(pred$pois)
-# Matrix de covariância completa e sem o alpha (marginal).
+# Preditos pela Quasi-Poisson.
+aux <- confint(glht(m1, linfct = X),
+ calpha = univariate_calpha())$confint
+colnames(aux)[1] <- "fit"
+pred$quasi <- cbind(pred$quasi, exp(aux))
+
+# Preditos pela Poisson Generalizada.
V <- vcov(m3)
V <- V[-1, -1]
-
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
@@ -559,6 +589,7 @@ pred$pgen <- cbind(pred$pgen,
exp(pred$pgen$eta + x)
}))
+# Junta o resultado dos 3 modelos.
pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, umid, K, modelo)
str(pred)
@@ -566,13 +597,14 @@ str(pred)
key <- list(type = "o", divide = 1,
lines = list(pch = 1:nlevels(pred$modelo),
lty = 1, col = 1),
- text = list(c("Poisson", "Poisson Generelizada")))
+ text = list(c("Poisson", "Quasi-Poisson",
+ "Poisson Generelizada")))
xyplot(fit ~ K | umid, data = pred,
layout = c(NA, 1), as.table = TRUE,
xlim = extendrange(range(pred$K), f = 0.1),
key = key, pch = pred$modelo,
- xlab = expression("Dose de potássion"~(mg~dm^{-3})),
+ xlab = expression("Dose de potássio"~(mg~dm^{-3})),
ylab = "Número de grãos por parcela",
ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
prepanel = prepanel.cbH,
@@ -584,10 +616,12 @@ xyplot(fit ~ K | umid, data = pred,
```{r}
#-----------------------------------------------------------------------
-# Número de grãos por vagem.
+# Número de grãos por vagem (offset).
xyplot(ngra/nvag ~ K | umid, data = soja, layout = c(NA, 1),
type = c("p", "smooth"),
+ xlab = expression("Dose de potássio"~(mg~dm^{-3})),
+ ylab = "Número de grãos por vagem",
strip = strip.custom(strip.names = TRUE, var.name = "Umidade"))
#-----------------------------------------------------------------------
@@ -595,22 +629,33 @@ xyplot(ngra/nvag ~ K | umid, data = soja, layout = c(NA, 1),
m0 <- glm(ngra ~ offset(log(nvag)) + bloc + umid * K,
data = soja, family = poisson)
+m1 <- update(m0, family = quasipoisson)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Medidas de decisão.
-anova(m0, test = "Chisq")
-summary(m0)
+# anova(m0, test = "Chisq")
+anova(m1, test = "F")
+summary(m1)
+
+# Declara um modelo aditivo.
+m0 <- glm(ngra ~ offset(log(nvag)) + bloc + umid + K,
+ data = soja, family = poisson)
+m1 <- update(m0, family = quasipoisson)
+anova(m1, test = "F")
#-----------------------------------------------------------------------
# Modelo Poisson Generalizado.
-L <- with(soja, list(y = ngra, offset = nvag, X = model.matrix(m0)))
+L <- with(soja,
+ list(y = ngra, offset = nvag, X = model.matrix(m0)))
-# Já que na verossimilhaça (1 + alpha * y) > -1, então o menor valor
-# possível para gamma para ter uma log-verossimilhança avaliável é
+# De acordo com a estimativa de phi da Quasi-Poisson, esse dado é
+# subdisperso. Já que na verossimilhaça (1 + alpha * y) > -1 quando
+# alpha < 0, então o menor valor possível para gamma para ter uma
+# log-verossimilhança avaliável é
-1/max(soja$ngra)
# Mesmo com esse lower bound, o valor chute para alpha foi definido por
@@ -639,10 +684,11 @@ cbind("PoissonGLM" = c(NA, coef(m0)),
# Perfil para o parâmetro de dispersão.
plot(profile(m3, which = "alpha"))
-abline(v = 0, lty = 2)
V <- cov2cor(vcov(m3))
-corrplot.mixed(V, upper = "ellipse", col = "gray50")
+corrplot.mixed(V, lower = "number", upper = "ellipse",
+ diag = "l", tl.pos = "lt", tl.col = "black",
+ tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()
# Tamanho das covariâncias com \alpha.
@@ -663,10 +709,84 @@ linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
coef. = coef(m3))
+
+#-----------------------------------------------------------------------
+# Predição com bandas de confiança.
+
+# Por causa do offset, não tem como usar a LSmatrix.
+pred <- unique(subset(soja, select = c("umid", "K")))
+str(pred)
+
+pred <- list(pois = pred, quasi = pred, pgen = pred)
+
+X <- model.matrix(formula(m0)[-2],
+ data = cbind(nvag = 1, bloc = soja$bloc[1], pred))
+i <- grep(x = colnames(X), pattern = "^bloc")
+X[, i] <- X[, i] * 0 + 1/(length(i) + 1)
+head(X)
+
+# Quantil normal.
+qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
+
+# Preditos pela Poisson.
+aux <- confint(glht(m0, linfct = X),
+ calpha = univariate_calpha())$confint
+colnames(aux)[1] <- "fit"
+pred$pois <- cbind(pred$pois, exp(aux))
+
+# Preditos pela Quasi-Poisson.
+aux <- confint(glht(m1, linfct = X),
+ calpha = univariate_calpha())$confint
+colnames(aux)[1] <- "fit"
+pred$quasi <- cbind(pred$quasi, exp(aux))
+
+# Preditos pela Poisson Generalizada.
+V <- vcov(m3)
+V <- V[-1, -1]
+U <- chol(V)
+aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
+ FUN = function(x) { sum(x^2) }))
+pred$pgen$eta <- c(X %*% coef(m3)[-1])
+pred$pgen <- cbind(pred$pgen,
+ apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
+ FUN = function(x) {
+ exp(pred$pgen$eta + x)
+ }))
+
+# Junta o resultado dos 3 modelos.
+pred <- ldply(pred, .id = "modelo")
+pred <- arrange(pred, umid, K, modelo)
+pred$K <- as.numeric(as.character(pred$K))
+
+key <- list(type = "o", divide = 1,
+ lines = list(pch = 1:nlevels(pred$modelo),
+ lty = 1, col = 1),
+ text = list(c("Poisson", "Quasi-Poisson",
+ "Poisson Generelizada")))
+
+xyplot(fit ~ K | umid, data = pred,
+ layout = c(NA, 1), as.table = TRUE,
+ xlim = extendrange(range(pred$K), f = 0.2),
+ key = key, pch = pred$modelo,
+ xlab = expression("Dose de potássio"~(mg~dm^{-3})),
+ ylab = "Número de grãos por parcela",
+ ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
+ prepanel = prepanel.cbH,
+ desloc = 8 * scale(as.integer(pred$modelo), scale = FALSE),
+ panel = panel.cbH)
```
## Número de Capulhos Produzidos em Algodão ##
+Experimento conduzido em casa de vegetação para avaliar o efeito da
+desfolha, em diferentes fases de desenvolvimento do algodão, sobre a
+produção da cultura. As parcelas foram vasos com duas plantas
+que tiveram a área das folhas removidas com uma tesoura, simulando o ataque de
+insetos desfolhadores, nos níveis de 0, 25, 50, 75 e 100% de remoção de
+área foliar. Em cada fase de desenvolvimento (de 5), 5 parcelas sofreram
+desfolha nos níveis mencionados. O número de capulhos ao final do
+experimento foi contado.
+
```{r}
#-----------------------------------------------------------------------
# Número de capulhos em função do nível de desfolha artificial e fase
@@ -689,17 +809,20 @@ p1
m0 <- glm(ncap ~ est * (des + I(des^2)),
data = capdesfo, family = poisson)
+m1 <- update(m0, family = quasipoisson)
par(mfrow = c(2, 2))
plot(m0); layout(1)
anova(m0, test = "Chisq")
-summary(m0)
+anova(m1, test = "F")
+summary(m1)
#-----------------------------------------------------------------------
# Modelo Poisson Generalizada.
-L <- with(capdesfo, list(y = ncap, offset = 1, X = model.matrix(m0)))
+L <- with(capdesfo,
+ list(y = ncap, offset = 1, X = model.matrix(m0)))
start <- c(alpha = log(1), coef(m0))
parnames(llpgnz) <- names(start)
@@ -725,7 +848,9 @@ cbind("PoissonGLM" = c(NA, coef(m0)),
"PGeneraliz" = coef(m3))
V <- cov2cor(vcov(m3))
-corrplot.mixed(V, upper = "ellipse", col = "gray50")
+corrplot.mixed(V, lower = "number", upper = "ellipse",
+ diag = "l", tl.pos = "lt", tl.col = "black",
+ tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()
# Tamanho das covariâncias com \alpha.
@@ -755,8 +880,7 @@ linearHypothesis(model = m0,
pred <- with(capdesfo, expand.grid(est = levels(est),
des = seq(0, 1, by = 0.025)))
X <- model.matrix(formula(m0)[-2], data = pred)
-
-pred <- list(pois = pred, pgen = pred)
+pred <- list(pois = pred, quasi = pred, pgen = pred)
# Quantil normal.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
@@ -766,11 +890,16 @@ aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
-str(pred$pois)
-# Matrix de covariância completa e sem o alpha.
+# Preditos pela Quasi-Poisson.
+aux <- confint(glht(m1, linfct = X),
+ calpha = univariate_calpha())$confint
+colnames(aux)[1] <- "fit"
+pred$quasi <- cbind(pred$quasi, exp(aux))
+
+# Preditos pela Poisson Generalizada.
V <- vcov(m3)
-V <- V[-1,-1]
+V <- V[-1, -1]
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
@@ -780,23 +909,22 @@ pred$pgen <- cbind(pred$pgen,
FUN = function(x) {
exp(pred$pgen$eta + x)
}))
-str(pred$pgen)
pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, est, des, modelo)
-str(pred)
-key <- list(lines = list(
- lty = 1,
- col = trellis.par.get("superpose.line")$col[1:2]),
- text = list(
- c("Poisson", "Poisson Generalizada")))
+key <- list(lines = list(lty = 1),
+ text = list(c("Poisson", "Quasi-Poisson",
+ "Poisson Generelizada")))
+key$lines$col <-
+ trellis.par.get("superpose.line")$col[1:nlevels(pred$modelo)]
p2 <- xyplot(fit ~ des | est, data = pred, groups = modelo,
layout = c(NA, 1), as.table = TRUE,
xlim = extendrange(range(pred$des), f = 0.1),
type = "l", key = key,
- ly = pred$lwr, uy = pred$upr, cty = "bands", alpha = 0.5,
+ ly = pred$lwr, uy = pred$upr,
+ cty = "bands", alpha = 0.35,
prepanel = prepanel.cbH,
panel.groups = panel.cbH,
panel = panel.superpose)
@@ -809,6 +937,171 @@ update(p1, type = "p", layout = c(NA, 1),
as.layer(p2, under = TRUE)
```
+## Número de Nematóides em Linhagens de Feijão
+
+```{r}
+data(nematoide, package = "MRDCr")
+str(nematoide)
+
+# Número de nematóides por grama de raíz.
+plot(nema ~ off, data = nematoide)
+
+# Média do número de nematóides por grama de raíz.
+mv <- aggregate(cbind(y = nema/off) ~ cult, data = nematoide,
+ FUN = function(x) c(m = mean(x), v = var(x)))
+
+xyplot(y[, "v"] ~ y[, "m"], data = mv,
+ xlab = "Média amostral",
+ ylab = "Variância amostral") +
+ layer(panel.abline(a = 0, b = 1, lty = 2))
+
+#-----------------------------------------------------------------------
+# Ajuste do Poisson.
+
+m0 <- glm(nema ~ offset(log(off)) + cult,
+ data = nematoide,
+ family = poisson)
+m1 <- update(m0, family = quasipoisson)
+
+# Diagnóstico.
+par(mfrow = c(2, 2))
+plot(m0); layout(1)
+
+# Estimativas dos parâmetros.
+summary(m1)
+
+# Quadro de deviance.
+# anova(m0, test = "Chisq")
+anova(m1, test = "F")
+
+#-----------------------------------------------------------------------
+# Ajuste da Poisson Generalizada.
+
+L <- with(nematoide,
+ list(y = nema, offset = off, X = model.matrix(m0)))
+
+start <- c(alpha = log(1), coef(m0))
+parnames(llpgnz) <- names(start)
+
+# Modelo Poisson também.
+m2 <- mle2(llpgnz, start = start, data = L,
+ fixed = list(alpha = 0), vecpar = TRUE)
+
+c(logLik(m2), logLik(m0))
+
+# Poisson Generalizada.
+m3 <- pgnz(formula(m0), data = nematoide)
+
+# Diferença de deviance.
+# 2 * diff(c(logLik(m0), logLik(m3)))
+anova(m3, m2)
+
+# Perfil de log-verossimilhança para o parâmetro de dispersão.
+plot(profile(m3, which = "alpha"))
+
+# Covariância.
+V <- cov2cor(vcov(m3))
+corrplot.mixed(V, lower = "number", upper = "ellipse",
+ diag = "l", tl.pos = "lt", tl.col = "black",
+ tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
+dev.off()
+
+# Tamanho das covariâncias com \alpha.
+each(sum, mean, max)(abs(V[1, -1]))
+
+# Gráfico das estimativas.
+pars <- data.frame(Pois = c(0, coef(m0)), PGen = coef(m3))
+xyplot(PGen ~ Pois, data = pars, aspect = "iso", grid = TRUE) +
+ layer(panel.abline(a = 0, b = 1, lty = 2))
+
+#-----------------------------------------------------------------------
+
+X <- model.matrix(m0)
+
+# # Predito do número de nematóides observado (considera o offset).
+# with(nematoide, {
+# cbind(y = nema,
+# Pois = nematoide$off * exp(X %*% coef(m0)),
+# PGen = nematoide$off * exp(X %*% coef(m1)[-1]))
+# })
+
+# Predito do número de nematóides por grama de raíz.
+pred <- with(nematoide, {
+ data.frame(y = nema/off,
+ Pois = c(exp(X %*% coef(m0))),
+ PGen = c(exp(X %*% coef(m3)[-1])))
+})
+str(pred)
+
+splom(pred) + layer(panel.abline(a = 0, b = 1))
+
+# Correlação predito x observado.
+cor(pred)
+
+# Média das observações de das estimativas por cultivar.
+predm <- aggregate(as.matrix(pred) ~ cult, data = nematoide, FUN = mean)
+cor(predm[, -1])
+
+#-----------------------------------------------------------------------
+# Predição com intervalos de confiança.
+
+pred <- unique(subset(nematoide, select = cult))
+X <- model.matrix(~cult, data = pred)
+
+pred <- list(pois = pred, quasi = pred, pgen = pred)
+
+# Quantil normal.
+qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
+
+# Preditos pela Poisson.
+aux <- confint(glht(m0, linfct = X),
+ calpha = univariate_calpha())$confint
+colnames(aux)[1] <- "fit"
+pred$pois <- cbind(pred$pois, exp(aux))
+
+# Preditos pela Quasi-Poisson.
+aux <- confint(glht(m1, linfct = X),
+ calpha = univariate_calpha())$confint
+colnames(aux)[1] <- "fit"
+pred$quasi <- cbind(pred$quasi, exp(aux))
+
+# Preditos pela Poisson Generalizada.
+V <- vcov(m3)
+V <- V[-1, -1]
+U <- chol(V)
+aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
+ FUN = function(x) { sum(x^2) }))
+pred$pgen$eta <- c(X %*% coef(m3)[-1])
+pred$pgen <- cbind(pred$pgen,
+ apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
+ FUN = function(x) {
+ exp(pred$pgen$eta + x)
+ }))
+
+pred <- ldply(pred, .id = "modelo")
+pred <- arrange(pred, cult, modelo)
+
+key <- list(type = "o", divide = 1,
+ lines = list(pch = 1:nlevels(pred$modelo),
+ lty = 1, col = 1),
+ text = list(c("Poisson", "Quasi-Poisson",
+ "Poisson Generelizada")))
+
+xyplot(nema/off ~ cult, data = nematoide,
+ key = key,
+ xlab = "Cultivar de feijão",
+ ylab = "Número de nematóides por grama de raíz") +
+ as.layer(
+ xyplot(fit ~ cult, data = pred,
+ pch = pred$modelo,
+ ly = pred$lwr, uy = pred$upr,
+ cty = "bars", length = 0,
+ prepanel = prepanel.cbH,
+ desloc = 0.25 * scale(as.integer(pred$modelo),
+ scale = FALSE),
+ panel = panel.cbH))
+```
+
## Pacote VGAM ##
```{r, eval=FALSE}