Calcular a Deviance com a Abordagem MLE

library(bellreg)

## Warning: pacote 'bellreg' foi compilado no R versão 4.4.1

1. Ajustar o modelo usando a abordagem MLE.
   

# Carregar os dados
data(faults)
glimpse(faults)

## Rows: 32
## Columns: 2
## $ nf    <int> 6, 4, 17, 9, 14, 8, 5, 7, 7, 7, 6, 8, 28, 4, 10, 4, 8, 9, 23, 9,…
## $ lroll <int> 551, 651, 832, 375, 715, 868, 271, 630, 491, 372, 645, 441, 895,…

# Ajuste de modelo usando MLE
fit <- bellreg(nf ~ lroll, data = faults, approach = "mle")
summary(fit)

## Call:
## bellreg(formula = nf ~ lroll, data = faults, approach = "mle")
## 
## Coefficients:
##               Estimate     StdErr z.value   p.value    
## (Intercept) 0.98525124 0.33219430  2.9659  0.003018 ** 
## lroll       0.00190933 0.00049003  3.8963 9.766e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## logLik = -88.96139   AIC = 181.9228

2. Extrair a log-verossimilhança do modelo ajustado.
   

# Obter a log-verossimilhança do modelo ajustado
loglik <- fit$logLik

3. Calculo a deviance.

deviance <- -2 * loglik; 
deviance

## [1] 177.9228

4. Calculo dos resíduos da deviance

# Obter os valores ajustados
fitted_values <- fitted(fit)

# Obter os valores observados
observed_values <- faults$nf

# Calcular os resíduos de deviance
deviance_residuals <- sign(observed_values - fitted_values) *
  sqrt(2 * (observed_values * log(observed_values / fitted_values) - 
              (observed_values - fitted_values)))

# Exibir os resíduos de deviance
deviance_residuals

##  [1] -0.62708707 -1.95737700  1.02521399  1.37445341  1.03041878 -1.75734470
##  [7]  0.23454378 -0.66780711  0.06108277  0.63593245 -1.12033896  0.68449792
## [13]  3.05256187 -1.02786186  0.28513094 -1.18263079  0.16087535 -1.27078567
## [19]  1.89221065  0.51611165 -0.48104724 -1.52503916 -0.12828355  0.15098899
## [25] -0.61142617  3.05787832  1.70680081 -0.36339684  0.04100604 -2.73821299
## [31] -2.02069290 -1.87714065

5. Calculo dos resíduos da deviance padronizada

# Calcular a variância dos resíduos (estimativa)
phi_hat <- sum((observed_values - fitted_values)^2 / fitted_values) /
  (length(observed_values) - length(coef(fit)))

# Calcular os resíduos de deviance padronizados
res_deviance_std <- deviance_residuals / sqrt(phi_hat)

# Exibir os resíduos de deviance padronizados
res_deviance_std

##  [1] -0.43051427 -1.34379858  0.70384045  0.94360388  0.70741369 -1.20647035
##  [7]  0.16102141 -0.45846980  0.04193517  0.43658688 -0.76914662  0.46992855
## [13]  2.09567616 -0.70565829  0.19575103 -0.81191185  0.11044580 -0.87243284
## [19]  1.29905991  0.35432627 -0.33025350 -1.04698556 -0.08807054  0.10365852
## [25] -0.41976258  2.09932606  1.17177044 -0.24948294  0.02815189 -1.87986613
## [31] -1.38726687 -1.28871390

6. Calculo dos resíduos de pearson padronizado

# Calcular os resíduos de Pearson
pearson_residuals <- (observed_values - fitted_values) / sqrt(fitted_values)

# Calcular a variância dos resíduos (estimativa)
phi_hat <- sum((observed_values - fitted_values)^2 / fitted_values) / 
  (length(observed_values) - length(coef(fit)))

# Calcular os resíduos de Pearson padronizados
res_pearson_std <- pearson_residuals / sqrt(phi_hat)

# Exibir os resíduos de Pearson padronizados
res_pearson_std

##  [1] -0.41394059 -1.19047226  0.73631863  1.03200443  0.74400513 -1.10795360
##  [7]  0.16396410 -0.44104380  0.04209810  0.45598893 -0.72011192  0.49097087
## [13]  2.35760291 -0.65613800  0.19880686 -0.74814420  0.11151811 -0.82027652
## [19]  1.40075857  0.36526086 -0.32027066 -0.88900427 -0.08745389  0.10500494
## [25] -0.40663295  2.52008664  1.26889842 -0.24492449  0.02822508 -1.59037987
## [31] -1.26667707 -1.10954855

par(mfrow = c(1, 2))
# Gráfico 1: Resíduos de Pearson Padronizados vs Valores ajustados
plot(fitted_values, res_pearson_std, 
     xlab = "Valor ajustado", 
     ylab = "Res. Pearson Std.", 
     main = "", 
     pch = 16)
abline(h = 0, col = "red")

# Gráfico 2: Resíduos de Deviance Padronizados vs Valores ajustados
plot(fitted_values, res_deviance_std, 
     xlab = "Valor ajustado", 
     ylab = "Res. Deviance Std.", 
     main = "", 
     pch = 16)
abline(h = 0, col = "red")