A análise de componentes principais (PCA) é amplamente utilizada para reduzir a dimensionalidade de grandes conjuntos de dados. No entanto, ele otimiza implicitamente uma função objetivo que é equivalente a uma verossimilhança gaussiana. Portanto, para dados como contagens discretas não negativas que não seguem a distribuição normal, o PCA pode ser inadequado.

Para remediar isso, Collins (2002) propôs generalizar o PCA para a família exponencial de uma maneira análoga à generalização da regressão linear para modelos lineares generalizados. Esse artigo baseia-se em ideias da família Exponencial, modelos lineares generalizados e distâncias de Bregman, para fornecer uma generalização do PCA para funções de perda que consideramos serem mais adequadas para outros tipos de dados.

Posteriormente fornecemos uma derivação detalhada do PCA generalizado (GLM-PCA) com foco na otimização usando o escore de Fisher. Também expandimos o modelo de Collins incorporando covariáveis e propomos transformações post hoc para melhorar a interpretabilidade de fatores latentes.



1. Introdução

A análise de componentes principais (PCA) é uma técnica de redução de dimensionalidade extremamente popular que tenta encontrar um subespaço de baixa dimensão passando perto de um determinado conjunto de pontos \(x_1,\cdots,x_n\in\mathbb{R}^d\).

Mais especificamente, no PCA, encontramos um subespaço de dimensão inferior que minimiza a soma dos quadrados das distâncias dos pontos de dados \(x\) às suas projeções \(\theta_i\) no subespaço, ou seja, \[\begin{equation} \tag{1} \sum_{i=1}^n || x_i-\theta_i||^2\cdot \end{equation}\]

Isso acaba sendo equivalente a escolher um subespaço que maximize a soma dos comprimentos quadrados das projeções \(\theta_i\), que é o mesmo que a variância empírica dessas projeções, se os dados estiverem centrados na origem, de modo que \(\sum_i x_i=0\).

O PCA também tem outra interpretação conveniente que talvez seja menos conhecida. Nesta interpretação probabilística, cada ponto \(x_i\) é pensado como um sorteio aleatório de alguma distribuição desconhecida \(P_{\theta_i}\), onde \(P_{\theta}\) denota a distribuição gaussiana com variância identidade e média \(\theta\in\mathbb{R}^d\).

O objetivo então do PCA é encontrar o conjunto de parâmetros \(\theta_1,\cdots,\theta_n\) que maximiza a verossimilhança dos dados, sujeito à condição de que todos esses parâmetros estejam em um subespaço de baixa dimensão. Em outras palavras, \(x_1,\cdots,x_n\) são considerados versões corrompidas por ruído de alguns pontos verdadeiros \(\theta_1,\cdots,\theta_n\) que se encontram em um subespaço; o objetivo é encontrar esses pontos verdadeiros e a suposição principal é que o ruído é gaussiano.

A equivalência dessa interpretação com as dadas acima decorre simplesmente do fato de que o negativo do logaritmo da verossimilhança sob esse modelo guassiano é igual, ignorando as constantes, à equação em (1).

Essa suposição gaussiana pode ser inadequada, por exemplo, se os dados tiverem valor binário, valor inteiro ou não forem negativos. Na verdade, a gaussiano é apenas uma das distribuições canônicas que compõem a família exponencial e é uma distribuição adaptada para dados de valores reais. A Poisson é mais adequada para dados inteiros e o Bernoulli para dados binários. Parece natural considerar variantes do PCA baseadas nessas outras distribuições no lugar da Gaussiana.

Estendemos o PCA para o restante da família exponencial. Seja \(\{P_{\theta}\}\) qualquer conjunto paramétrico de distribuições da família exponencial, onde \(\theta\) é o parâmetro natural de uma distribuição. Por exemplo, uma distribuição de Poisson unidimensional pode ser parametrizada por \[ P_\theta(n)=e^{n\theta-e^\theta}/n!=e^{-\lambda}\lambda^n/n!, \qquad n\in\{0,1,2,\cdots\}\cdot \] Dados os dados \(x_1,\cdots,x_n\in\mathbb{R}^d\), o objetivo agora é encontrar parâmetros \(\theta_1,\cdots,\theta_n\) que se encontram em um subespaço de baixa dimensão e para os quais a log-verossimilhança \(\sum_i \log P_{\theta_i}(x_i)\) é maximizada.

Uma abordagem unificada permite sem esforço esquemas de redução de dimensionalidade híbrida em que diferentes tipos de distribuições podem ser usados para diferentes atributos dos dados. Se os dados \(x_i\) tiverem alguns atributos binários e alguns atributos de valor inteiro, algumas coordenadas dos correspondentes \(\theta_i\) podem ser parâmetros de distribuições binomiais, enquanto outras são parâmetros de distribuições de Poisson. No entanto, para simplificar a apresentação, neste resumo assumimos que todas as distribuições são do mesmo tipo.

Os esquemas de redução de dimensionalidade para distribuições não gaussianas são substancialmente diferentes do PCA. Por exemplo, em PCA os parâmetros \(\theta_i\), que são médias de gaussianas, encontram-se em um espaço que coincide com o dos dados \(x_i\). Este não é o caso em geral e, portanto, embora os parâmetros \(\theta_i\) estejam em um subespaço linear, eles normalmente correspondem a uma superfície não linear no espaço dos dados.

A discrepância e interação entre o espaço dos parâmetros \(\theta\) e o espaço dos dados \(x\) é uma preocupação central no estudo de famílias exponenciais, modelos lineares generalizados (GLM’s) e distâncias de Bregman. Nossa exposição é inevitavelmente tecida em torno desses três assuntos intimamente relacionados. Em particular, mostramos que a maneira como generalizamos o PCA é exatamente análoga à maneira como a regressão é generalizada pelos GLMs.

