|
Modelos de Regressão para Dados de Contagem com o R
|
Análise de Contagens por Quase-Verossimilhança
Walmes M. Zeviani, Eduardo E. Ribeiro Jr & Cesar A. Taconeli
Os dados são referentes a um experimento realizado com o objetivo de investigar o efeito de uma intervenção, por parte do cuidador, no comportamento de ovelhas.
Para isso, foram consideradas ovelhas de duas linhagens distintas (pouco ou muito reativas), submetidas a dois tipos diferentes de intervenção (observação ou observação + intervenção).
A variável resposta aqui considerada é o número de mudanças na postura corporal do animal ao longo do período de observação (3 minutos).
# Pacotes requeridos.
library("lmtest")
library("boot")
library("car")
library("latticeExtra")
library("RColorBrewer")
library("sandwich")
library("hnp")
library("knitr")
library("MRDCr")Análise descritiva
data(postura)
str(postura)## 'data.frame': 38 obs. of 3 variables:
## $ trat : Factor w/ 2 levels "Observ","Observ + Interv": 1 1 1 1 2 2 2 2 2 2 ...
## $ linh : Factor w/ 2 levels "Pouco reativo",..: 2 2 2 2 1 1 1 1 1 1 ...
## $ npost: num 3 11 15 6 1 1 4 16 2 1 ...summary(postura)## trat linh npost
## Observ :18 Pouco reativo:21 Min. : 0.000
## Observ + Interv:20 Muito reativo:17 1st Qu.: 3.000
## Median : 6.000
## Mean : 7.895
## 3rd Qu.:11.750
## Max. :35.000bwplot(npost ~ linh | trat,
data = postura,
main = "Mudanças de postura vs tratamento e linhagem",
xlab = "Linhagem",
ylab = "Frequência",
pch = "|")
xyplot(npost ~ linh | trat,
data = postura,
main = "Mudanças de postura vs tratamento e linhagem",
xlab = "Linhagem",
ylab = "Frequência",
panel = panel.beeswarm, spread = 0.05)
# Variância e média amostrais por trat e linhagem.
mdp <- aggregate(npost ~ trat + linh,
data = postura,
FUN = function(x) {
c(mean = mean(x), var = var(x))
})
mdp## trat linh npost.mean npost.var
## 1 Observ Pouco reativo 9.600000 28.266667
## 2 Observ + Interv Pouco reativo 4.636364 20.454545
## 3 Observ Muito reativo 12.125000 101.267857
## 4 Observ + Interv Muito reativo 6.222222 38.444444Regressão Poisson com estimação por máxima verossimilhança
# Ajuste do modelo Poisson.
mP <- glm(npost ~ trat + linh,
data = postura,
family = poisson)
summary(mP)##
## Call:
## glm(formula = npost ~ trat + linh, family = poisson, data = postura)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.3579 -1.5592 -0.4320 0.7746 5.2884
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.25083 0.09286 24.238 < 2e-16 ***
## tratObserv + Interv -0.69665 0.12053 -5.780 7.48e-09 ***
## linhMuito reativo 0.25513 0.11549 2.209 0.0272 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 198.24 on 37 degrees of freedom
## Residual deviance: 158.46 on 35 degrees of freedom
## AIC: 298.65
##
## Number of Fisher Scoring iterations: 5# Avaliando possível efeito de interação.
mPi <- glm(npost ~ trat * linh,
data = postura,
family = poisson)
anova(mP, mPi, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: npost ~ trat + linh
## Model 2: npost ~ trat * linh
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 35 158.46
## 2 34 158.40 1 0.063318 0.8013exp(coef(mP)[2])## tratObserv + Interv
## 0.4982513# Estima-se que a taxa de variação de postura no grupo com intervenção
# seja aproximadamente a metade em relação ao grupo sem intervenção.
# Estimação do parâmetro de dispersão.
X2 <- sum(resid(mP, type = "pearson")^2)
phichap <- X2/df.residual(mP)
phichap## [1] 5.094182Diagnóstico do ajuste
# Diagnóstico do modelo - gráficos padrão do R.
par(mfrow = c(2, 2))
plot(mP)
envelope <- function(modelo) {
dados <- na.omit(modelo$data)
nsim <- 100
n <- modelo$df.null + 1
r1 <- sort(rstandard(modelo, type = "deviance"))
m1 <- matrix(0, nrow = n, ncol = nsim)
a2 <- simulate(modelo, nsim = nsim)
for (i in 1:nsim) {
dados$y <- a2[, i]
aj <- update(modelo, y ~ ., data = dados)
m1[, i] <- sort(rstandard(aj, type = "deviance"))
}
li <- apply(m1, 1, quantile, 0.025)
m <- apply(m1, 1, quantile, 0.5)
ls <- apply(m1, 1, quantile, 0.975)
quantis <- qnorm((1:n - 0.5)/n)
plot(rep(quantis, 2), c(li, ls), type = "n",
xlab = "Percentil da N(0,1)",
ylab = "Resíduos")
title("Gráfico Normal de Probabilidades")
lines(quantis, li, type = "l")
lines(quantis, m, type = "l", lty = 2)
lines(quantis, ls, type = "l")
points(quantis, r1, pch = 16, cex = 0.75)
}# Gráfico quantil-quantil com envelopes simulados.
layout(1)
envelope(mP)
Ajustando modelos por quase-verossimilhança
# Modelo quasi poisson (V(mu) = mu).
mQP0 <- glm(npost ~ trat + linh,
data = postura,
family = "quasipoisson")
# summary(mQP0)
# Forma alternativa de declarar o Modelo quase-poisson (V(mu) = mu).
mQP0 <- glm(npost ~ trat + linh,
data = postura,
family = quasi(variance = "mu", link = "log"))
# summary(mQP0)
# Modelo de quase-verossimilhança (V(mu) = mu^2).
mQP1 <- glm(npost ~ trat + linh,
data = postura,
family = quasi(variance = "mu^2", link = "log"))
summary(mQP1)##
## Call:
## glm(formula = npost ~ trat + linh, family = quasi(variance = "mu^2",
## link = "log"), data = postura)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5797 -0.6427 -0.1600 0.2356 1.6205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.2477 0.2337 9.616 2.34e-11 ***
## tratObserv + Interv -0.7007 0.2744 -2.554 0.0152 *
## linhMuito reativo 0.2655 0.2756 0.964 0.3419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasi family taken to be 0.7133735)
##
## Null deviance: 27.88 on 37 degrees of freedom
## Residual deviance: 22.67 on 35 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Gráficos de diagnóstico para o modelo de quase-verossimilhança.
par(mfrow = c(2,2))
plot(mQP1)
# Gráficos quantil-quantil para os resíduos dos modelos Poisson e de
# Quase-Verossimilhança.
par(mfrow = c(1, 2))
qqnorm(resid(mP, type = "deviance"),
pch = 20, main = "Poisson", las = 1)
qqline(resid(mP, type = "deviance"))
qqnorm(resid(mQP1, type = "deviance"),
pch = 20, main = "QL", las = 1)
qqline(resid(mQP1, type = "deviance"))
Usando estimação robusta e bootstrap
# Estimador sanduíche.
estrb <- coeftest(mP, vcov = sandwich)
estrb##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.25083 0.16483 13.6557 < 2.2e-16 ***
## tratObserv + Interv -0.69665 0.26594 -2.6196 0.008804 **
## linhMuito reativo 0.25513 0.25452 1.0024 0.316144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Usando bootstrap (R = 999 simulações)
mB <- Boot(mP)
# Resultados obtidos via bootstrap.
summary(mB)## R original bootBias bootSE bootMed
## (Intercept) 999 2.25083 -0.016707 0.18249 2.24898
## tratObserv + Interv 999 -0.69665 -0.006699 0.27398 -0.70538
## linhMuito reativo 999 0.25513 -0.016729 0.26200 0.24453Resumo geral dos resultados
erroz <- rbind(summary(mP)$coefficients[2,2:3],
summary(mQP0)$coefficients[2,2:3],
summary(mQP1)$coefficients[2,2:3],
estrb[2,2:3],
c(summary(mB)[2,4],
mean(mB$t[,2]/summary(mB)[2,4])))
ics <- rbind(confint.default(mP)[2,],
confint.default(mQP0)[2,],
confint.default(mQP1)[2,],
estrb[2,1] + c(-1.96,1.96) * estrb[2,2],
mean(mB$t[, 2]) +
c(-1.96, 1.96) * summary(mB)[2, 4])
quadres <- cbind(erroz, ics)
rownames(quadres) <- c("Poisson", "Quasi(mu)", "Quasi(mu^2)",
"Robusto (sanduiche)", "Bootstrap")
# Quadro resumo para as estimativas produzidas pelos quatro modelos para
# o efeito de intervenção.
kable(quadres,
format = "markdown",
caption = "Comparativo dos modelos ajustados")| Std. Error | z value | 2.5 % | 97.5 % | |
|---|---|---|---|---|
| Poisson | 0.1205303 | -5.779881 | -0.9328858 | -0.4604157 |
| Quasi(mu) | 0.2720405 | -2.560835 | -1.2298404 | -0.1634612 |
| Quasi(mu^2) | 0.2744137 | -2.553533 | -1.2385656 | -0.1628836 |
| Robusto (sanduiche) | 0.2659421 | -2.619558 | -1.2178973 | -0.1754043 |
| Bootstrap | 0.2739846 | -2.567114 | -1.2403596 | -0.1663399 |
Verificando efeito das observações com maiores resíduos na análise
postura.ex <- postura[-c(8, 18, 28), ]
# Ajustando o modelo Poisson sem as três observações.
mPe <- glm(npost ~ trat + linh,
data = postura.ex,
family = poisson)
# Estimativa do parâmetro de dispersão.
sum(resid(mPe, type = "pearson")^2)/mPe$df.residual## [1] 1.990952# Gráficos de diagnóstico.
par(mfrow = c(2, 2))
plot(mPe)
# Agora, o modelo quase-poisson.
mQPe <- glm(npost ~ trat + linh,
data = postura.ex,
family = quasi(variance = "mu", link = "log"))
# Estimativas produzidas pelo modelo quasipoisson com e sem as três
# observações.
c1 <- compareCoefs(mP,
mQP0,
mPe,
mQPe,
print = FALSE)
colnames(c1) <- c("Est. Poisson c/out", "E.P. Poisson c/out",
"Est. Quasi c/out", "E.P. Quasi c/out",
"Est. Poisson s/out", "E.P. Poisson s/out",
"Est. Quasi s/out", "E.P. Quasi s/out")
kable(round(c1, 3), align = "c")| Est. Poisson c/out | E.P. Poisson c/out | Est. Quasi c/out | E.P. Quasi c/out | Est. Poisson s/out | E.P. Poisson s/out | Est. Quasi s/out | E.P. Quasi s/out | |
|---|---|---|---|---|---|---|---|---|
| (Intercept) | 2.251 | 0.093 | 2.251 | 0.210 | 2.227 | 0.097 | 2.227 | 0.137 |
| tratObserv + Interv | -0.697 | 0.121 | -0.697 | 0.272 | -0.886 | 0.144 | -0.886 | 0.204 |
| linhMuito reativo | 0.255 | 0.115 | 0.255 | 0.261 | 0.005 | 0.134 | 0.005 | 0.190 |
# Efeito da intervenção desconsiderando as três observações.
exp(coef(mPe)[2])## tratObserv + Interv
## 0.4123804# O efeito de intervenção aumenta, e torna-se mais significativo,
# mediante exclusão dos outliers.Informações da sessão
## Atualizado em 21 de julho de 2016.
##
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04 LTS
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets
## [7] methods base
##
## other attached packages:
## [1] hnp_1.2 sandwich_2.3-4 car_2.1-2
## [4] boot_1.3-18 lmtest_0.9-34 zoo_1.7-13
## [7] multcomp_1.4-6 TH.data_1.0-7 MASS_7.3-45
## [10] survival_2.39-5 mvtnorm_1.0-5 doBy_4.5-15
## [13] rmarkdown_1.0 latticeExtra_0.6-28 RColorBrewer_1.1-2
## [16] lattice_0.20-33 knitr_1.13 MRDCr_0.0-2
## [19] gsubfn_0.6-6 proto_0.3-10 bbmle_1.0.18
## [22] devtools_1.12.0
##
## loaded via a namespace (and not attached):
## [1] splines_3.3.1 tcltk_3.3.1 testthat_1.0.2
## [4] htmltools_0.3.5 mgcv_1.8-12 yaml_2.1.13
## [7] nloptr_1.0.4 withr_1.0.2 stringr_1.0.0
## [10] MatrixModels_0.4-1 codetools_0.2-14 memoise_1.0.0
## [13] evaluate_0.9 SparseM_1.7 quantreg_5.26
## [16] pbkrtest_0.4-6 parallel_3.3.1 highr_0.6
## [19] Rcpp_0.12.5 formatR_1.4 lme4_1.1-12
## [22] digest_0.6.9 stringi_1.1.1 numDeriv_2014.2-1
## [25] grid_3.3.1 tools_3.3.1 magrittr_1.5
## [28] crayon_1.3.2 Matrix_1.2-6 minqa_1.2.4
## [31] roxygen2_5.0.1 R6_2.1.2 nnet_7.3-12
## [34] nlme_3.1-128 compiler_3.3.1
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 |