From a57f13f0cc907d8f6ad729aa0cde2381fa1a68b6 Mon Sep 17 00:00:00 2001 From: Walmes Zeviani Date: Thu, 5 May 2016 21:27:45 -0300 Subject: [PATCH] =?UTF-8?q?Retira=20as=20defini=C3=A7=C3=B5es=20de=20fun?= =?UTF-8?q?=C3=A7=C3=B5es.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- vignettes/v04_poisson_generelizada.Rmd | 337 +++---------------------- 1 file changed, 34 insertions(+), 303 deletions(-) diff --git a/vignettes/v04_poisson_generelizada.Rmd b/vignettes/v04_poisson_generelizada.Rmd index 35f3df9..2570a2e 100644 --- a/vignettes/v04_poisson_generelizada.Rmd +++ b/vignettes/v04_poisson_generelizada.Rmd @@ -56,7 +56,7 @@ $$ ```{r} # Função densidade na parametrização original. -dpg0 <- function(y, theta, gamma, m = 4) { +dpgnz0 <- function(y, theta, gamma, m = 4) { if (gamma < 0) { m <- max(c(m, floor(-theta/gamma))) if (gamma < max(c(-1, -theta/m))) { @@ -79,10 +79,9 @@ y <- 0:30 theta <- 10 gamma <- 0 -fy <- dpg0(y = y, theta = theta, gamma = gamma) -plot(fy ~ y, type = "h") -lines(y + 0.3, dpois(y, lambda = theta), type = "h", col = 2, - xlab = "y", ylab = "f(y)") +fy <- dpgnz0(y = y, theta = theta, gamma = gamma) +plot(fy ~ y, type = "h", xlab = "y", ylab = "f(y)") +lines(y + 0.3, dpois(y, lambda = theta), type = "h", col = 2) ``` ## Recursos interativos com o `rpanel` ## @@ -96,7 +95,7 @@ react <- function(panel){ from <- floor(max(c(0, m - 5 * s))) to <- ceiling(max(c(YMAX, m + 5 * s))) y <- from:to - py <- dpg0(y = y, theta = THETA, gamma = GAMMA) + py <- dpgnz0(y = y, theta = THETA, gamma = GAMMA) if (POIS) { pz <- dpois(y, lambda = m) } else { @@ -147,17 +146,7 @@ rp.do(panel = panel, action = react) ```{r, eval=FALSE} # Função densidade na parametrização de modelo de regressão. -dpg1 <- function(y, lambda, alpha) { - k <- lfactorial(y) - w <- 1 + alpha * y - z <- 1 + alpha * lambda - m <- alpha > pmax(-1/y, -1/lambda) - # fy <- (lambda/z)^(y) * w^(y - 1) * exp(-lambda * (w/z))/exp(k) - fy <- y * (log(lambda) - log(z)) + - (y - 1) * log(w) - lambda * (w/z) - k - fy[!m] <- 0 - return(m * exp(fy)) -} +MRDCr::dpgnz react <- function(panel){ with(panel, @@ -167,7 +156,7 @@ react <- function(panel){ from <- floor(max(c(0, m - 5 * s))) to <- ceiling(max(c(YMAX, m + 5 * s))) y <- from:to - py <- dpg1(y = y, lambda = LAMBDA, alpha = ALPHA) + py <- dpgnz(y = y, lambda = LAMBDA, alpha = ALPHA) if (POIS) { pz <- dpois(y, lambda = m) } else { @@ -218,14 +207,9 @@ rp.do(panel = panel, action = react) #----------------------------------------------------------------------- # Gráfico do espaço paramétrico de theta x gamma. -# debug(dpg1) -# dpg1(y = 0:10, lambda = 1, alpha = 0) -# dpg1(y = 0:10, lambda = 1, alpha = -0.1) -# undebug(dpg1) - fun <- Vectorize(vectorize.args = c("theta", "gamma"), FUN = function(theta, gamma) { - sum(dpg0(y = y, theta = theta, gamma = gamma)) + sum(dpgnz0(y = y, theta = theta, gamma = gamma)) }) grid <- list(theta = seq(1, 50, by = 1), @@ -250,7 +234,7 @@ levelplot(sum ~ theta + gamma, fun <- Vectorize(vectorize.args = c("lambda", "alpha"), FUN = function(lambda, alpha) { - sum(dpg1(y = y, lambda = lambda, alpha = alpha)) + sum(dpgnz(y = y, lambda = lambda, alpha = alpha)) }) grid <- list(lambda = seq(0.2, 50, by = 0.2), @@ -275,28 +259,7 @@ levelplot(sum ~ lambda + alpha, ```{r} # Função de log-Verossimilhança da Poisson Generalizada na # parametrização de modelo de regressão. -llpg <- function(theta, y, X, offset = NULL) { - # theta: vetor de parâmetros; - # theta[1]: parâmetro de dispersão (alpha); - # theta[-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; - #---------------------------------------- - if (is.null(offset)) { - offset <- 1L - } - lambda <- offset * exp(X %*% theta[-1]) - alpha <- theta[1] - w <- 1 + alpha * y - z <- 1 + alpha * lambda - fy <- y * (log(lambda) - log(z)) + - (y - 1) * log(w) - - lambda * (w/z) - - lfactorial(y) - # Negativo da log-likelihood. - -sum(fy) -} +MRDCr::llpgnz #----------------------------------------------------------------------- # Gerando uma amostra aleatória da Poisson. @@ -309,10 +272,10 @@ L <- list(y = y, X = cbind(rep(1, length(y)))) start <- c(alpha = 0, lambda = 1) -parnames(llpg) <- names(start) +parnames(llpgnz) <- names(start) # Como \alpha foi fixado em 1, essa ll corresponde à Poisson. -n0 <- mle2(minuslogl = llpg, +n0 <- mle2(minuslogl = llpgnz, start = start, data = L, fixed = list(alpha = 0), vecpar = TRUE) @@ -321,7 +284,7 @@ c(coef(n0)["lambda"], coef(glm(y ~ offset(log(L$offset)), family = poisson))) # Estimando o \alpha. -n1 <- mle2(llpg, start = start, data = L, vecpar = TRUE) +n1 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE) coef(n1) # Perfil de verossimilhança dos parâmetros. @@ -342,8 +305,6 @@ nesse experimento foram o número de vagens viáveis (e não viáveis) e o número total de sementes por parcela. ```{r} -library(lattice) - data(soja, package = "MRDCr") str(soja) @@ -376,10 +337,10 @@ L <- with(soja, list(y = nvag, offset = 1, X = model.matrix(m0))) # Usa as estimativas do Poisson como valore iniciais. start <- c(alpha = 0, coef(m0)) -parnames(llpg) <- names(start) +parnames(llpgnz) <- names(start) # Com alpha fixo em 0 corresponde à Poisson. -m2 <- mle2(llpg, start = start, data = L, +m2 <- mle2(llpgnz, start = start, data = L, fixed = list(alpha = 0), vecpar = TRUE) # Mesma medida de ajuste e estimativas. @@ -387,7 +348,7 @@ c(logLik(m2), logLik(m0)) cbind(coef(m2)[-1], coef(m0)) # Poisson Generalizada. -m3 <- mle2(llpg, start = start, data = L, vecpar = TRUE) +m3 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE) # Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0). anova(m3, m2) @@ -517,17 +478,17 @@ L <- with(soja, list(y = ngra, offset = 1, X = model.matrix(m0))) # Usa as estimativas do Poisson como valore iniciais. start <- c(alpha = 0, coef(m0)) -parnames(llpg) <- names(start) +parnames(llpgnz) <- names(start) # Com alpha fixo em 0 corresponde à Poisson. -m2 <- mle2(llpg, start = start, data = L, +m2 <- mle2(llpgnz, 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(llpg, start = start, data = L, vecpar = TRUE) +m3 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE) # Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0). anova(m3, m2) @@ -657,17 +618,17 @@ L <- with(soja, list(y = ngra, offset = nvag, X = model.matrix(m0))) # perfilhar encontra um mínimo melhor. Então, por tentativa erro # chegou-se no -0.0026. start <- c(alpha = -0.0026, coef(m0)) -parnames(llpg) <- names(start) +parnames(llpgnz) <- names(start) # Com alpha fixo em 0 corresponde à Poisson. -m2 <- mle2(llpg, start = start, data = L, +m2 <- mle2(llpgnz, 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(llpg, start = start, data = L, vecpar = TRUE) +m3 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE) # Teste para nulinidade do parâmetro de dispersão (H_0: alpha == 0). anova(m3, m2) @@ -741,16 +702,16 @@ summary(m0) L <- with(capdesfo, list(y = ncap, offset = 1, X = model.matrix(m0))) start <- c(alpha = log(1), coef(m0)) -parnames(llpg) <- names(start) +parnames(llpgnz) <- names(start) # Modelo Poisson também. -m2 <- mle2(llpg, start = start, data = L, +m2 <- mle2(llpgnz, start = start, data = L, fixed = list(alpha = 0), vecpar = TRUE) c(logLik(m2), logLik(m0)) # Modelo Poisson Generalizado. -m3 <- mle2(llpg, start = start, data = L, vecpar = TRUE) +m3 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE) logLik(m3) anova(m3, m2) @@ -848,249 +809,19 @@ update(p1, type = "p", layout = c(NA, 1), as.layer(p2, under = TRUE) ``` -## Mosca Branca ## - -```{r} -data(ninfas, package = "MRDCr") -str(ninfas) - -# Somente as cultivares que contém BRS na identificação -ninfas <- droplevels(subset(ninfas, - grepl("BRS.*RR", x = cult) & dias <= 22)) -ninfas$y <- ninfas$ntot -str(ninfas) - -xyplot(y ~ dias | cult, data = ninfas, - type = c("p", "spline"), - grid = TRUE, as.table = TRUE, layout = c(NA, 2), - xlab = "Dias", ylab = "Número total de ninfas") - -ninfas <- transform(ninfas, aval = factor(dias)) - -#----------------------------------------------------------------------- -# Modelo Poisson. - -m0 <- glm(y ~ bloco + cult + aval, - data = ninfas, family = poisson) - -par(mfrow = c(2, 2)) -plot(m0); layout(1) - -anova(m0, test = "Chisq") -summary(m0) - -#----------------------------------------------------------------------- -# Modelo Poisson Generalizada. - -L <- with(ninfas, list(y = y, offset = 1, X = model.matrix(m0))) - -start <- c(alpha = 0, coef(m0)) -parnames(llpg) <- names(start) - -# Modelo Poisson também. -m2 <- mle2(llpg, start = start, data = L, - fixed = list(alpha = 0), vecpar = TRUE) - -c(logLik(m2), logLik(m0)) - -# Modelo Poisson Generalizado. -m3 <- mle2(llpg, start = start, data = L, vecpar = TRUE) -logLik(m3) - -anova(m3, m2) - -summary(m3) - -plot(profile(m3, which = "alpha")) - -cbind("PoissonGLM" = c(NA, coef(m0)), - "PoissonML" = coef(m2), - "PGeneraliz" = coef(m3)) - -V <- cov2cor(vcov(m3)) -corrplot.mixed(V, upper = "ellipse", col = "gray50") -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(ninfas, expand.grid(bloco = factor(levels(bloco)[1], - levels = levels(bloco)), - cult = levels(cult), - aval = levels(aval), - KEEP.OUT.ATTRS = FALSE)) -X <- model.matrix(formula(m0)[-2], data = pred) - -bl <- attr(X, "assign") == 1 -X[, bl] <- X[, bl] * 0 + 1/(sum(bl) + 1) -head(X) - -pred <- list(pois = pred, pgen = 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)) -str(pred$pois) - -# Matrix de covariância completa e sem o alpha. -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$pgen$eta <- c(X %*% coef(m3)[-1]) -pred$pgen <- cbind(pred$pgen, - apply(outer(aux, qn, FUN = "*"), MARGIN = 2, - FUN = function(x) { - exp(pred$pgen$eta + x) - })) -str(pred$pgen) - -pred <- ldply(pred, .id = "modelo") -pred <- arrange(pred, cult, aval, modelo) -str(pred) - -key <- list(lines = list( - lty = 1, - col = trellis.par.get("superpose.line")$col[1:2]), - text = list( - c("Poisson", "Poisson Generalizada"))) - -xyplot(y ~ aval | cult, data = ninfas, - grid = TRUE, as.table = TRUE, layout = c(NA, 2), - xlab = "Dias", ylab = "Número total de ninfas", - key = key) + - as.layer(xyplot(fit ~ aval | cult, data = pred, - groups = modelo, pch = 19, type = "o", - ly = pred$lwr, uy = pred$upr, - cty = "bars", length = 0.05, - desloc = 0.2 * scale(as.integer(pred$modelo), - scale = FALSE), - prepanel = prepanel.cbH, - panel.groups = panel.cbH, - panel = panel.superpose), under = TRUE) - -#----------------------------------------------------------------------- - -fun <- Vectorize(vectorize.args = c("lambda", "alpha"), - FUN = function(lambda, alpha) { - sum(dpg1(y = y, lambda = lambda, alpha = alpha)) - }) - -# dpg1(y = 0:10, lambda = 5, alpha = -0) -# dpois(0:10, lambda = 5) - -head(sort(subset(pred, modelo = "pois")$fit, decreasing = TRUE)) -coef(m3)["alpha"] - -y <- 0:400 -grid <- list(lambda = seq(10, 200, by = 2), - alpha = seq(-0.05, 0.1, by = 0.001)) -grid$sum <- with(grid, outer(lambda, alpha, fun)) - -funcur - -y <- 0:800 -py <- dpg1(y = y, lambda = 190, alpha = 0.05) -plot(py ~ y, type = "h") -abline(v = 190) - -grid <- with(grid, - cbind(expand.grid(lambda = lambda, alpha = alpha), - data.frame(sum = c(sum)))) - -levelplot(sum ~ lambda + alpha, - data = subset(grid, round(sum, 3) == 1), - col.regions = gray.colors) + - layer(panel.abline(h = 0.05)) - -``` - -```{r} -funcur <- function(lambda, alpha, n = 10) { - m <- lambda - s <- sqrt(lambda) * (1 + alpha * lambda) - sy <- seq(max(c(0, m - 4 * s)), - ceiling(m + 4 * s), - length.out = n) - sy <- round(sy, 0) - fy <- dpg1(y = sy, lambda = lambda, alpha = alpha) - fy <- 0.8 * fy/max(fy) - return(cbind(lambda = lambda, alpha = alpha, y = sy, fy = fy)) -} - -den <- subset(pred, modelo == "pgen") - -L <- lapply(den$fit, FUN = funcur, alpha = 0.05, n = 50) -for (i in seq_along(L)) { - L[[i]] <- cbind(den[i, ], L[[i]]) - # L[[i]]$fy <- 0.75 * L[[i]]$fy/max(L[[i]]$fy) -} - -L <- do.call(rbind, L) - -xyplot(y ~ aval | cult, data = L, - py = L$fy, type = "l", groups = aval, col = 1, - panel = function(x, y, subscripts, py, ...) { - panel.xyplot(x = as.integer(x) + py[subscripts], - y = y, - subscripts = subscripts, ...) - }) + - as.layer(xyplot(ntot ~ aval | cult, data = ninfas)) + - as.layer(xyplot(fit ~ aval | cult, data = pred, - groups = modelo, pch = 19, type = "o", - ly = pred$lwr, uy = pred$upr, - cty = "bars", length = 0.05, - desloc = 0.2 * scale(as.integer(pred$modelo), - scale = FALSE), - prepanel = prepanel.cbH, - panel.groups = panel.cbH, - panel = panel.superpose), under = TRUE) - - -# Teria que estimar uma variância para cada avaliação. - -``` - ## Pacote VGAM ## ```{r, eval=FALSE} #----------------------------------------------------------------------- -# # http://finzi.psych.upenn.edu/library/VGAM/html/genpoisson.html -# library(VGAM) -# -# formula(m0) -# m1 <- vglm(formula(m0), data = cap, family = genpoisson, trace = TRUE) -# coef(m1, matrix = TRUE) -# summary(m1) -# -# logLik(m2) +# http://finzi.psych.upenn.edu/library/VGAM/html/genpoisson.html +library(VGAM) + +m1 <- vglm(ncap ~ est * (des + I(des^2)), + data = capdesfo, family = genpoisson, trace = TRUE) +coef(m1, matrix = TRUE) +summary(m1) + +logLik(m1) ``` -- GitLab