O problema de otimização derivado pode ser resolvido naturalmente por um algoritmo que minimiza alternadamente os componentes da análise e seus coeficientes; assim, o algoritmo é uma reminiscência dos procedimentos alternados de minização de Csiszáar and Tusnády (1984). Neste caso, cada lado da minimização é um programa convexo simples que pode ser interpretado como uma projeção em relação a uma distância de Bregman adequada; no entanto, o programa geral geralmente não é convexo.

No caso de distribuições gaussianas, o algoritmo coincide exatamente com o método de potência para calcular autovetores; neste sentido é uma generalização de um dos algoritmos mais antigos para PCA. Embora omitido por falta de espaço, pode-se mostrar que o procedimento converge em que qualquer ponto limite do espaço paramétrico dos coeficientes é um ponto estacionário da função de perda. Além disso, uma pequena modificação do critério de otimização garante a existência de pelo menos um ponto limite.

Alguns comentários sobre a notação: Todos os vetores neste artigo são vetores linha. Se \(M\) for uma matriz, denotamos sua \(i\)-ésima linha por \(m_i\) e seu \(i,j\)-ésimo elemento por \(m_{ij}\).


2. A família exponencial, MLG’s e distâncias de Bregman


2.1 Família exponencial e modelos lineares generalizados

Na família exponencial de distribuições, a probabilidade condicional de um valor \(x\) determinado valor de parâmetro \(\theta\) assume a seguinte forma: \[\begin{equation} \tag{2} \log P(x|\theta)=\log P_0(x) + x\theta - G(\theta)\cdot \end{equation}\]

Aqui \(\theta\) é o parâmetro natural da distribuição e geralmente pode assumir qualquer valor em reais. \(G(\theta)\) é uma função que garante que a soma ou integral de \(P(x|\theta)\) sobre o domínio de \(x\) é 1. Disto segue que \[ G(\theta)=\log \sum_{x\in\mathcal{X}} P_0(x)e^{x\theta}\cdot \]

Usamos \(\mathcal{X}\) para denotar o domínio de \(x\). A soma é substituída por uma integral no caso contínuo, onde \(P\) define uma densidade sobre \(\mathcal{X}\). \(P_0\) é um termo que depende apenas de \(x\) e geralmente pode ser ignorado como uma constante durante a estimação. A principal diferença entre os diferentes membros da família é a forma de \(G(\theta)\). Veremos que quase todos os conceitos dos algoritmos PCA neste artigo decorrem diretamente da definição de \(G\).

Um primeiro exemplo é uma distribuição normal, com média \(\mu\) e variância unitária, que tem uma densidade que geralmente é escrita como \[ \log P(x|\mu) = -\log\sqrt{2\pi} -(x-\mu)^2/2\cdot \] Pode-se verificar que este é um membro da família exponencial com \(\log P_0(x)=-\log\sqrt{2\pi}-x^2/2\), \(\theta=\mu\) e \(G(\theta)=\theta^2/2\). Outro caso comum é uma distribuição de Bernoulli para o caso de resultados binários. Nesse caso \(\mathcal{X}=\{0,1\}\). A probabilidade de \(x\in\mathcal{X}\) é usualmente escrita como \(P(x|p)=p^x(1-p)^{1-x}\) onde \(p\in[0,1]\) é o parâmetro. Esta distribuição é um membro da família exponencial com \(P_0(x)=1\), \(\theta=\log\big(p/1-p\big)\) e \(G(\theta)=\log\big(1+e^\theta\big)\).

