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

github.com/leg-ufpr/MRDCr
# Definições da sessão.
library(lattice)
library(latticeExtra)
library(grid)
library(plyr)
library(car)
library(doBy)
library(multcomp)
library(MRDCr)

Acomodando superdispersão com efeito aleatório

data(nematoide, package = "MRDCr")
str(nematoide)
## 'data.frame':    94 obs. of  5 variables:
##  $ cult: Factor w/ 19 levels "A","B","C","D",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ mfr : num  10.52 11.03 6.42 8.16 4.48 ...
##  $ vol : int  40 40 40 40 40 40 40 40 40 40 ...
##  $ nema: int  4 5 3 4 3 2 2 2 2 2 ...
##  $ off : num  0.263 0.276 0.161 0.204 0.112 ...
m0 <- glm(nema ~ offset(log(off)) + cult,
          data = nematoide, family = poisson)
m1 <- update(m0, family = quasipoisson)

library(lme4)
## Loading required package: Matrix
nematoide$ue <- 1:nrow(nematoide)

m2 <- glmer(nema ~ offset(log(off)) + (1 | cult) + (1 | ue),
            data = nematoide, family = poisson)
summary(m2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: nema ~ offset(log(off)) + (1 | cult) + (1 | ue)
##    Data: nematoide
## 
##      AIC      BIC   logLik deviance df.resid 
##    561.3    568.9   -277.6    555.3       91 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.38276 -0.28758 -0.00125  0.32557  2.16491 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ue     (Intercept) 0.2583   0.5082  
##  cult   (Intercept) 0.4259   0.6526  
## Number of obs: 94, groups:  ue, 94; cult, 19
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.1050     0.1665   24.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Número de Grãos em Soja

data(soja, package = "MRDCr")
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 ...
soja <- soja[-74, ]
soja$K <- factor(soja$K)
soja$ue <- 1:nrow(soja)
str(soja)
## 'data.frame':    74 obs. of  6 variables:
##  $ K   : Factor w/ 5 levels "0","30","60",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ 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 ...
##  $ ue  : int  1 2 3 4 5 6 7 8 9 10 ...
m0 <- glm(ngra ~ bloc + umid * K,
          data = soja, family = poisson)
m1 <- update(m0, family = quasipoisson)

m2 <- glmer(ngra ~ (1 | bloc) + umid * K,
            data = soja, family = poisson)
m3 <- update(m2, . ~ . + (1 | ue))


logLik(m0)
## 'log Lik.' -321.6692 (df=19)
logLik(m2)
## 'log Lik.' -327.9654 (df=16)
logLik(m3)
## 'log Lik.' -320.1019 (df=17)
anova(m2, m3)
## Data: soja
## Models:
## m2: ngra ~ (1 | bloc) + umid * K
## m3: ngra ~ (1 | bloc) + umid + K + (1 | ue) + umid:K
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m2 16 687.93 724.80 -327.97   655.93                             
## m3 17 674.20 713.37 -320.10   640.20 15.727      1  7.317e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: ngra ~ (1 | bloc) + umid + K + (1 | ue) + umid:K
##    Data: soja
## 
##      AIC      BIC   logLik deviance df.resid 
##    674.2    713.4   -320.1    640.2       57 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.52165 -0.53703 -0.00107  0.39545  1.82190 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.
##  ue     (Intercept) 0.0045576 0.06751 
##  bloc   (Intercept) 0.0009459 0.03076 
## Number of obs: 74, groups:  ue, 74; bloc, 5
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    4.74880    0.05319   89.27  < 2e-16 ***
## umid50         0.13376    0.07116    1.88  0.06014 .  
## umid62,5       0.18695    0.07061    2.65  0.00811 ** 
## K30            0.29857    0.06953    4.29 1.75e-05 ***
## K60            0.34414    0.06911    4.98 6.38e-07 ***
## K120           0.36651    0.06892    5.32 1.05e-07 ***
## K180           0.29494    0.06956    4.24 2.23e-05 ***
## umid50:K30     0.03898    0.09618    0.41  0.68525    
## umid62,5:K30  -0.13766    0.09650   -1.43  0.15373    
## umid50:K60     0.11420    0.09523    1.20  0.23048    
## umid62,5:K60   0.09077    0.09467    0.96  0.33766    
## umid50:K120    0.11716    0.09497    1.23  0.21731    
## umid62,5:K120  0.16415    0.09654    1.70  0.08906 .  
## umid50:K180    0.28731    0.09497    3.03  0.00248 ** 
## umid62,5:K180  0.21509    0.09464    2.27  0.02305 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.00247726 (tol = 0.001, component 1)
#-----------------------------------------------------------------------
# Predição com intervalos de confiança.

X <- LSmatrix(m0, effect = c("umid", "K"))
pred <- attr(X, "grid")
pred <- transform(pred,
                  K = as.integer(as.character(K)),
                  umid = factor(umid))
pred <- list(pois = pred, quasi = pred, pmis1 = pred, pmis2 = pred)

# Quantil normal para fazer um IC de 95%.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)

# Preditos pela Poisson.
aux <- confint(glht(m0, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))

aux <- confint(glht(m1, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))

# Removendo as colunas correspondentes ao blocos.
X <- X[, -grep(pattern = "^bloc", x = colnames(X))]

# Poisson Misto 1: ~ (1 | bloc)
aux <- confint(glht(m2, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pmis1 <- cbind(pred$pmis1, exp(aux))

# Poisson Misto 2: ~ (1 | bloc) + (1 | ue).
aux <- confint(glht(m3, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pmis2 <- cbind(pred$pmis2, exp(aux))

pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, umid, K, modelo)

key <- list(type = "o", divide = 1,
            lines = list(pch = 1:nlevels(pred$modelo),
                         lty = 1, col = 1),
            text = list(c("Poisson", "Quasi-Poisson",
                          "Poissin Misto 1", "Poissin Misto 2")))

xyplot(fit ~ K | umid, data = pred,
       layout = c(NA, 1), as.table = TRUE,
       xlim = extendrange(range(pred$K), f = 0.1),
       key = key, pch = pred$modelo,
       xlab = expression("Dose de potássio"~(mg~dm^{-3})),
       ylab = "Número de grãos por parcela",
       ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
       prepanel = prepanel.cbH,
       desloc = 6 * scale(as.integer(pred$modelo), scale = FALSE),
       panel = panel.cbH)

Resfriamento de Cobertura em Aviários na Mortalidade das Aves

#-----------------------------------------------------------------------

data(confterm, package = "MRDCr")
data(conftemp, package = "MRDCr")

xyplot(nap ~ idade | resfr, data = confterm,
       groups = galp, type = "o",
       xlab = "Idade das aves (dias)",
       ylab = "Número de aves perdidas por galpão",
       strip = strip.custom(factor.levels = c(
                                "Com sistema de resfriamento",
                                "Sem sistema de resfriamento")),
       auto.key = list(corner = c(0.05, 0.9)))

# Amplitude estendida das variáveis.
lim <- with(conftemp, apply(cbind(h, ctr, itgu), MARGIN = 2,
                            FUN = extendrange, f = 0.2))

# Anotação da eixo x do gráfico.
scales <- list(
    y = list(relation = "free"),
    x = list(at = seq(from = 1,
                      to = ceiling(max(conftemp$hr/24)) * 24,
                      by = 24)))
scales$x$labels <- seq_along(scales$x$at)

xyplot(h + ctr + itgu ~ hr, data = conftemp,
       outer = TRUE, type = "l", layout = c(1, NA),
       scales = scales, xlim = range(scales$x$at),
       xlab = "Dias",
       ylab = "Variáveis térmicas",
       panel = function(y, subscripts, ...) {
           wp <- which.packet()
           r <- lim[, wp[1]]
           panel.rect(10.5 + 24 * (scales$x$labels - 1), r[1],
                      20 + 24 * (scales$x$labels - 1), r[2],
                      col = "blue",
                      border = "transparent",
                      alpha = 0.25)
           panel.xyplot(y = y, subscripts = subscripts, ...)
       })

#-----------------------------------------------------------------------
# Juntando os datasets.

tempdia <- aggregate(cbind(hmax = h, cmax = ctr, imax = itgu) ~ idade,
                     data = conftemp, FUN = max)
splom(tempdia)

confterm <- merge(confterm, tempdia, by = "idade")
str(confterm)
## 'data.frame':    76 obs. of  7 variables:
##  $ idade: int  21 21 21 21 22 22 22 22 23 23 ...
##  $ resfr: Factor w/ 2 levels "C","S": 1 2 1 2 1 1 2 2 1 1 ...
##  $ galp : Factor w/ 4 levels "I","II","III",..: 1 4 2 3 2 1 4 3 1 2 ...
##  $ nap  : int  5 14 9 7 7 12 6 9 4 8 ...
##  $ hmax : num  81.4 81.4 81.4 81.4 78 ...
##  $ cmax : num  730 730 730 730 645 ...
##  $ imax : num  81.8 81.8 81.8 81.8 79.1 ...
summary(confterm)
##      idade    resfr   galp         nap              hmax      
##  Min.   :21   C:38   I  :19   Min.   :  4.00   Min.   :73.38  
##  1st Qu.:25   S:38   II :19   1st Qu.:  9.00   1st Qu.:77.56  
##  Median :30          III:19   Median : 18.50   Median :78.56  
##  Mean   :30          IV :19   Mean   : 28.09   Mean   :79.64  
##  3rd Qu.:35                   3rd Qu.: 36.25   3rd Qu.:82.03  
##  Max.   :39                   Max.   :186.00   Max.   :86.18  
##       cmax            imax      
##  Min.   :616.8   Min.   :76.20  
##  1st Qu.:650.3   1st Qu.:79.07  
##  Median :685.2   Median :80.43  
##  Mean   :686.9   Mean   :81.09  
##  3rd Qu.:718.0   3rd Qu.:83.36  
##  Max.   :758.2   Max.   :86.44
# Na escala original, ao ajustar o modelo de efeitos aleatórios,
# apareceu o aviso que segere diminuir a escala dos dados. Sendo assim,
# os dados serão padronizados.
#
# Warning messages:
# 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
#   Model failed to converge with max|grad| = 0.00183911 (tol = 0.001, component 1)
# 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
#   Model is nearly unidentifiable: very large eigenvalue
#  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
#  - Rescale variables?

confterm$ue <- 1:nrow(confterm)
confterm <- within(confterm, {
    idade <- idade - min(idade)
    idade <- idade/max(idade)
    hmax <- hmax - min(hmax)
    hmax <- hmax/max(hmax)
})

summary(confterm)
##      idade        resfr   galp         nap              hmax       
##  Min.   :0.0000   C:38   I  :19   Min.   :  4.00   Min.   :0.0000  
##  1st Qu.:0.2222   S:38   II :19   1st Qu.:  9.00   1st Qu.:0.3264  
##  Median :0.5000          III:19   Median : 18.50   Median :0.4048  
##  Mean   :0.5000          IV :19   Mean   : 28.09   Mean   :0.4886  
##  3rd Qu.:0.7778                   3rd Qu.: 36.25   3rd Qu.:0.6756  
##  Max.   :1.0000                   Max.   :186.00   Max.   :1.0000  
##       cmax            imax             ue       
##  Min.   :616.8   Min.   :76.20   Min.   : 1.00  
##  1st Qu.:650.3   1st Qu.:79.07   1st Qu.:19.75  
##  Median :685.2   Median :80.43   Median :38.50  
##  Mean   :686.9   Mean   :81.09   Mean   :38.50  
##  3rd Qu.:718.0   3rd Qu.:83.36   3rd Qu.:57.25  
##  Max.   :758.2   Max.   :86.44   Max.   :76.00
#-----------------------------------------------------------------------
# Ajuste dos modelos.

m0 <- glm(nap ~ galp + resfr * (idade + hmax),
          data = confterm, family = poisson)
m1 <- update(m0, family = quasipoisson)

anova(m0, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: nap
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                           75    1822.52             
## galp         3   243.04        72    1579.48  < 2e-16 ***
## resfr        0     0.00        72    1579.48             
## idade        1  1303.94        71     275.54  < 2e-16 ***
## hmax         1     1.15        70     274.39  0.28341    
## resfr:idade  1    68.01        69     206.38  < 2e-16 ***
## resfr:hmax   1     4.93        68     201.45  0.02643 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1, test = "F")
## Analysis of Deviance Table
## 
## Model: quasipoisson, link: log
## 
## Response: nap
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev        F    Pr(>F)    
## NULL                           75    1822.52                       
## galp         3   243.04        72    1579.48  24.6204 6.742e-11 ***
## resfr        0     0.00        72    1579.48                       
## idade        1  1303.94        71     275.54 396.2752 < 2.2e-16 ***
## hmax         1     1.15        70     274.39   0.3497    0.5562    
## resfr:idade  1    68.01        69     206.38  20.6684 2.307e-05 ***
## resfr:hmax   1     4.93        68     201.45   1.4975    0.2253    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m1)
## 
## Call:
## glm(formula = nap ~ galp + resfr * (idade + hmax), family = quasipoisson, 
##     data = confterm)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.5682  -0.9809  -0.2827   0.6872   5.8073  
## 
## Coefficients: (1 not defined because of singularities)
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.6980     0.2617   6.487 1.17e-08 ***
## galpII         0.4392     0.1341   3.276  0.00166 ** 
## galpIII       -0.3894     0.3596  -1.083  0.28265    
## galpIV        -0.7186     0.3618  -1.986  0.05106 .  
## resfrS             NA         NA      NA       NA    
## idade          1.9248     0.2514   7.656 9.24e-11 ***
## hmax          -0.1425     0.2791  -0.511  0.61126    
## resfrS:idade   1.6198     0.3434   4.717 1.23e-05 ***
## resfrS:hmax    0.4441     0.3630   1.223  0.22543    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 3.290501)
## 
##     Null deviance: 1822.52  on 75  degrees of freedom
## Residual deviance:  201.45  on 68  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m2 <- glmer(nap ~ (1 | galp) + resfr * (idade + hmax),
            data = confterm, family = poisson)

summary(m2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: nap ~ (1 | galp) + resfr * (idade + hmax)
##    Data: confterm
## 
##      AIC      BIC   logLik deviance df.resid 
##    592.3    608.7   -289.2    578.3       69 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3877 -0.9287 -0.3148  0.6909  6.3847 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  galp   (Intercept) 0.03528  0.1878  
## Number of obs: 76, groups:  galp, 4
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.9196     0.1910  10.051  < 2e-16 ***
## resfrS        -0.7753     0.2685  -2.888  0.00388 ** 
## idade          1.9248     0.1385  13.895  < 2e-16 ***
## hmax          -0.1425     0.1538  -0.927  0.35400    
## resfrS:idade   1.6198     0.1892   8.560  < 2e-16 ***
## resfrS:hmax    0.4441     0.2000   2.220  0.02641 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) resfrS idade  hmax   rsfrS:d
## resfrS      -0.711                             
## idade       -0.605  0.431                      
## hmax        -0.519  0.369  0.334               
## resfrS:idad  0.443 -0.623 -0.732 -0.245        
## resfrS:hmax  0.399 -0.512 -0.257 -0.769  0.368
anova(m2)
## Analysis of Variance Table
##             Df  Sum Sq Mean Sq   F value
## resfr        1   18.67   18.67   18.6716
## idade        1 1006.17 1006.17 1006.1733
## hmax         1    0.47    0.47    0.4731
## resfr:idade  1   69.28   69.28   69.2830
## resfr:hmax   1    4.92    4.92    4.9241
m3 <- update(m2, . ~ . + (1 | ue))

anova(m2, m3)
## Data: confterm
## Models:
## m2: nap ~ (1 | galp) + resfr * (idade + hmax)
## m3: nap ~ (1 | galp) + resfr + idade + hmax + (1 | ue) + resfr:idade + 
## m3:     resfr:hmax
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m2  7 592.35 608.66 -289.17   578.35                             
## m3  8 521.37 540.02 -252.69   505.37 72.973      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## nap ~ (1 | galp) + resfr + idade + hmax + (1 | ue) + resfr:idade +  
##     resfr:hmax
##    Data: confterm
## 
##      AIC      BIC   logLik deviance df.resid 
##    521.4    540.0   -252.7    505.4       68 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.63478 -0.39487 -0.09488  0.27327  3.10964 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ue     (Intercept) 0.06158  0.2482  
##  galp   (Intercept) 0.01373  0.1172  
## Number of obs: 76, groups:  ue, 76; galp, 4
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.8857     0.2045   9.221  < 2e-16 ***
## resfrS        -0.6865     0.2989  -2.297   0.0216 *  
## idade          1.9427     0.2000   9.713  < 2e-16 ***
## hmax          -0.1247     0.2269  -0.550   0.5826    
## resfrS:idade   1.4719     0.2891   5.090 3.57e-07 ***
## resfrS:hmax    0.4140     0.3171   1.306   0.1917    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) resfrS idade  hmax   rsfrS:d
## resfrS      -0.685                             
## idade       -0.707  0.484                      
## hmax        -0.674  0.461  0.258               
## resfrS:idad  0.491 -0.737 -0.692 -0.179        
## resfrS:hmax  0.483 -0.685 -0.185 -0.716  0.304 
## convergence code: 0
## Model failed to converge with max|grad| = 0.00111972 (tol = 0.001, component 1)
anova(m3)
## Analysis of Variance Table
##             Df Sum Sq Mean Sq  F value
## resfr        1  12.42   12.42  12.4170
## idade        1 366.81  366.81 366.8053
## hmax         1   0.07    0.07   0.0728
## resfr:idade  1  24.18   24.18  24.1849
## resfr:hmax   1   1.70    1.70   1.6980
#-----------------------------------------------------------------------
# Predição com bandas de confiança.

pred <- unique(subset(confterm, select = c(idade, resfr, hmax)))
X <- model.matrix(terms(m2), data = cbind(pred, nap = 1))
pred$nap <- NULL
str(pred)
## 'data.frame':    38 obs. of  3 variables:
##  $ idade: num  0 0 0.0556 0.0556 0.1111 ...
##  $ resfr: Factor w/ 2 levels "C","S": 1 2 1 2 1 2 1 2 1 2 ...
##  $ hmax : num  0.63 0.63 0.359 0.359 0.249 ...
# pred <- list(pois = pred, quasi = pred, pmis1 = pred, pmis2 = pred)
pred <- list(pmis1 = pred, pmis2 = pred)

# Quantil normal para fazer um IC de 95%.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)

# Poisson Misto 1: ~ (1 | galp)
aux <- confint(glht(m2, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pmis1 <- cbind(pred$pmis1, exp(aux))

# Poisson Misto 2: ~ (1 | galp) + (1 | ue).
aux <- confint(glht(m3, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pmis2 <- cbind(pred$pmis2, exp(aux))

pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, idade, resfr, modelo)

key <- list(type = "o", divide = 1,
            lines = list(pch = 1:nlevels(pred$modelo),
                         lty = 1, col = 1),
            text = list(c("Poisson", "Quasi-Poisson",
                          "Poissin Misto 1", "Poissin Misto 2")))

pred$idade <- 21 + (39 - 21) * pred$idade
confterm$idade <- 21 + (39 - 21) * confterm$idade

xyplot(fit ~ idade | modelo, groups = resfr, data = pred,
       layout = c(NA, 1), as.table = TRUE, type = "l",
       auto.key = list(lines = TRUE, points = FALSE,
                       text = c("Com resfr.", "Sem resfr.")),
       xlab = "Idade das aves (dias)",
       ylab = "Número de aves perdidas",
       strip = strip.custom(
           factor.levels = c("P. Misto 1", "P. Misto 2")),
       ly = pred$lwr, uy = pred$upr, cty = "bands", alpha = 0.25,
       prepanel = prepanel.cbH,
       panel.groups = panel.cbH,
       panel = panel.superpose) +
    as.layer(xyplot(nap ~ idade, groups = resfr, data = confterm))

Informações da sessão

## Atualizado em 21 de julho de 2016.
## 
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04 LTS
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    methods   grid      stats     graphics  grDevices
## [7] utils     datasets  base     
## 
## other attached packages:
##  [1] lme4_1.1-12         Matrix_1.2-6        MRDCr_0.0-2        
##  [4] bbmle_1.0.18        multcomp_1.4-6      TH.data_1.0-7      
##  [7] MASS_7.3-45         survival_2.39-5     mvtnorm_1.0-5      
## [10] doBy_4.5-15         car_2.1-2           plyr_1.8.4         
## [13] latticeExtra_0.6-28 RColorBrewer_1.1-2  lattice_0.20-33    
## [16] rmarkdown_1.0       knitr_1.13         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.5        formatR_1.4        nloptr_1.0.4      
##  [4] tools_3.3.1        digest_0.6.9       evaluate_0.9      
##  [7] nlme_3.1-128       mgcv_1.8-12        yaml_2.1.13       
## [10] parallel_3.3.1     SparseM_1.7        stringr_1.0.0     
## [13] MatrixModels_0.4-1 nnet_7.3-12        minqa_1.2.4       
## [16] magrittr_1.5       codetools_0.2-14   htmltools_0.3.5   
## [19] splines_3.3.1      pbkrtest_0.4-6     numDeriv_2014.2-1 
## [22] quantreg_5.26      sandwich_2.3-4     stringi_1.1.1     
## [25] zoo_1.7-13