Cria chunk apos dev.off. Resíduos de Pearson do PGen vs Pois.

704c977e · Walmes Marques Zeviani · 71508b1a · 704c977e · 704c977e · 704c977e
Commit 704c977e authored 9 years ago by Walmes Marques Zeviani
--- a/inst/doc/v04_poisson_generelizada.R
+++ b/inst/doc/v04_poisson_generelizada.R
@@ -322,6 +322,7 @@ dev.off()
 # Tamanho das covariâncias com \alpha.
 each(sum, mean, max)(abs(V[1, -1]))
+## -----------------------------------------------------------------
 #-----------------------------------------------------------------------
 # Testes de hipótese.
@@ -467,6 +468,7 @@ dev.off()
 # Tamanho das covariâncias com \alpha.
 each(sum, mean, max)(abs(V[1, -1]))
+## -----------------------------------------------------------------
 # Teste de Wald para a interação.
 a <- c(0, attr(model.matrix(m0), "assign"))
 ai <- a == max(a)
@@ -621,6 +623,7 @@ corrplot.mixed(V, lower = "number", upper = "ellipse",
               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]))
@@ -771,6 +774,7 @@ corrplot.mixed(V, lower = "number", upper = "ellipse",
               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]))
@@ -853,165 +857,204 @@ update(p1, type = "p", layout = c(NA, 1),
       key = key, spread = 0.07) +
    as.layer(p2, under = TRUE)
