Funções para ajuste dos modelos
- Log-verossimilhança
- Ajuste geral
Capulhos de algodão sob exposição à mosca-branca
Produção de soja sob efeito de umidade e adubação potássica
Capulhos de algodão sob efeito de desfolha
Ocorrência de ninfas de mosca-branca em variedades de soja

library(MRDCr)

Funções para ajuste dos modelos

Log-verossimilhança

Implentando a função de log-verossimilhança do modelo COM-Poisson, definida como:

\[ \sum_i^n y_i \log(\lambda_i) - \nu \sum_i^n \log(y_i!) - \sum_i^n \log(Z(\lambda_i, \nu)) \]

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

Detalhes computacionais

Reparametrização do parâmetro \(\nu\) para \(\phi = \log(\nu)\). Assim o espaço paramétrico do modelo são os reais \(\Re^n\).
Truncamento da série infinita \(Z(\lambda_i)\). sumto é tomado como argumento da função.
Para o cálculo de \(Z(\lambda_i)\) faz-se, minimizando problemas de overflow \[ \sum_{j=0}^\infty \lambda_i^j (j!)^{-\nu} = \sum_{j=0}^\infty \exp \left ( \log \left( \lambda_i^j (j!)^{-\nu} \right ) \right ) = \sum_{j=0}^\infty \exp(i \log(\lambda_i) - \nu \log(i!)) \]

llcmp

## function(params, y, X, sumto = ceiling(max(y)^1.2)){
##     ##-------------------------------------------
##     betas <- params[-1]
##     phi <- params[1]
##     nu <- exp(phi)
##     ##-------------------------------------------
##     Xb <- X %*% betas
##     ##-------------------------------------------
##     ## Obtendo a constante normatizadora Z.
##     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)
## }
## <environment: namespace:MRDCr>

Ajuste geral

Framework implementado em R que utiliza a forma de escrita de preditores no estilo de fórmulas, similar as funções lm, glm.

cmp

## function(formula, data, start = NULL, sumto = NULL, ...) {
##     ##-------------------------------------------
##     ## Constrói as matrizes do modelo
##     frame <- model.frame(formula, data)
##     terms <- attr(frame, "terms")
##     y <- model.response(frame)
##     X <- model.matrix(terms, frame)
##     if(!is.null(model.offset(frame))) {
##         stop("Este modelo ainda nao suporta offset")
##     }
##     if (is.null(sumto)) sumto <- ceiling(max(y)^1.2)
##     ##-------------------------------------------
##     ## Define os chutes iniciais
##     if (is.null(start)) {
##         m0 <- glm.fit(x = X, y = y, family = poisson())
##         start <- c("phi" = 0, m0$coefficients)
##     }
##     ##-------------------------------------------
##     ## Nomeia os parâmetros da função para métodos bbmle
##     bbmle::parnames(llcmp) <- names(start)
##     model <- bbmle::mle2(llcmp, start = start,
##                          data = list(y = y, 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, 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)

##-------------------------------------------
## O framework

## Constrói as matrizes para ajuste do modelo
frame <- model.frame(formula, data = da)
(X <- model.matrix(formula, data = da))

##   (Intercept) x I(x^2)
## 1           1 1      1
## 2           1 2      4
## 3           1 3      9
## 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 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, 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, sumto = 50),
     vecpar = TRUE)

## 
## Call:
## mle2(minuslogl = llcmp, start = start, data = list(y = y, X = X, 
##     sumto = 50), vecpar = TRUE)
## 
## Coefficients:
##         phi (Intercept)           x      I(x^2) 
##  -0.5769066  -4.4398609   4.0697971  -0.8354618 
## 
## Log-likelihood: -13.44

Capulhos de algodão sob exposição à mosca-branca

data(capmosca)
str(capmosca)

## 'data.frame':    60 obs. of  8 variables:
##  $ dexp   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ vaso   : Factor w/ 5 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
##  $ planta : Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
##  $ alt    : num  70.5 67.5 72 70 64 64.5 67.5 61 70 70 ...
##  $ pesocap: num  25.2 NA 28.5 NA 31.8 ...
##  $ nerep  : int  4 4 3 5 5 5 5 5 4 5 ...
##  $ ncap   : int  4 4 2 5 5 5 5 5 4 5 ...
##  $ nnos   : int  13 14 14 15 15 11 14 11 14 15 ...

## help(capmosca)

Experimento conduzido sob delineamento inteiramente casualizado na Universidade Federal da Grande Dourados, cujo objetivo foi avaliar os impactos da exposição de plantas de algodão à alta infestação da praga mosca-branca. No experimento avaliou-se duas plantas por vaso, nesta análise tomaremos como unidade amostral o vaso e o interesse será somente na variável número de capulhos produzidos.

capmosca <- aggregate(ncap ~ vaso + dexp, data = capmosca, FUN = sum)
str(capmosca)

## 'data.frame':    30 obs. of  3 variables:
##  $ vaso: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ dexp: int  0 0 0 0 0 1 1 1 1 1 ...
##  $ ncap: int  8 7 10 10 9 5 7 8 8 11 ...

Assim as variáveis consideradas são definidas como:

dexp: Dias de exposição à alta infestação de mosca-branca;
ncap: Número de capulhos de algodão produzidos ao final do experimento.

Análise Exploratória

## Experimento balanceado
xtabs(~dexp, data = capmosca)

## dexp
## 0 1 2 3 4 5 
## 5 5 5 5 5 5

library(lattice)
library(latticeExtra)

(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 ~ dexp, data = capmosca,
                 FUN = function(x) c(mean = mean(x), var = var(x))))

##   dexp ncap.mean ncap.var
## 1    0       8.8      1.7
## 2    1       7.8      4.7
## 3    2       6.8      5.2
## 4    3       6.8      1.7
## 5    4       7.4      2.3
## 6    5       7.6      2.3

Ajuste dos modelos

## Preditores considerados
f1 <- ncap ~ 1
f2 <- ncap ~ dexp
f3 <- ncap ~ dexp + I(dexp^2)

## Ajustando os modelos 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 = capmosca)
m2C <- cmp(f2, data = capmosca)
m3C <- cmp(f3, data = capmosca)

Comparação dos ajustes

## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(list(m1P, m2P, m3P), logLik),
      "COM-Poisson" = sapply(list(m1C, m2C, m3C), logLik))

##        Poisson COM-Poisson
## [1,] -63.76981   -58.77209
## [2,] -63.52399   -58.13540
## [3,] -62.93611   -56.49274

