Commit 71508b1a by Walmes Marques Zeviani

Remove chunck do pacote VGAM.

parent c99f458c
Pipeline #4374 failed with stage
This source diff could not be displayed because it is too large. You can view the blob instead.
 ## ----setup, include=FALSE----------------------------------------- source("_setup.R") ## ---- message=FALSE, error=FALSE, warning=FALSE------------------- # Definições da sessão. # devtools::load_all("../") library(lattice) library(latticeExtra) library(rpanel) library(bbmle) library(corrplot) library(plyr) library(car) library(doBy) library(multcomp) library(MRDCr) ## ----------------------------------------------------------------- # 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))) { ... ... @@ -15,8 +29,9 @@ dpg0 <- function(y, theta, gamma, m = 4) { m <- 1 } z <- theta + gamma * y k <- exp(lfactorial(y)) fy <- m * theta * z^(y - 1) * exp(-z)/k k <- lfactorial(y) # fy <- m * theta * z^(y - 1) * exp(-z)/exp(k) fy <- m * exp(log(theta) + (y - 1) * log(z) - z - k) return(fy) } ... ... @@ -25,13 +40,11 @@ y <- 0:30 theta <- 10 gamma <- 0 fy <- dpg0(y = y, theta = theta, gamma = gamma) plot(fy ~ y, type = "h") 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) ## ---- eval=FALSE-------------------------------------------------- # library(rpanel) # # react <- function(panel){ # with(panel, # { ... ... @@ -40,41 +53,37 @@ lines(y + 0.3, dpois(y, lambda = theta), type = "h", col = 2) # 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 { # pz <- 0 # } # if (any(!is.finite(py))) { # plot(x = NULL, y = NULL, # xlim = c(0, 1), ylim = c(0, 1), # axes = FALSE, ann = FALSE) # text(x = 0.5, y = 0.5, labels = "Valores não finitos") # } else { # 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) { # abline(v = m, col = 2) # } # if (POIS) { # lines(y + 0.3, pz, type = "h", col = 3) # } # # Colocar 0 para valores não finitos (-Inf, Inf e NaN) para # # fazer gráfico. # 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) { # abline(v = m, col = 2) # } # if (POIS) { # lines(y + 0.3, pz, type = "h", col = 3) # } # }) # panel # } # # panel <- rp.control(title = "Poisson Generalizada", # size = c(300, 100), YMAX = 30, m = 4) # size = c(300, 100), YMAX = 150, m = 4) # rp.slider(panel = panel, action = react, # variable = THETA, title = "theta", # from = 0.1, to = 30, # from = 0.1, to = 100, # initval = 5, resolution = 0.1, # showvalue = TRUE) # rp.slider(panel = panel, action = react, ... ... @@ -90,45 +99,9 @@ lines(y + 0.3, dpois(y, lambda = theta), type = "h", col = 2) # labels = "Adicionar a distribução Poisson?") # rp.do(panel = panel, action = react) ## ----------------------------------------------------------------- #----------------------------------------------------------------------- # Gráfico do espaço paramétrico de theta x gamma. library(latticeExtra) y <- 0:50 fun <- Vectorize(vectorize.args = c("theta", "gamma"), FUN = function(theta, gamma) { sum(dpg0(y = y, theta = theta, gamma = gamma)) }) grid <- list(theta = seq(0.1, 10, by = 0.1), gamma = seq(-0.98, 0.98, by = 0.02)) grid$sum <- with(grid, outer(theta, gamma, fun)) grid <- with(grid, cbind(expand.grid(theta = theta, gamma = gamma), data.frame(sum = c(sum)))) # levelplot(sum ~ theta + gamma, data = grid, # col.regions = heat.colors) + # layer(panel.abline(h = 0)) levelplot(sum ~ theta + gamma, data = subset(grid, round(sum, 3) == 1), col.regions = heat.colors) + layer(panel.abline(h = 0)) subset(grid, round(sum, 3) != 1 & theta > 6 & gamma < 0) ## ---- eval=FALSE-------------------------------------------------- # # Função densidade na parametrização de modelo de regressão. # dpg1 <- function(y, lambda, alpha) { # k <- exp(lfactorial(y)) # w <- 1 + alpha * y # z <- 1 + alpha * lambda # fy <- (lambda/z)^(y) * w^(y - 1) * exp(-lambda * (w/z))/k # return(fy) # } # MRDCr::dpgnz # # react <- function(panel){ # with(panel, ... ... @@ -138,46 +111,40 @@ subset(grid, round(sum, 3) != 1 & theta > 6 & gamma < 0) # 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 { # pz <- 0 # } # if (any(!is.finite(py))) { # plot(x = NULL, y = NULL, # xlim = c(0, 1), ylim = c(0, 1), # axes = FALSE, ann = FALSE) # text(x = 0.5, y = 0.5, labels = "Valores não finitos") # } else { # plot(py ~ y, type = "h", # ylim = c(0, max(c(py, pz))), # 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) { # abline(v = m, col = 2) # } # if (POIS) { # lines(y + 0.3, pz, type = "h", col = 3) # } # ylab = expression(f(y))) # mtext(side = 3, # text = substitute(sum(f(y)) == s, # list(s = round(sum(py), 5)))) # if (EX) { # abline(v = m, col = 2) # } # if (POIS) { # lines(y + 0.3, pz, type = "h", col = 3) # } # }) # panel # } # # panel <- rp.control(title = "Poisson Generalizada", # size = c(300, 100), YMAX = 30) # size = c(300, 100), YMAX = 150) # rp.slider(panel = panel, action = react, # variable = LAMBDA, title = "lambda", # from = 0.1, to = 30, # from = 0.1, to = 100, # initval = 5, resolution = 0.1, # showvalue = TRUE) # rp.slider(panel = panel, action = react, # variable = ALPHA, title = "alpha", # from = -0.2, to = 0.4, # from = -0.1, to = 0.4, # initval = 0, resolution = 0.01, # showvalue = TRUE) # rp.checkbox(panel = panel, ... ... @@ -189,36 +156,63 @@ subset(grid, round(sum, 3) != 1 & theta > 6 & gamma < 0) # rp.do(panel = panel, action = react) ## ----------------------------------------------------------------- library(bbmle) library(corrplot) #----------------------------------------------------------------------- # Gráfico do espaço paramétrico de theta x gamma. fun <- Vectorize(vectorize.args = c("theta", "gamma"), FUN = function(theta, gamma) { sum(dpgnz0(y = y, theta = theta, gamma = gamma)) }) grid <- list(theta = seq(1, 50, by = 1), gamma = seq(-0.5, 1, by = 0.05)) str(grid) y <- 0:500 my <- max(y) grid$sum <- with(grid, outer(theta, gamma, fun)) grid <- with(grid, cbind(expand.grid(theta = theta, gamma = gamma), data.frame(sum = c(sum)))) levelplot(sum ~ theta + gamma, data = subset(grid, round(sum, 3) == 1), col.regions = gray.colors) + layer(panel.abline(a = 0, b = -1/my)) + layer(panel.abline(h = 0, lty = 2)) #----------------------------------------------------------------------- # Gráfico do espaço paramétrico de lambda x alpha. fun <- Vectorize(vectorize.args = c("lambda", "alpha"), FUN = function(lambda, alpha) { sum(dpgnz(y = y, lambda = lambda, alpha = alpha)) }) grid <- list(lambda = seq(0.2, 50, by = 0.2), alpha = seq(-0.98, 0.98, by = 0.02)) grid$sum <- with(grid, outer(lambda, alpha, fun)) 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, lty = 2)) + layer(panel.curve(-1/x)) # Já que lambda * alpha > -1, então alpha = -1/lambda dá a fronteira. ## ----------------------------------------------------------------- # Função de log-Verossimilhança da Poisson Generalizada na # parametrização de modelo de regressão. PG_ll <- 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. # Gerando uma amostra aleatória da Poisson, para teste. # Offset = 2, lambda = 3. y <- rpois(100, lambda = 2 * 3) ... ... @@ -228,10 +222,10 @@ L <- list(y = y, X = cbind(rep(1, length(y)))) start <- c(alpha = 0, lambda = 1) parnames(PG_ll) <- names(start) parnames(llpgnz) <- names(start) # Como \alpha foi fixado em 1, essa ll corresponde à Poisson. n0 <- mle2(minuslogl = PG_ll, n0 <- mle2(minuslogl = llpgnz, start = start, data = L, fixed = list(alpha = 0), vecpar = TRUE) ... ... @@ -240,7 +234,7 @@ c(coef(n0)["lambda"], coef(glm(y ~ offset(log(L$offset)), family = poisson))) # Estimando o \alpha. n1 <- mle2(PG_ll, start = start, data = L, vecpar = TRUE) n1 <- mle2(llpgnz, start = start, data = L, vecpar = TRUE) coef(n1) # Perfil de verossimilhança dos parâmetros. ... ... @@ -248,11 +242,14 @@ plot(profile(n1)) # Covariância. V <- cov2cor(vcov(n1)) corrplot.mixed(V, upper = "ellipse", col = "gray50") 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() ## ----------------------------------------------------------------- library(lattice) #----------------------------------------------------------------------- # Carregando e explorando os dados. data(soja, package = "MRDCr") str(soja) ... ... @@ -262,6 +259,8 @@ 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)) ... ... @@ -270,26 +269,29 @@ 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") summary(m0) # 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))) L <- with(soja, list(y = nvag, offset = 1, X = model.matrix(m0))) # Usa as estimativas do Poisson como valore iniciais. # Usa as estimativas do Poisson como valores iniciais para a PGen. start <- c(alpha = 0, coef(m0)) parnames(PG_ll) <- names(start) parnames(llpgnz) <- names(start) # Com alpha fixo em 0 corresponde à Poisson. m2 <- mle2(PG_ll, start = start, data = L, m2 <- mle2(llpgnz, start = start, data = L, fixed = list(alpha = 0), vecpar = TRUE) # Mesma medida de ajuste e estimativas. ... ... @@ -297,23 +299,12 @@ c(logLik(m2), logLik(m0)) cbind(coef(m2)[-1], coef(m0)) # Poisson Generalizada. m3 <- mle2(PG_ll, 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) # est_stderr <- function(tb) { # sprintf("%s (%0.4f)", # formatC(tb[, 1], flag = " ", digits = 4, format = "f"), # tb[, 2]) # } # # L <- list("PoissonGLM" = rbind(NA, summary(m0)$coef), # "PoissonML*" = rbind(NA, summary(m2)@coef), # "PGeneraliz" = summary(m3)@coef) # # as.data.frame(sapply(L, est_stderr)) # Estimativas dos coeficientes. cbind("PoissonGLM" = c(NA, coef(m0)), "PoissonML" = coef(m2), "PGeneraliz" = coef(m3)) ... ... @@ -323,14 +314,17 @@ plot(profile(m3, which = "alpha")) abline(v = 0, lty = 2) V <- cov2cor(vcov(m3)) corrplot.mixed(V, upper = "ellipse", col = "gray50") 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() library(plyr) # 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) ... ... @@ -345,9 +339,8 @@ crossprod(L %*% coef(m3), solve(L %*% vcov(m3) %*% t(L), L %*% coef(m3))) library(car) # Teste de Wald para interação (poderia ser LRT, claro). # É necessário um objeto glm. # É necessário passar um objeto glm mesmo fornecendo o restante a parte. linearHypothesis(model = m0, hypothesis.matrix = L, vcov. = vcov(m3), ... ... @@ -356,9 +349,6 @@ linearHypothesis(model = m0, #----------------------------------------------------------------------- # Predição com bandas de confiança. library(doBy) library(multcomp) X <- LSmatrix(m0, effect = c("umid", "K")) pred <- attr(X, "grid") ... ... @@ -383,22 +373,18 @@ str(pred$pois) # 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$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, umid, K, modelo) str(pred) key <- list(type = "o", divide = 1, lines = list(pch = 1:nlevels(pred$modelo), ... ... @@ -409,6 +395,8 @@ 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), ... ... @@ -416,63 +404,64 @@ xyplot(fit ~ K | umid, data = pred, ## ----------------------------------------------------------------- #----------------------------------------------------------------------- # 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) 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(c(0:1), f = 0.15), xlab = "Nível de desfolhas artificial", ylab = "Número de capulho produzidos", spread = 0.035, panel = panel.beeswarm) p1 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")) #----------------------------------------------------------------------- # Modelo Poisson. m0 <- glm(ncap ~ est * (des + I(des^2)), data = capdesfo, family = poisson) m0 <- glm(ngra ~ bloc + umid * K, data = soja, family = poisson) m1 <- update(m0, family = quasipoisson) # Diagnóstico. par(mfrow = c(2, 2)) plot(m0); layout(1) anova(m0, test = "Chisq") summary(m0) # Medidas de decisão. # anova(m0, test = "Chisq") anova(m1, test = "Chisq") summary(m1) #----------------------------------------------------------------------- # Modelo Poisson Generalizada. # Modelo Poisson Generalizado. L <- with(capdesfo, list(y = ncap, offset = 1, X = model.matrix(m0))) L <- with(soja, list(y = ngra, offset = 1, X = model.matrix(m0))) start <- c(alpha = log(1), coef(m0)) parnames(PG_ll) <- names(start) # Usa as estimativas do Poisson como valores iniciais. start <- c(alpha = 0, coef(m0)) parnames(llpgnz) <- names(start) # Modelo Poisson também. m2 <- mle2(PG_ll, start = start, data = L, # Com alpha fixo em 0 corresponde à Poisson. m2 <- mle2(llpgnz, start = start, data = L, fixed = list(alpha = 0), vecpar = TRUE) # Mesma medida de ajuste e estimativas. c(logLik(m2), logLik(m0)) # Modelo Poisson Generalizado. m3 <- mle2(PG_ll, start = start, data = L, vecpar = TRUE) logLik(m3) # Poisson Generalizada. 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) summary(m3) plot(profile(m3, which = "alpha")) # Estimaitvas dos parâmetros. cbind("PoissonGLM" = c(NA, coef(m0)), "PoissonML" = coef(m2), "PGeneraliz" = coef(m3)) # 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, upper = "ellipse", col = "gray50") 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. ... ... @@ -484,16 +473,11 @@ ai <- a == max(a) L <- t(replicate(sum(ai), rbind(coef(m3) * 0), simplify = "matrix")) L[, ai] <- diag(sum(ai)) # Teste de Wald explicito. # t(L %*% coef(m3)) %*% # solve(L %*% vcov(m3) %*% t(L)) %*% # (L %*% coef(m3)) # 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 um objeto glm, mas necesse caso ele não usado para nada. linearHypothesis(model = m0, hypothesis.matrix = L, vcov. = vcov(m3), ... ... @@ -502,11 +486,13 @@ linearHypothesis(model = m0, #----------------------------------------------------------------------- # 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) X <- LSmatrix(m0, effect = c("umid", "K")) pred <- list(pois = pred, pgen = pred) pred <- attr(X, "grid") pred <- transform(pred, K = as.integer(K), umid = factor(umid)) pred <- list(pois = pred, quasi = pred, pgen = pred)