Capítulo 11: Análise Discriminante Linear (LDA) em R

DES4033: Análise Multivariada em Biologia

Bastão de Asclépio & Distribuição Normal

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))

# suppressMessages(library(MVLM, warn.conflicts=FALSE))
suppressMessages(library(apaTables, warn.conflicts=FALSE))
suppressMessages(library(aplpack, warn.conflicts=FALSE))
suppressMessages(library(bestNormalize, warn.conflicts=FALSE))
suppressMessages(library(calculus, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(cellWise, warn.conflicts=FALSE))
suppressMessages(library(cluster, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(diptest, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(EGAnet, warn.conflicts=FALSE))
suppressMessages(library(ellipse, warn.conflicts=FALSE))
suppressMessages(library(far, warn.conflicts=FALSE))
suppressMessages(library(fmsb, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(ggfortify, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(ggrepel, warn.conflicts=FALSE))
suppressMessages(library(HSAUR2, warn.conflicts=FALSE))
suppressMessages(library(irrCAC, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(latex2exp, warn.conflicts=FALSE))
suppressMessages(library(lava, warn.conflicts=FALSE))
suppressMessages(library(lavaan, warn.conflicts=FALSE))
suppressMessages(library(lawstat, warn.conflicts=FALSE))
suppressMessages(library(lfda, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(MatchIt, warn.conflicts=FALSE)) # distancia robusta
suppressMessages(library(matlib, warn.conflicts=FALSE))
suppressMessages(library(matrixcalc, warn.conflicts=FALSE))
suppressMessages(library(MatrixModels, warn.conflicts=FALSE))
suppressMessages(library(matrixStats, warn.conflicts=FALSE))
suppressMessages(library(MomTrunc, warn.conflicts=FALSE))
suppressMessages(library(MVA, warn.conflicts=FALSE)) # An Introduction to Applied Multivariate Analysis with R
suppressMessages(library(MVar.pt, warn.conflicts=FALSE))
suppressMessages(library(mvdalab, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(mvtnorm, warn.conflicts=FALSE))
suppressMessages(library(performance, warn.conflicts=FALSE))
suppressMessages(library(plotly, warn.conflicts=FALSE))
suppressMessages(library(polycor, warn.conflicts=FALSE))
suppressMessages(library(ppcor, warn.conflicts=FALSE))
suppressMessages(library(pracma, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(qcc, warn.conflicts=FALSE))
suppressMessages(library(rcompanion, warn.conflicts=FALSE))
suppressMessages(library(reticulate, warn.conflicts=FALSE))
suppressMessages(library(rgl, warn.conflicts=FALSE))
suppressMessages(library(scatterplot3d, warn.conflicts=FALSE))
source("agr2x2_G.R")

Incluir

• biotools::D2.disc: Discriminant Analysis Based on Mahalanobis Distance

Material

• HTML de Rmarkdown em RPubs
• Livro-texto:
  • JOHNSON, RA & WICHERN, DW (2007) Applied multivariate statistical analysis. 6^th ed. NJ: Prentice Hall.
  • JOHNSON, RA & WICHERN, DW (2007) Applied multivariate statistical analysis: Solutions Manual. 6^th ed. NJ: Prentice Hall.
• Arquivos de dados do livro-texto:
  • JW6Data
• Prof. José Siqueira: ResearchGate

Pensamento

“Não é paradoxo dizer que nos nossos momentos de inspiração mais teórica podemos estar o mais próximo possível de nossas aplicações mais práticas.”

WHITEHEAD, AN apud BOYER, CB (1974) História da matemática. São Paulo: Blücher/EDUSP, p. 419.

Sumário

1. Aspectos de análise multivariada
2. Álgebra matricial e vetor estocástico
3. Geometria amostral e amostragem aleatória
4. Distribuição normal multivariada
5. Inferência sobre vetor de média
6. Comparação de várias médias multivariadas
7. Modelo de regressão linear multivariada
8. Componentes principais
9. Análise de fatores e inferência sobre matriz de covariância
10. Análise de correlação canônica
11. Discriminação e classificação
12. Clusterização, métodos de distância e ordination

O que é Análise Discriminante Linear?

Análise Discriminante Linear (LDA) não tem hipótese nula. Portanto, LDA é análise descritiva exploratória multivarida.

Introdução

Discriminação e classificação são técnicas multivariadas relacionadas à separação de conjuntos distintos de objetos (ou observações) e à alocação de novos objetos (observações) em grupos previamente definidos. A análise discriminante é de natureza mais exploratória. Como procedimento de separação, é frequentemente empregada de forma pontual para investigar diferenças observadas quando as relações causais não são bem compreendidas. Os procedimentos de classificação são menos exploratórios, no sentido de que levam a regras bem definidas, que podem ser usadas para atribuir novos objetos. A classificação normalmente requer mais estrutura do problema do que a discriminação.

Assim, os objetivos imediatos da discriminação e da classificação, respectivamente, são os seguintes:

Objetivo 1. Descrever, seja graficamente (em três ou menos dimensões) ou algebricamente, as características diferenciais de objetos (observações) provenientes de várias coleções conhecidas (populações). Tentamos encontrar “discriminantes” cujos valores numéricos sejam tais que as coleções sejam separadas o máximo possível.

Objetivo 2. Classificar objetos (observações) em duas ou mais classes rotuladas. A ênfase está em derivar uma regra que possa ser usada para atribuir de forma ótima novos objetos às classes rotuladas.

Seguiremos a convenção e usaremos o termo discriminação para nos referirmos ao Objetivo 1. Essa terminologia foi introduzida por R.A. Fisher [9] no primeiro tratamento moderno de problemas de separação. Um termo mais descritivo para esse objetivo, contudo, é separação. Referir-nos-emos ao segundo objetivo como classificação ou alocação.

Uma função que separa objetos pode, às vezes, servir como um alocador e, inversamente, uma regra que aloca objetos pode sugerir um procedimento discriminatório. Na prática, os Objetivos 1 e 2 frequentemente se sobrepõem, e a distinção entre separação e alocação torna-se obscurecida.

Separação e classificação para duas populações

Para fixar as ideias, listemos situações em que se possa ter interesse em:

1. separar duas classes de objetos ou
   
2. atribuir um novo objeto a uma de duas classes (ou ambos).

É conveniente rotular as classes como \(\pi_1\) e \(\pi_2\). Os objetos são, em geral, separados ou classificados com base em medidas de, por exemplo, \(p\) variáveis aleatórias associadas \(\mathbf{X}^{\prime} = [X_1, X_2, \ldots, X_p]\).

Os valores observados de \(\mathbf{X}\) diferem, até certo ponto, de uma classe para outra. Se os valores de \(\mathbf{X}\) não fossem muito diferentes para objetos em \(\pi_1\) e \(\pi_2\), não haveria problema; isto é, as classes seriam indistinguíveis e novos objetos poderiam ser atribuídos indistintamente a qualquer uma das classes.

Podemos pensar o conjunto total de valores da primeira classe como sendo a população de valores de \(\mathbf{x}\) para \(\pi_1\) e os da segunda classe como a população de valores de \(\mathbf{x}\) para \(\pi_2\).

Essas duas populações podem então ser descritas por funções densidade de probabilidade \(f_1(\mathbf{x})\) e \(f_2(\mathbf{x})\) e, consequentemente, podemos falar em atribuir observações a populações ou objetos a classes de forma intercambiável.

Você talvez se lembre de que alguns dos exemplos das situações de separação-classificação a seguir foram introduzidos no Capítulo 1.

Exemplos de duas populações e variáveis medidas:

Populações \(\pi_1\) e \(\pi_2\)	Variáveis medidas \(\mathbf{X}\)
1. Companhias de seguro de responsabilidade patrimonial solventes e em dificuldade.	Ativos totais, custo de ações e títulos, valor de mercado de ações e títulos, despesas com sinistros, superávit, montante de prêmios emitidos.
2. Dispépticos não ulcerosos (aqueles com distúrbios gástricos) e controles (“normais”).	Medidas de ansiedade, dependência, culpa, perfeccionismo.
3. Federalist Papers escritos por James Madison e os escritos por Alexander Hamilton.	Frequências de diferentes palavras e comprimentos das sentenças.
4. Duas espécies de erva-de-passarinho (chickweed).	Comprimento do sépalo e da pétala, profundidade da fenda da pétala, comprimento da bráctea, comprimento do ápice escarioso, diâmetro do pólen.
5. Compradores de um novo produto e retardatários (os que demoram a comprar).	Escolaridade, renda, tamanho da família, quantidade de trocas prévias de marca.
6. Estudantes universitários bem-sucedidos ou malsucedidos (não se formam).	Notas no exame de admissão, média do ensino médio, número de atividades no ensino médio.
7. Homens e mulheres.	Medidas antropológicas, como circunferência e volume em crânios antigos.
8. Bons e maus riscos de crédito.	Renda, idade, número de cartões de crédito, tamanho da família.
9. Alcoólatras e não alcoólatras.	Atividade da enzima monoamina oxidase, atividade da enzima adenilato ciclase.

Vemos pelo item 5, por exemplo, que objetos (consumidores) devem ser separados em duas classes rotuladas (“compradores” e “retardatários”) com base nos valores observados de variáveis presumivelmente relevantes (escolaridade, renda e assim por diante).

Na terminologia de observação e população, queremos identificar uma observação da forma

\[ \mathbf{x}^{\prime} = [x_1(\text{escolaridade}),\; x_2(\text{renda}),\\ x_3(\text{tamanho da família}),\; x_4(\text{quantidade de troca de marca})] \]

como pertencente à população \(\pi_1\), compradores, ou à população \(\pi_2\), retardatários.

Neste ponto, concentraremos a atenção na classificação para duas populações, retornando à separação na Seção 11.3.

Regras de alocação ou classificação são, em geral, desenvolvidas a partir de amostras de “aprendizado”. Características medidas de objetos selecionados aleatoriamente, conhecidos como provenientes de cada uma das duas populações, são examinadas quanto a diferenças.

Essencialmente, o conjunto de todos os resultados amostrais possíveis é dividido em duas regiões, \(R_1\) e \(R_2\), de tal forma que, se uma nova observação cair em \(R_1\), ela é alocada à população \(\pi_1\), e, se cair em \(R_2\), é alocada à população \(\pi_2\).

Assim, um conjunto de valores observados favorece \(\pi_1\), enquanto o outro conjunto de valores favorece \(\pi_2\).

Você pode se perguntar, neste ponto, como sabemos que algumas observações pertencem a uma determinada população, mas estamos incertos sobre outras. (Isto, claro, é o que torna a classificação um problema!) Várias condições podem dar origem a essa aparente anomalia (ver [19]):

1. Conhecimento incompleto do desempenho futuro.

   Exemplos: No passado, valores extremos de certas variáveis financeiras foram observados 2 anos antes da falência subsequente de uma empresa. Classificar outra empresa como sólida ou em dificuldade com base nos valores observados desses indicadores antecipados pode permitir que os administradores tomem medidas corretivas, se necessário, antes que seja tarde demais.

   Um escritório de admissões de uma faculdade de medicina pode querer classificar um candidato como provável de tornar-se médico ou improvável de tornar-se médico com base em notas de testes e outros registros universitários. Aqui, a determinação real só pode ser feita ao final de vários anos de treinamento.

2. Informação “perfeita” requer destruir o objeto.

   Exemplo: A vida útil de uma bateria de calculadora é determinada usando-a até que ela falhe, e a resistência de uma peça de madeira é obtida carregando-a até que ela quebre. Produtos que falham não podem ser vendidos. Gostaríamos de classificar produtos como bons ou ruins (não atendendo às especificações) com base em certas medições preliminares.

3. Informação indisponível ou cara.

   Exemplos: Supõe-se que certos dos Federalist Papers foram escritos por James Madison ou Alexander Hamilton porque foram assinados por eles. Outros dos Papers, contudo, não foram assinados e é de interesse determinar qual dos dois homens escreveu os Papers não assinados. Claramente, não podemos perguntar a eles. Frequências de palavras e comprimentos de sentenças podem ajudar a classificar os Papers disputados.

   Muitos problemas médicos podem ser identificados conclusivamente apenas realizando uma operação cara. Geralmente, desejamos diagnosticar uma enfermidade a partir de sintomas externos facilmente observáveis, porém potencialmente falíveis. Essa abordagem ajuda a evitar operações desnecessárias — e caras.

\[\Diamond\]

Deve ficar claro, a partir desses exemplos, que regras de classificação normalmente não podem fornecer um método de atribuição livre de erros. Isso ocorre porque pode não haver uma distinção clara entre as características medidas das populações; isto é, os grupos podem se sobrepor.

Nesse caso, é possível, por exemplo, classificar incorretamente um objeto de \(\pi_2\) como pertencente a \(\pi_1\) ou um objeto de \(\pi_1\) como pertencente a \(\pi_2\).

Exemplo 11.1: Discriminando proprietários e não proprietários de cortadores de grama

Considere dois grupos em uma cidade: \(\pi_1\), proprietários de cortadores de grama, e \(\pi_2\), aqueles sem cortadores de grama — isto é, não proprietários. Para identificar os melhores clientes potenciais para uma campanha intensiva de vendas, um fabricante de cortadores de grama está interessado em classificar famílias como potenciais proprietários ou não proprietários com base em \(x_1 =\) renda e \(x_2 =\) tamanho do lote.

Amostras aleatórias de \(n_1 = 12\) proprietários atuais e \(n_2 = 12\) não proprietários atuais fornecem os valores mostrados na Tabela 11.1.

Tabela 11.1:

\(\pi_1\): Proprietários de cortador		\(\pi_2\): Não proprietários
\(x_1\) (Renda em milhares de dólares)	\(x_2\) (Área do lote em 1000 ft\(^2\))	\(x_1\) (Renda em milhares de dólares)	\(x_2\) (Área do lote em 1000 ft\(^2\))
60.0	18.4	75.0	19.6
85.5	16.8	52.8	20.8
64.8	21.6	64.8	17.2
61.5	20.8	43.2	20.4
87.0	23.6	84.0	17.6
110.1	19.2	49.2	17.6
108.0	17.6	59.4	16.0
82.8	22.4	66.0	18.4
69.0	20.0	47.4	16.4
93.0	20.8	33.0	18.8
51.0	22.0	51.0	14.0
81.0	20.0	63.0	14.8

Esses dados estão representados na Figura 11.1. Observamos que os proprietários de cortadores de grama tendem a ter maiores rendas e lotes maiores do que os não proprietários, embora a renda pareça ser um “discriminante” melhor do que o tamanho do lote. Por outro lado, há alguma sobreposição entre os dois grupos.

Se, por exemplo, alocássemos os valores de \((x_1, x_2)\) que caem na região \(R_1\) (como determinada pela linha contínua na figura) à \(\pi_1\), proprietários, e aqueles valores de \((x_1, x_2)\) que caem na região \(R_2\) à \(\pi_2\), não proprietários, cometeríamos alguns erros.

Alguns proprietários de cortadores seriam classificados incorretamente como não proprietários e, inversamente, alguns não proprietários como proprietários.

A ideia é criar uma regra (regiões \(R_1\) e \(R_2\)) que minimize as chances de cometer esses erros. (Ver Exercício 11.2.)

\[\Diamond\]

Um bom procedimento de classificação deve resultar em poucas classificações incorretas. Em outras palavras, as chances, ou probabilidades, de erro de classificação devem ser pequenas. Como veremos, existem características adicionais que uma regra de classificação “ótima” deve possuir.

Pode ocorrer que uma classe ou população tenha maior probabilidade de ocorrência do que outra, porque uma das duas populações é relativamente muito maior do que a outra. Por exemplo, tende a haver mais empresas financeiramente sólidas do que empresas falidas. Como outro exemplo, uma espécie de erva-de-passarinho pode ser mais prevalente do que outra.

Uma regra de classificação ótima deve levar em conta essas “probabilidades a priori de ocorrência”. Se realmente acreditarmos que a probabilidade (a priori) de uma empresa financeiramente em dificuldade e, no final, falida é muito pequena, então deve-se classificar uma empresa selecionada aleatoriamente como não falida, a menos que os dados favoreçam esmagadoramente a falência.

Outro aspecto da classificação é o custo. Suponha que classificar um objeto de \(\pi_1\) como pertencente a \(\pi_2\) represente um erro mais grave do que classificar um objeto de \(\pi_2\) como pertencente a \(\pi_1\). Então, deve-se ter cautela ao fazer a primeira atribuição. Como exemplo, deixar de diagnosticar uma doença potencialmente fatal é substancialmente mais “custoso” do que concluir que a doença está presente quando, na verdade, não está. Um procedimento de classificação ótimo deve, sempre que possível, levar em conta os custos associados à classificação incorreta.

Sejam \(f_1(\mathbf{x})\) e \(f_2(\mathbf{x})\) as funções densidade de probabilidade associadas ao vetor aleatório \(\mathbf{X}\) de dimensão \(p \times 1\) para as populações \(\pi_1\) e \(\pi_2\), respectivamente. Um objeto com medições associadas \(\mathbf{x}\) deve ser atribuído a \(\pi_1\) ou \(\pi_2\).

Seja \(\Omega\) o espaço amostral — isto é, o conjunto de todos os possíveis valores observados de \(\mathbf{x}\). Seja \(R_1\) o conjunto de valores de \(\mathbf{x}\) para os quais classificamos os objetos como pertencentes a \(\pi_1\) e \(R_2 = \Omega - R_1\) os valores restantes de \(\mathbf{x}\) para os quais classificamos os objetos como pertencentes a \(\pi_2\).

Como todo objeto deve ser atribuído a uma e somente uma das duas populações, os conjuntos \(R_1\) e \(R_2\) são mutuamente exclusivos e exaustivos.

Para \(p = 2\), podemos ter um caso como o ilustrado na Figura 11.2.

A probabilidade condicional \(P(2 \mid 1)\) de classificar um objeto como pertencente a \(\pi_2\) quando, na verdade, ele é de \(\pi_1\) é

\[ P(2 \mid 1) = P(\mathbf{X} \in R_2 \mid \pi_1) = \int_{R_2 = \Omega - R_1} f_1(\mathbf{x})\, d\mathbf{x} \tag{11-1} \]

Analogamente, a probabilidade condicional \(P(1 \mid 2)\) de classificar um objeto como pertencente a \(\pi_1\) quando ele é realmente de \(\pi_2\) é

\[ P(1 \mid 2) = P(\mathbf{X} \in R_1 \mid \pi_2) = \int_{R_1} f_2(\mathbf{x})\, d\mathbf{x} \tag{11-2} \]

O sinal de integral em (11-1) representa o volume formado pela função densidade \(f_1(\mathbf{x})\) sobre a região \(R_2\). Analogamente, o sinal de integral em (11-2) representa o volume formado por \(f_2(\mathbf{x})\) sobre a região \(R_1\). Isso é ilustrado na Figura 11.3 para o caso univariado, \(p = 1\).

Seja \(p_1\) a probabilidade a priori de \(\pi_1\) e \(p_2\) a probabilidade a priori de \(\pi_2\), onde \(p_1 + p_2 = 1\). Então, as probabilidades totais de classificar corretamente ou incorretamente objetos podem ser obtidas como o produto das probabilidades a priori pelas probabilidades condicionais de classificação:

\[ P(\text{observação corretamente classificada como } \pi_1) = \\ P(\text{observação vem de } \pi_1 \text{ e é corretamente classificada como } \pi_1) =\\ P(\mathbf{X} \in R_1 \mid \pi_1) P(\pi_1) = \\ P(1 \mid 1)\, p_1\\ \\ P(\text{observação incorretamente classificada como } \pi_1) = \\P(\text{observação vem de } \pi_2 \text{ e é incorretamente classificada como } \pi_1)=\\ P(\mathbf{X} \in R_1 \mid \pi_2) P(\pi_2)=\\ P(1 \mid 2)\, p_2 \\ P(\text{observação corretamente classificada como } \pi_2) = \\ P(\text{observação vem de } \pi_2 \text{ e é corretamente classificada como } \pi_2)=\\ P(\mathbf{X} \in R_2 \mid \pi_2) P(\pi_2)= \\ P(2 \mid 2)\, p_2 \\ P(\text{observação incorretamente classificada como } \pi_2) = \\ P(\text{observação vem de } \pi_1 \text{ e é incorretamente classificada como } \pi_2)=\\ P(\mathbf{X} \in R_2 \mid \pi_1) P(\pi_1)=\\ P(2 \mid 1)\, p_1 \tag{11-3} \]

Esquemas de classificação são frequentemente avaliados em termos de suas probabilidades de classificação incorreta (ver Seção 11.4), mas isso ignora o custo da classificação incorreta. Por exemplo, mesmo uma probabilidade aparentemente pequena, como \(0.06 = P(2 \mid 1)\), pode ser grande demais se o custo de fazer uma atribuição incorreta para \(\pi_2\) for extremamente alto. Uma regra que ignora custos pode causar problemas.

Os custos de classificação incorreta podem ser definidos por uma matriz de custos:

\[ \begin{array}{c|cc} \text{Classificar como:} & \pi_1 & \pi_2 \\ \hline \text{População verdadeira: } \pi_1 & 0 & c(2 \mid 1) \\ \pi_2 & c(1 \mid 2) & 0 \end{array} \tag{11-4} \]

Os custos são:

1. zero para classificação correta,
   
2. \(c(1 \mid 2)\) quando uma observação de \(\pi_2\) é incorretamente classificada como \(\pi_1\), e
   
3. \(c(2 \mid 1)\) quando uma observação de \(\pi_1\) é incorretamente classificada como \(\pi_2\).

Para qualquer regra, o custo médio, ou custo esperado de classificação incorreta (ECM), é obtido multiplicando os termos fora da diagonal em (11-4) por suas probabilidades de ocorrência, obtidas em (11-3). Consequentemente,

\[ \text{ECM} = c(2 \mid 1)\,P(2 \mid 1)\,p_1 + c(1 \mid 2)\,P(1 \mid 2)\,p_2 \tag{11-5} \]

Uma regra de classificação razoável deve ter um \(\text{ECM}\) tão pequeno, ou tão próximo de pequeno quanto possível.

Resultado 11.1. As regiões \(R_1\) e \(R_2\) que minimizam o \(\text{ECM}\) são definidas pelos valores de \(\mathbf{x}\) para os quais as seguintes desigualdades são satisfeitas:

$$ R_1:

$$

\[ (\text{razão de densidades}) \;\ge\; (\text{razão de custos}) \;(\text{razão das probabilidades a priori}) \tag{11-6} \]

\[ R_2:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} < \frac{c(1 \mid 2)}{c(2 \mid 1)} \frac{p_2}{p_1} \]

\[ (\text{razão de densidades}) \;<\; (\text{razão de custos}) \;(\text{razão das probabilidades a priori}) \]

Demonstração. Ver Exercício 11.3.

\[\Diamond\] Fica claro, a partir de (11-6), que a implementação da regra de \(\text{ECM}\) mínimo requer:

1. a razão das funções densidade avaliada em uma nova observação \(\mathbf{x}_0\),
   
2. a razão de custos e
   
3. a razão das probabilidades a priori.

O aparecimento de razões na definição das regiões de classificação ótimas é significativo. Frequentemente, é muito mais fácil especificar as razões do que suas partes componentes.

Por exemplo, pode ser difícil especificar os custos (em unidades apropriadas) de classificar um estudante como material de nível universitário quando, na verdade, ele ou ela não o é, e de classificar um estudante como não sendo material universitário quando, na verdade, ele ou ela é. O custo para os contribuintes de educar, por 2 anos, um estudante que abandonará a universidade, por exemplo, pode ser avaliado aproximadamente. Já o custo para a universidade e para a sociedade de não educar um estudante capaz é mais difícil de determinar.

Entretanto, pode ser que um número realista para a razão desses custos de classificação incorreta seja obtido. Quaisquer que sejam as unidades de medida, não admitir um possível formando universitário pode ser cinco vezes mais custoso, ao longo de um horizonte de tempo adequado, do que admitir um eventual evadido. Nesse caso, a razão de custos é cinco.

É interessante considerar as regiões de classificação definidas em (11-6) para alguns casos especiais.

Casos especiais das regiões de custo esperado mínimo:

1. \(p_2/p_1 = 1\) (probabilidades a priori iguais)

\[ R_1:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} \ge \frac{c(1 \mid 2)}{c(2 \mid 1)} \qquad\qquad R_2:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} < \frac{c(1 \mid 2)}{c(2 \mid 1)} \]

2. \(\dfrac{c(1 \mid 2)}{c(2 \mid 1)} = 1\) (custos de classificação incorreta iguais)

\[ R_1:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} \ge \frac{p_2}{p_1} \qquad\qquad R_2:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} < \frac{p_2}{p_1} \tag{11-7} \]

3. \(\dfrac{p_2}{p_1} = \dfrac{c(1 \mid 2)}{c(2 \mid 1)}= 1\) (probabilidades a priori iguais e custos de classificação incorreta iguais)

\[ R_1:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} \ge 1 \qquad\qquad R_2:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} < 1 \]

Quando as probabilidades a priori são desconhecidas, elas são frequentemente tomadas como iguais, e a regra do \(\text{ECM}\) mínimo envolve comparar a razão das densidades populacionais com a razão dos custos apropriados de classificação incorreta.

Se a razão dos custos de classificação incorreta for indeterminada, ela é geralmente tomada como unidade, e a razão das densidades populacionais é comparada com a razão das probabilidades a priori. (Note que as probabilidades a priori aparecem em ordem inversa às densidades.)

Finalmente, quando tanto a razão das probabilidades a priori quanto a razão dos custos de classificação incorreta são iguais a um, as regiões ótimas de classificação são determinadas simplesmente comparando os valores das funções densidade.

Nesse caso, se \(\mathbf{x}_0\) é uma nova observação e

\[ \frac{f_1(\mathbf{x}_0)}{f_2(\mathbf{x}_0)} \ge 1 \]

isto é, \(f_1(\mathbf{x}_0) \ge f_2(\mathbf{x}_0)\), atribuímos \(\mathbf{x}_0\) a \(\pi_1\).

Por outro lado, se

\[ \frac{f_1(\mathbf{x}_0)}{f_2(\mathbf{x}_0)} < 1 \]

ou \(f_1(\mathbf{x}_0) < f_2(\mathbf{x}_0)\), atribuímos \(\mathbf{x}_0\) a \(\pi_2\).

É prática comum usar arbitrariamente o caso (c) de (11-7) para classificação. Isso equivale a assumir probabilidades a priori iguais e custos de classificação incorreta iguais para a regra de \(\text{ECM}\) mínimo. Esta é a justificativa geralmente fornecida. Ela também é equivalente a assumir que a razão das probabilidades a priori é o recíproco da razão dos custos de classificação incorreta.

Exemplo 11.2: Classificando uma nova observação em uma das duas populações

Um pesquisador possui dados suficientes para estimar as funções densidade \(f_1(\mathbf{x})\) e \(f_2(\mathbf{x})\) associadas às populações \(\pi_1\) e \(\pi_2\), respectivamente. Suponha que
\(c(2 \mid 1) = 5\) unidades e \(c(1 \mid 2) = 10\) unidades.

Além disso, sabe-se que cerca de 20% de todos os objetos (para os quais as medidas \(\mathbf{x}\) podem ser registradas) pertencem a \(\pi_2\). Assim, as probabilidades a priori são
\(p_1 = 0.8\) e \(p_2 = 0.2\).

Dadas as probabilidades a priori e os custos de classificação incorreta, podemos usar (11-6) para obter as regiões de classificação \(R_1\) e \(R_2\). Especificamente, temos:

\[ R_1:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} \ge \frac{10}{5}\frac{0.2}{0.8} = 0.5 \]

\[ R_2:\quad \frac{f_1(\mathbf{x})}{f_2(\mathbf{x})} < \frac{10}{5}\frac{0.2}{0.8} = 0.5 \]

Suponha que as funções densidade avaliadas em uma nova observação \(\mathbf{x}_0\) forneçam
\(f_1(\mathbf{x}_0) = 0.3\) e \(f_2(\mathbf{x}_0) = 0.4\).

Devemos classificar a nova observação como pertencente a \(\pi_1\) ou a \(\pi_2\)?
Para responder a essa questão, formamos a razão:

\[ \frac{f_1(\mathbf{x}_0)}{f_2(\mathbf{x}_0)} = \frac{0.3}{0.4} = 0.75 \]

e a comparamos com \(0.5\) obtido anteriormente. Como

\[ \frac{f_1(\mathbf{x}_0)}{f_2(\mathbf{x}_0)} = 0.75 > \frac{c(1 \mid 2)}{c(2 \mid 1)}\frac{p_2}{p_1} = 0.5 \]

concluímos que \(\mathbf{x}_0 \in R_1\) e a classificamos como pertencente a \(\pi_1\).

\[\Diamond\] Critérios além do custo esperado de classificação incorreta podem ser usados para derivar procedimentos de classificação “ótimos”. Por exemplo, pode-se ignorar os custos de classificação incorreta e escolher \(R_1\) e \(R_2\) de modo a minimizar a probabilidade total de classificação incorreta (TPM):

\[ \begin{align} \text{TPM} &= P(\text{classificar incorretamente uma observação de } \pi_1\\ &\qquad\text{ou classificar incorretamente uma observação de } \pi_2)\\ &= P(\text{observação vem de } \pi_1 \text{ e é classificada incorretamente}) + \\ &\quad P(\text{observação vem de } \pi_2 \text{ e é classificada incorretamente})\\ \text{TPM}&= p_1 \int_{R_2} f_1(\mathbf{x})\, d\mathbf{x} + p_2 \int_{R_1} f_2(\mathbf{x})\, d\mathbf{x} \end{align} \tag{11-8} \]

Do ponto de vista matemático, este problema é equivalente a minimizar o custo esperado de classificação incorreta quando os custos de classificação incorreta são iguais. Consequentemente, as regiões ótimas, nesse caso, são dadas por (b) em (11-7).

Também podemos alocar uma nova observação \(\mathbf{x}_0\) à população com a maior probabilidade a posteriori \(P(\pi_i \mid \mathbf{x}_0)\). Pela regra de Bayes, as probabilidades a posteriori são:

\[ \begin{align} P(\pi_1 \mid \mathbf{x}_0) &= \frac{P(\pi_1 \text{ ocorre e observamos } \mathbf{x}_0)} {P(\text{observamos } \mathbf{x}_0)}\\ &= \frac{P(\mathbf{x}_0 \mid \pi_1)\,P(\pi_1)} {P(\mathbf{x}_0 \mid \pi_1)\,P(\pi_1) + P(\mathbf{x}_0 \mid \pi_2)\,P(\pi_2)}\\ P(\pi_1 \mid \mathbf{x}_0)&= \frac{p_1 f_1(\mathbf{x}_0)} {p_1 f_1(\mathbf{x}_0) + p_2 f_2(\mathbf{x}_0)}\\\\ P(\pi_2 \mid \mathbf{x}_0)&= 1 - P(\pi_1 \mid \mathbf{x}_0)\\ P(\pi_2 \mid \mathbf{x}_0)&= \frac{p_2 f_2(\mathbf{x}_0)} {p_1 f_1(\mathbf{x}_0) + p_2 f_2(\mathbf{x}_0)} \end{align} \tag{11-9} \]

Classificar uma observação \(\mathbf{x}_0\) como pertencente a \(\pi_1\) quando \(P(\pi_1 \mid \mathbf{x}_0) > P(\pi_2 \mid \mathbf{x}_0)\) é equivalente a usar a regra (b) para a probabilidade total de classificação incorreta em (11-7), pois os denominadores em (11-9) são os mesmos.

Entretanto, calcular as probabilidades a posteriori das populações \(\pi_1\) e \(\pi_2\) após observar \(\mathbf{x}_0\) (daí o nome probabilidade a posteriori) é frequentemente útil para identificar atribuições menos evidentes.

Classificação com duas populações normais multivariadas

Procedimentos de classificação baseados em populações normais predominam na prática estatística devido à sua simplicidade e à eficiência razoavelmente alta em uma ampla variedade de modelos populacionais. Passamos agora a supor que \(f_1(\mathbf{x})\) e \(f_2(\mathbf{x})\) são densidades normais multivariadas, a primeira com vetor média \(\boldsymbol{\mu}_1\) e matriz de covariâncias \(\boldsymbol{\Sigma}_1\) e a segunda com vetor média \(\boldsymbol{\mu}_2\) e matriz de covariâncias \(\boldsymbol{\Sigma}_2\).

O caso especial de matrizes de covariâncias iguais leva a uma estatística de classificação linear particularmente simples.

Classificação de populações normais quando \(\boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2 = \boldsymbol{\Sigma}\)

Suponha que as densidades conjuntas de \(\mathbf{X}^{\prime} = [X_1, X_2, \ldots, X_p]\) para as populações \(\pi_1\) e \(\pi_2\) sejam dadas por

\[ \begin{align} f_i(\mathbf{x}) &= \frac{1}{((2\pi)^{p}\,|\boldsymbol{\Sigma}|)^{1/2}} \exp\left[ -\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu}_i)^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}_i) \right]\\ i &= 1, 2 \end{align} \tag{11-10} \]

