Modelos de Regressão para Dados de Contagem com o R

github.com/leg-ufpr/MRDCr

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:

# Pacotes necessários.
library(lattice)
library(MASS)
library(effects)
library(knitr)
library(MRDCr)

Verificação do conteúdo e a estrutura dos dados

# Dez primeiras linhas da base.
head(seguros, 10)
##     tipo idade sexo    civil  valor expos nsinist
## 1  hatch    59    F   Casado 24.624  0.50       1
## 2   mono    45    M Solteiro 23.391  0.70       0
## 3    suv    42    M   Casado 86.639  0.79       0
## 4  sedan    63    F   Casado 77.545  0.01       0
## 5   mono    36    F   Casado 25.934  0.51       0
## 6  hatch    33    M   Casado 42.704  0.79       0
## 7  hatch    35    M   Casado 16.986  0.81       0
## 8  sedan    63    M   Casado 19.055  0.01       0
## 9    suv    54    M   Casado 52.374  0.76       0
## 10   suv    32    M   Casado 52.524  0.79       0
str(seguros)
## 'data.frame':    16483 obs. of  7 variables:
##  $ tipo   : Factor w/ 6 levels "hatch","mono",..: 1 2 5 4 2 1 1 4 5 5 ...
##  $ idade  : int  59 45 42 63 36 33 35 63 54 32 ...
##  $ sexo   : Factor w/ 2 levels "F","M": 1 2 2 1 1 2 2 2 2 2 ...
##  $ civil  : Factor w/ 4 levels "Casado","Divorciado",..: 1 3 1 1 1 1 1 1 1 1 ...
##  $ valor  : num  24.6 23.4 86.6 77.5 25.9 ...
##  $ expos  : num  0.5 0.7 0.79 0.01 0.51 0.79 0.81 0.01 0.76 0.79 ...
##  $ nsinist: int  1 0 0 0 0 0 0 0 0 0 ...
seguros$lexpo <- log(seguros$expos)

Análise descritiva da distribuição do número de sinistros

# Distribuição do números de sinistros.
tb <- table(seguros$nsinist)
tb
## 
##     0     1     2     3     4     5 
## 15689   602   166    22     3     1
barchart(tb, horizontal = FALSE)

# Taxa de sinistros na amostra.
taxageral <- sum(seguros$nsinist)/sum(seguros$expos)
taxageral
## [1] 0.09517499
tab <- aggregate(cbind(expos, nsinist) ~ sexo,
                 data = seguros, FUN = sum)
taxa <- with(tab, nsinist/expos)
tab <- cbind(tab, taxa)

# Distribuição do número de sinistros por sexo.
kable(tab, align = "c",
      caption = "**Taxa de sinistros segundo sexo**")
Taxa de sinistros segundo sexo
sexo expos nsinist taxa
F 4363.92 461 0.1056390
M 6321.66 556 0.0879516 
seguros$idadecat <- cut(seguros$idade,
                       breaks = c(18, 30, 60, 92),
                       include.lowest = TRUE)
tab <- aggregate(cbind(expos, nsinist) ~ idadecat,
                 data = seguros, FUN = sum)
taxa <- with(tab, nsinist/expos)
tab <- cbind(tab, taxa)

# Distribuição do número de sinistros por sexo.
kable(tab, align = "c",
      caption = "**Taxa de sinistros segundo idade**")
Taxa de sinistros segundo idade
idadecat expos nsinist taxa
[18,30] 462.91 46 0.0993714
(30,60] 6909.21 661 0.0956694
(60,92] 3311.64 310 0.0936092 
tabidsex <- aggregate(cbind(expos, nsinist) ~ sexo + idadecat,
                      data = seguros, FUN = sum)
taxa <- with(tabidsex, nsinist/expos)
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**")
Taxa de sinistros segundo sexo e idade
sexo idadecat expos nsinist taxa
F [18,30] 273.34 29 0.1060950
M [18,30] 189.57 17 0.0896766
F (30,60] 3035.82 318 0.1047493
M (30,60] 3873.39 343 0.0885529
F (60,92] 1054.74 114 0.1080835
M (60,92] 2256.90 196 0.0868448

Regressão usando o modelo log-linear Poisson

seguros <- na.omit(seguros)
mP <- glm(nsinist ~ sexo + idade + I(idade^2) + valor +
              offset(log(expos)),
          data = seguros, family = poisson)
