From 2d719ad85362e9466a54a2e8e13e9e7457bda8aa Mon Sep 17 00:00:00 2001 From: Walmes Zeviani Date: Thu, 12 May 2016 15:29:54 -0300 Subject: [PATCH] Adiciona o rpanel ao Suggests. --- DESCRIPTION | 3 +- vignettes/v04_poisson_generelizada.Rmd | 399 +++++++++++++++++++++---- 2 files changed, 348 insertions(+), 54 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 818371e..221860d 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 2570a2e..181eeca 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} -- GitLab