Regressao Logistica

Exemplos - Simulacao

sim_logistic_data = function(sample_size = 25, beta_0 = -2, beta_1 = 3) {
  x = rnorm(n = sample_size)
  eta = beta_0 + beta_1 * x
  p = 1 / (1 + exp(-eta))
  y = rbinom(n = sample_size, size = 1, prob = p)
  data.frame(y, x)
}

vamos simular dados do seguinte modelo:

\[ \log\left(\frac{p({\bf x})}{1 - p({\bf x})}\right) = -2 + 3 x \]

Outra maneira de escrevê-lo que melhor corresponde à função que estamos usando para simular os dados:

\[ Y_i \sim \text{Bern}(p_i) \\ p_i = p({\bf x_i}) = \frac{1}{1 + e^{-\eta({\bf x_i})}} \\ \eta({\bf x_i}) = -2 + 3 x_i \]

set.seed(1)

set.seed(1)
example_data = sim_logistic_data()
head(example_data)
##   y          x
## 1 0 -0.6264538
## 2 1  0.1836433
## 3 0 -0.8356286
## 4 1  1.5952808
## 5 0  0.3295078
## 6 0 -0.8204684

Após a simulacao, observa-se que a variável resposta “Y” usa apenas valores 0 e 1. Logo, confirma-se que a variável explicativa “X” é binária. Facamos dois tipos de regressao, linear simples.

# Regressao linear ordinaria
fit_lm  = lm(y ~ x, data = example_data)

# Regressao logistica
fit_glm = glm(y ~ x, data = example_data, family = binomial)

Observe que a sintaxe utilizada nos modelos propostos sao quase identicas. A diferenca do primeiro para o segundo modelo, esta no lm que passou a ser glm. O lm é uma versao mais especifica do glm. O glm utiliza como argumento padrao “familia gaussiana”, mas neste caso utilizamos “familia binomial”, como argumento.

# Modelo generalizado, familia gaussiana

glm(y ~ x, data = example_data)
## 
## Call:  glm(formula = y ~ x, data = example_data)
## 
## Coefficients:
## (Intercept)            x  
##      0.3090       0.3021  
## 
## Degrees of Freedom: 24 Total (i.e. Null);  23 Residual
## Null Deviance:       5.76 
## Residual Deviance: 3.782     AIC: 29.73

O argumento family para glm() especifica a distribuição e a função de link.

# Chamada mais detalhada para a regressao logistica - modelo generalizado familia  binomial

fit_glm = glm(y ~ x, data = example_data, family = binomial(link = "logit"))

Fazer previsões com um objeto do tipo ‘glm’ é um pouco diferente do que fazer previsões após o ajuste com lm (). No caso da regressão logística, com family = binomial, o type =“link” irá obter os log odds, enquanto type=“response” retornará \(P[Y=1\mid{\bf X}={\bf x}]\) para cada observação.

Com relacao ao gráfico abaixo, vale salientar que nao se trata de estimativa de Y. Estamos plotando a média estimada, pois, como a variável resposta Y recebe valores 0 e 1, resposta binária, a média estimada é uma estimativa de \(P[Y=1\mid{\bf X}={\bf x}]\). Logo, se a média é uma probabilidade, queremos probabilidades que estejam entre o intervalo [0,1].

plot(y ~ x, data = example_data, 
     pch = 20, ylab = "Estimated Probability", 
     main = "Ordinary vs Logistic Regression")
abline(fit_lm, col = "darkorange")
curve(predict(fit_glm, data.frame(x), type = "response"), 
      add = TRUE, col = "dodgerblue", lty = 2)
legend("topleft", c("Ordinary", "Logistic", "Data"), lty = c(1, 2, 0), 
       pch = c(NA, NA, 20), lwd = 2, col = c("darkorange", "dodgerblue", "black"))

Nota-se que a média está estimando valores menores que 0. Contudo a regressao linear ordinaria nao é adequada. É possivel ver que os dados não se ajustam a reta de regressao ajustada isso porque a mesma não é limitada ao intervalo [0,1], mas vai de [0,\(\infty\)]. Já na regressão logística, os dados parecem se ajustar melhor a curva de regressão logística estimada.

Como a saída da função logit inversa é restrita a estar entre 0 e 1, nossas estimativas fazem muito mais sentido como probabilidades. Observemos os coeficientes estimados.

round(coef(fit_glm), 1)
## (Intercept)           x 
##        -2.3         3.7

