diff --git a/vignettes/v02_quase_poisson.Rmd b/vignettes/v02_quase_poisson.Rmd index b5303b548c2872d63add0d3617ed9891de3e633e..eb8a1256445da371e29a376fcbfcdb922f0ad155 100644 --- a/vignettes/v02_quase_poisson.Rmd +++ b/vignettes/v02_quase_poisson.Rmd @@ -1,11 +1,11 @@ --- -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() +```