Suponha também que os parâmetros populacionais \(\boldsymbol{\mu}_1\), \(\boldsymbol{\mu}_2\) e \(\boldsymbol{\Sigma}\) sejam conhecidos. Então, após o cancelamento dos termos \((2\pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}\), as regiões de \(\text{ECM}\) mínimo em (11-6) tornam-se

\[ \begin{align} R_1:&\quad\exp\left[ -\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu}_1)^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}_1) +\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}_2) \right] \\ &\quad\ge \frac{c(1 \mid 2)}{c(2 \mid 1)} \frac{p_2}{p_1} \\[4pt] R_2: &\quad\exp\left[ -\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu}_1)^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}_1) +\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}_2) \right] \\ &\quad< \frac{c(1 \mid 2)}{c(2 \mid 1)} \frac{p_2}{p_1} \end{align} \tag{11-11} \] Dadas essas regiões \(R_1\) e \(R_2\), podemos construir a regra de classificação apresentada no seguinte resultado.

Resultado 11.2. Sejam as populações \(\pi_1\) e \(\pi_2\) descritas por densidades normais multivariadas da forma (11-10). Então, a regra de alocação que minimiza o \(\text{ECM}\) é a seguinte:

Aloque \(\mathbf{x}_0\) a \(\pi_1\) se