## Teste de razão de verossimilhanças
anova(m1P, m2P, m3P, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: ncap ~ 1
## Model 2: ncap ~ dexp
## Model 3: ncap ~ dexp + I(dexp^2)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        29     11.993                     
## 2        28     11.501  1  0.49164   0.4832
## 3        27     10.325  1  1.17575   0.2782

anova(m1C, m2C, m3C)

## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)
## Model 2: m2C, [llcmp]: phi+(Intercept)+dexp
## Model 3: m3C, [llcmp]: phi+(Intercept)+dexp+I(dexp^2)
##   Tot Df Deviance  Chisq Df Pr(>Chisq)  
## 1      2   117.54                       
## 2      3   116.27 1.2734  1     0.2591  
## 3      4   112.98 3.2853  1     0.0699 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Estimativas dos parâmetros
summary(m3P)

## 
## Call:
## glm(formula = f3, family = poisson, data = capmosca)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.72211  -0.32208   0.01768   0.34877   1.13445  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.17571    0.13842  15.718   <2e-16 ***
## dexp        -0.16923    0.13586  -1.246    0.213    
## I(dexp^2)    0.02870    0.02638   1.088    0.277    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 11.993  on 29  degrees of freedom
## Residual deviance: 10.325  on 27  degrees of freedom
## AIC: 131.87
## 
## Number of Fisher Scoring iterations: 4

summary(m3C)

## Maximum likelihood estimation
## 
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, 
##     X = X, sumto = sumto), vecpar = TRUE)
## 
## Coefficients:
##              Estimate Std. Error z value     Pr(z)    
## 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
## 
## -2 log L: 112.9855

Avaliando modelo proposto

## 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(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))
plot(m3P)

##-------------------------------------------
## Testando a nulidade do parâmetro phi

## Usando o ajuste Poisson
trv <- 2 * (logLik(m3C) - logLik(m3P))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)

## [1] 12.88675  0.00033

## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m3Cfixed <- cmp(f3, data = capmosca, fixed = list("phi" = 0))
anova(m3C, m3Cfixed)

## Likelihood Ratio Tests
## Model 1: m3C, [llcmp]: phi+(Intercept)+dexp+I(dexp^2)
## Model 2: m3Cfixed, [llcmp]: (Intercept)+dexp+I(dexp^2)
##   Tot Df Deviance  Chisq Df Pr(>Chisq)    
## 1      4   112.98                         
## 2      3   125.84 12.855  1  0.0003366 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Via perfil de log-verossimilhança
perf <- profile(m3C, which = 1)
confint(perf)

##     2.5 %    97.5 % 
## 0.5622126 1.5977357

plot(perf)

##-------------------------------------------
## Verificando a matriz ve variâncias e covariâncias
Vcov <- vcov(m3C)
Corr <- cov2cor(Vcov)

library(corrplot)
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)])

Predição

## Predição pontual/intervalar
pred <- with(capmosca,
             expand.grid(
                 dexp = seq(min(dexp), max(dexp), l = 50)
             ))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)

##-------------------------------------------
## Considerando a Poisson
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
f3; f3[-2]

## ncap ~ dexp + I(dexp^2)

## ~dexp + I(dexp^2)

X <- model.matrix(f3[-2], data = pred)

aux <- predict(m3C, newdata = X, type = "response",
               interval = "confidence")
aux <- data.frame(modelo = "COM-Poisson", aux)
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 ~ dexp,
                    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))

Produção de soja sob efeito de umidade e adubação potássica

data(soja)
str(soja)

## 'data.frame':    75 obs. of  5 variables:
##  $ K   : int  0 30 60 120 180 0 30 60 120 180 ...
##  $ umid: Factor w/ 3 levels "37,5","50","62,5": 1 1 1 1 1 2 2 2 2 2 ...
##  $ bloc: Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ngra: int  136 159 156 171 190 140 193 200 208 237 ...
##  $ nvag: int  56 62 66 68 82 63 86 94 86 97 ...

## help(soja)

Dados resultantes de um experimento fatorial 5 \(\times\) 3 conduzido em casa de vegetação, onde experimentou-se diferentes níveis de umidade do solo e níveis de adubação com potássio em parcelas que continham 2 plantas. O interesse deste estudo foi avaliar o impacto desses fatores na produção de soja, mensurada pelos componentes número grãos e número de vagegns viáveis. Para controle de variação local, as parcelas foram arranjadas em blocos. Como temos somente três níveis de umidade aplicados trabalharemos esta variável como categórica. Outra observação é quanto a observação 74, que foi considerada como outlier pelo pesquisador responsável e portanto retirada da análise.

## Observação 74 foi diagnosticada como outlier
soja <- soja[-74, ]
soja <- transform(soja, K = factor(K))

Abaixo definimos as variáveis que são utilizadas na modelagem.

bloc: Fator categórico para controle local no experimento. Foram controladas as variações de ambiente da casa de vegetação, como exposição ao sol, e mão de obra ao longo do experimento, como irrigação e capina dos vasos e controle de doenças.
umid: Fator de níveis categóricos. Faixas de umidade do solo. Procurou-se por meio de pesagens diárias dos vasos, manter a umidade do solo nas faixas de 35 a 40; 47,5 a 52,5; e 60 a 65 do volume total de poros.
K: fator de níveis métricos. Dose de potássio aplicado ao solo na forma de adubação, em mg dm\(^{-3}\).

Análise Exploratória

## Experimento (des)balanceado
xtabs(~K + umid, data = soja)

##      umid
## K     37,5 50 62,5
##   0      5  5    5
##   30     5  5    5
##   60     5  5    5
##   120    5  5    4
##   180    5  5    5

key <- list(
    type = "b", divide = 1,
    ## points = list(pch = 1, col = cols),
    lines = list(pch = 1, lty = 1, col = cols),
    text = list(c("Nº de grãos por parcela", "Nº de vagens viáveis")))

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)
                     c(mean = mean(x), var = var(x))))