-## ---- eval=FALSE--------------------------------------------------
+## -----------------------------------------------------------------
-#  data(nematoide, package = "MRDCr")
+#-----------------------------------------------------------------------
-#  str(nematoide)
-#  
+data(nematoide, package = "MRDCr")
-#  # Número de nematóides por grama de raíz.
+str(nematoide)
-#  plot(nema ~ off, data = nematoide)
-#  
+# Número de nematóides por grama de raíz.
-#  # Média do número de nematóides por grama de raíz.
+plot(nema ~ off, data = nematoide)
-#  mv <- aggregate(cbind(y = nema/off) ~ cult, data = nematoide,
-#                  FUN = function(x) c(m = mean(x), v = var(x)))
+# Média do número de nematóides por grama de raíz.
-#  
+mv <- aggregate(cbind(y = nema/off) ~ cult, data = nematoide,
-#  xyplot(y[, "v"] ~ y[, "m"], data = mv,
+                FUN = function(x) c(m = mean(x), v = var(x)))
-#         xlab = "Média amostral",
-#         ylab = "Variância amostral") +
+xyplot(y[, "v"] ~ y[, "m"], data = mv,
-#      layer(panel.abline(a = 0, b = 1, lty = 2))
+       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,
+# Ajuste do Poisson.
-#            data = nematoide,
-#            family = poisson)
+m0 <- glm(nema ~ offset(log(off)) + cult,
-#  m1 <- update(m0, family = quasipoisson)
+          data = nematoide,
-#  
+          family = poisson)
-#  # Diagnóstico.
+m1 <- update(m0, family = quasipoisson)
-#  par(mfrow = c(2, 2))
-#  plot(m0); layout(1)
+# Diagnóstico.
-#  
+par(mfrow = c(2, 2))
-#  # Estimativas dos parâmetros.
+plot(m0); layout(1)
-#  summary(m1)
-#  
+# Estimativas dos parâmetros.
-#  # Quadro de deviance.
+summary(m1)
-#  # anova(m0, test = "Chisq")
-#  anova(m1, test = "F")
+# Quadro de deviance.
-#  
+# anova(m0, test = "Chisq")
-#  #-----------------------------------------------------------------------
+anova(m1, test = "F")
-#  # Ajuste da Poisson Generalizada.
-#  
+#-----------------------------------------------------------------------
-#  L <- with(nematoide,
+# Ajuste da Poisson Generalizada.
-#            list(y = nema, offset = off, X = model.matrix(m0)))
-#  
+L <- with(nematoide,
-#  start <- c(alpha = log(1), coef(m0))
+          list(y = nema, offset = off, X = model.matrix(m0)))
-#  parnames(llpgnz) <- names(start)
-#  
+start <- c(alpha = log(1), coef(m0))
-#  # Modelo Poisson também.
+parnames(llpgnz) <- names(start)
-#  m2 <- mle2(llpgnz, start = start, data = L,
-#             fixed = list(alpha = 0), vecpar = TRUE)
+# Modelo Poisson também.
-#  
+m2 <- mle2(llpgnz, start = start, data = L,
-#  c(logLik(m2), logLik(m0))
+           fixed = list(alpha = 0), vecpar = TRUE)
-#  
-#  # Poisson Generalizada.
+c(logLik(m2), logLik(m0))
-#  m3 <- pgnz(formula(m0), data = nematoide)
-#  
+# Poisson Generalizada.
-#  # Diferença de deviance.
+m3 <- pgnz(formula(m0), data = nematoide)
-#  # 2 * diff(c(logLik(m0), logLik(m3)))
-#  anova(m3, m2)
+# Diferença de deviance.
-#  
+# 2 * diff(c(logLik(m0), logLik(m3)))
-#  # Perfil de log-verossimilhança para o parâmetro de dispersão.
+anova(m3, m2)
-#  plot(profile(m3, which = "alpha"))
-#  
+# Perfil de log-verossimilhança para o parâmetro de dispersão.
-#  # Covariância.
+plot(profile(m3, which = "alpha"))
-#  V <- cov2cor(vcov(m3))
-#  corrplot.mixed(V, lower = "number", upper = "ellipse",
+# Covariância.
-#                 diag = "l", tl.pos = "lt", tl.col = "black",
+V <- cov2cor(vcov(m3))
-#                 tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
+corrplot.mixed(V, lower = "number", upper = "ellipse",
-#  dev.off()
+               diag = "l", tl.pos = "lt", tl.col = "black",
-#  
+               tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
-#  # Tamanho das covariâncias com \alpha.
+dev.off()
-#  each(sum, mean, max)(abs(V[1, -1]))
-#  
+## -----------------------------------------------------------------
-#  # Gráfico das estimativas.
+# Tamanho das covariâncias com \alpha.
-#  pars <- data.frame(Pois = c(0, coef(m0)), PGen = coef(m3))
+each(sum, mean, max)(abs(V[1, -1]))
-#  xyplot(PGen ~ Pois, data = pars, aspect = "iso", grid = TRUE) +
-#      layer(panel.abline(a = 0, b = 1, lty = 2))
+# 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, {
+X <- model.matrix(m0)
-#  #     cbind(y = nema,
-#  #           Pois = nematoide$off * exp(X %*% coef(m0)),
+# # Predito do número de nematóides observado (considera o offset).
-#  #           PGen = nematoide$off * exp(X %*% coef(m1)[-1]))
+# with(nematoide, {
-#  # })
+#     cbind(y = nema,
-#  
+#           Pois = nematoide$off * exp(X %*% coef(m0)),
-#  # Predito do número de nematóides por grama de raíz.
+#           PGen = nematoide$off * exp(X %*% coef(m1)[-1]))
-#  pred <- with(nematoide, {
-#      data.frame(y = nema/off,
-#                 Pois = c(exp(X %*% coef(m0))),
-#                 PGen = c(exp(X %*% coef(m3)[-1])))
 # })
-#  str(pred)
-#  
+# Predito do número de nematóides por grama de raíz.
-#  splom(pred) + layer(panel.abline(a = 0, b = 1))
+pred <- with(nematoide, {
-#  
+    data.frame(y = nema/off,
-#  # Correlação predito x observado.
+               Pois = c(exp(X %*% coef(m0))),
-#  cor(pred)
+               PGen = c(exp(X %*% coef(m3)[-1])))
-#  
+})
-#  # Média das observações de das estimativas por cultivar.
+str(pred)
-#  predm <- aggregate(as.matrix(pred) ~ cult, data = nematoide, FUN = mean)
-#  cor(predm[, -1])
+splom(pred) + layer(panel.abline(a = 0, b = 1))
-#  
-#  #-----------------------------------------------------------------------
+# Correlação predito x observado.
-#  # Predição com intervalos de confiança.
+cor(pred)
-#  
-#  pred <- unique(subset(nematoide, select = cult))
+# Média das observações de das estimativas por cultivar.
-#  X <- model.matrix(~cult, data = pred)
+predm <- aggregate(as.matrix(pred) ~ cult, data = nematoide, FUN = mean)
-#  
+cor(predm[, -1])
-#  pred <- list(pois = pred, quasi = pred, pgen = pred)
-#  
+#-----------------------------------------------------------------------
-#  # Quantil normal.
+# Predição com intervalos de confiança.
-#  qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
-#  
+pred <- unique(subset(nematoide, select = cult))
-#  # Preditos pela Poisson.
+X <- model.matrix(~cult, data = pred)
-#  aux <- confint(glht(m0, linfct = X),
-#                 calpha = univariate_calpha())$confint
+pred <- list(pois = pred, quasi = pred, pgen = pred)
-#  colnames(aux)[1] <- "fit"
-#  pred$pois <- cbind(pred$pois, exp(aux))
+# Quantil normal.
-#  
+qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
-#  # Preditos pela Quasi-Poisson.
-#  aux <- confint(glht(m1, linfct = X),
+# Preditos pela Poisson.
-#                 calpha = univariate_calpha())$confint
+aux <- confint(glht(m0, linfct = X),
-#  colnames(aux)[1] <- "fit"
+               calpha = univariate_calpha())$confint
-#  pred$quasi <- cbind(pred$quasi, exp(aux))
+colnames(aux)[1] <- "fit"
-#  
+pred$pois <- cbind(pred$pois, exp(aux))
-#  # Preditos pela Poisson Generalizada.
-#  V <- vcov(m3)
+# Preditos pela Quasi-Poisson.
-#  V <- V[-1, -1]
+aux <- confint(glht(m1, linfct = X),
-#  U <- chol(V)
+               calpha = univariate_calpha())$confint
-#  aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
+colnames(aux)[1] <- "fit"
-#                    FUN = function(x) { sum(x^2) }))
+pred$quasi <- cbind(pred$quasi, exp(aux))
-#  pred$pgen$eta <- c(X %*% coef(m3)[-1])
-#  pred$pgen <- cbind(pred$pgen,
+# Preditos pela Poisson Generalizada.
-#                     apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
+V <- vcov(m3)
-#                           FUN = function(x) {
+V <- V[-1, -1]
-#                               exp(pred$pgen$eta + x)
+U <- chol(V)
-#                           }))
+aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
-#  
+                  FUN = function(x) { sum(x^2) }))
-#  pred <- ldply(pred, .id = "modelo")
+pred$pgen$eta <- c(X %*% coef(m3)[-1])
-#  pred <- arrange(pred, cult, modelo)
+pred$pgen <- cbind(pred$pgen,
-#  
+                   apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
-#  key <- list(type = "o", divide = 1,
+                         FUN = function(x) {
-#              lines = list(pch = 1:nlevels(pred$modelo),
+                             exp(pred$pgen$eta + x)
-#                           lty = 1, col = 1),
+                         }))
-#              text = list(c("Poisson", "Quasi-Poisson",
-#                            "Poisson Generelizada")))
+pred <- ldply(pred, .id = "modelo")
-#  
+pred <- arrange(pred, cult, modelo)
-#  xyplot(nema/off ~ cult, data = nematoide,
-#         key = key,
+key <- list(type = "o", divide = 1,
-#         xlab = "Cultivar de feijão",
+            lines = list(pch = 1:nlevels(pred$modelo),
-#         ylab = "Número de nematóides por grama de raíz") +
+                         lty = 1, col = 1),
-#      as.layer(
+            text = list(c("Poisson", "Quasi-Poisson",
-#          xyplot(fit ~ cult, data = pred,
+                          "Poisson Generelizada")))
-#                 pch = pred$modelo,
-#                 ly = pred$lwr, uy = pred$upr,
+xyplot(nema/off ~ cult, data = nematoide,
-#                 cty = "bars", length = 0,
+       key = key,
-#                 prepanel = prepanel.cbH,
+       xlab = "Cultivar de feijão",
-#                 desloc = 0.25 * scale(as.integer(pred$modelo),
+       ylab = "Número de nematóides por grama de raíz") +
-#                                      scale = FALSE),
+    as.layer(
-#                 panel = panel.cbH))
+        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))
+#-----------------------------------------------------------------------
+# Resíduos de Pearson.
+X <- model.matrix(m0)
+# # Resíduos de Pearson no Poisson.
+# with(nematoide,  {
+#     y <- nema
+#     # haty <- fitted(m0)
+#     haty <- nematoide$off * exp(X %*% coef(m0))
+#     sdy <- sqrt(haty)
+#     cbind((y - haty)/sdy,
+#           residuals(m0, type = "pearson"))
+# })
+# Resíduos de Pearson do Poisson Generalizado.
+rp <- with(nematoide,  {
+    y <- nema
+    alph <- coef(m3)["alpha"]
+    haty <- c(nematoide$off * exp(X %*% coef(m3)[-1]))
+    sdy <- sqrt(haty) * (1 + alph * haty)
+    (y - haty)/sdy
+})
+rp <- stack(data.frame(Pois = residuals(m0, type = "pearson"),
+                       PGen = rp))
+qqmath(~values | ind, data = rp,
+       xlab = "Quantis teóricos",
+       ylab = "Resíduos de Pearson",
+       panel = function(...) {
+           panel.qqmathline(...)
+           panel.qqmath(...)
+       })
--- a/inst/doc/v04_poisson_generelizada.html
+++ b/inst/doc/v04_poisson_generelizada.html
--- a/vignettes/v04_poisson_generelizada.Rmd
+++ b/vignettes/v04_poisson_generelizada.Rmd
@@ -384,7 +384,9 @@ dev.off()
 # Tamanho das covariâncias com \alpha.
 each(sum, mean, max)(abs(V[1, -1]))
