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

Merge branch 'master' into devel

Conflicts:
	.gitignore
parents ae210098 78010f14
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
---
title: "Análise de Contagens - Quase-Verossimilhança"
author: >
Walmes M. Zeviani,
Eduardo E. Ribeiro Jr &
Cesar A. Taconeli
vignette: >
%\VignetteIndexEntry{Análise de Contagens - Quase-Verossimilhança}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include=FALSE}
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.
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('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')
```
### Verificação do conteúdo e a estrutura dos dados.
```{r}
str(datapost4)
summary(datapost4)
```
### 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)))
mdp
```
### 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
phichap
```
### Diagnóstico do ajuste (gráficos).
```{r}
##### Diagnóstico do modelo - gráficos padrão do R.
par(mfrow=c(2,2))
plot(ajusteps, pch = 20, cex = 1.25)
```
```{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}
##### Gráfico quantil-quantil com envelopes simulados.
par(mfrow=c(1,1))
envelope(ajusteps)
```
### 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.
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'))
```
### Usando estimação robusta e bootstrap.
```{r}
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.
```
### 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])
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")
```
### 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.
```
---
title: "Análise de Contagens - Poisson e Binomial Negativa"
author: >
Walmes M. Zeviani,
Eduardo E. Ribeiro Jr &
Cesar A. Taconeli
vignette: >
%\VignetteIndexEntry{Análise de Contagens - Poisson e Binomial Negativa}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include=FALSE}
source("_setup.R")
```
Dados referentes ao número de sinistros registrados por 16483 clientes de
uma seguradora de automóveis ao longo de um ano, contemplando as
seguintes variáveis:
* **Sinistros**: Número de sinistros registrados;
* **Exposicao**: Período de cobertura do cliente durante o ano sob análise;
* **Idade**: Idade do cliente (em anos);
* **Sexo**: M para masculino e F para feminino;
* **Valor**: Valor do veículo segurado (em reais).
```{r, echo = FALSE, include=FALSE}
require(lattice)
require(MASS)
require(effects)
require(knitr)
dados <- read.csv('Baseauto2.csv', header=TRUE, sep=',', dec = ',')
options(scipen = 3) ### Era zero.
##### Preparando os dados.
names(dados) <- c('Bonus','Tipo','Estado','Idade','Sexo','Civil','Valor','Exposicao','Sinistros')
dados <- dados[,c('Idade','Sexo','Valor','Exposicao','Sinistros')]
dados <- dados[which(dados$Valor != 0),]
dados$lexpo <- log(dados$Exposicao)
dados$Valor <- dados$Valor/1000
```
### Verificação do conteúdo e a estrutura dos dados.
```{r}
head(dados, 10) ### Dez primeiras linhas da base.
str(dados)
```
### Análise descritiva da distribuição do número de sinistros.
```{r}
table(dados$Sinistros) ### Distribuição do números de sinistros.
taxageral <- sum(dados$Sinistros)/sum(dados$Exposicao); taxageral ### Taxa de sinistros na amostra.
tab <- aggregate(cbind(Exposicao, Sinistros) ~ Sexo, data = dados, sum)
taxa <- with(tab, Sinistros/Exposicao)
tab <- cbind(tab, taxa)
### Distribuição do número de sinistros por sexo.
kable(tab, align = 'c', caption = '**Taxa de sinistros segundo Sexo**')
dados$idadecat <- cut(dados$Idade, breaks=c(18,30,60, 92), include.lowest = T)
tab <- aggregate(cbind(Exposicao, Sinistros) ~ idadecat, data = dados, sum)
taxa <- with(tab, Sinistros/Exposicao)
tab <- cbind(tab, taxa)
### Distribuição do número de sinistros por sexo.
kable(tab, align = 'c', caption = '**Taxa de sinistros segundo Idade**')
tabidsex <- aggregate(cbind(Exposicao, Sinistros) ~ Sexo + idadecat, data = dados, sum)
taxa <- with(tabidsex, Sinistros/Exposicao)
tabidsex <- cbind(tabidsex, taxa)
### Distribuição do número de sinistros por idade e sexo.
kable(tabidsex, align = 'c', caption = '**Taxa de sinistros segundo Sexo e Idade**')
```
## Regressão usando o modelo log-linear poisson.
```{r}
dados <- na.omit(dados)
ajusteps <- glm(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+offset(log(Exposicao)), data = dados, family=poisson)
summary(ajusteps)
### Estimação do parâmetro de dispersão.
X2 <- sum(resid(ajusteps,type='pearson')**2)
phichap <- X2/ajusteps$df.residual
phichap ### Indicador de superdispersão.
```
```{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)
}
```
### Diagnóstico do ajuste (gráficos).
```{r}
##### Diagnóstico do modelo - gráficos.
par(mfrow=c(2,2))
plot(ajusteps)
par(mfrow=c(1,1))
envelope(ajusteps)
```
### Ajuste do modelo associando um parâmetro ao termo offset (log-exposicao)
```{r}
ajusteps2 <- glm(Sinistros ~ Sexo + Idade +I(Idade**2) + Valor + log(Exposicao), data = dados, family=poisson)
summary(ajusteps2)
anova(ajusteps, ajusteps2, test = 'Chisq')
```
## Regressão usando a distribuição binomial negativa.
```{r}
ajustenb2 <- glm.nb(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+log(Exposicao),data= dados)
summary(ajustenb2)
```
### Diagnóstico do ajuste (gráficos).
```{r}
##### Diagnóstico do modelo - gráficos.
par(mfrow=c(2,2))
plot(ajustenb2)
```
```{r, echo = FALSE}
dadosnb3 <- dados[,c('Sexo','Idade','Valor','Exposicao','Sinistros')]
dadosnb3$lexpo <- log(dados$Exposicao)
ajustenb2 <- glm.nb(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+lexpo,data= dadosnb3)
envelope=function(modelo){
dados=na.omit(dadosnb3)
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}
par(mfrow=c(1,1))
envelope(ajustenb2)
```
### Explorando os efeitos das covariáveis. Estimativas pontuais e ICs (95%)
```{r}
intervalos <- confint(ajustenb2)
estimat <- cbind(ajustenb2$coefficients, intervalos)
colnames(estimat)[1] <- 'Estimativa pontual'
### Quadro de estimativas
kable(round(estimat, 5), align = 'c', caption = '**Estimativas pontuais e intervalos de confiança - Modelo Binomial Negativo**')
```
### Gráficos de efeitos
```{r, echo = FALSE}
```
```{r}
efeitos <- allEffects(ajustenb2, given.values=c(lexpo=0))
trellis.par.set(list(axis.text = list(cex = 1.2)))
plot(efeitos[[2]], type='response',main=list(
label="Taxa de sinistros vs. Idade",
cex=1.5),
xlab=list(
label="Idade (anos)",
cex=1.5),
ylab=list(
label="Taxa de sinistros",
cex=1.5))
plot(efeitos[[1]], type='response',main=list(
label="Taxa de sinistros vs. Sexo",
cex=1.5),
xlab=list(
label="Sexo",
cex=1.5),
ylab=list(
label="Taxa de sinistros",
cex=1.5))
plot(efeitos[[4]], type='response',main=list(
label="Taxa de sinistros vs. Valor do automóvel",
cex=1.5),
xlab=list(
label="Valor (x1000 reais)",
cex=1.5),
ylab=list(
label="Taxa de sinistros",
cex=1.5))
```
```{r, echo = FALSE, include=F}
## Frequências ajustadas pelas duas distribuições, com e sem covariaveis.
### Poisson sem ajuste de covariáveis.
n <- nrow(dados)
mediasin <- mean(dados$Sinistros)
freqps <- round(n*dpois(0:10,mediasin))
### Poisson com covariaveis
pred1 <- predict(ajusteps2,type='response')
intervalo <- 0:10
matprob <- matrix(0, nrow=nrow(dados), ncol=length(intervalo))
probpois <- function(interv, taxa)
dpois(intervalo,taxa)
for(i in 1:nrow(dados))
matprob[i,] <- probpois(interv = intervalo, taxa = pred1[i])
pbarra <- colMeans(matprob)
freqpsaj <- round(n*pbarra)
### Binomial negativa sem covariaveis.
ajustenb <- glm.nb(Sinistros ~ 1,data=dados)
media <- exp(coefficients(ajustenb))
shape <- ajustenb$theta
freqbn <- round(n*dnbinom(0:10, mu = media, size = shape))
### Binomial negativa com covariaveis
pred2 <- predict(ajustenb2,type='response')
intervalo <- 0:10
matprob <- matrix(0,nrow=nrow(dados),ncol=length(intervalo))
probnb <- function(interv, media, shape)
dnbinom(intervalo, mu = media, size = shape)
for(i in 1:nrow(dados))
matprob[i,] <- probnb(interv = intervalo, media = pred2[i], shape = ajustenb2$theta)
pbarra <- colMeans(matprob)
frebnaj <- round(n*pbarra)
ams <- c(table(dados$Sinistros), rep(0,5))
matfreq <- rbind(ams, freqps, freqpsaj, freqbn, frebnaj)
colnames(matfreq) <- 0:10
rownames(matfreq) <- c('Amostra', 'Poisson não ajustada por covariáveis', 'Poisson ajustada por covariáveis',
'BN não ajustada por covariáveis', 'BN ajustada por covariáveis')
```
## Comparação dos ajustes
```{r, results = 'markup'}
kable(matfreq, format = "markdown", caption = "Frequências amostrais e frequências ajustadas para o número de sinistros")
```
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