##      K umid ngra.mean  ngra.var nvag.mean  nvag.var   ng2v.mean
## 1    0 37,5  115.8000  293.2000  48.80000  91.20000 2.395310245
## 2   30 37,5  156.0000  370.0000  64.20000  43.70000 2.427565924
## 3   60 37,5  163.4000  584.3000  66.40000 127.30000 2.470499952
## 4  120 37,5  166.8000   81.7000  67.80000  33.20000 2.466046370
## 5  180 37,5  155.6000  508.3000  63.00000 139.00000 2.484103301
## 6    0   50  132.2000   32.7000  55.80000  22.70000 2.377585272
## 7   30   50  186.0000 1423.0000  78.20000 327.70000 2.392742849
## 8   60   50  209.4000  801.8000  89.60000 180.30000 2.347156695
## 9  120   50  214.4000  294.3000  90.00000  15.50000 2.380056874
## 10 180   50  237.0000  940.5000  97.60000 120.30000 2.425452773
## 11   0 62,5  139.4000   88.8000  60.60000  50.80000 2.313648571
## 12  30 62,5  163.8000   92.2000  71.60000  48.80000 2.295677653
## 13  60 62,5  215.6000  508.8000  91.80000  91.70000 2.349956047
## 14 120 62,5  238.7500  566.9167  95.75000  94.91667 2.494250659
## 15 180 62,5  232.4000  500.8000  95.40000  86.30000 2.437854072
##       ng2v.var
## 1  0.033816078
## 2  0.013683073
## 3  0.013660997
## 4  0.010510982
## 5  0.011059128
## 6  0.020876812
## 7  0.011001292
## 8  0.022848924
## 9  0.012783676
## 10 0.012417259
## 11 0.025210453
## 12 0.013469841
## 13 0.006495655
## 14 0.004365371
## 15 0.010088987

##-------------------------------------------
## 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)
               })

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

Ajuste dos modelos

##-------------------------------------------
## 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)

Comparação dos ajustes

##-------------------------------------------
## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(
          list("nvagem adtv" = m1P.nv, "nvagem mult" = m2P.nv,
               "ngraos adtv" = m1P.ng, "ngraos mult" = m2P.ng),
          logLik),
      "COM-Poisson" = sapply(
          list(m1C.nv, m2C.nv, m1C.ng, m2C.ng),
          logLik))

##               Poisson COM-Poisson
## nvagem adtv -266.6918   -266.6018
## nvagem mult -259.6165   -259.3267
## ngraos adtv -343.1633   -326.6072
## ngraos mult -321.6692   -315.6448

##-------------------------------------------
## Teste de razão de verossimilhanças

## Dos modelos para o número de vagens viáveis
anova(m1P.nv, m2P.nv, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: nvag ~ bloc + umid + K
## Model 2: nvag ~ bloc + umid * K
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1        63     79.429                       
## 2        55     65.278  8    14.15  0.07793 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1C.nv, m2C.nv)

## Likelihood Ratio Tests
## Model 1: m1C.nv, [llcmp]: phi+(Intercept)+blocII+blocIII+
##           blocIV+blocV+umid50+umid62,5+K30+K60+K120+K180
## Model 2: m2C.nv, [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   533.20                      
## 2     20   518.65 14.55  8     0.0685 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Dos modelos para o número de grão por parcela
anova(m1P.ng, m2P.ng, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: ngra ~ bloc + umid + K
## Model 2: ngra ~ bloc + umid * K
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)    
## 1        63     167.71                         
## 2        55     124.72  8   42.988 8.83e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1C.ng, m2C.ng)

## Likelihood Ratio Tests
## Model 1: m1C.ng, [llcmp]: phi+(Intercept)+blocII+blocIII+
##           blocIV+blocV+umid50+umid62,5+K30+K60+K120+K180
## Model 2: m2C.ng, [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   653.21                        
## 2     20   631.29 21.925  8   0.005057 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##-------------------------------------------
## Estimativas dos parâmetros do modelo proposto para

## o número de vagens viáveis
summary(m2P.nv)

## 
## Call:
## glm(formula = f2.nv, family = poisson, data = soja)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9681  -0.6393  -0.0195   0.4854   2.3306  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.95369    0.06892  57.363  < 2e-16 ***
## blocII        -0.02928    0.04091  -0.716  0.47414    
## blocIII       -0.07265    0.04136  -1.756  0.07902 .  
## blocIV        -0.12544    0.04194  -2.991  0.00278 ** 
## blocV         -0.10795    0.04297  -2.512  0.01199 *  
## umid50         0.13404    0.08765   1.529  0.12619    
## umid62,5       0.21656    0.08602   2.518  0.01181 *  
## K30            0.27427    0.08493   3.229  0.00124 ** 
## K60            0.30797    0.08432   3.652  0.00026 ***
## K120           0.32883    0.08395   3.917 8.97e-05 ***
## K180           0.25540    0.08528   2.995  0.00275 ** 
## umid50:K30     0.06322    0.11557   0.547  0.58433    
## umid62,5:K30  -0.10747    0.11536  -0.932  0.35152    
## umid50:K60     0.16561    0.11370   1.457  0.14521    
## umid62,5:K60   0.10735    0.11220   0.957  0.33869    
## umid50:K120    0.14920    0.11338   1.316  0.18818    
## umid62,5:K120  0.11839    0.11404   1.038  0.29922    
## umid50:K180    0.30370    0.11361   2.673  0.00751 ** 
## umid62,5:K180  0.19838    0.11256   1.762  0.07800 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 322.684  on 73  degrees of freedom
## Residual deviance:  65.278  on 55  degrees of freedom
## AIC: 557.23
## 
## Number of Fisher Scoring iterations: 4

summary(m2C.nv)

## Maximum likelihood estimation
## 
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, 
##     X = X, sumto = sumto), vecpar = TRUE)
## 
## Coefficients:
##                Estimate Std. Error z value     Pr(z)    
## phi            0.128842   0.166355  0.7745 0.4386355    
## (Intercept)    4.498707   0.753608  5.9696 2.379e-09 ***
## blocII        -0.033261   0.043957 -0.7567 0.4492454    
## blocIII       -0.082558   0.046164 -1.7884 0.0737184 .  
## blocIV        -0.142556   0.050548 -2.8202 0.0047992 ** 
## blocV         -0.122674   0.050100 -2.4486 0.0143425 *  
## umid50         0.152300   0.096746  1.5742 0.1154352    
## umid62,5       0.246080   0.100275  2.4541 0.0141256 *  
## K30            0.311661   0.104125  2.9931 0.0027612 ** 
## K60            0.349970   0.106846  3.2755 0.0010548 ** 
## K120           0.373657   0.108688  3.4379 0.0005863 ***
## K180           0.290233   0.102751  2.8246 0.0047334 ** 
## umid50:K30     0.071910   0.123768  0.5810 0.5612386    
## umid62,5:K30  -0.122097   0.124603 -0.9799 0.3271400    
## umid50:K60     0.188280   0.125152  1.5044 0.1324761    
## umid62,5:K60   0.122069   0.121317  1.0062 0.3143177    
## umid50:K120    0.169663   0.124092  1.3672 0.1715518    
## umid62,5:K120  0.134661   0.123612  1.0894 0.2759844    
## umid50:K180    0.345239   0.133922  2.5779 0.0099397 ** 
## umid62,5:K180  0.225538   0.125683  1.7945 0.0727326 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 518.6533

## o número de grão por parcela
summary(m2P.ng)

## 
## Call:
## glm(formula = f2.ng, family = poisson, data = soja)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6297  -1.0707  -0.0873   0.7428   3.2810  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    4.80196    0.04477 107.254  < 2e-16 ***
## blocII        -0.01939    0.02655  -0.730 0.465260    
## blocIII       -0.03663    0.02667  -1.373 0.169673    
## blocIV        -0.10559    0.02715  -3.889 0.000101 ***
## blocV         -0.09313    0.02788  -3.340 0.000837 ***
## umid50         0.13245    0.05692   2.327 0.019968 *  
## umid62,5       0.18548    0.05623   3.299 0.000972 ***
## K30            0.29799    0.05486   5.432 5.56e-08 ***
## K60            0.34434    0.05432   6.339 2.32e-10 ***
## K120           0.36493    0.05409   6.746 1.52e-11 ***
## K180           0.29542    0.05489   5.383 7.35e-08 ***
## umid50:K30     0.04344    0.07482   0.581 0.561495    
## umid62,5:K30  -0.13669    0.07527  -1.816 0.069349 .  
## umid50:K60     0.11559    0.07361   1.570 0.116354    
## umid62,5:K60   0.09174    0.07289   1.259 0.208190    
## umid50:K120    0.11860    0.07329   1.618 0.105637    
## umid62,5:K120  0.16266    0.07367   2.208 0.027242 *  
## umid50:K180    0.28832    0.07327   3.935 8.33e-05 ***
## umid62,5:K180  0.21569    0.07285   2.961 0.003071 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 761.43  on 73  degrees of freedom
## Residual deviance: 124.72  on 55  degrees of freedom
## AIC: 681.34
## 
## Number of Fisher Scoring iterations: 4

summary(m2C.ng)

## Maximum likelihood estimation
## 
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, 
##     X = X, sumto = sumto), vecpar = TRUE)
## 
## Coefficients:
##                Estimate Std. Error z value     Pr(z)    
## phi           -0.517896   0.164787 -3.1428 0.0016732 ** 
## (Intercept)    2.859182   0.473115  6.0433 1.510e-09 ***
## blocII        -0.011573   0.020603 -0.5617 0.5743000    
## blocIII       -0.021861   0.020915 -1.0452 0.2959233    
## blocIV        -0.063028   0.023386 -2.6951 0.0070365 ** 
## blocV         -0.055592   0.023392 -2.3765 0.0174766 *  
## umid50         0.079129   0.045862  1.7254 0.0844633 .  
## umid62,5       0.110803   0.047094  2.3528 0.0186316 *  
## K30            0.177988   0.051448  3.4596 0.0005411 ***
## K60            0.205658   0.053825  3.8209 0.0001330 ***
## K120           0.217953   0.054973  3.9647 7.348e-05 ***
## K180           0.176455   0.051325  3.4380 0.0005861 ***
## umid50:K30     0.025872   0.057975  0.4463 0.6554086    
## umid62,5:K30  -0.081673   0.059681 -1.3685 0.1711582    
## umid50:K60     0.068923   0.058007  1.1882 0.2347640    
## umid62,5:K60   0.054660   0.057046  0.9582 0.3379736    
## umid50:K120    0.070709   0.057822  1.2229 0.2213771    
## umid62,5:K120  0.096967   0.059122  1.6401 0.1009777    
## umid50:K180    0.172004   0.063282  2.7180 0.0065670 ** 
## umid62,5:K180  0.128635   0.060140  2.1389 0.0324408 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 631.2897

Avaliando modelo proposto

## 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

## 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)

## 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)

## 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)

## [1] 0.57970 0.44643

## 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)

## Likelihood Ratio Tests
## Model 1: m2C.nv, [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
## Model 2: m2Cfixed.nv, [llcmp]: (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     20   518.65                     
## 2     19   519.23 0.5797  1     0.4464

##-------------------------------------------
## 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)

## [1] 12.04879  0.00052

## 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)

## Likelihood Ratio Tests
## Model 1: m2C.ng, [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
## Model 2: m2Cfixed.ng, [llcmp]: (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     20   631.29                         
## 2     19   643.34 12.049  1  0.0005183 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##      2.5 %     97.5 % 
## -0.2138722  0.4352045

##      2.5 %     97.5 % 
## -0.8644610 -0.2164388

## 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)])

Predição

## 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]

## ngra ~ bloc + umid * K

## ~bloc + umid * K

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)

