Skip to content
Snippets Groups Projects
Commit 1cc01f08 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Termina a vinheta da Gamma COunt.

parent ee11d675
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
inst/doc/v06_gamma_count.R 0 → 100644
+ 1286
− 0
View file @ 1cc01f08
## ----setup, include=FALSE-----------------------------------------
source("_setup.R")
## ---- message=FALSE, error=FALSE, warning=FALSE-------------------
# Definições da sessão.
library(lattice)
library(latticeExtra)
library(grid)
library(rpanel)
library(bbmle)
library(corrplot)
library(plyr)
library(car)
library(doBy)
library(multcomp)
library(MRDCr)
## -----------------------------------------------------------------
# Função densidade na parametrização original.
dgcnt
# 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",
panel = function(x, y, ...) {
m <- sum(x * y)
panel.xyplot(x, y, ...)
panel.abline(v = m, lty = 2)
# grid.text(label = sprintf("%0.1f", m),
# x = unit(0.95, "npc"),
# y = unit(0.9, "npc"),
# hjust = 1)
}),
strip = strip.custom(
strip.names = TRUE,
var.name = expression(lambda == ""),
sep = ""),
strip.left = strip.custom(
strip.names = TRUE,
var.name = expression(alpha == ""),
sep = ""))
## ---- eval=FALSE--------------------------------------------------
# # Função densidade na parametrização de modelo de regressão.
# MRDCr::dgcnt
#
# react <- function(panel){
# with(panel,
# {
# m <- LAMBDA
# s <- sqrt(LAMBDA/exp(ALPHA))
# y <- 0:max(c(YMAX, ceiling(m + 5 * s)))
# py <- dgcnt(y = y, lambda = LAMBDA, alpha = exp(ALPHA))
# if (POIS) {
# pz <- dpois(y, lambda = m)
# } else {
# pz <- 0
# }
# py[!is.finite(py)] <- 0
# plot(py ~ y, type = "h",
# ylim = c(0, max(c(py, pz))),
# xlab = expression(y),
# ylab = expression(f(y)))
# mtext(side = 3,
# text = substitute(sum(f(y)) == s,
# list(s = round(sum(py), 5))))
# if (EX) {
# m <- sum(y * py)
# # v <- sum((y - m)^2 * py)
# legend("topright", bty = "n",
# legend = substitute(E(Y) == mu,
# list(mu = m)))
# abline(v = m, col = 2)
# }
# if (POIS) {
# lines(y + 0.3, pz, type = "h", col = 3)
# }
# })
# panel
# }
#
# panel <- rp.control(title = "Gamma Count",
# size = c(300, 100), YMAX = 50)
# rp.slider(panel = panel, action = react,
# variable = LAMBDA, title = "lambda",
# from = 1, to = 0.75 * panel$YMAX,
# initval = 5, resolution = 0.1,
# showvalue = TRUE)
# rp.slider(panel = panel, action = react,
# variable = ALPHA, title = "alpha",
# from = -5, to = 5,
# initval = 0, resolution = 0.02,
# showvalue = TRUE)
# rp.checkbox(panel = panel,
# variable = EX, action = react, title = "E(Y)",
# labels = "Mostrar o valor esperado?")
# rp.checkbox(panel = panel,
# variable = POIS, action = react, title = "Poisson",
# labels = "Adicionar a distribução Poisson?")
# rp.do(panel = panel, action = react)
## -----------------------------------------------------------------
#-----------------------------------------------------------------------
# Função de log-Verossimilhança da Poisson Generalizada na
# parametrização de modelo de regressão.
MRDCr::llgcnt
#-----------------------------------------------------------------------
# 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)))
# Estimando o \alpha.
n1 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
coef(n1)
# 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()
# TODO documentar no pacote.
fitted.gcnt <- function(lambda, alpha, offset = NULL, ymax = 40) {
if (is.null(offset)) {
offset <- 1
}
pars <- cbind(lambda, alpha, offset)
y <- 1:ymax
py <- apply(pars, MARGIN = 1,
FUN = function(p) {
sum(pgamma(p[3],
shape = p[2] * y,
rate = p[2] * p[1]))
# sum(y * dgcnt(y = y, lambda = p[1], alpha = p[2]))
})
names(py) <- NULL
return(py)
}
mean(y)
fitted.gcnt(lambda = exp(coef(n1)[2]),
alpha = exp(coef(n1)[1]),
offset = 2)
#-----------------------------------------------------------------------
# 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]),
fitted.gcnt(lambda = exp(coef(n1)[2]), alpha = exp(coef(n1)[1])))
y <- 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]),
fitted.gcnt(lambda = exp(coef(n1)[2]),
alpha = exp(coef(n1)[1]), ymax = 100))
## -----------------------------------------------------------------
#-----------------------------------------------------------------------
# Carregando e explorando os dados.
data(soja, package = "MRDCr")
str(soja)
# 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")
summary(m1)
#-----------------------------------------------------------------------
# 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))
cbind(coef(m2)[-1], coef(m0))
# 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)
# Estimaitvas dos parâmetros.
c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0
splom(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()
## -----------------------------------------------------------------
# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))
#-----------------------------------------------------------------------
# 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)))
# 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))
#-----------------------------------------------------------------------
# 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)
# Quantil normal.
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))
# 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))
# Matrix de covariância completa e sem o alpha (marginal).
V <- vcov(m3)
V <- V[-1, -1]
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
pred$gcnt$eta <- c(X %*% coef(m3)[-1])
pred$gcnt <- cbind(pred$gcnt,
apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
FUN = function(x) {
# exp(pred$gcnt$eta + x)
fitted.gcnt(
lambda = exp(pred$gcnt$eta + x),
alpha = exp(coef(m3)["alpha"]),
ymax = 300)
}))
pred$gcnt$eta <- NULL
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)
## -----------------------------------------------------------------
#-----------------------------------------------------------------------
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))
#-----------------------------------------------------------------------
# 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")
summary(m1)
#-----------------------------------------------------------------------
# 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))
# 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)
# Estimaitvas dos parâmetros.
c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0
splom(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()
## -----------------------------------------------------------------
# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))
# 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)))
linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
coef. = coef(m3))
#-----------------------------------------------------------------------
# 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)
# Quantil normal.
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))
# 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))
# Matrix de covariância completa e sem o alpha (marginal).
V <- vcov(m3)
V <- V[-1, -1]
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
pred$gcnt$eta <- c(X %*% coef(m3)[-1])
pred$gcnt <- cbind(pred$gcnt,
apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
FUN = function(x) {
exp(pred$gcnt$eta + x)
# fitted.gcnt(
# lambda = exp(pred$gcnt$eta + x),
# alpha = exp(coef(m3)["alpha"]),
# ymax = 100)
}))
pred$gcnt$eta <- NULL
# Junta o resultado dos 3 modelos.
pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, umid, K, modelo)
str(pred)
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 (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")
summary(m1)
# 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")
#-----------------------------------------------------------------------
# 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))
# 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)
c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0
splom(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()
## -----------------------------------------------------------------
# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))
# 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)))
linearHypothesis(model = m0,
hypothesis.matrix = L,
vcov. = vcov(m3),
coef. = coef(m3))
#-----------------------------------------------------------------------
# X <- model.matrix(m0)
# soja$nvag
#
# cbind(y = soja$ngra,
# Pois = soja$nvag * exp(X %*% coef(m0)),
# Quas = soja$nvag * exp(X %*% coef(m1)),
# Gcnt = exp(X %*% coef(m3)[-1]))
#
# cbind(y = soja$ngra/soja$nvag,
# Pois = exp(X %*% coef(m0)),
# Quas = exp(X %*% coef(m1)),
# Gcnt = exp(X %*% coef(m3)[-1]))
#
# fitted.gcnt(lambda = 1, alpha = 1)
# fitted.gcnt(lambda = 2, alpha = 1)
# fitted.gcnt(lambda = 2, alpha = 0.4)
# fitted.gcnt(lambda = 2, alpha = 3.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)
pred <- list(pois = pred, quasi = pred, gcnt = pred)
# Quantil normal.
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))
# 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))
# # TODO
# m3 <- m2
# Matrix de covariância completa e sem o alpha (marginal).
V <- vcov(m3)
V <- V[-1, -1]
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
pred$gcnt$eta <- c(X %*% coef(m3)[-1])
# Comparando as estimativas de média para contagem.
c(Pois = pred$pois$fit[1],
Gcnt1 = exp(pred$gcnt$eta)[1],
GCnt2 = fitted.gcnt(
lambda = exp(pred$gcnt$eta)[1],
alpha = exp(coef(m3)["alpha"]),
ymax = 80))
pred$gcnt <- cbind(pred$gcnt,
apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
FUN = function(x) {
exp(pred$gcnt$eta + x)
# fitted.gcnt(
# lambda = exp(pred$gcnt$eta + x),
# alpha = exp(coef(m3)["alpha"]),
# ymax = 50)
}))
pred$gcnt$eta <- NULL
# 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 em função do nível de desfolha artificial e fase
# fenológica do algodoeiro.
data(capdesfo, package = "MRDCr")
str(capdesfo)
# 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")
summary(m1)
#-----------------------------------------------------------------------
# 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))
# Modelo Poisson Generalizado.
m3 <- mle2(llgcnt, start = start, data = L, vecpar = TRUE)
logLik(m3)
anova(m3, m2)
summary(m3)
c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0
splom(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()
## -----------------------------------------------------------------
# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))
# 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)))
# 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))
#-----------------------------------------------------------------------
# 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)
# Quantil normal.
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))
# 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 Poisson Generalizada.
V <- vcov(m3)
V <- V[-1, -1]
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
pred$gcnt$eta <- c(X %*% coef(m3)[-1])
# Comparando as estimativas de média para contagem.
c(Pois = pred$pois$fit[1],
Gcnt1 = exp(pred$gcnt$eta)[1],
GCnt2 = fitted.gcnt(
lambda = exp(pred$gcnt$eta)[1],
alpha = exp(coef(m3)["alpha"]),
ymax = 30))
pred$gcnt <- cbind(pred$gcnt,
apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
FUN = function(x) {
# exp(pred$gcnt$eta + x)
fitted.gcnt(
lambda = exp(pred$gcnt$eta + x),
alpha = exp(coef(m3)["alpha"]),
ymax = 50)
}))
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)
# p2
## ---- fig.width=7, fig.height=3.5---------------------------------
update(p1, type = "p", layout = c(NA, 1),
key = key, spread = 0.07) +
as.layer(p2, under = TRUE)
## -----------------------------------------------------------------
#-----------------------------------------------------------------------
data(nematoide, package = "MRDCr")
str(nematoide)
# 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)
# Quadro de deviance.
# anova(m0, test = "Chisq")
anova(m1, test = "F")
#-----------------------------------------------------------------------
# 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))
# Poisson Generalizada.
m3 <- gcnt(formula(m0), data = nematoide)
# Diferença de deviance.
# 2 * diff(c(logLik(m0), logLik(m3)))
anova(m3, m2)
c0 <- cbind("PoissonGLM" = c(NA, coef(m0)),
"PoissonML" = coef(m2),
"GCount" = coef(m3))
c0
splom(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()
## -----------------------------------------------------------------
# Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1]))
# Gráfico das estimativas.
pars <- data.frame(Pois = c(0, coef(m0)), Gcnt = coef(m3))
xyplot(Gcnt ~ Pois, data = pars, aspect = "iso", grid = TRUE) +
layer(panel.abline(a = 0, b = 1, lty = 2))
#-----------------------------------------------------------------------
# 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, gcnt = pred)
# Quantil normal.
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))
# 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 Poisson Generalizada.
V <- vcov(m3)
V <- V[-1, -1]
U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) }))
pred$gcnt$eta <- c(X %*% coef(m3)[-1])
# Comparando as estimativas de média para contagem.
c(Pois = pred$pois$fit[1],
Gcnt1 = exp(pred$gcnt$eta)[1],
GCnt2 = fitted.gcnt(
lambda = exp(pred$gcnt$eta)[1],
alpha = exp(coef(m3)["alpha"]),
ymax = 500))
pred$gcnt <- cbind(pred$gcnt,
apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
FUN = function(x) {
# exp(pred$gcnt$eta + x)
fitted.gcnt(
lambda = exp(pred$gcnt$eta + x),
alpha = exp(coef(m3)["alpha"]),
ymax = 500)
}))
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")))
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.25 * scale(as.integer(pred$modelo),
scale = FALSE),
panel = panel.cbH))
## -----------------------------------------------------------------
#-----------------------------------------------------------------------
# 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)
subset(db, home == "Fluminense")
# Níveis na mesma ordem?
all(levels(db$home) == levels(db$away))
# 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")]
print.table(t(Xha[1:5, ]), zero.print = ".")
#-----------------------------------------------------------------------
# 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)
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)
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))
c0
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)
# plot(profile(m1, which = "alpha"))
# plot(profile(m2, which = "alpha"))
# 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))))
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)
#-----------------------------------------------------------------------
# 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 = fitted.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 = fitted.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)
# # 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]
# 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]))
beta
exp(-beta)
# 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 = ".")
print.table(round(plagcnt, 2), zero.print = ".")
# Olhando só para os empates.
round(cbind(Pois = diag(plapois),
GCnt = diag(plagcnt)), 3)
# 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)
## ---- eval=FALSE, include=FALSE-----------------------------------
# opts_chunk$set(eval = FALSE)
inst/doc/v06_gamma_count.html 0 → 100644
+ 2361
− 0
View file @ 1cc01f08
Source diff could not be displayed: it is too large. Options to address this: view the blob.
vignettes/v06_gamma_count.Rmd
+ 267
− 2
View file @ 1cc01f08
...@@ -97,20 +97,54 @@ $$ ...@@ -97,20 +97,54 @@ $$
Como Como
$$ $$
F_n(T) = G(n\alpha, \beta) = F_n(T) = G(n\alpha, \beta T) =
\frac{1}{\Gamma(n\alpha)} \frac{1}{\Gamma(n\alpha)}
\int_{0}^{T} t^{n\alpha-1}\cdot\exp\{-\beta t\}\, \text{d}t, \int_{0}^{T} t^{n\alpha-1}\cdot\exp\{-\beta t\}\, \text{d}t,
$$ $$
podemos escrever $\Pr(N = n)$ como sendo a diferença de acumuladas da podemos escrever $\Pr(N = n)$ como sendo a diferença de acumuladas da
Gama, Gama,
$$ $$
\Pr(N = n) = G(n\alpha, \beta) - G((n+1)\alpha, \beta). \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 Em resumo, a distribuição para o número de eventos, quando se assume a
distribuição do intervalo entre eventos é Gamma, é resultado da distribuição do intervalo entre eventos é Gamma, é resultado da
diferença de duas acumuladas da distribuição Gamma. 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$.
```{r, message=FALSE, error=FALSE, warning=FALSE} ```{r, message=FALSE, error=FALSE, warning=FALSE}
# Definições da sessão. # Definições da sessão.
library(lattice) library(lattice)
...@@ -1215,6 +1249,237 @@ xyplot(nema/off ~ cult, data = nematoide, ...@@ -1215,6 +1249,237 @@ xyplot(nema/off ~ cult, data = nematoide,
## Número de Gols por Partida de Futebol ## ## 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.
```{r} ```{r}
#-----------------------------------------------------------------------
# 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)
subset(db, home == "Fluminense")
# Níveis na mesma ordem?
all(levels(db$home) == levels(db$away))
# 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")]
print.table(t(Xha[1:5, ]), zero.print = ".")
#-----------------------------------------------------------------------
# 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)
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)
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))
c0
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)
# plot(profile(m1, which = "alpha"))
# plot(profile(m2, which = "alpha"))
# 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))))
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)
#-----------------------------------------------------------------------
# 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 = fitted.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 = fitted.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)
# # 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]
# 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]))
beta
exp(-beta)
# 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 = ".")
print.table(round(plagcnt, 2), zero.print = ".")
# Olhando só para os empates.
round(cbind(Pois = diag(plapois),
GCnt = diag(plagcnt)), 3)
# 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)
```
```{r, eval=FALSE, include=FALSE}
opts_chunk$set(eval = FALSE) opts_chunk$set(eval = FALSE)
``` ```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment