VD categórica com mais de uma categoria

Uma extensão do modelo binomial se torna necessária quando temos mais do que duas categorias. No caso do plebiscito, por exemplo, a variável dependente original trazia a categoria dos indecisos. Repare que as três categorias, “sí”, “no” e “indeciso”, não possui uma ordem entre elas, o que caracateriza um modelo logit multinomial. Neste modelo calculamos qual a probabilidade do voto do eleitor \(i\) cair em alguma das categorias. Seja \(\pi_{ij}\) a probabilidade da \(iª\) observação cair na \(jª\) categoria, isto é, \(\pi_{ij} \equiv \text{Pr}(Y_{i} = {j})\) para \(j = 1,\dots,m\) categorias:

\[ \begin{align} \pi_{ij} &= \frac{\text{exp}(\gamma_{0j}+\gamma_{1j}X_{i1}+\dots+\gamma_{kj}X_{ik})}{ 1 + \sum_{l=1}^{m-1}\text{exp}(\gamma_{0j}+\gamma_{1j}X_{i1}+\dots+\gamma_{kj}X_{ik})} \text{ para } j = 1,\dots,m-1\\ \pi_{im} &= 1 - \sum_{l=1}^{m-1}\pi_{ij} \end{align} \]

Com alguma manipulação algébrica temos

\[ \text{log}\frac{\pi_{ij}}{\pi_{im}} = \gamma_{0j}+\gamma_{1j}X_{i1}+\dots+\gamma_{kj}X_{ik} \text{ para } j = 1,\dots,m-1\] Vemos que a última categoria se torna um tipo de base, e \(\pi_{im}\) entra no denominador do cálculo da probabilidade nas outras categorias. O coeficiente, portanto, nos dá o log da chance de se pertencer a categoria \(j\) contra a categoria base \(m\) [como fica a equação acima no caso de termos apenas duas categorias?]

O modelo multinomial é uma extensão da binomial onde em cada categoria temos um sequência de bernoullis. A distrbuição de probabilidade é dada por

\[ \begin{align} Pr(Y = y_i\mid \pi) &= \prod_{j}^{m} \prod_{i}^{n} \pi_{ij}^{w_{ij}} \\ \text{log}Pr(Y = y_i\mid \pi) &= \sum_{j}^{m} \sum_{i}^{n} {w_{ij} \text{log}(\pi_{ij}}) \\ \end{align} \] Onde \(w_{ij} = 1\) quando \(Y_i = j\) e\(w_{ij} = 0\) quando \(Y_i \neq j\) Substituindo a expressão para \(\pi\) na euqação acima obtmeos a verossimilhança e resolvendo numericamente obtemos os parâmetros que a maximiza.

[Exemplo Fox e Weisberg 1]

Muitas vezes temos mais do que duas categorias exaustivas e mutuamente exclusivas. Avaliação de governo em geral envolvem as categoras “Ótimo”, “Bom”, “Regular”, “Ruim” e “Péssimo”. Neste caso temos 5 categorias exaustivas e mutuamente exclusivas e que, ainda, mantém uma ordem entre si. Podemos pensar que essas cinco categorias são partições de um a’variável latente’ \(z\). Essa variável, não observada (ex: ‘conformismo’), é função linear de nossas variáveis independentes. Conforme \(z\) vai assumindo diferentes valores muda-se a categoria de \(Y\). Podemos estimar pontos latentes a partir do qual ocorre a transição entre categorias:

