# Nomes dos pacotes
packages <- c("foreign", "nnet", "ggplot2", "reshape2", "kableExtra")

# Instala os pacotes que ainda não foram instalados
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Carrega os pacotes
invisible(lapply(packages, library, character.only = TRUE))

Introdução

A regressão logística nominal é usada para modelar variáveis categóricas nominais.

Exemplo Os alunos ingressantes no ensino médio fazem escolhas de programa entre programa geral, programa vocacional e programa acadêmico. A escolha pode ser modelada usando sua pontuação de escrita e o status socioeconômico de cada aluno.

O banco de dados utilizado no modelo é o ‘hsbdemo’.

ml  <-  read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
head(ml) %>%
  kbl() %>%
  kable_styling("hover", full_width = F)
id female ses schtyp prog read write math science socst honors awards cid
45 female low public vocation 34 35 41 29 26 not enrolled 0 1
108 male middle public general 34 33 41 36 36 not enrolled 0 1
15 male high public vocation 39 39 44 26 42 not enrolled 0 1
67 male low public vocation 37 37 42 33 32 not enrolled 0 1
153 male middle public vocation 39 31 40 39 51 not enrolled 0 1
51 female high public general 42 36 42 31 39 not enrolled 0 1

Descrição dos dados

O conjunto de dados contém variáveis referente a 200 alunos. A variável resposta é prog, que é tipo de programa. As variáveis preditoras são status socioeconômico (ses) que é uma variável categórica com três níveis e pontuação de escrita (write) que é uma variável contínua.

Abaixo são apresentadas algumas estatísticas descritivas das variáveis de interesse, como frequência para cada grupo, média e desvio padrão.

kbl(with(ml, table(ses, prog))) %>%
  kable_styling("hover", full_width = F)
general academic vocation
low 16 19 12
middle 20 44 31
high 9 42 7 
kbl(with(ml, do.call(rbind, tapply(write, prog, function(x) c(Media = round(mean(x),2), DP = round(sd(x),2))))))  %>%
  kable_styling("hover", full_width = F)
Media DP
general 51.33 9.40
academic 56.26 7.94
vocation 46.76 9.32

Modelo ajustado

Utilizando a função multinom do pacote nnet, estimaremos um modelo de regressão logística nominal.

Primeiro, escolhemos a categoria de referência do desfecho, especificando com o uso da função relevel. Em seguida, calculamos o p-valor para os testes dos coeficientes da regressão, através do teste de Wald.

ml$prog2 <- relevel(ml$prog, ref = 'academic')
test <- multinom(prog2 ~ ses + write, data = ml) # nnet::multinom(prog2 ~ ses + write, data = ml)
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.982880
## final  value 179.981726 
## converged
summary(test) # resumo dos resultados obtidos pelo modelo ajustado
## Call:
## multinom(formula = prog2 ~ ses + write, data = ml)
## 
## Coefficients:
##          (Intercept)  sesmiddle    seshigh      write
## general     2.852198 -0.5332810 -1.1628226 -0.0579287
## vocation    5.218260  0.2913859 -0.9826649 -0.1136037
## 
## Std. Errors:
##          (Intercept) sesmiddle   seshigh      write
## general     1.166441 0.4437323 0.5142196 0.02141097
## vocation    1.163552 0.4763739 0.5955665 0.02221996
## 
## Residual Deviance: 359.9635 
## AIC: 375.9635

A saída de resumo do modelo possui um bloco de coeficientes e um bloco de erros padrão. Cada um desses blocos possui uma linha de valores correspondente a uma equação do modelo.

  • Cálculo para a interpretação de \(\beta_{13}\):

\(\beta_{13} = -0.0579287 \approx -0.06\)

**Interpretação: existe uma redução aproximada de 0,06 no logaritmo da chance da escolha do programa geral, em relação a acadêmica, dado o aumento em uma unidade no score de escrita.

  • Cálculo para a interpretação de \(\beta_{21}\):

\(\beta_{21} = 0.2914 \approx 0.30\)

**Interpretação: existe um aumento de aproximadamente 0,30 no logaritmo da chance da escolha do programa vocacional, em relação a acadêmico, dado a mudança do status social baixo para médio.

Teste de Wald

z <- summary(test)$coefficients/summary(test)$standard.errors # estatística do valor-z
kbl(round(z,1), caption = "Estatística do Teste de Wald") %>%
  kable_styling("hover", full_width = F)
Estatística do Teste de Wald
(Intercept) sesmiddle seshigh write
general 2.4 -1.2 -2.3 -2.7
vocation 4.5 0.6 -1.6 -5.1

O valor-z é uma estatística de teste para testes que medem a razão entre o coeficiente e seu erro padrão. Ele é utilizado para calcular o valor-p, auxiliando na tomada de uma decisão sobre a significância estatística dos termos.

