diff --git a/inst/doc/v05_compoisson.R b/inst/doc/v05_compoisson.R index 7a51b159d47ad0bc5ece60b6f1a835103e596eda..7c52b22bf334468de74a78b5de62909098e2d72c 100644 --- a/inst/doc/v05_compoisson.R +++ b/inst/doc/v05_compoisson.R @@ -19,13 +19,12 @@ cmp ## ------------------------------------------------------------------------ set.seed(2016) -x <- rep(1:3, 2) -t <- rnorm(6, 5) -y <- rpois(6, lambda = x*t) -(da <- data.frame(x, t, y)) +x <- rep(1:3, 3) +y <- rpois(9, lambda = x) +(da <- data.frame(x, y)) ## Definindo o prditor do modelo -formula <- y ~ x + I(x^2) + offset(log(t)) +formula <- y ~ x + I(x^2) ##------------------------------------------- ## O framework @@ -34,18 +33,17 @@ formula <- y ~ x + I(x^2) + offset(log(t)) frame <- model.frame(formula, data = da) (X <- model.matrix(formula, data = da)) (y <- model.response(frame)) -(o <- model.offset(frame)) ## Utiliza como valores iniciais as estimativas dos parametros de um ## modelo GLM-Poisson -m0 <- glm.fit(x = X, y = y, offset = o, family = poisson()) -start <- c(phi = 0, m0$coefficients) +m0 <- glm.fit(x = X, y = y, family = poisson()) +(start <- c(phi = 0, m0$coefficients)) ## Otimiza a função de log-verossimilhança via bbmle library(bbmle) parnames(llcmp) <- names(start) mle2(llcmp, start = start, - data = list(y = y, X = X, offset = o, sumto = 50), + data = list(y = y, X = X, sumto = 50), vecpar = TRUE) @@ -103,8 +101,6 @@ m3C <- cmp(f3, data = capmosca) ## ------------------------------------------------------------------------ -library(bbmle) - ## Verossimilhancas dos modelos ajustados cbind("Poisson" = sapply(list(m1P, m2P, m3P), logLik), "COM-Poisson" = sapply(list(m1C, m2C, m3C), logLik)) @@ -129,6 +125,9 @@ summary(m3C) ## necessária convergencez(m3C) +## Reajustando o modelo +logLik(m3C <- cmp(f3, data = capmosca, sumto = 30)) + ## ------------------------------------------------------------------------ @@ -211,7 +210,8 @@ aux <- c(loglambdas) + outer(se, qn, FUN = "*") aux <- apply(aux, MARGIN = 2, FUN = function(col) { sapply(col, FUN = function(p) { - y <- 0:30; py <- dcmp(y, p, phi, sumto = 50) + y <- 0:30 + py <- dcmp(y, exp(p), exp(phi), sumto = 50) sum(y*py) }) }) @@ -227,7 +227,7 @@ cols <- trellis.par.get("superpose.line")$col[1:2] key <- list( lines = list(lty = 1, col = cols), rect = list(col = cols, alpha = 0.1, lty = 3, border = NA), - text = list(c("Poisson", "COM-Poisson"))) + text = list(c("COM-Poisson", "Poisson"))) ## Gráfico dos valores preditos e IC de 95% para média update(xy, type = c("p", "g"), key = key, alpha = 0.7) + @@ -266,32 +266,18 @@ xtabs(~K + umid, data = soja) key <- list( type = "b", divide = 1, ## points = list(pch = 1, col = cols), - lines = list(pch = 1, lty = 1, - col = c(cols, trellis.par.get("dot.symbol")$col)), - text = list(c("Nº de grãos por parcela", "Nº de vagens viáveis", - "Nº de grãos por vagem"))) - -xy1 <- xyplot(ngra + nvag ~ K | umid, - data = soja, - xlab = "Nível de adubação potássica", - ylab = "Contagem", - type = c("p", "g", "smooth"), - key = key, - layout = c(NA, 1), - strip = strip.custom( - strip.names = TRUE, var.name = "Umidade")) - -xy2 <- xyplot(ngra/nvag ~ K | umid, - data = soja, - xlab = "Nível de adubação potássica", - ylab = "Contagem", - type = c("p", "g", "smooth"), - layout = c(NA, 1), - strip = strip.custom( - strip.names = TRUE, var.name = "Umidade")) + lines = list(pch = 1, lty = 1, col = cols), + text = list(c("Nº de grãos por parcela", "Nº de vagens viáveis"))) -print(xy1, split = c(1, 1, 1, 2), more = TRUE) -print(xy2, split = c(1, 2, 1, 2), more = FALSE) +xyplot(ngra + nvag ~ K | umid, + data = soja, + xlab = "Nível de adubação potássica", + ylab = "Contagem", + type = c("p", "g", "smooth"), + key = key, + layout = c(NA, 1), + strip = strip.custom( + strip.names = TRUE, var.name = "Umidade")) ## ------------------------------------------------------------------------ @@ -329,26 +315,11 @@ mvp2 <- xyplot(nvag[, "var"] ~ nvag[, "mean"], panel.abline(a = 0, b = 1, lty = 2) }) -##------------------------------------------- -## Para o número de grãos por vagem -xlim <- ylim <- extendrange(c(mv$ng2v), f = 0.05) -mvp3 <- xyplot(ng2v[, "var"] ~ ng2v[, "mean"], - data = mv, - xlim = xlim, ylim = ylim, - ylab = "Variância Amostral", - xlab = "Média Amostral", - main = "Número grãos por vagem", - panel = function(x, y) { - panel.xyplot(x, y, type = c("p", "r"), grid = TRUE) - panel.abline(a = 0, b = 1, lty = 2) - }) - -print(mvp1, split = c(1, 1, 3, 1), more = TRUE) -print(mvp2, split = c(2, 1, 3, 1), more = TRUE) -print(mvp3, split = c(3, 1, 3, 1), more = FALSE) +print(mvp1, split = c(1, 1, 2, 1), more = TRUE) +print(mvp2, split = c(2, 1, 2, 1), more = FALSE) -## ------------------------------------------------------------------------ +## ----ajuste2, cache = TRUE----------------------------------------------- ##------------------------------------------- ## Para o número de vagens viáveis por parcela @@ -380,21 +351,6 @@ m2P.ng <- glm(f2.ng, data = soja, family = poisson) m1C.ng <- cmp(f1.ng, data = soja, sumto = 700) m2C.ng <- cmp(f2.ng, data = soja, sumto = 700) -##------------------------------------------- -## Para o número de grãos produzidos por vagem - -## Preditores considerados -f1.gpv <- ngra ~ offset(log(nvag)) + bloc + umid + K -f2.gpv <- ngra ~ offset(log(nvag)) + bloc + umid * K - -## Ajustando os modelos Poisson -m1P.gpv <- glm(f1.gpv, data = soja, family = poisson) -m2P.gpv <- glm(f2.gpv, data = soja, family = poisson) - -## Ajustando os modelos COM-Poisson -m1C.gpv <- cmp(f1.gpv, data = soja, sumto = 350) -m2C.gpv <- cmp(f2.gpv, data = soja, sumto = 350) - ## ------------------------------------------------------------------------ @@ -402,11 +358,10 @@ m2C.gpv <- cmp(f2.gpv, data = soja, sumto = 350) ## Verossimilhancas dos modelos ajustados cbind("Poisson" = sapply( list("nvagem adtv" = m1P.nv, "nvagem mult" = m2P.nv, - "ngraos adtv" = m1P.ng, "ngraos mult" = m2P.ng, - "ngrapv adtv" = m1P.gpv, "ngrapv mult" = m2P.gpv), + "ngraos adtv" = m1P.ng, "ngraos mult" = m2P.ng), logLik), "COM-Poisson" = sapply( - list(m1C.nv, m2C.nv, m1C.ng, m2C.ng, m1C.gpv, m2C.gpv), + list(m1C.nv, m2C.nv, m1C.ng, m2C.ng), logLik)) ##------------------------------------------- @@ -420,10 +375,6 @@ anova(m1C.nv, m2C.nv) anova(m1P.ng, m2P.ng, test = "Chisq") anova(m1C.ng, m2C.ng) -## Dos modelos para o número de grão por vagem -anova(m1P.gpv, m2P.gpv, test = "Chisq") -anova(m1C.gpv, m2C.gpv) - ## ------------------------------------------------------------------------ ##------------------------------------------- @@ -437,10 +388,6 @@ summary(m2C.nv) summary(m2P.ng) summary(m2C.ng) -## o número de grão por vagem -summary(m1P.gpv) -summary(m1C.gpv) - ## ------------------------------------------------------------------------ @@ -449,17 +396,12 @@ summary(m1C.gpv) ## verificar se o ajuste proporcionou uma densidade válida se faz ## necessária -par(mfrow = c(1, 3)) - ## No modelo para o número de vagens viáveis convergencez(m2C.nv) ## No modelo para o número de grão por parcela convergencez(m2C.ng) -## No modelo para o número de grão por vagem -convergencez(m1C.gpv) - ## ------------------------------------------------------------------------ @@ -474,10 +416,6 @@ plot(m2P.nv) par(mfrow = c(2, 2)) plot(m2P.ng) -## Avaliação do modelo para o número de grão por vagem -par(mfrow = c(2, 2)) -plot(m1P.gpv) - ## ---- cache = TRUE------------------------------------------------------- @@ -509,33 +447,20 @@ m2Cfixed.ng <- cmp(f2.ng, data = soja, fixed = list("phi" = 0), sumto = 700) anova(m2C.ng, m2Cfixed.ng) -##------------------------------------------- -## No modelo para o número de grão por vagem -## Usando o ajuste Poisson -trv <- 2 * (logLik(m1C.gpv) - logLik(m1P.gpv)) -attributes(trv) <- NULL -round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5) - -## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1) -m1Cfixed.gpv <- cmp(f1.gpv, data = soja, fixed = list("phi" = 0), - sumto = 350) -anova(m1C.gpv, m1Cfixed.gpv) +## ----perf2, echo = FALSE, warnings = FALSE------------------------------- ##------------------------------------------- ## Via perfis de log-verossimilhança perf.nv <- profile(m2C.nv, which = "phi") perf.ng <- profile(m2C.ng, which = "phi") -perf.gpv <- profile(m1C.gpv, which = "phi") - -par(mfrow = c(1, 3)) -plot(perf.nv) -plot(perf.ng) -plot(perf.gpv) confint(perf.nv) confint(perf.ng) -confint(perf.gpv) + +par(mfrow = c(1, 2)) +plot(perf.nv) +plot(perf.ng) ## ------------------------------------------------------------------------ @@ -556,14 +481,6 @@ corrplot.mixed(Corr, lower = "number", upper = "ellipse", Vcov.ng <- vcov(m2C.ng) Corr <- cov2cor(Vcov.ng) -corrplot.mixed(Corr, lower = "number", upper = "ellipse", - diag = "l", tl.pos = "lt", tl.col = "black", - tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)]) - -## No modelo para o número de grão por vagem -Vcov.gpv <- vcov(m1C.gpv) -Corr <- cov2cor(Vcov.gpv) - corrplot.mixed(Corr, lower = "number", upper = "ellipse", diag = "l", tl.pos = "lt", tl.col = "black", tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)]) @@ -617,13 +534,6 @@ colnames(aux) <- c("fit", "lwr", "upr") aux <- data.frame(modelo = "Poisson", aux) predP.ng <- cbind(pred, aux) -## No modelo para o número de grãos por vagem -aux <- exp(confint(glht(m1P.gpv, linfct = X1), - calpha = univariate_calpha())$confint) -colnames(aux) <- c("fit", "lwr", "upr") -aux <- data.frame(modelo = "Poisson", aux) -predP.gpv <- cbind(pred, aux) - ##---------------------------------------------------------------------- ## Considerando a COM-Poisson @@ -649,7 +559,8 @@ aux <- c(loglambdas) + outer(se, qn, FUN = "*") aux <- apply(aux, MARGIN = 2, FUN = function(col) { sapply(col, FUN = function(p) { - y <- 0:200; py <- dcmp(y, p, phi, sumto = 300) + y <- 0:200 + py <- dcmp(y, exp(p), exp(phi), sumto = 300) sum(y*py) }) }) @@ -678,42 +589,14 @@ aux <- c(loglambdas) + outer(se, qn, FUN = "*") aux <- apply(aux, MARGIN = 2, FUN = function(col) { sapply(col, FUN = function(p) { - y <- 0:350; py <- dcmp(y, p, phi, sumto = 700) + y <- 0:350 + py <- dcmp(y, exp(p), exp(phi), sumto = 700) sum(y*py) }) }) aux <- data.frame(modelo = "COM-Poisson", aux) predC.ng <- cbind(pred, aux) -##------------------------------------------- -## No modelo para o número de grãos por vagem -## Obtendo os parâmetros da distribuição (lambdas e phi) -betas <- coef(m1C.gpv)[-1] -phi <- coef(m1C.gpv)[1] -loglambdas <- X1 %*% betas - -## Obtendo os erros padrão das estimativas -## Obs.: Deve-se usar a matriz de variâncias e covariâncias -## condicional, pois os parâmetros de locação (betas) e dispersão -## (phi) não são ortogonais. -Vc <- Vcov.gpv[-1, -1] - Vcov.gpv[-1, 1] %*% - solve(Vcov.gpv[1, 1]) %*% Vcov.gpv[1, -1] -U <- chol(Vc) -se <- sqrt(apply(X1 %*% t(U), MARGIN = 1, FUN = function(x) { - sum(x^2) -})) - -aux <- c(loglambdas) + outer(se, qn, FUN = "*") -aux <- apply(aux, MARGIN = 2, - FUN = function(col) { - sapply(col, FUN = function(p) { - y <- 0:30; py <- dcmp(y, p, phi, sumto = 50) - sum(y*py) - }) - }) -aux <- data.frame(modelo = "COM-Poisson", aux) -predC.gpv <- cbind(pred, aux) - ##------------------------------------------- ## Visualizando os valores preditos intervalares pelos dois modelos da.nv <- rbind(predP.nv, predC.nv) @@ -722,9 +605,6 @@ da.nv <- da.nv[with(da.nv, order(umid, K, modelo)), ] da.ng <- rbind(predP.ng, predC.ng) da.ng <- da.ng[with(da.ng, order(umid, K, modelo)), ] -da.gpv <- rbind(predP.gpv, predC.gpv) -da.gpv <- da.gpv[with(da.gpv, order(umid, K, modelo)), ] - source(paste0("https://gitlab.c3sl.ufpr.br/leg/legTools/raw/", "issue%2315/R/panel.segplot.by.R")) @@ -768,50 +648,10 @@ xy2 <- xyplot(ngra ~ K | umid, data = soja, panel = panel.segplot.by, f = 0.1, as.table = TRUE) ) -xy3 <- xyplot(ngra/nvag ~ K | umid, data = soja, - xlab = "Nível de adubação potássica", - ylab = "Contagem", - type = c("p", "g"), alpha = 0.7, - ## key = key, - layout = c(NA, 1), - as.table = TRUE, - strip = strip.custom( - strip.names = TRUE, var.name = "Umidade")) + - as.layer( - segplot( - K ~ lwr + upr | umid, - centers = fit, groups = modelo, data = da.gpv, - grid = TRUE, horizontal = FALSE, draw = FALSE, - lwd = 2, pch = 1:nlevels(da$modelo) + 3, - panel = panel.segplot.by, f = 0.1, as.table = TRUE) - ) - ## x11(width = 10, height = 50) -print(xy1, split = c(1, 1, 1, 3), more = TRUE) -print(xy2, split = c(1, 2, 1, 3), more = TRUE) -print(xy3, split = c(1, 3, 1, 3), more = FALSE) - - -## ---- eval = FALSE------------------------------------------------------- -## -## ## Testando o termo offset -## -## ##====================================================================== -## ## offset 2 e lambda = 5 -## trat <- rep(letters[1:3], 30) -## X <- model.matrix(~trat) -## betas <- c(-0.5, -2, 1.5) -## lambdas <- exp(X %*% betas) -## y <- rpois(nrow(X), lambda = lambdas * 2) -## -## kk <- data.frame(trat = trat, o = 2, y = y) -## kkm <- cmp(y ~ offset(log(o)) + trat, data = kk, sumto = 50) -## kkg <- glm(y ~ offset(log(o)) + trat, data = kk, family = poisson) -## -## cbind(coef(kkm), c(NA, coef(kkg))) -## -## ##====================================================================== -## +print(xy1, split = c(1, 1, 1, 2), more = TRUE) +print(xy2, split = c(1, 2, 1, 2), more = FALSE) + ## ------------------------------------------------------------------------ @@ -896,6 +736,9 @@ summary(m4C) ## necessária convergencez(m4C) +## Reajustando o modelo +logLik(m4C <- cmp(f4, data = capdesfo, sumto = 30)) + ## ------------------------------------------------------------------------ @@ -919,8 +762,11 @@ round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5) m4Cfixed <- cmp(f4, data = capdesfo, fixed = list("phi" = 0)) anova(m4C, m4Cfixed) + +## ----perf3, cache = TRUE------------------------------------------------- + ## Via perfil de log-verossimilhança -perf <- profile(m4C, which = 1) +perf <- profile(m4C, which = "phi") confint(perf) plot(perf) @@ -979,7 +825,8 @@ aux <- c(loglambdas) + outer(se, qn, FUN = "*") aux <- apply(aux, MARGIN = 2, FUN = function(col) { sapply(col, FUN = function(p) { - y <- 0:50; py <- dcmp(y, p, phi, sumto = 100) + y <- 0:50 + py <- dcmp(y, exp(p), exp(phi), sumto = 100) sum(y*py) }) }) @@ -990,6 +837,13 @@ predC <- cbind(pred, aux) ## Visualizando os valores preditos intervalares pelos dois modelos da <- rbind(predP, predC) +## Legenda +cols <- trellis.par.get("superpose.line")$col[1:2] +key <- list( + lines = list(lty = 1, col = cols), + rect = list(col = cols, alpha = 0.1, lty = 3, border = NA), + text = list(c("COM-Poisson", "Poisson"))) + ## Gráfico dos valores preditos e IC de 95% para média update(xy, type = c("p", "g"), key = key, alpha = 0.7) + as.layer(xyplot(fit ~ des | est, @@ -1005,6 +859,124 @@ update(xy, type = c("p", "g"), key = key, alpha = 0.7) + panel = panel.superpose)) +## ------------------------------------------------------------------------ + +## Reparametriza para a média (aproximada) +llcmp3 <- function (params, y, X, offset = NULL, sumto = 50) { + betas <- params[-1] + phi <- params[1] + nu <- exp(phi) + if (is.null(offset)) + offset <- 0 + ##------------------------------------------- + ## Reparametrização para a média + mu <- exp(X %*% betas) + Xb <- nu * log(mu + (nu - 1)/(2 * nu)) + ##------------------------------------------- + i <- 0:sumto + zs <- sapply(Xb, function(loglambda) sum(exp(i * loglambda - + nu * lfactorial(i)))) + Z <- sum(log(zs)) + ll <- sum(y * Xb - nu * lfactorial(y)) - Z + return(-ll) +} + + +## ------------------------------------------------------------------------ + +## Modifica framework para llcmp3 +cmp3 <- function (formula, data, start = NULL, sumto = NULL, ...) { + frame <- model.frame(formula, data) + terms <- attr(frame, "terms") + y <- model.response(frame) + X <- model.matrix(terms, frame) + ## off <- model.offset(frame) + ## if (is.null(sumto)) + ## sumto <- ceiling(max(y)^1.5) + if (is.null(start)) { + m0 <- glm.fit(x = X, y = y, family = poisson()) + start <- c(phi = 0, m0$coefficients) + } + bbmle::parnames(llcmp3) <- names(start) + model <- bbmle::mle2(llcmp3, start = start, data = list(y = y, + X = X), vecpar = TRUE, ...) + return(model) +} + + +## ------------------------------------------------------------------------ + +## Ajustando o modelo +m4C2 <- cmp3(f4, data = capdesfo) +convergencez(m4C2) + + +## ------------------------------------------------------------------------ + +## Perfis de verossimilhança para phi +perf <- profile(m4C2, which = "phi") +confint(perf) +plot(perf) + + +## ------------------------------------------------------------------------ + +##------------------------------------------- +## Verificando a matriz de variâncias e covariâncias +Vcov <- vcov(m4C2) +Corr <- cov2cor(Vcov) + +corrplot.mixed(Corr, lower = "number", upper = "ellipse", + diag = "l", tl.pos = "lt", tl.col = "black", + tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)]) + + +## ------------------------------------------------------------------------ + +X <- m4C2@data$X +mu <- c(exp(X %*% coef(m4C2)[-1])) + +##------------------------------------------- +## Considerando a COM-Poisson reparametrizada para a média +f4; f4[-2] +X <- model.matrix(f4[-2], data = pred) + +betas <- coef(m4C2)[-1] +phi <- coef(m4C2)[1] +logmu <- X %*% betas + +## Obtendo os erros padrão das estimativas +## Obs.: Deve-se usar a matriz de variâncias e covariâncias +## condicional, pois os parâmetros de locação (betas) e dispersão +## (phi) não são ortogonais. +Vc <- Vcov[-1, -1] - Vcov[-1, 1] %*% solve(Vcov[1, 1]) %*% Vcov[1, -1] +U <- chol(Vc) +se <- sqrt(apply(X %*% t(U), MARGIN = 1, FUN = function(x) { + sum(x^2) +})) + +aux <- exp(c(logmu) + outer(se, qn, FUN = "*")) +aux <- data.frame(modelo = "COM-Poisson", aux) +predC2 <- cbind(pred, aux) + +##------------------------------------------- +## Visualizando os valores preditos intervalares pelos dois modelos +da <- rbind(predP, predC2) + +update(xy, type = c("p", "g"), key = key, alpha = 0.7) + + as.layer(xyplot(fit ~ des | est, + groups = modelo, + data = da, + type = "l", + ly = da$lwr, + uy = da$upr, + cty = "bands", + alpha = 0.5, + prepanel = prepanel.cbH, + panel.groups = panel.cbH, + panel = panel.superpose)) + + ## ------------------------------------------------------------------------ data(ninfas) @@ -1060,8 +1032,8 @@ m1P <- glm(f1, data = ninfas, family = poisson) m2P <- glm(f2, data = ninfas, family = poisson) ## Ajustando os modelos COM-Poisson -m1C <- cmp(f1, data = ninfas, sumto = 600) -m2C <- cmp(f2, data = ninfas, sumto = 600) +m1C <- cmp(f1, data = ninfas, sumto = 800) +m2C <- cmp(f2, data = ninfas, sumto = 800) ## ------------------------------------------------------------------------ @@ -1088,7 +1060,7 @@ summary(m1C) ## constante de normalização Z. Assim uma visualização pós ajuste para ## verificar se o ajuste proporcionou uma densidade válida se faz ## necessária -convergencez(m1C, incremento = 100, maxit = 10) +convergencez(m1C) ## ------------------------------------------------------------------------ @@ -1113,8 +1085,11 @@ round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5) m1Cfixed <- cmp(f1, data = ninfas, fixed = list("phi" = 0)) anova(m1C, m1Cfixed) + +## ----perf4, cache = TRUE, warnings = FALSE------------------------------- + ## Via perfil de log-verossimilhança -perf <- profile(m1C, which = 1) +perf <- profile(m1C, which = "phi") confint(perf) plot(perf) @@ -1186,7 +1161,8 @@ aux <- c(loglambdas) + outer(se, qn, FUN = "*") aux <- apply(aux, MARGIN = 2, FUN = function(col) { sapply(col, FUN = function(p) { - y <- 0:400; py <- dcmp(y, p, phi, sumto = 600) + y <- 0:400 + py <- dcmp(y, exp(p), exp(phi), sumto = 600) sum(y*py) }) }) diff --git a/inst/doc/v05_compoisson.html b/inst/doc/v05_compoisson.html index 549848a193a760c71cd899a55aeab8f2c4b5de9c..210f2e2e8a8f2f396c530ffe6634afdd742e9deb 100644 --- a/inst/doc/v05_compoisson.html +++ b/inst/doc/v05_compoisson.html @@ -129,6 +129,11 @@ table.header {
• Comparação dos ajustes
• Avaliando modelo proposto
• Predição
• +
• Reparametrização para a média
• +
• Ajustando o modelo reparametrizado para a média
• Ocorrência de ninfas de mosca-branca em variedades de soja
• Análise Exploratória
• @@ -152,14 +157,12 @@ table.header { \]