summary(mP)
## 
## Call:
## glm(formula = nsinist ~ sexo + idade + I(idade^2) + valor + offset(log(expos)), 
##     family = poisson, data = seguros)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5778  -0.3935  -0.3646  -0.2971   5.2731  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.4583436  0.3809057  -3.829 0.000129 ***
## sexoM       -0.1796330  0.0638541  -2.813 0.004905 ** 
## idade       -0.0374658  0.0149912  -2.499 0.012448 *  
## I(idade^2)   0.0003276  0.0001416   2.314 0.020644 *  
## valor        0.0057360  0.0014905   3.848 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 6530.2  on 16478  degrees of freedom
## Residual deviance: 6502.0  on 16474  degrees of freedom
## AIC: 8229
## 
## Number of Fisher Scoring iterations: 6
# 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
## [1] 2.518473
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 modelo - gráficos.
par(mfrow = c(2, 2))
plot(mP)

par(mfrow = c(1, 1))
envelope(mP)

Ajuste do modelo associando um parâmetro ao termo offset (log-exposicao)

mPo <- glm(nsinist ~ sexo + idade + I(idade^2) + valor +
               log(expos),
           data = seguros, family = poisson)
summary(mPo)
## 
## Call:
## glm(formula = nsinist ~ sexo + idade + I(idade^2) + valor + log(expos), 
##     family = poisson, data = seguros)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5068  -0.3714  -0.3489  -0.3284   5.4639  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.8967565  0.3819147  -4.966 6.82e-07 ***
## sexoM       -0.1888165  0.0638419  -2.958  0.00310 ** 
## idade       -0.0327982  0.0149723  -2.191  0.02848 *  
## I(idade^2)   0.0003002  0.0001414   2.123  0.03376 *  
## valor        0.0045550  0.0014881   3.061  0.00221 ** 
## log(expos)   0.1865318  0.0410945   4.539 5.65e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 6319.8  on 16478  degrees of freedom
## Residual deviance: 6275.2  on 16473  degrees of freedom
## AIC: 8004.2
## 
## Number of Fisher Scoring iterations: 6
anova(mP, mPo, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: nsinist ~ sexo + idade + I(idade^2) + valor + offset(log(expos))
## Model 2: nsinist ~ sexo + idade + I(idade^2) + valor + log(expos)
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1     16474     6502.0                          
## 2     16473     6275.2  1   226.78 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regressão usando a distribuição binomial negativa

mBNo <- glm.nb(nsinist ~ sexo + idade + I(idade^2) + valor +
                   log(expos), data = seguros)
summary(mBNo)
## 
## Call:
## glm.nb(formula = nsinist ~ sexo + idade + I(idade^2) + valor + 
##     log(expos), data = seguros, init.theta = 0.1203144242, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4198  -0.3306  -0.3141  -0.2984   2.6758  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.8575829  0.4809945  -3.862 0.000112 ***
## sexoM       -0.1897891  0.0791377  -2.398 0.016475 *  
## idade       -0.0342408  0.0188090  -1.820 0.068690 .  
## I(idade^2)   0.0003144  0.0001774   1.772 0.076428 .  
## valor        0.0045735  0.0018819   2.430 0.015086 *  
## log(expos)   0.1957247  0.0480051   4.077 4.56e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.1203) family taken to be 1)
## 
##     Null deviance: 3201.7  on 16478  degrees of freedom
## Residual deviance: 3170.7  on 16473  degrees of freedom
## AIC: 7430.7
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.1203 
##           Std. Err.:  0.0105 
## 
##  2 x log-likelihood:  -7416.7460

Diagnóstico do ajuste

# Diagnóstico do modelo - gráficos.
par(mfrow = c(2, 2))
plot(mBNo)

dadosnb3 <-
    seguros[, c("sexo", "idade", "valor", "expos", "nsinist")]
dadosnb3$lexpo <- log(seguros$expos)

mBNo <- glm.nb(nsinist ~ 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)
}
par(mfrow = c(1, 1))
envelope(mBNo)

Explorando os efeitos das covariáveis

intervalos <- confint(mBNo)
## Waiting for profiling to be done...
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**"))
Estimativas pontuais e intervalos de confiança - Modelo Binomial Negativo
Estimativa pontual 2.5 % 97.5 %
(Intercept) -1.85758 -2.81686 -0.90742
sexoM -0.18979 -0.34453 -0.03493
idade -0.03424 -0.07148 0.00332
I(idade^2) 0.00031 -0.00004 0.00067
valor 0.00457 0.00078 0.00835
lexpo 0.19572 0.10164 0.29403

Gráficos de efeitos

efeitos <- allEffects(mBNo, given.values = c(lexpo = 0))

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")

