Skip to content
Snippets Groups Projects
Commit f8d8979b authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Padroniza a vinheta com dados de seguro.

parent 1b838b0a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
vignettes/v03_binomial_negativa.Rmd
+ 230
− 221
View file @ f8d8979b
---
title: "Análise de Contagens - Poisson e Binomial Negativa"
title: "Análise de Contagens com Distribuição Binomial Negativa"
author: >
Walmes M. Zeviani,
Eduardo E. Ribeiro Jr &
Cesar A. Taconeli
vignette: >
%\VignetteIndexEntry{Análise de Contagens - Poisson e Binomial Negativa}
%\VignetteIndexEntry{Análise de Contagens com Distribuição Binomial Negativa}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
......@@ -14,305 +14,314 @@ vignette: >
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
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('../data-raw/Baseauto2.csv', header=TRUE, sep=',', dec = ',')
* **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).
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
```{r, include=FALSE}
devtools::load_all()
```
```{r, results = "hide", message = FALSE}
# Pacotes necessários.
library(lattice)
library(MASS)
library(effects)
library(knitr)
```
### Verificação do conteúdo e a estrutura dos dados
## Verificação do conteúdo e a estrutura dos dados
```{r}
head(dados, 10) ### Dez primeiras linhas da base.
str(dados)
# Dez primeiras linhas da base.
head(seguro, 10)
str(seguro)
```
### Análise descritiva da distribuição do número de sinistros
## 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.
# Distribuição do números de sinistros.
tb <- table(seguro$Sinistros)
tb
barchart(tb, horizontal = FALSE)
tab <- aggregate(cbind(Exposicao, Sinistros) ~ Sexo, data = dados, sum)
# Taxa de sinistros na amostra.
taxageral <- sum(seguro$Sinistros)/sum(seguro$Exposicao)
taxageral
tab <- aggregate(cbind(Exposicao, Sinistros) ~ Sexo,
data = seguro, FUN = 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)
# Distribuição do número de sinistros por sexo.
kable(tab, align = "c",
caption = "**Taxa de sinistros segundo Sexo**")
seguro$idadecat <- cut(seguro$Idade,
breaks = c(18, 30, 60, 92),
include.lowest = TRUE)
tab <- aggregate(cbind(Exposicao, Sinistros) ~ idadecat,
data = seguro, FUN = 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)
# 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 = seguro, FUN = 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**')
# 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
## 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.
seguro <- na.omit(seguro)
mP <- glm(Sinistros ~ Sexo + Idade + I(Idade^2) + Valor +
offset(log(Exposicao)),
data = seguro, family = poisson)
summary(mP)
# Estimação do parâmetro de dispersão.
X2 <- sum(resid(mP, type = "pearson")^2)
phichap <- X2/mP$df.residual
# Indicador de superdispersão.
phichap
```
```{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)
}
```
### Diagnóstico do ajuste (gráficos)
## Diagnóstico do ajuste (gráficos)
```{r}
##### Diagnóstico do modelo - gráficos.
par(mfrow=c(2,2))
plot(ajusteps)
# Diagnóstico do modelo - gráficos.
par(mfrow = c(2, 2))
plot(mP)
par(mfrow=c(1,1))
envelope(ajusteps)
par(mfrow = c(1, 1))
envelope(mP)
```
### Ajuste do modelo associando um parâmetro ao termo offset (log-exposicao)
## 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')
mPo <- glm(Sinistros ~ Sexo + Idade + I(Idade^2) + Valor +
log(Exposicao),
data = seguro, family = poisson)
summary(mPo)
anova(mP, mPo, 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)
mBNo <- glm.nb(Sinistros ~ Sexo + Idade + I(Idade^2) + Valor +
log(Exposicao), data = seguro)
summary(mBNo)
```
### Diagnóstico do ajuste (gráficos)
## Diagnóstico do ajuste
```{r}
##### Diagnóstico do modelo - gráficos.
par(mfrow=c(2,2))
plot(ajustenb2)
# Diagnóstico do modelo - gráficos.
par(mfrow = c(2, 2))
plot(mBNo)
```
```{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}
dadosnb3 <-
seguro[, c("Sexo", "Idade", "Valor", "Exposicao", "Sinistros")]
dadosnb3$lexpo <- log(seguro$Exposicao)
mBNo <- 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)
par(mfrow = c(1, 1))
envelope(mBNo)
```
### Explorando os efeitos das covariáveis. Estimativas pontuais e ICs (95%)
## Explorando os efeitos das covariáveis
```{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**')
intervalos <- confint(mBNo)
estimat <- cbind(mBNo$coefficients, intervalos)
colnames(estimat)[1] <- "Estimativa pontual"
# Quadro de estimativas
kable(round(estimat, 5), align = "c",
caption = paste("**Estimativas pontuais e intervalos de",
"confiança - Modelo Binomial Negativo**"))
```
### Gráficos de efeitos
```{r, echo = FALSE}
```
## Gráficos de efeitos
```{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))
efeitos <- allEffects(mBNo, given.values = c(lexpo = 0))
trellis.par.set(list(axis.text = list(cex = 1.2)))
plot(efeitos[[2]],
type = "response",
main = "Taxa de sinistros vs. Idade",
xlab = "Idade (anos)",
ylab = "Taxa de sinistros")
plot(efeitos[[1]],
type = "response",
main = "Taxa de sinistros vs. Sexo",
xlab = "Sexo",
ylab = "Taxa de sinistros")
plot(efeitos[[4]],
type = "response",
main = "Taxa de sinistros vs. Valor do automóvel",
xlab = "Valor (x1000 reais)",
ylab = "Taxa de sinistros")
```
## Frequências ajustadas pelas duas distribuições
```{r}
# Poisson sem ajuste de covariáveis.
n <- nrow(seguro)
mediasin <- mean(seguro$Sinistros)
freqps <- round(n * dpois(0:10, mediasin))
```{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')
# Poisson com covariaveis
pred1 <- predict(mPo, 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])
matprob <- matrix(0, nrow = nrow(seguro), ncol = length(intervalo))
probpois <- function(interv, taxa) dpois(intervalo, taxa)
for (i in 1:nrow(seguro)) {
matprob[i, ] <- probpois(interv = intervalo, taxa = pred1[i])
}
pbarra <- colMeans(matprob)
freqpsaj <- round(n*pbarra)
freqpsaj <- round(n * pbarra)
### Binomial negativa sem covariaveis.
ajustenb <- glm.nb(Sinistros ~ 1,data=dados)
# Binomial negativa sem covariaveis.
ajustenb <- glm.nb(Sinistros ~ 1, data = seguro)
media <- exp(coefficients(ajustenb))
shape <- ajustenb$theta
freqbn <- round(n*dnbinom(0:10, mu = media, size = shape))
freqbn <- round(n * dnbinom(0:10, mu = media, size = shape))
### Binomial negativa com covariaveis
pred2 <- predict(ajustenb2,type='response')
# Binomial negativa com covariaveis
pred2 <- predict(mBNo, 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)
matprob <- matrix(0, nrow = nrow(seguro), ncol = length(intervalo))
probnb <- function(interv, media, shape) {
dnbinom(intervalo, mu = media,
size = shape)
}
for (i in 1:nrow(seguro)) {
matprob[i, ] <- probnb(interv = intervalo, media = pred2[i],
shape = mBNo$theta)
}
pbarra <- colMeans(matprob)
frebnaj <- round(n*pbarra)
ams <- c(table(dados$Sinistros), rep(0,5))
frebnaj <- round(n * pbarra)
ams <- c(table(seguro$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')
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")
```{r, results="markup"}
kable(matfreq, format = "markdown",
caption = paste("Frequências amostrais e frequências",
"ajustadas para o número de sinistros"))
```
## 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()
```
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