+```
+```{r}
 #-----------------------------------------------------------------------
 # Testes de hipótese.
@@ -534,7 +536,9 @@ dev.off()
 # Tamanho das covariâncias com \alpha.
 each(sum, mean, max)(abs(V[1, -1]))
+```
+```{r}
 # Teste de Wald para a interação.
 a <- c(0, attr(model.matrix(m0), "assign"))
 ai <- a == max(a)
@@ -691,7 +695,9 @@ 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()
+```
+```{r}
 # Tamanho das covariâncias com \alpha.
 each(sum, mean, max)(abs(V[1, -1]))
@@ -853,7 +859,9 @@ 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()
+```
+```{r}
 # Tamanho das covariâncias com \alpha.
 each(sum, mean, max)(abs(V[1, -1]))
@@ -940,7 +948,9 @@ update(p1, type = "p", layout = c(NA, 1),
 ## Número de Nematóides em Linhagens de Feijão
-```{r, eval=FALSE}
+```{r}
+#-----------------------------------------------------------------------
 data(nematoide, package = "MRDCr")
 str(nematoide)
@@ -1006,7 +1016,9 @@ 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()
+```
+```{r}
 # Tamanho das covariâncias com \alpha.
 each(sum, mean, max)(abs(V[1, -1]))
@@ -1101,4 +1113,40 @@ xyplot(nema/off ~ cult, data = nematoide,
               desloc = 0.25 * scale(as.integer(pred$modelo),
                                    scale = FALSE),
               panel = panel.cbH))
+#-----------------------------------------------------------------------
+# Resíduos de Pearson.
+X <- model.matrix(m0)
+# # Resíduos de Pearson no Poisson.
+# with(nematoide,  {
+#     y <- nema
+#     # haty <- fitted(m0)
+#     haty <- nematoide$off * exp(X %*% coef(m0))
+#     sdy <- sqrt(haty)
+#     cbind((y - haty)/sdy,
+#           residuals(m0, type = "pearson"))
+# })
+# Resíduos de Pearson do Poisson Generalizado.
+rp <- with(nematoide,  {
+    y <- nema
+    alph <- coef(m3)["alpha"]
+    haty <- c(nematoide$off * exp(X %*% coef(m3)[-1]))
+    sdy <- sqrt(haty) * (1 + alph * haty)
+    (y - haty)/sdy
+})
+rp <- stack(data.frame(Pois = residuals(m0, type = "pearson"),
+                       PGen = rp))
+qqmath(~values | ind, data = rp,
+       xlab = "Quantis teóricos",
+       ylab = "Resíduos de Pearson",
+       panel = function(...) {
+           panel.qqmathline(...)
+           panel.qqmath(...)
+       })
 ```