Atualiza conteúdo da inst/doc.

parent cbef5389
Pipeline #5002 passed with stage
This diff is collapsed.
This diff is collapsed.
......@@ -32,6 +32,8 @@ grid <- cbind(grid[rep(1:nrow(grid), each = length(y)), ],
useOuterStrips(xyplot(py ~ y | factor(lambda) + factor(alpha),
data = grid, type = "h",
xlab = expression(y),
ylab = expression(p(y)),
panel = function(x, y, ...) {
m <- sum(x * y)
panel.xyplot(x, y, ...)
......
This diff is collapsed.
## ----setup, include=FALSE-----------------------------------------
source("_setup.R")
## ---- message=FALSE, error=FALSE, warning=FALSE-------------------
# Definições da sessão.
library(lattice)
library(latticeExtra)
library(grid)
library(plyr)
library(car)
library(doBy)
library(multcomp)
library(MRDCr)
## ---- eval=FALSE, include=FALSE-----------------------------------
# opts_chunk$set(eval = FALSE)
## -----------------------------------------------------------------
data(nematoide, package = "MRDCr")
str(nematoide)
m0 <- glm(nema ~ offset(log(off)) + cult,
data = nematoide, family = poisson)
m1 <- update(m0, family = quasipoisson)
library(lme4)
nematoide$ue <- 1:nrow(nematoide)
m2 <- glmer(nema ~ offset(log(off)) + (1 | cult) + (1 | ue),
data = nematoide, family = poisson)
summary(m2)
# ???
## -----------------------------------------------------------------
data(soja, package = "MRDCr")
str(soja)
soja <- soja[-74, ]
soja$K <- factor(soja$K)
soja$ue <- 1:nrow(soja)
str(soja)
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)
logLik(m2)
logLik(m3)
anova(m2, m3)
summary(m3)
#-----------------------------------------------------------------------
# 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)
## -----------------------------------------------------------------
#-----------------------------------------------------------------------
# devtools::load_all()
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)
summary(confterm)
# 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)
#-----------------------------------------------------------------------
# Ajuste dos modelos.
m0 <- glm(nap ~ galp + resfr * (idade + hmax),
data = confterm, family = poisson)
m1 <- update(m0, family = quasipoisson)
anova(m0, test = "Chisq")
anova(m1, test = "F")
summary(m1)
m2 <- glmer(nap ~ (1 | galp) + resfr * (idade + hmax),
data = confterm, family = poisson)
summary(m2)
anova(m2)
m3 <- update(m2, . ~ . + (1 | ue))
anova(m2, m3)
summary(m3)
anova(m3)
#-----------------------------------------------------------------------
# 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)
# 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))
## ---- echo=FALSE, results="hold"----------------------------------
cat(format(Sys.time(),
format = "Atualizado em %d de %B de %Y.\n\n"))
sessionInfo()
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment