|
Modelos de Regressão para Dados de Contagem com o R
|
Análise de Contagens Modelando Excesso de Zeros
Walmes M. Zeviani, Eduardo E. Ribeiro Jr & Cesar A. Taconeli
library(MRDCr)## Loading required package: bbmle## Loading required package: stats4## ----------------------------------------------------------------------
## MRDCr: Modelos de Regressão para Dados de Contagem
##
## Para colaboração, suporte ou relato de bugs, visite:
## https://gitlab.c3sl.ufpr.br/leg/MRDCr
##
## MRDCr version 0.0-2 (feito em 2016-07-21) foi carregado.
## ----------------------------------------------------------------------Números de peixes capturados em um Parque Estadual
data(peixe)
peixe$campista <- as.factor(peixe$campista)
levels(peixe$campista) <- c("Não", "Sim")
str(peixe)## 'data.frame': 250 obs. of 4 variables:
## $ campista : Factor w/ 2 levels "Não","Sim": 1 2 1 2 1 2 1 1 2 2 ...
## $ ncriancas: int 0 0 0 1 0 2 1 3 2 0 ...
## $ npessoas : int 1 1 1 2 1 4 3 4 3 1 ...
## $ npeixes : int 0 0 0 0 1 0 0 0 0 1 ...## 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
## Estudo observacional
ftable(with(peixe, table(npessoas, ncriancas, campista)))## campista Não Sim
## npessoas ncriancas
## 1 0 22 35
## 1 0 0
## 2 0 0
## 3 0 0
## 2 0 17 18
## 1 12 23
## 2 0 0
## 3 0 0
## 3 0 8 13
## 1 8 11
## 2 5 12
## 3 0 0
## 4 0 6 13
## 1 10 11
## 2 11 5
## 3 4 6## Resumo das variáveis
summary(peixe)## campista ncriancas npessoas npeixes
## Não:103 Min. :0.000 Min. :1.000 Min. : 0.000
## Sim:147 1st Qu.:0.000 1st Qu.:2.000 1st Qu.: 0.000
## Median :0.000 Median :2.000 Median : 0.000
## Mean :0.684 Mean :2.528 Mean : 3.296
## 3rd Qu.:1.000 3rd Qu.:4.000 3rd Qu.: 2.000
## Max. :3.000 Max. :4.000 Max. :149.000## Resumo das variáveis, considerando somente as respostas não nulas
summary(subset(peixe, npeixes > 0))## campista ncriancas npessoas npeixes
## Não:34 Min. :0.0000 Min. :1.000 Min. : 1.00
## Sim:74 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.: 1.00
## Median :0.0000 Median :3.000 Median : 3.00
## Mean :0.3241 Mean :2.667 Mean : 7.63
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.: 6.00
## Max. :2.0000 Max. :4.000 Max. :149.00## 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)))## Proporção N. observ
## 0 0.568 142
## 1 0.124 31
## 2 0.080 20
## 3 0.048 12
## 4 0.024 6
## 5 0.040 10
## 6 0.016 4
## 7 0.012 3
## 8 0.008 2
## 9 0.008 2
## 10 0.004 1
## 11 0.004 1
## 13 0.004 1
## 14 0.004 1
## 15 0.008 2
## 16 0.004 1
## 21 0.008 2
## 22 0.004 1
## 29 0.004 1
## 30 0.004 1
## 31 0.004 1
## 32 0.008 2
## 38 0.004 1
## 65 0.004 1
## 149 0.004 1## 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]
})
## [1] "campista" "ncriancas" "npessoas"## 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
##======================================================================
## Ajuste de modelos hurdle
library(pscl)## Loading required package: MASS## Classes and Methods for R developed in the## Political Science Computational Laboratory## Department of Political Science## Stanford University## Simon Jackman## hurdle and zeroinfl functions by Achim Zeileis## 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
## 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)
)## Poisson HUPoisson ZIPoisson BinNeg HUBinNeg ZIBinNeg
## [1,] -837.0725 -751.6182 -752.7315 -405.2220 -395.1590 -395.5394
## [2,] -830.3639 -743.4591 -744.0157 -404.9135 -393.5848 -394.1446## Testes de razão de verossimilhanças para o efeito de interação
anova(m1BN, m2BN)## Likelihood ratio tests of Negative Binomial Models
##
## Response: npeixes
## Model theta Resid. df
## 1 campista + npessoas + ncriancas 0.4635288 246
## 2 campista * npessoas + ncriancas 0.4672296 245
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -810.444
## 2 -809.827 1 vs 2 1 0.6170103 0.4321604lmtest::lrtest(m1HBN, m2HBN)## Likelihood ratio test
##
## Model 1: npeixes ~ campista + npessoas + ncriancas
## Model 2: npeixes ~ campista * npessoas + ncriancas
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -395.16
## 2 11 -393.58 2 3.1484 0.2072lmtest::lrtest(m1ZBN, m2ZBN)## Likelihood ratio test
##
## Model 1: npeixes ~ campista + npessoas + ncriancas
## Model 2: npeixes ~ campista * npessoas + ncriancas
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -395.54
## 2 11 -394.14 2 2.7896 0.2479## Teste de Vuong para diferença entre os modelos BN e HUBN
vuong(m1BN, m1HBN)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -1.5078524 model2 > model1 0.065796
## AIC-corrected -0.9084884 model2 > model1 0.181810
## BIC-corrected 0.1468301 model1 > model2 0.441633## Teste de Vuong para diferença entre os modelos ZIBN e HUBN
vuong(m1ZBN, m1HBN)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -0.09895803 model2 > model1 0.46059
## AIC-corrected -0.09895803 model2 > model1 0.46059
## BIC-corrected -0.09895803 model2 > model1 0.46059## Estimativas dos parâmetros e testes de Wald
summary(m1BN)##
## Call:
## glm.nb(formula = f1, data = peixe, init.theta = 0.4635287626,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6673 -0.9599 -0.6590 -0.0319 4.9433
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6250 0.3304 -4.918 8.74e-07 ***
## campistaSim 0.6211 0.2348 2.645 0.00816 **
## npessoas 1.0608 0.1144 9.273 < 2e-16 ***
## ncriancas -1.7805 0.1850 -9.623 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4635) family taken to be 1)
##
## Null deviance: 394.25 on 249 degrees of freedom
## Residual deviance: 210.65 on 246 degrees of freedom
## AIC: 820.44
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4635
## Std. Err.: 0.0712
##
## 2 x log-likelihood: -810.4440summary(m1HBN)##
## Call:
## hurdle(formula = f1, data = peixe, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.72774 -0.51588 -0.25938 -0.03287 14.70973
##
## Count model coefficients (truncated negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6215 0.5960 -2.720 0.006521 **
## campistaSim 0.3746 0.3360 1.115 0.264934
## npessoas 1.0029 0.1551 6.465 1.01e-10 ***
## ncriancas -1.0945 0.3198 -3.422 0.000621 ***
## Log(theta) -1.0530 0.4974 -2.117 0.034266 *
## Zero hurdle model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3087 0.4612 -5.005 5.57e-07 ***
## campistaSim 1.0179 0.3246 3.136 0.00171 **
## npessoas 1.1104 0.1911 5.811 6.19e-09 ***
## ncriancas -2.1380 0.3107 -6.882 5.90e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta: count = 0.3489
## Number of iterations in BFGS optimization: 13
## Log-likelihood: -395.2 on 9 Dfsummary(m1ZBN)##
## Call:
## zeroinfl(formula = f1, data = peixe, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.71806 -0.56103 -0.38168 0.04398 16.16367
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6177 0.3202 -5.052 4.37e-07 ***
## campistaSim 0.3856 0.2461 1.567 0.117128
## npessoas 1.0901 0.1117 9.761 < 2e-16 ***
## ncriancas -1.2613 0.2473 -5.100 3.40e-07 ***
## Log(theta) -0.5929 0.1580 -3.753 0.000174 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.9920 64.4408 -0.186 0.852
## campistaSim -10.7704 64.3725 -0.167 0.867
## npessoas 0.2902 0.7314 0.397 0.692
## ncriancas 10.9517 64.3569 0.170 0.865
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 0.5527
## Number of iterations in BFGS optimization: 68
## Log-likelihood: -395.5 on 9 Df## 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)## Likelihood ratio test
##
## Model 1: npeixes ~ campista + npessoas + ncriancas
## Model 2: npeixes ~ npessoas + ncriancas | campista + npessoas + ncriancas
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -395.16
## 2 8 -395.75 -1 1.1744 0.2785## Refazendo o teste de Vuong para comparação de modelos BN e HUBN
vuong(m1BN, m3HBN)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -1.2007642 model2 > model1 0.11492
## AIC-corrected -0.8206083 model2 > model1 0.20593
## BIC-corrected -0.1512561 model2 > model1 0.43989Avaliação do modelo proposto
## Uma pequena avaliação dos resíduos
p1 <- xyplot(residuals(m3HBN) ~ fitted(m3HBN),
type = c("p", "g", "spline"))
p2 <- qqmath(~residuals(m3HBN, type = "pearson"),
type = c("p", "g"),
panel = function(x, ...) {
panel.qqmathline(x, ...)
panel.qqmath(x, ...)
})
qqsimul <- hnp::hnp(m3HBN, plot = FALSE)## Hurdle negative binomial modelp3 <- 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)##
## Call:
## hurdle(formula = npeixes ~ npessoas + ncriancas | campista +
## npessoas + ncriancas, data = peixe, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.72427 -0.49730 -0.25410 -0.01439 12.19296
##
## Count model coefficients (truncated negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4334 0.5719 -2.506 0.012196 *
## npessoas 1.0260 0.1552 6.612 3.8e-11 ***
## ncriancas -1.1642 0.3177 -3.665 0.000247 ***
## Log(theta) -1.1294 0.5185 -2.178 0.029387 *
## Zero hurdle model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3087 0.4612 -5.005 5.57e-07 ***
## campistaSim 1.0179 0.3246 3.136 0.00171 **
## npessoas 1.1104 0.1911 5.811 6.19e-09 ***
## ncriancas -2.1380 0.3107 -6.882 5.90e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta: count = 0.3232
## Number of iterations in BFGS optimization: 12
## Log-likelihood: -395.7 on 8 DfPredição
##----------------------------------------------------------------------
## 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"
))
## 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
## Modelo Hurdle (binomial e binomial negativa)
m3HBN$formula## npeixes ~ npessoas + ncriancas | campista + npessoas + ncriancas##-------------------------------------------
## 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"))## glm_binomial zeroinfl
## (Intercept) 2.308713 -2.308713
## campistaSim -1.017853 1.017853
## npessoas -1.110389 1.110389
## ncriancas 2.137983 -2.137983# ##-------------------------------------------
# ## 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)Número de sinistros em uma seguradora de veículos
Análise exploratória
data(seguros)
str(seguros)## 'data.frame': 16483 obs. of 7 variables:
## $ tipo : Factor w/ 6 levels "hatch","mono",..: 1 2 5 4 2 1 1 4 5 5 ...
## $ idade : int 59 45 42 63 36 33 35 63 54 32 ...
## $ sexo : Factor w/ 2 levels "F","M": 1 2 2 1 1 2 2 2 2 2 ...
## $ civil : Factor w/ 4 levels "Casado","Divorciado",..: 1 3 1 1 1 1 1 1 1 1 ...
## $ valor : num 24.6 23.4 86.6 77.5 25.9 ...
## $ expos : num 0.5 0.7 0.79 0.01 0.51 0.79 0.81 0.01 0.76 0.79 ...
## $ nsinist: int 1 0 0 0 0 0 0 0 0 0 ...## help(seguros)Dados referentes ao acompanhamento de clientes de uma seguradora de automóveis ao longo de um ano. Foram registrados as variáveis descritas abaixo para 16483 clientes.
tipo: Tipo de veículo segurado. Fator com seis níveis (hatch, mono, picape, sedan, wagon e suv).idade: Idade do cliente, em anos.sexo: Sexo do cliente. Fator com dois níveis, M para clientes do sexo masculino e F para feminino.civil: Estado civil do cliente. Fator com quatro níveis (Casado, Divorciado, Solteiro e Viuvo).valor: Valor do veículo segurado, em 1000 reais.expos: Período de cobertura do cliente durante o ano sob análise.nsinist: Número de sinistros registrados no período.
summary(seguros)## tipo idade sexo civil
## hatch :6080 Min. :18.00 F:6694 Casado :13327
## mono :1546 1st Qu.:41.00 M:9789 Divorciado: 729
## picape:1197 Median :52.00 Solteiro : 1896
## sedan :4696 Mean :52.13 Viuvo : 531
## suv :2499 3rd Qu.:62.00
## wagon : 465 Max. :96.00
## valor expos nsinist
## Min. : 1.60 Min. :0.0100 Min. :0.0000
## 1st Qu.: 22.07 1st Qu.:0.5900 1st Qu.:0.0000
## Median : 31.92 Median :0.7800 Median :0.0000
## Mean : 36.03 Mean :0.6483 Mean :0.0617
## 3rd Qu.: 45.14 3rd Qu.:0.8100 3rd Qu.:0.0000
## Max. :196.65 Max. :1.0000 Max. :5.0000Ajuste de modelos
library(pscl)
## Preditores
f0 <- nsinist ~ 1
f1 <- nsinist ~ sexo + valor + log(expos)
f2 <- nsinist ~ sexo * valor + log(expos)
## Poisson
m0P <- glm(f0, data = seguros, family = poisson)
m1P <- glm(f1, data = seguros, family = poisson)
m2P <- glm(f2, data = seguros, family = poisson)
## Hurdle Poisson
m0HP <- hurdle(f0, data = seguros, dist = "poisson")
m1HP <- hurdle(f1, data = seguros, dist = "poisson")
m2HP <- hurdle(f2, data = seguros, dist = "poisson")
## Zero Inflated Poisson
m0ZP <- zeroinfl(f0, data = seguros, dist = "poisson")
m1ZP <- zeroinfl(f1, data = seguros, dist = "poisson")
m2ZP <- zeroinfl(f2, data = seguros, dist = "poisson")
## Binomial Negativa
library(MASS)
m0BN <- glm.nb(f0, data = seguros)
m1BN <- glm.nb(f1, data = seguros)
m2BN <- glm.nb(f2, data = seguros)
## Hurdle Binomial Negativa
m0HBN <- hurdle(f0, data = seguros, dist = "negbin")
m1HBN <- hurdle(f1, data = seguros, dist = "negbin")
m2HBN <- hurdle(f2, data = seguros, dist = "negbin")
## Zero Inflated Poisson
m0ZBN <- zeroinfl(f0, data = seguros, dist = "negbin")
m1ZBN <- zeroinfl(f1, data = seguros, dist = "negbin")
m2ZBN <- zeroinfl(f2, data = seguros, dist = "negbin")Comparação dos ajustes
## Via log-verossimilhança
cbind("Poisson" = sapply(list(m0P, m1P, m2P), logLik),
"HUPoisson" = sapply(list(m0HP, m1HP, m2HP), logLik),
"ZIPoisson" = sapply(list(m0ZP, m1ZP, m2ZP), logLik),
"BinNeg" = sapply(list(m0BN, m1BN, m2BN), logLik),
"HUBinNeg" = sapply(list(m0HBN, m1HBN, m2HBN), logLik),
"ZIBinNeg" = sapply(list(m0ZBN, m1ZBN, m2ZBN), logLik)
)## Poisson HUPoisson ZIPoisson BinNeg HUBinNeg ZIBinNeg
## [1,] -4018.628 -3712.755 -3712.755 -3723.977 -3712.756 -3712.756
## [2,] -3998.623 -3696.327 -3696.251 -3710.146 -3696.327 -3696.251
## [3,] -3998.609 -3695.106 -3694.896 -3710.140 -3695.106 -3694.896## Testes de razão de verossimilhanças
lmtest::lrtest(m0HP, m1HP, m2HP)## Likelihood ratio test
##
## Model 1: nsinist ~ 1
## Model 2: nsinist ~ sexo + valor + log(expos)
## Model 3: nsinist ~ sexo * valor + log(expos)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 2 -3712.8
## 2 8 -3696.3 6 32.8576 1.117e-05 ***
## 3 10 -3695.1 2 2.4414 0.295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::lrtest(m0HBN, m1HBN, m2HBN)## Likelihood ratio test
##
## Model 1: nsinist ~ 1
## Model 2: nsinist ~ sexo + valor + log(expos)
## Model 3: nsinist ~ sexo * valor + log(expos)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -3712.8
## 2 9 -3696.3 6 32.8579 1.117e-05 ***
## 3 11 -3695.1 2 2.4406 0.2951
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Para a componente de contagens não negativas é necessário um modelo
## Binomial Negativa, considerando superdispersão dos dados?
vuong(m1HP, m1HBN)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw 0.2024475 model1 > model2 0.41978
## AIC-corrected 0.2024475 model1 > model2 0.41978
## BIC-corrected 0.2024475 model1 > model2 0.41978Avaliação do modelo proposto
## Estimativas do modelo proposto
summary(m1HP)##
## Call:
## hurdle(formula = f1, data = seguros, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.2981 -0.2157 -0.2059 -0.1962 13.9964
##
## Count model coefficients (truncated poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.463668 0.161143 -2.877 0.00401 **
## sexoM -0.207454 0.129131 -1.607 0.10815
## valor -0.001857 0.003340 -0.556 0.57822
## log(expos) 0.030484 0.097958 0.311 0.75565
## Zero hurdle model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.969721 0.086905 -34.172 < 2e-16 ***
## sexoM -0.148773 0.073399 -2.027 0.04267 *
## valor 0.005012 0.001721 2.913 0.00358 **
## log(expos) 0.186904 0.046857 3.989 6.64e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 17
## Log-likelihood: -3696 on 8 Df## Retira efeito de sexo na componente das contagens nulas
m3HP <- hurdle(nsinist ~ 1 | sexo + valor + log(expos),
data = seguros)
lmtest::lrtest(m1HP, m3HP)## Likelihood ratio test
##
## Model 1: nsinist ~ sexo + valor + log(expos)
## Model 2: nsinist ~ 1 | sexo + valor + log(expos)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 8 -3696.3
## 2 5 -3697.9 -3 3.095 0.3772summary(m3HP)##
## Call:
## hurdle(formula = nsinist ~ 1 | sexo + valor + log(expos),
## data = seguros)
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.2965 -0.2163 -0.2055 -0.1958 14.2575
##
## Count model coefficients (truncated poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.65911 0.06449 -10.22 <2e-16 ***
## Zero hurdle model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.969721 0.086905 -34.172 < 2e-16 ***
## sexoM -0.148773 0.073399 -2.027 0.04267 *
## valor 0.005012 0.001721 2.913 0.00358 **
## log(expos) 0.186904 0.046857 3.989 6.64e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 10
## Log-likelihood: -3698 on 5 Df## Avaliação de diferentes especificações para a componente de contagens
## nulas
m3HP.pois <- hurdle(nsinist ~ 1 | sexo + valor + log(expos),
data = seguros, zero.dist = "poisson")
m3HP.negb <- hurdle(nsinist ~ 1 | sexo + valor + log(expos),
data = seguros, zero.dist = "negbin")
m3HP.geom <- hurdle(nsinist ~ 1 | sexo + valor + log(expos),
data = seguros, zero.dist = "geometric")
## Comparação das funções de ligação
sapply(list("binomial" = m3HP, "poisson" = m3HP.pois,
"negbin" = m3HP.negb, "geometric" = m3HP.geom), logLik)## binomial poisson negbin geometric
## -3697.874 -3697.899 -3696.887 -3697.874Predição
## Região para predição
pred <- expand.grid(sexo = c("M", "F"),
valor = 1:200,
expos = c(0.1, 0.5, 1))
##-------------------------------------------
## Estimando as médias
aux <- predict(m3HP, newdata = pred, type = "response")
predmu <- cbind(pred, fit = aux)
xyplot(fit ~ valor | factor(expos),
layout = c(NA, 1),
groups = sexo,
data = predmu,
type = c("l", "g"),
strip = strip.custom(strip.names = TRUE,
var.name = "Exposição"),
auto.key = list(
columns = 2, cex.title = 1,
lines = TRUE, points = FALSE,
title = "Número de crianças"))
##-------------------------------------------
## Estimando probabilidades
pred <- expand.grid(sexo = c("M", "F"),
valor = c(1, 50),
expos = 0.5)
py <- c(t(predict(m3HP, type = "prob", at = 0:5, newdata = pred)))
carac <- with(pred, paste("Sexo:", sexo, "Valor:", valor))
da <- data.frame(y = 0:5, py = py,
sexo = rep(c("M", "F"), each = 6),
valor = rep(c(1, 50), each = 12))
useOuterStrips(
barchart(py ~ y | sexo + factor(valor), data = da,
stack = FALSE, horizontal = FALSE,
as.table = TRUE, between = list(y = 0.5),
scales = list( y = "free", x = list(labels = 0:5))
),
strip = strip.custom(strip.names = TRUE, var.name = "Sexo"),
strip.left = strip.custom(strip.names = TRUE, var.name = "Valor")
)
Modelos de Regressão para Dados de Contagem com o R |
|
| 61 RBRAS | Departamento de Estatística - UFPR |
| 23 a 25 de Maio de 2016 | LEG, PET Estatística |
| Salvador - Bahia | github.com/leg-ufpr/MRDCr |