v03_binomial_negativa.Rmd 8.46 KB
Newer Older
Cesar Taconeli's avatar
Cesar Taconeli committed
1
---
2
title: "Análise de Contagens com Distribuição Binomial Negativa"
Cesar Taconeli's avatar
Cesar Taconeli committed
3 4 5 6 7
author: >
  Walmes M. Zeviani,
  Eduardo E. Ribeiro Jr &
  Cesar A. Taconeli
vignette: >
8
  %\VignetteIndexEntry{Análise de Contagens com Distribuição Binomial Negativa}
Cesar Taconeli's avatar
Cesar Taconeli committed
9 10 11 12 13 14 15 16
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
source("_setup.R")
```

17 18
Dados referentes ao número de sinistros registrados por 16483 clientes
de uma seguradora de automóveis ao longo de um ano, contemplando as
Cesar Taconeli's avatar
Cesar Taconeli committed
19 20
seguintes variáveis:

21 22 23 24 25 26
  * **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).
Cesar Taconeli's avatar
Cesar Taconeli committed
27

28 29 30 31 32 33
```{r, results = "hide", message = FALSE}
# Pacotes necessários.
library(lattice)
library(MASS)
library(effects)
library(knitr)
34
library(MRDCr)
Cesar Taconeli's avatar
Cesar Taconeli committed
35 36
```

37
## Verificação do conteúdo e a estrutura dos dados
Cesar Taconeli's avatar
Cesar Taconeli committed
38 39

```{r}
40
# Dez primeiras linhas da base.
41 42 43 44
head(seguros, 10)
str(seguros)

seguros$lexpo <- log(seguros$expos)
Cesar Taconeli's avatar
Cesar Taconeli committed
45 46
```

47
## Análise descritiva da distribuição do número de sinistros
Cesar Taconeli's avatar
Cesar Taconeli committed
48 49

```{r}
50
# Distribuição do números de sinistros.
51
tb <- table(seguros$nsinist)
52 53 54
tb

barchart(tb, horizontal = FALSE)
Cesar Taconeli's avatar
Cesar Taconeli committed
55

56
# Taxa de sinistros na amostra.
57
taxageral <- sum(seguros$nsinist)/sum(seguros$expos)
58 59
taxageral

60 61 62
tab <- aggregate(cbind(expos, nsinist) ~ sexo,
                 data = seguros, FUN = sum)
taxa <- with(tab, nsinist/expos)
63
tab <- cbind(tab, taxa)
Cesar Taconeli's avatar
Cesar Taconeli committed
64

65 66
# Distribuição do número de sinistros por sexo.
kable(tab, align = "c",
67
      caption = "**Taxa de sinistros segundo sexo**")
68

69
seguros$idadecat <- cut(seguros$idade,
70 71
                       breaks = c(18, 30, 60, 92),
                       include.lowest = TRUE)
72 73 74
tab <- aggregate(cbind(expos, nsinist) ~ idadecat,
                 data = seguros, FUN = sum)
taxa <- with(tab, nsinist/expos)
75
tab <- cbind(tab, taxa)
Cesar Taconeli's avatar
Cesar Taconeli committed
76

77 78
# Distribuição do número de sinistros por sexo.
kable(tab, align = "c",
79
      caption = "**Taxa de sinistros segundo idade**")
80

81 82 83
tabidsex <- aggregate(cbind(expos, nsinist) ~ sexo + idadecat,
                      data = seguros, FUN = sum)
taxa <- with(tabidsex, nsinist/expos)
84
tabidsex <- cbind(tabidsex, taxa)
Cesar Taconeli's avatar
Cesar Taconeli committed
85

86 87
# Distribuição do número de sinistros por idade e sexo.
kable(tabidsex, align = "c",
88
      caption = "**Taxa de sinistros segundo sexo e idade**")
Cesar Taconeli's avatar
Cesar Taconeli committed
89 90
```

91
## Regressão usando o modelo log-linear Poisson
Cesar Taconeli's avatar
Cesar Taconeli committed
92 93

```{r}
94 95 96 97
seguros <- na.omit(seguros)
mP <- glm(nsinist ~ sexo + idade + I(idade^2) + valor +
              offset(log(expos)),
          data = seguros, family = poisson)
98 99 100 101 102 103 104 105
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
Cesar Taconeli's avatar
Cesar Taconeli committed
106 107
```

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
```{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)
Cesar Taconeli's avatar
Cesar Taconeli committed
137
}
138
```
Cesar Taconeli's avatar
Cesar Taconeli committed
139

140
## Diagnóstico do ajuste (gráficos)
141 142

```{r}
143 144 145
# Diagnóstico do modelo - gráficos.
par(mfrow = c(2, 2))
plot(mP)
146

147 148
par(mfrow = c(1, 1))
envelope(mP)
Cesar Taconeli's avatar
Cesar Taconeli committed
149 150
```

151
## Ajuste do modelo associando um parâmetro ao termo offset (log-exposicao)
Cesar Taconeli's avatar
Cesar Taconeli committed
152 153

```{r}
154 155 156
mPo <- glm(nsinist ~ sexo + idade + I(idade^2) + valor +
               log(expos),
           data = seguros, family = poisson)
157 158
summary(mPo)
anova(mP, mPo, test = "Chisq")
Cesar Taconeli's avatar
Cesar Taconeli committed
159 160
```

161
## Regressão usando a distribuição binomial negativa
Cesar Taconeli's avatar
Cesar Taconeli committed
162 163

```{r}
164 165
mBNo <- glm.nb(nsinist ~ sexo + idade + I(idade^2) + valor +
                   log(expos), data = seguros)
