Esse conjunto de scripts tem por finalidade mostrar o impacto de não levar em conta o efeito da superdispersão nas inferências de um MLG.
Vamos usar a base Ornstein, do pacote car. Vamos considerar essa amostra como uma população de referência e simular amostras a partir dela. O objetivo é estimar o efeito de uma das covariáveis (escolhi, arbitrariamente, log (assets)). Os passos são os seguintes:
1- Ajustar o MLG com todos os dados e esyimar o parâmetro de interesse (vamos considerar essa estimativa como o parâmetro a ser estimado).
2 - Seleção aleatória (e sem reposição) de uma amostra de tamanho n=70 da população de referência;
3 - Ajuste do MLG com distribuição Poisson para os dados amostrados e obtenção do intervalo de confiança (95%) para o efeito de log(assets);
4 - Verificar se o intervalo obtido no passo 3 contém ou não o parâmetro de interesse;
5 - Repetição dos passos 2-4 por 1000 vezes;
6- Estimação da taxa de cobertura do intervalo calculando a proporção de ICs que contém o parâmetro de interesse.
O algoritmo será repetido ajustando, no passo 3, modelos Binomial Negativo, quasi poisson e poisson com IC bootstrap.
## assets sector nation interlocks
## 1 147670 BNK CAN 87
## 2 133000 BNK CAN 107
## 3 113230 BNK CAN 94
## 4 85418 BNK CAN 48
## 5 75477 BNK CAN 66
## 6 40742 FIN CAN 69## assets sector nation interlocks
## Min. : 62 MIN :54 CAN:117 Min. : 0.0
## 1st Qu.: 519 MAN :48 OTH: 18 1st Qu.: 3.0
## Median : 1397 AGR :47 UK : 17 Median : 9.0
## Mean : 5978 FIN :22 US : 96 Mean : 13.6
## 3rd Qu.: 4326 MER :20 3rd Qu.: 18.0
## Max. :147670 WOD :19 Max. :107.0
## (Other):38ajuste1 <- glm(interlocks ~ log(assets) + nation + sector, family = poisson, data = Ornstein)
summary(ajuste1) ##
## Call:
## glm(formula = interlocks ~ log(assets) + nation + sector, family = poisson,
## data = Ornstein)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.711 -2.316 -0.459 1.282 6.285
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8394 0.1366 -6.14 8.1e-10 ***
## log(assets) 0.4514 0.0170 26.58 < 2e-16 ***
## nationOTH -0.1070 0.0744 -1.44 0.15030
## nationUK -0.3872 0.0895 -4.33 1.5e-05 ***
## nationUS -0.7724 0.0496 -15.56 < 2e-16 ***
## sectorBNK -0.1665 0.0958 -1.74 0.08204 .
## sectorCON -0.4893 0.2132 -2.29 0.02174 *
## sectorFIN -0.1116 0.0757 -1.47 0.14046
## sectorHLD -0.0149 0.1192 -0.13 0.90051
## sectorMAN 0.1219 0.0761 1.60 0.10949
## sectorMER 0.0616 0.0867 0.71 0.47760
## sectorMIN 0.2498 0.0689 3.63 0.00029 ***
## sectorTRN 0.1518 0.0789 1.92 0.05445 .
## sectorWOD 0.4983 0.0756 6.59 4.4e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3737.0 on 247 degrees of freedom
## Residual deviance: 1547.1 on 234 degrees of freedom
## AIC: 2473
##
## Number of Fisher Scoring iterations: 5## log(assets)
## 0.4514set.seed(98)
ICs <- matrix(0, nrow = 1000, ncol = 2)
for (i in 1:1000){
amostra <- sample(1:nrow(Ornstein), 70)
Ornstein_Novo <- Ornstein[amostra,]
ajuste_novo <- glm(interlocks ~ log(assets) + nation + sector, family = poisson, data = Ornstein_Novo)
ICs[i,] <- confint.default(ajuste_novo)[2,]
}
indica_cobert <- function(interval) ifelse(parametro > interval[1] & parametro < interval[2], 1, 0)cobert <- apply(ICs, 1, indica_cobert)
mean(cobert) # Taxa de cobertura estimada para o modelo Poisson.## [1] 0.58ICs <- matrix(0, nrow = 1000, ncol = 2)
for (i in 1:1000){
amostra <- sample(1:nrow(Ornstein), 70)
Ornstein_Novo <- Ornstein[amostra,]
ajuste_novo <- glm.nb(interlocks ~ log(assets) + nation + sector, data = Ornstein_Novo)
ICs[i,] <- confint.default(ajuste_novo)[2,]
}
cobert <- apply(ICs, 1, indica_cobert)
mean(cobert) # Taxa de cobertura estimada para o modelo binomial negativo.## [1] 0.954ICs <- matrix(0, nrow = 1000, ncol = 2)
for (i in 1:1000){
amostra <- sample(1:nrow(Ornstein), 70)
Ornstein_Novo <- Ornstein[amostra,]
ajuste_novo <- glm(interlocks ~ log(assets) + nation + sector, family = quasipoisson, data = Ornstein_Novo)
ICs[i,] <- confint.default(ajuste_novo)[2,]
}
cobert <- apply(ICs, 1, indica_cobert)
mean(cobert) # Taxa de cobertura estimada para o modelo quase Poisson.## [1] 0.942ICs <- matrix(0, nrow = 1000, ncol = 2)
for (i in 1:1000){
amostra <- sample(1:nrow(Ornstein), 70)
Ornstein_Novo <- Ornstein[amostra,]
ajuste_novo <- glm(interlocks ~ log(assets) + nation + sector, family = poisson, data = Ornstein_Novo)
boot_pois <- Boot(ajuste_novo)
ICs[i,] <- confint(boot_pois, type = 'perc')[2,]
}
cobert <- apply(ICs, 1, indica_cobert)
mean(cobert) # Taxa de cobertura = 0.98.## [1] 0.976