Commit 6bf9da25 authored by Eduardo E. R. Junior's avatar Eduardo E. R. Junior

Adiciona à vignette compoisson a análise dos dados 'soja'

parent aaa36c34
Pipeline #4086 passed with stage
## ----setup, include=FALSE-----------------------------------------
## ----setup, include=FALSE------------------------------------------------
source("_setup.R")
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
library(MRDCr)
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
llcmp
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
cmp
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
set.seed(2016)
x <- rep(1:3, 2)
......@@ -49,127 +49,120 @@ mle2(llcmp, start = start,
vecpar = TRUE)
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
data(capdesfo)
str(capdesfo)
## help(capdesfo)
data(capmosca)
str(capmosca)
## help(capmosca)
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
capmosca <- aggregate(ncap ~ vaso + dexp, data = capmosca, FUN = sum)
str(capmosca)
## ------------------------------------------------------------------------
## Experimento balanceado
xtabs(~est + des, data = capdesfo)
xtabs(~dexp, data = capmosca)
library(lattice)
library(latticeExtra)
(xy <- xyplot(ncap ~ des | est,
data = capdesfo,
xlab = "Nível de desfolha artificial",
ylab = "Número de capulhos produzidos",
type = c("p", "g", "smooth"),
panel = panel.beeswarm,
r = 0.05))
(xy <- xyplot(ncap ~ dexp,
data = capmosca,
xlab = "Dias de infestação",
ylab = "Número de capulhos produzidos",
type = c("p", "g", "smooth"),
panel = panel.beeswarm,
r = 0.05))
## Avaliando preliminarmente suposição de equidispersão
(mv <- aggregate(ncap ~ est + des, data = capdesfo,
(mv <- aggregate(ncap ~ dexp, data = capmosca,
FUN = function(x) c(mean = mean(x), var = var(x))))
xlim <- ylim <- extendrange(c(mv$ncap), f = 0.05)
xyplot(ncap[, "var"] ~ ncap[, "mean"],
data = mv,
xlim = xlim,
ylim = ylim,
ylab = "Variância Amostral",
xlab = "Média Amostral",
panel = function(x, y) {
panel.xyplot(x, y, type = c("p", "r"), grid = TRUE)
panel.abline(a = 0, b = 1, lty = 2)
})
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
## Preditores considerados
f1 <- ncap ~ 1
f2 <- ncap ~ des + I(des^2)
f3 <- ncap ~ est:des + I(des^2)
f4 <- ncap ~ est:(des + I(des^2))
f2 <- ncap ~ dexp
f3 <- ncap ~ dexp + I(dexp^2)
## Ajustando os modelos Poisson
m1P <- glm(f1, data = capdesfo, family = poisson)
m2P <- glm(f2, data = capdesfo, family = poisson)
m3P <- glm(f3, data = capdesfo, family = poisson)
m4P <- glm(f4, data = capdesfo, family = poisson)
m1P <- glm(f1, data = capmosca, family = poisson)
m2P <- glm(f2, data = capmosca, family = poisson)
m3P <- glm(f3, data = capmosca, family = poisson)
## Ajustando os modelos COM-Poisson
m1C <- cmp(f1, data = capdesfo)
m2C <- cmp(f2, data = capdesfo)
m3C <- cmp(f3, data = capdesfo)
m4C <- cmp(f4, data = capdesfo)
m1C <- cmp(f1, data = capmosca)
m2C <- cmp(f2, data = capmosca)
m3C <- cmp(f3, data = capmosca)
## ------------------------------------------------------------------------
## -----------------------------------------------------------------
library(bbmle)
## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(list(m1P, m2P, m3P, m4P), logLik),
"COM-Poisson" = sapply(list(m1C, m2C, m3C, m4C), logLik))
cbind("Poisson" = sapply(list(m1P, m2P, m3P), logLik),
"COM-Poisson" = sapply(list(m1C, m2C, m3C), logLik))
## Teste de razão de verossimilhanças
anova(m1P, m2P, m3P, m4P, test = "Chisq")
anova(m1C, m2C, m3C, m4C)
anova(m1P, m2P, m3P, test = "Chisq")
anova(m1C, m2C, m3C)
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
## Estimativas dos parâmetros
summary(m4P)
summary(m4C)
summary(m3P)
summary(m3C)
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
## Um dos problemas computacionais do modelo COM-Poisson é a obtenção da
## 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(m4C)
convergencez(m3C)
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
## 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))
plot(m4P)
plot(m3P)
## ---- cache = TRUE------------------------------------------------
## ---- cache = TRUE-------------------------------------------------------
##-------------------------------------------
## Testando a nulidade do parâmetro phi
## Usando o ajuste Poisson
trv <- 2 * (logLik(m4C) - logLik(m4P))
trv <- 2 * (logLik(m3C) - logLik(m3P))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)
## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m4Cfixed <- cmp(f4, data = capdesfo, fixed = list("phi" = 0))
anova(m4C, m4Cfixed)
m3Cfixed <- cmp(f3, data = capmosca, fixed = list("phi" = 0))
anova(m3C, m3Cfixed)
## Via perfil de log-verossimilhança
perf <- profile(m4C, which = 1)
perf <- profile(m3C, which = 1)
confint(perf)
plot(perf)
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
##-------------------------------------------
## Verificando a matriz ve variâncias e covariâncias
Vcov <- vcov(m4C)
Vcov <- vcov(m3C)
Corr <- cov2cor(Vcov)
library(corrplot)
......@@ -178,71 +171,477 @@ corrplot.mixed(Corr, lower = "number", upper = "ellipse",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
## -----------------------------------------------------------------
## ------------------------------------------------------------------------
## Predição pontual
pred <- with(capdesfo,
## Predição pontual/intervalar
pred <- with(capmosca,
expand.grid(
est = levels(est),
des = seq(min(des), max(des), l = 20)
dexp = seq(min(dexp), max(dexp), l = 50)
))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)
##-------------------------------------------
## Considerando a Poisson
mediasP <- exp(predict(m4P, newdata = pred))
aux <- data.frame(modelo = "Poisson", fit = mediasP)
aux <- predict(m3P, newdata = pred, se.fit = TRUE)
aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*")))
aux <- data.frame(modelo = "Poisson", aux)
predP <- cbind(pred, aux)
##-------------------------------------------
## Considerando a COM-Poisson
f4; f4[[2]] <- NULL; f4
X <- model.matrix(f4, data = pred)
f3; f3[-2]
X <- model.matrix(f3[-2], data = pred)
## Obtendo os parâmetros da distribuição (lambdas e phi)
betas <- coef(m4C)[-1]
phi <- coef(m4C)[1]
betas <- coef(m3C)[-1]
phi <- coef(m3C)[1]
loglambdas <- X %*% betas
## Aplicando a "inversa da função de ligação", ou seja, obtendo as
## contagens médias
mediasC <- sapply(loglambdas, FUN = function(p) {
y <- 0:50; py <- dcmp(y, p, phi, sumto = 100)
sum(y*py)
})
aux <- data.frame(modelo = "COM-Poisson", fit = mediasC)
## 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 <- 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 <- cbind(pred, aux)
##-------------------------------------------
## Visualizando os valores preditos pelos dois modelos
## Visualizando os valores preditos intervalares pelos dois modelos
da <- rbind(predP, predC)
update(xy, type = c("p", "g")) +
as.layer(xyplot(fit ~ des | est,
## 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("Poisson", "COM-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 ~ dexp,
groups = modelo,
data = da, type = "l"))
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(soja)
str(soja)
## help(soja)
## ------------------------------------------------------------------------
## Observação 74 foi diagnosticada como outlier
soja <- soja[-74, ]
soja <- transform(soja, K = factor(K))
## ---- fig.height = 12----------------------------------------------------
## Experimento (des)balanceado
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"))
## -----------------------------------------------------------------
print(xy1, split = c(1, 1, 1, 2), more = TRUE)
print(xy2, split = c(1, 2, 1, 2), more = FALSE)
## Predição intervalar
qn <- qnorm(0.975) * c(lwr = -1, upr = 1)
## ------------------------------------------------------------------------
## Avaliando preliminarmente suposição de equidispersão
(mv <- aggregate(cbind(ngra, nvag, ng2v = ngra/nvag) ~ K + umid,
data = soja, FUN = function(x)
c(mean = mean(x), var = var(x))))
##-------------------------------------------
## Considerando a Poisson
aux <- predict(m4P, newdata = pred, se.fit = TRUE)
aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*")))
predP <- cbind(predP, aux)
## Para o número de grãos
xlim <- ylim <- extendrange(c(mv$ngra), f = 0.05)
mvp1 <- xyplot(ngra[, "var"] ~ ngra[, "mean"],
data = mv,
xlim = xlim, ylim = ylim,
ylab = "Variância Amostral",
xlab = "Média Amostral",
main = "Número de grãos por parcela",
panel = function(x, y) {
panel.xyplot(x, y, type = c("p", "r"), grid = TRUE)
panel.abline(a = 0, b = 1, lty = 2)
})
##-------------------------------------------
## Para o número de vagens
xlim <- ylim <- extendrange(c(mv$nvag), f = 0.05)
mvp2 <- xyplot(nvag[, "var"] ~ nvag[, "mean"],
data = mv,
xlim = xlim, ylim = ylim,
ylab = "Variância Amostral",
xlab = "Média Amostral",
main = "Número de vagens viáveis",
panel = function(x, y) {
panel.xyplot(x, y, type = c("p", "r"), grid = TRUE)
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)
## ------------------------------------------------------------------------
##-------------------------------------------
## Para o número de vagens viáveis por parcela
## Preditores considerados
f1.nv <- nvag ~ bloc + umid + K
f2.nv <- nvag ~ bloc + umid * K
## Ajustando os modelos Poisson
m1P.nv <- glm(f1.nv, data = soja, family = poisson)
m2P.nv <- glm(f2.nv, data = soja, family = poisson)
## Ajustando os modelos COM-Poisson
m1C.nv <- cmp(f1.nv, data = soja, sumto = 300)
m2C.nv <- cmp(f2.nv, data = soja, sumto = 300)
##-------------------------------------------
## Para o número grão produzidos por parcela
## Preditores considerados
f1.ng <- ngra ~ bloc + umid + K
f2.ng <- ngra ~ bloc + umid * K
## Ajustando os modelos Poisson
m1P.ng <- glm(f1.ng, data = soja, family = poisson)
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)
## ------------------------------------------------------------------------
##-------------------------------------------
## 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),
logLik),
"COM-Poisson" = sapply(
list(m1C.nv, m2C.nv, m1C.ng, m2C.ng, m1C.gpv, m2C.gpv),
logLik))
##-------------------------------------------
## Teste de razão de verossimilhanças
## Dos modelos para o número de vagens viáveis
anova(m1P.nv, m2P.nv, test = "Chisq")
anova(m1C.nv, m2C.nv)
## Dos modelos para o número de grão por parcela
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)
## ------------------------------------------------------------------------
##-------------------------------------------
## Estimativas dos parâmetros do modelo proposto para
## o número de vagens viáveis
summary(m2P.nv)
summary(m2C.nv)
## o número de grão por parcela
summary(m2P.ng)
summary(m2C.ng)
## o número de grão por vagem
summary(m1P.gpv)
summary(m1C.gpv)
## ------------------------------------------------------------------------
## Um dos problemas computacionais do modelo COM-Poisson é a obtenção da
## 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
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)
## ------------------------------------------------------------------------
## 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
## Avaliação do modelo para o número de vagens viáveis
par(mfrow = c(2, 2))
plot(m2P.nv)
## Avaliação do modelo para o número de grão por parcela
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-------------------------------------------------------
## Testando a nulidade do parâmetro phi
##-------------------------------------------
## No modelo para o número de vagens viáveis
## Usando o ajuste Poisson
trv <- 2 * (logLik(m2C.nv) - logLik(m2P.nv))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)
## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m2Cfixed.nv <- cmp(f2.nv, data = soja, fixed = list("phi" = 0),
sumto = 300)
anova(m2C.nv, m2Cfixed.nv)
##-------------------------------------------
## No modelo para o número de grão por parcela
## Usando o ajuste Poisson
trv <- 2 * (logLik(m2C.ng) - logLik(m2P.ng))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)
## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
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)
##-------------------------------------------
## 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)
## ------------------------------------------------------------------------
## Verificando a matriz de variâncias e covariâncias
##-------------------------------------------
## No modelo para o número de vagens viáveis
Vcov.nv <- vcov(m2C.nv)
Corr <- cov2cor(Vcov.nv)
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 parcela
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)])
## ---- fig.height = 18----------------------------------------------------
## Predição pontual/intervalar
pred <- with(soja,
expand.grid(
bloc = factor(levels(bloc)[1], levels = levels(bloc)),
umid = levels(umid),
K = levels(K)
## nvag = 1
))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)
## Definindos as matrizes de delineamento
## Do modelo com interação
f2.ng; f2.ng[-2]
X2 <- model.matrix(f2.ng[-2], data = pred)
## Como não temos interesse na interpretação dos efeitos de blocos
## tomaremos a média desses efeitos para predição
bl2 <- attr(X2, "assign") == 1
X2[, bl2] <- X2[, bl2] * 0 + 1/(sum(bl2) + 1)
head(X2)
## Do modelo aditivo, apenas retiramos os termos referentes ao efeito de
## interação
X1 <- X2[, attr(X2, "assign") != 4]
head(X1)
library(multcomp)
##-------------------------------------------
## Considerando a Poisson
## No modelo para o número de vagens
aux <- exp(confint(glht(m2P.nv, linfct = X2),
calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Poisson", aux)
predP.nv <- cbind(pred, aux)
## No modelo para o número de grãos por parcela
aux <- exp(confint(glht(m2P.ng, linfct = X2),
calpha = univariate_calpha())$confint)
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)