\[ (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{x}_0 - \frac{1}{2} (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2) \ge \ln\left( \frac{c(1 \mid 2)}{c(2 \mid 1)} \frac{p_2}{p_1} \right) \tag{11-12} \]

Aloque \(\mathbf{x}_0\) a \(\pi_2\) em caso contrário.

Demonstração. Como as quantidades em (11-11) são não negativas para todo \(\mathbf{x}\), podemos tomar seus logaritmos naturais e preservar a ordem das desigualdades. Além disso (ver Exercício 11.5),

\[ -\frac{1}{2} (\mathbf{x} - \boldsymbol{\mu}_1)^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}_1) + \frac{1}{2} (\mathbf{x} - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}_2)=\\ (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{x} - \frac{1}{2} (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2) \tag{11-13} \]

e, consequentemente,

\[ \begin{align} R_1:\quad &(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{x} - \frac{1}{2} (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2)\\ &\ge \ln\left( \frac{c(1 \mid 2)}{c(2 \mid 1)} \frac{p_2}{p_1} \right) \\[4pt] R_2:\quad &(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{x} - \frac{1}{2} (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)^{\prime} \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2)\\ &< \ln\left( \frac{c(1 \mid 2)}{c(2 \mid 1)} \frac{p_2}{p_1} \right) \end{align} \tag{11-14} \]

Segue daí a regra de classificação de \(\text{ECM}\) mínimo.

\[\Diamond\]

Na maioria das situações práticas, as quantidades populacionais \(\boldsymbol{\mu}_1\), \(\boldsymbol{\mu}_2\) e \(\boldsymbol{\Sigma}\) são desconhecidas, de modo que a regra (11-12) deve ser modificada. Wald [31] e Anderson [2] sugeriram substituir os parâmetros populacionais por seus correspondentes amostrais.

Suponha, então, que tenhamos \(n_1\) observações da variável aleatória multivariada \(\mathbf{X}^{\prime} = [X_1, X_2, \ldots, X_p]\) provenientes de \(\pi_1\) e \(n_2\) medições dessa quantidade provenientes de \(\pi_2\), com \(n_1 + n_2 - 2 \ge p\). Então, as respectivas matrizes de dados são

\[ \underset{n_1 \times p}{\mathbf{X}_1} = \begin{bmatrix} \mathbf{X}_{11}^{\prime} \\ \mathbf{X}_{12}^{\prime} \\ \vdots \\ \mathbf{X}_{1 n_1}^{\prime} \end{bmatrix} \qquad\text{e}\qquad \underset{n_{2} \times p}{\mathbf{X}_2} = \begin{bmatrix} \mathbf{X}_{21}^{\prime} \\ \mathbf{X}_{22}^{\prime} \\ \vdots \\ \mathbf{X}_{2 n_2}^{\prime} \end{bmatrix} \tag{11-15} \]

A partir dessas matrizes de dados, os vetores médias amostrais e as matrizes de covariâncias amostrais são dados por

\[ \begin{align} \bar{\mathbf{x}}_1 &= \frac{1}{n_1} \sum_{j=1}^{n_1} \mathbf{x}_{1j} &\qquad \mathbf{s}_1 &= \frac{1}{n_1 - 1} \sum_{j=1}^{n_1} (\mathbf{x}_{1j} - \bar{\mathbf{x}}_1)(\mathbf{x}_{1j} - \bar{\mathbf{x}}_1)^{\prime} \\[6pt] \bar{\mathbf{x}}_2 &= \frac{1}{n_2} \sum_{j=1}^{n_2} \mathbf{x}_{2j} &\qquad \mathbf{s}_2 &= \frac{1}{n_2 - 1} \sum_{j=1}^{n_2} (\mathbf{x}_{2j} - \bar{\mathbf{x}}_2)(\mathbf{x}_{2j} - \bar{\mathbf{x}}_2)^{\prime} \end{align} \tag{11-16} \]

Como se assume que as populações de origem possuem a mesma matriz de covariância \(\boldsymbol{\Sigma}\), as matrizes de covariância amostrais \(\mathbf{s}_1\) e \(\mathbf{s}_2\) são combinadas (pooled) para obter uma única estimativa não viciada de \(\boldsymbol{\Sigma}\), como em (6-21). Em particular, a média ponderada

\[ \mathbf{s}_{\text{pooled}} = \frac{n_1 - 1}{n_1 + n_2 - 2} \, \mathbf{s}_1 + \frac{n_2 - 1}{n_1 + n_2 - 2} \, \mathbf{s}_2 \tag{11-17} \]

é uma estimativa não viciada de \(\boldsymbol{\Sigma}\) se as matrizes de dados \(\mathbf{X}_1\) e \(\mathbf{X}_2\) contiverem amostras aleatórias das populações \(\pi_1\) e \(\pi_2\), respectivamente.

Substituindo \(\bar{\mathbf{x}}_1\) por \(\boldsymbol{\mu}_1\), \(\bar{\mathbf{x}}_2\) por \(\boldsymbol{\mu}_2\) e \(\mathbf{s}_{\text{pooled}}\) por \(\boldsymbol{\Sigma}\) em (11-12), obtemos a regra de classificação “amostral”:

Regra estimada de \(\text{ECM}\) mínimo para duas populações normais:

Aloque \(\mathbf{x}_0\) a \(\pi_1\) se

\[ (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\mathbf{x}_0 -\frac{1}{2} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} (\bar{\mathbf{x}}_1 + \bar{\mathbf{x}}_2) \\ \ge \ln\left( \frac{c(1 \mid 2)}{c(2 \mid 1)}\frac{p_2}{p_1} \right) \tag{11-18} \]

Aloque \(\mathbf{x}_0\) a \(\pi_2\) caso contrário.

Se, em (11-18),

\[ \frac{c(1 \mid 2)}{c(2 \mid 1)}\frac{p_2}{p_1} = 1 \]

então \(\ln(1)=0\) e a regra estimada de \(\text{ECM}\) mínimo reduz-se a comparar a variável escalar

\[ \begin{align} \hat{y} &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\mathbf{x}\\ \hat{y}&= \hat{\mathbf{a}}^{\prime}\mathbf{x} \end{align} \tag{11-19} \]

avaliada em \(\mathbf{x}_0\), com o valor

\[ \begin{align} \hat{m} &= \frac{1}{2} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} (\bar{\mathbf{x}}_1 + \bar{\mathbf{x}}_2) \\ \hat{m}&= \frac{1}{2}(\bar{y}_1 + \bar{y}_2) \end{align} \tag{11-20} \]