Uma função crítica é a derivada \(G'(\theta)\), que denotaremos como \(g(\theta)\) ao longo deste artigo. Diferenciando \[ G(\theta)=\log_{x\in\mathcal{X}} P_0(x)e^{x\theta}, \] verifica-se que \(g(\theta)=\mbox{E}\big(x|\theta\big)\), a esperança de \(x\) sob \(P(x|\theta)\).

Na distribuição normal, \(\mbox{E}\big(x|\theta\big)=\mu\) e no caso de Bernoulli \(\mbox{E}\big(x|\theta\big)=p\). No caso geral, \(\mbox{E}\big(x|\theta\big)\) é referido como o “parâmetro de esperança” e \(g\) define uma função dos valores de parâmetros naturais para os valores de parâmetros na esperança.

A generalização do PCA é análoga à maneira pela qual os modelos lineares generalizados (GLMs) (McCullagh and J. A. Nelder, 1990) fornecem um tratamento unificado de regressão para a família exponencial, generalizando a regressão de mínimos quadrados para funções de perda que são mais apropriadas para outros membros dessa família. A configuração da regressão assume uma amostra de treinamento de pares \((x_i,y_i)\), onde \(x_i\in\mathbb{R}^d\) é \(y_i\in\mathbb{R}^d\) um vetor de atributos e é alguma variável de resposta.

Os parâmetros do modelo são um vetor \(\beta\in\mathbb{R}^d\). O produto escalar \(\beta\cdot x_i\) é considerado uma aproximação de \(y_i\). Na regressão de mínimos quadrados, os parâmetros ótimos \(\beta^*\) são definidos como \[ \beta^* =\arg\min_{\beta\in\mathbb{R}^d} \sum_i \big(y_i-\beta_i\cdot x_i \big)^2\cdot \]

Nos GLM’s, \(h(\beta\cdot x_i)\) é assumido para aproximar o parâmetro da esxperança do modelo exponencial, onde \(h\), é o inverso da “função de ligação”. Uma escolha natural é usar a “ligação canônica”, onde \(h=g\), \(g\) sendo a derivada \(G'(\theta)\). Neste caso, os parâmetros naturais são diretamente aproximados por \(\beta\cdot x_i\) e a log-verossimilhança \(\sum_i \log P(y_i|x_i)\) é simplesmente \[ \sum_i \log P_0(y_i)+y_i \beta\cdot x_i -G(\beta\cdot x_i)\cdot \]

No caso de uma distribuição normal com variância fixa, \(G(\theta)=\theta^2/2\) segue-se que o critério de máxima verossimilhança é equivalente ao critério de mínimos quadrados. Outro caso interessante é a regressão logística em que \(G(\theta)=\log\big(1+e^\theta\big)\) e o negativo da log-verossimilhança para os parâmetros \(\beta\) é \[ \sum_i \log\left( 1+e^{-y_i^*\beta\cdot x_i}\right), \] onde \(y_i^*=1\) se \(y_1=1\), \(y_i^*=-1\) se \(y_i=0\).


2.2 Distâncias de Bregman e a família exponencial

Seja \(F: \Delta\to \mathbb{R}\) ser uma função diferenciável e estritamente convexa definida em um conjunto convexo fechado \(\Delta\subseteq \mathbb{R}\). A distância de Bregman associada a \(F\) é definido para \(p,q\in\Delta\) como \[ B_F(p \, || \, q) = P(p)-F(q)-f(p)(p-q), \] onse \(f(x)=F'(x)\). Pode-se mostrar que, em geral, toda distância de Bregman é não negativa e é igual a zero se, e somente se, seus dois argumentos são iguais.

Para a família exponencial, a log-verossimilhança \(\log P(x|\theta)\) está diretamente relacionado a uma distância de Bregman. Especificamente, Katy S. Azoury and M. K. Warmuth (2001) e Jürgen Forster and Manfred Warmuth (2002) definem uma função “dual” \(F\) através de \(G\) e \(g\): \[\begin{equation} \tag{3} F\big( g(\theta)\big)+G(\theta)=g(\theta) \theta\cdot \end{equation}\]

Pode-se mostrar sob condições bastante gerais que \(f(x)=g^{-1}(x)\). A aplicação dessas identidades implica que o negativo da log-verossimilhança de um ponto pode ser expresso por meio de uma distância de Bregman: \[\begin{equation} \tag{4} -\log P(x|\theta) =-\log P_0(x) -F(x)+B_F(x \, || \, g(\theta))\cdot \end{equation}\]

Em outras palavras, o negativo da log-verossimilhança sempre pode ser escrita como uma distância de Bregman mais um termo constante em relação a \(\theta\) e que, portanto, pode ser ignorado. A Tabela 1 resume várias funções de interesse para exemplos da família exponencial.

\[ \begin{array}{cccc} & \mbox{Normal} & \mbox{Bernoulli} & \mbox{Poisson} \\\hline \mathcal{X} & \mathbb{R} & \{0,1\} & \{0,1,2,\cdots\} \\ G(\theta) & \theta^2/2 & \log\big(1+e^\theta\big) & e^\theta \\ g(\theta) & \theta & e^\theta/\big(1+e^\theta\big) & e^\theta \\ F(x) & x^2/2 & x\log(x)+(1-x)\log(1-x) & x\log(x)-x \\ f(x)=g^{-1}(x) & x & \log(x/(1-x)) & \log(x) \\ B_F(p \, || \, q) & (p-q)^2/2 & p\log(p/q)+(1-p)\log\big((1-p)/(1-q)\big) & p\log\big( p/q\big)+q-p \\ B_F(x \, || \, g(\theta)) & (x-\theta)^2/2 & \log\big(1+e^{x^*\cdot \beta}\big), \mbox{ onde } x^*=2x-1 & e^\theta-x\theta+x\log(x)-x \\\hline \end{array} \]

Tabela 1: Várias funções de interesse para três membros da família exponencial.

Acharemos útil estender a ideia de distâncias de Bregman para divergências entre vetores e matrizes. Se \(x\) e \(y\) são vetores e \(A\) e \(B\) são matrizes, sobrecarregamos a notação como \[ B_F(x \, || \, y)=\sum_i B_F(x_i \, || \, y_i) \] e \[ B_F(A \, || \, B)=\sum_i \sum_j B_F(a_{ij} \, || \, b_{ij})\cdot \]

A noção de distância de Bregman, bem como nossa generalização de PCA, podem ser estendidas a vetores de maneira mais geral; aqui, para simplificar, restringimos nossa atenção às distâncias de Bregman e aos problemas PCA dessa forma particular.


3. Análise de componentes principais para a família exponencial

Agora generalizamos o PCA para outros membros da família exponencial. Queremos encontrar \(\theta_i\)’s que estejam próximos dos \(x\)’s e que pertençam a um subespaço de dimensão inferior do espaço de parâmetros. Assim, nossa abordagem é encontrar uma base \(v_1,\cdots,v_\ell\in\mathbb{R}^d\) e representar cada \(\theta_i\) como a combinação linear desses elementos \[ \theta_i = \sum_k a_{ik} v_k, \] que esteja próximo de \(x_i\).

Seja \(X\) uma matriz \(n\times d\) com \(i\)-ésima linha \(x_i\). Seja \(V\) uma matriz \(\ell\times d\) de \(k\)-ésima linha \(v_k\) e seja \(A\) uma matriz \(n\times\ell\) de elementos \(a_{ik}\). Então \(\Theta=AV\) é a matriz \(n\times d\) de linhas \(\theta_i\) como acima. Esta é a matriz cujos valores são os parâmetros naturais que definem a probabilidade de cada ponto em \(\mathcal{X}\).

Seguindo a discussão na Seção 2, consideramos a função de perda assumindo a forma \[ L(V,A)=-\log P(X|A,V)=-\sum_i \sum_j \log P(x_{ij}|\theta_{ij})=C+\sum_i\sum_j \big( -x_{ij}\theta_{ij}+G(\theta_{ij})\big), \] onde \(C\) é um termo constante que será descartado daqui em diante.

A função de perda varia dependendo de qual membro da família exponencial é escolhido, o que simplesmente muda a forma de \(G\). Por exemplo, se \(X\) é uma matriz de valores reais e a distribuição normal é apropriada para os dados, então \(G=\theta²/2\) e o critério de perda é a perda quadrada usual para PCA. Para a distribuição de Bernoulli \(G(\theta)=\log (1+e^\theta)\). Se definirmos \(x_{ij}^*=2x_{ij}-1\), então \(L(V,A)=\sum_i\sum_j \log\big( 1+e^{-x_{ij}^*\theta_{ij}}\big)\).

A partir da relação entre log-verossimilhança e distâncias de Bregman, ver Eq. (4), a perda também pode ser escrita como \[ L(V,A)=\sum_i \sum_j B_F \big(x_{ij} \, || \, g(\theta_{ij})\big) = \sum_i B_F \big(x_i \, || \, g(\theta_i) \big), \] onde permitimos que \(g\) seja aplicado a vetores e matrizes de maneira pontual.

Uma vez que \(V\) e \(A\) foram encontrados para os pontos de dados, o \(i\)-ésimo ponto de dados \(x_i\in\mathbb{R}^d\) pode ser representado como o vetor \(a_i\) no espaço de dimensão inferior \(\mathbb{R}^\ell\).

Então \(a_i\) são os coeficientes que definem uma projeção de Bregman do vetor \(x_i\): \[\begin{equation} \tag{5} a_i = \arg\min_{a\in\mathbb{R}^\ell} B_F \big( x_i \, || \, g(aV)\big)\cdot \end{equation}\] A forma generalizada de PCA também pode ser considerada como a busca por uma base de baixa dimensão, matriz \(V\), que define uma superfície próxima a todos os pontos de dados \(x_i\). Definimos o conjunto de pontos \(\mathcal{Q}(V)\) como sendo \[ \mathcal{Q}(V)= \left\{ g(aV) \, | \, a\in\mathbb{R}^\ell\right\}\cdot \]

O valor ideal para \(V\) então minimiza a soma das distâncias de projeção: \[ V^*=\arg \min_V \sum_i \min_{q\in\mathcal{Q}(V)} B_F(x_i \, | \, q)\cdot \] Observe que para a distribuição normal \(g(\theta)=\theta\) e a distância de Bregman é a distância euclidiana, de modo que a operação de projeção na Eq. (5) é uma projeção linear simples \(a_i=x_i V^\top\). \(\mathcal{Q}(V)\) também é simplificado no caso normal, sendo simplesmente o hiperplano cuja base é \(V\).

Para resumir, uma vez membro da família exponencial e por implicação, escolhendo \(G(\theta)\) uma função convexa, o PCA regular é generalizado da seguinte maneira:

  1. A função de perda é \(-\log P(x|\theta)=-x\theta+G(\theta)+C\), sendo \(C\) uma constante.
  2. A matriz \(\Theta=AV\) é considerada uma matriz de valores de parâmetros naturais.
  3. A derivada \(g(\theta)\) de \(G(\theta)\) define uma matriz de esperanças dos parâmetros \(g(AV)\).
  4. Uma função \(F\) é derivado de \(G\) e \(g\). Uma distância de Bregman \(B_F\) é derivada de \(F\).
  5. A perda é uma soma das distâncias de Bregman dos elementos \(x_{ij}\) para os valores \(g(\theta_{ij})\).
  6. O PCA também pode ser pensado como uma busca por uma matriz \(V\) que define a superfície \(\mathcal{Q}(V)\) que está “próxima” de todos os pontos de dados.

A distribuição normal é um caso simples porque \(g(\theta)=\theta\) e a divergência é a distância euclidiana. A operação de projeção é uma operação linear e \(\mathcal{Q}(V)\) é o hiperplano que tem \(V\) como base.


4. Algoritmos genéricos para minimizar a função de perda

Descrevemos agora um algoritmo genérico para minimização da função de perda. Primeiro, nos concentramos no caso mais simples em que há apenas um único componente, de modo que \(\ell=1\). Eliminamos o subscrito \(k\) de \(a_{ik}\) e \(v_{jk}\). O método é iterativo, com uma escolha inicial aleatória para o valor de \(V\). Sejam \(V^{(t)}\), \(A^{(t)}\), etc. denotando os valores na \(t\)-ésima iteração, e seja \(V^{(0)}\) a escolha aleatória inicial.

Collins et. all. (2001) proporam as atualizações iterativas \[ A^{(t)}=\arg\min_A L(V^{(t-1)},A) \qquad \mbox{e} \qquad V^{(t)}=\arg\min_V L(V,A^{(t)})\cdot \] Assim, \(L\) é minimizado alternadamente com respecto aos dois argumentos, cada vez otimizando um argumento enquanto mantém o outro fixo, reminiscente dos procedimentos alternados de minização de Csiszár and Tusnády (1984).

É útil escrever esses problemas de minimização da seguinte forma: \[ \mbox{Para } i=1,\cdots,n, \qquad a_i^{(t)}=\arg\min_{a\in\mathbb{R}}\sum_j B_F\Big( x_{ij} || g(av_j^{(t-1)})\Big)\\ \mbox{Para } j=1,\cdots,d, \qquad v_j^{(t)}=\arg\min_{v\in\mathbb{R}}\sum_i B_F\Big(x_{ij} || g\big(a^{(t)}v\big) \Big)\\ \]

Podemos então ver que existem \(n+d\) problemas de otimização e que cada um é essencialmente idêntico a um problema de regressão GLM muito simples, onde há um único parâmetro sendo otimizado. Esses subproblemas são facilmente resolvidos, pois as funções são convexas no argumento que está sendo otimizado e a ampla literatura sobre estimação de máxima verossimilhança em GLMs pode ser aplicada diretamente ao problema.

Essas atualizações assumem forma simples para a distribuição normal: \[ A^{(t)}=X\big(V^{(t-1)}\big)^\top\big/|| V^{(t-1)}||^2 \qquad \mbox{e} \qquad V^{(t)}=\big(A^{(t)}\big)^\top X\big/ ||A^{(t)}||^2\cdot \] Segue que \(V^{(t)}=V^{(t-1)} X^\top X/C\), onde \(C\) é um escalar. O método é então equivalente ao método de potência, ver Jolliffe (1986), para encontrar o autovetor de \(X^\top X\) com o maior autovalor, que é a melhor solução de componente único para \(V\). Assim, o algoritmo genérico generaliza um dos algoritmos mais antigos para resolver o problema PCA regular.

A função perda é convexa em qualquer um de seus argumentos com o outro fixo, mas em geral não é convexa nos dois argumentos juntos. Isso torna muito difícil provar a convergência para o mínimo global. A distribuição normal é um caso especial interessante a esse respeito, o método da potência é conhecido por convergir para a solução ótima, apesar da natureza não convexa da superfície de perda. Uma prova simples disso vem das propriedades dos autovetores (Jolliffe, 1986). Também pode ser explicado pela análise do Hessiano \(H\): para qualquer ponto estacionário que não é mínimo global, \(H\) não é semi-definido positivo. Assim, esses pontos estacionários são pontos de sela em vez de mínimos locais.

O Hessiano para a função de perda generalizada é mais complexo; permanece um problema em aberto se também não é semidefinido positivo em pontos estacionários diferentes do mínimo global. Também está em aberto para determinar sob quais condições este algoritmo genérico irá convergir para um mínimo global. Em estudos numéricos preliminares, o algoritmo se comportou bem nesse aspecto. Além disso, qualquer ponto limite da sequência \(\Theta^{(t)}=A^{(t)}V^{(t)}\) será um ponto estacionário.

No entanto, é possível que essa sequência diverja, pois o ótimo pode estar no infinito. Para evitar tais escolhas degeneradas de \(\Theta\), podemos usar uma perda modificada \[ \sum_i \sum_j \left( B_F\big( x_{ij} \, || \, g(\theta_{ij})\big)+\epsilon B_F\big(\mu_0 \, || \, g(\theta_{ij})\big) \right), \] onde \(\epsilon\) é uma pequena constante positiva e \(\mu_0\) é qualquer valor no intervalo de \(g\) e, portanto, para o qual \(g^{-1}(\mu_0)\) é finito.

Isso é aproximadamente equivalente a adicionar um conjugado a priori e encontrar a solução máxima a posteriori. Pode-se provar, para esta perda modificada, que a sequência \(\Theta^{(t)}\) permanece em uma região limitada e, portanto, sempre tem pelo menos um ponto limite que deve ser um ponto estacionário.


5. Exemplo

Utilizaremos uma base da dados no UCI Machine Learning Repository do Center for Machine Learning and Intelligent Systems. O UCI Machine Learning Repository é uma coleção de bancos de dados, domínio de teorias e geradores de dados que são usados pela comunidade de aprendizado de máquina para a análise empírica de algoritmos de aprendizado de máquina. O arquivo foi criado como um arquivo ftp em 1987 por David Aha e outros alunos de pós-graduação da University of California Irvine. Desde então, tem sido amplamente utilizado por estudantes, educadores e pesquisadores em todo o mundo como fonte primária de conjuntos de dados de aprendizado de máquina.

Como indicação do impacto do arquivo, ele foi citado mais de 1.000 vezes, tornando-se um dos 100 “artigos” mais citados em toda a ciência da computação. A versão atual do site foi projetada em 2007 por Arthur Asuncion e David Newman e este projeto é em colaboração com Rexa.info da Universidade de Massachusetts Amherst. O apoio financeiro da National Science Foundation é reconhecido com gratidão.

Muitas pessoas merecem agradecimentos por tornar o repositório um sucesso. Em primeiro lugar entre eles estão os doadores e criadores dos bancos de dados e geradores de dados. Agradecimentos especiais também vão para os antigos bibliotecários do repositório: David Aha, Patrick Murphy, Christopher Merz, Eamonn Keogh, Cathy Blake, Seth Hettich e David Newman. 

parkinsons=read.csv("http://archive.ics.uci.edu/ml/machine-learning-databases/parkinsons/parkinsons.data", 
                      header = TRUE)
names(parkinsons)
##  [1] "name"             "MDVP.Fo.Hz."      "MDVP.Fhi.Hz."     "MDVP.Flo.Hz."    
##  [5] "MDVP.Jitter..."   "MDVP.Jitter.Abs." "MDVP.RAP"         "MDVP.PPQ"        
##  [9] "Jitter.DDP"       "MDVP.Shimmer"     "MDVP.Shimmer.dB." "Shimmer.APQ3"    
## [13] "Shimmer.APQ5"     "MDVP.APQ"         "Shimmer.DDA"      "NHR"             
## [17] "HNR"              "status"           "RPDE"             "DFA"             
## [21] "spread1"          "spread2"          "D2"               "PPE"


Em particular o conjunto de dados apareceu no artigo ‘Exploiting Nonlinear Recurrence and Fractal Scaling Properties for Voice Disorder Detection’, Little MA, McSharry PE, Roberts SJ, Costello DAE, Moroz IM. BioMedical Engineering OnLine 2007, 6:23 (26 June 2007).

As variáveis consideradas são as seguintes:

  1. name - nome e número da gravação
  2. MDVP.Fo.Hz. - Frequência vocal fundamental média
  3. MDVP.Fhi.Hz. - Frequência vocal fundamental máxima
  4. MDVP.Flo.Hz. - Frequência vocal fundamental mínima
  5. MDVP.Jitter…
    MDVP.Jitter.Abs.
    MDVP.RAP
    MDVP.PPQ
    Jitter.DDP - Várias medidas de variação na frequência fundamental
  6. MDVP.Shimmer
    MDVP.Shimmer.dB.
    Shimmer.APQ3
    Shimmer.APQ5
    MDVP.APQ
    Shimmer.DDA - Várias medidas de variação em amplitude
  7. NHR
    HNR - Duas medidas de proporção de ruído para componentes tonais na voz
  8. status - Estado de saúde do sujeito (um) - Parkinson, (zero) - saudável
  9. RPDE
    D2 - Duas medidas de complexidade dinâmica não linear
  10. DFA - Expoente de escala fractal de sinal
  11. spread1
    spread2
    PPE - Três medidas não lineares de variação de frequência fundamental

Vamos utilizar o gráfico das coordenadas paralelas para tentarmos identificar quais variáveis eventualmente exercem maior influência na resposta. Gostariamos de identificar quais, dentre estas variáveis explicativas, identificam os indivíduos como saudáveis (zero) ou com Parkinsons (um). 

library(plotly)
fig <- parkinsons %>% plot_ly(type = 'parcoords',
          line = list(color = ~ status, colorscale = list(c(0,'cyan'),c(1,'black'))),
          dimensions = list(
            list(range = c(80,270), label = '', values = ~ MDVP.Fo.Hz.),
            list(range = c(100,600), label = '', values = ~ MDVP.Fhi.Hz.),
            list(range = c(60,240), label = '', values = ~ MDVP.Flo.Hz.),
            list(range = c(0,0.035), label = '', values = ~ MDVP.Jitter...),
            list(range = c(0,0.0003), label = '', values = ~ MDVP.Jitter.Abs.),
            list(range = c(0,0.025), label = '', values = ~ MDVP.RAP),
            list(range = c(0,0.02), label = '', values = ~ MDVP.PPQ),
            list(range = c(0,0.065), label = '', values = ~ Jitter.DDP),
            list(range = c(0,0.12), label = '', values = ~ MDVP.Shimmer),
            list(range = c(0,1.35), label = '', values = ~ MDVP.Shimmer.dB.),
            list(range = c(0,0.06), label = '', values = ~ Shimmer.APQ3),
            list(range = c(0,0.08), label = '', values = ~ Shimmer.APQ5),
            list(range = c(0,0.14), label = '', values = ~ MDVP.APQ),
            list(range = c(0,0.17), label = '', values = ~ Shimmer.DDA),
            list(range = c(0,0.32), label = '', values = ~ NHR),
            list(range = c(8,34), label = '', values = ~ HNR),
            list(range = c(0.2,0.7), label = '', values = ~ RPDE),
            list(range = c(1.4,3.7), label = '', values = ~ D2),
            list(range = c(0.5,0.85), label = '', values = ~ DFA),
            list(range = c(-8,-2), label = '', values = ~ spread1),
            list(range = c(0,0.5), label = '', values = ~ spread2),
            list(range = c(0,0.55), label = '', values = ~ PPE)
            )
          )
fig

Na figura, as linhas representam os indivíduos e em preto os diagnosticados com Parkinsons. Podemos observar sobreposição entre as linhas, não sendo possível identificar um padrão claro de distinção enre indivíduos sadios e doentes. 

set.seed(150) 
N = nrow(parkinsons)
N
## [1] 195
index = sample(1:N, 0.7*N) 
table(parkinsons$status)
## 
##   0   1 
##  48 147

A tabela acima mostra que nosso conjunto de dados contêm 147 (75%) indivíduos diagnosticados com Parkinsons e 48 (25%) individuos sadios; não sendo considerada uma base de dados muito desbalanceada. Identifiquemos agora as componentes principais. 

training = parkinsons[index,] # cria conjunto de dados de treino 
testing  <- parkinsons[-index,]
Y = training$status
nn = table(Y)
nn
## Y
##   0   1 
##  33 103
X = as.matrix(training[,c(2:17,19:24)])

Nas linas acima primeiro, com a função set.seed fixamos a semente de geração de dados simulados para criarmos um arquivo de dados de treino training, com 70% dos dados e um outro de teste testing, com os restantes 30% dos dados. Ainda criamos o vetor Y, com as respostas selecionadas no arquivo de treino e a matriz X com as covariáveis ou variáveis explicativas, também com somente com as informações selecionadas no arquivo de treino. 

scaled.X = scale(X)
scaled.dados = data.frame(Y,scaled.X)
names(scaled.dados)
##  [1] "Y"                "MDVP.Fo.Hz."      "MDVP.Fhi.Hz."     "MDVP.Flo.Hz."    
##  [5] "MDVP.Jitter..."   "MDVP.Jitter.Abs." "MDVP.RAP"         "MDVP.PPQ"        
##  [9] "Jitter.DDP"       "MDVP.Shimmer"     "MDVP.Shimmer.dB." "Shimmer.APQ3"    
## [13] "Shimmer.APQ5"     "MDVP.APQ"         "Shimmer.DDA"      "NHR"             
## [17] "HNR"              "RPDE"             "DFA"              "spread1"         
## [21] "spread2"          "D2"               "PPE"
library(plotly)
scaled.fig <- scaled.dados %>% plot_ly(type = 'parcoords',
          line = list(color = ~ Y, colorscale = list(c(0,'cyan'),c(1,'black'))),
          dimensions = list(
            list(label = '', values = ~ MDVP.Fo.Hz.),
            list(label = '', values = ~ MDVP.Fhi.Hz.),
            list(label = '', values = ~ MDVP.Flo.Hz.),
            list(label = '', values = ~ MDVP.Jitter...),
            list(label = '', values = ~ MDVP.Jitter.Abs.),
            list(label = '', values = ~ MDVP.RAP),
            list(label = '', values = ~ MDVP.PPQ),
            list(label = '', values = ~ Jitter.DDP),
            list(label = '', values = ~ MDVP.Shimmer),
            list(label = '', values = ~ MDVP.Shimmer.dB.),
            list(label = '', values = ~ Shimmer.APQ3),
            list(label = '', values = ~ Shimmer.APQ5),
            list(label = '', values = ~ MDVP.APQ),
            list(label = '', values = ~ Shimmer.DDA),
            list(label = '', values = ~ NHR),
            list(label = '', values = ~ HNR),
            list(label = '', values = ~ RPDE),
            list(label = '', values = ~ D2),
            list(label = '', values = ~ DFA),
            list(label = '', values = ~ spread1),
            list(label = '', values = ~ spread2),
            list(label = '', values = ~ PPE)
            )
          )
scaled.fig
correla = cor(X)
library(corrplot)
par(cex = 0.7)
corrplot(correla[1:6,], method = "number", type = "upper")

corrplot(correla[7:12,], method = "number", type = "upper")

corrplot(correla[13:18,], method = "number", type = "upper")

corrplot(correla[19:22,], method = "number", type = "upper")

scaled.dados1 = scaled.dados[,-c(3,4,19)]
names(scaled.dados1)
##  [1] "Y"                "MDVP.Fo.Hz."      "MDVP.Jitter..."   "MDVP.Jitter.Abs."
##  [5] "MDVP.RAP"         "MDVP.PPQ"         "Jitter.DDP"       "MDVP.Shimmer"    
##  [9] "MDVP.Shimmer.dB." "Shimmer.APQ3"     "Shimmer.APQ5"     "MDVP.APQ"        
## [13] "Shimmer.DDA"      "NHR"              "HNR"              "RPDE"            
## [17] "spread1"          "spread2"          "D2"               "PPE"
library(Compositional)
Y = as.matrix(scaled.dados1[,1])
scaled.X = as.matrix(scaled.dados1[,-1])
ajuste1 = glm.pcr(Y,scaled.X,k=5)
pred0 = scaled.X%*%ajuste1$vec
PC1 = pred0[,1]; PC2 = pred0[,2];  PC3 = pred0[,3];  PC4 = pred0[,4];  PC5 = pred0[,5]

Utilizamos a função glm.pcr, no pacote Compositional para encontrarmos as 5 componentes principais que eventualemnte podem explicar a resposta. O vetor Y e a matriz X foram criados acima por exigência do comando glm.pcr.

Como resultado do ajuste anteior, criamos as novas variáveis componentes principais PC1, PC2, PC3, PC4 e PC2, as quais serão utilizadas como novas variáveis explicativas das resposta. 

fig1 <- scaled.dados %>% plot_ly(type = 'parcoords',
          line = list(color = ~ Y, colorscale = list(c(0,'cyan'),c(1,'black'))),
          dimensions = list(
            list(label = '', values = ~ PC1),
            list(label = '', values = ~ PC2),
            list(label = '', values = ~ PC3),
            list(label = '', values = ~ PC4),
            list(label = '', values = ~ PC5)
            )
          )
fig1

Ainda não e possível distinguir alguma componente como identificadora da saúde do indivíduo. 

total.scaled = nn[1]+nn[2]
pesos = 1:total.scaled
pesos = ifelse(scaled.dados$Y == 1,nn[1],nn[2])/total.scaled 
dados = data.frame(Y,PC1,PC2,PC3,PC4,PC5,V3=scaled.dados[,3],V4=scaled.dados[,4],V19=scaled.dados[,19])
ajuste2 = glm(Y ~ 0 + ., data = dados, family = binomial("logit"), weights = pesos)
ajuste3 = update(ajuste2, family = binomial("probit"))
ajuste4 = update(ajuste2, family = binomial("cloglog"))
ajuste5 = update(ajuste2, family = binomial("cauchit"))
AIC(ajuste2); AIC(ajuste3); AIC(ajuste4); AIC(ajuste5)
## [1] 29.04672
## [1] 28.98579
## [1] 31.67049
## [1] 29.891

Escolhemos como mais adequado o modelo no objeto ajuste3, ou seja, escolhemos como modelo mais adequado (menor AIC) o modelo de regressão logística com ligação complementar log-log. 

summary(ajuste3,correlation = TRUE)
## 
## Call:
## glm(formula = Y ~ 0 + ., family = binomial("probit"), data = dados, 
##     weights = pesos)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.03051   0.02417   0.35077   0.63141   1.10189  
## 
## Coefficients:
##     Estimate Std. Error z value Pr(>|z|)   
## PC1 -0.24958    0.09226  -2.705  0.00683 **
## PC2  0.35543    0.22032   1.613  0.10668   
## PC3 -0.25698    0.23146  -1.110  0.26689   
## PC4  0.28300    0.25140   1.126  0.26028   
## PC5  0.20633    0.40313   0.512  0.60878   
## V3  -0.08856    0.24465  -0.362  0.71737   
## V4  -0.27862    0.32088  -0.868  0.38524   
## V19  0.14647    0.36733   0.399  0.69009   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 69.294  on 136  degrees of freedom
## Residual deviance: 41.659  on 128  degrees of freedom
## AIC: 28.986
## 
## Number of Fisher Scoring iterations: 7
## 
## Correlation of Coefficients:
##     PC1   PC2   PC3   PC4   PC5   V3    V4   
## PC2 -0.19                                    
## PC3 -0.05 -0.01                              
## PC4 -0.04 -0.03 -0.13                        
## PC5 -0.25  0.20  0.24 -0.19                  
## V3   0.10  0.14  0.16 -0.25  0.01            
## V4  -0.20  0.34  0.20 -0.15  0.10  0.00      
## V19  0.19 -0.23 -0.35  0.30 -0.68  0.10 -0.21
library(glmulti)
# Using AICc as criterion. Could use crit = "bic" or "aic" instead.
# Using "good ~ ." to include all variables from data (other than "good")
ajuste7 <- glmulti(y = Y ~ 0 + ., data = dados, fitfunction = "glm", plotty = FALSE, weights = pesos, 
             level = 1, method = "h", crit = "aicc", family = binomial(link = "cloglog"))
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
## 
## After 50 models:
## Best model: Y~1+PC1
## Crit= 28.7343822088385
## Mean crit= 39.3308853393314
## 
## After 100 models:
## Best model: Y~1+PC1
## Crit= 28.7343822088385
## Mean crit= 39.4021855678481
## 
## After 150 models:
## Best model: Y~1+PC1
## Crit= 28.7343822088385
## Mean crit= 35.8343317171347
## 
## After 200 models:
## Best model: Y~1+PC1
## Crit= 28.7343822088385
## Mean crit= 33.4163145020537
## 
## After 250 models:
## Best model: Y~1+PC1
## Crit= 28.7343822088385
## Mean crit= 32.9310210432757
## Completed.
summary(ajuste7)$bestmodel
## [1] "Y ~ 1 + PC1"

Observemos que, pela regressão de todos os subconjuntos, somente apresenta-se significativamente influente na resposta a variável PC1 considerando ainda o intercepto no modelo. 

ajuste8 = glm(Y ~ PC1, data = dados, family = binomial("cloglog"), weights = pesos)
summary(ajuste8)
## 
## Call:
## glm(formula = Y ~ PC1, family = binomial("cloglog"), data = dados, 
##     weights = pesos)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.60728   0.00000   0.07193   0.57629   1.11190  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.5545     0.3603   1.539 0.123790    
## PC1          -0.7054     0.1935  -3.646 0.000266 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 69.294  on 135  degrees of freedom
## Residual deviance: 41.469  on 134  degrees of freedom
## AIC: 28.644
## 
## Number of Fisher Scoring iterations: 8

Vejamos agora o resultado do ajuste deste modelo. 

par(mar=c(4,4,3,1), pch=19)
min(ajuste8$linear.predictors)
## [1] -3.20384
max(ajuste8$linear.predictors)
## [1] 13.30574
preditor = seq(-1,2,by=0.01)
est = family(ajuste8)$linkinv(preditor)
plot(preditor, est, xlab = expression(paste("Preditor linear ", eta)), 
     ylab = "Probabilidades estimadas", type = "l", ylim = c(0,1))
points(ajuste8$linear.predictor[scaled.dados$Y==1],scaled.dados$Y[scaled.dados$Y==1], col = "red")
text(0,0.9,"Pacientes com Parkinsons", col = "red")
points(ajuste8$linear.predictor[dados$Y==0],dados$Y[dados$Y==0], col = "black")
text(0,0.1,"Pacientes sadios")
grid()

A área sob a curva ROC, chamado de AUC é uma medida de desempenho para os problemas de classificação em várias configurações. Curva ROC é uma curva de probabilidade e AUC representa o grau ou medida de separabilidade. Diz o quanto o modelo é capaz de distinguir entre as classes. Quanto maior a AUC, melhor será o modelo em prever 0 classes como 0 e 1 classes como 1. Por analogia, quanto maior a AUC, melhor será a capacidade do modelo em distinguir entre pacientes com a doença e sem doença. 

library(ROCit)
class <- ajuste8$y
score <- ajuste8$fitted.values
roc_empirical <- rocit(score = score, class = class)
roc_empirical$AUC
## [1] 0.8855546
plot(roc_empirical, type = "l")

Sensibilidade (Sensitivity): mede a proporção de respostas positivas que são corretamente identificadas como positivas pelo classificador. Em outras palavras, é a taxa de verdadeiro positivo e pode ser calculada diretamente das entradas da matriz de confusão. Por sua vez, a especificidade (Specificity): a especificidade mede a proporção das respostas negativas que são corretamente identificadas como negativas pelo classificador. Em outras palavras, é a verdadeira taxa negativa e pode ser calculada diretamente das entradas da matriz de confusão. 

preditos = ifelse(family(ajuste8)$linkinv(ajuste8$fitted.values)>0.8,1,0)
x = table(ajuste8$y,preditos)
x
##    preditos
##      0  1
##   0 28  5
##   1 32 71
prop.table(x, margin = 1)
##    preditos
##             0         1
##   0 0.8484848 0.1515152
##   1 0.3106796 0.6893204

Percebemos que nosso modelo, quando considerada resposta 1 para a probabilidade estimada maior do que 0.8, consegue prever 68% dos indivídos com a doença de Parkinsons, errando 31% das vezes; caso contrário, identifica corretamente 85% dos indivíduos saudáveis.

Vejamos o comportamento do modelo no arquivo de teste. 

n.Y = testing$status
n.X = as.matrix(testing[,c(2,5:17,20:24)])
n.pred0 = n.X%*%ajuste1$vec
n.PC1 = n.pred0[,1]
n.dados = data.frame(Y=n.Y,PC1=n.PC1)
novos.preditos = predict(ajuste8, newdata = n.dados, type = "response")
novos.Y = ifelse(family(ajuste8)$linkinv(novos.preditos)>0.8,1,0)
table(n.Y,novos.Y)
##    novos.Y
## n.Y  0
##   0 15
##   1 44

Observe que nosso modelo não erra na identificação das escolhidas como novas observações.


6. Referências