##   (Intercept) blocII blocIII blocIV blocV umid50 umid62,5 K30 K60
## 1           1    0.2     0.2    0.2   0.2      0        0   0   0
## 2           1    0.2     0.2    0.2   0.2      1        0   0   0
## 3           1    0.2     0.2    0.2   0.2      0        1   0   0
## 4           1    0.2     0.2    0.2   0.2      0        0   1   0
## 5           1    0.2     0.2    0.2   0.2      1        0   1   0
## 6           1    0.2     0.2    0.2   0.2      0        1   1   0
##   K120 K180 umid50:K30 umid62,5:K30 umid50:K60 umid62,5:K60
## 1    0    0          0            0          0            0
## 2    0    0          0            0          0            0
## 3    0    0          0            0          0            0
## 4    0    0          0            0          0            0
## 5    0    0          1            0          0            0
## 6    0    0          0            1          0            0
##   umid50:K120 umid62,5:K120 umid50:K180 umid62,5:K180
## 1           0             0           0             0
## 2           0             0           0             0
## 3           0             0           0             0
## 4           0             0           0             0
## 5           0             0           0             0
## 6           0             0           0             0

## Do modelo aditivo, apenas retiramos os termos referentes ao efeito de
## interação
X1 <- X2[, attr(X2, "assign") != 4]
head(X1)

##   (Intercept) blocII blocIII blocIV blocV umid50 umid62,5 K30 K60
## 1           1    0.2     0.2    0.2   0.2      0        0   0   0
## 2           1    0.2     0.2    0.2   0.2      1        0   0   0
## 3           1    0.2     0.2    0.2   0.2      0        1   0   0
## 4           1    0.2     0.2    0.2   0.2      0        0   1   0
## 5           1    0.2     0.2    0.2   0.2      1        0   1   0
## 6           1    0.2     0.2    0.2   0.2      0        1   1   0
##   K120 K180
## 1    0    0
## 2    0    0
## 3    0    0
## 4    0    0
## 5    0    0
## 6    0    0

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)

##----------------------------------------------------------------------
## Considerando a COM-Poisson

## No modelo para o número de vagens
aux <- predict(m2C.nv, newdata = X2, type = "response",
               interval = "confidence")
aux <- data.frame(modelo = "COM-Poisson", aux)
predC.nv <- cbind(pred, aux)

## No modelo para o número de grãos por parcela
aux <- predict(m2C.ng, newdata = X2, type = "response",
               interval = "confidence")
aux <- data.frame(modelo = "COM-Poisson", aux)
predC.ng <- cbind(pred, aux)

##-------------------------------------------
## Visualizando os valores preditos intervalares pelos dois modelos
da.nv <- rbind(predP.nv, predC.nv)
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)), ]

source(paste0("https://gitlab.c3sl.ufpr.br/leg/legTools/raw/",
              "issue%2315/R/panel.segplot.by.R"))

key <- list(type = "o", divide = 1,
            lines = list(pch = 1:nlevels(da.nv$modelo) + 3, lty = 1),
            text = list(c("Poisson", "COM-Poisson")))

xy1 <- xyplot(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.nv,
            grid = TRUE, horizontal = FALSE, draw = FALSE,
            lwd = 2, pch = 1:nlevels(da$modelo) + 3,
            panel = panel.segplot.by, f = 0.1, as.table = TRUE)
    )

xy2 <- xyplot(ngra ~ 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.ng,
            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, 2), more = TRUE)
print(xy2, split = c(1, 2, 1, 2), more = FALSE)

Capulhos de algodão sob efeito de desfolha

data(capdesfo)
str(capdesfo)

## 'data.frame':    125 obs. of  4 variables:
##  $ est : Factor w/ 5 levels "vegetativo","botão floral",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ des : num  0 0 0 0 0 0.25 0.25 0.25 0.25 0.25 ...
##  $ rept: int  1 2 3 4 5 1 2 3 4 5 ...
##  $ ncap: int  10 9 8 8 10 11 9 10 10 10 ...

## help(capdesfo)

Experimento conduzido sob delineamento inteiramente casualizado em casa de vegetação onde avaliou-se o número de capulhos produzidos por plantas de algodão submetidas à níveis de desfolha artificial de remoção foliar em combinação com o estágio fenológico no qual a desfolha foi aplicada.

est: Estágio fenológico com cinco níveis (vegetativo, botão floral, florecimento, maça, capulho);
des: Nível de desfolha artificial de remoção foliar (0, 25, 50, 75, 100%);
ncap: Número de capulhos de algodão produzidos ao final da ciclo cultura.

Análise Exploratória

## Experimento balanceado
xtabs(~est + des, data = capdesfo)

##               des
## est            0 0.25 0.5 0.75 1
##   vegetativo   5    5   5    5 5
##   botão floral 5    5   5    5 5
##   florecimento 5    5   5    5 5
##   maça         5    5   5    5 5
##   capulho      5    5   5    5 5

(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))

## Avaliando preliminarmente suposição de equidispersão
(mv <- aggregate(ncap ~ est + des, data = capdesfo,
                 FUN = function(x) c(mean = mean(x), var = var(x))))

##             est  des ncap.mean ncap.var
## 1    vegetativo 0.00       9.0      1.0
## 2  botão floral 0.00       8.4      1.3
## 3  florecimento 0.00       9.4      2.8
## 4          maça 0.00       8.6      1.8
## 5       capulho 0.00      10.0      4.0
## 6    vegetativo 0.25      10.0      0.5
## 7  botão floral 0.25       9.4      3.3
## 8  florecimento 0.25       5.8      1.2
## 9          maça 0.25       7.0      2.0
## 10      capulho 0.25      10.0      3.5
## 11   vegetativo 0.50       8.6      0.8
## 12 botão floral 0.50       8.8      0.7
## 13 florecimento 0.50       5.8      1.7
## 14         maça 0.50       8.8      2.2
## 15      capulho 0.50       8.2      1.2
## 16   vegetativo 0.75       8.0      1.0
## 17 botão floral 0.75       8.8      2.7
## 18 florecimento 0.75       6.0      2.5
## 19         maça 0.75       6.2      0.2
## 20      capulho 0.75       8.8      1.2
## 21   vegetativo 1.00       6.2      1.2
## 22 botão floral 1.00       7.2      0.2
## 23 florecimento 1.00       4.6      0.3
## 24         maça 1.00       3.0      0.5
## 25      capulho 1.00       9.0      1.0

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)
       })

Ajuste dos modelos

## 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))

## 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)

## 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)

Comparação dos ajustes

## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(list(m1P, m2P, m3P, m4P), logLik),
      "COM-Poisson" = sapply(list(m1C, m2C, m3C, m4C), logLik))

##        Poisson COM-Poisson
## [1,] -279.9328   -272.4794
## [2,] -271.3543   -256.0893
## [3,] -258.6741   -220.1977
## [4,] -255.8031   -208.2500

