rm(list=ls(all=TRUE))
setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/AULAS/LOGÍSTICA")
dados <- read.table("MichelinNY.txt", sep = ",", header = TRUE)Regressão Logística - Regressão 2
Regressão Logística
Realizando alguns exemplos da regressão logística com dados diferentes
Exemplo 1
A revista é o guia mais famoso de hoteis e restaurantes de Nova York. A inclusão de um restaurante no guia é baseado em um processo de avaliação meticuloso e confidencial no qual inspetores do guia visitam anonimamente restaurantes e hoteis de Nova York e fazem sua avaliação de local.
Já o guia é baseado na submissão de avaliações por clientes que responderam por e-mail ou carta.
Queremos avaliar a probabilidade de um restaurante francês ser incluido no , baseado na avaliação dos clientes no .
A variável resposta InMichelin:
\[y_i = 1, \quad ser \hspace{0.1cm} incluído \hspace{0.1cm} no \hspace{0.1cm} guia \hspace{0.1cm} Michelin \]
\[y_i = 0, \quad não \hspace{0.1cm} ser \hspace{0.1cm} incluído \hspace{0.1cm} no \hspace{0.1cm} guia \hspace{0.1cm} Michelin\] As outras variáveis são:
| Variável \(X_i\) | Descrição |
|---|---|
| Food | Avaliação da comida (nota de 0 a 30) |
| Décor | Avaliação da decoração (nota de 0 a 30) |
| Service | Avaliação do serviço (nota de 0 a 30) |
| Price | Preço (em dólares) do jantar |
Análise Descritiva dos dados
DescTools::Desc(dados, digits = 2, plotit = F)------------------------------------------------------------------------------
Describe dados (data.frame):
data frame: 76 obs. of 6 variables
76 complete cases (100.0%)
Nr ColName Class NAs Levels
1 InMichelin integer .
2 Adress character .
3 Food integer .
4 Decor integer .
5 Service integer .
6 Price integer .
------------------------------------------------------------------------------
1 - InMichelin (integer - dichotomous)
length n NAs unique
76 76 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
0 42 55.26% 44.10% 65.92%
1 34 44.74% 34.08% 55.90%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
2 - Adress (character)
length n NAs unique levels dupes
76 76 0 76 76 n
100.0% 0.0%
level freq perc cumfreq cumperc
1 14 Wall Street 1 1.32% 1 1.32%
2 212th 1 1.32% 2 2.63%
3 26 Seats 1 1.32% 3 3.95%
4 44 1 1.32% 4 5.26%
5 A 1 1.32% 5 6.58%
6 A.O.C. 1 1.32% 6 7.89%
7 A.O.C. Bedford 1 1.32% 7 9.21%
8 Aix 1 1.32% 8 10.53%
9 Alain Ducasse 1 1.32% 9 11.84%
10 Alouette 1 1.32% 10 13.16%
11 Arabelle 1 1.32% 11 14.47%
12 Artisanal 1 1.32% 12 15.79%
... etc.
[list output truncated]
------------------------------------------------------------------------------
3 - Food (integer)
length n NAs unique 0s mean meanCI'
76 76 0 12 0 21.47 20.85
100.0% 0.0% 0.0% 22.10
.05 .10 .25 median .75 .90 .95
17.75 18.00 19.00 21.00 23.00 25.00 27.00
range sd vcoef mad IQR skew kurt
11.00 2.72 0.13 2.97 4.00 0.43 -0.60
value freq perc cumfreq cumperc
1 17 4 5.3% 4 5.3%
2 18 5 6.6% 9 11.8%
3 19 11 14.5% 20 26.3%
4 20 13 17.1% 33 43.4%
5 21 9 11.8% 42 55.3%
6 22 6 7.9% 48 63.2%
7 23 13 17.1% 61 80.3%
8 24 2 2.6% 63 82.9%
9 25 7 9.2% 70 92.1%
10 26 1 1.3% 71 93.4%
11 27 4 5.3% 75 98.7%
12 28 1 1.3% 76 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
4 - Decor (integer)
length n NAs unique 0s mean meanCI'
76 76 0 16 0 19.04 18.18
100.0% 0.0% 0.0% 19.89
.05 .10 .25 median .75 .90 .95
14.00 15.00 16.75 19.00 21.00 25.00 26.25
range sd vcoef mad IQR skew kurt
15.00 3.74 0.20 2.97 4.25 0.50 -0.52
lowest : 12, 13, 14 (4), 15 (10), 16 (3)
highest: 23 (4), 24, 25 (2), 26 (3), 27 (4)
heap(?): remarkable frequency (18.4%) for the mode(s) (= 17)
' 95%-CI (classic)
------------------------------------------------------------------------------
5 - Service (integer)
length n NAs unique 0s mean meanCI'
76 76 0 15 0 19.80 19.03
100.0% 0.0% 0.0% 20.57
.05 .10 .25 median .75 .90 .95
14.75 16.00 18.00 19.00 22.00 23.50 27.00
range sd vcoef mad IQR skew kurt
15.00 3.38 0.17 2.97 4.00 0.43 -0.13
lowest : 13, 14 (3), 15 (2), 16 (7), 17 (5)
highest: 23 (8), 24, 25, 27 (5), 28
heap(?): remarkable frequency (15.8%) for the mode(s) (= 19)
' 95%-CI (classic)
------------------------------------------------------------------------------
6 - Price (integer)
length n NAs unique 0s mean meanCI'
76 76 0 38 0 50.49 45.49
100.0% 0.0% 0.0% 55.49
.05 .10 .25 median .75 .90 .95
27.75 34.50 40.00 46.00 53.00 71.50 89.75
range sd vcoef mad IQR skew kurt
166.00 21.89 0.43 8.90 13.00 2.97 13.92
lowest : 13, 22, 24, 27, 28 (2)
highest: 82, 88, 95 (2), 100, 179
heap(?): remarkable frequency (11.8%) for the mode(s) (= 50)
' 95%-CI (classic)str(dados)'data.frame': 76 obs. of 6 variables:
$ InMichelin: int 0 0 0 1 0 0 1 1 1 0 ...
$ Adress : chr "14 Wall Street" "212th" "26 Seats" "44" ...
$ Food : int 19 17 23 19 23 18 24 23 27 20 ...
$ Decor : int 20 17 17 23 12 17 21 22 27 17 ...
$ Service : int 19 16 21 16 19 17 22 21 27 19 ...
$ Price : int 50 43 35 52 24 36 51 61 179 42 ...Ajuste do modelo
modelo <- glm(InMichelin~Food+Decor+Service+Price, family = binomial(link="logit"), data = dados)Verificando aceitação do modelo:
Fixamos \(\alpha = 0.05\)
phi1 <- summary(modelo)$dispersion
desvio1 <- summary(modelo)$deviance/phi1
q.quadr1 <- qchisq(0.95, desvio1)Temos do resultado acima que: O modelo em análise foi aceito.
Regressão Stepwise
menor <- glm(InMichelin~1, family = binomial(link="logit"), data = dados)
maior <- glm(InMichelin~Food+Decor+Service+Price, family = binomial(link="logit"), data = dados)
step(menor, scope = list(upper = maior), direction = "both")Start: AIC=106.51
InMichelin ~ 1
Df Deviance AIC
+ Price 1 81.806 85.806
+ Decor 1 83.804 87.804
+ Food 1 84.523 88.523
+ Service 1 92.519 96.519
<none> 104.515 106.515
Step: AIC=85.81
InMichelin ~ Price
Df Deviance AIC
+ Food 1 76.719 82.719
<none> 81.806 85.806
+ Decor 1 79.963 85.963
+ Service 1 81.801 87.801
- Price 1 104.515 106.515
Step: AIC=82.72
InMichelin ~ Price + Food
Df Deviance AIC
+ Service 1 73.497 81.497
<none> 76.719 82.719
+ Decor 1 75.770 83.770
- Food 1 81.806 85.806
- Price 1 84.523 88.523
Step: AIC=81.5
InMichelin ~ Price + Food + Service
Df Deviance AIC
<none> 73.497 81.497
+ Decor 1 72.357 82.357
- Service 1 76.719 82.719
- Food 1 81.801 87.801
- Price 1 84.510 90.510
Call: glm(formula = InMichelin ~ Price + Food + Service, family = binomial(link = "logit"),
data = dados)
Coefficients:
(Intercept) Price Food Service
-10.0103 0.1003 0.4921 -0.2812
Degrees of Freedom: 75 Total (i.e. Null); 72 Residual
Null Deviance: 104.5
Residual Deviance: 73.5 AIC: 81.5Como resultado, observamos que o modelo final foi dado por:
\[log \left ( \frac{\pi_i}{1-\pi_i} \right ) = \beta_1 + \beta_2. preço_i + \beta_3. comida_i + \beta_4. serviço_i\]
fit <- glm(InMichelin ~ Price + Food + Service, family = binomial(link = "logit"),data = dados)
phi11 <- summary(fit)$dispersion
desvio11 <- summary(fit)$deviance/phi11
q.quadr11 <- qchisq(0.95, desvio11)| Modelo | |
|---|---|
| Desvio/\(\phi\) | 73.4972023 |
| \(\chi^2\) | 94.5103411 |
Temos do resultado acima que: O modelo em análise foi aceito.
ajuste<-summary(fit)
knitr::kable(ajuste$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -10.0102569 | 3.0925493 | -3.236895 | 0.0012084 |
| Price | 0.1002969 | 0.0338208 | 2.965540 | 0.0030215 |
| Food | 0.4920919 | 0.1741115 | 2.826303 | 0.0047089 |
| Service | -0.2811900 | 0.1594775 | -1.763196 | 0.0778675 |
\[log \left ( \frac{\pi_i}{1-\pi_i} \right ) = -10.01 + 0.1. preço_i + 0.492. comida_i -0.281. serviço_i\]
Análise de Diagnósticos
Teste de Normalidade
rd <- boot::glm.diag(fit)$rd
teste1 <- nortest::lillie.test(rd)
teste2 <- shapiro.test(rd)| Teste 1 | Teste 2 | |
|---|---|---|
| p-valor | 2.3922984^{-4} | 0.0033125 |
q1 <- ggplot2::ggplot(dados, ggplot2::aes(sample=rd))+
ggplot2::stat_qq(size=4, color="#e7ad52") + ggplot2::ylim(-4, 4) + ggplot2::xlim(-4, 4) +
ggplot2::stat_qq_line(colour = '#001132') +
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
q1
É observado através da análise gráfica e dos testes que os resíduos não são distribuidos normalmente.
Ligação
va <- fitted(fit)
g1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va,y=InMichelin))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,1) + ggplot2::xlim(0,1) + ggplot2::ggtitle("Ligação")
G1 <- g1+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
G1
Variância
q1 <- ggplot2::ggplot(dados, ggplot2::aes(sample = rd))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Variância")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
q1
Pontos aberrantes
indice=as.vector(1:length(dados$InMichelin))
abe <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=rd)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = -2, slope = 0, color='#001132')+ ggplot2::geom_abline(intercept = 2, slope = 0, color='#001132') + ggplot2::ylim(-4,4) + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
abeConforme vemos a análise gráfica, temos 2 pontos aberrantes, que são os pontos 11 e 14. Sendo as observações:
dados[which(rd>2 | rd< -2),] InMichelin Adress Food Decor Service Price
11 0 Arabelle 25 26 27 71
14 0 Atelier 27 25 27 95Pontos de Alavanca
h<- boot::glm.diag(fit)$h
h3 <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=h)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = 2*(length(coef(fit))/length(dados$InMichelin)), slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
h3
Conforme o gráfico vemos que tem alguns pontos de alavanca, mais precisamente 4 pontos. Que são as observações de número:11, 31, 44 e 46.
dados[which(h>2*length(coef(fit))/length(dados$InMichelin)),] InMichelin Adress Food Decor Service Price
11 0 Arabelle 25 26 27 71
31 0 La Bergamote 26 13 14 13
44 1 Le Bilboquet 20 15 15 50
46 0 Le Charlot 19 15 15 50Pontos de influência
DC<-boot::glm.diag(fit)$cook
dc <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=DC)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = qchisq(0.05,length(coef(fit)))/2, slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
dc
É observado que temos dois pontosde influência, que são as observações de número: 14 e 31. Sendo os valores observados:
dados[which(DC>qchisq(0.05,length(coef(fit)))/2),] InMichelin Adress Food Decor Service Price
14 0 Atelier 27 25 27 95
31 0 La Bergamote 26 13 14 13Interpretação de cada parâmetro
knitr::kable(exp(fit$coefficients))| x | |
|---|---|
| (Intercept) | 0.0000449 |
| Price | 1.1054991 |
| Food | 1.6357344 |
| Service | 0.7548849 |
Curva ROC
fitt <- fitted(fit)
pred <- ROCR::prediction(fitt, dados$InMichelin)
perf <- ROCR::performance(pred, "tpr", "fpr") # Escolha do ponto de corte
area <- ROCR::performance(pred, "auc") # Calcula a area sob a curva ROC
plot(perf) # Constroi o gráfico da curva ROC
Outra forma
Epi::ROC(form = InMichelin~Food+Service+Price, data = dados, plot = "ROC", MX = TRUE)
Conforme o gráfico é observado que o ponto de corte é de 0.283
Construção da matriz Confusão
threshold = 0.283
y.hat <- ifelse(fitt>threshold,1,0)
mat.conf <- caret::confusionMatrix(as.factor(y.hat), as.factor(dados$InMichelin), positive = "1")
knitr::kable(mat.conf$table)| 0 | 1 | |
|---|---|---|
| 0 | 31 | 2 |
| 1 | 11 | 32 |
Exemplo 2
Como ilustração, vamos considerar os dados de um experimento desenvolvido para avaliar a influência da quantidade de ar inspirado na ocorrência de vaso-constrição na pele dos dedos da mão (Finney, 1978; Pregibon, 1981). Os dados do experimento estão no arquivo pregibon.txt.
A resposta, nesse exemplo, é a ocorrência (\(Y = 1\)) ou ausência (\(Y=0\)) de compressão de vasos e as covariáveis são o logaritmo do volume e o logaritmo da razão de ar inspirado. Vamos supor para a \(i\)-ésima unidade experimental que \(Y_i \sim Be(\pi_i)\), em que:
\[log \left ( \frac{\pi_i}{1-\pi_i} \right ) = \beta_1 + \beta_2. volume_i + \beta_3. razão_i\]
rm(list=ls(all=TRUE))
setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/AULAS/LOGÍSTICA")
dados <- read.table("pregibon.txt", sep = "", header = TRUE)Análise Descritiva dos dados
DescTools::Desc(dados, digits = 2, plotit = F)------------------------------------------------------------------------------
Describe dados (data.frame):
data frame: 39 obs. of 3 variables
39 complete cases (100.0%)
Nr ColName Class NAs Levels
1 Resposta integer .
2 Volume numeric .
3 Razao numeric .
------------------------------------------------------------------------------
1 - Resposta (integer - dichotomous)
length n NAs unique
39 39 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
0 19 48.72% 33.87% 63.80%
1 20 51.28% 36.20% 66.13%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
2 - Volume (numeric)
length n NAs unique 0s mean meanCI'
39 39 0 26 0 1.36 1.10
100.0% 0.0% 0.0% 1.62
.05 .10 .25 median .75 .90 .95
0.59 0.60 0.80 1.10 1.65 2.42 3.23
range sd vcoef mad IQR skew kurt
3.30 0.82 0.60 0.52 0.85 1.35 1.12
lowest : 0.40, 0.55, 0.60 (3), 0.70, 0.75 (3)
highest: 2.35, 2.70, 3.20, 3.50, 3.70
' 95%-CI (classic)
------------------------------------------------------------------------------
3 - Razao (numeric)
length n NAs unique 0s mean meanCI'
39 39 0 31 0 1.69 1.40
100.0% 0.0% 0.0% 1.97
.05 .10 .25 median .75 .90 .95
0.45 0.71 1.08 1.62 2.00 3.04 3.35
range sd vcoef mad IQR skew kurt
3.72 0.88 0.52 0.79 0.92 0.49 -0.27
lowest : 0.03, 0.40, 0.45, 0.57, 0.75 (3)
highest: 3.00, 3.20, 3.33, 3.50, 3.75
' 95%-CI (classic)Ajuste do modelo
modelo <- glm(Resposta~Volume+Razao, family = binomial(link="logit"), data = dados)Verificando aceitação do modelo:
Fixamos \(\alpha = 0.05\)
phi1 <- summary(modelo)$dispersion
desvio1 <- summary(modelo)$deviance/phi1
q.quadr1 <- qchisq(0.95, desvio1)Temos do resultado acima que: O modelo em análise foi aceito.
| Modelo | |
|---|---|
| Desvio/\(\phi\) | 29.7723045 |
| \(\chi^2\) | 43.4964192 |
ajuste<-summary(modelo)
knitr::kable(ajuste$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -9.529586 | 3.2331906 | -2.947425 | 0.0032043 |
| Volume | 3.882150 | 1.4286087 | 2.717434 | 0.0065790 |
| Razao | 2.649118 | 0.9142176 | 2.897689 | 0.0037592 |
\[log \left ( \frac{\pi_i}{1-\pi_i} \right ) = -9.53 + 3.882. volume_i + 2.649. razão_i\]
exp(modelo$coefficients) (Intercept) Volume Razao
7.266973e-05 4.852846e+01 1.414156e+01 Análise de Diagnósticos
Teste de Normalidade
rd <- boot::glm.diag(modelo)$rd
teste1 <- nortest::lillie.test(rd)
teste2 <- shapiro.test(rd)| Teste 1 | Teste 2 | |
|---|---|---|
| p-valor | 0.5890915 | 0.2011729 |
q1 <- ggplot2::ggplot(dados, ggplot2::aes(sample=rd))+
ggplot2::stat_qq(size=4, color="#e7ad52") + ggplot2::ylim(-4, 4) + ggplot2::xlim(-4, 4) +
ggplot2::stat_qq_line(colour = '#001132') +
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
q1
É observado através da análise gráfica e dos testes que os resíduos são distribuidos normalmente.
Ligação
va <- fitted(modelo)
g1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va,y=Resposta))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,1) + ggplot2::xlim(0,1) + ggplot2::ggtitle("Ligação")
G1 <- g1+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
G1
Variância
q1 <- ggplot2::ggplot(dados, ggplot2::aes(sample = rd))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Variância")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
q1
Pontos aberrantes
indice=as.vector(1:length(dados$Resposta))
abe <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=rd)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = -2, slope = 0, color='#001132')+ ggplot2::geom_abline(intercept = 2, slope = 0, color='#001132') + ggplot2::ylim(-4,4) + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
abe
Conforme vemos a análise gráfica, temos 2 pontos aberrantes, que são os pontos 4 e 18. Sendo as observações:
dados[which(rd>2 | rd< -2),] Resposta Volume Razao
4 1 0.75 1.500
18 1 0.85 1.415Pontos de Alavanca
h<- boot::glm.diag(modelo)$h
h3 <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=h)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = 2*(length(coef(modelo))/length(dados$Resposta)), slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
h3
Conforme o gráfico vemos que tem alguns pontos de alavanca, mais precisamente 3 pontos. Que são as observações de número:12, 13 e 32.
dados[which(h>2*length(coef(modelo))/length(dados$Resposta)),] Resposta Volume Razao
12 0 0.55 2.75
13 0 0.60 3.00
32 0 2.35 0.03Pontos de influência
DC<-boot::glm.diag(modelo)$cook
dc <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=DC)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = qchisq(0.05,length(coef(modelo)))/2, slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
dc
É observado que temos três pontosde influência, que são as observações de número: 4,18 e 32. Sendo os valores observados:
dados[which(DC>qchisq(0.05,length(coef(modelo)))/2),] Resposta Volume Razao
4 1 0.75 1.500
18 1 0.85 1.415
32 0 2.35 0.030Curva ROC
Epi::ROC(form = Resposta~Volume+Razao, data = dados, plot = "ROC", MX = TRUE)
Conforme o gráfico é observado que o ponto de corte é de 0.474
Construção da matriz Confusão
fitt <- fitted(modelo)
threshold = 0.474
y.hat <- ifelse(fitt>threshold,1,0)
mat.conf <- caret::confusionMatrix(as.factor(y.hat), as.factor(dados$Resposta), positive = "1")
knitr::kable(mat.conf$table)| 0 | 1 | |
|---|---|---|
| 0 | 17 | 3 |
| 1 | 2 | 17 |
Conforme é observado temos uma boa precisão no modelo ajustado.
Exemplo 3
Agora vamos ver a aplicação com resposta binária para analisar parte dos dados descritos no arquivo prefauto.dat sobre a preferência de consumidores americanos com relação a automóveis. Uma amostra aleatória de 263 consumidores foi considerada. As seguintes variáveis foram observadas para cada comprador: preferência do tipo de automóvel (1:americano,0:japonês), idade (em anos), sexo (0: masculino; 1: feminino) e estado civil (0:casado, 1: solteiro). Para maiores detalhes ver Foster, Stine e Waterman (1998, pgs. 338-339). Os dados estão no arquivo
A resposta, nesse exemplo, \(Y_i\) é a preferência com relação ao tipo do automóvel pelo \(i\)-ésimo comprador (1: americano e 0:japonês). O modelo seria \(Y_i \sim Be(\pi_i)\), em que:
\[log \left ( \frac{\pi_i}{1-\pi_i} \right ) = \beta_1 + \beta_2. idade_i + \beta_3. sexo_i + \beta_4.ecivil_i\]
com \(\pi_i\) denotando a probabilidade de preferência de automóvel americano.
rm(list=ls(all=TRUE))
setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/AULAS/LOGÍSTICA")
dados <- read.table("prefauto.txt", sep = "", header = TRUE)
dados$Sexo <- as.factor(dados$Sexo)
dados$EstCivil <- as.factor(dados$EstCivil)Análise Descritiva dos dados
DescTools::Desc(dados, digits = 2, plotit = F)------------------------------------------------------------------------------
Describe dados (data.frame):
data frame: 263 obs. of 4 variables
263 complete cases (100.0%)
Nr ColName Class NAs Levels
1 Pref integer .
2 Idade integer .
3 Sexo factor . (2): 1-0, 2-1
4 EstCivil factor . (2): 1-0, 2-1
------------------------------------------------------------------------------
1 - Pref (integer - dichotomous)
length n NAs unique
263 263 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
0 148 56.27% 50.23% 62.14%
1 115 43.73% 37.86% 49.77%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
2 - Idade (integer)
length n NAs unique 0s mean meanCI'
263 263 0 30 0 30.70 29.97
100.0% 0.0% 0.0% 31.44
.05 .10 .25 median .75 .90 .95
23.00 24.00 26.00 30.00 34.00 38.00 40.00
range sd vcoef mad IQR skew kurt
42.00 6.05 0.20 5.93 8.00 0.96 2.09
lowest : 18, 20, 21 (4), 22 (5), 23 (13)
highest: 45, 46 (2), 51, 54, 60
heap(?): remarkable frequency (8.7%) for the mode(s) (= 25)
' 95%-CI (classic)
------------------------------------------------------------------------------
3 - Sexo (factor - dichotomous)
length n NAs unique
263 263 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
0 144 54.75% 48.71% 60.66%
1 119 45.25% 39.34% 51.29%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
4 - EstCivil (factor - dichotomous)
length n NAs unique
263 263 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
0 170 64.64% 58.69% 70.17%
1 93 35.36% 29.83% 41.31%
' 95%-CI (Wilson)Ajuste do modelo
modelo <- glm(Pref~Idade+Sexo+EstCivil, family = binomial(link="logit"), data = dados)Verificando aceitação do modelo:
Fixamos \(\alpha = 0.05\)
phi1 <- summary(modelo)$dispersion
desvio1 <- summary(modelo)$deviance/phi1
q.quadr1 <- qchisq(0.95, desvio1)Temos do resultado acima que: O modelo em análise foi aceito.
| Modelo | |
|---|---|
| Desvio/\(\phi\) | 349.6956659 |
| \(\chi^2\) | 394.3024885 |
ajuste<-summary(modelo)
knitr::kable(ajuste$coefficients)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -1.6528404 | 0.7077086 | -2.3354815 | 0.0195183 |
| Idade | 0.0498201 | 0.0215768 | 2.3089669 | 0.0209454 |
| Sexo1 | 0.0941828 | 0.2555684 | 0.3685228 | 0.7124834 |
| EstCivil1 | -0.5178675 | 0.2725333 | -1.9001988 | 0.0574070 |
\[log \left ( \frac{\pi_i}{1-\pi_i} \right ) = -1.653 + 0.05. idade_i + 0.094. sexo_i -0.518. estcivil_i\]
exp(modelo$coefficients)(Intercept) Idade Sexo1 EstCivil1
0.1915052 1.0510820 1.0987606 0.5957897 Análise de Diagnósticos
Teste de Normalidade
rd <- boot::glm.diag(modelo)$rd
teste1 <- nortest::lillie.test(rd)
teste2 <- shapiro.test(rd)| Teste 1 | Teste 2 | |
|---|---|---|
| p-valor | 2.1951339^{-66} | 3.2303722^{-19} |
q1 <- ggplot2::ggplot(dados, ggplot2::aes(sample=rd))+
ggplot2::stat_qq(size=4, color="#e7ad52") + ggplot2::ylim(-4, 4) + ggplot2::xlim(-4, 4) +
ggplot2::stat_qq_line(colour = '#001132') +
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
q1Warning: Removed 2 rows containing missing values (`geom_path()`).
É observado através da análise gráfica e dos testes que os resíduos não são distribuidos normalmente.
Ligação
va <- fitted(modelo)
g1 <- ggplot2::ggplot(dados, ggplot2::aes(x=va,y=Pref))+
ggplot2::geom_point(size=4, color = "#e7ad52")+
ggplot2::ylim(0,1) + ggplot2::xlim(0,1) + ggplot2::ggtitle("Ligação")
G1 <- g1+ ggplot2::geom_abline(intercept = 0, slope = 1, color = "#001132")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
G1
Variância
q1 <- ggplot2::ggplot(dados, ggplot2::aes(sample = rd))+
ggplot2::stat_qq(size=4, color = "#e7ad52")+
ggplot2::ylim(-4,4) + ggplot2::xlim(-4,4) +
ggplot2::stat_qq_line(colour = "#001132") +
ggplot2::ggtitle("Variância")+
ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "#001132"),
axis.text.x = ggplot2::element_text(colour = "#001132"),
axis.title.y = ggplot2::element_text(colour = "#001132"),
axis.text.y = ggplot2::element_text(colour = "#001132"))
q1Warning: Removed 2 rows containing missing values (`geom_path()`).
Pontos aberrantes
indice=as.vector(1:length(dados$Pref))
abe <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=rd)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = -2, slope = 0, color='#001132')+ ggplot2::geom_abline(intercept = 2, slope = 0, color='#001132') + ggplot2::ylim(-4,4) + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
abe
Não foi diagnosticado nenhum ponto aberrante
Pontos de Alavanca
h<- boot::glm.diag(modelo)$h
h3 <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=h)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = 2*(length(coef(modelo))/length(dados$Pref)), slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))
h3
Conforme o gráfico vemos que tem alguns pontos de alavanca, mais precisamente 7 pontos. Que são as observações de número:99, 106, 168, 220, 223, 235 e 250.
dados[which(h>2*length(coef(modelo))/length(dados$Pref)),] Pref Idade Sexo EstCivil
99 0 60 1 1
106 1 41 1 1
168 1 45 1 1
220 1 46 0 1
223 1 54 0 1
235 1 51 1 0
250 0 46 1 1Pontos de influência
DC<-boot::glm.diag(modelo)$cook
dc <- ggplot2::ggplot(dados, ggplot2::aes(x=indice, y=DC)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = qchisq(0.05,length(coef(modelo)))/2, slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132'))+ggplot2::ylim(0,0.5)
dc
Conforme a análise gráfica não há nenhum ponto de influência
Curva ROC
Epi::ROC(form = Pref~Idade+Sexo+EstCivil, data = dados, plot = "ROC", MX = TRUE)
Conforme o gráfico é observado que o ponto de corte é de 0.448
Construção da matriz Confusão
fitt <- fitted(modelo)
threshold = 0.448
y.hat <- ifelse(fitt>threshold,1,0)
mat.conf <- caret::confusionMatrix(as.factor(y.hat), as.factor(dados$Pref), positive = "1")
knitr::kable(mat.conf$table)| 0 | 1 | |
|---|---|---|
| 0 | 93 | 48 |
| 1 | 55 | 67 |
O modelo não ficou com um bom valor preditivo. Sendo a sua acurácia baixa.