em que $$Z(\lambda_i, \nu) = \sum_{j=0}^\infty \lambda_i^j (j!)^{-\nu}$$ e $$\lambda_i = \exp(X_i\beta)$$

llcmp
-
## function (params, y, X, offset = NULL, sumto = ceiling(max(y)^1.5))
+## function (params, y, X, sumto = ceiling(max(y)^1.2))
## {
##     betas <- params[-1]
##     phi <- params[1]
##     nu <- exp(phi)
-##     if (is.null(offset))
-##         offset <- 0
-##     Xb <- X %*% betas + offset
+##     Xb <- X %*% betas
##     i <- 0:sumto
##     zs <- sapply(Xb, function(loglambda) sum(exp(i * loglambda -
##         nu * lfactorial(i))))
@@ -169,9 +172,7 @@ table.header {
## }
## <environment: namespace:MRDCr>
## attr(,"parnames")
-##  [1] "phi"         "(Intercept)" "blocII"      "blocIII"
-##  [5] "blocIV"      "blocV"       "umid50"      "umid62,5"
-##  [9] "K30"         "K60"         "K120"        "K180"
+## [1] "phi"         "(Intercept)" "x"           "I(x^2)"

Detalhes computacionais

• Reparametrização do parâmetro $$\nu$$ para $$\phi = \exp(\nu)$$. Assim o espaçõ paramétrico do modelo são os reais $$\Re^n$$.

• @@ -195,34 +196,38 @@ table.header { ## terms <- attr(frame, "terms") ## y <- model.response(frame) ## X <- model.matrix(terms, frame) -## off <- model.offset(frame) +## if (!is.null(model.offset(frame))) { +## stop("Este modelo ainda nao suporta offset") +## } ## if (is.null(sumto)) -## sumto <- ceiling(max(y)^1.5) +## sumto <- ceiling(max(y)^1.2) ## if (is.null(start)) { -## m0 <- glm.fit(x = X, y = y, offset = off, family = poisson()) +## m0 <- glm.fit(x = X, y = y, family = poisson()) ## start <- c(phi = 0, m0$coefficients) ## } ## bbmle::parnames(llcmp) <- names(start) ## model <- bbmle::mle2(llcmp, start = start, data = list(y = y, -## X = X, offset = off, sumto = sumto), vecpar = TRUE, ...) +## X = X, sumto = sumto), vecpar = TRUE, ...) ## return(model) ## } ## <environment: namespace:MRDCr> Um exemplo de como são construídas as matrizes, definidos os chutes iniciais e ajustados os modelos na função: set.seed(2016) -x <- rep(1:3, 2) -t <- rnorm(6, 5) -y <- rpois(6, lambda = x*t) -(da <- data.frame(x, t, y)) - ## x t y -## 1 1 4.085258 2 -## 2 2 6.001248 16 -## 3 3 4.943577 13 -## 4 1 5.296645 3 -## 5 2 2.208529 5 -## 6 3 4.717260 9 +x <- rep(1:3, 3) +y <- rpois(9, lambda = x) +(da <- data.frame(x, y)) + ## x y +## 1 1 0 +## 2 2 1 +## 3 3 5 +## 4 1 0 +## 5 2 2 +## 6 3 1 +## 7 1 1 +## 8 2 4 +## 9 3 0 ## Definindo o prditor do modelo -formula <- y ~ x + I(x^2) + offset(log(t)) +formula <- y ~ x + I(x^2) ##------------------------------------------- ## O framework @@ -237,34 +242,36 @@ frame <- model.frame(formula, data = da) ## 4 1 1 1 ## 5 1 2 4 ## 6 1 3 9 +## 7 1 1 1 +## 8 1 2 4 +## 9 1 3 9 ## attr(,"assign") ## [1] 0 1 2 (y <- model.response(frame)) - ## 1 2 3 4 5 6 -## 2 16 13 3 5 9 - (o <- model.offset(frame)) - ## [1] 1.4073849 1.7919674 1.5980892 1.6670736 0.7923267 1.5512280 + ## 1 2 3 4 5 6 7 8 9 +## 0 1 5 0 2 1 1 4 0 ## Utiliza como valores iniciais as estimativas dos parametros de um ## modelo GLM-Poisson -m0 <- glm.fit(x = X, y = y, offset = o, family = poisson()) -start <- c(phi = 0, m0$coefficients)
-
-## Otimiza a função de log-verossimilhança via bbmle
+m0 <- glm.fit(x = X, y = y, family = poisson())
+(start <- c(phi = 0, m0$coefficients)) + ## phi (Intercept) x I(x^2) +## 0.000000 -5.144583 5.096001 -1.050030 + ## Otimiza a função de log-verossimilhança via bbmle library(bbmle) parnames(llcmp) <- names(start) mle2(llcmp, start = start, - data = list(y = y, X = X, offset = o, sumto = 50), + data = list(y = y, X = X, sumto = 50), vecpar = TRUE) ## ## Call: ## mle2(minuslogl = llcmp, start = start, data = list(y = y, X = X, -## offset = o, sumto = 50), vecpar = TRUE) +## sumto = 50), vecpar = TRUE) ## ## Coefficients: ## phi (Intercept) x I(x^2) -## 0.3338908 -4.5188932 5.4556688 -1.1173870 +## -0.5769066 -4.4398609 4.0697971 -0.8354618 ## -## Log-likelihood: -11.13 +## Log-likelihood: -13.44 @@ -343,9 +350,7 @@ m3C <- cmp(f3, data = capmosca) ## Comparação dos ajustes - library(bbmle) - -## Verossimilhancas dos modelos ajustados +## Verossimilhancas dos modelos ajustados cbind("Poisson" = sapply(list(m1P, m2P, m3P), logLik), "COM-Poisson" = sapply(list(m1C, m2C, m3C), logLik)) ## Poisson COM-Poisson @@ -404,14 +409,14 @@ summary(m3P) ## ## Call: ## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, -## X = X, offset = off, sumto = sumto), vecpar = TRUE) +## X = X, sumto = sumto), vecpar = TRUE) ## ## Coefficients: ## Estimate Std. Error z value Pr(z) -## phi 1.125858 0.261906 4.2987 1.718e-05 *** -## (Intercept) 6.824427 1.817389 3.7551 0.0001733 *** -## dexp -0.499138 0.266245 -1.8747 0.0608297 . -## I(dexp^2) 0.084627 0.050231 1.6848 0.0920312 . +## phi 1.125801 0.261944 4.2979 1.724e-05 *** +## (Intercept) 6.824074 1.817549 3.7545 0.0001737 *** +## dexp -0.499139 0.266247 -1.8747 0.0608317 . +## I(dexp^2) 0.084627 0.050231 1.6848 0.0920341 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -424,7 +429,10 @@ summary(m3P) ## verificar se o ajuste proporcionou uma densidade válida se faz ## necessária convergencez(m3C) - + + ## Reajustando o modelo +logLik(m3C <- cmp(f3, data = capmosca, sumto = 30)) + ## 'log Lik.' -56.49274 (df=4) ## Dado que o modelo COM-Poisson leva as mesmas estimativas pontuais que ## o modelo Poisson a análise de diagnóstico padrão pode ser utilizada par(mfrow = c(2, 2)) @@ -509,7 +517,8 @@ aux <- c(loglambdas) + outer(se, qn, FUN = "*") aux <- apply(aux, MARGIN = 2, FUN = function(col) { sapply(col, FUN = function(p) { - y <- 0:30; py <- dcmp(y, p, phi, sumto = 50) + y <- 0:30 + py <- dcmp(y, exp(p), exp(phi), sumto = 50) sum(y*py) }) }) @@ -525,7 +534,7 @@ cols <- trellis.par.get("superpose.line")$col[1:2]
key <- list(
lines = list(lty = 1, col = cols),
rect = list(col = cols, alpha = 0.1, lty = 3, border = NA),
-    text = list(c("Poisson", "COM-Poisson")))
+    text = list(c("COM-Poisson", "Poisson")))

## Gráfico dos valores preditos e IC de 95% para média
update(xy, type = c("p", "g"), key = key, alpha = 0.7) +
@@ -540,7 +549,7 @@ update(xy, type = c("p", "g"), key = key, alpha = 0.7) +
prepanel = prepanel.cbH,
panel.groups = panel.cbH,
panel = panel.superpose))
-

