MLG - Atividade: teste de hipóteses
Mauro Citro Lalucci
2024-03-06
Para os seguintes exercícios, foram utilizados os dados “sinistros.csv” apresentados em aula.
Pacotes:
require(statmod)## Carregando pacotes exigidos: statmod## Warning: package 'statmod' was built under R version 4.2.3require(ISLR)## Carregando pacotes exigidos: ISLR## Warning: package 'ISLR' was built under R version 4.2.3require(car)## Carregando pacotes exigidos: car## Warning: package 'car' was built under R version 4.2.3## Carregando pacotes exigidos: carData## Warning: package 'carData' was built under R version 4.2.3require(hnp)## Carregando pacotes exigidos: hnp## Warning: package 'hnp' was built under R version 4.2.3## Carregando pacotes exigidos: MASSrequire(pscl)## Carregando pacotes exigidos: pscl## Warning: package 'pscl' was built under R version 4.2.3## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.require(coefplot)## Carregando pacotes exigidos: coefplot## Warning: package 'coefplot' was built under R version 4.2.3## Carregando pacotes exigidos: ggplot2require(effects)## Carregando pacotes exigidos: effects## Warning: package 'effects' was built under R version 4.2.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.require(sandwich)## Carregando pacotes exigidos: sandwich## Warning: package 'sandwich' was built under R version 4.2.3require(lmtest)## Carregando pacotes exigidos: lmtest## Warning: package 'lmtest' was built under R version 4.2.3## Carregando pacotes exigidos: zoo## Warning: package 'zoo' was built under R version 4.2.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericSobre o dados:
Regressão para dados de contagens. Informações referentes a 500 portadores de seguro de automóvel de uma seguradora particular. Foram filtrados apenas os segurados cujas apólices têm mais de cinco anos. As variáveis são as seguintes:
idade: idade do segurado (em anos);
sexo: Masc para masculino e
Fem para feminino;
usop: uso principal do veículo. Cidade para uso
urbano e Estrada para uso rodoviário;
anosest: escolaridade do
segurado, em anos de estudo;
claims: número de sinistros produzidos
pelo segurado nos últimos cinco anos.
Dados sinistros.csv
dados <- read.csv2("C:/Users/crist/Downloads/sinistros.csv")[,-1]
head(dados, 10) ### Dez primeiras linhas da base.## idade sexo usop anosest claims
## 1 41 Fem Estrada 10 3
## 2 39 Masc Cidade 16 0
## 3 46 Masc Estrada 5 1
## 4 45 Fem Cidade 9 0
## 5 41 Fem Cidade 12 1
## 6 33 Masc Estrada 8 7
## 7 50 Masc Estrada 6 1
## 8 47 Masc Estrada 4 0
## 9 44 Masc Estrada 13 3
## 10 41 Masc Estrada 13 1summary(dados) ### Algumas descritivas dos dados.## idade sexo usop anosest
## Min. :16.00 Length:500 Length:500 Min. : 4.00
## 1st Qu.:35.00 Class :character Class :character 1st Qu.: 7.00
## Median :40.00 Mode :character Mode :character Median :10.00
## Mean :39.85 Mean :10.31
## 3rd Qu.:45.00 3rd Qu.:13.00
## Max. :61.00 Max. :17.00
## claims
## Min. :0.000
## 1st Qu.:1.000
## Median :1.000
## Mean :1.698
## 3rd Qu.:2.000
## Max. :9.000Exercício 1 - Aplicar o teste de razão de verossimilhanças para testar a hipótese: hipótese nula H0: B_sexoMasc = B_anosest = 0.
ajH0 <- glm(claims ~ idade + usop , family = poisson(link = 'log'),
data = dados)
ajuste2 <- glm(claims ~ idade + sexo + usop + anosest, family = poisson(link = 'log'),
data = dados)
l0 <- as.numeric(logLik(ajH0)); l0 #log-verossimilhança maximizada## [1] -784.3994l1 <- as.numeric(logLik(ajuste2)); l1 #log-verossimilhança maximizada## [1] -784.2589Lambda <- -2*(l0-l1); Lambda #estatastíca do teste da razão de verossimilhanças## [1] 0.2811146Estamos testando 2 parâmetros, a distribuição de referência para o teste é a qui-quadrado com 2 grau de liberdade. Vamos obter o valor crítico para um nível de significância de 5%:
qchisq(0.05, df = 2, lower.tail = FALSE)## [1] 5.991465Como o valor da estatística do teste não excede o valor crítico, a hipótese nula (não efeito de sexo e escolaridade) não é rejeitada ao nível de significância de 5%. Vamos calcular o p-valor do teste.
pchisq(Lambda, df = 2, lower.tail = FALSE)## [1] 0.86887390.8688739 < 5.991465 Não rejeitamos a hipótese nula
Exercício 2 - Aplicar o método de Wald no teste da hipótese nula: H0: B_sexoMasc = B_idade = 0.
#Começamos extraindo estimativas pontuais e a parte da matriz de variâncias e covariâncias referentes aos dois parâmetros sob teste.
B0 <- matrix(coef(ajuste2)[c(2,3)]); B0## [,1]
## [1,] -0.05870844
## [2,] -0.01638627VarB0 <- vcov(ajuste2)[c(2,3),c(2,3)]; VarB0## idade sexoMasc
## idade 2.527460e-05 -1.217083e-05
## sexoMasc -1.217083e-05 6.051997e-03Estatística do teste de Wald
W <- t(B0) %*% solve(VarB0) %*% B0
W## [,1]
## [1,] 136.6992qchisq(0.05, df = 2, lower.tail = FALSE) #valor crítico para um nível de significância de 5%, considerando 2 graus de liberdade## [1] 5.991465pchisq(W, df = 2, lower.tail = FALSE) #valor p-valor do teste## [,1]
## [1,] 2.070826e-30ajH0 <- update(ajuste2, claims ~.-sexo - idade)Usando a função waldtest do pacote lmtest:
ajH0 <- update(ajuste2, claims ~.-sexo - idade)
waldtest(ajH0, ajuste2, test = 'Chisq')## Wald test
##
## Model 1: claims ~ usop + anosest
## Model 2: claims ~ idade + sexo + usop + anosest
## Res.Df Df Chisq Pr(>Chisq)
## 1 497
## 2 495 2 136.7 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como o valor de probabilidade foi quase nulo (2.2e-16), rejeitamos a hipótese nula.
Exercício 3 - Construir o gráfico do perfil da verossimilhança e obter os intervalos de confiança com níveis de 90 e 95% para B_sexoM. (Dica: converta a variável sexo para numérica para construir o gráfico do perfil)
confint(ajuste2, parm = "sexoMasc", level = 0.95)ajuste <- glm(claims ~ offset(-0.016 * sexo) + usop + anosest + sexo, family = poisson(link = 'log'),data=dados)Ao tentar ajustar o modelo de regressão para a variável sexo, o R retornou um erro. Pois a variável sexo é categórica e não numérica. Para corrigir isso, é necessário transformar a variável sexo em uma variável dummy.
Transformando a variável sexo em uma variável dummy
dados$sexo <- ifelse(dados$sexo == "Masc", 1, 0)Agora podemos ajustar o modelo de regressão para a variável sexo
ajuste <- glm(claims ~ offset(-0.016 * sexo) + usop + anosest + sexo, family = poisson(link = 'log'),data=dados)
B_sexo_grid <- seq(-1, 1, length.out = 20)vet_logLik <- numeric()
for(i in 1:length(B_sexo_grid)){
B_sexo <- B_sexo_grid[i]
ajuste <- glm(claims ~ offset(B_sexo*sexo) + usop + anosest + idade, family = poisson(link = 'log'),
data = dados)
vet_logLik[i] <- logLik(ajuste)
}Estatística do teste da razão de verossimilhanças para cada valor de B_idade_grid.
par(cex = 1.4, las = 1)
plot(B_sexo_grid, vet_logLik, type = 'b', pch = 20, xlab = 'B_idade', ylab = 'Log-verossimilhan?a')
abline(h = logLik(ajuste2) - 3.84/2)
abline(v = B_sexo_grid[which.max(vet_logLik)], col = 'red', lwd = 2, lty = 2)
Calculando o IC baseado no perfil da verossimilhança usando a função confint:
confint(ajuste2, level = 0.95) #95%## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 2.28162415 3.11523675
## idade -0.06855816 -0.04885103
## sexoMasc -0.16817436 0.13685523
## usopEstrada 0.10861796 0.40834260
## anosest -0.02280212 0.01383107confint(ajuste2, level = 0.99) #99%## Waiting for profiling to be done...## 0.5 % 99.5 %
## (Intercept) 2.14968891 3.24532857
## idade -0.07165167 -0.04575198
## sexoMasc -0.21560746 0.18534609
## usopEstrada 0.06188893 0.45586413
## anosest -0.02857790 0.01957455Plotando os limites de confiança sobre o gráfico.
plot(B_sexo_grid, vet_logLik, type = 'b', pch = 20, xlab = 'B_idade', ylab = 'Log-verossimilhanca')
abline(v = c(-0.16817436, -0.13685523), lty = 2, col = 'red')