Regressão Logística
by Edneide Ramalho
30 Abril, 2021
1 O modelo
A estimação dos coeficientes da regressão logística ao contrário da regressão múltipla que utiliza o método dos mínimos quadrados, é efetuada pelo uso da máxima verossimilhança.
A qualidade do ajuste do modelo é avaliada pelo “pseudo” \(R^2\) e pelo exame da precisão preditiva (matriz de confusão).
# Bibliotecas
library(readr)
library(caret)
library(pROC)
library(mfx)
library(ResourceSelection)
library(modEvA)
library(foreign)
library(stargazer)
library(dplyr)2 Regressão Logística Simples
Neste exemplo vamos utilizar somente uma variável independente, neste caso numérica.
A variável dependente é a ocorrência ou não (1 ou 0) de doença coronária cardíaca (CHD), associando-se com a idade (AGE) dos indivíduos, criando assim um modelo de regressão logística.
require(readr)
chd <- read_delim("https://github.com/Smolski/livroavancado/raw/master/cdh.csv",
";", escape_double = FALSE, col_types = cols(CHD = col_factor(levels = c())),
trim_ws = TRUE)
summary(chd)## AGE AGRP CHD
## Min. :20.00 Min. :1.00 0:57
## 1st Qu.:34.75 1st Qu.:2.75 1:43
## Median :44.00 Median :4.00
## Mean :44.38 Mean :4.48
## 3rd Qu.:55.00 3rd Qu.:7.00
## Max. :69.00 Max. :8.00require(ggplot2)
ggplot(chd, aes(x=AGE, y=CHD)) +
geom_point() +
stat_smooth(method="glm", method.args=list(family="binomial"), se=FALSE)## `geom_smooth()` using formula 'y ~ x'## Warning: Computation failed in `stat_smooth()`:
## y values must be 0 <= y <= 1
Pelo gráfico, parece que uma relação entre CHD e idade.
Construindo o modelo:
m1 = glm(CHD ~ AGE, family = binomial(link = "logit"),
data = chd)
summary(m1)##
## Call:
## glm(formula = CHD ~ AGE, family = binomial(link = "logit"), data = chd)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9718 -0.8456 -0.4576 0.8253 2.2859
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.30945 1.13365 -4.683 2.82e-06 ***
## AGE 0.11092 0.02406 4.610 4.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 136.66 on 99 degrees of freedom
## Residual deviance: 107.35 on 98 degrees of freedom
## AIC: 111.35
##
## Number of Fisher Scoring iterations: 4- Função de ligação estimada do modelo:
\[ \ln\left(\frac{prob_{CHD}}{1 - prob_{CHD}}\right) = -5.309 + 0.1109 \times AGE \]
O coeficiente de 0,1109 para idade, pelo fato de ser positivo informa que quando a idade (AGE) se eleva, elevam-se as chances de ocorrência de CHD (p < 0.001).
Vamos usar o modelo para prever as probabilidades de CHD para todas as idades da base de dados.
# Filtrnado a idade dos indivíduos
idade <- chd[, 1]# Criando predição
chd$PRED <- predict(m1, newdata = idade, type = "response")# Plotando a probabilidade predita pelo modelo
require(ggplot2)
ggplot(chd, aes(x = AGE, y = PRED)) +
geom_point()
2.1 Estimando a Razão de Chances (OR)
- O log das chances da variável idade no modelo é 0,1109.
require(mfx)
logitor(CHD ~ AGE, data = chd)## Call:
## logitor(formula = CHD ~ AGE, data = chd)
##
## Odds Ratio:
## OddsRatio Std. Err. z P>|z|
## AGE 1.117307 0.026882 4.6102 4.022e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1razão de chances da variável AGE foi de 1.1173. Interpretação: para cada variação unitária na idade (AGE), as chances de ocorrência de CHD aumentam 1,1173 vezes.
Dito de outra forma, para cada variação unitária em AGE, aumentam-se 11,73% ((1,1173-1)*100) as chances da ocorrência de CHD.
2.2 Determinando o Intervalo de Confiança
exp(cbind(OR=coef(m1), confint(m1)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 0.004944629 0.0004412621 0.0389236
## AGE 1.117306795 1.0692223156 1.17586812.3 Predição de Probabilidades
- Usamos a função
predict()para calcular a probabilidade de CHD para uma certa idade.
media = data.frame(AGE = mean(chd$AGE))
media## AGE
## 1 44.38media$pred.prob = predict(m1, newdata = media, type = "response")
media## AGE pred.prob
## 1 44.38 0.4044944- A probabilidade de uma pessoa de 44 anos ter CHD é de 40.45%.
2.4 Matriz de confusão
- Primeiro precisa-se determinar um ponto de corte para classificar em 0 ou 1.
Precisão: representa a proporção das predições corretas do modelo sobre o total.
\[ \text{ACC} = \frac{VP + VN}{P + N} \]
P representa o ttotal de eventos positivos (Y = 1) e N é o total de “não-eventos” (Y = 0, ou negativo)
Sensibilidade: representa a proporção de verdadeiros positivos, ou seja, a capacidade do modelo em avaliar o evento como \(\hat{Y} = 1\) (estimado) dado que ele é evento real Y = 1.
\[ \text{SENS} = \frac{VP}{FN} \]
Especificidade: a proporção apresentada dos verdadeiros negativos, ou seja, o poder de predição do modelo em avaliar como “não evento” \(\hat{Y} = 0\) sendo que ele não é evento Y=0:
\[ ESPEC = \frac{VN}{VN + VP} \]
Verdadeiro Preditivo Positivo: se caracteriza como proporção de verdadeiros positivos com relação ao total de predições positivas, ou seja, se o evento é real Y=1 dada a classificação do modelo \(\hat{Y}=1\).
\[ VPP = \frac{VP}{VN + FP} \]
Verdadeiro Preditivo Negativo: se caracteriza pela proporção de verdadeiros negativos comparando-se com o total de predições negativas, ou seja, o indivíduo não ser evento Y=0Y=0 dada classificação do modelo como “não evento” ^Y=0Y^=0:
\[ VPN = \frac{VN}{VN + FN} \]
require(caret)
chd$pdata <- as.factor(
ifelse(
predict(m1,
newdata = chd,
type = "response")
> 0.5, "1", "0"
)
)confusionMatrix(chd$pdata, chd$CHD, positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 45 14
## 1 12 29
##
## Accuracy : 0.74
## 95% CI : (0.6427, 0.8226)
## No Information Rate : 0.57
## P-Value [Acc > NIR] : 0.0003187
##
## Kappa : 0.4666
##
## Mcnemar's Test P-Value : 0.8445193
##
## Sensitivity : 0.6744
## Specificity : 0.7895
## Pos Pred Value : 0.7073
## Neg Pred Value : 0.7627
## Prevalence : 0.4300
## Detection Rate : 0.2900
## Detection Prevalence : 0.4100
## Balanced Accuracy : 0.7319
##
## 'Positive' Class : 1
## 2.5 Curva ROC
- A curva ROC é produzida bi-dimensionalmente pela obtenção da relação entre a taxa dos verdadeiros positivos do modelo e da taxa dos falsos positivos preditos.
- Hosmer e Lemeschow (2000) sugere a utilização de AUC acima de 0,7 como aceitável.
require(pROC)
roc1 = plot.roc(chd$CHD, fitted(m1))## Setting levels: control = 0, case = 1## Setting direction: controls < cases
plot(roc1,
print.auc=TRUE,
auc.polygon=TRUE,
grud=c(0.1,0.2),
grid.col=c("green","red"),
max.auc.polygon=TRUE,
auc.polygon.col="lightgreen",
print.thres=TRUE)
2.6 O teste de Hosmer e Lemeshow
O teste de Hosmer e Lemeshow é utilizado para demonstrar a qualidade do ajuste do modelo, ou seja, se o modelo pode explicar os dados observados.
require(ResourceSelection) hl=hoslem.test(chd$CHD,fitted(m1),g=10)## Warning in Ops.factor(1, y): '-' not meaningful for factorshl## ## Hosmer and Lemeshow goodness of fit (GOF) test ## ## data: chd$CHD, fitted(m1) ## X-squared = 100, df = 8, p-value < 2.2e-16A hipótese nula \(H_0\) do qui-quadrado (
p=0,05) deste teste é a de que as proporções observadas e esperadas são as mesmas ao longo da amostra.O teste rejeita a hipótese nula.
2.7 Pseudo \(R^2\)
require(modEvA)
RsqGLM(m1)## $CoxSnell
## [1] 0.2540516
##
## $Nagelkerke
## [1] 0.3409928
##
## $McFadden
## [1] 0.2144684
##
## $Tjur
## [1] 0.2705749
##
## $sqPearson
## [1] 0.27255183 Regressão Logística Múltipla
require(haven)## Loading required package: havenmy_data <- read_dta("http://dss.princeton.edu/training/Panel101.dta")
glimpse(my_data)## Rows: 70
## Columns: 9
## $ country <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ year <dbl> 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 19…
## $ y <dbl> 1342787840, -1899660544, -11234363, 2645775360, 3008334848, 32…
## $ y_bin <dbl> 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,…
## $ x1 <dbl> 0.27790365, 0.32068470, 0.36346573, 0.24614404, 0.42462304, 0.…
## $ x2 <dbl> -1.1079559, -0.9487200, -0.7894840, -0.8855330, -0.7297683, -0…
## $ x3 <dbl> 0.28255358, 0.49253848, 0.70252335, -0.09439092, 0.94613063, 1…
## $ opinion <dbl+lbl> 1, 3, 3, 3, 3, 1, 3, 1, 3, 4, 2, 1, 2, 4, 3, 4, 1, 3, 2, 4…
## $ op <dbl> 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1,…- Vamos criari um modelo de regressão logística em que nossa variável resposta é y, e nossas variáveis independentes são $x_1, x_2, x_3$.
logit = glm(y_bin ~ x1 + x2 + x3, data = my_data, family = binomial(link = "logit"))
summary(logit)##
## Call:
## glm(formula = y_bin ~ x1 + x2 + x3, family = binomial(link = "logit"),
## data = my_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0277 0.2347 0.5542 0.7016 1.0839
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4262 0.6390 0.667 0.5048
## x1 0.8618 0.7840 1.099 0.2717
## x2 0.3665 0.3082 1.189 0.2343
## x3 0.7512 0.4548 1.652 0.0986 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.056 on 69 degrees of freedom
## Residual deviance: 65.512 on 66 degrees of freedom
## AIC: 73.512
##
## Number of Fisher Scoring iterations: 5Quando x3 se eleva de uma unidade, o log das chances altera-se em 0.7512.
As 3 variáveis independentes possuem efeitos positivos para a determinação das chances do preditor ser igual a 1.
A variável x3 revelou significância estatística a 10%, no entando o valor ususal para considerá-la significante é de 5%.
require(stargazer)
stargazer(logit, title = "Resultados", type = "text")##
## Resultados
## =============================================
## Dependent variable:
## ---------------------------
## y_bin
## ---------------------------------------------
## x1 0.862
## (0.784)
##
## x2 0.367
## (0.308)
##
## x3 0.751*
## (0.455)
##
## Constant 0.426
## (0.639)
##
## ---------------------------------------------
## Observations 70
## Log Likelihood -32.756
## Akaike Inf. Crit. 73.512
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01- Transformações necessárias para encontrar o OR:
require(mfx)
logitor(y_bin ~ x1 + x2 + x3, data = my_data)## Call:
## logitor(formula = y_bin ~ x1 + x2 + x3, data = my_data)
##
## Odds Ratio:
## OddsRatio Std. Err. z P>|z|
## x1 2.36735 1.85600 1.0992 0.27168
## x2 1.44273 0.44459 1.1894 0.23427
## x3 2.11957 0.96405 1.6516 0.09861 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- O resultado acima evidencia que para uma alteração em 1 (uma) unidade em x3, a chance de que y seja igual a 1 aumenta em 112% ((2,12-1)*100). Dito de outra forma, a chance de y=1 é 2,12 vezes maior quando x3 aumenta em uma unidade (sendo que aqui mantêm-se as demais variáveis independentes constantes).
exp(coef(logit))## (Intercept) x1 x2 x3
## 1.531417 2.367352 1.442727 2.119566- Desta forma, para cada incremento unitário em x2 e mantendo as demais variáveis constantes, conclui-se que é 1,443 vezes provável que y seja igual a 1 em oposição a não ser (igual a zero), ou seja, as chances aumentam em 44,30%.
Intervalo de confinça
exp(cbind(OR = coef(logit), confint(logit)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 1.531417 0.4387468 5.625299
## x1 2.367352 0.5129380 11.674641
## x2 1.442727 0.8041221 2.737965
## x3 2.119566 1.0038973 5.718637Predições
allmean = data.frame(x1 = mean(my_data$x1),
x2 = mean(my_data$x2),
x3 = mean(my_data$x3))
allmean## x1 x2 x3
## 1 0.6480006 0.1338694 0.761851Predição do modelo
allmean$pred.prob = predict(logit, newdata = allmean,
type = "response")
allmean## x1 x2 x3 pred.prob
## 1 0.6480006 0.1338694 0.761851 0.83285554 Método Stepwise
direções “both”, “backward”, “forward”
AIC - Akaike Infromation Criterion
Quanto menor o AIC, melhor o ajuste
\[ \text{AIC} = -2\log(L_p) + 2[(p + 1) + 1] \]
- \(L_p\) é a função da máxima verossimilhança e p é o número de variáveis explicativas do modelo.
step(logit, direction = 'both')## Start: AIC=73.51
## y_bin ~ x1 + x2 + x3
##
## Df Deviance AIC
## - x1 1 66.736 72.736
## - x2 1 66.996 72.996
## <none> 65.512 73.512
## - x3 1 69.402 75.402
##
## Step: AIC=72.74
## y_bin ~ x2 + x3
##
## Df Deviance AIC
## - x2 1 67.330 71.330
## <none> 66.736 72.736
## + x1 1 65.512 73.512
## - x3 1 70.032 74.032
##
## Step: AIC=71.33
## y_bin ~ x3
##
## Df Deviance AIC
## <none> 67.330 71.330
## - x3 1 70.056 72.056
## + x2 1 66.736 72.736
## + x1 1 66.996 72.996##
## Call: glm(formula = y_bin ~ x3, family = binomial(link = "logit"),
## data = my_data)
##
## Coefficients:
## (Intercept) x3
## 1.1339 0.4866
##
## Degrees of Freedom: 69 Total (i.e. Null); 68 Residual
## Null Deviance: 70.06
## Residual Deviance: 67.33 AIC: 71.334.1 VIF - Variance Inflation Factor
Multicolinearidade
VIF é um in’dice que não deve ficar abaixo de 10 para representar baixo problema de multicolinearidade.
require(faraway)## Loading required package: faraway##
## Attaching package: 'faraway'## The following object is masked from 'package:lattice':
##
## melanomavif(logit)## x1 x2 x3
## 9.291872 12.317808 29.859879- x2 e x3 têm VIF maiores que 10.
5 Regressão Logística Múltipla com variável categórica
Queremos descobrir a probabildade de um aluno ser admitido em uma universidade com base
no rank da escola em que estudou no ensino médio
com as notas nas provas (GPA - Grade Point Average), e
GRE (GRaduate Record Examinations)
# Carregando o arquivo
library(readr)
binary <- read_csv("http://www.karlin.mff.cuni.cz/~pesta/prednasky/NMFM404/Data/binary.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## admit = col_double(),
## gre = col_double(),
## gpa = col_double(),
## rank = col_double()
## )glimpse(binary)## Rows: 400
## Columns: 4
## $ admit <dbl> 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1…
## $ gre <dbl> 380, 660, 800, 640, 520, 760, 560, 400, 540, 700, 800, 440, 760,…
## $ gpa <dbl> 3.61, 3.67, 4.00, 3.19, 2.93, 3.00, 2.98, 3.08, 3.39, 3.92, 4.00…
## $ rank <dbl> 3, 3, 1, 4, 4, 2, 1, 2, 3, 2, 4, 1, 1, 2, 1, 3, 4, 3, 2, 1, 3, 2…- Transformando a variável rank em categórica:
binary$rank <- factor(binary$rank)- Determinando a regressão:
mylogit <- glm(admit ~ gre + gpa + rank, data = binary,
family = binomial(link = "logit"))summary(mylogit)##
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "logit"),
## data = binary)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6268 -0.8662 -0.6388 1.1490 2.0790
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.989979 1.139951 -3.500 0.000465 ***
## gre 0.002264 0.001094 2.070 0.038465 *
## gpa 0.804038 0.331819 2.423 0.015388 *
## rank2 -0.675443 0.316490 -2.134 0.032829 *
## rank3 -1.340204 0.345306 -3.881 0.000104 ***
## rank4 -1.551464 0.417832 -3.713 0.000205 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.52 on 394 degrees of freedom
## AIC: 470.52
##
## Number of Fisher Scoring iterations: 4Utilizando o teste de análise de variância pode-se obsevar quea variância para o modelo nulo é 500.
Quando as variáveis explicativas foram incluídas, estas reduziram a variânica do modelo para 459, contribuindo para o modelo.
anova(mylogit, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: admit
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 399 499.98
## gre 1 13.9204 398 486.06 0.0001907 ***
## gpa 1 5.7122 397 480.34 0.0168478 *
## rank 3 21.8265 394 458.52 7.088e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- A função update permite criar uma nova regressão com base na regressão anterior. Vamos usar a mesma regressão excluindo a variável gre por exemplo.
mylogit2 = update(mylogit, ~. - gre)anova(mylogit, mylogit2, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: admit ~ gre + gpa + rank
## Model 2: admit ~ gpa + rank
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 394 458.52
## 2 395 462.88 -1 -4.3578 0.03684 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A variância para este novo modelo foi de 463, piorando o mesmo em relação ao modelo anterior.
Voltando para o modelo anterior:
\[ \ln(P = Y) = -3.98998 + 0.00226 \times gre + 0.80404 * gpa - 0.67544 \times rank2 - 1.34020 \times rank3 - 1.55146 \times rank4 \]
- OR para cada variável:
exp(cbind(OR = coef(mylogit), confint(mylogit)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 0.0185001 0.001889165 0.1665354
## gre 1.0022670 1.000137602 1.0044457
## gpa 2.2345448 1.173858216 4.3238349
## rank2 0.5089310 0.272289674 0.9448343
## rank3 0.2617923 0.131641717 0.5115181
## rank4 0.2119375 0.090715546 0.4706961Interpretando os resultados:
Para variáveis numéricas, a cada 1 unidade aumentada gera-se uma elevação de chances na probabilidade de ser admitido.
Para a variável gre, para cada aumento unitário, mantendo-se as outras variáveis constantes, elevem-se em 0.23% ((OR - 1) * 100) (1.023 - 1) * 100, as chances de que o aluno seja aprovado.
Já para variáveis teóricas como o rank da escola de ensino médio, diminuem-se as chances em 49,11% ((0.5089-1)*100) de que o aluno seja aprovado se ele for de uma escola com ranking 2, e em 73.82% se ele for de uma escola de ranking 3. Já, se ele for de uma escola de ranking 4, suas chances diminuem em 78.81%.
Predição de probabilidade
Qual a probabilidade de um aluno que tirou gre de 700 e 3.67 no gpa e estudou numa escola de ranking 1 ser admitido?
pred = data.frame(gre = 700,
gpa = 3.67,
rank = factor(1))pred$prob = predict(mylogit, newdata = pred, type = "response")
pred## gre gpa rank prob
## 1 700 3.67 1 0.6331924Comparando com um aluno que tirou as mesmas notas no gre e gpa, mas que estudou em uma escola ranking 4:
pred = data.frame(gre = 700, gpa = 3.67, rank = factor(4)) pred$prob2 = predict(mylogit, newdata=pred, type = "response") pred## gre gpa rank prob2 ## 1 700 3.67 4 0.2678562- comparando os ranks das escolas:
novosdados = with(binary,
data.frame(gre = mean(gre),
gpa = mean(gpa),
rank = factor(1:4)))novosdados = cbind(novosdados, predict(mylogit,
newdata = novosdados, type = "response", se.fit = TRUE))
novosdados## gre gpa rank fit se.fit residual.scale
## 1 587.7 3.3899 1 0.5166016 0.06631526 1
## 2 587.7 3.3899 2 0.3522846 0.03978481 1
## 3 587.7 3.3899 3 0.2186120 0.03825061 1
## 4 587.7 3.3899 4 0.1846684 0.04863624 1# Renomeando as variáveis
names(novosdados)[names(novosdados)=='fit'] = 'prob'
names(novosdados)[names(novosdados) == 'se.fit'] = 'se.prob'novosdados## gre gpa rank prob se.prob residual.scale
## 1 587.7 3.3899 1 0.5166016 0.06631526 1
## 2 587.7 3.3899 2 0.3522846 0.03978481 1
## 3 587.7 3.3899 3 0.2186120 0.03825061 1
## 4 587.7 3.3899 4 0.1846684 0.04863624 1# Estimando os intervalos de confiança
novosdados$LL=novosdados$prob-1.96*novosdados$se.prob
novosdados$UL=novosdados$prob+1.96*novosdados$se.probrequire(ggplot2)
ggplot(novosdados, aes(x=rank,y=prob))+
geom_errorbar(aes(ymin=LL, ymax=UL), width=0.2,lty=1,lwd=1,col="red")+
geom_point(shape=18, size=5, fill="black")+
scale_x_discrete(limits=c("1","2","3","4"))+
labs(title="Probabilidades preditas", x="Ranking",y="Pr(y=1)")