Commit a37a977e authored by Eduardo E. R. Junior's avatar Eduardo E. R. Junior

Adiciona vignette sobre excesso de zeros com análise de 'peixe'

parent 3a226d48
Pipeline #4488 failed with stage
## ----setup, include=FALSE------------------------------------------------
source("_setup.R")
## ------------------------------------------------------------------------
library(MRDCr)
## ------------------------------------------------------------------------
data(peixe)
peixe$campista <- as.factor(peixe$campista)
levels(peixe$campista) <- c("Não", "Sim")
str(peixe)
## help(peixe)
## ------------------------------------------------------------------------
## Estudo observacional
ftable(with(peixe, table(npessoas, ncriancas, campista)))
## Resumo das variáveis
summary(peixe)
## Resumo das variáveis, considerando somente as respostas não nulas
summary(subset(peixe, npeixes > 0))
## Contagens (marginal aos efeitos das covariáveis)
p1 <- histogram(~npeixes, data = peixe, nint = 50, grid = TRUE)
p2 <- histogram(~npeixes, data = subset(peixe, npeixes > 0),
nint = 50)
print(p1, split = c(1, 1, 2, 1), more = TRUE)
print(p2, split = c(2, 1, 2, 1))
## Proporção dos valores observados
(proptb <- cbind("Proporção" = prop.table(table(peixe$npeixes)),
"N. observ" = table(peixe$npeixes)))
## Disposição das covariáveis
par(mfrow = c(1, 3))
sapply(1:3, function(x) {
barplot(table(peixe[, x]))
title(main = names(peixe)[x])
names(peixe)[x]
})
## Contagens com relação as covariáveis
peixe$lnpeixes <- log(peixe$npeixes + 0.5)
xyplot(lnpeixes ~ ncriancas + npessoas,
groups = campista,
auto.key = list(
columns = 2, cex.title = 1, lines = TRUE,
title = "Presença de campista"),
data = peixe,
jitter.x = TRUE,
type = c("p", "g", "spline"))
## ------------------------------------------------------------------------
##======================================================================
## Ajuste de modelos hurdle
library(pscl)
## Preditores
f1 <- npeixes ~ campista + npessoas + ncriancas
f2 <- npeixes ~ campista * npessoas + ncriancas
## Poisson
m1P <- glm(f1, data = peixe, family = poisson)
m2P <- glm(f2, data = peixe, family = poisson)
## Hurdle Poisson
m1HP <- hurdle(f1, data = peixe, dist = "poisson")
m2HP <- hurdle(f2, data = peixe, dist = "poisson")
## Zero Inflated Poisson
m1ZP <- zeroinfl(f1, data = peixe, dist = "poisson")
m2ZP <- zeroinfl(f2, data = peixe, dist = "poisson")
## Binomial Negativa
library(MASS)
m1BN <- glm.nb(f1, data = peixe)
m2BN <- glm.nb(f2, data = peixe)
## Hurdle Binomial Negativa
m1HBN <- hurdle(f1, data = peixe, dist = "negbin")
m2HBN <- hurdle(f2, data = peixe, dist = "negbin")
## Zero Inflated Binomial Negativa
m1ZBN <- zeroinfl(f1, data = peixe, dist = "negbin")
m2ZBN <- zeroinfl(f2, data = peixe, dist = "negbin")
## ------------------------------------------------------------------------
## Via log-verossimilhança
cbind("Poisson" = sapply(list(m1P, m2P), logLik),
"HUPoisson" = sapply(list(m1HP, m2HP), logLik),
"ZIPoisson" = sapply(list(m1ZP, m2ZP), logLik),
"BinNeg" = sapply(list(m1BN, m2BN), logLik),
"HUBinNeg" = sapply(list(m1HBN, m2HBN), logLik),
"ZIBinNeg" = sapply(list(m1ZBN, m2ZBN), logLik)
)
## ------------------------------------------------------------------------
## Testes de razão de verossimilhanças para o efeito de interação
anova(m1BN, m2BN)
lmtest::lrtest(m1HBN, m2HBN)
lmtest::lrtest(m1ZBN, m2ZBN)
## ------------------------------------------------------------------------
## Teste de Vuong para diferença entre os modelos BN e HUBN
vuong(m1BN, m1HBN)
## Teste de Vuong para diferença entre os modelos ZIBN e HUBN
vuong(m1ZBN, m1HBN)
## ------------------------------------------------------------------------
## Estimativas dos parâmetros e testes de Wald
summary(m1BN)
summary(m1HBN)
summary(m1ZBN)
## ------------------------------------------------------------------------
## Ajuste de preditor do modelo proposto
## Retira o efeito de campista no preditor para as contagens não nula
m3HBN <- hurdle(npeixes ~ npessoas + ncriancas |
campista + npessoas + ncriancas,
data = peixe, dist = "negbin")
lmtest::lrtest(m1HBN, m3HBN)
## Refazendo o teste de Vuong para comparação de modelos BN e HUBN
vuong(m1BN, m3HBN)
## ----diag, cache = TRUE--------------------------------------------------
## Uma pequena avaliação dos resíduos
p1 <- xyplot(residuals(m3HBN) ~ fitted(m3HBN),
type = c("p", "g", "spline"))
p2 <- qqmath(~residuals(m3HBN, type = "pearson"),
aspect = "xy",
type = c("p", "g"),
panel = function(x, ...) {
panel.qqmathline(x, ...)
panel.qqmath(x, ...)
})
qqsimul <- hnp::hnp(m3HBN, plot = FALSE)
p3 <- with(qqsimul,
xyplot(residuals ~ x, pch = 20) +
as.layer(
xyplot(median + lower + upper ~ x,
type = "l", col = 1, lty = c(2, 1, 1)))
)
print(p1, split = c(1, 1, 3, 1), more = TRUE)
print(p2, split = c(2, 1, 3, 1), more = TRUE)
print(p3, split = c(3, 1, 3, 1), more = FALSE)
## ------------------------------------------------------------------------
## Estimativas dos parâmetros e testes de Wald
summary(m3HBN)
## ------------------------------------------------------------------------
##----------------------------------------------------------------------
## Região para predição
pred <- expand.grid(campista = c("Não", "Sim"),
ncriancas = 0:3, npessoas = 1:4)
##-------------------------------------------
## Estimando as médias
aux <- predict(m3HBN, newdata = pred, type = "response")
predmu <- cbind(pred, fit = aux)
xyplot(fit ~ npessoas | campista,
groups = ncriancas, data = predmu,
type = c("l", "g"),
auto.key = list(
columns = 2, cex.title = 1,
lines = TRUE, points = FALSE,
title = "Número de crianças"),
strip = strip.custom(
strip.names = TRUE, var.name = "campista"
))
##-------------------------------------------
## Estimando probabilidades
pred <- expand.grid(campista = "Sim",
ncriancas = 0:3, npessoas = 4)
py <- c(t(predict(m3HBN, type = "prob", at = 0:20, newdata = pred)))
da <- data.frame(y = 0:20, py = py, ncriancas = rep(0:3, each = 21))
barchart(py ~ y | factor(ncriancas), data = da,
stack = FALSE, horizontal = FALSE,
as.table = TRUE, between = list(y = 0.5),
scales = list( y = "free", x = list(labels = 0:20)),
strip = strip.custom(
strip.names = TRUE, var.name = "nº de crianças"
))
## ----boot, cache = TRUE--------------------------------------------------
## Intervalos de confiança para predição
##-------------------------------------------
## Via reamostragem
n <- nrow(peixe)
formula <- m3HBN$formula
start <- list(zero = coef(m3HBN, "zero"), count = coef(m3HBN, "count"))
boots <- replicate(100, {
id <- sample(1:n, replace = TRUE)
model <- hurdle(formula, data = peixe[id, ], start = start)
yhat <- predict(model, newdata = predmu, type = "response")
})
pred2 <- cbind(predmu, t(apply(boots, 1, function(x) {
quantile(x, probs = c(0.025, 0.975))})))
names(pred2)[5:6] <- c("lwr", "upr")
## Viasualizando graficamente
xyplot(fit ~ npessoas | campista,
groups = ncriancas, data = pred2,
type = c("l", "g"), cty = "bands",
ly = pred2$lwr, uy = pred2$upr,
prepanel = prepanel.cbH, alpha = 0.5,
panel = panel.superpose,
panel.groups = panel.cbH,
auto.key = list(
columns = 2, cex.title = 1,
lines = TRUE, points = FALSE,
title = "Número de crianças"),
strip = strip.custom(
strip.names = TRUE, var.name = "campista"
))
## ------------------------------------------------------------------------
## Modelo Hurdle (binomial e binomial negativa)
m3HBN$formula
##-------------------------------------------
## Componente da contagem nula (f_zero)
indica <- with(peixe, ifelse(npeixes == 0, 1, 0))
mzero <- glm(indica ~ campista + npessoas + ncriancas,
family = binomial, data = peixe)
## Comparando os coeficientes
cbind("glm_binomial" = coef(mzero),
"zeroinfl" = coef(m3HBN, model = "zero"))
##-------------------------------------------
## Componente da contagem nula (f_count)
library(VGAM)
countp <- subset(peixe, npeixes > 0)
mcount <- vglm(npeixes ~ npessoas + ncriancas,
family = posnegbinomial, data = countp)
## Comparando os coeficientes (betas e theta (da BN))
cbind("vglm_posnegbin" = coef(mcount)[-2],
"zeroinfl" = coef(m3HBN, model = "count"))
cbind("vglm_posnegbin" = exp(coef(mcount)[2]),
"zeroinfl" = m3HBN$theta)
This diff is collapsed.
---
title: "Análise de Contagens Modelando Excesso de Zeros"
author: >
Walmes M. Zeviani,
Eduardo E. Ribeiro Jr &
Cesar A. Taconeli
vignette: >
%\VignetteIndexEntry{Análise de Contagens Modelando Excesso de Zeros}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include=FALSE}
source("_setup.R")
```
```{r}
library(MRDCr)
```
# Números de peixes capturados em um Parque Estadual #
```{r}
data(peixe)
peixe$campista <- as.factor(peixe$campista)
levels(peixe$campista) <- c("Não", "Sim")
str(peixe)
## help(peixe)
```
Dados sobre 250 grupos que foram ao Parque Estadual para pescar. As
informações coletadas foram refentes a presenção ou não de um campista,
ao número de crianças no grupo e ao número de indivíduos no grupo. Assim
as variáveis definidas são:
`campista`: Fator com dois níveis que representa a presença (_Sim_) ou
ausência (_Não_) de um campista no grupo.
`ncriancas`: Número de crianças no grupo.
`npessoas`: Número total de pessoas no grupo.
`npeixes`: Número de peixes capturados pelo grupo.
## Análise exploratória ##
```{r}
## Estudo observacional
ftable(with(peixe, table(npessoas, ncriancas, campista)))
## Resumo das variáveis
summary(peixe)
## Resumo das variáveis, considerando somente as respostas não nulas
summary(subset(peixe, npeixes > 0))
## Contagens (marginal aos efeitos das covariáveis)
p1 <- histogram(~npeixes, data = peixe, nint = 50, grid = TRUE)
p2 <- histogram(~npeixes, data = subset(peixe, npeixes > 0),
nint = 50)
print(p1, split = c(1, 1, 2, 1), more = TRUE)
print(p2, split = c(2, 1, 2, 1))
## Proporção dos valores observados
(proptb <- cbind("Proporção" = prop.table(table(peixe$npeixes)),
"N. observ" = table(peixe$npeixes)))
## Disposição das covariáveis
par(mfrow = c(1, 3))
sapply(1:3, function(x) {
barplot(table(peixe[, x]))
title(main = names(peixe)[x])
names(peixe)[x]
})
## Contagens com relação as covariáveis
peixe$lnpeixes <- log(peixe$npeixes + 0.5)
xyplot(lnpeixes ~ ncriancas + npessoas,
groups = campista,
auto.key = list(
columns = 2, cex.title = 1, lines = TRUE,
title = "Presença de campista"),
data = peixe,
jitter.x = TRUE,
type = c("p", "g", "spline"))
```
## Ajuste de modelos ##
```{r}
##======================================================================
## Ajuste de modelos hurdle
library(pscl)
## Preditores
f1 <- npeixes ~ campista + npessoas + ncriancas
f2 <- npeixes ~ campista * npessoas + ncriancas
## Poisson
m1P <- glm(f1, data = peixe, family = poisson)
m2P <- glm(f2, data = peixe, family = poisson)
## Hurdle Poisson
m1HP <- hurdle(f1, data = peixe, dist = "poisson")
m2HP <- hurdle(f2, data = peixe, dist = "poisson")
## Zero Inflated Poisson
m1ZP <- zeroinfl(f1, data = peixe, dist = "poisson")
m2ZP <- zeroinfl(f2, data = peixe, dist = "poisson")
## Binomial Negativa
library(MASS)
m1BN <- glm.nb(f1, data = peixe)
m2BN <- glm.nb(f2, data = peixe)
## Hurdle Binomial Negativa
m1HBN <- hurdle(f1, data = peixe, dist = "negbin")
m2HBN <- hurdle(f2, data = peixe, dist = "negbin")
## Zero Inflated Binomial Negativa
m1ZBN <- zeroinfl(f1, data = peixe, dist = "negbin")
m2ZBN <- zeroinfl(f2, data = peixe, dist = "negbin")
```
## Comparação dos ajustes ##
```{r}
## Via log-verossimilhança
cbind("Poisson" = sapply(list(m1P, m2P), logLik),
"HUPoisson" = sapply(list(m1HP, m2HP), logLik),
"ZIPoisson" = sapply(list(m1ZP, m2ZP), logLik),
"BinNeg" = sapply(list(m1BN, m2BN), logLik),
"HUBinNeg" = sapply(list(m1HBN, m2HBN), logLik),
"ZIBinNeg" = sapply(list(m1ZBN, m2ZBN), logLik)
)
```
```{r}
## Testes de razão de verossimilhanças para o efeito de interação
anova(m1BN, m2BN)
lmtest::lrtest(m1HBN, m2HBN)
lmtest::lrtest(m1ZBN, m2ZBN)
```
```{r}
## Teste de Vuong para diferença entre os modelos BN e HUBN
vuong(m1BN, m1HBN)
## Teste de Vuong para diferença entre os modelos ZIBN e HUBN
vuong(m1ZBN, m1HBN)
```
```{r}
## Estimativas dos parâmetros e testes de Wald
summary(m1BN)
summary(m1HBN)
summary(m1ZBN)
```
```{r}
## Ajuste de preditor do modelo proposto
## Retira o efeito de campista no preditor para as contagens não nula
m3HBN <- hurdle(npeixes ~ npessoas + ncriancas |
campista + npessoas + ncriancas,
data = peixe, dist = "negbin")
lmtest::lrtest(m1HBN, m3HBN)
## Refazendo o teste de Vuong para comparação de modelos BN e HUBN
vuong(m1BN, m3HBN)
```
## Avaliação do modelo proposto ##
```{r diag, cache = TRUE}
## Uma pequena avaliação dos resíduos
p1 <- xyplot(residuals(m3HBN) ~ fitted(m3HBN),
type = c("p", "g", "spline"))
p2 <- qqmath(~residuals(m3HBN, type = "pearson"),
aspect = "xy",
type = c("p", "g"),
panel = function(x, ...) {
panel.qqmathline(x, ...)
panel.qqmath(x, ...)
})
qqsimul <- hnp::hnp(m3HBN, plot = FALSE)
p3 <- with(qqsimul,
xyplot(residuals ~ x, pch = 20) +
as.layer(
xyplot(median + lower + upper ~ x,
type = "l", col = 1, lty = c(2, 1, 1)))
)
print(p1, split = c(1, 1, 3, 1), more = TRUE)
print(p2, split = c(2, 1, 3, 1), more = TRUE)
print(p3, split = c(3, 1, 3, 1), more = FALSE)
```
```{r}
## Estimativas dos parâmetros e testes de Wald
summary(m3HBN)
```
## Predição ##
```{r}
##----------------------------------------------------------------------
## Região para predição
pred <- expand.grid(campista = c("Não", "Sim"),
ncriancas = 0:3, npessoas = 1:4)
##-------------------------------------------
## Estimando as médias
aux <- predict(m3HBN, newdata = pred, type = "response")
predmu <- cbind(pred, fit = aux)
xyplot(fit ~ npessoas | campista,
groups = ncriancas, data = predmu,
type = c("l", "g"),
auto.key = list(
columns = 2, cex.title = 1,
lines = TRUE, points = FALSE,
title = "Número de crianças"),
strip = strip.custom(
strip.names = TRUE, var.name = "campista"
))
##-------------------------------------------
## Estimando probabilidades
pred <- expand.grid(campista = "Sim",
ncriancas = 0:3, npessoas = 4)
py <- c(t(predict(m3HBN, type = "prob", at = 0:20, newdata = pred)))
da <- data.frame(y = 0:20, py = py, ncriancas = rep(0:3, each = 21))
barchart(py ~ y | factor(ncriancas), data = da,
stack = FALSE, horizontal = FALSE,
as.table = TRUE, between = list(y = 0.5),
scales = list( y = "free", x = list(labels = 0:20)),
strip = strip.custom(
strip.names = TRUE, var.name = "nº de crianças"
))
```
```{r boot, cache = TRUE}
## Intervalos de confiança para predição
##-------------------------------------------
## Via reamostragem
n <- nrow(peixe)
formula <- m3HBN$formula
start <- list(zero = coef(m3HBN, "zero"), count = coef(m3HBN, "count"))
boots <- replicate(100, {
id <- sample(1:n, replace = TRUE)
model <- hurdle(formula, data = peixe[id, ], start = start)
yhat <- predict(model, newdata = predmu, type = "response")
})
pred2 <- cbind(predmu, t(apply(boots, 1, function(x) {
quantile(x, probs = c(0.025, 0.975))})))
names(pred2)[5:6] <- c("lwr", "upr")
## Viasualizando graficamente
xyplot(fit ~ npessoas | campista,
groups = ncriancas, data = pred2,
type = c("l", "g"), cty = "bands",
ly = pred2$lwr, uy = pred2$upr,
prepanel = prepanel.cbH, alpha = 0.5,
panel = panel.superpose,
panel.groups = panel.cbH,
auto.key = list(
columns = 2, cex.title = 1,
lines = TRUE, points = FALSE,
title = "Número de crianças"),
strip = strip.custom(
strip.names = TRUE, var.name = "campista"
))
```
## Ajuste em separado das componentes do modelo Hurdle ##
```{r}
## Modelo Hurdle (binomial e binomial negativa)
m3HBN$formula
##-------------------------------------------
## Componente da contagem nula (f_zero)
indica <- with(peixe, ifelse(npeixes == 0, 1, 0))
mzero <- glm(indica ~ campista + npessoas + ncriancas,
family = binomial, data = peixe)
## Comparando os coeficientes
cbind("glm_binomial" = coef(mzero),
"zeroinfl" = coef(m3HBN, model = "zero"))
##-------------------------------------------
## Componente da contagem nula (f_count)
library(VGAM)
countp <- subset(peixe, npeixes > 0)
mcount <- vglm(npeixes ~ npessoas + ncriancas,
family = posnegbinomial, data = countp)
## Comparando os coeficientes (betas e theta (da BN))
cbind("vglm_posnegbin" = coef(mcount)[-2],
"zeroinfl" = coef(m3HBN, model = "count"))
cbind("vglm_posnegbin" = exp(coef(mcount)[2]),
"zeroinfl" = m3HBN$theta)
```
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment