Commit 5b0599f2 authored by Eduardo E. R. Junior's avatar Eduardo E. R. Junior

Merge branch 'walmes' into devel

Conflicts:
	DESCRIPTION
	NAMESPACE
	R/methods.R
	R/panel.beeswarm.R
	inst/doc/v01_poisson.html
	inst/slides/slides-mrdcr.Rnw
	man/panel.beeswarm.Rd
parents 0ab8cc65 706885dd
......@@ -17,3 +17,4 @@ ci-c3sl.R
playground.R
style.css
plano.*
_vignettes/
......@@ -38,6 +38,7 @@ Imports:
Suggests:
MASS,
car,
corrplot,
plyr,
ggplot2,
reshape,
......@@ -45,6 +46,8 @@ Suggests:
knitr,
rmarkdown,
shiny,
corrplot
rpanel,
testthat,
methods
VignetteBuilder: knitr
RoxygenNote: 5.0.1
......@@ -8,8 +8,15 @@ export(cmp)
export(convergencez)
export(dcmp)
export(llcmp)
export(dgcnt)
export(dpgnz)
export(gcnt)
export(llgcnt)
export(llpgnz)
export(panel.beeswarm)
export(panel.cbH)
export(panel.groups.segplot)
export(pgnz)
export(prepanel.cbH)
import(bbmle)
import(doBy)
......@@ -23,5 +30,7 @@ importFrom(stats,model.frame)
importFrom(stats,model.matrix)
importFrom(stats,model.offset)
importFrom(stats,model.response)
importFrom(stats,pgamma)
importFrom(stats,poisson)
importFrom(stats,qnorm)
importFrom(utils,combn)
......@@ -29,6 +29,7 @@
#'
#' @seealso \code{\link[doBy]{LSmatrix}}.
#' @import doBy multcomp
#' @importFrom utils combn
#' @examples
#'
#' X <- diag(3)
......@@ -79,7 +80,7 @@ apc <- function(lfm, lev = NULL) {
lev <- as.character(1:nlev)
}
}
cbn <- combn(x = seq_along(lev), m = 2)
cbn <- utils::combn(x = seq_along(lev), m = 2)
M <- lfm[cbn[1, ], ] - lfm[cbn[2, ], ]
if (is.vector(M)) {
dim(M) <- c(1, length(M))
......
#' @name cambras
#' @title Resultados do Campeonato Brasileiro de 2010
#' @description Dados sobre o número de gols dos times mandantes e
#' desafiantes de todas as partidas do Campeonato Brasileiro de
#' 2010.
#' @format Um \code{data.frame} com 380 observações e 5 variáveis.
#' \describe{
#'
#' \item{\code{rod}}{Inteiro que identifica o número da rodada.}
#'
#' \item{\code{home}}{Fator que identifica o time mandante da partida,
#' aquele que jogou em casa.}
#'
#' \item{\code{away}}{Fator que identifica o time desafiante da partida,
#' aquele que jogou fora de casa.}
#'
#' \item{\code{h}}{Número de gols feitos pelo time mandante na partida.}
#'
#' \item{\code{a}}{Número de gols feitos pelo time desafiante na
#' partida.}
#'
#' }
#' @examples
#'
#' # Resultados finais na internet.
#' # browseURL(paste0("http://esportes.terra.com.br/futebol/",
#' # "brasileiro/2010/noticias/0,,OI4339585-EI15413",
#' # ",00-Classificacao+e+Jogos+Serie+A.html"))
#'
#' Xh <- model.matrix(~-1 + home, cambras)
#' Xa <- model.matrix(~-1 + away, cambras)
#'
#' Xha <- Xh - Xa
#' cambras[1:5, c("home", "away")]
#' print(as.table(t(Xha[1:5, ])), zero.print = ".")
#'
#' gsgc <- cbind(scored = t(Xh) %*% cambras$h + t(Xa) %*% cambras$a,
#' conceded = t(Xa) %*% cambras$h + t(Xh) %*% cambras$a)
#' colnames(gsgc) <- c("gScored", "gConced")
#' gsgc
#'
#' # Pontos em cada partida (empate = 1, vitória = 3).
#' pts <- with(cambras, {
#' d <- h - a
#' draw <- d == 0
#' winH <- d > 0
#' winA <- !(draw | winH)
#' x <- cbind(h = winH * 3 + draw, a = winA * 3 + draw)
#' return(x)
#' })
#'
#' tableIn <- function(x) {
#' tb <- table(x)
#' f <- rep(0, 3)
#' names(f) <- c(0, 1, 3)
#' f[names(tb)] <- tb
#' return(f)
#' }
#'
#' # Derrotas, empates e vitórias.
#' ldw <-
#' do.call(rbind, lapply(tapply(pts[, 1], cambras$home,
#' FUN = tableIn),
#' FUN = as.vector)) +
#' do.call(rbind, lapply(tapply(pts[, 2], cambras$away,
#' FUN = tableIn),
#' FUN = as.vector))
#' colnames(ldw) <- c("Lose", "Draw", "Win")
#' ldw
#'
#' # Tabela completa de 2010.
#' pl10 <- t(Xh) %*% pts[, "h"] + t(Xa) %*% pts[, "a"]
#' pl10 <- cbind(team = levels(cambras$home), data.frame(pts = pl10))
#' pl10 <- cbind(pl10, ldw[, 3:1], gsgc)
#' pl10 <- plyr::arrange(pl10, -pts)
#' pl10$pos <- rank(-pl10$pts)
#' pl10
#'
NULL
......@@ -9,8 +9,16 @@
#' aves por calor. Nesse experimento, o sistema só foi utilizado a
#' partir dos 21 dias de idade. A cada dia foi contado o número de
#' aves encontradas mortas no aviário.
#' @format Um \code{data.frame} com 176 observações e 4 variáveis, em
#' que
#'
#' Fora dos galpões, um sistema de monitoramento das variáveis
#' ambientais registrou, em intervalos de 1 hora dos 21 aos 39 dias
#' de idade, as variáveis para que fossem determinados: a entalpia
#' específica do ar (H), a carga térmica de radiação (CTR) e o
#' índice de temperatura de globo negro e umidade (ITGU). Essas
#' variáveis tem a finalidade de explicar a variação da mortalidade
#' das aves nos sistemas de resfriamento ao longo dos dias.
#' @format \code{confterm} é um \code{data.frame} com 176 observações e
#' 4 variáveis, em que
#'
#' \describe{
#'
......@@ -27,6 +35,26 @@
#' \item{\code{nap}}{Número de aves perdidas (ou mortas) por dia.}
#'
#' }
#'
#' \code{conftemp} é um \code{data.frame} com 456 observações e 6
#' variáveis, em que
#'
#' \describe{
#'
#' \item{\code{hora}}{As horas em cada dia, retomando do 0 em cada novo
#' dia.}
#'
#' \item{\code{hr}}{As horas a partir o primeiro dia continuamente.}
#'
#' \item{\code{idade}}{A idade dos animais, em dias.}
#'
#' \item{\code{h}}{Entalpia específica do ar.}
#'
#' \item{\code{ctr}}{Carga térmica de radiação.}
#'
#' \item{\code{itgu}}{Índice de temperatura de globo negro e umidade.}
#'
#' }
#' @source MACHADO, N. S.; TINÔCO, I. D. F. F.; ZOLNIER, S.; MOGAMI,
#' C. A.; DAMASCENO, F. A.; ZEVIANI, W. M. Resfriamento da cobertura
#' de aviários e seus efeitos na mortalidade e nos índices de
......@@ -35,9 +63,14 @@
#' \url{http://www.ufv.br/dea/ambiagro/gallery/publicações/Artigo5.pdf}.
#' @examples
#'
#' #-----------------------------------------
#' # Gráfico da mortalidade das aves.
#'
#' library(lattice)
#' library(latticeExtra)
#'
#' str(confterm)
#' summary(confterm)
#'
#' xtabs(~idade + resfr, data = confterm)
#'
......@@ -50,4 +83,48 @@
#' "Sem sistema de resfriamento")),
#' auto.key = list(corner = c(0.05, 0.9)))
#'
#' #-----------------------------------------
#' # Gráfico das variáveis térmicas.
#'
#' # Amplitude estendida das variáveis.
#' lim <- with(conftemp, apply(cbind(h, ctr, itgu), MARGIN = 2,
#' FUN = extendrange, f = 0.2))
#'
#' # Anotação da eixo x do gráfico.
#' scales <- list(
#' y = list(relation = "free"),
#' x = list(at = seq(from = 1,
#' to = ceiling(max(conftemp$hr/24)) * 24,
#' by = 24)))
#' scales$x$labels <- seq_along(scales$x$at)
#'
#' xyplot(h + ctr + itgu ~ hr, data = conftemp,
#' outer = TRUE, type = "l", layout = c(1, NA),
#' scales = scales, xlim = range(scales$x$at),
#' xlab = "Dias",
#' ylab = "Variáveis térmicas",
#' panel = function(y, subscripts, ...) {
#' wp <- which.packet()
#' r <- lim[, wp[1]]
#' panel.rect(10.5 + 24 * (scales$x$labels - 1), r[1],
#' 20 + 24 * (scales$x$labels - 1), r[2],
#' col = "blue",
#' border = "transparent",
#' alpha = 0.25)
#' panel.xyplot(y = y, subscripts = subscripts, ...)
#' })
#'
#' # Valores máximos do dia.
#' tempdia <- aggregate(cbind(hm = h, cm = ctr, im = itgu) ~ idade,
#' data = conftemp, FUN = max)
#'
#' splom(tempdia[, -1])
#'
#' confterm <- merge(confterm, tempdia, by = "idade")
#' str(confterm)
#'
NULL
#' @name conftemp
#' @rdname confterm
NULL
This diff is collapsed.
#' @name nematoide
#' @title Número de Nematóides em Raízes de Feijoeiro
#' @description Resultados de um experimento em casa de vegetação que
#' estudou a reprodução de nematóides em cultivares/linhagens de
#' feijoeiro. O solo dos vasos foi inicialmente contaminado com
#' namatóides e as parcelas tiveram duas plantas. Ao final do
#' experimento, as raízes das duas plantas por parcela foram
#' lavadas, trituradas, peneiradas e diluídas para fazer a contagem
#' dos nematóides em aliquotas dessa solução.
#' @format Um \code{data.frame} com 94 observações e 4 variáveis.
#'
#' \describe{
#'
#' \item{\code{cult}}{Fator categórico que indica a linhagem de
#' feijoeiro semeada em vasos com solo contaminado com nematóide.}
#'
#' \item{\code{mfr}}{Massa fresca de raízes (g) produzida por parcela
#' (duas plantas) que foi lavada, triturada, peneirada e diluída
#' para fazer a contagem dos nematóides.}
#'
#' \item{\code{vol}}{Volume (ml) usado para diluir a massa fresca de
#' raízes. Esse volume foi agitado para homogeneização e depois uma
#' alíquota de 1 ml foi extraída e colocada em uma lâmina de
#' contagem.}
#'
#' \item{\code{nema}}{Número de nematóides na alíquota de 1 ml,
#' determinado por contagem direta na lâmina.}
#'
#' \item{\code{off}}{É o offset da contagem, o equivalente em massa de
#' fresca de raízes de uma aliquota de 1 ml, ou seja, é
#' \code{off = mfr/vol}.}
#'
#' }
#' @source Cedido para fins acadêmicos por Andressa Cristina Zamboni
#' Machado, pesquisadora do Instituto Agronômico do Paraná (IAPAR),
#' e pelo técnico agrícola do IAPAR Santino Aleandro da Silva.
#'
#' O nome das cultivares, a espécie do nematóide e outras informações
#' não foram dadas para preservar a originalidade da Pesquisa.
#' @examples
#'
#' m0 <- glm(nema ~ offset(log(off)) + cult,
#' data = nematoide,
#' family = poisson)
#'
#' # Diagnóstico.
#' par(mfrow = c(2, 2))
#' plot(m0); layout(1)
#'
#' # Estimativas dos parâmetros.
#' summary(m0)
#'
#' # Quadro de deviance.
#' anova(m0, test = "Chisq")
#'
#' library(bbmle)
#'
#' # Poisson Generalizada.
#' m1 <- pgnz(formula(m0), data = nematoide)
#'
#' # Diferença de deviance.
#' 2 * diff(c(logLik(m0), logLik(m1)))
#'
#' # Estimativas dos parâmetros.
#' summary(m1)
#'
NULL
......@@ -37,13 +37,17 @@
#'
#' @references Suekane, R., Degrande, P. E., de Lima Junior, I. S., de
#' Queiroz, M. V. B. M., & Rigoni, E. R. (2013). Danos da
#' Mosca-Branca Bemisia Tabaci e distribuição vertical das ninfas em
#' cultivares de soja em casa de vegetação. Arquivos do Instituto
#' Biológico, 80(2), 151-158.
#' Mosca-Branca \emph{Bemisia Tabaci} e distribuição vertical das
#' ninfas em cultivares de soja em casa de vegetação. Arquivos do
#' Instituto Biológico, 80(2), 151-158.
#' \href{http://200.129.209.183/arquivos/arquivos/78/MESTRADO-DOUTORADO-AGRONOMIA/Dissertação Renato Suekane.pdf}{Dissertação de Renato Suekane.}
#' @examples
#'
#' data(ninfas)
#' str(ninfas)
#'
#' xtabs(~data + cult, data = ninfas)
#'
#' library(lattice)
#'
#' xyplot(ntot ~ dias | cult,
......@@ -63,4 +67,3 @@
#' as.table = TRUE,
#' layout = c(NA, 2))
NULL
#' @name panel.beeswarm
#' @aliases panel.beeswarm
#' @export
#' @author Walmes Zeviani
#' @title Função para Gráfico de Dispersão BeesWarm com a Lattice
#' @author Walmes Zeviani baseado no pacote
#' \href{http://www.cbs.dtu.dk/~eklund/beeswarm/}{beeswarm}.
#' @title Diagrama de Dispersão com Aranjo dos Pontos como Colmeia
#'
#' @description Essa função permite a construção de gráficos de
#' dispersão unidimensionais similares a função
#' \code{\link{stripchart}}. Porém com pontos de mesmas coordenadas
#' exibidos lado a lado e não sobrepostos.
#' @description Used to make scatter plot of discrete variables with no
#' overlapping points. Observations with the same y value are spread.
#'
#' @param x,y,subscripts,... Argumentos passados para a
#' função \code{\link[lattice]{xyplot}}.
#' @param spread Um escalar numérico a distância entre os pontos. Esse
#' valor é obtido por tentativa erro e toda vez que mudar as
#' dimensões do gráfico, eles precisam ser novamente fornecidos, no
#' entanto são valores na escala do eixo \code{x} e por isso são
#' baseados nas distâncias entre os níveis do fator representado
#' neste eixo. Como sugestão, abra sempre a janela gráfica
#' (\code{x11()}) ou faça a exportação (\code{png()}, \code{pdf()},
#' etc) com dimensões conhecidas e calibre o \code{spred} para que
#' seja exibido adequadamente.
#'
#' @param r Valor para espaçamento entre os pontos de mesmas
#' coordenadas.
#' @param x,y,subscripts,... Argumentos passados para a
#' \code{\link[lattice]{panel.xyplot}}.
#'
#' @return Painel \code{xyplot} padrão, com as coordenadas \code{x}
#' devidamente corrigidas para agrupamento lado a lado.
#' @return A função passa conteúdo para o argumento \code{panel}.
#'
#' @seealso \code{\link[lattice]{xyplot}}.
#' @import lattice
#'
#' @references Aron Eklund (2016). beeswarm: The Bee Swarm Plot, an
#' Alternative to Stripchart. R package version
#' 0.2.3. \url{https://CRAN.R-project.org/package=beeswarm}
#' @examples
#'
#' data(capdesfo)
#' str(capdesfo)
#'
#' library(lattice)
#'
#' set.seed(2016)
#' da <- data.frame(x = rep(letters[1:5], 10), y = rpois(50, 5))
#' # x11(width = 7, height = 2.8)
#' xyplot(ncap ~ des | est, data = capdesfo,
#' layout = c(5, 1), as.table = TRUE,
#' type = c("p", "smooth"), col = 1, col.line = "gray50",
#' xlim = extendrange(c(0:1), f = 0.15),
#' xlab = "Nível de desfolha artificial",
#' ylab = "Número de capulho produzidos",
#' spread = 0.07, panel = panel.beeswarm)
#'
#' xyplot(y ~ x, data = da, jitter.x = TRUE)
#' xyplot(y ~ x, data = da, panel = panel.beeswarm, r = 0.1)
#' xyplot(y ~ x, data = da, panel = panel.beeswarm, r = 0.05)
#'
#' xyplot(drat ~ carb | am, data = mtcars, jitter.x = TRUE)
#' xyplot(drat ~ carb | am, data = mtcars, panel = panel.beeswarm, r = 0.2)
#' # x11(width = 7, height = 2.8)
#' xyplot(ncap ~ est | factor(des), data = capdesfo,
#' layout = c(5, 1), as.table = TRUE,
#' type = c("p", "smooth"), col = 1, col.line = "gray50",
#' xlab = "Fase de desenvolvimento da planta",
#' ylab = "Número de capulhos produzidos",
#' scales = list(x = list(
#' at = 1:nlevels(capdesfo$est),
#' labels = substr(levels(capdesfo$est),
#' start = 1, stop = 3))),
#' spread = 0.35, panel = panel.beeswarm)
#'
panel.beeswarm <- function(x, y, subscripts, r, ...) {
panel.beeswarm <- function(x, y, subscripts, spread, ...){
xx <- x
yy <- y
aux <- by(cbind(yy, xx, subscripts), xx, function(i) {
or <- order(i[, 1])
ys <- i[or, 1]
yt <- table(ys)
dv <- sapply(unlist(yt), function(j) {
seq(1, j, l = j) - (j + 1)/2
})
if (!is.list(dv)) {
dv <- as.list(dv)
}
xs <- i[or, 2] + r * do.call(c, dv)
cbind(x = xs, y = ys, subscripts[or])
})
aux <- by(cbind(yy, xx, subscripts),
INDICES = xx,
FUN = function(i) {
or <- order(i[, 1])
ys <- i[or, 1]
yt <- table(ys)
dv <- sapply(unlist(yt),
FUN = function(j) {
seq(from = 1,
to = j,
length.out = j) -
(j + 1)/2
})
if (!is.list(dv)) {
dv <- as.list(dv)
}
xs <- i[or, 2] + spread * do.call(c, dv)
cbind(x = xs, y = ys, subscripts = subscripts[or])
})
aux <- do.call(rbind, aux)
panel.xyplot(aux[, 1], aux[, 2], subscripts = aux[, 3], ...)
}
......@@ -7,8 +7,6 @@
#' @description Essa função permite representar intervalos de confiança
#' e bandas de confiança em gráficos do pacote lattice.
#'
#' @param y Valor central ou estimativa pontual.
#'
#' @param ly Limite inferior do intervalo/banda de confiança.
#'
#' @param uy Limite superior do intervalo/banda de confiança.
......@@ -21,7 +19,7 @@
#' quantidades a somar/subtrair dos valores de \code{x} para não
#' sobrepor intervalos. Com esse argumento pode-se representar mais
#' de um intervalor por valor de \code{x}. Não é usado quando
#' \code{cty = "bands"}
#' \code{cty = "bands"}.
#'
#' @param fill Uma representação de cor para preencher o interior das
#' bandas de confiança. Não é usado quando \code{cty = "bars"}.
......@@ -34,8 +32,9 @@
#' formam o "T". O valor default é 0.05. Não é usado quando
#' \code{cty = "bands"}.
#'
#' @param x,subscripts,col.line,lwd,... Argumentos passados para a
#' função \code{\link[lattice]{xyplot}}.
#' @param x,y,subscripts,col.line,lwd,... Argumentos passados para a
#' \code{\link[lattice]{panel.xyplot}} pela função
#' \code{\link[lattice]{xyplot}}.
#'
#' @return São usadas dentro de funções do pacote \pkg{lattice}.
#'
......@@ -127,7 +126,6 @@ panel.cbH <- function(x, y, ly, uy,
#' @name prepanel.cbH
#' @rdname panel.cbH
#' @author Walmes Zeviani baseado na lista de discussão R-help.
#' @export
prepanel.cbH <- function(y, ly, uy, subscripts){
ly <- as.numeric(ly[subscripts])
......
#' @name panel.groups.segplot
#' @export
#' @author Walmes Zeviani
#' @title Painel para Fazer Intervalos para Grupos no \code{segplot()}
#'
#' @description Essa função permite fazer intervalos de confiança (ou
#' com barras de erro) em gráficos da \code{latticeExtra::segplot()}
#' para grupos (\code{groups =}) de tal forma que eles não fiquem
#' sobrepostos.
#'
#' @param x,y,z,centers,data,subscripts,... Argumentos passados para a
#' \code{\link[latticeExtra]{segplot}}.
#'
#' @param groups A variável (\code{factor}) de agrupamento, de mesmo
#' comprimento de \code{lwr} e \code{upr}.
#'
#' @param gap Escalar que representa a distância entre os segmentos. O
#' valor default é 0,1. Como um fator com \eqn{k} níveis é
#' representado pelos números inteiros \eqn{1, 2, \ldots, k}, então
#' \eqn{0 \leq \textrm{gap} < 1/k}. Se for usado um valor negativo,
#' os intervalos serão apresentados em ordem inversa.
#'
#' @return A função é passada para o argumento \code{panel} e retorna
#' elementos necessários para a \code{\link[latticeExtra]{segplot}}.
#'
#' @seealso \code{\link[latticeExtra]{segplot}}.
#' @import latticeExtra
#' @examples
#'
#' library(latticeExtra)
#'
#' m0 <- lm(log(breaks) ~ wool * tension, data = warpbreaks)
#'
#' pred <- with(warpbreaks, expand.grid(KEEP.OUT.ATTRS = TRUE,
#' wool = levels(wool),
#' tension = levels(tension)))
#'
#' pred <- cbind(pred,
#' predict(m0, newdata = pred, interval = "confidence"))
#' str(pred)
#'
#' segplot(wool ~ lwr + upr, centers = fit, data = pred,
#' draw = FALSE, horizontal = FALSE)
#'
#' segplot(wool ~ lwr + upr, centers = fit, data = pred,
#' draw = FALSE, horizontal = FALSE,
#' groups = tension, gap = NULL,
#' panel = panel.groups.segplot)
#'
panel.groups.segplot <- function(x, y, z, centers,
groups, gap = NULL,
data, subscripts, ...) {
if (!missing(data)) {
data <- eval(data, envir = parent.frame())
groups <- data[, deparse(substitute(groups))]
}
stopifnot(is.factor(groups))
stopifnot(length(groups) == length(z))
if (is.null(gap)) {
gap <- 0.5/nlevels(groups)
}
d <- 2 * ((as.numeric(groups) - 1)/(nlevels(groups) - 1)) - 1
z <- as.numeric(z) + gap * d
panel.segplot(x, y, z, centers = centers,
subscripts = subscripts, ...)
}
#' @title Log-Verossimilhança do Modelo Poisson Generalizada
#' @author Walmes Zeviani, \email{walmes@@ufpr.br}
#' @export
#' @description Calcula a log-verossimilhança de um modelo de regressão
#' para a média com distribuição Poisson Generalizada para a
#' resposta de contagem (y).
#' @details A função de log-verossimilhança da Poisson Generalizada, na
#' parametrização de modelo de regressão é:
#'
#' \deqn{\ell(\lambda,\alpha,y) =
#' y (\log(\lambda) - \log(1 + \alpha\lambda)) +
#' (y - 1) \log(1 + \alpha y) -
#' \lambda\left(\frac{1 + \alpha y}{1 + \alpha\lambda}\right) -
#' \log(y!), }
#'
#' em que \eqn{\alpha} é o parâmetro de dispersão e \eqn{\lambda > 0} é
#' a média \eqn{E(Y) = \lambda} e \eqn{y = 0,1,ldots} é vetor
#' observado da variável de contagem. Nessa parametrização,
#' \eqn{V(Y) = \lambda (1 + \alpha\lambda)^2}. A Poisson
#' Generalizada em a Poisson como caso particular quando \eqn{\alpha
#' = 0}.
#'
#' Para o modelo de regressão, um preditor linear é ligado à média pela
#' função de ligação log, \eqn{\log(\lambda) = X\beta}, tal como é
#' para o modelo Poisson com link log.
#'
#' O espaço paramétrico de \eqn{\alpha} não limitado para o lado direito
#' do zero (\eqn{\alpha} positivo) mas para o lado esquerdo
#' (\eqn{\alpha} negativo) o limite inferior é dependente do
#' parâmetro \eqn{\lambda} e dos valores observados de
#' \eqn{y}. Valores não finitos podem ser retornados durante a
#' estimação quando \eqn{1 + \alpha\lambda} ou \eqn{1 + \alpha y}
#' não forem maiores que zero.
#' @param params Um vetor de (estimativas dos) parâmetros do modelo de
#' regressão. O primeiro elemento desse vetor é o parâmetro de
#' dispersão do modelo e os restantes são parâmetros de
#' locação. Essa função retorna o negativo da log-verossimilhança
#' pois foi construída para ser usada na \code{\link[bbmle]{mle2}}.
#' @param y Um vetor com variável dependente do modelo, resposta do tipo
#' contagem.
#' @param X A matriz de delineamento correspondente ao modelo linear
#' ligado à média pela função de ligação log. A matriz do modelo
#' pode ser construída com a função
#' \code{\link[stats]{model.matrix}}.
#' @param offset Um vetor, de mesmo comprimento de \code{y}, com valores
#' que correspondem aos offsets (ou exposição) para cada valor
#' observado. Se \code{NULL}, é usado 1 como offset.
#' @return O negativo da log-verossimilhança do modelo Poisson
#' Generalizado para os parâmetros e dados informados.
#' @seealso \code{\link[bbmle]{mle2}}.
#' @examples
#'
#' set.seed(123)
#' y <- rpois(10, lambda = 5)
#'
#' # Log-verossimilhança pela Poisson.
#' sum(dpois(y, lambda = 5, log = TRUE))
#'
#' # Log-verossimilhança pela PGNZ usando alpha = 0
#' llpgnz(params = c(0, log(5)), y = y, X = cbind(y * 0 + 1))
#'
#' set.seed(121)
#' y <- rpois(100, lambda = exp(1))
#' X <- cbind(0 * y + 1)
#'
#' grid <- expand.grid(alpha = seq(-0.1, 0.4, by = 0.01),
#' lambda = seq(0.1, 2.1, by = 0.025))
#'
#' grid$ll <- apply(grid, MARGIN = 1,
#' FUN = function(vec) {
#' llpgnz(params = vec, y = y, X = X, offset = NULL)
#' })
#'
#' library(latticeExtra)
#'
#' levelplot(ll ~ alpha + lambda, data = grid) +
#' layer(panel.abline(v = 0, h = 1, lty = 2))
#'
llpgnz <- function(params, y, X, offset = NULL) {
# params: vetor de parâmetros;
# params[1]: parâmetro de dispersão (alpha);
# params[-1]: parâmetro de locação (lambda);
# y: variável resposta (contagem);
# X: matriz do modelo linear;
# offset: tamanho do domínio onde y foi medido (exposição);
#----------------------------------------
if (is.null(offset)) {
offset <- 1L
}
alpha <- params[1]
lambda <- offset * exp(X %*% params[-1])
z <- 1 + alpha * lambda
w <- 1 + alpha * y
ll <- y * (log(lambda) - log(z)) +
(y - 1) * log(w) -
lambda * (w/z) -
lfactorial(y)
# Negativo da log-likelihood.
-sum(ll)
}
#' @title Probabilidades do Modelo Poisson Generalizado
#' @author Walmes Zeviani, \email{walmes@@ufpr.br}
#' @export
#' @description Calcula as probabilidades para uma variável aleatória
#' com distribuição Poisson Generalizada na parametrização
#' \eqn{\lambda-\alpha}:
#'
#' \deqn{p(y,\lambda,\alpha) =
#' \left(\frac{\lambda}{1+\alpha\lambda}\right)^{y}
#' \,\frac{(1+\alpha y)^{y-1}}{y!}
#' \exp\left\{-\lambda\left(
#' \frac{1+\alpha y}{1+\alpha\lambda}\right)\right\},
#' }
#'
#' em que \eqn{\lambda > 0} é a média da variável aleatória e
#' \eqn{\alpha} é o parâmetro de dispersão, sendo que \eqn{V(Y) =
#' \lambda (1 + \alpha\lambda)^2}. O espaço paramétrico de
#' \eqn{\alpha} depende de \eqn{\lambda} e \eqn{y} pois
#' \eqn{1+\alpha\lambda > 0} e \eqn{1+\alpha y > 0}.
#' @param y Valor da variável de contagem.
#' @param lambda Valor do parâmetro \eqn{\lambda} que é a média da
#' distribuição .
#' @param alpha Valor do parâmetro \eqn{\alpha} que é o parâmetro de
#' dispersão.
#' @return Retorna uma probabilidade, ou seja \eqn{\Pr(Y = y | \lambda,
#' \alpha) = p(y, \lambda, \alpha)}.
#' @examples
#'
#' dpois(5, lambda = 5)
#' dpgnz(5, lambda = 5, alpha = 0)
#'
#' probs <- data.frame(y = 0:30)
#' within(probs, {
#' py0 <- dpois(y, lambda = 15)
#' py1 <- dpgnz(y, lambda = 15, alpha = 0)
#' py2 <- dpgnz(y, lambda = 15, alpha = 1/30)
#' py3 <- dpgnz(y, lambda = 15, alpha = -1/30)
#' plot(py0 ~ y, type = "h",
#' ylim = c(0, max(c(py0, py2, py3))),
#' ylab = expression(Pr(Y == y)))
#' points(y + 0.1, py1, type = "h", col = 2)
#' points(y - 0.3, py2, type = "h", col = 3)
#' points(y + 0.3, py3, type = "h", col = 4)
#' legend("topleft", bty = "n",
#' col = c(1:4), lty = 1,
#' legend = expression(
#' Poisson(lambda == 15),
#' PG(lambda == 15, alpha == 0),
#' PG(lambda == 15, alpha == 1/30),
#' PG(lambda == 15, alpha == -1/30)))
#' })
#'
dpgnz <- function(y, lambda, alpha) {
k <- lfactorial(y)
w <- 1 + alpha * y
z <- 1 + alpha * lambda
m <- alpha > pmax(-1/y, -1/lambda)
fy <- y * (log(lambda) - log(z)) +
(y - 1) * log(w) - lambda * (w/z) - k
fy[!m] <- 0
return(m * exp(fy))
}
#' @title Ajuste do Modelo Poisson Generalizado
#' @author Walmes Zeviani, \email{walmes@@ufpr.br}
#' @export
#' @description Estima os parâmetros de um modelo Poisson Generalizado
#' pela otimização da função de log-verossimilhança definida em
#' \code{\link{llpgnz}}. A sintaxe assemelha-se com a função
#' \code{\link[stats]{glm}} (Generalized Linear Models).
#' @param formula Um objeto da classe \code{\link{formula}}. Se
#' necessária a inclusão de \emph{offset} deve-se indicá-lo como
#' \code{\link{offset}}.
#' @param data Um objeto de classe \code{data.frame} que contém as
#' variáveis descritas na \code{formula}.
#' @param start Um vetor com os valores iniciais para os parâmetros do