## Teste de razão de verossimilhanças
anova(m1P, m2P, m3P, m4P, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: ncap ~ 1
## Model 2: ncap ~ des + I(des^2)
## Model 3: ncap ~ est:des + I(des^2)
## Model 4: ncap ~ est:(des + I(des^2))
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       124     75.514                          
## 2       122     58.357  2   17.157 0.0001881 ***
## 3       118     32.997  4   25.360 4.258e-05 ***
## 4       114     27.255  4    5.742 0.2192611    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1C, m2C, m3C, m4C)

## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)
## Model 2: m2C, [llcmp]: phi+(Intercept)+des+I(des^2)
## Model 3: m3C, [llcmp]: phi+(Intercept)+I(des^2)+
##           estvegetativo:des+estbotão floral:des+
##           estflorecimento:des+estmaça:des+estcapulho:des
## Model 4: m4C, [llcmp]: phi+(Intercept)+estvegetativo:des+
##           estbotão floral:des+estflorecimento:des+estmaça:des+
##           estcapulho:des+estvegetativo:I(des^2)+estbotão
##           floral:I(des^2)+estflorecimento:I(des^2)+
##           estmaça:I(des^2)+estcapulho:I(des^2)
##   Tot Df Deviance  Chisq Df Pr(>Chisq)    
## 1      2   544.96                         
## 2      4   512.18 32.780  2  7.619e-08 ***
## 3      8   440.40 71.783  4  9.537e-15 ***
## 4     12   416.50 23.895  4  8.382e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Estimativas dos parâmetros
summary(m4P)

## 
## Call:
## glm(formula = f4, family = poisson, data = capdesfo)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.06976  -0.31729  -0.02347   0.34779   1.27069  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.189560   0.063333  34.572   <2e-16 ***
## estvegetativo:des         0.436859   0.515560   0.847   0.3968    
## estbotão floral:des       0.289715   0.507751   0.571   0.5683    
## estflorecimento:des      -1.242481   0.603706  -2.058   0.0396 *  
## estmaça:des               0.364871   0.565803   0.645   0.5190    
## estcapulho:des            0.008949   0.503760   0.018   0.9858    
## estvegetativo:I(des^2)   -0.805215   0.583915  -1.379   0.1679    
## estbotão floral:I(des^2) -0.487850   0.566394  -0.861   0.3891    
## estflorecimento:I(des^2)  0.672837   0.680195   0.989   0.3226    
## estmaça:I(des^2)         -1.310346   0.672780  -1.948   0.0515 .  
## estcapulho:I(des^2)      -0.019970   0.552505  -0.036   0.9712    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 75.514  on 124  degrees of freedom
## Residual deviance: 27.255  on 114  degrees of freedom
## AIC: 533.61
## 
## Number of Fisher Scoring iterations: 4

summary(m4C)

## Maximum likelihood estimation
## 
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, 
##     X = X, sumto = sumto), vecpar = TRUE)
## 
## Coefficients:
##                           Estimate Std. Error z value     Pr(z)    
## phi                       1.584701   0.127619 12.4174 < 2.2e-16 ***
## (Intercept)              10.896190   1.404117  7.7602 8.481e-15 ***
## estvegetativo:des         2.018067   1.140738  1.7691 0.0768791 .  
## estbotão floral:des       1.342978   1.109141  1.2108 0.2259619    
## estflorecimento:des      -5.749747   1.479754 -3.8856 0.0001021 ***
## estmaça:des               1.595068   1.229222  1.2976 0.1944168    
## estcapulho:des            0.037665   1.087991  0.0346 0.9723840    
## estvegetativo:I(des^2)   -3.723683   1.341879 -2.7750 0.0055206 ** 
## estbotão floral:I(des^2) -2.264620   1.254617 -1.8050 0.0710702 .  
## estflorecimento:I(des^2)  3.133964   1.504322  2.0833 0.0372233 *  
## estmaça:I(des^2)         -5.894220   1.611865 -3.6568 0.0002554 ***
## estcapulho:I(des^2)      -0.090159   1.193309 -0.0756 0.9397741    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 416.4999

Avaliando modelo proposto

## 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)

## Reajustando o modelo
logLik(m4C <- cmp(f4, data = capdesfo, sumto = 30))

## 'log Lik.' -208.25 (df=12)

## 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)

##-------------------------------------------
## Testando a nulidade do parâmetro phi

## Usando o ajuste Poisson
trv <- 2 * (logLik(m4C) - logLik(m4P))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)

## [1] 95.10632  0.00000

## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m4Cfixed <- cmp(f4, data = capdesfo, fixed = list("phi" = 0))
anova(m4C, m4Cfixed)

## Likelihood Ratio Tests
## Model 1: m4C, [llcmp]: phi+(Intercept)+estvegetativo:des+
##           estbotão floral:des+estflorecimento:des+estmaça:des+
##           estcapulho:des+estvegetativo:I(des^2)+estbotão
##           floral:I(des^2)+estflorecimento:I(des^2)+
##           estmaça:I(des^2)+estcapulho:I(des^2)
## Model 2: m4Cfixed, [llcmp]: (Intercept)+estvegetativo:des+
##           estbotão floral:des+estflorecimento:des+estmaça:des+
##           estcapulho:des+estvegetativo:I(des^2)+estbotão
##           floral:I(des^2)+estflorecimento:I(des^2)+
##           estmaça:I(des^2)+estcapulho:I(des^2)
##   Tot Df Deviance  Chisq Df Pr(>Chisq)    
## 1     12    416.5                         
## 2     11    511.6 95.096  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Via perfil de log-verossimilhança
perf <- profile(m4C, which = "phi")
confint(perf)

##    2.5 %   97.5 % 
## 1.323327 1.824517

plot(perf)

##-------------------------------------------
## Verificando a matriz ve variâncias e covariâncias
Vcov <- vcov(m4C)
Corr <- cov2cor(Vcov)

library(corrplot)
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)])

Predição

## Predição pontual/intervalar
pred <- with(capdesfo,
             expand.grid(
                 est = levels(est),
                 des = seq(min(des), max(des), l = 20)
             ))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)

##-------------------------------------------
## Considerando a Poisson
aux <- predict(m4P, 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]

## ncap ~ est:(des + I(des^2))

## ~est:(des + I(des^2))

X <- model.matrix(f4[-2], data = pred)

aux <- predict(m4C, newdata = X, type = "response",
               interval = "confidence")
aux <- data.frame(modelo = "COM-Poisson", aux)
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,
                    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))

Reparametrização para a média

Definindo a log-verossimilhança

## 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)
}

Definindo a função de ajuste

## 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 reparametrizado para a média

## Ajustando o modelo
m4C2 <- cmp3(f4, data = capdesfo)
convergencez(m4C2)

## Perfis de verossimilhança para phi
perf <- profile(m4C2, which = "phi")
confint(perf)

##    2.5 %   97.5 % 
## 1.320391 1.821609

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]

## ncap ~ est:(des + I(des^2))

## ~est:(des + I(des^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))

Ocorrência de ninfas de mosca-branca em variedades de soja

data(ninfas)
str(ninfas)

## 'data.frame':    240 obs. of  8 variables:
##  $ data : Date, format: "2009-12-11" ...
##  $ dias : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ cult : Factor w/ 10 levels "BRS 239","BRS 243 RR",..: 3 3 3 3 2 2 2 2 4 4 ...
##  $ bloco: Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ nsup : int  168 30 28 6 31 4 3 48 100 54 ...
##  $ nmed : int  71 46 36 1 22 6 12 76 50 15 ...
##  $ ninf : int  39 22 7 31 28 10 27 13 34 7 ...
##  $ ntot : int  278 98 71 38 81 20 42 137 184 76 ...

## help(ninfas)

Experimento conduzido em casa de vegetação sob o delineamento de blocos casualizados. No experimento foram avaliadas plantas de diferentes cultivares de soja contabilizando o número de ninfas de mosca-branca nos folíolos dos terços superior, médio e inferior das plantas. As avaliações ocorreram em 6 datas dentre os 38 dias do estudo.

Nesta análise serão consideradas somente as cultivares com prefixo , sendo o número total de ninfas de mosca-branca nos folíolos a variável de interesse.

## Somente as cultivares que contém BRS na identificação
ninfas <- droplevels(subset(ninfas, grepl("BRS", x = cult)))

## Categorizando a variável dias em aval
ninfas$aval <- factor(ninfas$dias)

str(ninfas)