+

@@ -580,33 +589,19 @@ xtabs(~K + umid, data = soja)
key <- list(
type = "b", divide = 1,
## points = list(pch = 1, col = cols),
-    lines = list(pch = 1, lty = 1,
-                 col = c(cols, trellis.par.get("dot.symbol")$col)), - text = list(c("Nº de grãos por parcela", "Nº de vagens viáveis", - "Nº de grãos por vagem"))) - -xy1 <- xyplot(ngra + nvag ~ K | umid, - data = soja, - xlab = "Nível de adubação potássica", - ylab = "Contagem", - type = c("p", "g", "smooth"), - key = key, - layout = c(NA, 1), - strip = strip.custom( - strip.names = TRUE, var.name = "Umidade")) - -xy2 <- xyplot(ngra/nvag ~ K | umid, - data = soja, - xlab = "Nível de adubação potássica", - ylab = "Contagem", - type = c("p", "g", "smooth"), - layout = c(NA, 1), - strip = strip.custom( - strip.names = TRUE, var.name = "Umidade")) + lines = list(pch = 1, lty = 1, col = cols), + text = list(c("Nº de grãos por parcela", "Nº de vagens viáveis"))) -print(xy1, split = c(1, 1, 1, 2), more = TRUE) -print(xy2, split = c(1, 2, 1, 2), more = FALSE) - +xyplot(ngra + nvag ~ K | umid, + data = soja, + xlab = "Nível de adubação potássica", + ylab = "Contagem", + type = c("p", "g", "smooth"), + key = key, + layout = c(NA, 1), + strip = strip.custom( + strip.names = TRUE, var.name = "Umidade")) + ## Avaliando preliminarmente suposição de equidispersão (mv <- aggregate(cbind(ngra, nvag, ng2v = ngra/nvag) ~ K + umid, data = soja, FUN = function(x) @@ -671,24 +666,9 @@ mvp2 <- xyplot(nvag[, "var"] ~ nvag[, "mean"], panel.abline(a = 0, b = 1, lty = 2) }) -##------------------------------------------- -## Para o número de grãos por vagem -xlim <- ylim <- extendrange(c(mv$ng2v), f = 0.05)
-mvp3 <- xyplot(ng2v[, "var"] ~ ng2v[, "mean"],
-               data = mv,
-               xlim = xlim, ylim = ylim,
-               ylab = "Variância Amostral",
-               xlab = "Média Amostral",
-               main = "Número grãos por vagem",
-               panel = function(x, y) {
-                   panel.xyplot(x, y, type = c("p", "r"), grid = TRUE)
-                   panel.abline(a = 0, b = 1, lty = 2)
-               })
-
-print(mvp1, split = c(1, 1, 3, 1), more = TRUE)
-print(mvp2, split = c(2, 1, 3, 1), more = TRUE)
-print(mvp3, split = c(3, 1, 3, 1), more = FALSE)
-

