From f8d8979b0f0979766327f8a37a931e1b06a1fbd4 Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Thu, 21 Jul 2016 15:20:28 -0300
Subject: [PATCH] Padroniza a vinheta com dados de seguro.
---
vignettes/v03_binomial_negativa.Rmd | 451 ++++++++++++++--------------
1 file changed, 230 insertions(+), 221 deletions(-)
diff --git a/vignettes/v03_binomial_negativa.Rmd b/vignettes/v03_binomial_negativa.Rmd
index 0316b85..f986500 100644
--- a/vignettes/v03_binomial_negativa.Rmd
+++ b/vignettes/v03_binomial_negativa.Rmd
@@ -1,11 +1,11 @@
---
-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()
```
--
GitLab