diff --git a/NAMESPACE b/NAMESPACE index 26d8bebfc1e5522d510b6753cfc68bab989fecfb..eb464be636a32c83c64561f92a6c78b642062ec5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,8 +1,11 @@ # Generated by roxygen2: do not edit by hand export(apc) +export(cmp) +export(convergencez) export(panel.cbH) export(prepanel.cbH) +import(bbmle) import(doBy) import(lattice) import(latticeExtra) diff --git a/R/cmp.R b/R/cmp.R new file mode 100644 index 0000000000000000000000000000000000000000..e0b4e9d5cc9bee0a6f783001222826cca515adeb --- /dev/null +++ b/R/cmp.R @@ -0,0 +1,157 @@ +#' @title Avaliação da Convergência da Constante Normalizadora +#' @description Avalia a convergência da constante de normalização de +#' um modelo COM-Poisson definida por: \deqn{Z = \sum +#' \frac{\lambda^i}{(i!)^\nu}}, em que o parâmetro \eqn{\nu} é +#' tomado como \eqn{\exp{\phi}}. +#' @param model Objeto resultante da função \code{\link[MRDCr]{cmp}}. +#' @param tol Critério de parada do algoritmo, representa o valor +#' tolerado para a diferença do valor de \eqn{Z(\lambda, \phi)} +#' entre duas iterações. O valor padrão é 1e-4 +#' @param incremento Número de incrementos da soma a serem considerados +#' a cada iteração. Padrão definido como 10, ou seja, a cada +#' iteração 10 passos incrementos são somados a Z. +#' @param maxit Número máximo de iterações a serem realizadas pelo +#' algoritmo. Se este número for atingido e o critério de tolerância +#' não for atendido, uma mensagem de aviso será exibida. +#' @param plot Argumento lógico. Se \code{TRUE} (padrão) os gráficos dos +#' incrementos da constante são exibidos. +#' @return Uma lista com os incrementos das constantes Z, +#' \eqn{Z(\lambda, \phi)} da distribuição COM-Poisson, calculados +#' para cada observação. +#' @author Eduardo E. R. Junior, \email{edujrrib@gmail.com} +#' @export + +convergencez <- function(model, tol = 1e-4, incremento = 10, + maxit = 50, plot = TRUE) { + ##------------------------------------------- + ## Calcula Z para um c(lambda, phi) + calcula <- function(loglambda, phi) { + nu <- exp(phi) + zg <- vector("list", maxit) + t <- incremento + i <- 1:t + j <- 1 + zg[[j]] <- exp(i * loglambda - nu * lfactorial(i)) + dif <- abs(zg[[j]][t-1] - zg[[j]][t]) + while (dif > tol && is.finite(dif) && j <= maxit) { + i <- (i[t] + 1):(i[t] + t) + j <- j + 1 + zg[[j]] <- exp(i * loglambda - nu * lfactorial(i)) + dif <- abs(zg[[j]][t-1] - zg[[j]][t]) + } + if (j == maxit && dif > tol || !is.finite(dif)) { + print(c(loglambda, nu)) + warning("Soma n\\u00e3o convergiu") + } + return(z = unlist(zg)) + } + ##------------------------------------------- + X <- model@data$X + y <- model@data$y + sumto <- model@data$sumto + betas <- model@coef[-1] + phi <- model@coef[1] + loglambdas <- X %*% betas + zs <- sapply(loglambdas, calcula, phi, simplify = FALSE) + ##------------------------------------------- + if (plot) { + mx <- max(sapply(zs, max)) + lx <- max(sapply(zs, length)) + plot(zs[[1]], type = "l", xlim = c(0, lx), ylim = c(0, mx)) + abline(v = sumto) + for (i in 2:length(zs)) lines(zs[[i]], type = "l") + } + invisible(zs) +} + +#' @title Log-Verossimilhança do Modelo Conway-Maxwell-Poisson +#' @description Calcula a log-verossimilhança de um modelo COM-Poisson +#' considerando os dados e as estimativas dos parâmetros informadas. +#' @details A função de log-verossimilhança toma a forma: \deqn{-Z - y * +#' \lambda - \nu \log{y!}}, onde \eqn{Z = \sum +#' \frac{\lambda^i}{(i!)^\nu}} e \eqn{\nu = \exp{\phi}}. +#' @param params Um vetor de parâmetros do modelo COM-Poisson. É +#' necessário que seja informado como primeiro elemento do vetor, o +#' valor de \eqn{\phi}. Os demais elementos são os \eqn{\beta}'s do +#' preditor linear \eqn{\lambda}. +#' @param y Um vetor de contagens, considerado como variável resposta. +#' @param X A matriz de delineamento do modelo. +#' @param offset Um vetor de valores a serem adicionados ao preditor +#' linear. +#' @param sumto Número de incrementos a serem considerados para a soma +#' da constante normalizadora. Como padrão, o número de incrementos +#' é o valor inteiro de \eqn{(\max y)^1.2}. +#' @return O valor da log-verossimilhança do modelo +#' Conway-Maxwell-Poisson com os parâmetros e dados informados. +#' @author Eduardo E. R. Junior, \email{edujrrib@gmail.com} +#' @seealso \code{\link[bbmle]{mle2}} + +llcmp <- function(params, y, X, offset = NULL, + sumto = ceiling(max(y)^1.2)){ + ##------------------------------------------- + betas <- params[-1] + phi <- params[1] + nu <- exp(phi) + ##------------------------------------------- + if (is.null(offset)) offset <- 0 + Xb <- X %*% betas + offset + ##------------------------------------------- + ## Obtendo a constante normatizadora Z. + i <- 1:sumto + zs <- sapply(Xb, function(loglambda) + sum(exp(i * loglambda - nu * lfactorial(i)))) + Z <- sum(log(zs + 1)) + ##------------------------------------------- + ll <- sum(y * Xb - nu * lfactorial(y)) - Z + return(-ll) +} + +#' @title Ajuste do Modelo Conway-Maxwell-Poisson +#' @description Estima os parâmetros de um modelo COM-Poisson sob a +#' otimização da função de log-verossimilhança. A sintaxe +#' assemelha-se com a função \code{\link{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}, cujo contém as +#' variáveis descritas na \code{formula}. +#' @param start Um vetor de parâmetros do modelo tomados como valores +#' iniciais para o algoritmo iterativo. Se \code{NULL} os parâmetros +#' de um modelo log-linear Poisson e \eqn{\phi = 0} são tomados como +#' valores iniciais para o preditor e para o parâmetro de dispersão. +#' @param sumto Número de incrementos a serem considerados para a soma +#' da constante normalizadora. Como padrão, \code{NULL} o número de +#' incrementos é o valor inteiro de \eqn{(\max y)^1.2}. +#' @param ... Argumentos opcionais do framework de maximização numérica +#' \code{\link[bbmle]{mle2}}. +#' @return Um objeto de classe \code{mle2}, retornado da função de +#' \code{\link[bbmle]{mle2}}, usada para ajuste de modelos por +#' máxima verossimilhança. +#' @author Eduardo E. R. Junior, \email{edujrrib@gmail.com} +#' @import bbmle +#' @export + +cmp <- function(formula, data, start = NULL, sumto = NULL, ...) { + ##------------------------------------------- + ## Constrói as matrizes do modelo + frame <- model.frame(formula, data) + terms <- attr(frame, "terms") + y <- model.response(frame) + X <- model.matrix(terms, frame) + off <- model.offset(frame) + if (is.null(sumto)) sumto <- ceiling(max(y)^1.2) + ##------------------------------------------- + ## Define os chutes iniciais + if (is.null(start)) { + m0 <- glm.fit(x = X, y = y, offset = off, family = poisson()) + start <- c("phi" = 0, m0$coefficients) + } + ##------------------------------------------- + ## Nomeia os parâmetros da função para métodos bbmle + parnames(llcmp) <- names(start) + model <- bbmle::mle2(llcmp, start = start, + data = list(y = y, X = X, sumto = sumto), + vecpar = TRUE) + return(model) +} diff --git a/man/cmp.Rd b/man/cmp.Rd new file mode 100644 index 0000000000000000000000000000000000000000..d1832a220ea81ac1b580b3f0ea35b65ace6d80c8 --- /dev/null +++ b/man/cmp.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cmp.R +\name{cmp} +\alias{cmp} +\title{Ajuste do Modelo Conway-Maxwell-Poisson} +\usage{ +cmp(formula, data, start = NULL, sumto = NULL, ...) +} +\arguments{ +\item{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}}.} + +\item{data}{Um objeto de classe \code{data.frame}, cujo contém as +variáveis descritas na \code{formula}.} + +\item{start}{Um vetor de parâmetros do modelo tomados como valores +iniciais para o algoritmo iterativo. Se \code{NULL} os parâmetros +de um modelo log-linear Poisson e \eqn{\phi = 0} são tomados como +valores iniciais para o preditor e para o parâmetro de dispersão.} + +\item{sumto}{Número de incrementos a serem considerados para a soma +da constante normalizadora. Como padrão, \code{NULL} o número de +incrementos é o valor inteiro de \eqn{(\max y)^1.2}.} + +\item{...}{Argumentos opcionais do framework de maximização numérica +\code{\link[bbmle]{mle2}}.} +} +\value{ +Um objeto de classe \code{mle2}, retornado da função de + \code{\link[bbmle]{mle2}}, usada para ajuste de modelos por + máxima verossimilhança. +} +\description{ +Estima os parâmetros de um modelo COM-Poisson sob a + otimização da função de log-verossimilhança. A sintaxe + assemelha-se com a função \code{\link{glm}} (Generalized Linear + Models). +} +\author{ +Eduardo E. R. Junior, \email{edujrrib@gmail.com} +} + diff --git a/man/convergencez.Rd b/man/convergencez.Rd new file mode 100644 index 0000000000000000000000000000000000000000..ac4101be193ac1b6c87f5c43ff67424ff46333ed --- /dev/null +++ b/man/convergencez.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cmp.R +\name{convergencez} +\alias{convergencez} +\title{Avaliação da Convergência da Constante Normalizadora} +\usage{ +convergencez(model, tol = 1e-04, incremento = 10, maxit = 50, + plot = TRUE) +} +\arguments{ +\item{model}{Objeto resultante da função \code{\link[MRDCr]{cmp}}.} + +\item{tol}{Critério de parada do algoritmo, representa o valor +tolerado para a diferença do valor de \eqn{Z(\lambda, \phi)} +entre duas iterações. O valor padrão é 1e-4} + +\item{incremento}{Número de incrementos da soma a serem considerados +a cada iteração. Padrão definido como 10, ou seja, a cada +iteração 10 passos incrementos são somados a Z.} + +\item{maxit}{Número máximo de iterações a serem realizadas pelo +algoritmo. Se este número for atingido e o critério de tolerância +não for atendido, uma mensagem de aviso será exibida.} + +\item{plot}{Argumento lógico. Se \code{TRUE} (padrão) os gráficos dos +incrementos da constante são exibidos.} +} +\value{ +Uma lista com os incrementos das constantes Z, + \eqn{Z(\lambda, \phi)} da distribuição COM-Poisson, calculados + para cada observação. +} +\description{ +Avalia a convergência da constante de normalização de + um modelo COM-Poisson definida por: \deqn{Z = \sum + \frac{\lambda^i}{(i!)^\nu}}, em que o parâmetro \eqn{\nu} é + tomado como \eqn{\exp{\phi}}. +} +\author{ +Eduardo E. R. Junior, \email{edujrrib@gmail.com} +} + diff --git a/man/llcmp.Rd b/man/llcmp.Rd new file mode 100644 index 0000000000000000000000000000000000000000..afb7409c2869d31e67afd925a745098fb5bd1d3f --- /dev/null +++ b/man/llcmp.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cmp.R +\name{llcmp} +\alias{llcmp} +\title{Log-Verossimilhança do Modelo Conway-Maxwell-Poisson} +\usage{ +llcmp(params, y, X, offset = NULL, sumto = ceiling(max(y)^1.2)) +} +\arguments{ +\item{params}{Um vetor de parâmetros do modelo COM-Poisson. É +necessário que seja informado como primeiro elemento do vetor, o +valor de \eqn{\phi}. Os demais elementos são os \eqn{\beta}'s do +preditor linear \eqn{\lambda}.} + +\item{y}{Um vetor de contagens, considerado como variável resposta.} + +\item{X}{A matriz de delineamento do modelo.} + +\item{offset}{Um vetor de valores a serem adicionados ao preditor +linear.} + +\item{sumto}{Número de incrementos a serem considerados para a soma +da constante normalizadora. Como padrão, o número de incrementos +é o valor inteiro de \eqn{(\max y)^1.2}.} +} +\value{ +O valor da log-verossimilhança do modelo + Conway-Maxwell-Poisson com os parâmetros e dados informados. +} +\description{ +Calcula a log-verossimilhança de um modelo COM-Poisson + considerando os dados e as estimativas dos parâmetros informadas. +} +\details{ +A função de log-verossimilhança toma a forma: \deqn{-Z - y * + \lambda - \nu \log{y!}}, onde \eqn{Z = \sum + \frac{\lambda^i}{(i!)^\nu}} e \eqn{\nu = \exp{\phi}}. +} +\author{ +Eduardo E. R. Junior, \email{edujrrib@gmail.com} +} +\seealso{ +\code{\link[bbmle]{mle2}} +} +