## 'data.frame':    96 obs. of  9 variables:
##  $ data : Date, format: "2009-12-11" ...
##  $ dias : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ cult : Factor w/ 4 levels "BRS 239","BRS 243 RR",..: 3 3 3 3 2 2 2 2 4 4 ...
##  $ bloco: Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ nsup : int  168 30 28 6 31 4 3 48 100 54 ...
##  $ nmed : int  71 46 36 1 22 6 12 76 50 15 ...
##  $ ninf : int  39 22 7 31 28 10 27 13 34 7 ...
##  $ ntot : int  278 98 71 38 81 20 42 137 184 76 ...
##  $ aval : Factor w/ 6 levels "0","8","13","22",..: 1 1 1 1 1 1 1 1 1 1 ...

Assim as variáveis consideradas são definidas como:

bloco: Fator com 4 níveis que representam os blocos utilizados para controle de variação local.
cult: Fator com a identificação da cultivar de soja. Foram 4 cultivares (com prefixo BRS).
aval: Fator que indica o número de dias, após o início do experimento em que realizou-se as avaliações. Houve avaliações em 6 datas.

Análise Exploratória

## Experimento balanceado
xtabs(~aval + cult, data = ninfas)

##     cult
## aval BRS 239 BRS 243 RR BRS 245 RR BRS 246 RR
##   0        4          4          4          4
##   8        4          4          4          4
##   13       4          4          4          4
##   22       4          4          4          4
##   31       4          4          4          4
##   38       4          4          4          4

(xy <- xyplot(ntot ~ aval | cult,
              data = ninfas,
              type = c("p", "g", "smooth"),
              jitter.x = TRUE))

## Avaliando preliminarmente suposição de equidispersão
(mv <- aggregate(ntot ~ data + cult, data = ninfas,
                 FUN = function(x) c(mean = mean(x), var = var(x))))

##          data       cult   ntot.mean    ntot.var
## 1  2009-12-11    BRS 239   111.50000   475.66667
## 2  2009-12-19    BRS 239   149.00000  5262.66667
## 3  2009-12-24    BRS 239   170.25000  2889.58333
## 4  2010-01-02    BRS 239    54.75000  1634.91667
## 5  2010-01-11    BRS 239     6.25000    28.91667
## 6  2010-01-18    BRS 239     4.75000    40.91667
## 7  2009-12-11 BRS 243 RR    70.00000  2631.33333
## 8  2009-12-19 BRS 243 RR    76.25000  8435.58333
## 9  2009-12-24 BRS 243 RR    62.50000  2367.00000
## 10 2010-01-02 BRS 243 RR    30.00000   680.66667
## 11 2010-01-11 BRS 243 RR     4.50000     9.00000
## 12 2010-01-18 BRS 243 RR     3.50000    40.33333
## 13 2009-12-11 BRS 245 RR   121.25000 11522.25000
## 14 2009-12-19 BRS 245 RR   127.00000  9622.00000
## 15 2009-12-24 BRS 245 RR   113.50000  7052.33333
## 16 2010-01-02 BRS 245 RR    32.00000   568.66667
## 17 2010-01-11 BRS 245 RR     9.50000    51.00000
## 18 2010-01-18 BRS 245 RR     9.00000   134.66667
## 19 2009-12-11 BRS 246 RR    86.75000  4691.58333
## 20 2009-12-19 BRS 246 RR   124.00000  4790.66667
## 21 2009-12-24 BRS 246 RR   133.25000  5726.91667
## 22 2010-01-02 BRS 246 RR    38.75000   648.91667
## 23 2010-01-11 BRS 246 RR     7.00000    72.66667
## 24 2010-01-18 BRS 246 RR     4.50000    69.66667

xlim <- ylim <- extendrange(c(mv$ntot), f = 0.05)
xyplot(ntot[, "var"] ~ ntot[, "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)
       })

Ajuste dos modelos

## Preditores considerados
f1 <- ntot ~ bloco + cult + aval
f2 <- ntot ~ bloco + cult * aval

## Ajustando os modelos Poisson
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 = 800)
m2C <- cmp(f2, data = ninfas, sumto = 800)

Comparação dos ajustes

## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(list(m1P, m2P), logLik),
      "COM-Poisson" = sapply(list(m1C, m2C), logLik))

##        Poisson COM-Poisson
## [1,] -922.9782   -410.4441
## [2,] -879.2287   -407.1332

## Teste de razão de verossimilhanças
anova(m1P, m2P, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: ntot ~ bloco + cult + aval
## Model 2: ntot ~ bloco + cult * aval
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        84     1371.3                          
## 2        69     1283.8 15   87.499 2.897e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1C, m2C)

## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)+blocoII+blocoIII+
##           blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
##           aval8+aval13+aval22+aval31+aval38
## Model 2: m2C, [llcmp]: phi+(Intercept)+blocoII+blocoIII+
##           blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
##           aval8+aval13+aval22+aval31+aval38+cultBRS 243
##           RR:aval8+cultBRS 245 RR:aval8+cultBRS 246 RR:aval8+
##           cultBRS 243 RR:aval13+cultBRS 245 RR:aval13+cultBRS 246
##           RR:aval13+cultBRS 243 RR:aval22+cultBRS 245 RR:aval22+
##           cultBRS 246 RR:aval22+cultBRS 243 RR:aval31+cultBRS 245
##           RR:aval31+cultBRS 246 RR:aval31+cultBRS 243 RR:aval38+
##           cultBRS 245 RR:aval38+cultBRS 246 RR:aval38
##   Tot Df Deviance  Chisq Df Pr(>Chisq)
## 1     13   820.89                     
## 2     28   814.27 6.6218 15     0.9673

## Estimativas dos parâmetros
summary(m1P)

## 
## Call:
## glm(formula = f1, family = poisson, data = ninfas)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -13.2911   -3.0119   -0.7386    1.7743   12.3560  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     5.33397    0.03480 153.262  < 2e-16 ***
## blocoII        -0.52019    0.03228 -16.114  < 2e-16 ***
## blocoIII       -0.81268    0.03555 -22.857  < 2e-16 ***
## blocoIV        -0.99360    0.03792 -26.204  < 2e-16 ***
## cultBRS 243 RR -0.69921    0.03894 -17.954  < 2e-16 ***
## cultBRS 245 RR -0.18595    0.03332  -5.582 2.38e-08 ***
## cultBRS 246 RR -0.23060    0.03373  -6.837 8.10e-12 ***
## aval8           0.20108    0.03416   5.887 3.94e-09 ***
## aval13          0.20788    0.03411   6.095 1.09e-09 ***
## aval22         -0.91822    0.04743 -19.360  < 2e-16 ***
## aval31         -2.65981    0.09908 -26.846  < 2e-16 ***
## aval38         -2.88525    0.11016 -26.191  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 7106.3  on 95  degrees of freedom
## Residual deviance: 1371.3  on 84  degrees of freedom
## AIC: 1870
## 
## Number of Fisher Scoring iterations: 5

summary(m1C)

## Maximum likelihood estimation
## 
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y, 
##     X = X, sumto = sumto), vecpar = TRUE)
## 
## Coefficients:
##                  Estimate Std. Error  z value     Pr(z)    
## phi            -3.0828771  0.2122292 -14.5262 < 2.2e-16 ***
## (Intercept)     0.2434861  0.0526161   4.6276 3.699e-06 ***
## blocoII        -0.0268210  0.0089188  -3.0072 0.0026364 ** 
## blocoIII       -0.0428707  0.0113906  -3.7637 0.0001674 ***
## blocoIV        -0.0533302  0.0130955  -4.0724 4.653e-05 ***
## cultBRS 243 RR -0.0380972  0.0113820  -3.3472 0.0008165 ***
## cultBRS 245 RR -0.0096905  0.0078164  -1.2398 0.2150606    
## cultBRS 246 RR -0.0120552  0.0080308  -1.5011 0.1333264    
## aval8           0.0103003  0.0079612   1.2938 0.1957316    
## aval13          0.0106428  0.0079638   1.3364 0.1814228    
## aval22         -0.0523132  0.0143943  -3.6343 0.0002787 ***
## aval31         -0.2305570  0.0454114  -5.0771 3.833e-07 ***
## aval38         -0.2708058  0.0532561  -5.0850 3.677e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 820.8882

Avaliando modelo proposto

## 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(m1C)

## 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(m1P)

##-------------------------------------------
## Testando a nulidade do parâmetro phi

## Usando o ajuste Poisson
trv <- 2 * (logLik(m1C) - logLik(m1P))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)