+print(mvp1, split = c(1, 1, 2, 1), more = TRUE) +print(mvp2, split = c(2, 1, 2, 1), more = FALSE) +

## Ajuste dos modelos

@@ -720,22 +700,7 @@ m2P.ng <- glm(f2.ng, data = soja, family = poisson) ## Ajustando os modelos COM-Poisson m1C.ng <- cmp(f1.ng, data = soja, sumto = 700) -m2C.ng <- cmp(f2.ng, data = soja, sumto = 700) - -##------------------------------------------- -## Para o número de grãos produzidos por vagem - -## Preditores considerados -f1.gpv <- ngra ~ offset(log(nvag)) + bloc + umid + K -f2.gpv <- ngra ~ offset(log(nvag)) + bloc + umid * K - -## Ajustando os modelos Poisson -m1P.gpv <- glm(f1.gpv, data = soja, family = poisson) -m2P.gpv <- glm(f2.gpv, data = soja, family = poisson) - -## Ajustando os modelos COM-Poisson -m1C.gpv <- cmp(f1.gpv, data = soja, sumto = 350) -m2C.gpv <- cmp(f2.gpv, data = soja, sumto = 350) +m2C.ng <- cmp(f2.ng, data = soja, sumto = 700)

## Comparação dos ajustes

@@ -743,19 +708,16 @@ m2C.gpv <- cmp(f2.gpv, data = soja, sumto = 350) ## Verossimilhancas dos modelos ajustados cbind("Poisson" = sapply( list("nvagem adtv" = m1P.nv, "nvagem mult" = m2P.nv, - "ngraos adtv" = m1P.ng, "ngraos mult" = m2P.ng, - "ngrapv adtv" = m1P.gpv, "ngrapv mult" = m2P.gpv), + "ngraos adtv" = m1P.ng, "ngraos mult" = m2P.ng), logLik), "COM-Poisson" = sapply( - list(m1C.nv, m2C.nv, m1C.ng, m2C.ng, m1C.gpv, m2C.gpv), + list(m1C.nv, m2C.nv, m1C.ng, m2C.ng), logLik))
##               Poisson COM-Poisson
## nvagem mult -259.6165   -259.3267
-## ngraos mult -321.6692   -315.6448
-## ngrapv mult -271.3018   -261.3122
+## ngraos mult -321.6692 -315.6448
##-------------------------------------------
## Teste de razão de verossimilhanças

