cmp.R 6.96 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#' @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}}
88
#' @export
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

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)
}