sendo que

\[ \bar{y}_1 = (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\bar{\mathbf{x}}_1 = \hat{\mathbf{a}}^{\prime}\bar{\mathbf{x}}_1\\ \bar{y}_2 = (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\bar{\mathbf{x}}_2 = \hat{\mathbf{a}}^{\prime}\bar{\mathbf{x}}_2 \]

Isto é, a regra estimada de \(\text{ECM}\) mínimo para duas populações normais equivale a construir duas populações univariadas para os valores \(y\), tomando uma combinação linear apropriada das observações de \(\pi_1\) e \(\pi_2\), e então classificar uma nova observação \(\mathbf{x}_0\) em \(\pi_1\) ou \(\pi_2\), dependendo se \(\hat{y}_0 = \hat{\mathbf{a}}^{\prime}\mathbf{x}_0\) cai à direita ou à esquerda do ponto médio entre as duas médias univariadas \(\bar{y}_1\) e \(\bar{y}_2\).

Uma vez inseridas estimativas dos parâmetros populacionais correspondentes, não há garantia de que a regra resultante minimiza o custo esperado de classificação incorreta em uma aplicação particular. Isso ocorre porque a regra ótima em (11-12) foi derivada assumindo que as densidades normais multivariadas \(f_1(\mathbf{x})\) e \(f_2(\mathbf{x})\) eram completamente conhecidas. A expressão (11-18) é apenas uma estimativa da regra ótima. Ainda assim, é razoável esperar que ela tenha bom desempenho quando os tamanhos amostrais forem grandes. À medida que os tamanhos amostrais aumentam, \(\bar{\mathbf{x}}_1\), \(\bar{\mathbf{x}}_2\) e \(\mathbf{S}_{\text{pooled}}\) tornam-se, com probabilidade tendendo a 1, indistinguíveis de \(\boldsymbol{\mu}_1\), \(\boldsymbol{\mu}_2\) e \(\boldsymbol{\Sigma}\), respectivamente [ver (4-26) e (4-27)].

Em resumo, se os dados parecerem ser normais multivariados, a estatística de classificação do lado esquerdo da desigualdade em (11-18) pode ser calculada para cada nova observação \(\mathbf{x}_0\). Pelo menos, as distribuições de frequência marginais das observações em cada variável devem ser verificadas quanto à normalidade. Isso deve ser feito para as amostras em ambas as populações. Frequentemente, algumas variáveis precisam ser transformadas para torná-las mais “normais”. (Ver Seções 4.6 e 4.8.)

Essas observações são então classificadas comparando-se os valores dessa estatística com o valor de \(\ln\left(\frac{c(1 \mid 2)}{c(2 \mid 1)}\frac{p_2}{p_1}\right)\).

Exemplo 11.3: Classificação com duas populações normais: \(\boldsymbol{\Sigma}\) comum e custos iguais)

Este exemplo é adaptado de um estudo [4] relacionado à detecção de portadoras de hemofilia A. (Ver também Exercício 11.32.)

Para construir um procedimento para detectar potenciais portadoras de hemofilia A, amostras de sangue foram analisadas para dois grupos de mulheres, e medidas em duas variáveis foram registradas:

\[ \begin{align} X_1 &= \log_{10}(\text{atividade de AHF}) \\ X_2 &= \log_{10}(\text{antígeno semelhante ao AHF}) \end{align} \]

sendo que AHF denota fator anti-hemofílico.

O primeiro grupo, com \(n_1 = 30\) mulheres, foi selecionado de uma população de mulheres que não carregavam o gene da hemofilia. Esse grupo foi denominado grupo normal.

O segundo grupo, com \(n_2 = 45\) mulheres, foi selecionado entre portadoras conhecidas de hemofilia A (filhas de hemofílicos, mães com mais de um filho hemofílico, e mães com um filho hemofílico e outros parentes hemofílicos). Esse grupo foi denominado portadoras obrigatórias.

Os pares de observações \((x_1, x_2)\) para os dois grupos estão apresentados na Figura 11.4. Também são mostradas as curvas de nível estimadas contendo 50% e 95% da probabilidade para distribuições normais bivariadas centradas em \(\bar{\mathbf{x}}_1\) e \(\bar{\mathbf{x}}_2\), respectivamente. A matriz de covariância comum foi tomada como a matriz de covariância amostral combinada \(\mathbf{s}_{\text{pooled}}\). Neste exemplo, distribuições normais bivariadas parecem ajustar-se razoavelmente bem aos dados.

Os investigadores [4] fornecem a seguinte informação:

\[ \bar{\mathbf{x}}_1 = \begin{bmatrix} -0.1349 \\ -0.0778 \end{bmatrix} \qquad \bar{\mathbf{x}}_2 = \begin{bmatrix} -0.3079 \\ -0.0060 \end{bmatrix} \]

e

\[ \mathbf{s}_{\text{pooled}}^{-1} = \begin{bmatrix} 86.090 & -61.487 \\ -61.487 & 90.199 \end{bmatrix} \]

Portanto, a função discriminante com custos iguais e probabilidades a priori iguais (ver (11-19)) é

\[ \begin{align} \hat{y} &= \hat{\mathbf{a}}^{\prime}\mathbf{x}\\ &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\mathbf{x} \\[4pt] &= [0.173 \;-0.0719] \begin{bmatrix} 86.090 & -61.487 \\ -61.487 & 90.199 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\[4pt] \hat{y}&= 19.319 \,x_1 - 17.124\,x_2 \end{align} \]

Além disso,

\[ \begin{align} \bar{y}_1 &= \hat{\mathbf{a}}^{\prime}\bar{\mathbf{x}}_1 = [\,19.319 \;\; -17.124\,] \begin{bmatrix} -0.1349 \\ -0.0778 \end{bmatrix} = -1.27 \\[6pt] \bar{y}_2 &= \hat{\mathbf{a}}^{\prime}\bar{\mathbf{x}}_2 = [\,19.319 \;\; -17.124\,] \begin{bmatrix} -0.3079 \\ -0.0060 \end{bmatrix} = -5.85 \end{align} \]

O ponto médio entre essas médias (ver (11-20)) é

\[ \hat{m} = \frac{1}{2}(\bar{y}_1 + \bar{y}_2) = \frac{1}{2}(-1.27 - 5.85) = -3.56 \]

As medições de atividade de AHF e de antígeno semelhante ao AHF em uma mulher que pode ser portadora de hemofilia A fornecem
\(x_1 = -0.210\) e \(x_2 = -0.044\).

Devemos classificar essa mulher em \(\pi_1\) (normal) ou em \(\pi_2\) (portadora obrigatória)?

Usando (11-18) com custos iguais e probabilidades a priori iguais, de modo que \(\ln(1)=0\), obtemos:

\[ \text{Alocar } \mathbf{x}_0 \text{ a } \pi_1 \text{ se } \hat{y}_0 = \hat{\mathbf{a}}^{\prime}\mathbf{x}_0 \ge \hat{m} = -3.56 \]

\[ \text{Alocar } \mathbf{x}_0 \text{ a } \pi_2 \text{ se } \hat{y}_0 = \hat{\mathbf{a}}^{\prime}\mathbf{x}_0 < \hat{m} = -3.56 \]

sendo que \(\mathbf{x}_0^{\prime} = [-0.210,\,-0.044]\).

Como

\[ \hat{y}_0 = [\,19.319 \;\; -17.124\,] \begin{bmatrix} -0.210 \\ -0.044 \end{bmatrix} = -3.30 > -3.56 \]

classificamos a mulher como pertencente a \(\pi_1\), isto é, uma não portadora.

Suponha agora que as probabilidades a priori de pertencimento aos grupos sejam conhecidas.

Por exemplo, para uma prima em primeiro grau materna de um hemofílico, temos \(p_1 = 0.75\) e \(p_2 = 0.25\).

Assumindo custos de erro iguais e usando a estatística

\[ \begin{align} \hat{w} &= \hat{\mathbf{a}}^{\prime}\mathbf{x}_0 - \hat{m} \end{align} \]

com \(\hat{\mathbf{a}}^{\prime}\mathbf{x}_0 = -3.30\) e \(\hat{m} = -3.56\), obtemos

\[ \hat{w} = -3.30 - (-3.56) = 0.26 \]

Aplicando (11-18):

\[ \hat{w} = 0.26 > \ln\left(\frac{p_2}{p_1}\right) = \ln\left(\frac{0.25}{0.75}\right) = \ln\left(\frac{1}{3}\right) = -1.10 \]

Portanto, mesmo com probabilidades a priori desiguais, a classificação é novamente \(\pi_1\), isto é, não portadora.

\[\Diamond\] * JW6_Example11.3.R

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
))
Dados <- read.csv("JW6Data/T11-8.dat", header = FALSE, sep = "")
colnames(Dados) <- c("Group", "AHF_act", "AHF_antig")

Dados$Group <- factor(Dados$Group,
                      labels = c("noncarriers", "oblig.carriers"))

print(psych::describeBy(AHF_antig+AHF_act~Group, 
                        mat=1,
                        digits=2,
                        data = Dados))

           item         group1 vars  n  mean   sd median trimmed  mad   min
AHF_antig1    1    noncarriers    1 30 -0.08 0.13  -0.07   -0.07 0.13 -0.48
AHF_antig2    2 oblig.carriers    1 45 -0.01 0.16   0.00    0.00 0.17 -0.34
AHF_act1      3    noncarriers    2 30 -0.13 0.14  -0.14   -0.13 0.08 -0.53
AHF_act2      4 oblig.carriers    2 45 -0.31 0.15  -0.34   -0.31 0.17 -0.69
             max range  skew kurtosis   se
AHF_antig1  0.21  0.69 -0.59     1.09 0.02
AHF_antig2  0.29  0.63 -0.24    -0.83 0.02
AHF_act1    0.15  0.68 -0.61     0.93 0.03
AHF_act2   -0.01  0.68  0.00    -0.53 0.02

suppressMessages(suppressWarnings(
  invisible(
    print(GGally::ggpairs(data = Dados,
                          mapping = ggplot2::aes(color = Group,
                                                 alpha = .5),
                          progress = FALSE)+ 
            ggplot2::theme_bw()+
            ggplot2::theme(panel.grid = 
                             ggplot2::element_blank())))))

out <- MVN::mvn(data = Dados,
                subset = "Group",
                univariate_test = "SW",
                desc = FALSE)
print(out$multivariate_normality, digits=3)

           Group          Test Statistic p.value      MVN
1    noncarriers Henze-Zirkler     0.427   0.490 ✓ Normal
2 oblig.carriers Henze-Zirkler     0.529   0.364 ✓ Normal

res <- heplots::boxM(Y=Dados[,2:3],
                     group=Dados$Group)
print(res)

    Box's M-test for Homogeneity of Covariance Matrices

data:  Dados[, 2:3]
Chi-Sq (approx.) = 5.3383, df = 3, p-value = 0.1486

plot(res)

car::scatterplot(data=Dados,
                 AHF_antig~AHF_act|Group,
                 regLine=FALSE, 
                 ellipse=list(levels=c(.68), 
                              robust=TRUE, fill=FALSE), 
                 grid=FALSE,
                 col=c("black", "grey30"),
                 smooth=FALSE)

# As colunas já estão em log10 no arquivo, então só renomeio:
Dados$X1 <- Dados$AHF_act
Dados$X2 <- Dados$AHF_antig

# Separar grupos
X1mat <- as.matrix(subset(Dados, Group == "noncarriers")[, c("X1", "X2")])
X2mat <- as.matrix(subset(Dados, Group == "oblig.carriers")[, c("X1", "X2")])

n1 <- nrow(X1mat)
n2 <- nrow(X2mat)

# Vetores média amostral
xbar1 <- colMeans(X1mat)
xbar2 <- colMeans(X2mat)

xbar1

         X1          X2 
-0.13487000 -0.07785667 

xbar2

          X1           X2 
-0.307946667 -0.005991111 

# Matrizes de covariância amostral
S1 <- cov(X1mat)
S2 <- cov(X2mat)

# Matriz de covariância combinada (pooled)
S_pooled <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n1 + n2 - 2)
S_pooled_inv <- solve(S_pooled)

S_pooled_inv

          X1        X2
X1  86.09006 -61.48727
X2 -61.48727  90.19932

# Vetor discriminante â = (x̄1 - x̄2)' S_pooled^{-1}
a_hat <- t(xbar1 - xbar2) %*% S_pooled_inv   # 1 x 2
a_hat   

         X1        X2
[1,] 19.319 -17.12424

# Função discriminante  ŷ = â' x
y_hat <- function(x) {
  as.numeric(a_hat %*% x)
}

# Médias discriminantes dos grupos
ybar1 <- y_hat(xbar1)  
ybar2 <- y_hat(xbar2)   

ybar1; ybar2

[1] -1.272317

[1] -5.846628

# Ponto médio m̂
m_hat <- 0.5 * (ybar1 + ybar2)  
m_hat

[1] -3.559472

# Nova observação do exemplo
x0 <- c(-0.210, -0.044)

y0 <- y_hat(x0)   # ŷ0
y0                

[1] -3.303523

## 1) Custos iguais, priors iguais
class_equal <- if (y0 >= m_hat) "pi1 (noncarriers)" else "pi2 (oblig.carriers)"
class_equal   # deve dar "pi2 (oblig.carriers)"

[1] "pi1 (noncarriers)"

## 2) Priors p1 = 0.75, p2 = 0.25, custos iguais

p1 <- 0.75
p2 <- 0.25

w_hat <- function(x) {
  term1 <- (xbar1 - xbar2) %*% S_pooled_inv %*% x
  term2 <- 0.5 * (xbar1 - xbar2) %*% S_pooled_inv %*% (xbar1 + xbar2)
  as.numeric(term1 - term2)
}

w0 <- w_hat(x0)          
w0

[1] 0.2559494

threshold <- log(p2 / p1) 
threshold

[1] -1.098612

class_priors <- if (w0 >= threshold) "pi1 (noncarriers)" else "pi2 (oblig.carriers)"
class_priors   # também "pi2 (oblig.carriers)"

[1] "pi1 (noncarriers)"

Escalonamento (Scaling)

O vetor de coeficientes discriminantes é