166
summary(mBNo)
Cesar Taconeli's avatar
Cesar Taconeli committed
167 168
```

169
## Diagnóstico do ajuste
Cesar Taconeli's avatar
Cesar Taconeli committed
170 171

```{r}
172 173 174
# Diagnóstico do modelo - gráficos.
par(mfrow = c(2, 2))
plot(mBNo)
Cesar Taconeli's avatar
Cesar Taconeli committed
175 176
```

177 178
```{r}
dadosnb3 <-
179 180
    seguros[, c("sexo", "idade", "valor", "expos", "nsinist")]
dadosnb3$lexpo <- log(seguros$expos)
181

182 183
mBNo <- glm.nb(nsinist ~ sexo + idade + I(idade^2) +
                   valor + lexpo,
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
               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)
Cesar Taconeli's avatar
Cesar Taconeli committed
214 215 216 217
}
```

```{r}
218 219
par(mfrow = c(1, 1))
envelope(mBNo)
Cesar Taconeli's avatar
Cesar Taconeli committed
220 221
```

222
## Explorando os efeitos das covariáveis
Cesar Taconeli's avatar
Cesar Taconeli committed
223 224

```{r}
225 226 227 228 229 230 231 232
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**"))
Cesar Taconeli's avatar
Cesar Taconeli committed
233 234
```

235
## Gráficos de efeitos
Cesar Taconeli's avatar
Cesar Taconeli committed
236

237
```{r}
238 239 240 241
efeitos <- allEffects(mBNo, given.values = c(lexpo = 0))

plot(efeitos[[2]],
     type = "response",
242
     main = "Taxa de sinistros vs. idade",
243 244 245 246 247
     xlab = "Idade (anos)",
     ylab = "Taxa de sinistros")

plot(efeitos[[1]],
     type = "response",
248
     main = "Taxa de sinistros vs. sexo",
249 250 251 252 253
     xlab = "Sexo",
     ylab = "Taxa de  sinistros")

plot(efeitos[[4]],
     type = "response",
254
     main = "Taxa de sinistros vs. valor do automóvel",
255 256
     xlab = "Valor (x1000 reais)",
     ylab = "Taxa de sinistros")
257 258

dev.off()
Cesar Taconeli's avatar
Cesar Taconeli committed
259 260
```

261
## Frequências ajustadas pelas duas distribuições
Cesar Taconeli's avatar
Cesar Taconeli committed
262

263
```{r}
Cesar Taconeli's avatar
Cesar Taconeli committed
264

265
# Poisson sem ajuste de covariáveis.
266 267
n <- nrow(seguros)
mediasin <- mean(seguros$nsinist)
268
freqps <- round(n * dpois(0:10, mediasin))
Cesar Taconeli's avatar
Cesar Taconeli committed
269

270 271
# Poisson com covariaveis
pred1 <- predict(mPo, type = "response")
Cesar Taconeli's avatar
Cesar Taconeli committed
272
intervalo <- 0:10
273
matprob <- matrix(0, nrow = nrow(seguros), ncol = length(intervalo))
274
probpois <- function(interv, taxa) dpois(intervalo, taxa)
275
for (i in 1:nrow(seguros)) {
276 277
    matprob[i, ] <- probpois(interv = intervalo, taxa = pred1[i])
}
Cesar Taconeli's avatar
Cesar Taconeli committed
278
pbarra <- colMeans(matprob)
279
freqpsaj <- round(n * pbarra)
Cesar Taconeli's avatar
Cesar Taconeli committed
280

281
# Binomial negativa sem covariaveis.
282
ajustenb <- glm.nb(nsinist ~ 1, data = seguros)
Cesar Taconeli's avatar
Cesar Taconeli committed
283 284 285

media <- exp(coefficients(ajustenb))
shape <- ajustenb$theta
286
freqbn <- round(n * dnbinom(0:10, mu = media, size = shape))
Cesar Taconeli's avatar
Cesar Taconeli committed
287

288 289
# Binomial negativa com covariaveis
pred2 <- predict(mBNo, type = "response")
Cesar Taconeli's avatar
Cesar Taconeli committed
290 291

intervalo <- 0:10
292
matprob <- matrix(0, nrow = nrow(seguros), ncol = length(intervalo))
293 294 295 296
probnb <- function(interv, media, shape) {
    dnbinom(intervalo, mu = media,
            size = shape)
}
297
for (i in 1:nrow(seguros)) {
298 299 300
    matprob[i, ] <- probnb(interv = intervalo, media = pred2[i],
                           shape = mBNo$theta)
}
Cesar Taconeli's avatar
Cesar Taconeli committed
301
pbarra <- colMeans(matprob)
302
frebnaj <- round(n * pbarra)
303
ams <- c(table(seguros$nsinist), rep(0, 5))
Cesar Taconeli's avatar
Cesar Taconeli committed
304 305
matfreq <- rbind(ams, freqps, freqpsaj, freqbn, frebnaj)
colnames(matfreq) <- 0:10
306 307 308 309 310
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")
Cesar Taconeli's avatar
Cesar Taconeli committed
311 312 313 314
```

## Comparação dos ajustes

315 316 317 318 319 320 321 322 323 324 325 326
```{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()
Cesar Taconeli's avatar
Cesar Taconeli committed
327
```
328 329 330 331

```{r, include=FALSE}
detach("package:effects")
```