\[ \begin{align} \nonumber y_{i} = \left\{ \begin{array}{l l} 1 \text{ se } z_{i} < c_{1.5}\\ 2 \text{ se } z_{i} \in(c_{1.5},c_{2.5})\\ 3 \text{ se } z_{i} \in(c_{2.5},c_{3.5})\\ 4 \text{ se } z_{i} \in(c_{3.5},c_{4.5})\\ 5 \text{ se } z_{i} > c_{4.5}\\ \end{array} \right. \\ \\ z_{i} \sim \text{logistic}(x_{i}, \sigma^2) \end{align} \]

Onde \(c_{(.)}\) é o valor dos pontos de transição na escala de \(x\) e \(\sigma\) nos dá o quão forte é a transição. A distribuição de probabilidade acumulada de \(y\) é dada por:

\[ \begin{align} \text{Pr}(y \leq j \mid \boldsymbol{x}) &= \text{Pr}(z \leq c_{j} \mid \boldsymbol{x}) \\ &= \text{Pr}(\alpha + \beta_{1}x_{1} + \dots + \beta_{k}x_{k} + \epsilon \leq c_{j} \mid \boldsymbol{x}) \\ &= \text{Pr}(\epsilon \leq c_{j} - \alpha - \beta_{1}x_{1} - \dots - \beta_{k}x_{k} \mid \boldsymbol{x}) \end{align} \]

Para \(j = 1,2 \dots, q-1\), isto é, a quebra de \(z\) em \(q\) categorias acrescenta \(q-1\) parâmetros ao modelo. Assumindo que \(\epsilon\) segue uma distribuição logística então temos o modelo logit ordinal:

\[ \begin{align} \text{logit}[\text{Pr}(y >j\mid \boldsymbol{x})] &= \text{log} \frac{\text{Pr}(y >j\mid \boldsymbol{x})}{\text{Pr}(y \leq j\mid \boldsymbol{x})} \\ &= (\alpha - c_{j}) - \beta_{1}x_{1} - \dots - \beta_{k}x_{k} \end{align} \]

Repare que, na equação acima, a única difernça entre os diferentes \(j\)’s está no intercepto \((\alpha - c_{j})\), portanto as curvas em um gráfico desta função serão paralelas para cada \(j\). Isso quer dizer que os logits, ou as chances, são proporcionais. Por isso alguns chama esse modelo de modelo de chance proporcinal ou, em inglês, proportional-odds logistic model. Ao utilizar um programa para estima-lo ele retornará os coeficientes de \(x_{i}\) e os pontos de transição \(c_{j}\).

[Exemplo Fox e Weisberg 2]

Alternativamente podemos estimar esse modelo com logits separados por categoria, se forem poucas, ou por MQO, se forem muitas. O problema do MQO é que ele considera os intervalos entre as categorias como sendo do mesmo tamanho, o que, por se tratar de variáveis qualitativas nem sempre é verdade.

Modelos não paramétricos de classificação

Os modelos que vimos até agora foram estimados utilizando métodos paramétricos. Nós estabelecemos uma forma funcional para nosso modelo que tenta se aproximar da verdadeira forma da distribuição \(f(Y \mid X)\) e então estimamos parâmetros, isto é, valores que sintetizam a relação entre \(Y\) e \(X\). Se a forma funcional que estabelecemos para nosso modelo for uma distribuição normal, por exemplo, estimamos a média e variância e temos assim uma noção de como \(Y\) varia com \(X\), do mesmo modo se nossa distribuição for da família da binomial estimando \(\pi\) obtemos uma decrição de \(f(Y\mid X)\).

Um possíve problema do método paramétrico é que o modelo que usamos nem sempre reflete corretamente a verdadeira forma de \(f(Y \mid X)\). Um outro método para conhecermos a distribuição de \(f(Y\mid X)\) não parte de uma forma funcional, mas tenta apreender esta forma diretamente dos dados. Estes métodos tentam se aproximar o máximo possível de todos os pontos dos dados evitando, ao mesmo tempo, variar muito.

Um exemplo de método não paramétrico é, como vimos, o uso de lowess para estimar o voto no plebiscito no Chile. Existem vários outros métodos não paramétricos ou semi-paramétricos para VD’s tanto quantitativas quanto qualitativas. Um deles é o K-Nearest Neighbors. Por este método atribuimos a observação \(i\) à categoria A ou B com uma probabilidade dada pela quantidade de vizinhos que pertencem a uma ou outra categoria.

\[ \text{Pr}(Y = j\mid X = x) = \frac{1}{K} \sum_{i \in \mathcal{N(x)}} I(y_i = j)\]

Onde \(j\) é uma das categorias e \(\mathcal{N}\) é a vizinhança de \(i\). O KNN atribui \(x\) à categoria \(j\) com a maior probabilidade de ocorrência entre os vizinhos.

Observem que o KNN é uma abordagem completamente não paramétrica e não nos diz quais variáveis são importantes, não obtemos nada parecido com uma tabela de coeficientes como vimos quando usamos outros modelos, isto é, um modelo pouco interpretável. Esse tipo de modelo é usado quando estamos interessados mais em predição do que em inferência. Neste último caso estamos interessados na relação entre a variévale dependente e as variáveis explicativas, isto é, de que forma um determinado preditor, como a adesão ao status quo no caso do plebiscito, impacta na VD, o voto no sim ou não. Observamos os coeficientes para testar nossas teorias sobre essas relações, nossos modelos teóricos. Na predição estamos interessados em estimar futuros eventos com base nos dados e algumas suposições sobre como aqueles dados surgiram. No KNN, por exemplo, com base no conhecimento que adquirirmos ao observar dados da relação entre \(Y\) e \(X\), conseguims predizer, com algum erro, as categorias a que pertencerão um novo \(\tilde X\). Vamos falar mais de predição na próxima aula.

Um método menos radical de classificação é o Latent Discriminant Analysis (LDA). Neste método supomos que uma variável independente se distribui de determinada forma em cada categoria da variável dependente, isto é, os valores que \(X\) assume é condicinal à categoria que \(Y\) assume, ou \(f_{k}(X) = Pr(X = x\mid Y = k)\). Sendo \(\pi_{k}\) uma probabilidade a priori de que uma observação aleatória venha da \(kª\) categoria da variável \(Y\), então, pelo teorema de Bayes temos:

\[ \text{Pr}(Y = k \mid X = x) = \frac{\pi_{k}f_{k}(X)}{\sum_{l=1}^{K}\pi_{l}f_{l}(X)}\] Isso sugere que no lugar de calcular \(\text{Pr}(Y = k \mid X)\) como fizemos quando estimamos o modelo logit, podemos apenas estimar \(\pi_{k}\) e \(f_{k}(X)\), substituir na equação e classificar na categoria \(k\) a observação que tiver a maior \(\text{Pr}(Y = k \mid X)\). Para estimar \(f_{k}(X)\) assume-se que sua distribuição é normal. Com apenas um preditor a densidade da normal tem a forma

\[ f_{k}(X) = \frac{1}{\sqrt{2 \pi \sigma_{k}}} \text{exp} \left(-\frac{1}{\sigma_{k}^{2}}(x - \mu_{k})^2) \right)\]