\[ \hat{\mathbf{a}} = \mathbf{s}_{\text{pooled}}^{-1}(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2) \]

e é único apenas até uma constante multiplicativa. Assim, para qualquer \(c \neq 0\), o vetor \(c\,\hat{\mathbf{a}}\) também serve como vetor de coeficientes discriminantes.

Com frequência, o vetor \(\hat{\mathbf{a}}\) é “escalonado” ou “normalizado” para facilitar a interpretação de seus elementos. Duas das normalizações mais utilizadas são:

1. Definir

\[ \hat{\mathbf{a}}^{*} = \frac{\hat{\mathbf{a}}}{\sqrt{\hat{\mathbf{a}}^{\prime}\hat{\mathbf{a}}}} \quad (11\text{-}21) \]

de modo que \(\hat{\mathbf{a}}^{*}\) tenha norma unitária.

2. Definir

\[ \hat{\mathbf{a}}^{*} = \frac{\hat{\mathbf{a}}}{\hat{a}_1} \quad (11\text{-}22) \]

de modo que o primeiro elemento do novo vetor de coeficientes \(\hat{\mathbf{a}}^{*}\) seja igual a 1.

Em ambos os casos, \(\hat{\mathbf{a}}^{*}\) é da forma \(c\,\hat{\mathbf{a}}\).

Para a normalização (1):

\[ c = (\hat{\mathbf{a}}^{\prime}\hat{\mathbf{a}})^{-1/2} \]

Para a normalização (2):

\[ c = \hat{a}_1^{-1} \]

As magnitudes de \(\hat{a}^{*}_{1}, \hat{a}^{*}_{2}, \ldots, \hat{a}^{*}_{p}\) em (11-21) pertencem todas ao intervalo \([-1,1]\).
Em (11-22), \(\hat{a}^{*}_{1} = 1\) e \(\hat{a}^{*}_{2}, \ldots, \hat{a}^{*}_{p}\) são expressos como múltiplos de \(\hat{a}^{*}_{1}\).

Restringir os \(\hat{a}^{*}_{i}\) ao intervalo \([-1,1]\) geralmente facilita uma comparação visual dos coeficientes.

De modo semelhante, expressar os coeficientes como múltiplos de \(\hat{a}^{*}_{1}\) permite avaliar facilmente a importância relativa (em relação à \(X_1\)) das variáveis \(X_2, \ldots, X_p\) como discriminantes.

A normalização dos \(\hat{a}_i\) é recomendada apenas se as variáveis \(X\) tiverem sido padronizadas.

Se esse não for o caso, é necessário extremo cuidado na interpretação dos resultados.

• JW6_Example11.3_Scaling.R

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
))
Dados <- read.csv("JW6Data/T11-8.dat", header = FALSE, sep = "")
colnames(Dados) <- c("Group", "AHF_act", "AHF_antig")

Dados$Group <- factor(Dados$Group,
                      labels = c("noncarriers", "oblig.carriers"))

# As colunas já estão em log10 no arquivo, então só renomeio:
Dados$X1 <- Dados$AHF_act
Dados$X2 <- Dados$AHF_antig

# Separar grupos
X1mat <- as.matrix(subset(Dados, Group == "noncarriers")[, c("X1", "X2")])
X2mat <- as.matrix(subset(Dados, Group == "oblig.carriers")[, c("X1", "X2")])

n1 <- nrow(X1mat)
n2 <- nrow(X2mat)

# Vetores média amostral
xbar1 <- colMeans(X1mat)
xbar2 <- colMeans(X2mat)

# Matrizes de covariância amostral
S1 <- cov(X1mat)
S2 <- cov(X2mat)

# Matriz de covariância combinada (pooled)
S_pooled <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n1 + n2 - 2)

# Nova observação do exemplo
x0 <- data.frame(X1 = -0.210,
                 X2 = -0.044)

############################################
## (1) LDA com custos iguais e priors iguais
############################################

fit_eq <- MASS::lda(Group ~ X1 + X2,
                    data  = Dados,
                    prior = c(0.5, 0.5))   # priors iguais

print(fit_eq)

Call:
lda(Group ~ X1 + X2, data = Dados, prior = c(0.5, 0.5))

Prior probabilities of groups:
   noncarriers oblig.carriers 
           0.5            0.5 

Group means:
                       X1           X2
noncarriers    -0.1348700 -0.077856667
oblig.carriers -0.3079467 -0.005991111

Coefficients of linear discriminants:
         LD1
X1 -9.032787
X2  8.006605

# Médias projetadas (escala lda)
ybar <- as.numeric(as.matrix(fit_eq$means) %*% fit_eq$scaling)
names(ybar) <- rownames(fit_eq$means)
ybar

   noncarriers oblig.carriers 
     0.5948844      2.7336482 

# ponto de corte m_hat na escala lda (priors iguais)
m_hat_eq <- 0.5 * (ybar["noncarriers"] + ybar["oblig.carriers"])
m_hat_eq

noncarriers 
   1.664266 

# classificação da nova observação com priors iguais
pred_eq <- predict(fit_eq, newdata = x0)
pred_eq$class

[1] noncarriers
Levels: noncarriers oblig.carriers

pred_eq$posterior

  noncarriers oblig.carriers
1   0.5636403      0.4363597

pred_eq$x        # escore discriminante da nova observação

         LD1
1 -0.1196716

a_manual <- drop(t(xbar1 - xbar2) %*% solve(S_pooled))
a_manual

       X1        X2 
 19.31900 -17.12424 

a_mass  <- drop(fit_eq$scaling)  # vetor 2x1
a_mass

       X1        X2 
-9.032787  8.006605 

# fator de reescala:
k <- a_manual[1] / a_mass[1]
as.numeric(k)

[1] -2.138764

# reescalar os coeficientes de lda para obter exatamente o vetor do livro
a_recuperado <- k * a_mass
a_recuperado

       X1        X2 
 19.31900 -17.12424 

############################################
## (2) LDA com priors p1 = 0.75, p2 = 0.25
############################################

fit_pr <- MASS::lda(Group ~ X1 + X2,
                    data  = Dados,
                    prior = c(0.75, 0.25))  # p1 = 0.75, p2 = 0.25

print(fit_pr)

Call:
lda(Group ~ X1 + X2, data = Dados, prior = c(0.75, 0.25))

Prior probabilities of groups:
   noncarriers oblig.carriers 
          0.75           0.25 

Group means:
                       X1           X2
noncarriers    -0.1348700 -0.077856667
oblig.carriers -0.3079467 -0.005991111

Coefficients of linear discriminants:
         LD1
X1 -9.032787
X2  8.006605

# classificação da nova observação com priors desiguais
pred_pr <- predict(fit_pr, newdata = x0)
pred_pr$class

[1] noncarriers
Levels: noncarriers oblig.carriers

pred_pr$posterior

  noncarriers oblig.carriers
1   0.7948744      0.2051256

pred_pr$x        # escore discriminante na mesma direção (apenas o cutoff muda)

        LD1
1 0.4150193

Abordagem de Fisher para classificação com duas populações

Fisher chegou à estatística de classificação linear (11-19) usando um argumento totalmente diferente.

A ideia de Fisher era transformar as observações multivariadas \(\mathbf{x}\) em observações univariadas \(y\), de modo que os valores de \(y\) derivados das populações \(\pi_1\) e \(\pi_2\) fossem o mais separados possível.

Fisher sugeriu tomar combinações lineares de \(\mathbf{x}\) para criar \(y\), pois elas são funções suficientemente simples de \(\mathbf{x}\) para serem tratadas facilmente.

Sua abordagem não assume que as populações sejam normais. No entanto, assume implicitamente que as matrizes de covariância das populações sejam iguais, pois utiliza-se uma estimativa combinada da matriz de covariância comum.

Uma combinação linear fixa dos \(\mathbf{x}\) assume os valores
\(y_{11}, y_{12}, \ldots, y_{1n_1}\) para as observações da primeira população, e \(y_{21}, y_{22}, \ldots, y_{2n_2}\) para as observações da segunda população.

A separação desses dois conjuntos de \(y\) univariados é avaliada em termos da diferença entre \(\bar{y}_1\) e \(\bar{y}_2\), expressa em unidades de desvio padrão. Isto é,

\[ \text{separação} = \frac{|\bar{y}_1 - \bar{y}_2|}{s_y} \]

sendo que

\[ s_y^2 = \frac{\sum_{j=1}^{n_1}(y_{1j} - \bar{y}_1)^2 + \sum_{j=1}^{n_2}(y_{2j} - \bar{y}_2)^2} {n_1 + n_2 - 2} \]

é a estimativa combinada da variância.

O objetivo é selecionar a combinação linear de \(\mathbf{x}\) que maximize a separação entre as médias amostrais \(\bar{y}_1\) e \(\bar{y}_2\).

Resultado 11.3. A combinação linear

\[ \hat{y} = \hat{\mathbf{a}}^{\prime}\mathbf{x} = (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\mathbf{x} \]

maximiza a razão

\[ \begin{align} \frac{\text{distância^2 entre médias amostrais de } y} {\text{variância amostral de } y} &= \frac{(\bar{y}_1 - \bar{y}_2)^2}{s_y^2}\\ &= \frac{(\hat{\mathbf{a}}^{\prime}\bar{\mathbf{x}}_1 - \hat{\mathbf{a}}^{\prime}\bar{\mathbf{x}}_2)^2} {\hat{\mathbf{a}}^{\prime}\mathbf{s}_{\text{pooled}}\hat{\mathbf{a}}}\\ &= \frac{(\hat{\mathbf{a}}^{\prime}\mathbf{d})^2} {\hat{\mathbf{a}}^{\prime}\mathbf{s}_{\text{pooled}}\hat{\mathbf{a}}} \end{align} \tag{11-23} \]

sobre todos os vetores de coeficientes possíveis \(\hat{\mathbf{a}}\), sendo que \(\mathbf{d} = \bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2\).

O máximo da razão (11-23) é

\[ d^2 = (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2) \] Demonstração. O máximo da razão em (11-23) é obtido aplicando-se diretamente (2-50).

Assim, definindo \(\mathbf{d} = (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)\), temos

\[ \begin{align} \max_{\hat{\mathbf{a}}} \frac{(\hat{\mathbf{a}}^{\prime}\mathbf{d})^2} {\hat{\mathbf{a}}^{\prime}\mathbf{s}_{\text{pooled}}\hat{\mathbf{a}}} &= \mathbf{d}^{\prime}\mathbf{s}_{\text{pooled}}^{-1}\mathbf{d} \\ &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)\\ &= d^2 \end{align} \]

sendo que \(d^2\) é a distância quadrática amostral entre as duas médias.

\[\Diamond\]

Observe que \(s_y^2\) em (11-33) pode ser calculado como

\[ s_y^2 = \frac{ \displaystyle \sum_{j=1}^{n_1} (y_{1j} - \bar{y}_1)^2 + \sum_{j=1}^{n_2} (y_{2j} - \bar{y}_2)^2 }{ n_1 + n_2 - 2 } \tag{11-24} \]