@@ -807,26 +769,6 @@ anova(m1P.ng, m2P.ng, test = "Chisq")
## 2 20 631.29 21.925 8 0.005057 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -
## Dos modelos para o número de grão por vagem
-anova(m1P.gpv, m2P.gpv, test = "Chisq")
-
## Analysis of Deviance Table
-##
-## Model 1: ngra ~ offset(log(nvag)) + bloc + umid + K
-## Model 2: ngra ~ offset(log(nvag)) + bloc + umid * K
-##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1        63     26.634
-## 2        55     23.985  8    2.649   0.9544
-
anova(m1C.gpv, m2C.gpv)
-
## Likelihood Ratio Tests
-## Model 1: m1C.gpv, [llcmp]: phi+(Intercept)+blocII+blocIII+
-##           blocIV+blocV+umid50+umid62,5+K30+K60+K120+K180
-## Model 2: m2C.gpv, [llcmp]: phi+(Intercept)+blocII+blocIII+
-##           blocIV+blocV+umid50+umid62,5+K30+K60+K120+K180+
-##           umid50:K30+umid62,5:K30+umid50:K60+umid62,5:K60+
-##           umid50:K120+umid62,5:K120+umid50:K180+umid62,5:K180
-##   Tot Df Deviance  Chisq Df Pr(>Chisq)
-## 1     12   532.19
-## 2     20   522.62 9.5648  8     0.2969
##-------------------------------------------
## Estimativas dos parâmetros do modelo proposto para