Logo o modelo estimado é:

\[ \log\left(\frac{\hat{p}({\bf x})}{1 - \hat{p}({\bf x})}\right) = -2.3 + 3.7 x \]

Devemos ter cautela na interpretação de \((\hat{\beta}_1 = 3.7)\), já que não estamos estimando diretamente a média, mas sim uma função da média. A interpretação adequada é:

  • Para um aumento de uma unidade em x, as odds de log mudam (nesse caso, aumentam) em 3,7. Além disso, como \((\hat{\beta}_1 )\) é positivo, à medida que aumentamos x, também aumentamos \(\hat{p}({\bf x})\).

Partindo para outro exemplo, consideremos o modelo abaixo.

\[ \log\left(\frac{p({\bf x})}{1 - p({\bf x})}\right) = 1 + -4 x. \]

podemos reescrevê-lo para melhor correspondencia à função que estamos usando para simular os dados:

\[ Y_i \sim \text{Bern}(p_i) \\ p_i = p({\bf x_i}) = \frac{1}{1 + e^{-\eta({\bf x_i})}} \\ \eta({\bf x_i}) = 1 + -4 x_i \]

Nesse modelo, à medida que x aumenta, as log-verossimilhanças diminuem.

set.seed(1)
example_data = sim_logistic_data(sample_size = 50, beta_0 = 1, beta_1 = -4)

Novamente simulamos algumas observações deste modelo, depois ajustamos a regressão logística.

fit_glm = glm(y ~ x, data = example_data, family = binomial)
plot(y ~ x, data = example_data, 
     pch = 20, ylab = "Estimated Probability", 
     main = "Logistic Regression, Decreasing Probability")
curve(predict(fit_glm, data.frame(x), type = "response"), 
      add = TRUE, col = "dodgerblue", lty = 2)
legend("bottomleft", c("Estimated Probability", "Data"), lty = c(2, 0), 
       pch = c(NA, 20), lwd = 2, col = c("dodgerblue", "black"))

Note que, desta vez, quando x aumenta, \(\hat{p}({\bf x})\) diminui.

Agora vamos ver um exemplo em que a probabilidade estimada nem sempre aumenta ou diminui simplesmente. Assim como a regressão linear comum, a combinação linear de preditores pode conter transformações de preditores (neste caso, um termo quadrático) e interações.

sim_quadratic_logistic_data = function(sample_size = 25) {
  x = rnorm(n = sample_size)
  eta = -1.5 + x ^ 2
  p = 1 / (1 + exp(-eta))
  y = rbinom(n = sample_size, size = 1, prob = p)
  data.frame(y, x)
}

\[ \log\left(\frac{p({\bf x})}{1 - p({\bf x})}\right) = -1.5 + x^2. \]

Novamente, podemos reescrever para melhor corresponder à função que estamos usando para simular os dados:

\[ Y_i \sim \text{Bern}(p_i) \\ p_i = p({\bf x_i}) = \frac{1}{1 + e^{-\eta({\bf x_i})}} \\ \eta({\bf x_i}) = -1.5 + x_i^2 \]

set.seed(42)
example_data = sim_quadratic_logistic_data(sample_size = 50)
fit_glm = glm(y ~ x + I(x^2), data = example_data, family = binomial)
plot(y ~ x, data = example_data, 
     pch = 20, ylab = "Estimated Probability", 
     main = "Logistic Regression, Quadratic Relationship")
curve(predict(fit_glm, data.frame(x), type = "response"), 
      add = TRUE, col = "dodgerblue", lty = 2)
legend("bottomleft", c("Prob", "Data"), lty = c(2, 0), 
       pch = c(NA, 20), lwd = 2, col = c("dodgerblue", "black"))

Exemplo com dados reais (MichelinFood)

# Carregando os dados

MichelinFood <- read.table("MichelinFood.txt", header=TRUE)
attach(MichelinFood)
plot(Food,proportion,ylab="Sample proportion",xlab="Zagat Food Rating")

Figura 1. Gráfico da proporção da amostra de “sucessos” em relação às classificações de alimentos.

Observa-se que os dados plotados no gráfico tem formato de “S”. Com isso sabe-se que a variável “Y” recebe como resposta valores 0 e 1, sabe-se também que a variável explicativa “X” é binária e que utilizaremos a regressao logistica para análise do modelo.

# Modelo
m1 <- glm(cbind(InMichelin,NotInMichelin)~Food,family=binomial)

