Padronização da vinheta quasi-poisson.

parent b71bd572
---
title: "Análise de Contagens - Quase-Verossimilhança"
title: "Análise de Contagens por Quase-Verossimilhança"
author: >
Walmes M. Zeviani,
Eduardo E. Ribeiro Jr &
Cesar A. Taconeli
vignette: >
%\VignetteIndexEntry{Análise de Contagens - Quase-Verossimilhança}
%\VignetteIndexEntry{Análise de Contagens por Quase-Verossimilhança}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
......@@ -14,252 +14,262 @@ vignette: >
source("_setup.R")
```
Dados referentes a um experimento realizado com o objetivo de
investigar o efeito de uma intervenção, por parte do cuidador, no
comportamento de ovelhas.
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).
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).
```{r, echo = FALSE, include=FALSE}
##### Carregamento e tratamento inicial dos dados
dados <- read.csv('../data-raw/Dadoscomp.csv',sep=',')
dados$tratamento <- factor(dados$tratamento)
levels(dados$tratamento) <- c('Observ', 'Observ + Interv')
dados$linhagem <- factor(dados$linhagem)
levels(dados$linhagem) <- c('Pouco reativo', 'Muito reativo')
dados2 <- dados[1:38,c(1:5,30:53)]
dados3 <- dados[1:38,30:53] ### Somente as mudanças de postura durante a intervenção.
r1 <- rowSums(dados3[,1:3])-1
table(r1)
r2 <- rowSums(dados3[,4:6])-1
table(r2)
r4 <- rowSums(dados3[,12:16])-1
table(r4)
d2 <- rowSums(dados[1:38,57:59])-1
table(d2)
datapost1 <- data.frame(r1, dados2[ ,c('tratamento', 'linhagem')])
names(datapost1)[1] <- 'Nposturas'
datapost2 <- data.frame(r2, dados2[ ,c('tratamento', 'linhagem')])
names(datapost2)[1] <- 'Nposturas'
datapost4 <- data.frame(r4, dados2[ ,c('tratamento', 'linhagem')])
names(datapost4)[1] <- 'Nposturas'
##### Pacotes requeridos
require('lmtest')
require('boot')
require('car')
require('RColorBrewer')
require('sandwich')
require('hnp')
require('knitr')
```{r, echo=FALSE, include=FALSE}
devtools::load_all()
```
### Verificação do conteúdo e a estrutura dos dados.
```{r}
str(datapost4)
summary(datapost4)
```{r, results="hide", message=FALSE}
# Pacotes requeridos.
library("lmtest")
library("boot")
library("car")
library("latticeExtra")
library("RColorBrewer")
library("sandwich")
library("hnp")
library("knitr")
```
### Análise descritiva
## Análise descritiva
```{r}
tab <- data.frame(with(datapost4, table(tratamento, Nposturas)))
myColours <- brewer.pal(2,"Blues")
my.settings <- list(
superpose.polygon=list(col=myColours[2:5], border="transparent"),
strip.background=list(col=myColours[6]),
strip.border=list(col="black")
)
bwplot(Nposturas ~ linhagem | tratamento,
data=datapost4,
main="Mudanças de postura vs tratamento e linhagem",
xlab="Linhagem", ylab="Frequência")
### Variância e média amostrais por tratamento e linhagem.
mdp <- aggregate(Nposturas ~ tratamento:linhagem, datapost4, function(x) c(mean = mean(x), var = var(x)))
str(postura)
summary(postura)
# names(postura) <- c("npost", "trat", "linh")
# postura <- postura[, c(2, 3, 1)]
# use_data(postura, overwrite = TRUE)
# tab <- xtabs(~trat + npost, data = postura)
bwplot(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
```
### Regressão poisson com estimação por máxima verossimilhança.
## Regressão poisson com estimação por máxima verossimilhança
```{r}
ajusteps <- glm(Nposturas ~ tratamento + linhagem, data=datapost4, family=poisson)
summary(ajusteps)
### Avaliando possível efeito de interação.
ajustepsint <- glm(Nposturas ~ tratamento * linhagem, data=datapost4, family=poisson)
anova(ajusteps,ajustepsint, test = 'Chisq')
exp(coef(ajusteps)[2])
### 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(ajusteps,type='pearson')**2)
phichap <- X2/ajusteps$df.residual
# Ajuste do modelo Poisson.
mP <- glm(npost ~ trat + linh,
data = postura,
family = poisson)
summary(mP)
# Avaliando possível efeito de interação.
mPi <- glm(npost ~ trat * linh,
data = postura,
family = poisson)
anova(mP, mPi, test = "Chisq")
exp(coef(mP)[2])
# 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
```
### Diagnóstico do ajuste (gráficos).
## Diagnóstico do ajuste
```{r}
##### Diagnóstico do modelo - gráficos padrão do R.
par(mfrow=c(2,2))
plot(ajusteps, pch = 20, cex = 1.25)
# Diagnóstico do modelo - gráficos padrão do R.
par(mfrow = c(2, 2))
plot(mP)
```
```{r, echo = FALSE}
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)
```{r}
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)
}
```
```{r}
##### Gráfico quantil-quantil com envelopes simulados.
par(mfrow=c(1,1))
envelope(ajusteps)
# Gráfico quantil-quantil com envelopes simulados.
layout(1)
envelope(mP)
```
### Ajustando modelos por quase-verossimilhança.
## Ajustando modelos por quase-verossimilhança
```{r}
### Modelo quasi poisson (V(mu) = mu).
ajuste12 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = 'quasipoisson')
# summary(ajuste12)
### Forma alternativa de declarar o Modelo quase-poisson (V(mu) = mu).
ajuste13 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = quasi(variance = 'mu', link='log'))
# summary(ajuste13)
### Modelo de quase-verossimilhança (V(mu) = mu^2).
ajuste14 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = quasi(variance = 'mu^2', link='log'))
summary(ajuste14)
### Gráficos de diagnóstico para o modelo de 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)
# Gráficos de diagnóstico para o modelo de quase-verossimilhança.
par(mfrow = c(2,2))
plot(ajuste14, pch = 20, cex = 1.25)
### Gráficos quantil-quantil para os resíduos dos modelos Poisson e de Quase-Verossimilhança.
par(mfrow=c(1,2))
qqnorm(resid(ajusteps,type='deviance'), pch = 20, main = 'Poisson', las = 1)
qqline(resid(ajusteps,type='deviance'))
qqnorm(resid(ajuste14,type='deviance'), pch = 20, main = 'QL', las = 1)
qqline(resid(ajuste14,type='deviance'))
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.
## Usando estimação robusta e bootstrap
```{r}
# Estimador sanduíche.
estrb <- coeftest(mP, vcov = sandwich)
estrb
estrb <- coeftest(ajusteps, vcov=sandwich) ### Estimador sanduíche
estrb ### Estimação robusta.
### Usando bootstrap (R=999 simulações)
ajusteboot <- Boot(ajusteps)
summary(ajusteboot) ### Resultados obtidos via bootstrap.
# Usando bootstrap (R = 999 simulações)
mB <- Boot(mP)
# Resultados obtidos via bootstrap.
summary(mB)
```
### Resumo geral dos resultados.
## Resumo geral dos resultados
```{r}
erroz <- rbind(summary(ajusteps)$coefficients[2,2:3], summary(ajuste13)$coefficients[2,2:3],
summary(ajuste14)$coefficients[2,2:3], estrb[2,2:3], c(summary(ajusteboot)[2,4],
mean(ajusteboot$t[,2]/summary(ajusteboot)[2,4])))
ics <- rbind(confint.default(ajusteps)[2,],confint.default(ajuste13)[2,], confint.default(ajuste14)[2,],
estrb[2,1] + c(-1.96,1.96) * estrb[2,2], mean(ajusteboot$t[,2])+c(-1.96,1.96)*summary(ajusteboot)[2,4])
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")
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")
```
### Verificando efeito das observações com maiores resíduos na análise.
## Verificando efeito das observações com maiores resíduos na análise
```{r}
dadosexclud <- datapost4[-c(8,18,28),]
### Ajustando o modelo Poisson sem as três observações.
ajusteexcludpoi <- glm(Nposturas ~ tratamento+linhagem, data=dadosexclud, family = poisson)
### Estimativa do parâmetro de dispersão.
sum(resid(ajusteexcludpoi, type = 'pearson')**2)/ajusteexcludpoi$df.residual
### Gráficos de diagnóstico.
par(mfrow=c(2,2))
plot(ajusteexcludpoi)
### Agora, o modelo quase-poisson
ajusteexcludquasi <- glm(Nposturas ~ tratamento+linhagem, data=dadosexclud, family = quasi(variance = 'mu', link='log'))
### Estimativas produzidas pelo modelo quasipoisson com e sem as três observações.
c1 <- compareCoefs(ajusteps, ajuste13, ajusteexcludpoi, ajusteexcludquasi, 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')
### Efeito da intervenção desconsiderando as três observações.
exp(coef(ajusteexcludpoi)[2])
### O efeito de intervenção aumenta, e torna-se mais significativo, mediante exclusão dos outliers.
```
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
# 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")
# Efeito da intervenção desconsiderando as três observações.
exp(coef(mPe)[2])
# O efeito de intervenção aumenta, e torna-se mais significativo,
# mediante exclusão dos outliers.
```
## Informações da sessão
```{r, echo=FALSE, results="hold"}
cat(format(Sys.time(),
format = "Atualizado em %d de %B de %Y.\n\n"))
sessionInfo()
```
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