@@ -876,7 +818,7 @@ summary(m2P.nv)
## ## Call: ## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, -## X = X, offset = off, sumto = sumto), vecpar = TRUE) +## X = X, sumto = sumto), vecpar = TRUE) ## ## Coefficients: ## Estimate Std. Error z value Pr(z) @@ -950,7 +892,7 @@ summary(m2P.ng) ## ## Call: ## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, -## X = X, offset = off, sumto = sumto), vecpar = TRUE) +## X = X, sumto = sumto), vecpar = TRUE) ## ## Coefficients: ## Estimate Std. Error z value Pr(z) @@ -978,64 +920,6 @@ summary(m2P.ng) ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## -2 log L: 631.2897 -
## o número de grão por vagem
-summary(m1P.gpv)
-
##
-## Call:
-## glm(formula = f1.gpv, family = poisson, data = soja)
-##
-## Deviance Residuals:
-##      Min        1Q    Median        3Q       Max
-## -1.60825  -0.36054   0.04036   0.41699   1.25361
-##
-## Coefficients:
-##              Estimate Std. Error z value Pr(>|z|)
-## (Intercept)  0.855552   0.031241  27.385   <2e-16 ***
-## blocII       0.011126   0.026566   0.419    0.675
-## blocIII      0.036245   0.026683   1.358    0.174
-## blocIV       0.018896   0.027174   0.695    0.487
-## blocV        0.012659   0.027784   0.456    0.649
-## umid50      -0.026632   0.021695  -1.228    0.220
-## umid62,5    -0.026250   0.021865  -1.201    0.230
-## K30          0.007448   0.030216   0.247    0.805
-## K60          0.012552   0.029286   0.429    0.668
-## K120         0.039757   0.029467   1.349    0.177
-## K180         0.041632   0.028944   1.438    0.150
-## ---
-## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## (Dispersion parameter for poisson family taken to be 1)
-##
-##     Null deviance: 34.238  on 73  degrees of freedom
-## Residual deviance: 26.634  on 63  degrees of freedom
-## AIC: 567.25
-##
-## Number of Fisher Scoring iterations: 3
-
summary(m1C.gpv)
-
## Maximum likelihood estimation
-##
-## Call:
-## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y,
-##     X = X, offset = off, sumto = sumto), vecpar = TRUE)
-##
-## Coefficients:
-##               Estimate Std. Error z value     Pr(z)
-## phi          0.2687369  0.0740709  3.6281 0.0002855 ***
-## (Intercept)  2.3181295  0.4609397  5.0291 4.927e-07 ***
-## blocII       0.0060163  0.0304179  0.1978 0.8432099
-## blocIII      0.0253652  0.0306999  0.8262 0.4086741
-## blocIV      -0.0115586  0.0325485 -0.3551 0.7224998
-## blocV       -0.0158195  0.0330418 -0.4788 0.6321009
-## umid50       0.0510988  0.0348517  1.4662 0.1425998
-## umid62,5     0.0524112  0.0352500  1.4868 0.1370560
-## K30          0.0894600  0.0430913  2.0761 0.0378886 *
-## K60          0.1403132  0.0523141  2.6821 0.0073155 **
-## K120         0.1802570  0.0556233  3.2407 0.0011925 **
-## K180         0.1879750  0.0566969  3.3154 0.0009150 ***
-## ---
-## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## -2 log L: 532.1893