Os odds a favor do “sucesso” são definidas como a razão da probabilidade de que “sucesso” ocorrerá, para a probabilidade de que “sucesso” não ocorrerá. Seja x a classificação de alimentos Zagat para um determinado restaurante francês e \(\theta\)(x) denota a probabilidade de que este restaurante esteja incluído no guia Michelin. Então nosso modelo de regressão logística para a resposta, \(\theta\) baseado na variável preditora x é dado por:

\(\theta\) = 1/ [1 + exp(- {\(\hat{\beta}_0\) + \(\hat{\beta}_1\)X})]

summary(m1)
## 
## Call:
## glm(formula = cbind(InMichelin, NotInMichelin) ~ Food, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4850  -0.7987  -0.1679   0.5913   1.5889  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -10.84154    1.86236  -5.821 5.84e-09 ***
## Food          0.50124    0.08768   5.717 1.08e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 61.427  on 13  degrees of freedom
## Residual deviance: 11.368  on 12  degrees of freedom
## AIC: 41.491
## 
## Number of Fisher Scoring iterations: 4

O rearranjo da equação do modelo ajustado fornece o log(odds) ou logit

log[\(\hat{\theta}\)(x) / (1 - \(\hat{\theta}\)(x))] = {\(\hat{\beta}_0\) + \(\hat{\beta}_1\)X} = -10.842 + 0.501x

Observe que o log(odds) ou logit é uma função linear de x. Os odds estimados sendo incluídos no guia Michelin são dadas por:

(\(\hat{\theta}\)(x) / 1 - \(\hat{\theta}\)(x)) = exp(- {\(\hat{\beta}_0\) + \(\hat{\beta}_1\)X}) = exp(-10.842 + 0.501x)

Temos que \((\hat{\beta}_0 = -10.84154)\) e que \((\hat{\beta}_1 = 0.50124)\). Nesse caso pode-se interpretar a estimativa dos betas da seguinte forma:

  • Para um aumento de uma unidade em x, as odds de log aumentam em 0,5. E pelo fato de \((\hat{\beta}_1 )\) ser positivo, à medida que aumentamos x, aumentamos a funcao \(\hat{\theta}({\bf x})\).
#Figura 2 
x <- seq(15,28,0.05)
y <- 1/(1+exp(-1*(m1$coeff[1] + m1$coeff[2]*x)))
plot(Food,proportion,ylab="Probability of inclusion in the Michelin Guide",xlab="Zagat Food Rating")
lines(x,y)

Figura 2. Gráfico da proporção da amostra de “sucessos” em relação às classificações de alimentos.

O gráfico acima mostra o modelo de regressao logistica ajustado nos dados da figura 1 e é marcado com uma curva suave. Por exemplo, se x, a classificação de alimentos de Zagat é aumentada uma unidade, as chances de ser incluído no guia Michelin aumentam exp (0,501) = 1,7 e se é aumentada em cinco unidades, as chances de ser incluído no guia Michelin aumentam exp (5 × 0,501) = 12,2.

#Table 2
thetahat <- m1$fitted.values
odds_ratio <- m1$fitted.values/(1-m1$fitted.values)
cbind(Food,round(thetahat,3),round(odds_ratio,3))
##    Food             
## 1    15 0.035  0.036
## 2    16 0.056  0.060
## 3    17 0.089  0.098
## 4    18 0.140  0.162
## 5    19 0.211  0.268
## 6    20 0.306  0.442
## 7    21 0.422  0.729
## 8    22 0.546  1.204
## 9    23 0.665  1.988
## 10   24 0.766  3.281
## 11   25 0.844  5.416
## 12   26 0.899  8.941
## 13   27 0.937 14.759
## 14   28 0.961 24.364

Tabela 2. Probabilidades estimadas e odds obtidas a partir do modelo logístico.

A Tabela 2 apresenta as probabilidades estimadas e as probabilidades obtidas a partir do modelo logístico. Tomando a proporção de entradas sucessivas na última coluna da Tabela 2. Exemplo:

0.060 / 0.036 = 0.098 / 0.060 =. = 24.364 / 14.759 = 1.7) reproduz o resultado que o aumento de x (Zagat food rating) por uma unidade aumenta as chances de ser incluído no guia Michelin por 1.7.

Observe na Tabela 2 que os odds são maiores que 1 quando a probabilidade é maior que 0.5. Nestas circunstâncias, a probabilidade de “sucesso” é maior que a probabilidade de “fracasso”.