## [1] 1025.068    0.000

## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m1Cfixed <- cmp(f1, data = ninfas, fixed = list("phi" = 0))
anova(m1C, m1Cfixed)

## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)+blocoII+blocoIII+
##           blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
##           aval8+aval13+aval22+aval31+aval38
## Model 2: m1Cfixed, [llcmp]: (Intercept)+blocoII+blocoIII+
##           blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
##           aval8+aval13+aval22+aval31+aval38
##   Tot Df Deviance  Chisq Df Pr(>Chisq)    
## 1     13   820.89                         
## 2     12  1845.96 1025.1  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Via perfil de log-verossimilhança
perf <- profile(m1C, which = "phi")
confint(perf)

##     2.5 %    97.5 % 
## -3.550089 -2.700390

plot(perf)

##-------------------------------------------
## Verificando a matriz ve variâncias e covariâncias
Vcov <- vcov(m1C)
Corr <- cov2cor(Vcov)

library(corrplot)
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)])

Predição

## Predição pontual/intervalar
pred <- with(ninfas,
             expand.grid(
                 bloco = factor(levels(bloco)[1],
                                levels = levels(bloco)),
                 cult = levels(cult),
                 aval = levels(aval)
             ))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)

f1; f1[-2]

## ntot ~ bloco + cult + aval

## ~bloco + cult + aval

X <- model.matrix(f1[-2], data = pred)

## Como não temos interesse na interpretação dos efeitos de blocos
## tomaremos a média desses efeitos para predição

bl <- attr(X, "assign") == 1
X[, bl] <- X[, bl] * 0 + 1/(sum(bl) + 1)
head(X)

##   (Intercept) blocoII blocoIII blocoIV cultBRS 243 RR
## 1           1    0.25     0.25    0.25              0
## 2           1    0.25     0.25    0.25              1
## 3           1    0.25     0.25    0.25              0
## 4           1    0.25     0.25    0.25              0
## 5           1    0.25     0.25    0.25              0
## 6           1    0.25     0.25    0.25              1
##   cultBRS 245 RR cultBRS 246 RR aval8 aval13 aval22 aval31 aval38
## 1              0              0     0      0      0      0      0
## 2              0              0     0      0      0      0      0
## 3              1              0     0      0      0      0      0
## 4              0              1     0      0      0      0      0
## 5              0              0     1      0      0      0      0
## 6              0              0     1      0      0      0      0

library(multcomp)

##-------------------------------------------
## Considerando a Poisson
aux <- exp(confint(glht(m1P, linfct = X),
               calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Poisson", aux)
predP <- cbind(pred, aux)

##-------------------------------------------
## Considerando a COM-Poisson
aux <- predict(m1C, newdata = X, type = "response",
               interval = "confidence")
aux <- data.frame(modelo = "COM-Poisson", aux)
predC <- cbind(pred, aux)

##-------------------------------------------
## Visualizando os valores preditos intervalares pelos dois modelos
da <- rbind(predP, predC)
da <- da[order(da$cult, da$aval, da$modelo), ]

source(paste0("https://gitlab.c3sl.ufpr.br/leg/legTools/raw/",
              "issue%2315/R/panel.segplot.by.R"))

key <- list(type = "o", divide = 1,
            lines = list(pch = 1:nlevels(da$modelo) + 3, lty = 1),
            text = list(c("Poisson", "COM-Poisson")))

update(xy, type = c("p", "g"), key = key, alpha = 0.7) +
    as.layer(segplot(
        aval ~ lwr + upr | cult,
        centers = fit, groups = modelo, data = da,
        grid = TRUE, horizontal = FALSE, draw = FALSE,
        lwd = 2, pch = 1:nlevels(da$modelo) + 3,
        panel = panel.segplot.by, f = 0.1))

Modelos de Regressão para Dados de Contagem com o R

`github.com/leg-ufpr/MRDCr`

Análise de Contagens com o Modelo COM-Poisson

Walmes M. Zeviani, Eduardo E. Ribeiro Jr & Cesar A. Taconeli

Funções para ajuste dos modelos

Log-verossimilhança

Ajuste geral

Capulhos de algodão sob exposição à mosca-branca

Análise Exploratória

Ajuste dos modelos

Comparação dos ajustes

Avaliando modelo proposto

Predição

Produção de soja sob efeito de umidade e adubação potássica

Análise Exploratória

Ajuste dos modelos

Comparação dos ajustes

Avaliando modelo proposto

Predição

Capulhos de algodão sob efeito de desfolha

Análise Exploratória

Ajuste dos modelos

Comparação dos ajustes

Avaliando modelo proposto

Predição

Reparametrização para a média

Definindo a log-verossimilhança

Definindo a função de ajuste

Ajustando o modelo reparametrizado para a média

Ocorrência de ninfas de mosca-branca em variedades de soja

Análise Exploratória

Ajuste dos modelos

Comparação dos ajustes

Avaliando modelo proposto

Predição

Modelos de Regressão para Dados de Contagem com o R