com \(y_{1j} = \hat{\mathbf{a}}^{\prime}\mathbf{x}_{1j}\) e \(y_{2j} = \hat{\mathbf{a}}^{\prime}\mathbf{x}_{2j}\).

Exemplo 11.4: Função discriminante linear de Fisher para os dados de hemofilia

Considere a detecção de portadoras de hemofilia A introduzida no Exemplo 11.3.

Lembre que a função discriminante linear de Fisher, com custos iguais e probabilidades a priori iguais, era

\[ \begin{align} \hat{y} &= \hat{\mathbf{a}}^{\prime}\mathbf{x}\\ &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\mathbf{x} \\[4pt] &= [0.173 \;-0.072] \begin{bmatrix} 86.090 & -61.487 \\ -61.487 & 90.199 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\[4pt] \hat{y}&= 19.319 \,x_1 - 17.124\,x_2 \end{align} \]

Essa função discriminante linear é a função linear de Fisher, que maximiza a separação entre as duas populações, e a separação máxima nas amostras é

\[ \begin{align} d^2 &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2) \\[6pt] &= [0.173 \;-0.072] \begin{bmatrix} 86.090 & -61.487 \\ -61.487 & 90.199 \end{bmatrix} \begin{bmatrix} 0.173 \\ -0.072 \end{bmatrix} \\[6pt] d^2&= 4.574 \end{align} \] \[\Diamond\]

A solução de Fisher para o problema de separação também pode ser usada para classificar novas observações.

• JW6_Example11.4.R

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
))
Dados <- read.csv("JW6Data/T11-8.dat", header = FALSE, sep = "")
colnames(Dados) <- c("Group", "AHF_act", "AHF_antig")

Dados$Group <- factor(Dados$Group,
                      labels = c("noncarriers", "oblig.carriers"))

# As colunas já estão em log10 no arquivo, então só renomeio:
Dados$X1 <- Dados$AHF_act
Dados$X2 <- Dados$AHF_antig

# Separar grupos
X1mat <- as.matrix(subset(Dados, Group == "noncarriers")[, c("X1", "X2")])
X2mat <- as.matrix(subset(Dados, Group == "oblig.carriers")[, c("X1", "X2")])

n1 <- nrow(X1mat)
n2 <- nrow(X2mat)

# Vetores média amostral
xbar1 <- colMeans(X1mat)
xbar2 <- colMeans(X2mat)

# Matrizes de covariância amostral
S1 <- cov(X1mat)
S2 <- cov(X2mat)

# Matriz de covariância combinada (pooled)
S_pooled <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n1 + n2 - 2)
S_pooled_inv <- solve(S_pooled)

# Vetor discriminante â = (x̄1 - x̄2)' S_pooled^{-1}
a_hat <- t(xbar1 - xbar2) %*% S_pooled_inv   # 1 x 2
a_hat

         X1        X2
[1,] 19.319 -17.12424

d2 <- t(xbar1 - xbar2) %*% S_pooled_inv %*% (xbar1 - xbar2)
d2

        [,1]
[1,] 4.57431

Regra de Alocação Baseada na Função Discriminante de Fisher

Devemos ter \((n_1 + n_2 - 2) \ge p\), caso contrário \(\mathbf{s}_{\text{pooled}}\) é singular, e a inversa usual, \(\mathbf{s}_{\text{pooled}}^{-1}\), não existe.

Aloque \(\mathbf{x}_0\) em \(\pi_1\) se

\[ \begin{align} \hat{y}_0 &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\mathbf{x}_0 \\ &\ge \hat{m} = \dfrac{1}{2} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} (\bar{\mathbf{x}}_1 + \bar{\mathbf{x}}_2) \end{align} \tag{11-25} \]

ou, de modo equivalente,

\[ \hat{y}_0 - \hat{m} \ge 0 \]

Aloque \(\mathbf{x}_0\) em \(\pi_2\) se

\[ \hat{y}_0 < \hat{m} \]

ou, de modo equivalente,

\[ \hat{y}_0 - \hat{m} < 0 \]

O procedimento da equação (11-23) é ilustrado esquematicamente, para \(p=2\), na Figura 11.5. Todos os pontos nos gráficos de dispersão são projetados em uma reta na direção de \(\hat{\mathbf{a}}\), e essa direção é variada até que as amostras estejam maximamente separadas.

A função discriminante linear de Fisher em (11-25) foi desenvolvida sob a suposição de que as duas populações, quaisquer que sejam suas formas, possuem uma matriz de covariância comum. Consequentemente, não é surpreendente que o método de Fisher corresponda a um caso particular da regra de mínima esperança de custo de classificação incorreta. O primeiro termo,

\[ \hat{y} = (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{S}_{\text{pooled}}^{-1}\mathbf{x}, \]

na regra de classificação (11-18), é a função linear obtida por Fisher que maximiza a variabilidade “entre grupos” (between) relativa à variabilidade “dentro dos grupos” (within). [Veja (11-23).]

A expressão completa é

\[ \begin{align} \hat{w} &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1}\mathbf{x} - \dfrac{1}{2} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} (\bar{\mathbf{x}}_1 + \bar{\mathbf{x}}_2)\\ \hat{w} &= (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \mathbf{s}_{\text{pooled}}^{-1} \left( \mathbf{x} - \dfrac{1}{2}(\bar{\mathbf{x}}_1 + \bar{\mathbf{x}}_2) \right) \end{align} \tag{11-26} \]

A quantidade \(\hat{w}\) é frequentemente chamada estatística de classificação de Anderson. Mais uma vez, se

\[ \dfrac{c(1\mid 2)}{c(2\mid 1)} \dfrac{p_2}{p_1} = 1 \]

então

\[ \ln\left( \dfrac{c(1\mid 2)}{c(2\mid 1)} \dfrac{p_2}{p_1} \right) = 0 \]

e a Regra (11-18) é comparável à Regra (11-26), baseada na função discriminante linear de Fisher.

Assim, desde que as duas populações normais tenham a mesma matriz de covariâncias, a regra de classificação de Fisher é equivalente à regra de mínimo ECM com probabilidades a priori iguais e custos iguais de misclassificação.

A Classificação é uma Boa Ideia?

Para duas populações, a separação relativa máxima que pode ser obtida considerando combinações lineares das observações multivariadas é igual à distância \(D^2\).

Isso é conveniente porque \(D^2\) pode ser usado, em certas situações, para testar se as médias populacionais \(\mu_1\) e \(\mu_2\) diferem significativamente.

Consequentemente, um teste para diferenças entre vetores de médias pode ser interpretado como um teste da “significância” da separação que pode ser alcançada.

Suponha que as populações \(\pi_1\) e \(\pi_2\) sejam normais multivariadas com uma matriz de covariância comum \(\Sigma\).

Então, como na Seção 6.3, um teste de \(H_0:\boldsymbol{\mu}_1=\boldsymbol{\mu}_2\) versus \(H_1:\boldsymbol{\mu}_1 \ne \boldsymbol{\mu}_2\) é realizado referindo

\[ \frac{n_1+n_2-p-1}{(n_1+n_2-2)p} \frac{n_1 n_2}{n_1 + n_2} D^2 \sim F_{p,\,n_1 + n_2 - p - 1} \]

Se \(H_0\) é rejeitada, podemos concluir que a separação entre as duas populações \(\pi_1\) e \(\pi_2\) é significante.

Comentário. Separação significativa não implica necessariamente boa classificação. Como veremos na Seção 11.4, a eficácia de um procedimento de classificação pode ser avaliada independentemente de qualquer teste de separação. Por outro lado, se a separação não é significativa, a busca por uma regra de classificação útil provavelmente será infrutífera.

Classificação de populações normais quando \(\boldsymbol{\Sigma}_1 \neq \boldsymbol{\Sigma}_2\)

Como era de se esperar, as regras de classificação tornam-se mais complicadas quando as matrizes de covariância populacionais são desiguais.

Considere as densidades normais multivariadas em (11-10), com \(\boldsymbol{\Sigma}_i\), \(i = 1,2\), substituindo \(\boldsymbol{\Sigma}\). Assim, as matrizes de covariância, bem como os vetores de médias, são diferentes para as duas populações. Como vimos, as regiões de mínimo ECM e de mínima probabilidade total de classificação incorreta (TPM) dependem da razão das densidades, \(f_1(\mathbf{x}) / f_2(\mathbf{x})\), ou, de forma equivalente, do logaritmo natural dessa razão, \(\ln\{ f_1(\mathbf{x}) / f_2(\mathbf{x}) \} = \ln f_1(\mathbf{x}) - \ln f_2(\mathbf{x})\). Quando as densidades normais multivariadas têm estruturas de covariância diferentes, os termos na razão de densidades que envolvem \(\lvert \boldsymbol{\Sigma}_i \rvert^{1/2}\) não se cancelam como ocorre quando \(\boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\). Além disso, as formas quadráticas nos expoentes de \(f_1(\mathbf{x})\) e \(f_2(\mathbf{x})\) não se combinam para produzir a expressão simples obtida em (11-13).

Substituindo densidades normais multivariadas com diferentes matrizes de covariância em (11-6) e, em seguida, tomando logaritmos naturais e simplificando (ver Exercício 11.15), obtemos as regiões de classificação

\[ \begin{align} R_1:\;& -\frac12\,\mathbf{x}^{\prime} \left(\boldsymbol{\Sigma}_1^{-1} - \boldsymbol{\Sigma}_2^{-1}\right)\mathbf{x} + \left( \boldsymbol{\mu}_1^{\prime}\boldsymbol{\Sigma}_1^{-1} - \boldsymbol{\mu}_2^{\prime}\boldsymbol{\Sigma}_2^{-1} \right)\mathbf{x}-k\\ &\ge \ln\!\left[ \frac{c(1\mid 2)}{c(2\mid 1)} \frac{p_2}{p_1} \right] \\[4pt] R_2:\;& -\frac12\,\mathbf{x}^{\prime} \left(\boldsymbol{\Sigma}_1^{-1} - \boldsymbol{\Sigma}_2^{-1}\right)\mathbf{x} + \left( \boldsymbol{\mu}_1^{\prime}\boldsymbol{\Sigma}_1^{-1} - \boldsymbol{\mu}_2^{\prime}\boldsymbol{\Sigma}_2^{-1} \right)\mathbf{x}-k\\ &< \ln\left[ \frac{c(1\mid 2)}{c(2\mid 1)} \frac{p_2}{p_1} \right] \end{align} \tag{11-27} \]

sendo que

\[ \begin{align} k &= \frac12\, \left(\ln\left( \frac{\lvert\boldsymbol{\Sigma}_1\rvert} {\lvert\boldsymbol{\Sigma}_2\rvert} \right) + \left( \boldsymbol{\mu}_1^{\prime}\boldsymbol{\Sigma}_1^{-1}\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2^{\prime}\boldsymbol{\Sigma}_2^{-1}\boldsymbol{\mu}_2 \right)\right) \end{align} \tag{11-28} \]

As regiões de classificação são definidas por funções quadráticas em \(\mathbf{x}\). Quando \(\boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\), o termo quadrático \(-\tfrac12\,\mathbf{x}^{\prime}(\boldsymbol{\Sigma}_1^{-1} - \boldsymbol{\Sigma}_2^{-1})\mathbf{x}\) desaparece, e as regiões definidas por (11-27) se reduzem àquelas dadas por (11-14).