substituindo na equação anterior temos

\[ \text{Pr}(Y = k \mid X = x) = \frac{\pi_{k}\frac{1}{\sqrt{2 \pi \sigma}} \text{exp}\left(-\frac{1}{\sigma^{2}}(x - \mu_{k})^2\right) }{\sum_{l=1}^{K}\pi_{l}\frac{1}{\sqrt{2 \pi \sigma}} \text{exp} \left(-\frac{1}{\sigma^{2}}(x - \mu_{l})^2) \right)}\]

Tirando o log e rearranjando temos

\[ \delta_{k}(x) = x \cdot \frac{\mu_{k}}{\sigma^2} - \frac{\mu_{k}^2}{2\sigma^2} + log(\pi_{k})\]

O método do Linear Discriminant Analysis (LDA) estima \(\pi_{k}\), \(\mu_{k}\) e \(\sigma^2\) a partir das observações:

\[ \begin{align} \tilde\mu_{k} &= \frac{1}{n_{k}} \sum_{i:y_{i} = k} x_{i} \\ \sigma^2 &= \frac{1}{n + K} \sum_{k = 1}^{k} \sum_{i:y_{i}=k}(x_i - \tilde\mu_{k}^2) \\ \pi_{k} &= \frac{n_{k}}{n} \end{align} \]

O LDA assume, além da normalidade da distribuição de \(X\) por categoria, que a variância \(\sigma^2\) é constante. Quando relaxamos essa premissa temos o método da Quadratic Discriminant Analysis (QDA). Esses métodos não partem de um modelo a priori da relação entre \(Y\) e \(X\). A partir das observações, calculamos \(\tilde\mu_{k}\), \(\sigma^2\) e \(\pi_{k}\)e com eles encontramos a função discriminante \(\delta_{k}(x)\). Classificamos \(x\) na categoria em que o discrimante apresentar o maior valor..