dev.off()
## null device 
##           1

Frequências ajustadas pelas duas distribuições

# Poisson sem ajuste de covariáveis.
n <- nrow(seguros)
mediasin <- mean(seguros$nsinist)
freqps <- round(n * dpois(0:10, mediasin))

# Poisson com covariaveis
pred1 <- predict(mPo, type = "response")
intervalo <- 0:10
matprob <- matrix(0, nrow = nrow(seguros), ncol = length(intervalo))
probpois <- function(interv, taxa) dpois(intervalo, taxa)
for (i in 1:nrow(seguros)) {
    matprob[i, ] <- probpois(interv = intervalo, taxa = pred1[i])
}
pbarra <- colMeans(matprob)
freqpsaj <- round(n * pbarra)

# Binomial negativa sem covariaveis.
ajustenb <- glm.nb(nsinist ~ 1, data = seguros)

media <- exp(coefficients(ajustenb))
shape <- ajustenb$theta
freqbn <- round(n * dnbinom(0:10, mu = media, size = shape))

# Binomial negativa com covariaveis
pred2 <- predict(mBNo, type = "response")

intervalo <- 0:10
matprob <- matrix(0, nrow = nrow(seguros), ncol = length(intervalo))
probnb <- function(interv, media, shape) {
    dnbinom(intervalo, mu = media,
            size = shape)
}
for (i in 1:nrow(seguros)) {
    matprob[i, ] <- probnb(interv = intervalo, media = pred2[i],
                           shape = mBNo$theta)
}
pbarra <- colMeans(matprob)
frebnaj <- round(n * pbarra)
ams <- c(table(seguros$nsinist), 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

kable(matfreq, format = "markdown",
      caption = paste("Frequências amostrais e frequências",
                      "ajustadas para o número de sinistros"))
0 1 2 3 4 5 6 7 8 9 10
Amostra 15685 602 166 22 3 1 0 0 0 0 0
Poisson não ajustada por covariáveis 15493 956 30 1 0 0 0 0 0 0 0
Poisson ajustada por covariáveis 15494 954 31 1 0 0 0 0 0 0 0
BN não ajustada por covariáveis 15683 632 122 30 8 2 1 0 0 0 0
BN ajustada por covariáveis 15683 633 121 30 8 2 1 0 0 0 0

Informações da sessão

## Atualizado em 21 de julho de 2016.
## 
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04 LTS
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets 
## [7] methods   base     
## 
## other attached packages:
##  [1] effects_3.1-1       hnp_1.2             sandwich_2.3-4     
##  [4] car_2.1-2           boot_1.3-18         lmtest_0.9-34      
##  [7] zoo_1.7-13          multcomp_1.4-6      TH.data_1.0-7      
## [10] MASS_7.3-45         survival_2.39-5     mvtnorm_1.0-5      
## [13] doBy_4.5-15         rmarkdown_1.0       latticeExtra_0.6-28
## [16] RColorBrewer_1.1-2  lattice_0.20-33     knitr_1.13         
## [19] MRDCr_0.0-2         gsubfn_0.6-6        proto_0.3-10       
## [22] bbmle_1.0.18        devtools_1.12.0    
## 
## loaded via a namespace (and not attached):
##  [1] splines_3.3.1      tcltk_3.3.1        colorspace_1.2-6  
##  [4] testthat_1.0.2     htmltools_0.3.5    mgcv_1.8-12       
##  [7] yaml_2.1.13        nloptr_1.0.4       withr_1.0.2       
## [10] stringr_1.0.0      MatrixModels_0.4-1 codetools_0.2-14  
## [13] memoise_1.0.0      evaluate_0.9       SparseM_1.7       
## [16] quantreg_5.26      pbkrtest_0.4-6     parallel_3.3.1    
## [19] highr_0.6          Rcpp_0.12.5        formatR_1.4       
## [22] lme4_1.1-12        digest_0.6.9       stringi_1.1.1     
## [25] numDeriv_2014.2-1  grid_3.3.1         tools_3.3.1       
## [28] magrittr_1.5       crayon_1.3.2       Matrix_1.2-6      
## [31] minqa_1.2.4        roxygen2_5.0.1     R6_2.1.2          
## [34] nnet_7.3-12        nlme_3.1-128       compiler_3.3.1

Modelos de Regressão para Dados de Contagem com o R

61 RBRAS Departamento de Estatística - UFPR
23 a 25 de Maio de 2016 LEG, PET Estatística
Salvador - Bahia github.com/leg-ufpr/MRDCr