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

github.com/leg-ufpr/MRDCr
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)\)

llcmp
## function (params, y, X, offset = NULL, sumto = ceiling(max(y)^1.5)) 
## {
##     betas <- params[-1]
##     phi <- params[1]
##     nu <- exp(phi)
##     if (is.null(offset)) 
##         offset <- 0
##     Xb <- X %*% betas + offset
##     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>

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\).

  • Inclusão do termo offset. Somado diretamente ao preditor \(X_i \beta\), pois \(X_i \beta\) representa o parâmetro \(\lambda\) de locação, da distribuição COM-Poisson.

  • Truncamento da série infinita \(Z(\lambda_i)\). sumto é tomado como argumento da função, que por padrão assume \(\max(\underline{y})^{1.2}\).

  • 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!)) \]

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, ...) 
## {
##     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, offset = off, 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, ...)
##     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
## Definindo o prditor do modelo
formula <- y ~ x + I(x^2) + offset(log(t))

##-------------------------------------------
## 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
## 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
## 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
library(bbmle)
## Loading required package: stats4
parnames(llcmp) <- names(start)
mle2(llcmp, start = start,
     data = list(y = y, X = X, offset = o, sumto = 50),
     vecpar = TRUE)
## 
## Call:
## mle2(minuslogl = llcmp, start = start, data = list(y = y, X = X, 
##     offset = o, sumto = 50), vecpar = TRUE)
## 
## Coefficients:
##         phi (Intercept)           x      I(x^2) 
##   0.3338908  -4.5188932   5.4556688  -1.1173870 
## 
## Log-likelihood: -11.13

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.

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

## 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, offset = off, sumto = sumto), vecpar = TRUE)
## 
## Coefficients:
##                           Estimate Std. Error z value     Pr(z)    
## phi                       1.584744   0.127631 12.4166 < 2.2e-16 ***
## (Intercept)              10.896662   1.404311  7.7594 8.531e-15 ***
## estvegetativo:des         2.018744   1.140780  1.7696 0.0767909 .  
## estbotão floral:des       1.343071   1.109168  1.2109 0.2259408    
## estflorecimento:des      -5.750451   1.479861 -3.8858 0.0001020 ***
## estmaça:des               1.594963   1.229254  1.2975 0.1944575    
## estcapulho:des            0.037655   1.088012  0.0346 0.9723918    
## estvegetativo:I(des^2)   -3.724525   1.341954 -2.7754 0.0055125 ** 
## estbotão floral:I(des^2) -2.264721   1.254650 -1.8051 0.0710650 .  
## estflorecimento:I(des^2)  3.134664   1.504387  2.0837 0.0371891 *  
## estmaça:I(des^2)         -5.894333   1.611937 -3.6567 0.0002555 ***
## estcapulho:I(des^2)      -0.090110   1.193331 -0.0755 0.9398080    
## ---
## 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)

## 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.50                         
## 2     11   511.61 95.106  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 = 1)
confint(perf)
##    2.5 %   97.5 % 
## 1.323326 1.824516
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
pred <- with(capdesfo,
             expand.grid(
                 est = levels(est),
                 des = seq(min(des), max(des), l = 20)
             ))

##-------------------------------------------
## Considerando a Poisson
mediasP <- exp(predict(m4P, newdata = pred))
aux <- data.frame(modelo = "Poisson", fit = mediasP)
predP <- cbind(pred, aux)

##-------------------------------------------
## Considerando a COM-Poisson
f4; f4[[2]] <- NULL; f4
## ncap ~ est:(des + I(des^2))
## ~est:(des + I(des^2))
X <- model.matrix(f4, data = pred)

## Obtendo os parâmetros da distribuição (lambdas e phi)
betas <- coef(m4C)[-1]
phi <- coef(m4C)[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)
predC <- cbind(pred, aux)

##-------------------------------------------
## Visualizando os valores preditos pelos dois modelos
da <- rbind(predP, predC)

update(xy, type = c("p", "g")) +
    as.layer(xyplot(fit ~ des | est,
                    groups = modelo,
                    data = da, type = "l"))