P-valor

p <- (1 - pnorm(abs(z),0,1)) * 2 # teste de Wald bilateral
kbl(round(p,4)) %>%
  kable_styling("hover", full_width = F)
(Intercept) sesmiddle seshigh write
general 0.0145 0.2294 0.0237 0.0068
vocation 0.0000 0.5408 0.0989 0.0000

O p-valor é uma medida de probabilidade que mede a evidência contra a hipótese nula. As estimativas mais próximas a zero (valor p menor que 0,05, supondo 5% de significância) fornecem evidências mais fortes contra a hipótese nula.

Razão de chances (RC)

As razões de chances são as medidas de associação produzidas por modelos de regressão logística. Neste modelo elas são obtidas exponenciando os coeficientes de regressão.

kbl(round(exp(coef(test)),4)) %>%
  kable_styling("hover", full_width = F)
(Intercept) sesmiddle seshigh write
general 17.3258 0.5867 0.3126 0.9437
vocation 184.6126 1.3383 0.3743 0.8926
  • Cálculo para a interpretação da \(RC= 0.9437\):

\(0.9437 \approx 0.94\)

\(100\% - 94\% = 6\%\)

**Interpretação: o aumento em uma unidade na variável write acarreta uma redução aproximada de 6% na chance de se escolher o programa geral ao invés do acadêmico.

  • Cálculo para a interpretação da \(RC= 0.3126\):

\(0.3126 \approx 0.31\)

\(100\% - 31\% = 69\%\)

**Interpretação: mudar o status da classe social baixa para alta reduz aproximadamente 69% a chance de se escolher o programa geral ao invés do acadêmico.

  • Cálculo para a interpretação da \(RC= 1.3383\):

\(1.3383 \approx 1.34\)

\(1.34 - 1 = 0.34\)

**Interpretação: mudar o status da classe social baixa para média aumenta aproximadamente 34% a chance de se escolher o programa vocacional ao invés do acadêmico (esta razão de chances é não significativa).

IC95% para as Razões de Chances

kbl(round(exp(confint(test)),2)) %>%
  kable_styling("hover", full_width = F)
2.5 %.general 97.5 %.general 2.5 %.vocation 97.5 %.vocation
(Intercept) 1.76 170.44 18.87 1805.84
sesmiddle 0.25 1.40 0.53 3.40
seshigh 0.11 0.86 0.12 1.20
write 0.90 0.98 0.85 0.93

Estes intervalos de confiança (IC95%) permitem ao pesquisador conhecer a precisão das estimativas das razões de chances.

Ajuste das Probabilidades Preditas

É possível estimar as probabilidades preditas para os sujeitos da amostra de acordo com seu status social e seu score de escrita, usando a função fitted.

kbl(round(head(pp <- fitted(test)),4)) %>%
  kable_styling("hover", full_width = F)
academic general vocation
0.1483 0.3382 0.5135
0.1202 0.1806 0.6992
0.4187 0.2368 0.3445
0.1727 0.3508 0.4765
0.1001 0.1689 0.7309
0.3534 0.2378 0.4088

Foram criados pequenos conjuntos de dados, variando uma variável enquanto mantemos a outra constante. Primeiro mantendo fixa a média de write e examinando as probabilidades preditas para cada nível de ses.

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
kbl(round(predict(test, newdata = dses, "probs"),4)) %>%
  kable_styling("hover", full_width = F)
academic general vocation
0.4397 0.3582 0.2021
0.4777 0.2283 0.2939
0.7009 0.1785 0.1206

Outra maneira de entender o modelo usando as probabilidades preditas é olhar para as probabilidades médias preditas para os diferentes valores da variável preditora contínua write dentro de cada nível de ses.

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
    3))

## armazena as probabilidades preditas para cada valor de ses e write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))

## calcula as probabilidades médias dentro de cada nível de ses
by(pp.write[, 3:5], pp.write$ses, colMeans)
## pp.write$ses: high
##  academic   general  vocation 
## 0.6164315 0.1808037 0.2027648 
## ------------------------------------------------------------ 
## pp.write$ses: low
##  academic   general  vocation 
## 0.3972977 0.3278174 0.2748849 
## ------------------------------------------------------------ 
## pp.write$ses: middle
##  academic   general  vocation 
## 0.4256198 0.2010864 0.3732938

Análise gráfica

lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")

ggplot(lpp, aes(x = write, y = probability, colour = ses)) + 
  geom_line() + 
  facet_grid(variable ~ ., scales = "free") + 
  theme_minimal()

Com a análise gráfica acima podemos evidenciar que o aluno que possui uma classe social mais alta está relacionado ao programa acadêmico. Além disso, na escolha do programa vocacional ao acadêmico, o status social baixo muda para alto.