|
Modelos de Regressão para Dados de Contagem com o R
|
Análise de Contagens com Modelo Gamma Count
Walmes M. Zeviani, Eduardo E. Ribeiro Jr & Cesar A. Taconeli
A Distribuição Gamma Count
Se uma variável aleatória \(Y\) tem distribuição de probabilidades Gamma Count, então sua função de probabilidade é \[ p(y,\lambda,\alpha) = \left(\int_{0}^{1} \frac{(\alpha\lambda)^{y\alpha}}{\Gamma(y\alpha)}\, u^{y\alpha-1} \exp\{-\alpha\lambda u\}\, \textrm{d}u \right) - \left(\int_{0}^{1} \frac{(\alpha\lambda)^{y\alpha}}{\Gamma((y+1)\alpha)}\, u^{(y+1)\alpha-1} \exp\{-\alpha\lambda u\}\, \textrm{d}u \right), \] em que \(\lambda\) é o parâmetro de locação e \(\alpha\) o parâmetro de dispersão.
Essa função de probabilidade resulta da relação que existe entre intervalo de tempo entre eventos e do número de eventos dentro de um intervalo.
Considere que \(\tau\) seja a variável aleatória tempo entre eventos que acontecem ao longo de um domínio com distribuição \(\tau \sim \text{Gama}(\alpha, \beta)\). Sendo assim, a função desidade de \(\tau\) é
\[ f(\tau, \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} \cdot \tau^{\alpha-1}\cdot \exp\{-\beta\tau\}, \] cuja média é \(\text{E}(\tau) = \frac{\alpha}{\beta}\) e a variância é \(\text{V}(\tau) = \frac{\alpha}{\beta^2}.\)
Ainda, o tempo decorrido até a ocorrência o n-ésimo evento é portanto,
\[ \vartheta_n = \tau_1+\cdots+\tau_n ~ \sim \text{Gama}(n\alpha, \beta), \] pois a soma de variáveis aleatórias Gama tem distribuição Gama, cuja média é \(\text{E}(\vartheta) = \frac{n\alpha}{\beta}\) e a variância é, \(\text{V}(\vartheta) = \frac{n\alpha}{\beta^2}\). A função densidade de \(\vartheta_n\) então é
\[ f_n(\vartheta, \alpha, \beta) = \frac{\beta^{n\alpha}}{\Gamma(n\alpha)}\cdot \vartheta^{n\alpha-1}\cdot \exp\{-\beta\vartheta\}, \] e a função acumulada é
\[ F_n(T) = \Pr(\vartheta_n \leq T) = \int_{0}^{T} \frac{\beta^{n\alpha}}{\Gamma(n\alpha)}\cdot t^{n\alpha-1}\cdot \exp\{-\beta t\}\,\text{d}t. \]
Decorre que, se \([0, T)\) é um intervalo e \(N\) é o número de eventos ocorridos dentro desse intervalo, então \(N < n\) se e somente se \(\vartheta_n \geq T\), ou seja \[ \Pr(N < n) = \Pr(\vartheta_n \geq T) = 1-F_n(T). \]
Como \[ F_n(T) = G(n\alpha, \beta T) = \int_{0}^{T} \frac{\beta^{n\alpha}}{\Gamma(n\alpha)} t^{n\alpha-1}\cdot\exp\{-\beta t\}\, \text{d}t, \] podemos escrever \(\Pr(N = n)\) como sendo a diferença de acumuladas da Gama, \[ \Pr(N = n) = G(n\alpha, \beta T) - G((n+1)\alpha, \beta T). \]
Em resumo, a distribuição para o número de eventos, quando se assume a distribuição do intervalo entre eventos é Gamma, é resultado da diferença de duas acumuladas da distribuição Gamma.
A média da variável de contagem é resultado de \[ \begin{align*} E(N) &= \sum_{i=0}^{\infty} i\cdot \Pr(i) \\ &= \sum_{i=1}^{\infty} i\cdot \Pr(i)\\ &= \sum_{i=1}^{\infty} G(i\alpha, \beta T).\\ \end{align*} \]
Para um \(T\) cada vez maior, tem-se que \[ N(T)\; \dot{\sim}\; \text{Normal}\left( \frac{\beta}{\alpha}, \frac{\beta}{\alpha^2}\right). \]
Considere que \[ \frac{\beta}{\alpha} = \exp\{x^{\top}\theta\} \Rightarrow \beta = \alpha \exp\{x^{\top}\theta\}. \]
Essa parametrização produz um modelo de regressão para a média do tempo entre eventos definida por \[ \text{E}(\tau|x) = \frac{\alpha}{\beta} = \exp\{-x^{\top}\theta\}. \]
O modelo de regressão é para o tempo entre eventos (\(\tau\)) e não diretamente para contagem porque, a menos que \(\alpha = 1\), não é certo que \(\text{E}(N_i|x_i) = [\text{E}(\tau_i|x_i)]^{-1}\). Dessa maneira, \(-\theta\) mede a variação percentual do tempo médio entre eventos para uma unidade de \(x\).
# Definições da sessão.
library(lattice)
library(latticeExtra)
library(grid)
library(bbmle)
library(corrplot)
library(plyr)
library(car)
library(doBy)
library(multcomp)
library(MRDCr)# Função densidade na parametrização original.
dgcnt## function (y, lambda, alpha)
## {
## p <- pgamma(q = 1, shape = y * alpha, rate = alpha * lambda) -
## pgamma(q = 1, shape = (y + 1) * alpha, rate = alpha *
## lambda)
## return(p)
## }
## <environment: namespace:MRDCr># Caso Poisson (gamma == 0).
grid <- expand.grid(lambda = c(2, 8, 15),
alpha = c(0.5, 1, 2.5))
y <- 0:30
py <- mapply(FUN = dgcnt,
lambda = grid$lambda,
alpha = grid$alpha,
MoreArgs = list(y = y), SIMPLIFY = FALSE)
grid <- cbind(grid[rep(1:nrow(grid), each = length(y)), ],
y = y,
py = unlist(py))
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, ...)
panel.abline(v = m, lty = 2)
}),
strip = strip.custom(
strip.names = TRUE,
var.name = expression(lambda == ""),
sep = ""),
strip.left = strip.custom(
strip.names = TRUE,
var.name = expression(alpha == ""),
sep = ""))
A Média da Gamma-Count
#-----------------------------------------------------------------------
# A média da Gamma-Count.
y <- rpois(100, lambda = 5)
L <- list(y = y, X = cbind(rep(1, length(y))))
start <- c(alpha = 0, lambda = 1)
parnames(llgcnt) <- names(start)
n1 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
c(mean(y),
exp(coef(n1)[2]),
calc_mean_gcnt(lambda = exp(coef(n1)[2]),
alpha = exp(coef(n1)[1])))## lambda
## 4.860000 4.840264 4.860013y <- rpois(500, lambda = 50)
L <- list(y = y, X = cbind(rep(1, length(y))))
start <- c(alpha = 0, lambda = 1)
parnames(llgcnt) <- names(start)
n1 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
c(mean(y),
exp(coef(n1)[2]),
calc_mean_gcnt(lambda = exp(coef(n1)[2]),
alpha = exp(coef(n1)[1])))## lambda
## 49.77000 49.76293 49.77000#-----------------------------------------------------------------------
# A E(y) por sum(y * p(py)) e por lambda como função de alpha.
grid <- expand.grid(lambda = exp(seq(-3, 3, length.out = 30)),
alpha = exp(seq(-2, 2, length.out = 11)),
KEEP.OUT.ATTRS = FALSE)
summary(grid)## lambda alpha
## Min. : 0.04979 Min. :0.1353
## 1st Qu.: 0.21188 1st Qu.:0.3012
## Median : 1.00536 Median :1.0000
## Mean : 3.57508 Mean :2.0125
## 3rd Qu.: 4.71960 3rd Qu.:3.3201
## Max. :20.08554 Max. :7.3891y <- 0:500
py <- mapply(FUN = dgcnt,
lambda = grid$lambda,
alpha = grid$alpha,
MoreArgs = list(y = y), SIMPLIFY = FALSE)
grid$m <- sapply(py, FUN = function(py) sum(y * py))
grid$alpha <- round(grid$alpha, 3)
str(grid)## 'data.frame': 330 obs. of 3 variables:
## $ lambda: num 0.0498 0.0612 0.0753 0.0926 0.1139 ...
## $ alpha : num 0.135 0.135 0.135 0.135 0.135 0.135 0.135 0.135 0.135 0.135 ...
## $ m : num 1.12 1.19 1.26 1.34 1.43 ...cl <- brewer.pal(n = 10, name = "Spectral")
xyplot(m ~ lambda, groups = alpha, data = grid, type = "l",
aspect = "iso", grid = TRUE, lwd = 2,
xlab = expression(lambda),
ylab = expression(sum(y %.% f(y))),
auto.key = list(space = "right",
title = expression(alpha),
points = FALSE, lines = TRUE),
par.settings = list(superpose.line = list(col = cl))) +
layer(panel.abline(a = 0, b = 1, lty = 2)) +
layer(panel.rect(xleft = 0, ybottom = 0,
xright = 5, ytop = 5, lty = 2)) +
layer(panel.arrows(7, 13, 13, 7, length = 0.1)) +
layer(panel.text(x = 13, y = 7,
labels = expression(alpha), pos = 4))
xyplot(m ~ lambda, groups = alpha, type = "l",
data = subset(grid, findInterval(lambda, vec = c(0, 5)) == 1),
xlab = expression(lambda),
ylab = expression(sum(y %.% f(y))),
aspect = "iso", grid = TRUE, lwd = 2,
auto.key = list(space = "right",
title = expression(alpha),
points = FALSE, lines = TRUE),
par.settings = list(superpose.line = list(col = cl))) +
layer(panel.abline(a = 0, b = 1, lty = 2)) +
layer(panel.abline(a = -0.5, b = 1, lty = 2)) +
layer(panel.abline(a = 1, b = 1, lty = 2))
Função de Risco
h <- function(...) {
dgamma(...)/(1 - pgamma(...))
}
shape <- seq(0.5, 1.5, by = 0.1)
col <- brewer.pal(n = 5, name = "Spectral")
col <- colorRampPalette(colors = col)(length(shape))
curve(dgamma(x, shape = shape[1], rate = 1),
from = 0, to = 5, col = col[1], lwd = 2,
xlab = expression(tau),
ylab = expression(f(tau)))
for (s in 2:length(shape)) {
curve(dgamma(x, shape = shape[s], rate = 1),
add = TRUE, col = col[s], lwd = 2)
}
legend("topright", legend = sprintf("%0.3f", shape),
col = col, lty = 1, lwd = 2, bty = "n",
title = expression(alpha))
curve(h(x, shape = shape[1], rate = 1),
from = 0, to = 10, col = col[1], lwd = 2,
ylim = c(0, 2.5),
xlab = expression(tau),
ylab = expression(f(tau)/(1 - F(tau))))
for (s in 2:length(shape)) {
curve(h(x, shape = shape[s], rate = 1), add = TRUE,
col = col[s], lwd = 2)
}
legend("topright", legend = sprintf("%0.3f", shape),
col = col, lty = 1, lwd = 2, bty = "n",
title = expression(alpha))
Verossimilhança e Estimação
#-----------------------------------------------------------------------
# Função de log-Verossimilhança da Poisson Generalizada na
# parametrização de modelo de regressão.
MRDCr::llgcnt## function (params, y, X, offset = NULL)
## {
## if (is.null(offset)) {
## offset <- 1L
## }
## alpha <- exp(params[1])
## eXb <- exp(X %*% params[-1])
## alpha.eXb <- alpha * eXb
## alpha.y <- alpha * y
## ll <- -sum(log(pgamma(offset, shape = alpha.y, rate = alpha.eXb) -
## pgamma(offset, shape = alpha.y + alpha, rate = alpha.eXb)))
## return(ll)
## }
## <environment: namespace:MRDCr>#-----------------------------------------------------------------------
# Gerando uma amostra aleatória da Poisson, para teste.
# Offset = 2, lambda = 3.
y <- rpois(100, lambda = 2 * 3)
L <- list(y = y,
offset = rep(2, length(y)),
X = cbind(rep(1, length(y))))
start <- c(alpha = 0, lambda = 1)
parnames(llgcnt) <- names(start)
# Como \alpha foi fixado em 1, essa ll corresponde à Poisson.
n0 <- mle2(minuslogl = llgcnt,
start = start, data = L,
fixed = list(alpha = 0), vecpar = TRUE)
# Para conferir.
c(coef(n0)["lambda"],
coef(glm(y ~ offset(log(L$offset)), family = poisson)))## lambda (Intercept)
## 1.131402 1.131402# Estimando o \alpha.
n1 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
coef(n1)## alpha lambda
## 0.006955925 1.131961388# Perfil de verossimilhança dos parâmetros.
plot(profile(n1))
# Covariância.
V <- cov2cor(vcov(n1))
corrplot.mixed(V, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()## null device
## 1Número de Vagens Produzidas em Soja
Resultados de um experimento fatorial (3 x 5), em delineamento de blocos casualizados, que estudou a influência de níveis de potássio na adubação de soja em combinação com irrigação em casa de vegetação. As variáveis de contagem registradas nesse experimento foram o número de vagens viáveis (e não viáveis) e o número total de sementes por parcela com duas plantas de soja.
#-----------------------------------------------------------------------
# Carregando e explorando os dados.
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 ...# help(soja, package = "MRDCr", help_type = "html")
# A observação 74 é um outlier.
soja <- soja[-74, ]
xyplot(nvag ~ K | umid, data = soja, layout = c(NA, 1),
type = c("p", "smooth"),
ylab = "Número de vagens por parcela",
xlab = expression("Dose de potássio aplicada"~(mg ~ dm^{3})),
strip = strip.custom(strip.names = TRUE, var.name = "Umidade"))
soja <- transform(soja, K = factor(K))
#-----------------------------------------------------------------------
# Modelo Poisson.
m0 <- glm(nvag ~ bloc + umid * K, data = soja, family = poisson)
m1 <- update(m0, family = quasipoisson)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Medidas de decisão.
# anova(m0, test = "Chisq")
anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: nvag
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 73 322.68
## bloc 4 14.284 69 308.40 2.9785 0.02684 *
## umid 2 92.911 67 215.49 38.7485 3.157e-11 ***
## K 4 136.060 63 79.43 28.3719 8.316e-13 ***
## umid:K 8 14.150 55 65.28 1.4754 0.18750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m1)##
## Call:
## glm(formula = nvag ~ bloc + umid * K, family = quasipoisson,
## data = soja)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9681 -0.6393 -0.0195 0.4854 2.3306
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.95369 0.07547 52.389 < 2e-16 ***
## blocII -0.02928 0.04479 -0.654 0.516036
## blocIII -0.07265 0.04529 -1.604 0.114420
## blocIV -0.12544 0.04592 -2.732 0.008457 **
## blocV -0.10795 0.04705 -2.294 0.025606 *
## umid50 0.13404 0.09597 1.397 0.168117
## umid62,5 0.21656 0.09418 2.299 0.025303 *
## K30 0.27427 0.09300 2.949 0.004670 **
## K60 0.30797 0.09233 3.336 0.001530 **
## K120 0.32883 0.09192 3.577 0.000734 ***
## K180 0.25540 0.09338 2.735 0.008376 **
## umid50:K30 0.06322 0.12654 0.500 0.619324
## umid62,5:K30 -0.10747 0.12631 -0.851 0.398532
## umid50:K60 0.16561 0.12449 1.330 0.188897
## umid62,5:K60 0.10735 0.12286 0.874 0.386028
## umid50:K120 0.14920 0.12414 1.202 0.234563
## umid62,5:K120 0.11839 0.12487 0.948 0.347229
## umid50:K180 0.30370 0.12439 2.441 0.017873 *
## umid62,5:K180 0.19838 0.12325 1.610 0.113208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.198902)
##
## Null deviance: 322.684 on 73 degrees of freedom
## Residual deviance: 65.278 on 55 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4#-----------------------------------------------------------------------
# Modelo Poisson Generalizado.
L <- with(soja,
list(y = nvag, offset = 1, X = model.matrix(m0)))
# Usa as estimativas do Poisson como valores iniciais para a PGen.
start <- c(alpha = 0, coef(m0))
parnames(llgcnt) <- names(start)
# Com alpha fixo em 0 corresponde à Poisson.
m2 <- mle2(llgcnt, start = start, data = L,
fixed = list(alpha = 0), vecpar = TRUE)
# Mesma medida de ajuste e estimativas.
c(logLik(m2), logLik(m0))## [1] -259.6165 -259.6165cbind(coef(m2)[-1], coef(m0))## [,1] [,2]
## (Intercept) 3.95368723 3.95368724
## blocII -0.02927855 -0.02927854
## blocIII -0.07265048 -0.07265048
## blocIV -0.12543798 -0.12543798
## blocV -0.10794857 -0.10794857
## umid50 0.13404356 0.13404356
## umid62,5 0.21656458 0.21656458
## K30 0.27427290 0.27427290
## K60 0.30796674 0.30796674
## K120 0.32883188 0.32883188
## K180 0.25540441 0.25540441
## umid50:K30 0.06322288 0.06322288
## umid62,5:K30 -0.10747272 -0.10747272
## umid50:K60 0.16561471 0.16561471
## umid62,5:K60 0.10735066 0.10735066
## umid50:K120 0.14920392 0.14920392
## umid62,5:K120 0.11838578 0.11838578
## umid50:K180 0.30369921 0.30369921
## umid62,5:K180 0.19837927 0.19837927# Gamma-Count estimando o alpha.
m3 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
# Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0).
anova(m3, m2)## Likelihood Ratio Tests
## Model 1: m3, [llgcnt]: alpha+(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: m2, [llgcnt]: (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.5802 1 0.4462# Estimaitvas dos parâmetros.
c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0## PoissonGLM PoissonML GCount
## NA 0.00000000 0.12875777
## (Intercept) 3.95368724 3.95368723 3.95487435
## blocII -0.02927854 -0.02927855 -0.02925363
## blocIII -0.07265048 -0.07265048 -0.07259276
## blocIV -0.12543798 -0.12543798 -0.12534177
## blocV -0.10794857 -0.10794857 -0.10785674
## umid50 0.13404356 0.13404356 0.13388357
## umid62,5 0.21656458 0.21656458 0.21632039
## K30 0.27427290 0.27427290 0.27397250
## K60 0.30796674 0.30796674 0.30763762
## K120 0.32883188 0.32883188 0.32848139
## K180 0.25540441 0.25540441 0.25512492
## umid50:K30 0.06322288 0.06322288 0.06321988
## umid62,5:K30 -0.10747272 -0.10747272 -0.10732539
## umid50:K60 0.16561471 0.16561471 0.16553912
## umid62,5:K60 0.10735066 0.10735066 0.10734196
## umid50:K120 0.14920392 0.14920392 0.14914452
## umid62,5:K120 0.11838578 0.11838578 0.11838048
## umid50:K180 0.30369921 0.30369921 0.30351780
## umid62,5:K180 0.19837927 0.19837927 0.19829595splom(c0[-(1:2), ]) + layer(panel.abline(a = 0, b = 1, lty = 2))
# Perfil para o parâmetro de dispersão.
plot(profile(m3, which = "alpha"))
abline(v = 0, lty = 2)
V <- cov2cor(vcov(m3))
corrplot.mixed(V, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()## null device
## 1# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))## sum mean max
## 0.061875318 0.003256596 0.022121891#-----------------------------------------------------------------------
# Testes de hipótese.
# Teste de Wald para a interação.
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a == max(a)
L <- t(replicate(sum(ai), rbind(coef(m3) * 0), simplify = "matrix"))
L[, ai] <- diag(sum(ai))
# Cáclculo da estatística Chi-quadrado.
crossprod(L %*% coef(m3),
solve(L %*% vcov(m3) %*% t(L),
L %*% coef(m3)))## [,1]
## [1,] 15.96558# Teste de Wald para interação (poderia ser LRT, claro).
# É necessário passar um objeto glm mesmo fornecendo o restante a parte.
linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
coef. = coef(m3))## Linear hypothesis test
##
## Hypothesis:
## umid50:K30 = 0
## umid62,5:K30 = 0
## umid50:K60 = 0
## umid62,5:K60 = 0
## umid50:K120 = 0
## umid62,5:K120 = 0
## umid50:K180 = 0
## umid62,5:K180 = 0
##
## Model 1: restricted model
## Model 2: nvag ~ bloc + umid * K
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 63
## 2 55 8 15.966 0.04288 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Predição com bandas de confiança.
X <- LSmatrix(m0, effect = c("umid", "K"))
pred <- attr(X, "grid")
pred <- transform(pred,
K = as.integer(K),
umid = factor(umid))
pred <- list(pois = pred, quasi = pred, gcnt = pred)
# Preditos pela Poisson.
aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
# Preditos pela Quasi-Poisson.
aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))
# Preditos pela Gamma-Count.
aux <- predict(m3, newdata = X,
interval = "confidence", type = "link")
pred$gcnt <- cbind(pred$gcnt, exp(aux[, c(2, 1, 3)]))
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",
"Gamma-Count")))
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 vagens por parcela",
ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
prepanel = prepanel.cbH,
desloc = 8 * scale(as.integer(pred$modelo), scale = FALSE),
panel = panel.cbH)
Número de Grãos Produzidas em Soja
Análise do número de grãos por pacela do experimento com soja.
#-----------------------------------------------------------------------
xyplot(ngra ~ K | umid, data = soja, layout = c(NA, 1),
type = c("p", "smooth"),
ylab = "Número de grãos por parcela",
xlab = expression("Dose de potássio aplicada"~(mg ~ dm^{3})),
strip = strip.custom(strip.names = TRUE, var.name = "Umidade"))
# NOTE: Essa contagem é alta e uma análise preliminar não retornou
# Hessiana para o modelo Gamma-Count ajustado com a mle2. A suspeita que
# é seja pela ordem de magnitude dos dados. Sendo assim, vamos
# considerar um offset artifical de 10 apenas para correr as análises.
#
# Warning message:
# In mle2(llgcnt, start = start, data = L, fixed = list(alpha = 0), :
# couldn't invert Hessian
#
# Warning message:
# In mle2(llgcnt, start = start, data = L, vecpar = TRUE) :
# couldn't invert Hessian
soja$off <- 10
fivenum(with(soja, ngra/off))## [1] 9.20 14.70 17.05 21.60 27.10#-----------------------------------------------------------------------
# Modelo Poisson.
m0 <- glm(ngra ~ offset(log(off)) + bloc + umid * K,
data = soja, family = poisson)
m1 <- update(m0, family = quasipoisson)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Medidas de decisão.
# anova(m0, test = "Chisq")
anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: ngra
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 73 761.43
## bloc 4 28.34 69 733.09 3.0962 0.02273 *
## umid 2 185.06 67 548.03 40.4318 1.584e-11 ***
## K 4 380.32 63 167.71 41.5466 5.186e-16 ***
## umid:K 8 42.99 55 124.72 2.3480 0.03004 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m1)##
## Call:
## glm(formula = ngra ~ offset(log(off)) + bloc + umid * K, family = quasipoisson,
## data = soja)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6297 -1.0707 -0.0873 0.7428 3.2810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.49937 0.06773 36.902 < 2e-16 ***
## blocII -0.01939 0.04017 -0.483 0.631230
## blocIII -0.03663 0.04035 -0.908 0.367964
## blocIV -0.10559 0.04107 -2.571 0.012884 *
## blocV -0.09313 0.04218 -2.208 0.031433 *
## umid50 0.13245 0.08611 1.538 0.129739
## umid62,5 0.18548 0.08506 2.180 0.033517 *
## K30 0.29799 0.08299 3.591 0.000703 ***
## K60 0.34434 0.08218 4.190 0.000102 ***
## K120 0.36493 0.08183 4.459 4.1e-05 ***
## K180 0.29542 0.08303 3.558 0.000778 ***
## umid50:K30 0.04344 0.11318 0.384 0.702600
## umid62,5:K30 -0.13669 0.11386 -1.201 0.235081
## umid50:K60 0.11559 0.11136 1.038 0.303817
## umid62,5:K60 0.09174 0.11027 0.832 0.409039
## umid50:K120 0.11860 0.11088 1.070 0.289461
## umid62,5:K120 0.16266 0.11145 1.460 0.150096
## umid50:K180 0.28832 0.11085 2.601 0.011917 *
## umid62,5:K180 0.21569 0.11021 1.957 0.055430 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 2.28854)
##
## Null deviance: 761.43 on 73 degrees of freedom
## Residual deviance: 124.72 on 55 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4#-----------------------------------------------------------------------
# Modelo Poisson Generalizado.
L <- with(soja,
list(y = ngra, offset = off, X = model.matrix(m0)))
# Usa as estimativas do Poisson como valores iniciais.
start <- c(alpha = 0, coef(m0))
parnames(llgcnt) <- names(start)
# Com alpha fixo em 0 corresponde à Poisson.
m2 <- mle2(llgcnt, start = start, data = L,
fixed = list(alpha = 0), vecpar = TRUE)
# Mesma medida de ajuste e estimativas.
c(logLik(m2), logLik(m0))## [1] -321.6692 -321.6692# Poisson Generalizada.
m3 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
# Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0).
anova(m3, m2)## Likelihood Ratio Tests
## Model 1: m3, [llgcnt]: alpha+(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: m2, [llgcnt]: (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.051 1 0.0005176 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimaitvas dos parâmetros.
c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0## PoissonGLM PoissonML GCount
## NA 0.00000000 -0.52305059
## (Intercept) 2.49937463 2.49937463 2.49649191
## blocII -0.01939070 -0.01939070 -0.01943043
## blocIII -0.03662632 -0.03662632 -0.03669577
## blocIV -0.10559332 -0.10559332 -0.10578690
## blocV -0.09313375 -0.09313375 -0.09331619
## umid50 0.13245136 0.13245136 0.13282522
## umid62,5 0.18548293 0.18548293 0.18598885
## K30 0.29799144 0.29799144 0.29876306
## K60 0.34433662 0.34433662 0.34520871
## K120 0.36493092 0.36493092 0.36585027
## K180 0.29542405 0.29542405 0.29619069
## umid50:K30 0.04343930 0.04343930 0.04341282
## umid62,5:K30 -0.13669277 -0.13669277 -0.13709575
## umid50:K60 0.11559375 0.11559375 0.11568092
## umid62,5:K60 0.09174072 0.09174072 0.09174645
## umid50:K120 0.11859658 0.11859658 0.11867938
## umid62,5:K120 0.16266223 0.16266223 0.16275434
## umid50:K180 0.28832017 0.28832017 0.28870557
## umid62,5:K180 0.21568848 0.21568848 0.21591249splom(c0[-(1:2), ]) + layer(panel.abline(a = 0, b = 1, lty = 2))
# Perfil para o parâmetro de dispersão.
plot(profile(m3, which = "alpha"))
abline(v = 0, lty = 2)
V <- cov2cor(vcov(m3))
corrplot.mixed(V, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()## null device
## 1# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))## sum mean max
## 0.054983277 0.002893857 0.020065055# Teste de Wald para a interação.
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a == max(a)
L <- t(replicate(sum(ai), rbind(coef(m3) * 0), simplify = "matrix"))
L[, ai] <- diag(sum(ai))
# Cáclculo da estatística Chi-quadrado.
crossprod(L %*% coef(m3),
solve(L %*% vcov(m3) %*% t(L),
L %*% coef(m3)))## [,1]
## [1,] 25.28744linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
coef. = coef(m3))## Linear hypothesis test
##
## Hypothesis:
## umid50:K30 = 0
## umid62,5:K30 = 0
## umid50:K60 = 0
## umid62,5:K60 = 0
## umid50:K120 = 0
## umid62,5:K120 = 0
## umid50:K180 = 0
## umid62,5:K180 = 0
##
## Model 1: restricted model
## Model 2: ngra ~ offset(log(off)) + bloc + umid * K
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 63
## 2 55 8 25.287 0.001389 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Predição com bandas de confiança.
# Por causa do offset, não tem como usar a LSmatrix.
pred <- unique(subset(soja, select = c("umid", "K")))
X <- model.matrix(formula(m0)[-2],
data = cbind(off = 1, bloc = soja$bloc[1], pred))
i <- grep(x = colnames(X), pattern = "^bloc")
X[, i] <- X[, i] * 0 + 1/(length(i) + 1)
pred <- transform(pred,
K = as.numeric(as.character(K)),
umid = factor(umid))
pred <- list(pois = pred, quasi = pred, gcnt = pred)
# Preditos pela Poisson.
aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
# Preditos pela Quasi-Poisson.
aux <- confint(glht(m1, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))
# Preditos pela Gamma-Count.
aux <- predict(m3, newdata = X,
interval = "confidence", type = "link")
pred$gcnt <- cbind(pred$gcnt, exp(aux[, c(2, 1, 3)]))
# Junta o resultado dos 3 modelos.
pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, umid, K, modelo)
str(pred)## 'data.frame': 45 obs. of 6 variables:
## $ modelo: Factor w/ 3 levels "pois","quasi",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ umid : Factor w/ 3 levels "37,5","50","62,5": 1 1 1 1 1 1 1 1 1 1 ...
## $ K : num 0 0 0 30 30 30 60 60 60 120 ...
## $ fit : num 11.6 11.6 11.5 15.6 15.6 ...
## $ lwr : num 10.7 10.2 10.4 14.5 14 ...
## $ upr : num 12.6 13.1 12.8 16.7 17.3 ...key <- list(type = "o", divide = 1,
lines = list(pch = 1:nlevels(pred$modelo),
lty = 1, col = 1),
text = list(c("Poisson", "Quasi-Poisson",
"Gamma-Count")))
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 = 8 * scale(as.integer(pred$modelo), scale = FALSE),
panel = panel.cbH)
Número de Grãos por Vagem
#-----------------------------------------------------------------------
# Número de grãos por vagem (offset).
xyplot(ngra/nvag ~ K | umid, data = soja, layout = c(NA, 1),
type = c("p", "smooth"),
xlab = expression("Dose de potássio"~(mg~dm^{-3})),
ylab = "Número de grãos por vagem",
strip = strip.custom(strip.names = TRUE, var.name = "Umidade"))
#-----------------------------------------------------------------------
# Modelo Poisson.
m0 <- glm(ngra ~ offset(log(nvag)) + bloc + umid * K,
data = soja, family = poisson)
m1 <- update(m0, family = quasipoisson)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Medidas de decisão.
# anova(m0, test = "Chisq")
anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: ngra
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 73 34.238
## bloc 4 2.0264 69 32.212 1.1637 0.33686
## umid 2 1.7457 67 30.466 2.0050 0.14439
## K 4 3.8324 63 26.634 2.2009 0.08082 .
## umid:K 8 2.6490 55 23.984 0.7606 0.63838
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m1)##
## Call:
## glm(formula = ngra ~ offset(log(nvag)) + bloc + umid * K, family = quasipoisson,
## data = soja)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.52741 -0.46025 0.06782 0.45274 1.05490
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.848160 0.029477 28.774 <2e-16 ***
## blocII 0.011313 0.017535 0.645 0.5215
## blocIII 0.036250 0.017615 2.058 0.0444 *
## blocIV 0.019384 0.017940 1.080 0.2846
## blocV 0.016100 0.018448 0.873 0.3866
## umid50 -0.001908 0.037574 -0.051 0.9597
## umid62,5 -0.031173 0.037124 -0.840 0.4047
## K30 0.022867 0.036200 0.632 0.5302
## K60 0.034499 0.035862 0.962 0.3403
## K120 0.035322 0.035703 0.989 0.3268
## K180 0.040586 0.036221 1.121 0.2674
## umid50:K30 -0.018692 0.049385 -0.378 0.7065
## umid62,5:K30 -0.028459 0.049680 -0.573 0.5691
## umid50:K60 -0.048182 0.048583 -0.992 0.3257
## umid62,5:K60 -0.014035 0.048103 -0.292 0.7716
## umid50:K120 -0.030438 0.048367 -0.629 0.5317
## umid62,5:K120 0.045524 0.048714 0.935 0.3541
## umid50:K180 -0.016302 0.048354 -0.337 0.7373
## umid62,5:K180 0.016911 0.048076 0.352 0.7264
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.4353299)
##
## Null deviance: 34.238 on 73 degrees of freedom
## Residual deviance: 23.984 on 55 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 3# Declara um modelo aditivo.
m0 <- glm(ngra ~ offset(log(nvag)) + bloc + umid + K,
data = soja, family = poisson)
m1 <- update(m0, family = quasipoisson)
anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: ngra
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 73 34.238
## bloc 4 2.0264 69 32.212 1.2021 0.31880
## umid 2 1.7457 67 30.466 2.0711 0.13455
## K 4 3.8324 63 26.634 2.2734 0.07114 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Modelo Poisson Generalizado.
L <- with(soja,
list(y = ngra, offset = nvag, X = model.matrix(m0)))
start <- c(alpha = 0, coef(m0))
parnames(llgcnt) <- names(start)
# Com alpha fixo em 0 corresponde à Poisson.
m2 <- mle2(llgcnt, start = start, data = L,
fixed = list(alpha = 0), vecpar = TRUE)
# Mesma medida de ajuste e estimativas.
c(logLik(m2), logLik(m0))## [1] -272.6263 -272.6263# Poisson Generalizada.
m3 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
# Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0).
anova(m3, m2)## Likelihood Ratio Tests
## Model 1: m3, [llgcnt]: alpha+(Intercept)+blocII+blocIII+blocIV+
## blocV+umid50+umid62,5+K30+K60+K120+K180
## Model 2: m2, [llgcnt]: (Intercept)+blocII+blocIII+blocIV+blocV+
## umid50+umid62,5+K30+K60+K120+K180
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 12 517.26
## 2 11 545.25 27.989 1 1.22e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0## PoissonGLM PoissonML GCount
## NA 0.000000000 1.01872022
## (Intercept) 0.855552235 0.855552234 0.85824271
## blocII 0.011126118 0.011126118 0.01116767
## blocIII 0.036245340 0.036245340 0.03631263
## blocIV 0.018895822 0.018895822 0.01904826
## blocV 0.012659290 0.012659290 0.01284491
## umid50 -0.026631767 -0.026631767 -0.02708107
## umid62,5 -0.026250165 -0.026250165 -0.02670737
## K30 0.007448488 0.007448488 0.00687900
## K60 0.012551823 0.012551823 0.01172824
## K120 0.039757260 0.039757260 0.03886651
## K180 0.041632187 0.041632187 0.04072724splom(c0[-(1:2), ]) + layer(panel.abline(a = 0, b = 1, lty = 2))
# Perfil para o parâmetro de dispersão.
plot(profile(m3, which = "alpha"))
V <- cov2cor(vcov(m3))
corrplot.mixed(V, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()## null device
## 1# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))## sum mean max
## 0.038796736 0.003526976 0.013298824# Teste de Wald para a interação.
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a == max(a)
L <- t(replicate(sum(ai), rbind(coef(m3) * 0), simplify = "matrix"))
L[, ai] <- diag(sum(ai))
# Cáclculo da estatística Chi-quadrado.
crossprod(L %*% coef(m3),
solve(L %*% vcov(m3) %*% t(L),
L %*% coef(m3)))## [,1]
## [1,] 10.33342linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
coef. = coef(m3))## Linear hypothesis test
##
## Hypothesis:
## K30 = 0
## K60 = 0
## K120 = 0
## K180 = 0
##
## Model 1: restricted model
## Model 2: ngra ~ offset(log(nvag)) + bloc + umid + K
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 67
## 2 63 4 10.333 0.03517 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Predição com bandas de confiança.
# Por causa do offset, não tem como usar a LSmatrix.
pred <- unique(subset(soja, select = c("umid", "K")))
X <- model.matrix(formula(m0)[-2],
data = cbind(nvag = 1, bloc = soja$bloc[1], pred))
i <- grep(x = colnames(X), pattern = "^bloc")
X[, i] <- X[, i] * 0 + 1/(length(i) + 1)
head(X)## (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 0 0 1 0
## 3 1 0.2 0.2 0.2 0.2 0 0 0 1
## 4 1 0.2 0.2 0.2 0.2 0 0 0 0
## 5 1 0.2 0.2 0.2 0.2 0 0 0 0
## 6 1 0.2 0.2 0.2 0.2 1 0 0 0
## K120 K180
## 1 0 0
## 2 0 0
## 3 0 0
## 4 1 0
## 5 0 1
## 6 0 0pred <- list(pois = pred, quasi = pred, gcnt = pred)
# Preditos pela Poisson.
aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
# Preditos pela Quasi-Poisson.
aux <- confint(glht(m1, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))
# Comparando as estimativas de média para contagem.
cbind(Pois = pred$pois$fit,
Gcnt1 = c(exp(predict(m3, newdata = X))),
GCnt2 = c(predict(m3, newdata = X, type = "response")))## Pois Gcnt1 GCnt2
## [1,] 2.390106 2.396759 2.077286
## [2,] 2.407975 2.413303 2.093831
## [3,] 2.420295 2.425034 2.105562
## [4,] 2.487044 2.491747 2.172275
## [5,] 2.491711 2.496387 2.176916
## [6,] 2.327293 2.332723 2.013249
## [7,] 2.344693 2.348825 2.029352
## [8,] 2.356689 2.360243 2.040770
## [9,] 2.421684 2.425173 2.105701
## [10,] 2.426228 2.429690 2.110218
## [11,] 2.328181 2.333595 2.014121
## [12,] 2.345587 2.349703 2.030229
## [13,] 2.357588 2.361125 2.041652
## [14,] 2.422608 2.426079 2.106607
## [15,] 2.427154 2.430598 2.111126# Preditos pela Gamma-Count.
aux <- predict(m3, newdata = X,
interval = "confidence", type = "link")
pred$gcnt <- cbind(pred$gcnt, exp(aux[, c(2, 1, 3)]))
# Junta o resultado dos 3 modelos.
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",
"Gamma-Count")))
xyplot(ngra/nvag ~ K | umid, data = soja, layout = c(NA, 1),
type = c("p", "smooth"),
xlim = extendrange(range(as.numeric(pred$K)), f = 0.2),
key = key,
xlab = expression("Dose de potássio"~(mg~dm^{-3})),
ylab = "Número de grãos por vagem",
strip = strip.custom(strip.names = TRUE, var.name = "Umidade")) +
as.layer(
xyplot(fit ~ K | umid, data = pred,
layout = c(NA, 1), as.table = TRUE,
pch = pred$modelo,
ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
prepanel = prepanel.cbH,
desloc = 0.15 * scale(as.integer(pred$modelo),
scale = FALSE),
panel = panel.cbH))
Número de Capulhos Produzidos em Algodão
Experimento conduzido em casa de vegetação para avaliar o efeito da desfolha, em diferentes fases de desenvolvimento do algodão, sobre a produção da cultura. As parcelas foram vasos com duas plantas que tiveram a área das folhas removida com uma tesoura, simulando o ataque de insetos desfolhadores, nos níveis de 0, 25, 50, 75 e 100% de remoção de área foliar. Em cada fase de desenvolvimento (de 5), 5 parcelas sofreram desfolha nos níveis mencionados. Trata-se então de um experimento em arranjo fatorial 5 desfolhas \(\times\) 5 fases com 5 repetições por cela. O número de capulhos ao final do experimento foi contado.
#-----------------------------------------------------------------------
# Número de capulhos em função do nível de desfolha artificial e fase
# fenológica do algodoeiro.
data(capdesfo, package = "MRDCr")
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, package = "MRDCr", help_type = "html")
p1 <- xyplot(ncap ~ des | est, data = capdesfo,
col = 1, type = c("p", "smooth"), col.line = "gray50",
layout = c(2, 3), as.table = TRUE,
xlim = extendrange(capdesfo$des, f = 0.15),
xlab = "Nível de desfolhas artificial",
ylab = "Número de capulho produzidos",
spread = 0.035, panel = panel.beeswarm)
p1
#-----------------------------------------------------------------------
# Modelo Poisson.
m0 <- glm(ncap ~ est * (des + I(des^2)),
data = capdesfo, family = poisson)
m1 <- update(m0, family = quasipoisson)
par(mfrow = c(2, 2))
plot(m0); layout(1)
# anova(m0, test = "Chisq")
anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: ncap
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 124 75.514
## est 4 19.9582 120 55.556 21.5543 3.775e-13 ***
## des 1 15.8643 119 39.692 68.5317 3.269e-13 ***
## I(des^2) 1 1.2926 118 38.399 5.5839 0.01988 *
## est:des 4 6.7085 114 31.691 7.2450 3.236e-05 ***
## est:I(des^2) 4 6.3592 110 25.331 6.8677 5.674e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m1)##
## Call:
## glm(formula = ncap ~ est * (des + I(des^2)), family = quasipoisson,
## data = capdesfo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.07706 -0.30981 -0.02283 0.27044 1.16650
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.21424 0.06706 33.018 < 2e-16 ***
## estbotão floral -0.08003 0.09655 -0.829 0.40898
## estflorecimento -0.02723 0.09626 -0.283 0.77781
## estmaça -0.14861 0.09867 -1.506 0.13489
## estcapulho 0.11291 0.09247 1.221 0.22465
## des 0.34859 0.32741 1.065 0.28935
## I(des^2) -0.73843 0.32397 -2.279 0.02458 *
## estbotão floral:des 0.13638 0.46401 0.294 0.76937
## estflorecimento:des -1.58189 0.49140 -3.219 0.00169 **
## estmaça:des 0.47548 0.49046 0.969 0.33444
## estcapulho:des -0.82100 0.45204 -1.816 0.07206 .
## estbotão floral:I(des^2) 0.10435 0.45451 0.230 0.81884
## estflorecimento:I(des^2) 1.40435 0.49033 2.864 0.00501 **
## estmaça:I(des^2) -0.92938 0.49667 -1.871 0.06397 .
## estcapulho:I(des^2) 1.07565 0.44314 2.427 0.01683 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.2314882)
##
## Null deviance: 75.514 on 124 degrees of freedom
## Residual deviance: 25.331 on 110 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4#-----------------------------------------------------------------------
# Modelo Poisson Generalizada.
L <- with(capdesfo,
list(y = ncap, offset = 1, X = model.matrix(m0)))
start <- c(alpha = log(1), coef(m0))
parnames(llgcnt) <- names(start)
# Modelo Poisson também.
m2 <- mle2(llgcnt, start = start, data = L,
fixed = list(alpha = 0), vecpar = TRUE)
c(logLik(m2), logLik(m0))## [1] -254.8414 -254.8414# Modelo Poisson Generalizado.
m3 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
logLik(m3)## 'log Lik.' -203.7501 (df=16)anova(m3, m2)## Likelihood Ratio Tests
## Model 1: m3, [llgcnt]: alpha+(Intercept)+estbotão floral+
## estflorecimento+estmaça+estcapulho+des+I(des^2)+
## estbotão floral:des+estflorecimento:des+estmaça:des+
## estcapulho:des+estbotão floral:I(des^2)+
## estflorecimento:I(des^2)+estmaça:I(des^2)+
## estcapulho:I(des^2)
## Model 2: m2, [llgcnt]: (Intercept)+estbotão floral+
## estflorecimento+estmaça+estcapulho+des+I(des^2)+
## estbotão floral:des+estflorecimento:des+estmaça:des+
## estcapulho:des+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 16 407.50
## 2 15 509.68 102.18 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m3)## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = llgcnt, start = start, data = L, vecpar = TRUE)
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## alpha 1.710979 0.135195 12.6556 < 2.2e-16 ***
## (Intercept) 2.258038 0.059339 38.0532 < 2.2e-16 ***
## estbotão floral -0.076524 0.085405 -0.8960 0.3702456
## estflorecimento -0.025313 0.085113 -0.2974 0.7661593
## estmaça -0.139845 0.087219 -1.6034 0.1088485
## estcapulho 0.108438 0.081763 1.3263 0.1847562
## des 0.329398 0.289643 1.1373 0.2554329
## I(des^2) -0.699663 0.286619 -2.4411 0.0146430 *
## estbotão floral:des 0.133744 0.410476 0.3258 0.7445554
## estflorecimento:des -1.501986 0.434453 -3.4572 0.0005459 ***
## estmaça:des 0.421813 0.433578 0.9729 0.3306201
## estcapulho:des -0.782023 0.399797 -1.9561 0.0504593 .
## estbotão floral:I(des^2) 0.094333 0.402116 0.2346 0.8145248
## estflorecimento:I(des^2) 1.338205 0.433533 3.0867 0.0020236 **
## estmaça:I(des^2) -0.833309 0.439023 -1.8981 0.0576829 .
## estcapulho:I(des^2) 1.022210 0.391972 2.6079 0.0091108 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 407.5003c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0## PoissonGLM PoissonML GCount
## NA 0.00000000 1.71097890
## (Intercept) 2.21424242 2.21424242 2.25803759
## estbotão floral -0.08002756 -0.08002756 -0.07652439
## estflorecimento -0.02722855 -0.02722855 -0.02531284
## estmaça -0.14861263 -0.14861263 -0.13984549
## estcapulho 0.11291268 0.11291268 0.10843781
## des 0.34858564 0.34858564 0.32939763
## I(des^2) -0.73843427 -0.73843427 -0.69966320
## estbotão floral:des 0.13638395 0.13638395 0.13374406
## estflorecimento:des -1.58188650 -1.58188650 -1.50198575
## estmaça:des 0.47548300 0.47548300 0.42181332
## estcapulho:des -0.82100307 -0.82100307 -0.78202287
## estbotão floral:I(des^2) 0.10435096 0.10435096 0.09433350
## estflorecimento:I(des^2) 1.40435178 1.40435178 1.33820484
## estmaça:I(des^2) -0.92937720 -0.92937720 -0.83330941
## estcapulho:I(des^2) 1.07565430 1.07565430 1.02220995splom(c0[-(1:2), ]) + layer(panel.abline(a = 0, b = 1, lty = 2))
# Perfil para o parâmetro de dispersão.
plot(profile(m3, which = "alpha"))
V <- cov2cor(vcov(m3))
corrplot.mixed(V, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()## null device
## 1# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))## sum mean max
## 0.070733515 0.004715568 0.022155377# Teste de Wald para a interação.
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a == max(a)
L <- t(replicate(sum(ai), rbind(coef(m3) * 0), simplify = "matrix"))
L[, ai] <- diag(sum(ai))
# Teste de Wald explicito.
crossprod(L %*% coef(m3),
solve(L %*% vcov(m3) %*% t(L),
L %*% coef(m3)))## [,1]
## [1,] 30.51573# Teste de Wald para interação (poderia ser LRT, claro).
# É necessário um objeto glm, mas necesse caso ele não usado para nada.
linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
coef. = coef(m3))## Linear hypothesis test
##
## Hypothesis:
## estbotão floral:I(des^2) = 0
## estflorecimento:I(des^2) = 0
## estmaça:I(des^2) = 0
## estcapulho:I(des^2) = 0
##
## Model 1: restricted model
## Model 2: ncap ~ est * (des + I(des^2))
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 114
## 2 110 4 30.516 3.843e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Predição com bandas de confiança.
pred <- with(capdesfo, expand.grid(est = levels(est),
des = seq(0, 1, by = 0.025)))
X <- model.matrix(formula(m0)[-2], data = pred)
pred <- list(pois = pred, quasi = pred, gcnt = pred)
# Preditos pela Poisson.
aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
# Preditos pela Quasi-Poisson.
aux <- confint(glht(m1, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))
# Comparando as estimativas de média para contagem.
cbind(Pois = pred$pois$fit,
Gcnt1 = c(exp(predict(m3, newdata = X))),
GCnt2 = c(predict(m3, newdata = X, type = "response")))[1:10, ]## Pois Gcnt1 GCnt2
## [1,] 9.154471 9.564302 9.154646
## [2,] 8.450409 8.859703 8.450047
## [3,] 8.908571 9.325240 8.915585
## [4,] 7.890265 8.316089 7.906434
## [5,] 10.248743 10.659755 10.250099
## [6,] 9.230337 9.639172 9.229517
## [7,] 8.550098 8.959491 8.549835
## [8,] 8.641685 9.059457 8.649801
## [9,] 8.046113 8.465629 8.055973
## [10,] 10.130548 10.541938 10.132282# Preditos pela Gamma-Count.
aux <- predict(m3, newdata = X,
interval = "confidence", type = "response")
pred$gcnt <- cbind(pred$gcnt, aux[, c(2, 1, 3)])
pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, est, des, modelo)
key <- list(lines = list(lty = 1),
text = list(c("Poisson", "Quasi-Poisson",
"Gamma-Count")))
key$lines$col <-
trellis.par.get("superpose.line")$col[1:nlevels(pred$modelo)]
p2 <- xyplot(fit ~ des | est, data = pred, groups = modelo,
layout = c(NA, 1), as.table = TRUE,
xlim = extendrange(range(pred$des), f = 0.1),
type = "l", key = key,
ly = pred$lwr, uy = pred$upr,
cty = "bands", alpha = 0.35,
prepanel = prepanel.cbH,
panel.groups = panel.cbH,
panel = panel.superpose)
# p2update(p1, type = "p", layout = c(NA, 1),
key = key, spread = 0.07) +
as.layer(p2, under = TRUE)
Número de Nematoides em Linhagens de Feijão
#-----------------------------------------------------------------------
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 ...# help(nematoide, package = "MRDCr", help_type = "html")
# Número de nematóides por grama de raíz.
plot(nema ~ off, data = nematoide)
xyplot(nema/off ~ cult, data = nematoide,
xlab = "Linhagens de feijão",
ylab = "Número de nematoides por grama de raíz")
#-----------------------------------------------------------------------
# Ajuste do Poisson.
m0 <- glm(nema ~ offset(log(off)) + cult,
data = nematoide,
family = poisson)
m1 <- update(m0, family = quasipoisson)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Estimativas dos parâmetros.
summary(m1)##
## Call:
## glm(formula = nema ~ offset(log(off)) + cult, family = quasipoisson,
## data = nematoide)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.9455 -0.6693 0.1100 1.1519 4.7603
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.9293 0.4496 6.515 7.42e-09 ***
## cultB 0.1662 0.7657 0.217 0.828703
## cultC 0.3417 0.6769 0.505 0.615154
## cultD 0.9128 0.6359 1.436 0.155278
## cultE 0.9692 0.6076 1.595 0.114879
## cultF 0.3973 0.6903 0.576 0.566621
## cultG -0.4507 0.7931 -0.568 0.571522
## cultH 0.3218 0.8665 0.371 0.711433
## cultI 0.8808 0.6018 1.464 0.147507
## cultJ 1.1982 0.6359 1.884 0.063390 .
## cultK 1.4511 0.5965 2.433 0.017369 *
## cultL 1.4299 0.5362 2.667 0.009381 **
## cultM 1.6138 0.5507 2.931 0.004482 **
## cultN 1.7743 0.4958 3.579 0.000610 ***
## cultO 1.5776 0.5268 2.995 0.003718 **
## cultP 1.6719 0.5103 3.277 0.001593 **
## cultQ 2.2105 0.4913 4.499 2.45e-05 ***
## cultR 2.2080 0.5103 4.327 4.60e-05 ***
## cultS 1.9155 0.4958 3.863 0.000236 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.841169)
##
## Null deviance: 556.73 on 93 degrees of freedom
## Residual deviance: 212.54 on 75 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6# Quadro de deviance.
# anova(m0, test = "Chisq")
anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: nema
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 93 556.73
## cult 18 344.19 75 212.54 4.978 3.475e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Ajuste da Poisson Generalizada.
L <- with(nematoide,
list(y = nema, offset = off, X = model.matrix(m0)))
start <- c(alpha = log(1), coef(m0))
parnames(llgcnt) <- names(start)
# Modelo Poisson também.
m2 <- mle2(llgcnt, start = start, data = L,
fixed = list(alpha = 0), vecpar = TRUE)
c(logLik(m2), logLik(m0))## [1] -274.5007 -274.5007# Poisson Generalizada.
m3 <- gcnt(formula(m0), data = nematoide)
# Diferença de deviance.
# 2 * diff(c(logLik(m0), logLik(m3)))
anova(m3, m2)## Likelihood Ratio Tests
## Model 1: m3, [llgcnt]: alpha+(Intercept)+cultB+cultC+cultD+
## cultE+cultF+cultG+cultH+cultI+cultJ+cultK+cultL+
## cultM+cultN+cultO+cultP+cultQ+cultR+cultS
## Model 2: m2, [llgcnt]: (Intercept)+cultB+cultC+cultD+cultE+
## cultF+cultG+cultH+cultI+cultJ+cultK+cultL+cultM+
## cultN+cultO+cultP+cultQ+cultR+cultS
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 20 488.72
## 2 19 549.00 60.278 1 8.236e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0## PoissonGLM PoissonML GCount
## NA 0.0000000 -1.14112640
## (Intercept) 2.9293041 2.9293041 2.64609924
## cultB 0.1662485 0.1662485 -0.04559607
## cultC 0.3417426 0.3417426 0.35670405
## cultD 0.9128492 0.9128492 0.91203676
## cultE 0.9692234 0.9692234 1.01528664
## cultF 0.3973326 0.3973326 0.31963344
## cultG -0.4507108 -0.4507108 -0.63094615
## cultH 0.3217814 0.3217814 -0.11142064
## cultI 0.8808018 0.8808018 0.94473617
## cultJ 1.1982163 1.1982163 1.19639146
## cultK 1.4511238 1.4511238 1.51617919
## cultL 1.4298744 1.4298744 1.59550343
## cultM 1.6137977 1.6137977 1.75476839
## cultN 1.7743062 1.7743062 1.99708212
## cultO 1.5776191 1.5776191 1.75826104
## cultP 1.6718841 1.6718841 1.87564022
## cultQ 2.2104990 2.2104990 2.43923574
## cultR 2.2080215 2.2080215 2.41168192
## cultS 1.9154798 1.9154798 2.14065393splom(c0[-(1:2), ]) + layer(panel.abline(a = 0, b = 1, lty = 2))
# Perfil de log-verossimilhança para o parâmetro de dispersão.
plot(profile(m3, which = "alpha"))
# Covariância.
V <- cov2cor(vcov(m3))
corrplot.mixed(V, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
dev.off()## null device
## 1# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))## sum mean max
## 1.5793846 0.0831255 0.1879248#-----------------------------------------------------------------------
# Predição com intervalos de confiança.
pred <- unique(subset(nematoide, select = cult))
X <- model.matrix(~cult, data = pred)
pred <- list(pois = pred, quasi = pred, gcnt1 = pred, gcnt2 = pred)
# Preditos pela Poisson.
aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
# Preditos pela Quasi-Poisson.
aux <- confint(glht(m1, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))
# Preditos pela Gamma-Count.
aux <- predict(m3, newdata = X,
interval = "confidence", type = "link")
pred$gcnt1 <- cbind(pred$gcnt1, exp(aux[, c(2, 1, 3)]))
aux <- predict(m3, newdata = X,
interval = "confidence", type = "response")
pred$gcnt2 <- cbind(pred$gcnt2, aux[, c(2, 1, 3)])
pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, cult, modelo)
key <- list(type = "o", divide = 1,
lines = list(pch = 1:nlevels(pred$modelo),
lty = 1, col = 1),
text = list(c("Poisson", "Quasi-Poisson",
"Gamma-Count 1", "Gamma-Count 2")))
xyplot(nema/off ~ cult, data = nematoide,
key = key,
xlab = "Cultivar de feijão",
ylab = "Número de nematóides por grama de raíz") +
as.layer(
xyplot(fit ~ cult, data = pred,
pch = pred$modelo,
ly = pred$lwr, uy = pred$upr,
cty = "bars", length = 0,
prepanel = prepanel.cbH,
desloc = 0.15 * scale(as.integer(pred$modelo),
scale = FALSE),
panel = panel.cbH))
Número de Gols por Partida de Futebol
Análise do número de gols feitos pelos times mandantes e desafiantes no Campeonato Brasileiro de 2010.
Para essa análise assume-se que o número de gols feito pelo time mandante (\(y_h\)) e pelo time desafiante (\(y_a\)) na mesma partida são variáveis de contagem independentes. O preditor correspondente ao número de gols do time mandante é \[ \eta_h = \gamma_h - \delta_a + \omega, \] em que \(\gamma_h\) é a força de ataque do time mandante, \(\delta_a\) é a força de defesa do time desafiante e \(\omega\) é vantagem pro time mandante devido a jogar em casa. Para o número de gols do time desafiante é análogo \[ \eta_a = \gamma_a - \delta_h, \] exceto pelo fato de não adicionarmos parâmetro devido ao mando de campo.
Considerando a independência entre os gols dos times que se enfrentam, temos então que \[ \Pr(Y_h = y_h \cap Y_a = y_a) = \Pr(Y_h = y_h) \times \Pr(Y_a = y_a). \]
Dessa forma, a verossimilhança de uma partida é
\[ L(y_h, y_a; \gamma_h, \gamma_a, \delta_h, \delta_a, \omega) = \Pr(y_h | \eta_h = \gamma_h - \delta_a + \omega) \times \Pr(y_a | \eta_a = \gamma_a - \delta_h). \]
Ao serem considerados 20 times, serão 41 parâmetros estimados presentes nos nestes preditores, além dos eventuais parâmetros complementares da distribuição que determina \(\Pr(Y)\), como parâmetros de dispersão.
#-----------------------------------------------------------------------
# Log-verossimilhança para o número de gols dos times em uma partida.
llgols <- function(theta, gh, ga, Xh, Xa){
# theta: Vetor de parâmetros.
# gh, ga: Gols do mandante e do visitante.
# Xh, Xa: Matrizes de incidência das partidas.
gamma <- 1:20 # Parâmetros de ataque.
delta <- 21:40 # Parâmetros de defesa.
omega <- 41 # Vantagem do mando de campo.
alpha <- 42 # Dispersão da Gamma-Count. Se 0 então Poisson.
#----------------------------------------
# Preditor dos times mandantes e desafiante.
eta.h <- (Xh %*% theta[gamma] - Xa %*% theta[delta]) + theta[omega]
eta.a <- (Xa %*% theta[gamma] - Xh %*% theta[delta])
# Parâmetro de dispersão.
alpha <- exp(theta[alpha])
#----------------------------------------
# Quantidades auxiliares.
alpha.gh <- alpha * gh
alpha.eXb.h <- alpha * exp(eta.h)
alpha.ga <- alpha * ga
alpha.eXb.a <- alpha * exp(eta.a)
offset <- 1
#-------------------------------------------------------------------
# Verossimilhanças.
ll.h <- sum(log(pgamma(offset,
shape = alpha.gh,
rate = alpha.eXb.h) -
pgamma(offset,
shape = alpha.gh + alpha,
rate = alpha.eXb.h)))
ll.a <- sum(log(pgamma(offset,
shape = alpha.ga,
rate = alpha.eXb.a) -
pgamma(offset,
shape = alpha.ga + alpha,
rate = alpha.eXb.a)))
# Verossimilhança total.
ll <- sum(ll.h + ll.a)
return(-ll)
}
#-----------------------------------------------------------------------
# Análise dos jogos do Campeonado Brasileiro.
# Tomando as primeiras rodadas do campeonato.
db <- subset(cambras, rod <= 10)
# Quantos jogos cada time fez em casa e fora de casa.
addmargins(cbind(home = xtabs(~home, db),
away = xtabs(~away, db)), margin = 2)## home away Sum
## Atlético/GO 6 4 10
## Atlético/MG 6 4 10
## Atlético/PR 5 5 10
## Avaí 5 5 10
## Botafogo 5 5 10
## Ceará 5 5 10
## Corinthians 5 5 10
## Cruzeiro 4 6 10
## Flamengo 6 4 10
## Fluminense 5 5 10
## Goiás 4 6 10
## Grêmio 5 5 10
## Grêmio Prudente 4 6 10
## Guarani 6 4 10
## Internacional 5 5 10
## Palmeiras 6 4 10
## Santos 4 6 10
## São Paulo 5 5 10
## Vasco 4 6 10
## Vitória 5 5 10subset(db, home == "Fluminense")## rod home away h a
## 2545 2 Fluminense Atlético/GO 1 0
## 2567 4 Fluminense Flamengo 2 1
## 2586 6 Fluminense Vitória 2 1
## 2614 8 Fluminense Grêmio Prudente 1 1
## 2632 10 Fluminense Cruzeiro 1 0# Níveis na mesma ordem?
all(levels(db$home) == levels(db$away))## [1] TRUE# Matrizes de delineamento.
Xh <- model.matrix(~ -1 + home, db) # h: home.
Xa <- model.matrix(~ -1 + away, db) # a: away.
# Verificação.
Xha <- Xh - Xa
db[1:5, c("home", "away")]## home away
## 2535 Botafogo Santos
## 2536 Atlético/GO Grêmio
## 2537 Palmeiras Vitória
## 2538 Flamengo São Paulo
## 2539 Atlético/MG Vascoprint.table(t(Xha[1:5, ]), zero.print = ".")## 2535 2536 2537 2538 2539
## homeAtlético/GO . 1 . . .
## homeAtlético/MG . . . . 1
## homeAtlético/PR . . . . .
## homeAvaí . . . . .
## homeBotafogo 1 . . . .
## homeCeará . . . . .
## homeCorinthians . . . . .
## homeCruzeiro . . . . .
## homeFlamengo . . . 1 .
## homeFluminense . . . . .
## homeGoiás . . . . .
## homeGrêmio . -1 . . .
## homeGrêmio Prudente . . . . .
## homeGuarani . . . . .
## homeInternacional . . . . .
## homePalmeiras . . 1 . .
## homeSantos -1 . . . .
## homeSão Paulo . . . -1 .
## homeVasco . . . . -1
## homeVitória . . -1 . .#-----------------------------------------------------------------------
# Ajuste do modelos.
# Valores iniciais para o modelo Poisson.
start <- c(rep(1, 40), 0.5, 0)
names(start) <- c(paste0("att", 1:20),
paste0("def", 1:20),
"hfa", "alpha")
parnames(llgols) <- names(start)
# Ajuste da Poisson, pois log(alpha) fixo em 0.
m0 <- mle2(minuslogl = llgols,
start = start,
data = list(gh = db$h, ga = db$a, Xh = Xh, Xa = Xa),
fixed = list(alpha = 0),
vecpar = TRUE)
# Valores iniciais para o modelo Gamma-Count.
start <- coef(m0)
parnames(llgols) <- names(start)
# Ajuste da Gamma-Count, alpha será estimado.
m2 <- mle2(minuslogl = llgols,
start = start,
data = list(gh = db$h, ga = db$a, Xh = Xh, Xa = Xa),
vecpar = TRUE)
#-----------------------------------------------------------------------
# Comparando resultados.
# Estimativas dos parâmetros.
c0 <- data.frame(Pois = coef(m0),
GCnt = coef(m2))
xyplot(GCnt ~ Pois, data = c0, aspect = "iso",
groups = gsub(x = rownames(c0), "\\d", ""),
auto.key = TRUE, grid = TRUE) +
layer(panel.abline(a = 0, b = 1))
# Teste para o parâmetro de dispersão.
anova(m0, m2)## Likelihood Ratio Tests
## Model 1: m0, [llgols]: att1+att2+att3+att4+att5+att6+att7+
## att8+att9+att10+att11+att12+att13+att14+att15+
## att16+att17+att18+att19+att20+def1+def2+def3+
## def4+def5+def6+def7+def8+def9+def10+def11+def12+
## def13+def14+def15+def16+def17+def18+def19+def20+
## hfa
## Model 2: m2, [llgols]: att1+att2+att3+att4+att5+att6+att7+
## att8+att9+att10+att11+att12+att13+att14+att15+
## att16+att17+att18+att19+att20+def1+def2+def3+
## def4+def5+def6+def7+def8+def9+def10+def11+def12+
## def13+def14+def15+def16+def17+def18+def19+def20+
## hfa+alpha
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 41 522.30
## 2 42 509.64 12.654 1 0.0003748 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Função de risco.
al <- exp(coef(m2)[42])
curve(dgamma(x, al, 1)/(1 - pgamma(x, al, 1)),
from = 0, to = 2,
xlab = "t",
ylab = expression(h(t) == f(t)/(1 - F(t))))
#-----------------------------------------------------------------------
# Estimativas intervalares.
ci <- rbind(cbind(1, coef(m0)[-42], confint(m0, method = "quad")),
cbind(2, coef(m2), confint(m2, method = "quad")))
ci <- as.data.frame(ci)
names(ci) <- c("model", "est", "lwr", "upr")
ci$model <- factor(ci$model, labels = c("Pois", "GCnt"))
ci$par <- factor(sub(pattern = "\\d+",
replacement = "",
x = rownames(ci)))
ci$team <- factor(levels(db$home)[
as.numeric(sub(pattern = "\\D+",
replacement = "",
x = rownames(ci)))])
sb <- subset(ci, par == "att" & model == "GCnt")
ci$team <- factor(ci$team, levels = sb$team[order(sb$est)])
segplot(team ~ lwr + upr | model,
centers = est, data = ci,
draw = FALSE, groups = par, gap = 0.2,
pch = as.integer(ci$par),
key = list(title = "Parâmetro", cex.title = 1.1,
type = "o", divide = 1,
text = list(c("Ataque", "Defesa")),
lines = list(pch = 2:3)),
ylab = "Times (ordenados pelo ataque)",
xlab = "Estimativa do parâmetro",
panel = panel.groups.segplot)
ci <- arrange(ci, par, team, model)
segplot(team ~ lwr + upr | par,
centers = est, data = subset(ci, par %in% c("att", "def")),
draw = FALSE, groups = model, gap = 0.2,
pch = as.integer(ci$model),
key = list(title = "Modelo", cex.title = 1.1,
type = "o", divide = 1,
text = list(c("Poisson", "Gamma-Count")),
lines = list(pch = 2:1)),
ylab = "Times (ordenados pelo ataque)",
xlab = "Estimativa do parâmetro",
panel = panel.groups.segplot)
#-----------------------------------------------------------------------
# Gols observados e preditos dos times mandantes e desafiantes nestas
# rodadas.
gg <-
rbind(
cbind(x = 1,
y = db$h,
Pois = c(exp(cbind(Xh, -Xa) %*%
coef(m0)[1:40] + coef(m0)["hfa"])),
GCnt1 = c(exp(cbind(Xh, -Xa) %*%
coef(m2)[1:40] + coef(m2)["hfa"])),
GCnt2 = calc_mean_gcnt(lambda = exp(cbind(Xh, -Xa) %*%
coef(m2)[1:40] +
coef(m2)["hfa"]),
alpha = exp(coef(m2)["alpha"]))),
cbind(x = 2,
y = db$a,
Pois = c(exp(cbind(Xa, -Xh) %*% coef(m0)[1:40])),
GCnt1 = c(exp(cbind(Xa, -Xh) %*% coef(m2)[1:40])),
GCnt2 = calc_mean_gcnt(lambda = exp(cbind(Xa, -Xh) %*%
coef(m2)[1:40]),
alpha = exp(coef(m2)["alpha"]))))
# Comparação de observado com predito.
splom(gg[, -1], groups = gg[, 1]) +
layer(panel.abline(a = 0, b = 1))
#-----------------------------------------------------------------------
# Fazendo previsão dos resultados da rodada a seguir.
levels(db$home)## [1] "Atlético/GO" "Atlético/MG" "Atlético/PR"
## [4] "Avaí" "Botafogo" "Ceará"
## [7] "Corinthians" "Cruzeiro" "Flamengo"
## [10] "Fluminense" "Goiás" "Grêmio"
## [13] "Grêmio Prudente" "Guarani" "Internacional"
## [16] "Palmeiras" "Santos" "São Paulo"
## [19] "Vasco" "Vitória"# # Na rodada 11 deu Cruzeiro 2x2 Grêmio.
# i <- which(levels(db$home) %in% c("Cruzeiro", "Grêmio"))
# levels(db$home)[i]
# Na rodada 11 dei Ceará 0x0 Palmeiras.
i <- which(levels(db$home) %in% c("Ceará", "Palmeiras"))
levels(db$home)[i]## [1] "Ceará" "Palmeiras"# Dois times jogando em território neutro.
# HomeAttack - AwayDefense
# AwayAttack - HomeDefense
alp <- exp(coef(m2)[42])
beta <- exp(coef(m2)[i] - rev(coef(m2)[20 + i]))
names(beta) <- levels(db$home)[i]
beta## Ceará Palmeiras
## 0.7866764 0.3602727# Tempo médio entre Gols.
exp(-beta)## Ceará Palmeiras
## 0.4553557 0.6974861# Probabilidade dos placares de 0x0, 0x1, ..., 6x6.
y <- 0:6
dnn <- lapply(substr(levels(db$home)[i], 0, 3), FUN = paste, y)
plapois <- outer(dgcnt(y = y, lambda = beta[1], alpha = 1),
dgcnt(y = y, lambda = beta[2], alpha = 1),
FUN = "*")
plagcnt <- outer(dgcnt(y = y, lambda = beta[1], alpha = alp),
dgcnt(y = y, lambda = beta[2], alpha = alp),
FUN = "*")
dimnames(plapois) <- dnn
dimnames(plagcnt) <- dnn
print.table(round(plapois, 2), zero.print = ".")## Pal 0 Pal 1 Pal 2 Pal 3 Pal 4 Pal 5 Pal 6
## Cea 0 0.32 0.11 0.02 . . . .
## Cea 1 0.25 0.09 0.02 . . . .
## Cea 2 0.10 0.04 0.01 . . . .
## Cea 3 0.03 0.01 . . . . .
## Cea 4 0.01 . . . . . .
## Cea 5 . . . . . . .
## Cea 6 . . . . . . .print.table(round(plagcnt, 2), zero.print = ".")## Pal 0 Pal 1 Pal 2 Pal 3 Pal 4 Pal 5 Pal 6
## Cea 0 0.43 0.09 . . . . .
## Cea 1 0.32 0.07 . . . . .
## Cea 2 0.07 0.01 . . . . .
## Cea 3 0.01 . . . . . .
## Cea 4 . . . . . . .
## Cea 5 . . . . . . .
## Cea 6 . . . . . . .# Olhando só para os empates.
round(cbind(Pois = diag(plapois),
GCnt = diag(plagcnt)), 3)## Pois GCnt
## [1,] 0.318 0.428
## [2,] 0.090 0.067
## [3,] 0.006 0.001
## [4,] 0.000 0.000
## [5,] 0.000 0.000
## [6,] 0.000 0.000
## [7,] 0.000 0.000# Draw - Win - Lose.
dwl <- rbind(Pois = c(sum(diag(plapois)),
sum(plapois[lower.tri(plapois)]),
sum(plapois[upper.tri(plapois)])),
GCnt = c(sum(diag(plagcnt)),
sum(plagcnt[lower.tri(plagcnt)]),
sum(plagcnt[upper.tri(plagcnt)])))
colnames(dwl) <- sprintf(paste(levels(db$home)[i], collapse = " %s "),
c("=", ">", "<"))
t(dwl)## Pois GCnt
## Ceará = Palmeiras 0.4142014 0.49552479
## Ceará > Palmeiras 0.4288175 0.40618126
## Ceará < Palmeiras 0.1569623 0.09829394Informaçõ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] methods stats4 grid stats graphics grDevices
## [7] utils datasets base
##
## other attached packages:
## [1] MRDCr_0.0-2 multcomp_1.4-6 TH.data_1.0-7
## [4] MASS_7.3-45 survival_2.39-5 mvtnorm_1.0-5
## [7] doBy_4.5-15 car_2.1-2 plyr_1.8.4
## [10] corrplot_0.77 bbmle_1.0.18 latticeExtra_0.6-28
## [13] RColorBrewer_1.1-2 lattice_0.20-33 rmarkdown_1.0
## [16] 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 lme4_1.1-12
## [7] evaluate_0.9 nlme_3.1-128 mgcv_1.8-12
## [10] Matrix_1.2-6 yaml_2.1.13 parallel_3.3.1
## [13] SparseM_1.7 stringr_1.0.0 MatrixModels_0.4-1
## [16] nnet_7.3-12 minqa_1.2.4 magrittr_1.5
## [19] codetools_0.2-14 htmltools_0.3.5 splines_3.3.1
## [22] pbkrtest_0.4-6 numDeriv_2014.2-1 quantreg_5.26
## [25] sandwich_2.3-4 stringi_1.1.1 zoo_1.7-13
Modelos de Regressão para Dados de Contagem com o R |
|
| 61 RBRAS | Departamento de Estatística - UFPR |
| 23 a 25 de Maio de 2016 | LEG, PET Estatística |
| Salvador - Bahia | github.com/leg-ufpr/MRDCr |