Exemplo (http://www.lac.inpe.br/~rafael.santos/Docs/CAP394/WholeStory-Iris.html):

O LDA, ou o QDA, tem a vantagem de não pressupor um modelo e, ao mesmo tempo, ser um pouco mais interpretável do que o KNN.

Predição

Ao ajustar um modelo estamos interessados em quão bem os valores preditos correspondem aos valores observados. Calculamos a qualidade do ajuste por meio do Erro Médio Quadrado, no caso de VD quantitativa ou regressão, e da Taxa de Erro para VDs qualitativas ou classificação:

\[ \text{EMQ} = \frac{1}{n} \sum_{i=1}^{n}(y_{i} - \hat{f}(x_{i}))^2\]

e

\[ \text{TE} = \frac{1}{n} \sum_{i=1}^{n} I(y_{i} \neq \hat{y_{i}})\] Como vimos acima podemos tornar nossos modelos flexíveis o bastante para se ajustar quase que perfeitamente aos dados. O problema é que, quando trabalhamos com predição, não nos interessa tanto que nosso modelo se ajuste às observações nas quais ele foi treinado, mas sim a observações ainda não vistas. Não nos interessa, por exemplo, descrever os fatores que levaram a determinados golpes de estado, mas sim conseguir prever se haverá um golpe ou não. Também não nos interessa apenas descrever os determinantes do desempenho passado de alunos, mas sim qual será o desempenho conforme variações nestes determinantes.

Se tivermos dados o suficiente podemos usar nossa estimativa \(\hat f\) para calacular o MSE ou o TE de observações teste e escolher o método que nos dá o menor erro no teste.

[exemplo viés e variância]

Observe que, a partir de certa quantidade de flexibilidade, o MSE do teste começa a subir apesar do MSE do treino diminuir continuamente. Quando o MSE treino é pequeno, mas o MSE teste é grande temos o que se chama de overfitting. A partir de certo ponto nosso modelo começa a capturar padrões nos dados que ocorrem ao acaso e não tem nada a ver com o processo ‘verdadeiro’ que gerou aqueles dados. Quando ajustamos esse modelo ao teste aqueles padrões aleatórios não estão presentes e o ajuste não é tão bom. Esse padrão de “U” no MSE de teste ocorre sempre em qualquer base de dados e em qualquer modelo.

Portanto devemos procurar ajustar nosso modelo a novos dados para evitar overffiting e garantir validade externa a eles. No entanto, muitas vezes não temos dados o suficiente para separar uma base de teste. Quando isso ocorre podemos usar um conjunto de técnicas que permitem extrair várias bases teste de um mesmo conjunto de dados. Essas técnicas recebem o nome de cross validation ou validação cruzada.

Na validação simples fazemos como no exemplo que acabamos de ver no R: sorteamos aleatoriamente algumas observações de nossos dados e usamos para treino, o resto será a base teste. Neste proceso acabamos com aenas uma base de teste. Na validação cruzada tentamos maximizar o uso dos dados para validação. No limite temos a Leave-One-Out Cross-Validation (LOOCV) onde separamos uma observação para teste e usamos todas as outras para treinar o modelo, repetindo o procedimento para todas as observações. Esse processo exige muita computação. Um métod mais eficiente é o K-Fold Cross-Validation.

No k-fold CV dividimos, aleatoriamente, nossa base em k grupos (folds). Um grupo servirá como teste e os outros para treino, repetindo para todos os grupos. Para cada teste obtemos um MSE.Tiramos a média de todos esses k’s MSE e obtemos o erro quadrado médio do k-fold CV.

No R:

library(ISLR)
data("Auto")
library(boot)
data(Auto)
set.seed(17)
cv.erro.10 <- rep(0,10)
for(i in 1:10){
  glm.fit <- glm(mpg~poly(horsepower, i), data = Auto)
  cv.erro.10[i] <- cv.glm(Auto, glm.fit, K = 10)$delta[1]
}
cv.erro.10
##  [1] 24.27207 19.26909 19.34805 19.29496 19.03198 18.89781 19.12061 19.14666
##  [9] 18.87013 20.95520