## Predição intervalar
qn <- qnorm(0.975) * c(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 = "*")))
predP <- cbind(predP, aux)

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

## 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:50; py <- dcmp(y, p, phi, sumto = 100)
                     sum(y*py)
                 })
             })
predC <- cbind(predC, aux)

##-------------------------------------------
## Visualizando os valores preditos intervalares pelos dois modelos
da <- rbind(predP, predC)

update(xy, type = c("p", "g")) +
    as.layer(xyplot(fit ~ des | est,
                    groups = modelo,
                    data = da,
                    type = "l",
                    ly = da$lwr,
                    uy = da$upr,
                    cty = "bands",
                    alpha = 0.3,
                    prepanel = prepanel.cbH,
                    panel.groups = panel.cbH,
                    panel = panel.superpose))

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:

Análise Exploratória

## Experimento balanceado
xtabs(~dexp, data = capmosca)
## dexp
## 0 1 2 3 4 5 
## 5 5 5 5 5 5
(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, offset = off, 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 .  
## ---
## 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)

## 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.87 12.887  1  0.0003309 ***
## ---
## 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]] <- NULL; f3
## ncap ~ dexp + I(dexp^2)
## ~dexp + I(dexp^2)
X <- model.matrix(f3, data = pred)

## Obtendo os parâmetros da distribuição (lambdas e phi)
betas <- coef(m3C)[-1]
phi <- coef(m3C)[1]
loglambdas <- 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 <- 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 intervalares pelos dois modelos
da <- rbind(predP, predC)

update(xy, type = c("p", "g")) +
    as.layer(xyplot(fit ~ dexp,
                    groups = modelo,
                    data = da,
                    type = "l",
                    ly = da$lwr,
                    uy = da$upr,
                    cty = "bands",
                    alpha = 0.3,
                    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:

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 = 600)
m2C <- cmp(f2, data = ninfas, sumto = 600)
## Warning in bbmle::mle2(llcmp, start = start, data = list(y = y, X
## = X, offset = off, : convergence failure: code=1 (iteration limit
## 'maxit' reached)

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.4440
## [2,] -879.2287   -407.1331
## 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.6219 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, offset = off, sumto = sumto), vecpar = TRUE)
## 
## Coefficients:
##                  Estimate Std. Error  z value     Pr(z)    
## phi            -3.0829138  0.2092160 -14.7356 < 2.2e-16 ***
## (Intercept)     0.2434784  0.0518954   4.6917 2.709e-06 ***
## blocoII        -0.0268210  0.0088711  -3.0234 0.0024994 ** 
## blocoIII       -0.0428702  0.0113061  -3.7918 0.0001496 ***
## blocoIV        -0.0533294  0.0129892  -4.1057 4.031e-05 ***
## cultBRS 243 RR -0.0380971  0.0113435  -3.3585 0.0007837 ***
## cultBRS 245 RR -0.0096912  0.0078249  -1.2385 0.2155260    
## cultBRS 246 RR -0.0120559  0.0080379  -1.4999 0.1336482    
## aval8           0.0103005  0.0079581   1.2943 0.1955470    
## aval13          0.0106430  0.0079610   1.3369 0.1812566    
## aval22         -0.0523118  0.0143159  -3.6541 0.0002581 ***
## aval31         -0.2305532  0.0452463  -5.0955 3.478e-07 ***
## aval38         -0.2708014  0.0530935  -5.1005 3.388e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 820.888

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, incremento = 100, maxit = 10)

## 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 = 1)
confint(perf)
##     2.5 %    97.5 % 
## -3.550614 -2.700415
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]] <- NULL; f1
## ntot ~ bloco + cult + aval
## ~bloco + cult + aval
X <- model.matrix(f1, 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

## Obtendo os parâmetros da distribuição (lambdas e phi)
betas <- coef(m1C)[-1]
phi <- coef(m1C)[1]
loglambdas <- 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 <- 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)
                     sum(y*py)
                 })
             })
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"))

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