A regra de classificação para populações normais multivariadas gerais decorre diretamente de (11-27).

Resultado 11.4. Sejam as populações \(\pi_1\) e \(\pi_2\) descritas por densidades normais multivariadas com vetores de médias e matrizes de covariância \(\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1\) e \(\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2\), respectivamente.

A regra de alocação que minimiza o custo esperado de classificação incorreta é dada por:

Alocar \(\mathbf{x}_0\) em \(\pi_1\) se

\[ -\frac12\,\mathbf{x}_0^{\prime} \left( \boldsymbol{\Sigma}_1^{-1} - \boldsymbol{\Sigma}_2^{-1} \right) \mathbf{x}_0 + \left( \boldsymbol{\mu}_1^{\prime}\boldsymbol{\Sigma}_1^{-1} - \boldsymbol{\mu}_2^{\prime}\boldsymbol{\Sigma}_2^{-1} \right) \mathbf{x}_0 - k \;\ge\; \ln\!\left[ \frac{c(1\mid 2)}{c(2\mid 1)} \frac{p_2}{p_1} \right] \]

e alocar \(\mathbf{x}_0\) em \(\pi_2\) caso contrário.

O valor de \(k\) encontra-se definido em (11-28).

\[\Diamond\]

Na prática, a regra de classificação em Result 11.5 é implementada substituindo as quantidades amostrais \(\mathbf{\bar{x}}_1\), \(\mathbf{\bar{x}}_2\), \(\mathbf{s}_1\) e \(\mathbf{s}_2\)
(ver (11-16)) para \(\boldsymbol{\mu}_1\), \(\boldsymbol{\mu}_2\), \(\boldsymbol{\Sigma}_1\) e \(\boldsymbol{\Sigma}_2\), respectivamente.

As desigualdades \(n_1 > p\) e \(n_2 > p\) devem ambas ser satisfeitas para que \(\mathbf{s}_1^{-1}\) e \(\mathbf{s}_2^{-1}\) existam. Essas quantidades são usadas no lugar de \(\boldsymbol{\Sigma}_1^{-1}\) e \(\boldsymbol{\Sigma}_2^{-1}\), respectivamente, no análogo amostral (11-29).

A classificação com funções quadráticas é bastante trabalhosa em mais do que duas dimensões e pode levar a alguns resultados estranhos. Isso é particularmente verdadeiro quando os dados não são (essencialmente) normais multivariados.

A Figura 11.6(a) mostra as regras com custos iguais e probabilidades a priori iguais, baseadas no caso idealizado de duas distribuições normais com variâncias diferentes. Esta regra quadrática leva a uma região \(R_1\) consistindo de dois subconjuntos disjuntos de pontos.

Figura 11.6: Regras quadráticas para (a) duas distribuições normais com variâncias desiguais e (b) duas distribuições, sendo uma delas não normal — regra não apropriada.

Em muitas aplicações, a cauda inferior da distribuição \(\pi_1\) será menor do que a prescrita por uma distribuição normal. Nesse caso, como mostrado na Figura 11.6(b), a parte inferior da região \(R_1\), produzida pelo procedimento quadrático, não se ajusta bem às distribuições populacionais e pode levar a taxas de erro elevadas. Um aspecto sério da regra quadrática é que ela é sensível a desvios de normalidade.

Se os dados não forem normais multivariados, duas opções estão disponíveis.

Primeiro, os dados não normais podem ser transformados para ficarem mais próximos da normalidade, e pode-se realizar um teste de igualdade das matrizes de covariância (ver Seção 6.6) para decidir se a regra linear (11-18) ou a regra quadrática (11-29) é apropriada. Transformações são discutidas no Capítulo 4. (Os testes usuais de homogeneidade de covariância são muito afetados pela não normalidade. A transformação dos dados não normais para dados normais deve ser feita antes desses testes.)

Segundo, pode-se usar uma regra linear (ou quadrática) sem se preocupar com a forma das populações geradoras e esperar que ela funcione razoavelmente bem. Estudos (ver [22] e [23]) mostraram, entretanto, que existem casos não normais nos quais uma função discriminante linear tem desempenho ruim, mesmo quando as matrizes de covariância populacionais são iguais. A moral da história é sempre verificar o desempenho de qualquer procedimento de classificação. No mínimo, isso deve ser feito com os conjuntos de dados usados para construir o classificador. Idealmente, deve haver dados suficientes disponíveis para fornecer amostras de “treinamento” e amostras de “validação”. As amostras de treinamento podem ser usadas para desenvolver a função de classificação, e as amostras de validação podem ser usadas para avaliar seu desempenho.

Exemplo: ## Classificação LDA e QDA com iris

Do ponto de vista teórico, a LDA assume covariâncias iguais nos grupos, enquanto a QDA permite matrizes de covariância diferentes.

A seguir, aplicamos LDA e QDA ao arquivo iris.

• JW6_Example11_iris.R

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
))

# Dados: excluir versicolor, ficar só com setosa e virginica
Dados <- subset(datasets::iris,
                Species != "versicolor",
                select = c(Sepal.Length, Sepal.Width, Species))
Dados$Species <- droplevels(Dados$Species)
print(summary(Dados), digits = 2)

  Sepal.Length    Sepal.Width         Species  
 Min.   :4.300   Min.   :2.200   setosa   :50  
 1st Qu.:5.000   1st Qu.:3.000   virginica:50  
 Median :5.700   Median :3.200                 
 Mean   :5.797   Mean   :3.201                 
 3rd Qu.:6.500   3rd Qu.:3.425                 
 Max.   :7.900   Max.   :4.400                 

# Ajuste LDA e QDA 
lda_fit <- MASS::lda(Species ~ Sepal.Length + Sepal.Width,
                     data = Dados)

qda_fit <- MASS::qda(Species ~ Sepal.Length + Sepal.Width,
                     data = Dados)

print(lda_fit)

Call:
lda(Species ~ Sepal.Length + Sepal.Width, data = Dados)

Prior probabilities of groups:
   setosa virginica 
      0.5       0.5 

Group means:
          Sepal.Length Sepal.Width
setosa           5.006       3.428
virginica        6.588       2.974

Coefficients of linear discriminants:
                   LD1
Sepal.Length  2.208596
Sepal.Width  -2.511742

print(qda_fit)

Call:
qda(Species ~ Sepal.Length + Sepal.Width, data = Dados)

Prior probabilities of groups:
   setosa virginica 
      0.5       0.5 

Group means:
          Sepal.Length Sepal.Width
setosa           5.006       3.428
virginica        6.588       2.974

# MATRIZES DE CONFUSÃO 
# Predições no próprio conjunto de dados (erro aparente)

# LDA
pred_lda_all <- predict(lda_fit, newdata = Dados)
tab_lda <- table(Observado = Dados$Species,
                 Predito = pred_lda_all$class)
cat("\nMatriz de confusão LDA (erro aparente):\n")

Matriz de confusão LDA (erro aparente):

print(tab_lda)

           Predito
Observado   setosa virginica
  setosa        50         0
  virginica      1        49

erro_lda <- 1 - sum(diag(tab_lda)) / sum(tab_lda)
cat("Taxa de erro aparente LDA:", round(erro_lda, 4), "\n")

Taxa de erro aparente LDA: 0.01 

# QDA
pred_qda_all <- predict(qda_fit, newdata = Dados)
tab_qda <- table(Observado = Dados$Species,
                 Predito = pred_qda_all$class)
cat("\nMatriz de confusão QDA (erro aparente):\n")

Matriz de confusão QDA (erro aparente):

print(tab_qda)

           Predito
Observado   setosa virginica
  setosa        50         0
  virginica      0        50

erro_qda <- 1 - sum(diag(tab_qda)) / sum(tab_qda)
cat("Taxa de erro aparente QDA:", round(erro_qda, 4), "\n")

Taxa de erro aparente QDA: 0 

# Grade no plano (Sepal.Length, Sepal.Width)
x_seq <- seq(min(Dados$Sepal.Length) - 0.2,
             max(Dados$Sepal.Length) + 0.2,
             length.out = 200)
y_seq <- seq(min(Dados$Sepal.Width) - 0.2,
             max(Dados$Sepal.Width) + 0.2,
             length.out = 200)

grid_xy <- expand.grid(Sepal.Length = x_seq,
                       Sepal.Width  = y_seq)

# Predições na grade para LDA e QDA
pred_lda <- predict(lda_fit, newdata = grid_xy)
pred_qda <- predict(qda_fit, newdata = grid_xy)

# Diferença de probabilidades a posteriori (setosa - virginica)
grid_xy$lda_diff <- pred_lda$posterior[, "setosa"] -
  pred_lda$posterior[, "virginica"]

grid_xy$qda_diff <- pred_qda$posterior[, "setosa"] -
  pred_qda$posterior[, "virginica"]

# Gráfico: fronteiras discriminantes LDA (linear) e QDA (quadrática)
g <- ggplot2::ggplot(
  data = Dados,
  mapping = ggplot2::aes(x = Sepal.Length,
                         y = Sepal.Width,
                         colour = Species,
                         shape  = Species)
) +
  ggplot2::geom_point(size = 2) +
  # Fronteira LDA (linha onde posterior_setosa - posterior_virginica = 0)
  ggplot2::stat_contour(
    data = grid_xy,
    mapping = ggplot2::aes(x = Sepal.Length,
                           y = Sepal.Width,
                           z = lda_diff),
    breaks     = 0,
    linewidth  = 0.7,
    colour     = "black",
    linetype   = 2,
    inherit.aes = FALSE
  ) +
  # Fronteira QDA (linha onde posterior_setosa - posterior_virginica = 0)
  ggplot2::stat_contour(
    data = grid_xy,
    mapping = ggplot2::aes(x = Sepal.Length,
                           y = Sepal.Width,
                           z = qda_diff),
    breaks     = 0,
    linewidth  = 0.7,
    colour     = "red",
    linetype   = 1,
    inherit.aes = FALSE
  ) +
  ggplot2::labs(
    x = "Comprimento da sépala",
    y = "Largura da sépala",
    title = "Fronteiras discriminantes\n linear (LDA, tracejada) e quadrática (QDA)"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(panel.grid = ggplot2::element_blank())

print(g)

Questões

• 11.1, 11.2, 11.6, 11.7, 11.10, 11.11, 11.12, 11.13, 11.14, 11.23, 11.27, 11.28, 11.32, 11.33, 11.34, 11.35

Referências

• Livro-texto:
  • JOHNSON, RA & WICHERN, DW (2007) Applied multivariate statistical analysis. 6^th ed. NJ: Prentice Hall.
  • JOHNSON, RA & WICHERN, DW (2007) Applied multivariate statistical analysis: Solutions Manual. 6^th ed. NJ: Prentice Hall.
• Silveira, PSP & Siqueira, JO (2022) Better to be in agreement than in bad company: A critical analysis of many kappa-like tests. Behavior Research Methods 55(7): 3326-47. DOI: doi: 10.3758/s13428-022-01950-0