## Avaliando modelo proposto

@@ -1044,17 +928,12 @@ summary(m1P.gpv) ## verificar se o ajuste proporcionou uma densidade válida se faz ## necessária -par(mfrow = c(1, 3)) - ## No modelo para o número de vagens viáveis -convergencez(m2C.nv) - -## No modelo para o número de grão por parcela -convergencez(m2C.ng) - -## No modelo para o número de grão por vagem -convergencez(m1C.gpv) -

+convergencez(m2C.nv) +

+
## No modelo para o número de grão por parcela
+convergencez(m2C.ng)
+

## Dado que o modelo COM-Poisson leva as mesmas estimativas pontuais que
## o modelo Poisson a análise de diagnóstico padrão pode ser utilizada

@@ -1066,10 +945,6 @@ plot(m2P.nv)
par(mfrow = c(2, 2)) plot(m2P.ng)

-
## Avaliação do modelo para o número de grão por vagem
-par(mfrow = c(2, 2))
-plot(m1P.gpv)
-

## Testando a nulidade do parâmetro phi

##-------------------------------------------
@@ -1122,59 +997,19 @@ anova(m2C.ng, m2Cfixed.ng)
## 2 19 643.34 12.049 1 0.0005183 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 -
##-------------------------------------------
-## No modelo para o número de grão por vagem
-
-## Usando o ajuste Poisson
-trv <- 2 * (logLik(m1C.gpv) - logLik(m1P.gpv))
-attributes(trv) <- NULL
-round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)
-
## [1] 13.06333  0.00030
-
## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
-m1Cfixed.gpv <- cmp(f1.gpv, data = soja, fixed = list("phi" = 0),
-                    sumto = 350)
-anova(m1C.gpv, m1Cfixed.gpv)
-
## Likelihood Ratio Tests
-## Model 1: m1C.gpv, [llcmp]: phi+(Intercept)+blocII+blocIII+
-##           blocIV+blocV+umid50+umid62,5+K30+K60+K120+K180
-## Model 2: m1Cfixed.gpv, [llcmp]: (Intercept)+blocII+blocIII+
-##           blocIV+blocV+umid50+umid62,5+K30+K60+K120+K180
-##   Tot Df Deviance  Chisq Df Pr(>Chisq)
-## 1     12   532.19
-## 2     11   545.25 13.063  1  0.0003011 ***
-## ---
-## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-
##-------------------------------------------
-## Via perfis de log-verossimilhança
-perf.nv <- profile(m2C.nv, which = "phi")
## Warning in onestep(step): Error encountered in profile: Error in optim(par = structure(c(4.49870697260928, -0.0332608749877799,  :
##   initial value in 'vmmin' is not finite
## Warning in onestep(step - 1 + dstep): Error encountered in profile: Error in optim(par = structure(c(4.49870697260928, -0.0332608749877799,  :
##   initial value in 'vmmin' is not finite
-
perf.ng <- profile(m2C.ng, which = "phi")
## Warning in onestep(step): Error encountered in profile: Error in optim(par = structure(c(2.85918209769605, -0.0115732234418501,  :
##   initial value in 'vmmin' is not finite
## Warning in onestep(step - 1 + dstep): Error encountered in profile: Error in optim(par = structure(c(2.85918209769605, -0.0115732234418501,  :
##   initial value in 'vmmin' is not finite
-
perf.gpv <- profile(m1C.gpv, which = "phi")
-
## Warning in onestep(step): Error encountered in profile: Error in optim(par = structure(c(2.31812952440581, 0.00601634110167321,  :
-##   initial value in 'vmmin' is not finite
-
## Warning in onestep(step - 1 + dstep): Error encountered in profile: Error in optim(par = structure(c(2.31812952440581, 0.00601634110167321,  :
-##   initial value in 'vmmin' is not finite
-
par(mfrow = c(1, 3))
-plot(perf.nv)
-plot(perf.ng)
-plot(perf.gpv)
-

-
confint(perf.nv)
##      2.5 %     97.5 %
## -0.2138722  0.4352045
-
confint(perf.ng)
##      2.5 %     97.5 %
## -0.8644610 -0.2164388
-
confint(perf.gpv)
-
##     2.5 %    97.5 %
-## 0.1235600 0.4136233
+

## Verificando a matriz de variâncias e covariâncias

##-------------------------------------------
@@ -1195,14 +1030,6 @@ corrplot.mixed(Corr, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])