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