Construção de redes usando estatística clássica e Bayesiana-uma

Construção de redes usando estatísticaclássica e Bayesiana - uma comparação

Lina Dornelas Thomas

Dissertação apresentadaao

Instituto de Matemática e Estatísticada

Universidade de São Paulopara

obtenção do títulode

Mestre em Ciências

Programa: Estatística

Orientador: Prof. Dr. Anatoli Iambartsev

Durante o desenvolvimento deste trabalho o autor recebeu auxílionanceiro do CNPq

São Paulo, maio de 2012

Construção de redes usando estatísticaclássica e Bayesiana - uma comparação

Esta tese/dissertação contém as correções e alteraçõessugeridas pela Comissão Julgadora durante a defesarealizada por Lina Dornelas Thomas em 13/03/2012.

O original encontra-se disponível no Instituto deMatemática e Estatística da Universidade de São Paulo.

Comissão Julgadora:

• Prof. Dr. Anatoly Yambartsev (orientador) - IME-USP

• Profa. Dra. Florência Graciela Leonardi - IME-USP

• Prof. Dr. Iouri Soukhov - DPMMS/CMS-University ofCambridge

Agradecimentos

Gostaria de agradecer especialmente ao meu orientador, Anatoli Iambartsev,por todos os ensinamentos, paciência, dedicação e interesse pelo trabalho.

Agradeço também a minha mãe, Regina Andrade Dornelas, por seu in-centivo emocional e nanceiro e aos meus familiares que sempre me deramapoio incondicional.

Durante esses dois anos no instituto, contei com vários amigos para es-tudar para provas, fazer trabalhos e tirar dúvidas. Agradeço a cada umadessas pessoas. Em particular, gostaria de agradecer a Victor Fossaluza, umcolega e amigo que forneceu ideias e bastante suporte no estudo de Esta-tística Bayesiana, e a Renata Stella Khouri, pelas guras do Capítulo 7. Epor último, a Manuel González, que sempre demonstrou curiosidade pelapesquisa e, dessa forma, ajudou a manter minha motivação!

i

ii

Resumo

Construção de redes usando estatística clássica e Bayesiana -uma comparação

Nesta pesquisa, estudamos e comparamos duas maneiras de se construirredes. O principal objetivo do nosso estudo é encontrar uma forma efetivade se construir redes, especialmente quando temos menos observações do quevariáveis. A construção das redes é realizada através da estimação do coe-ciente de correlação parcial com base na estatística clássica (inverse method)e na Bayesiana (priori conjugada Normal - Wishart invertida). No presentetrabalho, para resolver o problema de se ter menos observações do que variá-veis, propomos uma nova metodologia, a qual chamamos correlação parciallocal, que consiste em selecionar, para cada par de variáveis, as demais variá-veis que apresentam maior coeciente de correlação com o par. Aplicamosessas metodologias em dados simulados e as comparamos traçando curvasROC. O resultado mais atrativo foi que, mesmo com custo computacionalalto, usar inferência Bayesiana é melhor quando temos menos observaçõesdo que variáveis. Em outros casos, ambas abordagens apresentam resultadossatisfatórios.

Palavras-chave: redes, correlação parcial, método inverso, estatística Baye-siana.

iii

iv

Abstract

Building complex networks through classical and Bayesianstatistics - a comparison

This research is about studying and comparing two dierent ways of buil-ding complex networks. The main goal of our study is to nd an eective wayto build networks, particularly when we have fewer observations than vari-ables. We construct networks estimating the partial correlation coecienton Classic Statistics (Inverse Method) and on Bayesian Statistics (Normal- Invese Wishart conjugate prior). In this current work, in order to solvethe problem of having less observations than variables, we propose a newmethodology called local partial correlation, which consists of selecting, foreach pair of variables, the other variables most correlated to the pair. Weapplied these methods on simulated data and compared them through ROCcurves. The most atractive result is that, even though it has high computati-onal costs, to use Bayesian inference is better when we have less observationsthan variables. In other cases, both approaches present satisfactory results.

Keywords: complex networks, partial correlation, inverse method, Bayesianstatistics.

v

vi

Sumário

1 Introdução 1

1.1 Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Organização do texto . . . . . . . . . . . . . . . . . . . . . . . 2

2 Conceitos Básicos 3

2.1 Coeciente de Correlaçao . . . . . . . . . . . . . . . . . . . . 32.2 Independência Condicional . . . . . . . . . . . . . . . . . . . . 42.3 Redes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.4 Dados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.5 Teorema de Bayes . . . . . . . . . . . . . . . . . . . . . . . . 72.6 Matriz Positiva Denida . . . . . . . . . . . . . . . . . . . . . 7

3 Distribuição Normal Multivariada e Wishart invertida 9

3.1 Multivariada Normal . . . . . . . . . . . . . . . . . . . . . . . 93.2 Distribuição Wishart . . . . . . . . . . . . . . . . . . . . . . . 113.3 Distribuição Wishart-invertida . . . . . . . . . . . . . . . . . 12

4 Independência Condicional 15

4.1 Correlação Parcial . . . . . . . . . . . . . . . . . . . . . . . . 154.1.1 Regressão Linear . . . . . . . . . . . . . . . . . . . . . 164.1.2 Método Inverso . . . . . . . . . . . . . . . . . . . . . . 174.1.3 Covariância Parcial . . . . . . . . . . . . . . . . . . . . 17

4.2 Correlação Condicional . . . . . . . . . . . . . . . . . . . . . . 18

5 Estatística Bayesiana 19

5.1 Distribuição a Priori . . . . . . . . . . . . . . . . . . . . . . . 195.2 Distribuição a Posteriori . . . . . . . . . . . . . . . . . . . . . 195.3 Distribuições a Priori Conjugadas . . . . . . . . . . . . . . . . 20

5.3.1 Multivariada Normal- Wishart Invertida . . . . . . . . 205.4 Full Bayesian Statistics Test . . . . . . . . . . . . . . . . . . . 21

vii

viii SUMÁRIO

6 Medidas de Diagnóstico 23

6.1 Teste de Hipóteses . . . . . . . . . . . . . . . . . . . . . . . . 236.2 Sensibilidade X Especicidade . . . . . . . . . . . . . . . . . . 236.3 Curva ROC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

7 Simulações 27

7.1 Algoritmos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 297.1.1 Árvores . . . . . . . . . . . . . . . . . . . . . . . . . . 297.1.2 Redes Mistas . . . . . . . . . . . . . . . . . . . . . . . 29

7.2 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 327.2.1 Estatística Clássica . . . . . . . . . . . . . . . . . . . . 327.2.2 Estatística Bayesiana . . . . . . . . . . . . . . . . . . . 327.2.3 Correlação Parcial Local . . . . . . . . . . . . . . . . . 34

7.3 Estudo das curvas ROC . . . . . . . . . . . . . . . . . . . . . 347.3.1 Comparação das Curvas ROC da Estatística clássica

com a Bayesiana . . . . . . . . . . . . . . . . . . . . . 37

A Demonstrações 39

A.1 Demonstração da fórmula (4.1) no caso unidimensional . . . . 39A.2 Equivalência do método inverso com a denição de correlação

parcial no caso unidimencional . . . . . . . . . . . . . . . . . 40A.3 Equivalência do método da covariância parcial com a denição

de correlação parcial no caso unidimencional . . . . . . . . . . 41

B Grácos 43

C Scripts em R 47

C.1 Gerador de árvores . . . . . . . . . . . . . . . . . . . . . . . . 47C.2 Gerador de redes mistas . . . . . . . . . . . . . . . . . . . . . 49C.3 Gerador de redes mistas piloto . . . . . . . . . . . . . . . . . 52C.4 Correlação parcial - método inverso . . . . . . . . . . . . . . . 53C.5 Correlação parcial - Bayesiana . . . . . . . . . . . . . . . . . . 55C.6 Correlação parcial local - método inverso . . . . . . . . . . . . 59C.7 Correlação parcial local- Bayesiana . . . . . . . . . . . . . . . 62C.8 Média das curvas ROC . . . . . . . . . . . . . . . . . . . . . . 70

Lista de Figuras

2.1 Representação de um grafo. Os cículos desempenham o papelde vértices, enquanto as linhas retas entre os vértices repre-sentam as ligações entre os mesmos. À esquerda apresentamosum grafo não direcionado, enquanto à direita temos um grafodirecionado e com peso. . . . . . . . . . . . . . . . . . . . . . 5

6.1 Exemplo de curva ROC . . . . . . . . . . . . . . . . . . . . . 25

7.1 Exemplo de árvores . . . . . . . . . . . . . . . . . . . . . . . . 287.2 Exemplo de redes mistas . . . . . . . . . . . . . . . . . . . . . 287.3 Gráco do número máximo de variáveis no cálculo da corre-

lação parcial local pela área das curvas ROC de 5 árvores (àesquerda) e 5 redes (à direita). . . . . . . . . . . . . . . . . . 34

7.4 Comparação da média (à esquerda) e mediana(à direita) dascurvas ROC das árvores com 50, 250, 500 e 1000 observaçõesna visão clássica. . . . . . . . . . . . . . . . . . . . . . . . . . 35

7.5 Comparação da média (à esquerda) e mediana(à direita) dascurvas ROC das árvores com 50, 250, 500 e 1000 observaçõesna visão Bayesiana. . . . . . . . . . . . . . . . . . . . . . . . . 36

7.6 Comparação da média (à esquerda) e mediana(à direita) dascurvas ROC das redes mistas com 50, 250, 500 e 1000 obser-vações na visão clássica. . . . . . . . . . . . . . . . . . . . . . 36

7.7 Comparação da média (à esquerda) e mediana(à direita) dascurvas ROC das redes mistas com 50, 250 e 500 observaçõesna visão Bayesiana. . . . . . . . . . . . . . . . . . . . . . . . . 37

7.8 Comparação entre as curvas ROC provenientes da estatísticaclássica e bayesiana para árvores com 50 (à esq), 250, 500 e1000 obs (à dir) . . . . . . . . . . . . . . . . . . . . . . . . . . 38

7.9 Comparação entre as curvas ROC provenientes da estatísticaclássica e bayesiana para redes mistas com 50, 250 e 500 obs) 38

B.1 Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das árvores com 50, 250, 500 e 1000 observações na visãoclássica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

ix

x LISTA DE FIGURAS

B.2 Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das árvores com 50, 250, 500 e 1000 observações na visãoBayesiana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

B.3 Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das redes mistas com 50, 250, 500 e 1000 observações navisão clássica. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

B.4 Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das redes mistas com 50, 250 e 500 observações na visãoBayesiana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Capítulo 1

Introdução

1.1 Motivação

A construção de redes é um tema cada vez mais presente entre os pesquisa-dores. Elas podem representar inúmeras situações dentre os mais variadostemas. Como exemplo podemos citar a internet, redes sociais de parentescoou outra característica que conecta indivíduos, redes organizacionais ou derelações entre empresas, redes neurais, cadeia alimentar, logística, redes dereferências entre artigos.

A maior motivação deste trabalho foi o artigo [11], onde redes regulatóriasforam construídas. Os dados foram tirados de uma base de um estudo caso-controle sobre câncer colo-uterino. A m de construir uma rede de expressãogênica, os autores usaram o método inverso para calcular o coeciente decorrelação parcial.

Contudo, nos estudos de microarray, é comum termos milhares de genesversus dezenas de observações, como acontece em [11]. Isso gera problemade identicabilidade, já que não temos informação suciente para fazer infe-rência sobre os dados. Estatisticamente, quanto mais conhecimento temos arespeito dos dados, mais realista são os resultados do estudo. Para estudardados com mais variáveis do que observações, é necessário a elaboração denovas metodologias.

Embora nosso objetivo neste estudo não seja encontrar a melhor metodo-logia dentre as já existentes, vimos que há diferentes maneiras para enfrentaressa situação. Em [13] a correlação parcial foi estimada através de regressãolinear em conjuntos esparsos. Na literatura, encontramos outras referênciascom aplicações em dados oriundos de microarray. Em [8], o coeciente decorrelação foi calculado sobre clusters de genes. Já nos artigos [20], [21],[22] e [18], modelos grácos Gaussianos (Graphical Gaussian Models) foramusados para construção de redes regulatórias. Em [20] também foi usadaanálise de cluster, enquanto que em [18] foi usado o lasso ponderado nosmodelos grácos Gaussianos. Conjuntos de genes foram usados em [7] para

1

2 CAPÍTULO 1. INTRODUÇÃO

a obtençao de redes de expressão gênica.

1.2 Objetivos

Os objetivos desse trabalho, utilizando dados simulados, são:

1. Estudar e comparar a construção de redes usando o coeciente de cor-relação parcial através de duas metodologias. A primeira é por meioda estatística clássica. Para a estimação do coeciente de correlaçãoparcial, utilizaremos o método inverso, como em [11], e, para o testede hipóteses, o valor-p calculado pela transformada-z de Fisher. A se-gunda metodologia é embasada na escola Bayesiana. Estimaremos amatriz de covariância parcial e usaremos as entradas dessa matriz parao cálculo do coeciente de correlação parcial. Testaremos as hipótesescom o valor-e vindo do FBST [14].

2. Estudar e comparar a construção de redes quando temos mais variáveisdo que observações utilizando o que chamamos de correlação parciallocal. Dizemos local porque, para cada par de variáveis, escolhemosapenas um subconjunto de variáveis correlacionadas com o par paracalcular o coeciente de correlação parcial. Esse estudo também éfeito usando as duas metodologias descritas no item anterior.

1.3 Organização do texto

No Capítulo 2 recordaremos alguns conceitos básicos que serão utilizados nodecorrer do texto. As denições, densidades e propriedades das distribuiçõesempregadas nas estimações encontram-se no Capítulo 3. Já o Capítulo 4 falasobre a denição de coeciente de correlação parcial e outras duas possíveismaneiras de obtê-lo. Noções básicas de estatística Bayesiana são introduzidasno Capítulo 5. As metodologias são comparadas através da curva ROC, aqual foi especicada no Capítulo 6. Finalmente, no Capítulo 7, estudaremosas simulações e discutiremos seus resultados.

Capítulo 2

Conceitos Básicos

Neste capítulo vamos introduzir algumas noções básicas, as quais serão usa-das ao longo desse trabalho. Iniciaremos com os conceitos relacionados coma construção de grafos: coeciente de correlação de Pearson, independênciacondicional, grafos e dados. Depois introduziremos dois teoremas bastanteconhecidos (lei da probabilidade total e teorema de Bayes) que serão usa-dos no Capítulo 5, onde falaremos um pouco sobre Estatística Bayesiana.Também apresentaremos a denição e algumas propriedades de matrizes po-sitivas denidas, essenciais nas distribuições Wishart e Wishart invertida, asquais serão tratadas no Capítulo 3.

2.1 Coeciente de Correlaçao

A correlação, também chamada coeciente de correlação, entre duas variáveisaleatórias, X e Y , é uma das medidas de dependência mais usadas e bemconhecidas na literaratura cientíca. Ela indica a força e a direção da relaçãolinear entre essas variáveis. A fórmula de maior uso comum é o coecientede correlação de Pearson. Uma de suas vantagens é que esse coeciente éadimensional e pode ser usado para comparar variáveis, pois é padronizadopelos desvios padrões.

Denição 2.1.1 (Coeciente de Correlaçao de Pearson). Sejam X e Y duasvariáveis aleatórias unidimensionais. O coeciente de correlação de Pearsonρ(X,Y ) entre X e Y é denida por

ρ(X,Y ) =Cov(X,Y )√V ar(X)V ar(Y )

.

Pela desigualdade de Cauchy-Schwartz podemos ver que |ρ(X,Y )| ≤ 1.Note que se X e Y forem não correlacionadas, então ρ(X,Y ) = 0, e que|ρ(X,Y )| = 1 se e somente se Y = βX, β 6= 0. Outra observação que mostracomo a correlação mede a relação linear entre as variáveis é que, quando

3

4 CAPÍTULO 2. CONCEITOS BÁSICOS

ρ(X,Y ) ∈ (−1, 1), podemos escrever Y = βX + ε, em que V ar(ε) 6= 0 eCov(Y, ε) = 0.

Denição 2.1.2 (Matriz de Correlação). Sejam X1, · · · , Xp variáveis ale-atórias unidimensionais. A matriz de correlação é simplesmente a matrizcujos elementos são os coecientes de correlação entre cada par de variáveisaleatórias (Xi, Xj), i, j = 1 · · · p, ou seja, -

ρ(X1, · · · , Xp) =

ρ(X1, X1) · · · ρ(X1, Xp)...

......

ρ(Xp, X1) · · · ρ(Xp, Xp)

,em que ρ(Xi, Xi) = 1,∀i = 1, · · · , p

Para saber mais sobre aspectos históricos do coeciente de correlação esobre medidas de dependência consulte [17] e [19] respectivamente.

2.2 Independência Condicional

SejamX e Y variáveis aleatórias denidas no espaço de probabilidade (Ω,A,P),com função densidade de probabilidade (f.d.p.) fX(x) e fY (y) respectiva-mente.

Denição 2.2.1 (Independência). Dizemos que as variáveis aleatóriasX e Ysão independentes e denotamos por X ⊥⊥ Y se sua f.d.p. conjunta fX,Y (x, y)pode ser escrita como o produto das densidades marginais:

fX,Y (x, y) = fX(x)fY (y),

ou, analogamente, se a densidade de Y dado X = x não depende de x, ouseja,

fY |X(y | x) = fY (y).

Denição 2.2.2 (Independência Condicional). Agora considere duas variá-veis aleatórias X, Y e um vetor aleatório Z = (Z1, · · · , Zp). Seja fX|Y,Z(x |y, z), z = (z1, · · · , zp), a f.d.p.de X condicional a Y e a Z, e fX|Z(x | z) af.d.p.de X condicional somente a Z. Dizemos que X e Y são condicional-mente independentes dado Z e denotamos por X ⊥⊥ Y | Z se, para cada valorde z, a seguinte equação é válida:

fX|Y,Z(x | y, z) = fX|Z(x | z).

Outra forma de denir a independência condicional, encontrada em [3],é se podemos fatorar a densidade conjunta fX,Y,Z(x, y, z) como um produtode duas funções, uma que não envolve x e uma não envolve y, isto é,

fX,Y,Z(x, y, z) = h(x, z)k(y, z).

2.3. REDES 5

Figura 2.1: Representação de um grafo. Os cículos desempenham o papelde vértices, enquanto as linhas retas entre os vértices representam as liga-ções entre os mesmos. À esquerda apresentamos um grafo não direcionado,enquanto à direita temos um grafo direcionado e com peso.

2.3 Redes

Vamos introduzir algumas denições relacionadas a redes, as quais utilizare-mos mais adiante.

Denição 2.3.1 (Grafo). Uma rede (ou grafo) é um conjunto nito V devértices, ou nós, que conectam-se entre si através de ligações ou elos. Oconjunto de elos em uma rede é denotado por E. Denotamos uma rede,mais conhecida como grafo no meio matemático, por G = (V,E). UsaremosX acompanhado de um índice para denotar os vértices, por exemplo, Xi,i = 1, · · · , p, p ∈ N, em que p é o número total de vértices em um grafo. Asligações entre Xi e Xj serão escritas das seguintes formas: [XiXj ] ou [XjXi].A Figura 2.1 mostra um diagrama exemplicando alguns tipos de redes.

Os elos podem ter pesos, por exemplo, a correlação entre dois genes, etambém estar apontados para uma única direção, isto é, direcionados. Grafoscontendo elos direcionados são chamados grafos direcionados, ou digrafos(digraphs). Um diagrama de um grafo direcionado e com peso pode serencontrado na Figura 2.1.

Denição 2.3.2 (Adjacência). Dizemos que Xi, Xj ∈ V são adjacentes,Xi ∼ Xj , se existe uma ligação entre eles, isto é, [XiXj ] ∈ E. Representamostodos os pares adjacentes através da matriz de adjacência Ap×p, em que seuselementos aij recebem 1 quando existe a ligação entre Xi e Xj e 0, caso


contrário. No caso da rede direcionada, aij = 1 indica que a ligação sai deXi e vai em direção a Xj , ou seja, a seta aponta para Xj . Se a direção forapenas em um sentido, de Xi para Xj , então aji = 0.

Denição 2.3.3 (Subgrafo). Qualquer subconjunto de vértices u ⊆ V geraum subgrafo de G. Denotaremos Gu = (u,H), sendo H as ligações em E asquais ambos os vértices de cada uma das ligações pertencem a u.

2.4 Dados

Os conjuntos de dados estudados neste trabalho são da seguinte forma:

Xn×p =

x11 x12 · · · x1p

x21 x22 · · · x2p...

......

...xn1 xn2 · · · xnp

,que são amostras de um vetor multidimensional e cada entrada da matrixxij , i = 1 · · ·n, j = 1 · · · p é a i-ésima observação da variável Xj . Essasobservações podem ser altura ou peso de um indivíduo, expressão gênica,disparos neuronais (spikes), etc.

Assumimos que as variáveis X1, X2, · · ·Xp seguem uma distribuição comfunção densidade de probabilidade conjunta

fθ(x1, x2, · · · , xp),

em que θ simboliza parâmetros desconhecidos. Os dados são obtidos de umexperimento no qual n observações são amostradas independentemente se-guindo fθ.

Nosso objetivo é trabalhar com modelos através de redes. Como descritoem ??, construimos um grafo G = (V,E), em que V é o conjunto de variáveisno modelo e E são os elos entre os pares de variáveis. Mas como estabeleceras ligações? Para cada par de variáveis Xi e Xj , se Xi ⊥⊥ Xj | X−ij , em queX−ij aloca todas as demais variáveis do modelo, então não há ligação entreelas. Entretanto, se as variáveis não forem condicionalmente independentes,dado X−ij , então existe um elo entre elas. Essa construção de grafo é co-nhecida como propriedade de Markov por par (pairwise Markov property).

Se a condição u ⊥⊥ v | w é válida para quaisquer conjuntos u, v e w ⊆ V,então w separa u e v no grafo de independência. Portanto, todas as inde-pendências condicionais e marginais podem ser diretamente obtidas do grafo(perfeição de Markov).

Outras propriedades de grafos podem ser encontradas em [5] e [12].

2.5. TEOREMA DE BAYES 7

2.5 Teorema de Bayes

Teorema 2.5.1 (Lei da Probabilidade Total). Suponha que os eventos A1,

A2, · · · , Ap em (Ω,F, P ) formam uma partição de Ω e P (Ai) > 0, i =1, 2, · · · , n. Então, para qualquer evento B

P (B) =n∑i=1

P (B | Ai)P (Ai)

Teorema 2.5.2 (Teorema de Bayes). Suponha que os eventos A1, A2, · · · , Anestão em (Ω,F, P ), formam uma partição de Ω e P (Ai) > 0, i = 1, 2, · · · , n.Seja A um evento qualquer com P (B) > 0. Então, ∀i = 1, 2, · · · , n, temos

P (Ai | B) =P (B | Ai)P (Ai)∑nj=1 P (B | Aj)P (Aj)

. (2.1)

2.6 Matriz Positiva Denida

Denição 2.6.1 (Matriz Positiva Denida). Uma matriz A é dita positivadenida se uTAu > 0 para todos os vetores u não nulos com entradas reais(u ∈ Rn), em que uT denota o vetor transposto de u.

Além da denição, existem várias condições que caracterizam quandouma matriz A é positiva denida (A > 0) ou semipositiva denida (A ≥ 0).Listaremos algumas delas. Para isso, considere λi os auto-valores de A e A∗

a matriz transposta conjugada de A.

1. A ≥ 0 se e somente se A for Hermitiana, ou seja, A = A∗, e λi ≥ 0,∀i.A > 0 se e somente se λi > 0,∀i

2. A ≥ 0 se e somente se A for Hermitiana e todos os seus menoresprincipais forem não negativos. A > 0 se e somente se todos os seusmenores principais forem positivos.

3. A ≥ 0 se e somente se A = B ∗ B para alguma matriz B. A > 0 se esomente se B for não singular.

4. A ≥ 0 se e somente se A = T ∗ T para alguma matriz triangularsuperior. Mais ainda, a diagonal de T pode ser formada por númerosnão negativos. Se A > 0, então T é única. Essa é a decomposição deCholesky de A. A > 0 se e somente se T for não singular.

5. A ≥ 0 se e somente se B2 para alguma matriz B. Essa matriz B éúnica. A > 0 se e somente se B > 0.

Um estudo mais detalhado sobre matrizes positivas denidas pode serencontrado em [2].


Capítulo 3

Distribuição Normal

Multivariada e Wishart

invertida

Neste capítulo falaremos um pouco sobre as distribuições consideradas nestetrabalho e suas propriedades. Assumiremos que os dados seguem uma distri-buição Multivariada Normal. A Wishart-invertida é utilizada nas distribui-ções a priori e posteriori, na Estatística Bayesiana (Capítulo 5). Citaremos adistribuição Wishart já que ela dá origem à distribuição Wishart-invertida.

3.1 Multivariada Normal

Nos estudos realizados em análise multivariada é comum que o pesquisadorassuma que os vetores aleatórios sigam uma distribuição normal multivari-dada. Abaixo listamos alguns motivos para tal procedimento.

• A distribuição Normal univariada pode ser facilmente generalizadaatravés da distribuição Normal multivariada, além de a análise destaser semelhante à do caso univariado.

• A distribuição multivariada normal é inteiramente denida pelo pri-meiro e segundo momento, o que implica uma menor quantidade deparâmetros quando comparada com outras distribuições, diminuindo,assim, a diculdade dos problemas de estimação.

• Correlação nula nesse tipo de distribuição implica independência. Ecorrelação nula dois a dois implica independência total. A prova dissopode ser encontrada em [10].

• Mesmo quando os dados originais não são normais, podemos apelarpara o Teorema do Limite Central, se temos um número grande dedados.

9

10CAPÍTULO 3. DISTRIBUIÇÃO NORMALMULTIVARIADA EWISHART INVERTIDA

Denição 3.1.1. Dizemos que X = (X1, · · · , Xp) segue uma distribuiçãoNormal p-variada se e somente se, ∀a = (a1, · · · , ap) ∈ Rp, a variável aTX éNormal univariada. A distribuição Normal p-variada é denotada da seguinteforma:

X ∼ Np(µ,Σ),

em que µp×1 é a média do vetor X e Σp×p é a matriz de covariância de X.

Densidade

Sua f.d.p., quando Σ é uma matriz positiva denida (Σ > 0), tem aseguinte forma:

f(X) = (2π)−p/2|Σ|−1/2e−12

(x−µ)T Σ−1(x−µ) (3.1)

Como já mencionamos anteriormente,

Esperança

E(X) = µ

Variância

V ar(X) = Σ.

Características e Propriedades

Teorema 3.1.1. Particionemos X em X = (X(1)T , X(2)T ) ∼ Np(µ,Σ), emque Σ é particionada da seguinte forma:

Σ =

[Σ11 Σ12

Σ21 Σ22

]A distribuição condicional de X(1)T dado X(2)T é

X(1)T | X(2)T ∼ Nr(µ2 + Σ21Σ−111 (X − µ1),Σ22.1). (3.2)

em que r é o comprimento do vetor X2 e

Σ22.1 = Σ22 − Σ21Σ−111 Σ12. (3.3)

A prova desse teorema pode ser encontrada em [10].

3.2. DISTRIBUIÇÃO WISHART 11

3.2 Distribuição Wishart

A distribuição Wishart é uma generalização matricial da distribuição qui-quadrado. É uma família de distribuições denida sobre matrizes aleatóriassimétricas e positivas denidas. Essas distribuições são de grande impor-tância na estimação de matrizes de covariância na estatística multivariada.Na inferência Bayesiana, assumindo que os dados seguem uma distribuiçãoNormal multivariada, ela é a priori conjugada da matriz de precisão (matrizde covariância inversa).

Denição 3.2.1. Suponha que X é uma matriz k×p, com cada linha vindoindependentemente de uma distribuição Normal p-variada com média 0 ematriz de covariância S:

Xi = (Xi1, Xi2, · · ·Xip) ∼ Np(0, S), i = 1, · · · , k

A distribuição Wishart é a distribuição de probabilidade de uma matrixaleatória M = XTX, também conhecida como matriz de dispersão(scattermatrix). Notacionalmente:

M ∼Wp(S, k)

em que S é a matriz de escala, S > 0 e k representa os graus de liberdade,k > 0.

Densidade

Seja M uma matriz simétrica p× p, aleatória e positiva denida. Seja Suma matriz constante também positiva denida de tamanho p× p. Quandok ≥ p (vide Observação 3.2.1), M tem uma distribuição Wishart com matrizde escala S e k graus de liberdade se tiver uma f.d.p. dada por

f(M) =|M |(k−p−1)/2 exp[−1

2 tr(S−1M)]

2kp/2|S|k/2Γp(k2

) (3.4)

em que |M | e |S| são os determinantes de M e S respectivamente e Γp(·) éa função gamma multivariada denida como

Γp

(k

2

)= πp(p−1)/a

p∏j=1

Γ

[n

2+

(1− j)2

](3.5)

Observação 3.2.1. Lembre que o posto (rank) de uma matriz Xk×p é omínimo entre k e p, ou seja, r(X) ≤ min(k, p). Lembre também quer(XTX) = min(r(XT ), r(X)) ≤ min(k, p). Se k < p, então r(M) < p.Logo, a matriz M não terá posto completo e, portanto, será singular, isto é,será não invertível.


Esperança

E(M) = E(XTX) =k∑i=1

E(XiXTi ) = kS.

Variância

V ar(M)ij = k(s2ij + siisjj),

com i, j = 1, · · · , p e sij representando os elementos ij da matriz S.Características e Propriedades

Segundo [10], quando p = 1, a distribuição W1(s2, k) é dada por XTX,em que os elementos de Xk×1 são variáveis iid N(0, s2), isto é, a distribuiçãoWishart W1(s2, k) é a mesma que a distribuição s2χ2

k.

A matriz de escala S tem a mesma função na distribuição Wishart ques2 tem na distribuição s2χ2

k.

3.3 Distribuição Wishart-invertida

A distribuição Wishart-invertida (inverse-Wishart) é a generalização multi-variada da distribuição Gamma-invertida (veja [16] para mais informações).Se M é uma matriz que segue a distribuição Wishart, então sua inversa,M−1 segue a distribuição Wishart-invertida. Essa distribuição é a prioriconjugada natural para a matriz de covariância quando os dados seguemuma distribuição multivariada Normal.

Densidade

Se M ∼ Wp(S, k), em que S > 0 e k ≥ p (vide Observação 3.2.1), entãoU = M−1 segue a distribuição Wishart-invertida IWp(Ψ, k), sendo Ψ = S−1.Sua f.d.p. é dada por:

g(U) =|U |−(k+p+1)/2 exp[−1

2 tr(ΨU−1)]

2kp/2|Ψ|−k/2Γp(k2

) (3.6)

Esperança

E[U ] =Ψ

k − p− 1(3.7)

Variância

V ar(σii) =2ψ2

ii

(k − p− 1)2(k − p− 3), k − p > 3, (3.8)

3.3. DISTRIBUIÇÃO WISHART-INVERTIDA 13

V ar(σij) =ψiiψjj + k−p+1

k−p−1ψ2ij

(k − p)(k − p− 1)(k − p− 3)(3.9)

em que ψij é o elemento ij da matriz Ψ.Características e Propriedades

Teorema 3.3.1. Seja U ∼ IWp(Σ, k),Σ > 0, p < k entao

U22.1 ∼ IWq(Σ22.1, k − q), (3.10)

em que U22.1 = U22 − U21U−111 U12.

Demonstração.

U ∼ IWp(Σ, k)⇒ U ∼Wp(Σ−1, k)

Agora, considere U particionada da seguinte forma:

U =

(U11 U12

U21 U22

),

em que U11 é uma matriz p× p, U12, p× q, U12 : q × p, U22 : q × q.

Como podemos ver em [15] U−1 pode ser escrita da seguinte forma:

M = U−1 =

(U−1

11.2 −U−111 U12U22.1−1

−U−122 U21U11.2−1 U−1

22.1

), (3.11)

em que U−111.2 = (U11 − U12U22−1U21)−1

Por teorema visto em [10], as matrizes diagonais de U−1 também seguemuma distribuição Wishart, logo

(U−1)22 ∼W ((Σ−1)22, k)

ou seja,

U−122.1 ∼W (Σ−1

22.1, k)⇒ U22.1 ∼ IW ((Σ22.1, k)


Capítulo 4

Independência Condicional

Como já vimos no Capítulo 2, independência condicional é fundamental namodelagem gráca, pois denimos um elo entre dois vértices se e somente seas variáveis do modelo não forem condicionalmente independentes. Entre-tanto, apesar de ser baseada na teoria das probabilidades, seu cálculo podeser bastante complexo. Por isso, na prática, é razoável substituir indepen-dência condicional por correlação parcial igual a 0 ou correlação condicionaligual a 0. Optamos por trabalhar com a correlação parcial nula, já que estaé mais fácil de ser calculada do que correlação condicional, a qual dependeda forma da distribuição.

Em [1] prova-se que correlação parcial é idêntica a correlação condicio-nal se a correlação condicional é independente da condição e também se aesperança condicional é linear.

Neste capítulo deniremos correlação parcial e apresentaremos duas ou-tras formas mais simples de calculá-la. Provaremos que essas formas sãoequivalentes à denição original. Introduziremos também a noção de corre-lação condicional.

4.1 Correlação Parcial

Estamos interessados na independência condicional de Xi, n×1, e Xj , n×1,dado Y = X−ij , i, j = 1, · · · p, i 6= j, especialmente através da correlação par-cial. Podemos calcular a correlação parcial por três caminhos que acabamgerando os mesmos resultados. O primeiro caminho é através da deniçãooriginal usando regrassão linear. O primeiro caminho é o jeito mais tra-balhoso, o segundo é pelo Método Inverso e o terceiro utiliza a covariânciaparcial.

15

16 CAPÍTULO 4. INDEPENDÊNCIA CONDICIONAL

4.1.1 Regressão Linear

Denição 4.1.1 (Correlação Parcial). Seja Y um vetor linha aleatório decomprimento p − 2 correlacionado com Xi e Xj . Se calcularmos o coe-ciente de correlação entre Xi e Xj , provavelmente concluiremos que existeuma certa correlação entre elas. Contudo, não sabemos se essa correlaçãovem da correlação dessas duas variáveis com Y, ou se elas são realmentecorrelacionadas independente de Y. Se mantivermos Y constante, será queuma alteração de Xi inuenciará Xj? A ideia principal dessa abordagem éjustamente retirar o efeito de Y, examinando as correlações entre os resíduosdas projeções de Xi e Xj no espaço linear gerado por Y, ou seja, subtraindoaquela parte da relação linear devido a Y. Matematicamente temos:

Xi = αi + Yβi + εi

Xj = αj + Yβj + εj ,

em que βi = (βi1, · · · , βi(p−2))T .

Através da regressão pelo método dos mínimos quadrados, obtemos Xi

e Xj , que são os estimadores de Xi e Xj respectivamente, e seus resíduosdados por:

Resi = Xi − Xi

Resj = Xj − Xj .

Note que, Resi eResj são ortogonais aY e portanto a correlação ρ(Resi, Resj)é a correlação entre os componentes de Xi e Xj que não apresentam depen-dência linear com Y. Portanto, a correlação parcial será dada por

ρij.Y = ρ(Resi, Resj) = ρ(Xi − Xi, Xj − Xj)

Se ρij.Y ≈ 0, então Xi e Xj são não correlacionados se não levarmos Yem consideração.

Lemma 4.1.1. Mantendo as notações da Denição 4.1.1,

ρ(Resi, Resj) =ρXiXj − ρXiYρXjY√1− ρ2

XiY

√1− ρ2

XjY

,

ou seja, o cáculo da correlação parcial pode vir direto da fórmula:

ρij.Y =ρ(Xi, Xj)− ρ(Xi,Y)ρ(Xj ,Y)√1− ρ(Xi,Y)2

√1− ρ(Xj ,Y)2

, (4.1)

4.1. CORRELAÇÃO PARCIAL 17

A prova do Lema 4.1.1, considerando Y uma variável unidimensional,pode ser encontrada no Apêndice A.1. [19] prova para o caso multidimensi-onal. Um exemplo desse caso é quando

Y = (X1, · · · , Xi−1, Xi+1, · · · , Xj−1, Xj+1, · · · , Xp),

como foi denotado na Seção 2.4.

Observação 4.1.1. Note que, em geral, a correlação parcial não é o mesmoque correlação condicional, embora na literatura exista um certo grau deconfusão sobre estes dois conceitos. Isto se deve à existência de diferen-tes denições de correlação condicional e também de casos para o qual acorrelação parcial e a condicional são equivalentes. Notadamente, no casogaussiano, em que a correlação parcial nula indica independência condicio-nal, obtém-se os mesmos valores para as duas correlações. Uma discussãobastante interessante sobre as diferenças e condições de equivalência entreos dois conceitos pode ser encontrada em [1]. Para nalizar esta pequenaconsideração, em geral, a palavra condicional é reservada para quantidadesem que elas mesmas são variáveis aleatórias, o que não é o caso da correlaçãoparcial que é sempre um valor não aleatório. [19]

4.1.2 Método Inverso

Existe outra abordagem que nos permite calcular todas as correlações parci-ais entre qualquer par de variávies aleatórias Xi e Xj , dado Y = Y−ij . Sejaa matriz de correlação R, em que rij = ρ(Xi, Xj) são os elementos da matrizR. Se R é invertível, então denimos P = R−1. A correlação parcial é dadapor:

ρij.Y = − pij√piipjj

(4.2)

Note que, se o número de observações (linhas) n da matriz X for menorque o número de variáveis aleatórias (colunas) p, então R será singular. Paraobtermos R−1 é necessário ao menos um método para fazer a matriz R naosingular, positiva denida.

A demonstração da equivalência do método inverso com a denição ori-ginal de correlação parcial pode ser encontrado no Apêndice A.2

4.1.3 Covariância Parcial

A matrix de covariância parcial de Xij = (Xi, Xj) é denotada por

Σij.Y =

[σii.Y σij.Yσji.Y σjj.Y

],

a qual pode ser calculada através da fórmula

18 CAPÍTULO 4. INDEPENDÊNCIA CONDICIONAL

Σij.Y = ΣXij − ΣXijYΣ−1Y ΣYXij , (4.3)

particionando a matriz de covariância de (Xij ,Y) em

Cov

([Xij

Y

])=

[ΣXij ΣXijY

ΣYXij ΣY

],

em que ΣXij é 2 × 2, ΣXijY é 2 × (p − 2), ΣYXij é (p − 2) × 2 e ΣY é(p− 2)× (p− 2).

A correlação parcial é dada por

ρij.Y =σij.Y√

σii.Yσjj.Y. (4.4)

Observe que a fórmula (4.3) é igual à fórmula da variância da distribuiçãocondicional da normal multivariada (3.2).

Através do Apêndice A.3 podemos ver que essa última abordagem é equi-valente à abordagem original, quando as variáveis Xi, i = 1, · · · , p têm médianula.

4.2 Correlação Condicional

Podemos calcular a correlação condicional de forma análoga à correlação par-cial utilizando a covariância condicional. A matriz de covariância condicionalde Xi e Xj dado Y é denida por

Cov(Xi, Xj | Y) = E((Xi − E(Xi | Y)

)(Xj − E(Xj | Y)

)| Y)

= σij|Y.

A covariância condicional será denotada da seguinte forma:

ΣXX|Y =

[σii|Y σij|Yσji|Y σjj|Y

]Calculamos a correlação condicional como segue.

ρij|Y =σij|Y

√σii|Yσjj|Y

Essa é uma fórmula comumente considerada para averiguar a indepen-dência condicional, entretanto, em [1], prova-se que a correlação parcial iguala 0 ou a correlação condicional igual a 0 não implica em independência con-dicional, exceto no caso da distribuição normal.

Capítulo 5

Estatística Bayesiana

Seja X = (X1, X2, · · · , Xp) um vetor aleatório denido em (Ω,A,Pθ). Soba perspectiva baysiana, a incerteza sobre o parâmetro desconhecido θ é des-crito em forma de probabilidade. Deste modo, atribui-se uma distribuiçãode probabilidade para θ, usualmente chamada de distribuição a priori. De-pois de observarmos os dados, obtemos uma distribuição de θ condicionalà informação trazida pela amostra, também conhecida como distribuição aposteriori.

5.1 Distribuição a Priori

Considere os problemas em que as observações vêm de uma distribuiçãof(x | θ), sendo θ os parâmetros com valores desconhecidos. Dessa forma,consideraremos θ um vetor aleatório denido em (Θ,B, π). Muitas vezes,antes que qualquer observação de f(x | θ) esteja disponível, o experimen-tador ou o estatístico conseguirá resumir seu conhecimento e informaçãoanteriores sobre θ ∈ Θ por meio de uma distribuição de probabilidade paraθ. Em outras palavras, antes que qualquer dado seja coletado ou observado,a experiência e conhecimento do experimentador no passado vai levá-lo aacreditar que é mais provável que θ se encontre em certas regiões de Θ doque em outras.

5.2 Distribuição a Posteriori

A distribuição de θ, condicional aos dados observados é chamada distribuiçãoa posteriori de θ. O teorema de Bayes (2.5) nos revela como calcular a f.d.pa posteriori de θ após observarmos os dados. Vamos rever o teorema deBayes usando a notação especíca das distribuições a priori e a posterioridos parâmetros.

Suponha queX1, · · · ,Xn sejam amostras aleatórias independentes de umvetor aleatório X = (X1, · · · , Xp), que, condicionadas ao conhecimento de

19

20 CAPÍTULO 5. ESTATÍSTICA BAYESIANA

θ, são i.i.d. A f.d.p. conjunta será

f(x | θ) = f(x1, · · · ,xn | θ) = fp(x1 | θ) · · · fp(xn | θ) (5.1)

.Seja Θ o espaço paramétrico contendo θ. Agora suponha que θ siga a

distribuição com f.d.p. ξ(θ), a qual podemos chamar de distribuição a prioride θ. Assim, usando a Lei da Probabilidade Total (2.5.1), podemos calculara f.d.p. conjunta de (X1, · · · ,Xn) da seguinte forma:

g(x) =

∫Θf(x | θ)ξ(θ)dθ (5.2)

A distribuição a posteriori de θ é simplesmente a distribuição de θ apósos valores de X1, · · · ,Xn serem observados, ou seja, a f.d.p. condicional deθ dado X1 = x1, · · · ,Xn = xn, denotada por ξ(θ | x). Utilizando o Teoremade Bayes (2.5.2), podemos escrever ξ(θ | x) como segue:

ξ(θ | x) =f(x | θ)ξ(θ)

g(x)para θ ∈ Θ (5.3)

Como o denominador da equação (5.3) não é uma função de θ, quandoconhecermos os valores de x1, · · · ,xn, g(x) será uma constante. Logo, po-demos reescrever (5.3) através da seguinte relação:

ξ(θ | x) ∝ f(x | θ)ξ(θ). (5.4)

Observe que f(x | θ), quando em função apenas de θ, é chamada defunção de verossimilhança e denotada por Lx(θ).

5.3 Distribuições a Priori Conjugadas

Quando a distribuição a posteriori pertencer à mesma família de distribuiçõesque a distribuição a priori, então chamados essa família de família conjugada.Normalmente, escolher uma priori que pertença a uma família conjugada(priori conjugada), simplica os cálculos da distribuição a posteriori.

5.3.1 Multivariada Normal- Wishart Invertida

Neste estudo, trabalharemos com a distribuição multivariada normal commédia, µ, e matriz de covariância, Σ, desconhecidas.

Suponhamos, a priori, que Σ siga uma distribuição Wishart Invertidacom parâmetros k0 e Ψ0 e também que µ | Σ siga uma distribuição normalcom média λ0, variância Σ

v0, em que v0 ∈ R é um parâmetro de escala, Ψ0 é

uma matriz positiva denida p× p e λ0 é um vetor com tamanho p. Usandonotação:

5.4. FULL BAYESIAN STATISTICS TEST 21

Σ ∼ IW (k0,Ψ0)

µ | Σ ∼ Normal(λ0,1

v0Σ)

Fazendo os cálculos, é possível ver que a posteriori teremos

Σ ∼ IW (k1,Ψ1)

µ | Σ ∼ Normal(λ0,1

v0Σ),

em que

λ1 =nx+ v0λ0

n+ v0,

v1 = n+ v0,

k1 = n+ k0 e

Ψ1 = Ψ0 +

n∑i=1

(xi − x)(xi − x)T +nv0

n+ v0(x− λ0)(x− λ0)T ,

sendo x a média amostral.

Os leitores interessados podem obter mais observações em [4].

5.4 Full Bayesian Statistics Test

Para denir esse teste, vamos focar apenas no espaço de probabilidade aposteriori (Θ,B, ξx), em que Θ é o espaço paramétrico e, nesse caso, o espaçoamostral. B é a σ-álgebra gerada por Θ e ξx é a medida de probabilidadecondicional em B dado o ponto amostral observado x, ou seja, ξx = ξ(θ | x).

Agora denimos ξ(θ∗ | x) como o máximo da densidade a posteriori sobrea hipótese nula, Θ0, em que

θ∗ ∈ argmaxθ∈Θ0ξ(θ | x).

Agora, considere o subconjunto de Θ em que ξ(θ | x) > ξ(θ∗ | x),

T ∗ = θ ∈ Θ | ξ(θ | x) > ξ(θ∗ | x)

A credibilidade de T ∗ é sua probabilidade a posteriori,

K∗ =

∫T ∗ξ(θ | x)dθ,

em que f∗(θ | x) = f(θ | x) se f(θ | x) > f(θ∗ | x) e 0 caso contrário.A medida de evidência, ou e-valor, da hipótese nula é o complemento da

probabilidade do conjunto T ∗, ou seja,

22 CAPÍTULO 5. ESTATÍSTICA BAYESIANA

Ev(H) = 1−K∗.

Note que, quando K∗ for grande, o espaço paramétrico nulo encontra-seem uma reagião de baixa probabilidade e a evidência é alta, ou seja, contraa hipótese nula. Por outro lado, quando K∗ for pequena, então o espaçoparamétrico nulo está em uma reagião de alta probabilidade e, portanto,evidência baixa, isto é, a favor da hipótese nula.

Capítulo 6

Medidas de Diagnóstico

Agora que já conhecemos todos os métodos, como escolher o melhor? Comodenir o modelo mais adequado aos dados? As Medidas de Diagnósticoservem justamente para essa nalidade, pois medem a ecácia do modelo.Nesse capítulo, apresentamos as medidas usadas neste trabalho.

6.1 Teste de Hipóteses

Na modelagem gráca segundo [5], como já mencionamos nos Capítulos 2 e4, um elo entre duas variáveis só existe se essas variáveis não são condici-onalmente independentes. Devido a maior facilidade dos cálculos, optamospor usar a correlação parcial para substituir a independência condicional.Para isso, zemos o seguinte teste de hipóteses para cada par de variáveis(Xi, Xj):

H0 : ρij.Y = 0

H1 : ρij.Y 6= 0

Se rejeitarmos H0, então traçamos um elo entre as variáveis.

6.2 Sensibilidade X Especicidade

Primeiramente vamos introduzir alguns conceitos para o nosso estudo.

• Verdadeiros Positivos (VP): Pares corretamente classicados como cor-relacionados.

• Falsos Positivos (FP): Pares incorretamente classicados como corre-lacionados.

• Verdadeiros Negativos (VN): Pares corretamente classicados comonão correlacionados.

23

24 CAPÍTULO 6. MEDIDAS DE DIAGNÓSTICO

• Falsos Negativos (FN): Pares incorretamente classicados como nãocorrelacionados.

Calculamos a sensibilidade da seguinte forma:

Sensitividade =#verdadeiros positivos

#verdadeiros positivos + #falsos negativos

A sensitividade serve para relacionar a habilidade do teste de identicarresultados positivos. Modelos com sensibilidade alta apresentam erro tipo 2pequeno.

Já a especicidade tem a seguinte fórmula:

Especificidade =#verdadeiros negativos

#verdadeiros negativos + #falsos positivos

A especicidade serve para relacionar a habilidade do teste de identicarresultados negativos. Modelos com especicidade alta apresentam erro tipo1 pequeno.

Podemos também calcular a acurácia do modelo:

Acuracia =V P + V N

V P + V N + FP + FN

6.3 Curva ROC

A Curva ROC é o gráco da sensibilidade × (1 − especicidade). Tambémé conhecido como curva Receiving Operating Characteristics (ROC) porqueé uma comparação de duas características operacionais (Sensibilidade e Es-pecicidade) a medida que o critério aumenta. Como critério, usaremos ovalor-p da estatística clássica ou o valor-e da estatística Bayesiana. Veja umexemplo na Figura 6.1.

6.3. CURVA ROC 25

Figura 6.1: Exemplo de curva ROC

26 CAPÍTULO 6. MEDIDAS DE DIAGNÓSTICO

Capítulo 7

Simulações

Nessa parte do trabalho, a m de conrmar a ecácia da teoria descrita an-teriormente e de constatarmos o melhor método, simulamos 2 tipos de redes,as quais chamamos árvores e redes mistas. Para cada tipo de rede, zemos500 simulações com 100 variáveis, calculamos a correlação parcial, o p-valor eo e-valor. No caso de termos mais variáveis do que observações, computamosa correlação parcial local. Para averiguar a qualidade do modelo, traçamosa curva ROC nas visões clássica e Bayesiana para cada árvore, calculamosa curva ROC média e depois comparamos a área dessas curvas médias. Ainterpretação é que quanto maior a área, melhor o modelo. Isso foi feito paranúmeros de observações diferentes: 50, 250, 500 e 1000.

Árvores são redes em que cada nó possui apenas um pai"e zero ou maislhos,exceto o primeiro nó, o qual não possui pai. Um exemplo de umaárvore pode ser encontrado na Figura 7.1. A distribuição Lognormal(1,1)foi escolhida para simular a quantidade de lhos (distribuição de saída).Já as redes mistas que apresentam distribuição de entrada (distribuição donúmero de pais) e distribuição de saída. Como exemplo, temos a Figura7.2. No nosso caso vamos trabalhar com redes mistas com distribuiçaõ desaída igual a das árvores, ou seja, Lognormal(1,1) e distribuição de entradaGeométrica(0.6).

A escolha dessas distribuições de entrada e saída foram baseadas nosresultados do trabalho [11], que concluem que as redes genéticas possuemdistribuição de entrada próxima a uma Geométrica(0.6) e distribuição desaída com cauda pesada apresentando média e variância próximas à distri-buição Lognormal(1,1).

Mas como denimos esse parentesco de pais e lhos? Simples. Os lhossão uma combinação linear dos pais com um acréscimo de um ruído branco,ou seja, uma perturbação seguindo uma distribuição normal com média 0 evariância bem pequena. No nosso caso, utilizamos variância igual a 0, 01.

27

28 CAPÍTULO 7. SIMULAÇÕES

Figura 7.1: Exemplo de árvores

Figura 7.2: Exemplo de redes mistas

7.1. ALGORITMOS 29

7.1 Algoritmos

7.1.1 Árvores

Como já descrevemos anteriormente, o número de lhos de cada nó segueuma distribuição lognomal(1,1).

A construção das árvores foi feita da seguinte maneira: Considerandoque estamos trabalhando com p vértices,

Algorithm 7.1.1 (Árvores). 1. A árvore começa com um vértice na pri-

meira geração.

2. Geramos o númerod n1 de lhos de X1 com seguindo uma distribuição

Lognormal(1, 1), com a condiçao que o primeiro vértice tenha ao menos

um descendente.

3. Para cada vértice da segunda geração, (X2, · · · , X1+n1), geramos o nú-

mero de lhos n2i, i = 1, · · · , n1, seguindo uma distribuição Lognor-

mal(1,1), com a condição que a segunda geração tenha ao menos um

descendente e 1 + n1 + n2 ≤ p, n2 = n21 + · · ·+ n2n1.

4. Para cada geração seguinte, repetimos o item anterior até obtermos pvariáveis.

Agora que já temos a estrutura de parentesco, precisamos simular os

dados.

1. Geramos a primeira variável X1 seguindo uma distribuição Normal(0,1).

2. A relação entre pais e lhos vem de uma combinação linear. Gera-

mos os p− 1 coecientes de p− 1 combinações lineares seguindo uma

distribuição Uniform(0, 7; 2). Os sinais positivo e negativo dos coeci-

entes seguem uma distribuição Bernoulli(23), em que sucesso é o sinal

positivo.

7.1.2 Redes Mistas

Como já descrevemos anteriormente, o número de lhos de cada nó segue umadistribuição Lognormal(1,1) e o número de pais segue uma Geométrica(0.6).

A construção das árvores foi feita da seguinte maneira: Considerandoque estamos trabalhando com p vértices,

Algorithm 7.1.2 (Redes Mistas). 1. Geramos o número de variáveis na

primeira geração, n1, com uma distribuição Lognormal(1,1)


2. Geramos o número de lhos dessas variáveis seguindo a distribuição

Lognormal(1,1), com a restrição que a geração tenha no mínimo 1

lho.

3. Para o primeiro vértice da segunda geração, geramos o número de pais

seguindo uma Geométrica(0.6), com a restrição que esse vértice te-

nha pelo menos um pai e no máximo a quantidade total de vértices da

geração anterior que tenha ao menos um lho.

4. Para especicar quem são os pais desse primeiro vértice da segunda

geração, usamos uma distribuição em que as probabilidades são os pesos

do número de lhos de cada vértice da geração anterior sobre o número

total de lhos dessa mesma geração.

7.1. ALGORITMOS 31

5. Como cada lho da geração anterior conecta com um vértice da geração

em questão, vericamos se há ainda lhos desconectados. Se houver,

criamos mais um vértice nessa geração e geramos o número de pais

como no item anterior.

6. Para especicar quem são os pais desse novo vértice, usamos uma dis-

tribuição em que as probabilidades são os pesos do número de lhos

desconectados de cada vértice da geração anterior sobre o número total

de lhos desconectados dessa mesma geração.

7. Seguimos assim sucessivamente até conectarmos todos os lhos a um

vértice da segunda geração.

Agora que já temos a estrutura de parentesco, precisamos simular os

dados.

1. Geramos as variáveis da primeira geração com uma distribuição p1-

variada Normal(0, Ip1).

2. A relação entre pais e lhos vem de uma combinação linear. Para

cada ligação entre pai e lho, geramos um coeciente seguindo uma


distribuição uniforme entre 0, 7 e 2. Os sinais positivo e negativo dos

coecientes seguem uma distribuição Bernoulli(23), em que sucesso é o

sinal positivo.

3. Para cada vértice a partir da segunda geração, fazemos as combinações

lineares do vértice com seus pais e adicionamos um ruído branco.

7.2 Metodologia

7.2.1 Estatística Clássica

O primeiro passo foi simular árvores e redes mistas como descrito nos al-goritmos 7.1.1 e 7.1.1. Com os dados simulados, estimaremos a correlaçãoparcial através do método inverso, equação (4.2).

Como estamos procurando correlações parciais não nulas ou nulas paradenirmos se há ou não ligação entre cada par de vértices, precisamos testara hipótese nula da correlação parcial de Xi e Xj dado Y ser igual a 0, ouseja, H0 : ρij.Y = 0 contra a alternativa HA : ρij.Y 6= 0. Para isso, podemosusar a transformada-z de Fisher:

z(ρij.Y ) =1

2ln(

1 + ρij.Y1− ρij.Y

).

Rejeitamos H0 com nível de signicância α se:√n− |Y | − 3 · |z(ρij.Y )| > Φ−1(1− α/2), (7.1)

em que Φ(·) é a função de distribuição acumulada da Normal(0, σ2) e n é onúmero de observações.

Para mais informações veja [9] e [6].

7.2.2 Estatística Bayesiana

Escolha da priori

Assim como já citamos no Capítulo 5, assumimos que os dados seguem umadistribuição Normal com média µ e matriz de covariância Σ. Devido à maiorfacilidade dos cálculos, optamos por trabalhar com a priori conjugada, ouseja, Σ ∼ IW (k0,Ψ0) e µ | Σ ∼ (λ0,

1v0

Σ).

Mas, se não temos nenhuma informação sobre os parâmetros, como de-vemos escolher os parâmetros da distribuição a priori? Nessas situações, umpesquisador pode trabalhar com uma amostra piloto que vem da mesma po-pulação que os dados do estudo e, a partir dela, fazer suas suposições sobre os

7.2. METODOLOGIA 33

dados para a escolha da priori. Nesse trabalho, optamos por trabalhar comamostras piloto com o mesmo número de observações que os dados de estudo.

Com essa nalidade, simulamos duas árvores com a mesma matriz deadjacência e mesma dependência linear. Uma funciona como amostra pilotoe a outra é vista como os dados do estudo. Fizemos o mesmo para redesmistas. Denotaremos a matriz piloto por Z e seus elementos, por zij

Deveríamos escolher uma priori com média próxima ao que esperamos,mas com variância grande. Lembre que, no caso da Wishart invertida,

E(Σ) =Ψ0

k0 − p− 1, k0 − p > 1.

Se repararmos as fórmulas das variâncias (3.8) e (3.9), percebemos quequanto menor for k0, maior será a variância. Então, como, para E(Σ) existir,k0 > p+ 1, escolhemos

k0 = p+ 3.

Baseado na esperança, poderíamos escolher Ψ0 = (k0 − p − 1) ∗ Q, tal queQ seja uma matriz próxima do que esperamos da matriz Σ. Por essa razão,optamos por trabalhar com a matriz de covariância estimada de Z, ou seja,

Ψ0 = (k0 − p− 1) ∗ ΣZ .

Usamos o estimador de máxima verossimilhança da covariância, ou seja,

σij =1

n

n∑i=1

(zij − zj)(zik − zk),

em que σij são as entradas da matriz ΣZ .Utilizaremos a fórmula (4.4), que vem da covariância parcial (4.3) para

o cálculo da correlação parcial.Para testar a hipótese nula de correlção parcial igual a 0, usaremos o

teste FBST que foi explicado na seção 5.4. Note que precisamos da distri-buição da correlação parcial. Porém, o mais perto que conseguimos chegarfoi encontrar a distribuição da covariância parcial, que também segue umaWishart invertida. Veja o Teorema 3.3.1.

Como não temos a distribuição da correlação parcial, simulamos a co-variância parcial 1000 vezes e estimamos a densidade da correlação parcialatravés da função density do programa estatístico R, possibilitando o cálculodo e-valor. Para calcularmos a área da densidade que nos fornecerá o e-valor,utilizamos a aproximação pela área dos retângulos.


7.2.3 Correlação Parcial Local

Já se sabe que, em modelos multivariados, quando o número de observaçõesé menor que a quantidade de variáveis no modelo, as estatísticas envolvendotodas as variáveis são viesadas, e também não encontramos na literaturacomo se calcula o valor-p para essas situações.

Como alternativa, utilizamos o que chamamos de correlação parcial lo-cal. Nesse cálculo, para cada par de vértices correlacionados, começamosconsiderando apenas os vértices correlacionados com ao menos um dos ele-mentos do par, os quais denominamos vizinhança do par. Se essa vizinhançacontiver mais vértices do que observações, então selecionamos os n

10 vérticescom as maiores correlações no caso de árvores e n

2 vértices no caso de redesmistas. Esses valores foram selecionados de acordo com o número de vérticesque geravam, para cada tipo de rede, a maior área abaixo da curva ROC, Se-ção 6.3. Devido ao alto custo computacional da abordagem Bayesiana, esseestudo foi feito apenas na abordagem clássica. Por conseguinte, utilizamoso mesmo número de vértices para o cáculo da correlação parcial local emambas as abordagens.

Figura 7.3: Gráco do número máximo de variáveis no cálculo da correlaçãoparcial local pela área das curvas ROC de 5 árvores (à esquerda) e 5 redes(à direita).

7.3 Estudo das curvas ROC

A m de fazer uma análise de diagnóstico, zemos um estudo das curvasROC. Simulamos 500 redes de cada tipo e calculamos a sensibilidade e espe-cicidade de cada rede com o nível de signicância α variando de 0 a 1. Em

7.3. ESTUDO DAS CURVAS ROC 35

seguida, para cada valor de α, computamos as médias, medianas, primeiroquartil e terceiro quartil da sensibilidade e da especicidade e calculamosuma curva ROC com cada uma dessas medições. No caso Bayesiano, comoo custo computacional é alto, optamos por trabalhar com apenas 50 redespara calcular a média das curvas ROC, com exceção da correlação parciallocal das redes mistas, em que foram usadas apenas 20 redes. Nesse últimocaso, também devido ao custo computacional alto, apenas zemos o estudopara 50, 250 e 500 obs.

Se repararmos as Figuras 7.4, 7.5, 7.6 e 7.7, podemos ver que as áreasdas curvas aumentam à medida que o número de observações aumenta, comojá era esperado. O que chama mais atenção é que a área da curva referenteà 50 observações está consideravelmente próxima à da curva de 250 obser-vações, especialmente no gráco de árvores na visão Bayesiana, Figura 7.5.Vale ressaltar que estamos trabalhando com redes compostas de 100 variá-veis, ou seja, as curvas das redes com menos observações que parâmetrosse encontram relativamente próximas às curvas com mais observações queparâmetros. A m de conrmar que as curvas ROC apresentam pequenadispersão, o leitor encontrará grácos com a média, mediana, primeiro eterceiro quartis no Apêndice B.

Figura 7.4: Comparação da média (à esquerda) e mediana(à direita) das cur-vas ROC das árvores com 50, 250, 500 e 1000 observações na visão clássica.


Figura 7.5: Comparação da média (à esquerda) e mediana(à direita) dascurvas ROC das árvores com 50, 250, 500 e 1000 observações na visão Baye-siana.

Figura 7.6: Comparação da média (à esquerda) e mediana(à direita) dascurvas ROC das redes mistas com 50, 250, 500 e 1000 observações na visãoclássica.

7.3. ESTUDO DAS CURVAS ROC 37

Figura 7.7: Comparação da média (à esquerda) e mediana(à direita) das cur-vas ROC das redes mistas com 50, 250 e 500 observações na visão Bayesiana.

7.3.1 Comparação das Curvas ROC da Estatística clássicacom a Bayesiana

Como podemos ver na Figura 7.8 à direita, quando temos mais observações(250, 500, 1000) do que variáveis (100), as médias das curvas ROC na abor-dagem clássica e na Bayesiana cam sobrepostas, o que indica que ambasapresentam resultados tão bons quanto. Em contrapartida, a curva ROCna visão Bayesiana para árvores com apenas 50 observações apresenta umaárea bem maior que a curva ROC clássica. Já no caso das redes mistaspara 50, 250 e 500 observações, as médias das curvas ROC na abordagemclássica e na Bayesiana cam sobrepostas. Concluimos que, em geral, é me-lhor usarmos a estatística clássica, pois esta apresenta resultados tão bonsquanto da estatística Bayesiana, mas possui um custo computacional maisbaixo. Entretanto, no caso de árvores com problema de identicabilidade,apesar do custo computacional ser maior na estatítica Bayesiana, o seu usoé aconselhável.


Figura 7.8: Comparação entre as curvas ROC provenientes da estatísticaclássica e bayesiana para árvores com 50 (à esq), 250, 500 e 1000 obs (à dir)

Figura 7.9: Comparação entre as curvas ROC provenientes da estatísticaclássica e bayesiana para redes mistas com 50, 250 e 500 obs)

Apêndice A

Demonstrações

A.1 Demonstração da fórmula (4.1) no caso unidi-mensional

Demonstração. Sejam Xi, Xj , Y variáveis aleatórias e εi, εj os erros dasregressões de Xi em Y e de Xj em Y respectivamente. Por denição,

ρij.Y = ρ(εi, εj)

=Cov(εi, εj)√V ar(εi)V ar(εj)

(A.1)

Observe que

Cov(εi, εj) = Cov(εi, Xj − Y βj − αj)= Cov(εi, Xj)

= Cov(Xi − Y βi − αi, Xj)

= Cov(Xi, Xj)− βiCov(Xj , Y )) (A.2)

Por outro lado,

Cov(Xi, Y ) = Cov(αi + βiY + εi, Y )

= βiV ar(Y )⇒⇒ βi = Cov(Xi, Y )V ar(Y )−1 (A.3)

Logo,

Cov(εi, εj) = Cov(Xi, Xj)−Cov(Xi, Y )Cov(Xj , Y )

V ar(Y )) (A.4)

39

40 APÊNDICE A. DEMONSTRAÇÕES

Obtemos a variância de εi da seguinte maneira:

V ar(εi) = Cov(εi, εi)

= Cov(εi, Xi − βiY − αi)= Cov(εi, Xi)

= E(εiXi)− E(εi)E(Xi) (A.5)

Como, neste trabalho, Xi tem média 0,

V ar(εi) = E((Xi − βiY − αi)Xi)

= E(X2i )− βiE(Y Xi)

= V ar(Xi)−Cov(Xi, Y )2

V ar(Y )(A.6)

A variância de εj segue de forma análoga. Substituindo (A.6) e (A.4) em(A.1) temos

ρij.Y =Cov(Xi, Xj)− Cov(Xi,Y )Cov(Xj ,Y )

V ar(Y )√V ar(Xi)− Cov(Xi,Y )2

V ar(Y )

√V ar(Xj)− Cov(Xj ,Y )2

V ar(Y )

=Cov(Xi, Xj)− Cov(Xi,Y )Cov(Xj ,Y )√

V ar(Y )√V ar(Y )√

V ar(Xi)(

1− Cov(Xi,Y )2

V ar(Xi)V ar(Y )

)√V ar(Xj)

(1− Cov(Xj ,Y )2

V ar(Xj)V ar(Y )

)

=

Cov(Xi,Xj)√V ar(Xi)V ar(Xj)

− Cov(Xi,Y )Cov(Xj ,Y )√V ar(Xi)V ar(Y )

√V ar(Xj)V ar(Y )√(

1− Cov(Xi,Y )2

V ar(Xi)V ar(Y )

)√(1− Cov(Xj ,Y )2

V ar(Xj)V ar(Y )

)=

ρ(Xi, Xj)− ρ(Xi, Y )ρ(Xj , Y )√(1− ρ(Xi, Y )2) (1− ρ(Xj , Y )2)

A.2 Equivalência do método inverso com a deni-ção de correlação parcial no caso unidimencio-nal

Demonstração. Sejam Xi, Xj e Y variáveis aleatórias e R a matriz de cor-relação de Xi, Xj e Y , ou seja,

A.3. EQUIVALÊNCIA DOMÉTODODA COVARIÂNCIA PARCIAL COMADEFINIÇÃODE CORRELAÇÃO PARCIAL NO CASO UNIDIMENCIONAL41

Y =

1 ρ(Xi, Xj) ρ(Xi, Y )ρ(Xj , Xi) 1 ρ(Xj , Y )ρ(Y,Xi) ρ(Y,Xj) 1

=

1 ρ(Xi, Xj) ρ(Xi, Y )ρ(Xi, Xj) 1 ρ(Xj , Y )ρ(Xi, Y ) ρ(Xj , Y ) 1

Uma das formas de obtermos a matriz inversa é através da matriz adjuntaadj(R), que é a transposta da matriz dos cofatores. Lembre que o cofatorRij é o determinante da matriz resultante de retirar de R a linha i e a colunaj (isso é, o determinante menor) multiplicado por (−1)i+j .

P = R−1 = |R|−1adj(R) (A.7)

Pela fórmula (4.2), vimos que ρij.Y = − pij√piipjj

, em que, de acordo com

(A.7)

pii = 1− ρ(Xj , Y )2

pjj = 1− ρ(Xi, Y )2

pij = −(ρ(Xi, Xj)− ρ(Xi, Y )ρ(Xj , Y )

)

Portanto, substituindo, temos que

ρij.Y = − ρ(Xi, Xj)− ρ(Xi, Y )ρ(Xj , Y )√(1− ρ(Xi, Y )2)(1− ρ(Xj , Y )2)

A.3 Equivalência do método da covariância parcialcom a denição de correlação parcial no casounidimencional

Iniciaremos demonstrando a fómula da covariância parcial.

Demonstração. Sejam Xi, i = 1, · · · , p variáveis aleatórias unidimensionais eY um vetor linha aleatório (p-2)-dimensional.

Primeiro vamos provar que a variância parcial de Xi.Y é

V ar(Xi.Y ) = V ar(Xi)− Cov(Xi, Y )V ar(Y )−1Cov(Y,Xi)

Seja εi o erro da regressão de Xi em Y . Assim,

42 APÊNDICE A. DEMONSTRAÇÕES

Xi = αi + Y βi + εi ⇒ εi = Xi − βiY − αi (A.8)

. Logo,Lembre de (A.6) que

V ar(Xi.Y ) = V ar(εi)

= E(εiXi)− E(εi)E(Xi)

= E(X2i )− βiE(Y Xi) (A.9)

Como vimos em (A.3):

βi = Cov(Xi, Y )V ar(Y )−1 (A.10)

Substituindo (A.10) em (A.9) temos que

V ar(Xi.Y ) = E(X2i )− Cov(Xi, Y )V ar(Y )−1E(Y Xi) (A.11)

Portanto,

V ar(Xi · Y ) = V ar(Xi)− Cov(Xi, Y )V ar(Y )−1Cov(Y,Xi) (A.12)

O cálculo da covariância parcial é análogo ao da variância parcial:

Cov(Xi, Xj .Y ) = V ar(Xi : Xj)−Cov(Xi : Xj , Y )V ar(Y )−1Cov(Y,Xi : Xj)

Segundo (A.12),

V ar(Xi : Xj · Y ) = V ar(Xi : Xj)− Cov(Xi : Xj , Y )V ar(Y )−1Cov(Y,Xi : Xj)

=

[σ2i σij

σji σ2j

]−[σiY (σ2

Y )−1σY i σiY (σ2Y )−1σY j

σjY (σ2Y )−1σY i σjY (σ2

Y )−1σY j

]A continuação da prova da equivalência do método da covariância parcial,

Seção 4.1.3, com a denição, Seção 4.1.1, para Y unidimensional é análogaao Apêndice A.1

Apêndice B

Grácos

Figura B.1: Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das árvores com 50, 250, 500 e 1000 observações na visão clássica.

43

44 APÊNDICE B. GRÁFICOS

Figura B.2: Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das árvores com 50, 250, 500 e 1000 observações na visão Bayesiana.

45

Figura B.3: Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das redes mistas com 50, 250, 500 e 1000 observações na visão clássica.

46 APÊNDICE B. GRÁFICOS

Figura B.4: Quartil 1 (0.25), média, mediana e quartil 3 (0.75) das curvasROC das redes mistas com 50, 250 e 500 observações na visão Bayesiana.

Apêndice C

Scripts em R

C.1 Gerador de árvores

# Di s t r i b u i ç ão dos f i l h o sd i s t=function (n)

x=f loor (rlnorm (n , 1 , 1 ) )return ( x )

# p é o número t o t a l de v é r t i c e sfor (p in c (100 ) )for ( i 1 in c ( 1 ) )

# Número de observaçõesn=5000

# Calculando os c o e f i c i e n t e s

coe f 1=c (round( runif (p , 0 . 7 , 2 ) , 2 ) )s i n a l=sample (c (1 ,−1) ,p , replace=TRUE, c (2/3 ,1/3) )coe f 1=coe f 1∗ s i n a l

# Matriz de ad jacênc iaA=matrix (0 , ncol=p ,nrow=p)

# Matriz com os dadosX=matrix ( , ncol=p ,nrow=n)X[ ,1 ]=rnorm(n)

# Lis ta com número de f i l h o s de cada ramo por geraçãoger=l i s t ( )

47

48 APÊNDICE C. SCRIPTS EM R

# Vetor com número de f i l h o s de cada geraçãof i l h o s=vector ( )

# Gerando número de f i l h o s na 1a geração ex i g indo que ao# menos um f i l h o e re spe i t ando o número de v a r i á v e i s p

repeat

ger [ [ 1 ] ] = d i s t (1 )i f ( ger [ [ 1 ] ] !=0 && ger [ [ 1 ] ] < p)break

f i l h o s [1 ]=sum( ger [ [ 1 ] ] )

# Calculando as observações da 1a geração e as entradas da# matr iz de ad jacênc ia

i=1cont=iwhile ( i<=ger [ [ 1 ] ] )

i=i+1X[ , i ]= coe f 1 [ i ] ∗X[ ,1 ]+rnorm(n , 0 , 0 . 0 1 )A[ 1 , i ]=1A[ i ,1 ]=1

k=2cont2=1while ( i<p)

# Gerando , para cada geração k , o número de f i l h o s# de cada ramo da geração anter ior , r e spe i t ando o# número de v a r i á v e i s p e ex i g indo ao menos um f i l h o# em cada geração k

repeat

ger [ [ k ] ]= d i s t (sum( ger [ [ k−1 ] ] ) )f i l h o s [ [ k ] ]=sum( ger [ [ k ] ] )i f (sum( f i l h o s [ 1 : k])>=p)

soma=sum( f i l h o s [ 1 : k−1])l=0while ( soma<p)

l=l+1

C.2. GERADOR DE REDES MISTAS 49

soma=soma+ger [ [ k ] ] [ l ]ger [ [ k ] ] [ l ]= ger [ [ k ] ] [ l ]−(soma−(p−1))ger [ [ k ] ] [ l +1: length ( ger [ [ k ] ] ) ]= 0f i l h o s [ [ k ] ]=sum( ger [ [ k ] ] )break

elsei f (sum( ger [ [ k ] ] ) !=0) break

# Calculando , para cada geração k , as observações e# as entradas da matr iz ad jacênc ia o r i g i n a l

for ( j in 1 : length ( ger [ [ k ] ] ) )i f ( ger [ [ k ] ] [ j ] !=0)

cont=iwhile ( i<cont+ger [ [ k ] ] [ j ] )

i=i+1X[ , i ]= coe f 1 [ i ] ∗X[ , cont2+j ]+

rnorm(n , 0 , 0 . 0 1 )A[ cont2+j , i ]=1A[ i , cont2+j ]=1

cont2=cont2+sum( ger [ [ k−1 ] ] )k=k+1

write . table (X, f i l e=paste ("arvore" , i1 , "_" ,p , ".txt" ,

sep="" ) )write . table (A, f i l e=paste ("original" , i1 , "_" ,p , ".txt" ,

sep="" ) )write . table ( coe f1 , f i l e=paste ("coeficientes" , i1 , "_" ,p , ".txt" ,

sep="" ) )

C.2 Gerador de redes mistas

require (mnormt)

# Di s t r i b u i ç ão do número de f i l h o sd i s t . out=function (n)


x=f loor (rlnorm (n , 1 , 1 ) )return ( x )

# Di s t r i b u i ç ão do número de pa i sd i s t . in=function (n)

x=rgeom(n , 0 . 6 )return ( x )

# Número de observaçõesn=10000# Número de v a r i á v e i sp=100

for ( i 3 in 1 : 500 )

X=matrix (ncol=p ,nrow=n)A=matrix (0 , ncol=p ,nrow=p)

# Número de f i l h o snf=l i s t ( )

# Número de pa i snp=l i s t ( )pa i . todos=l i s t ( )

# pai srepeat

n1=d i s t . out (1 )i f ( n1 !=0) break

pai . todos [ 1 : n1]=0

pai=vector ( length=n1 )

X1=rmnorm(n , mean = rep (0 , n1 ) , diag ( n1 ) )X[ , 1 : n1 ] = X1

# Usando como d i s t r i b u i ç ã o de cauda pesada uma# log normal

# Número de f i l h o sn2=0k=0

C.2. GERADOR DE REDES MISTAS 51

while ( n1<p)

k=k+1repeatnf [ [ k ] ]= d i s t . out ( length ( pa i ) )i f (sum( nf [ [ k ] ] ) !=0) break

temp=nf [ [ k ] ]

np [ [ k ] ]=vector ( )pa i=l i s t ( )

i=0while (sum( temp) !=0)

i=i+1nmax=length ( temp [ temp !=0 ] )

np [ [ k ] ] [ i ]= min((1+ d i s t . in ( 1 ) ) , nmax)prob=vector ( )for ( j in 1 : length ( temp ) )prob [ j ]=temp [ j ] /sum( temp)

pai [ [ i ] ]= sample (c ( 1 : length ( temp ) ) ,np [ [ k ] ] [ i ] , replace=FALSE,prob )

temp [ pai [ [ i ] ] ]= temp [ pai [ [ i ] ] ] −1

i f ( i+n1<=p)

A[ i+n1 , n2+pai [ [ i ] ] ]=1A[ n2+pai [ [ i ] ] , i+n1]=1pai . todos [ [ i+n1 ] ]=n2+pai [ [ i ] ]

n2=n1n1=n1+length ( pa i )

# Determinado os c o e f i c i e n t e scoef=l i s t ( )for ( i in 1 : p )

coef [ [ i ] ]= c (round( runif ( length (


pai . todos [ [ i ] ] ) , 0 . 7 , 2 ) , 2 ) )s i n a l=sample (c (1 ,−1) , length ( pa i . todos [ [ i ] ] ) ,

replace=TRUE, c (2/3 ,1/3) )coef [ [ i ] ]= coef [ [ i ] ] ∗ s i n a l

# Combinações Linearesfor ( i in 1 : p )

i f (sum( pa i . todos [ [ i ] ] )==0) X[ , i ]=X[ , i ] elsei f ( length ( pa i . todos [ [ i ] ] )==1)

X[ , i ]= X[ , pa i . todos [ [ i ] ] ] ∗coef [ [ i ] ]+rnorm(n , 0 , 0 . 0 1 ) else

X[ , i ]=X[ , pa i . todos [ [ i ] ] ]%∗%as .matrix (coef [ [ i ] ] )+rnorm(n , 0 , 0 . 0 1 )

max. np=max( sapply (np , function ( x ) max( x ) ) )

Coef=matrix (0 , p ,max. np )for ( j in 1 : p ) Coef [ j , 1 : length ( pa i . todos [ [ j ] ] ) ] =

coef [ [ j ] ]

write . table (X, f i l e=paste ("arvore_inout_" , i3 , "_" ,p ,".txt" , sep="" ) )

write . table (A, f i l e=paste ("original_inout_" , i3 , "_" ,p ,".txt" , sep="" ) )

write . table ( Coef , f i l e=paste ("coeficientes_inout_" ,i3 , "_" ,p , ".txt" , sep="" ) )

C.3 Gerador de redes mistas piloto

require (mnormt)

#Número de observaçõesn=10000

for (p in 100)for ( i 3 in 2 : 500 )

A=read . table (paste ("original_inout_" , i3 , "_" ,p ,

".txt" , sep="" ) )c o e f i c i e n t e s=read . table (paste ("coeficientes_

inout_" , i3 , "_" ,p , ".txt" , sep="" ) )

pa i . todos=l i s t ( )coef=l i s t ( )

C.4. CORRELAÇÃO PARCIAL - MÉTODO INVERSO 53

for ( k in 1 : ncol (A) )

pai . todos [ [ k ] ] = which(A[ , k ] [ row(A)<k]==1)coef [ [ k ] ]= t ( as .matrix ( c o e f i c i e n t e s [ k ,

1 : length ( pa i . todos [ [ k ] ] ) ] ) )

n1=0for ( k in 1 : length ( pa i . todos ) )

i f (sum( ! i s .na( pa i . todos [ [ k ] ] ) )==1) breakn1=n1+1

X=matrix (nrow=n , ncol=p)

X1=rmnorm(n , mean = rep (0 , n1 ) , diag ( n1 ) )X[ , 1 : n1 ] = X1

for ( i in 1 : p )i f (sum( pa i . todos [ [ i ] ] )==0) X[ , i ]=X[ , i ] elsei f ( length ( pa i . todos [ [ i ] ] )==1) X[ , i ]=X[ , pa i . todos [ [

i ] ] ] ∗coef [ [ i ] ]+rnorm(n , 0 , 0 . 0 1 ) elseX[ , i ]=X[ , pa i . todos [ [ i ] ] ]%∗%as .matrix ( coef [ [

i ] ] )+rnorm(n , 0 , 0 . 0 1 )

write . table (X, paste ("arvore_piloto_inout_" , i3 , "_" ,p , "_.txt" , sep="" ) )

C.4 Correlação parcial - método inverso

require ( corpcor )require ( igraph )

# Função que c a l c u l a a cor r e l a cao e o pva lo rcor . prob <− function (X, d f r = nrow(X) − 2)

R <− cor (X)above <− row(R) < col (R)r2 <− R[ above ]^2Fstat <− r2 ∗ d f r / (1 − r2 )R[ above ] <− 1 − pf ( Fstat , 1 , d f r )R


# Função que c a l c u l a a matr iz com pva l o r e s da cor r e l a ção# pa r c i a l o b t i da a t ravé s do método inve r so

ginvmet2=function (X, alpha1 )

p=ncol (X)

# Matriz de cova r i a c i acovX=cov (X)

# Matriz das co r r e l a co e scorX=cor (X)

# P−va l o r é a par te t r i n g u l a r de cima da matr izpva lcor=cor . prob (X)

# Selec ionando os pares s i gn i f i c an t emen t e# cor re l a c i onadosA=matrix (99 , ncol=ncol (X) ,nrow=ncol (X) )for ( i in 1 : ( ncol (X)−1))for ( j in ( i +1) :(ncol (X) ) )

A[ i , j ]= i f e l s e ( pva l cor [ i , j ]<alpha1 , corX [ i , j ] ,0) A[ j , i ]=A[ i , j ]

diag (A)=0

# Calculando a cor re l a cao p a r c i a l a t r ave s do metodo# inver soco r r . pos<−make . p o s i t i v e . d e f i n i t e ( covX)co r r . inv<−cor2pcor ( co r r . pos )

# Calculando o pva lo r da cor re l a cao p a r c i a lco r r . inv . p<−array (99 ,dim=c (p , p ) )

for ( i in 1 : ( p−1))for ( j in ( i +1):p )

i f (A[ i , j ] !=0)co r r . inv . p [ i , j ]<−1−pnorm(abs ( sqrt (n−

(p−2)−3)∗log ((1+ cor r . inv [ i ,j ] ) /(1− co r r . inv [ i , j ] ) ) /2) )

l i s t a=l i s t (A, co r r . inv . p)return ( l i s t a )

p=100for (n in c (250 ,500 ,1000) )

C.5. CORRELAÇÃO PARCIAL - BAYESIANA 55

for ( i 1 in c ( 1 : 5 0 0 ) )

# Lendo os dadosX=as .matrix ( read . table ( f i l e=paste ("arvore_inout_" ,

i1 , "_" ,p , ".txt" , sep="" ) ) ) [ 1 : n , ]

colnames (X)=c ( seq ( 1 : p ) )

# Chamando a matr iz de pva l o r do método inve r so para# cor parc . Esta função devo l v e uma l i s t a . 1) matr iz# de ad jacênc ia com as co r r e l a çõ e s s i g n i f i c a n t e s . 2)# matr iz com pva l o r e s

A=ginvmet2 (X, 0 . 0 1 )

write . table (A [ [ 2 ] ] , f i l e=paste ("pvalor_arvore_inout_", i1 , "_" ,p , "_" ,n , ".txt" , sep="" ) )

C.5 Correlação parcial - Bayesiana

require (mnormt)require (MCMCpack)require ( corpcor )require ( igraph )require ( zoo )

# Função rwish para casos em que S nao é p o s i t i v a d e f i n i d a no# método cho l

r iw i sh2=function (v , S )

return ( solve ( rwish2 (v , solve (S ) ) ) )

rwish2=function (v , S )

i f ( ! i s .matrix (S ) )S <− matrix (S)

i f (nrow(S) != ncol (S ) ) stop ( message = "S not square in rwish().\n" )

i f ( v < nrow(S ) )

stop ( message = "v is less than the dimension of S in

rwish().\n" )


p <− nrow(S)CC <− chol (S , p ivot=TRUE)Z <− matrix (0 , p , p )diag (Z) <− sqrt ( rchisq (p , v : ( v − p + 1 ) ) )i f (p > 1)

pseq <− 1 : ( p − 1)Z [ rep (p ∗ pseq , pseq ) + unlist ( lapply ( pseq , seq ) ) ] =

rnorm(p ∗ (p − 1)/2)return ( crossprod (Z %∗% CC))

# Função f b s t

f b s t . area=function (X, c t e )

Y <− density ( subset (X, subset=(X<1)))

i f (min(Y$x ) >0| |max(Y$x)<0) return (0 )else

x=min(Y$x [Y$x>cte ] )i f ( x==In f ) x=max(Y$x [Y$x<cte ] )

xt <− Y$x [5]−Y$x [ 4 ]yt <− Y$y [Y$y<Y$y [Y$x==x ] ]

return (sum( xt∗yt ) )

# Função que c a l c u l a cor re l a ção e pva lo r

cor . prob <− function (X, d f r = nrow(X) − 2) R <− cor (X)above <− row(R) < col (R)r2 <− R[ above ]^2Fstat <− r2 ∗ d f r / (1 − r2 )R[ above ] <− 1 − pf ( Fstat , 1 , d f r )R

p i l o t o=1

for (p in c (100 ) )for (n in c (1000) )for ( i 1 in c ( 1 : 5 0 ) )

# Lendo os dados

C.5. CORRELAÇÃO PARCIAL - BAYESIANA 57

X=as .matrix ( read . table ( f i l e=paste ("arvore_inout_" , i1 ,"_" ,p , ".txt" , sep="" ) ) ) [ 1 : n , ]

X. p=as .matrix ( read . table ( f i l e=paste ("arvore_piloto_ inout_" , i1 , "_" ,p , "_.txt" , sep="" ) ) ) [ 1 : n , ]

colnames (X)=c ( seq ( 1 : p ) )

Xor i g ina l=XXor i g ina l . p=X. p

covX=cov (X. p)

corX=cor (X)pva lcor=cor . prob (X)

alpha=0.001

A=matrix ( , ncol=ncol (X) ,nrow=ncol (X) )for ( i in 1 : ( ncol (X)−1))for ( j in ( i +1) :(ncol (X) ) )

A[ i , j ]= i f e l s e ( pva l cor [ i , j ]<alpha , corX [ i , j ] , 0 )A[ j , i ]=A[ i , j ]

diag (A)=0

# Bayesiana

# Construindo a rede usando rwish2

# A p r i o r i

mediaX=colMeans (X. p , na .rm = FALSE, dims = 1)

lambda0=mediaXk0=p+3ps i 0=(k0−p−1)∗covXv0=1

# A po s t e r i o r i

k1=k0+nsoma=0for ( i in 1 : n )soma=soma + (X[ i , ]−mediaX)%∗%t (X[ i , ]−mediaX)ps i 1=ps i 0+soma+(n∗v0/ (n+v0 ) )∗ (mediaX−lambda0 )%∗%

t (mediaX−lambda0 )

covpar=matrix ( )


corpar=vector ( )

f b s t=matrix ( rep (0 , p^2) ,ncol=p ,nrow=p)

for ( i in 1 : ( p−1))for ( j in ( i +1) :(p ) )i f (A[ i , j ] !=0)

ps i 1 . i j=ps i 1 [ c ( i , j ) , c ( i , j )]− ps i 1 [ c ( i , j ) ,−c ( i , j ) ]%∗%solve ( p s i 1 [−c ( i , j ) ,−c ( i ,j ) ] )%∗%ps i 1 [−c ( i , j ) , c ( i , j ) ]

for ( k2 in 1 :1000)

i f ( i s . p o s i t i v e . d e f i n i t e ( solve (p s i 1 . i j ) , method="chol")==

TRUE)

covar=r iw i sh ( k1 , p s i 1 . i j ) else

invcovar=rwish2 ( k1 , p s i 1 . i j )

while ( i s . p o s i t i v e . d e f i n i t e (invcovar)==FALSE)invcovar=rwish2 ( k1 ,

p s i 1 . i j )

covar=solve ( invcovar )

covpar=covar

corpar [ k2 ]=( covpar [ 1 , 2 ] ) /sqrt (covpar [ 1 , 1 ] ∗covpar [ 2 , 2 ] )

f b s t [ i , j ]= f b s t . area ( corpar , 0 )f b s t [ j , i ]= f b s t [ i , j ]

write . table ( fbs t , f i l e=paste ("fbst_global_arvore_ piloto_inout_" , i1 , "_" ,p , "_" ,n , ".txt" , sep="" ) )

C.6. CORRELAÇÃO PARCIAL LOCAL - MÉTODO INVERSO 59

C.6 Correlação parcial local - método inverso

require (MCMCpack)require ( corpcor )require ( igraph )


#ve to r com opcoes de número de observações

########lendo os dados

for (p in c (100 ) )for (n in 50)for ( i 3 in 1 : 50 )

i f (p>=n)

#X=as . matrix ( read . t a b l e ( f i l e=pas t e (" arvore_" , i3 ,"_",p , ".txt" , sep="" ) ) ) [ 1 : n , ]

X=as .matrix ( read . table ( f i l e=paste ("arvore_inout_" ,i3 , "_" ,p , ".txt" , sep="" ) ) ) [ 1 : n , ]

p=ncol (X)colnames (X)=c ( seq ( 1 : ncol (X) ) )

Xor i g ina l=X

covX=cov (X)corX=cor (X)

pva lcor=cor . prob (X)alpha=0.01

A=matrix ( , ncol=p ,nrow=p)colnames (A)=c ( 1 : p )for ( i in 1 : ( p−1))for ( j in ( i +1):p )A[ i , j ]= i f e l s e ( pva l cor [ i , j ]<alpha , corX [ i , j ] , 0 )


A[ j , i ]=A[ i , j ]diag (A)=0

g1=graph . adjacency (A,mode=("undirected" ) ,weighted=TRUE)

V( g1 )$name <− c ( seq ( 1 : p ) )

# Acessando os e l o s

edges=get . a d j l i s t ( g1 )

Ap=matrix (99 , ncol=ncol (X) ,nrow=ncol (X) )Af=matrix ( rep (0 , p^2) ,ncol=ncol (X) ,nrow=ncol (X) )

for ( k in 1 : ncol (X) )

v e r t i c e s=edges [ [ k ] ]

i f (sum( ! i s .na( v e r t i c e s ) ) !=0)for ( j in 1 : length ( v e r t i c e s ) )i f ( v e r t i c e s [ j ]+1>k && A[ k , v e r t i c e s [ j ]+

1 ] !=0)

i 1=as . character ( k )j 1=as . character ( v e r t i c e s [ j ]+

1)v i z l i s t=graph . neighborhood (g1 , 1 , nodes=c ( i1 , j 1 ) ,mode=c ("all" ) )

v i z=sort ( as . integer (unique (c (V( v i z l i s t [ [ 1 ] ] ) $name ,V(v i z l i s t [ [ 2 ] ] ) $name ) ) ,d e c r ea s ing=TRUE) )

#Diminuindo a quant idade de v i z inhos , caso tenhamos# mais v i z i n ho s do que observações

i f ( length ( v i z)>=n)

pos i cao . i=vector ( )pos i cao . j=vector ( )

# Pegando as n maiores co r r e l a co e s de cada v é r t i c e do# par

cor . maiores . i=sort (A[viz , c ( k ) ] [−k ] ,d e c r ea s ing=TRUE) [ 1 : n ]cor . maiores . j=sort (A[

C.6. CORRELAÇÃO PARCIAL LOCAL - MÉTODO INVERSO 61

viz , c ( v e r t i c e s [ j ]+1)] [−( v e r t i c e s [ j ]+1 ) ] ,d e c r ea s ing=TRUE) [ 1 : n ]

# Pegando as n−4 maiores co r r e l a ç õ e s g e r a i scor . maiores=sort (abs (

unique (c (cor . maiores . j ,cor . maiores . i ) ) ) ,d e c r ea s ing=TRUE) [1 : ( ( n/ 2 ) ) ]

cor . maiores . i=cor . maiores . i [abs ( cor . maiores . i )%in%cor . maiores ]

cor . maiores . j=cor . maiores . j [abs ( cor . maiores . j )

%in%cor . maiores ]

# Pegando os v é r t i c e s correspondentes as n maiores# co r r e l a çõ e s em cada v é r t i c e do par

for ( i 2 in 1 : length (cor . maiores . i ) )

pos i cao . i [ i 2 ]=which(A[ ,c ( k)]==cor . maiores . i [i 2 ] , useNames = TRUE)

for ( i 2 in 1 : length (cor . maiores . j ) )

pos i cao . j [ i 2 ]=which(A[ ,c ( v e r t i c e s [ j ]+1)]==cor . maiores . j [ i 2 ] ,useNames = TRUE)

# Pegando os v e r t i c e s v i z i n ho s mais co r r e l a c i onados# do par , assim como o propr io par

v i z=unique ( sort (c (k ,v e r t i c e s [ j ]+1 ,unique (c ( pos i cao . i ,pos i cao . j ) ) ) ) )

l c l<−length ( v i z )

Csub=covX [ viz , v i z ]c o r r . inv<−cor2pcor (Csub )

rownames( co r r . inv)=v i zcolnames ( co r r . inv)=v i z


co r r . inv . p<−1−pnorm(abs ( sqrt (n−( l c l −2)−3)∗log ((1+ cor r . inv [i1 , j 1 ] ) /(1− co r r . inv [ i1 , j 1 ] ) ) /2) )

Ap[ k , ( v e r t i c e s [ j ]+1)]=co r r . inv . p

Ap [ ( v e r t i c e s [ j ]+1) , k]=Ap[ k , ( v e r t i c e s [ j ]+1) ]

write . table (Ap, f i l e=paste (" pvalor_arvore_inout_local_" , i3 , "_" ,p , "_" ,

n , "_2.txt" , sep="" ) )

C.7 Correlação parcial local- Bayesiana

#### carregando b i b l i o t e c a s

require (MCMCpack)require ( corpcor )require ( igraph )require ( zoo )

##### funcoes

r iw i sh2=function (v , S )

return ( solve ( rwish2 (v , solve (S ) ) ) )

rwish2=function (v , S )

i f ( ! i s .matrix (S ) )S <− matrix (S)

i f (nrow(S) != ncol (S ) ) stop ( message = "S not square in rwish().\n" )

i f ( v < nrow(S ) )

C.7. CORRELAÇÃO PARCIAL LOCAL- BAYESIANA 63

stop ( message = "v is less than the dimension

of S in rwish().\n" )p <− nrow(S)CC <− chol (S , p ivot=TRUE)Z <− matrix (0 , p , p )diag (Z) <− sqrt ( rchisq (p , v : ( v − p + 1 ) ) )i f (p > 1)

pseq <− 1 : ( p − 1)Z [ rep (p ∗ pseq , pseq ) + unlist ( lapply ( pseq , seq ) )

] = rnorm(p ∗ (p − 1)/2)return ( crossprod (Z %∗% CC))

### funcao f b s t . area

f b s t . area=function (X, c t e )Y <− density ( subset (X, subset=(X<1)))

i f (min(Y$x ) >0| |max(Y$x)<0) return (0 )elsex=min(Y$x [Y$x>cte ] )i f ( x==In f ) x=max(Y$x [Y$x<cte ] )

#xt <− d i f f (Y$x [Y$y<Y$y [Y$x==x ] ] )#yt <− ro l lmean (Y$y [Y$y<Y$y [Y$x==x ] ] , 2 )xt <− Y$x [5]−Y$x [ 4 ]yt <− Y$y [Y$y<Y$y [Y$x==x ] ]

return (sum( xt∗yt ) )


#### ca lcu lando f b s t


p i l o t o=1

########lendo os dados

for (p in c (100 ) )for (n in c ( 50 ) )for ( i 3 in c ( 4 7 : 2 2 ) )#marca a hora do i n i c i oi n i c i o=Sys . time ( )

#lendo os dados da arvoreX=as .matrix ( read . table ( f i l e=paste ("arvore_inout_" , i3 ,

"_" ,p , ".txt" , sep="" ) ) ) [ 1 : n , ]

## se f o r t r a b a l h a r com a amostra p i l o t o para a p r i o r i

i f ( p i l o t o == 1)X. p=as .matrix ( read . table ( f i l e=paste (" arvore_piloto_inout_" , i3 , "_" ,p , "_.txt" , sep=""

) ) ) [ 1 : n , ]

#dando nomes as co lunas na arvorecolnames (X)=c ( seq ( 1 : p ) )

Xor i g i na l=XX. p o r i g i n a l=X. p##################

i f ( p i l o t o==1) covX=cov (X. p) else covX=cov (X)corX=cor (X)

#Pvalor no t r i a n gu l o supe r i o r da matrix e cor re l a cao# no t r i an gu l o i n f e r i o rpva lcor=cor . prob (X)

alpha=0.01

#Marcando os pares que sao s i gn i ca t emente co r r e l a c i onadosA=matrix ( , ncol=ncol (X) ,nrow=ncol (X) )for ( i in 1 : ( ncol (X)−1))for ( j in ( i +1) :(ncol (X) ) )A[ i , j ]= i f e l s e ( pva l cor [ i , j ]<alpha , corX [ i , j ] , 0 )A[ j , i ]=A[ i , j ]


diag (A)=0

#gra fo com as co r r e l a co e sg1=graph . adjacency (A,mode=("undirected" ) ,weighted=TRUE)#nomeando os v e r t i c e s do gra foV( g1 )$name <− c ( seq ( 1 : p ) )

############# Acessando os e l o s

##### bayes iana# Pegando , para cada v é r t i c e , os demais v e r t i c e s que tem# l i g a c ao com esse v é r t i c eedges=get . a d j l i s t ( g1 )

f b s t=matrix ( , ncol=p ,nrow=p)

Ap=matrix ( , ncol=p ,nrow=p)Af=matrix ( rep (0 , p^2) ,ncol=p ,nrow=p)

for ( k in 1 : p)

#v e r t i c e s que fazem l i g a c ao com o v e r t i c e kv e r t i c e s=edges [ [ k ] ]

i f ( length ( v e r t i c e s ) !=0)

i f ( length ( v e r t i c e s )==1)

f b s t [ k , ( v e r t i c e s +1)]= pva lcor [ k ,( v e r t i c e s +1)]

f b s t [ ( v e r t i c e s +1) ,k]= f b s t [ k , (v e r t i c e s +1)]

elsefor ( j in 1 : length ( v e r t i c e s ) )

# Condicao que c e r t i f i c a que so vamos pegar cada# par apenas uma única vez e so vamos c a l c u l a r a# corre l a ção p a r c i a l para pares s i gn i f i c an t emen t e# cor re l a c i onados

i f ( v e r t i c e s [ j ]+1>k && A[ k ,v e r t i c e s [ j ]+1] !=0)

X=Xor i g ina lX. p=X. p o r i g i n a l

i 1=as . character ( k )j 1=as . character ( v e r t i c e s [ j ]+1)


#Lis ta com todos os v i z i n ho s de cada v e r t i c e do par em#ques tao

v i z l i s t=graph . neighborhood (g1 , 1 , nodes=c ( i1 , j 1 ) ,mode=c ("all" ) )

#todos os v i z i n ho s dos do i s v e r t i c e s do par ,#inc l u indo e l e s .

v i z=sort ( as . integer (unique (c (V( v i z l i s t [ [ 1 ] ] ) $name ,V( v i z l i s t [ [ 2 ] ] ) $name ) ) ,d e c r ea s ing=TRUE) )

pos i cao . i=vector ( )pos i cao . j=vector ( )

#diminuindo a quant idade de#v i z inhos , caso tenhamos mais#v i z i n ho s do que observacoesi f ( length ( v i z)>=n)

#pegando as n maiores#co r r e l a co e s de cada#v e r t i c e do parcor . maiores . i=sort (A[

viz , c ( k ) ] [−k ] ,d e c r ea s ing=TRUE) [1 : n ]

cor . maiores . j=sort (A[viz , c ( v e r t i c e s [ j ]+1) ] [ − ( v e r t i c e s [ j ]+1 ) ] ,d e c r ea s ing=TRUE) [ 1 : n ]

#pegando as n−4 maiores#co r r e l a co e s g e r a i scor . maiores=sort (abs (

unique (c (cor . maiores . j ,cor . maiores . i ) ) ) ,d e c r ea s ing=TRUE) [1 : ( n/ 2 ) ]

cor . maiores . i=cor . maiores . i [ abs (cor . maiores . i ) %in%cor . maiores ]

cor . maiores . j=cor . maiores . j [ abs (cor . maiores . j ) %in%cor . maiores ]


# pegando os v e r t i c e s#correspondentes as n#maiores co r r e l a co e s em#cada v e r t i c e do parfor ( i 2 in 1 : length (

cor . maiores . i ) )pos i cao . i [ i 2 ]=

which(A[ , c ( k)]==cor . maiores . i [ i 2 ] ,useNames = TRUE)

for ( i 2 in 1 : length (cor . maiores . j ) )pos i cao . j [ i 2 ]=

which(A[ , c (v e r t i c e s [ j ]+1) ]==cor . maiores . j [i 2 ] , useNames =TRUE)

#Pegando os v e r t i c e s v i z i n ho s mais co r r e l a c i onados do#par , assim como o propr io par

v i z=unique ( sort (c (k ,v e r t i c e s [ j ]+1 ,unique (c ( pos i cao . i ,pos i cao . j ) ) ) ) )

#Pegando os dados s imulados apenas da v i z inhança#mais cor re l ac ionada do par

X=X[ , v i z ]X. p=X. p [ , v i z ]

l c l<−ncol (X)# A pr i o r i

mediaX=colMeans (X. p ,na .rm = FALSE, dims = 1)

lambda0=mediaXk0=l c l +3ps i 0=(k0−l c l −1)∗cov (X. p)v0=1

# A po s t e r i o r imediaX=colMeans (X,

na .rm = FALSE, dims = 1)k1=k0+nsoma=0for ( i in 1 : n )

soma=soma + (X[ i , ]−mediaX)%∗%t (X[ i , ]−


mediaX)ps i 1=ps i 0+soma+(n∗v0/ (n+v0 ) )∗

(mediaX−lambda0 )%∗%t (mediaX−lambda0 )

# Estimando cor re l a cao p a r c i a l usando l i s t a s

covpar=matrix ( )corpar=vector ( )

# Posição do v e r t i c e k na v i z inhança mais# corre l ac ionada

i 2=which( v i z==k)

# Posicao do outro v e r t i c e do par da i t e r a cao na# vi z inhanca mias cor re l ac ionada

j 2=which( v i z==v e r t i c e s [ j ]+1)

# Calculando a o paramentro da d i s t r i b u i c a o in v e r s e# wishar t da covar ianc ia p a r c i a l

p s i . i j=ps i 1 [ c ( i2 , j 2 ) , c ( i2 , j 2 ) ]−ps i 1 [ c ( i2 , j 2 ) ,−c ( i2 , j 2 )]%∗%solve ( p s i 1 [−c ( i2 , j 2 ) ,−c (i2 , j 2 ) ] )%∗%ps i 1 [−c ( i2 , j 2 ) ,c ( i2 , j 2 ) ]

i f ( i s . p o s i t i v e . d e f i n i t e (p s i . i j , method="chol")==FALSE)

p s i . i j=make . p o s i t i v e . d e f i n i t e (p s i . i j , t o l=1e−20)

covpar=matrix ( )

for ( k2 in 1 :1000)

i f ( i s . p o s i t i v e . d e f i n i t e (solve ( p s i . i j , t o l=1e−20) ,method="chol")==TRUE) covar=solve (rwish ( k1 , solve ( p s i . i j ,t o l=1e−20))) else

invcovar=rwish2 ( k1 ,

p s i . i j )while (i s . p o s i t i v e . d e f i n i t e (invcovar)==FALSE)invcovar=rwish2 ( k1 ,p s i . i j )

covar=solve ( invcovar )

C.8. MÉDIA DAS CURVAS ROC 69

#covpar=covar [ c ( i2 , j2 ) ,#c ( i2 , j2 )]− covar [ c ( i2 , j 2 ) ,#−c ( i2 , j2 ) ]%∗%s o l v e ( covar [#−c ( i2 , j2 ),−c ( i2 , j2 ) ] )%∗%#covar [−c ( i2 , j2 ) , c ( i2 , j 2 ) ]covpar=covari f ( covpar [1 ,1 ] <0) covpar [ 1 , 1

]=−covpar [ 1 , 1 ]i f ( covpar [2 ,2 ] <0) covpar [ 2 , 2

]=−covpar [ 2 , 2 ]

corpar [ k2 ]=( covpar [ 1 , 2 ] ) /sqrt ( covpar [ 1 , 1 ] ∗covpar [ 2 , 2 ] )

corpar=subset ( corpar , subset=(corpar !="NaN" ) )

#f b s t para H0: cor r e l a cao pa r c i a l = 0f b s t [ k , ( v e r t i c e s [ j ]+1)]=f b s t . area ( corpar , 0 )

f b s t [ ( v e r t i c e s [ j ]+1) , k]=f b s t [ k , ( v e r t i c e s [ j ]+1) ]

write . table ( fbs t , f i l e=paste ("fbst_local_arvore_inout_" ,i3 , "_" ,p , "_" ,n , "_2.txt" , sep="" ) )

f im=Sys . time ( )fim− i n i c i o

op ; l

C.8 Média das curvas ROC


####################### comeca aqui

require ( igraph )

# funcao que c a l c u l a a cor r e l a cao e o pva lo r na# mesma matr izcor . prob <− function (X, d f r = nrow(X) − 2)

R <− cor (X)above <− row(R) < col (R)r2 <− R[ above ]^2Fstat <− r2 ∗ d f r / (1 − r2 )R[ above ] <− 1 − pf ( Fstat , 1 , d f r )R

i n i c i o=Sys . time ( )

p=100for (n in c (1000) )s e n s i t i v i t y=l i s t ( )s p e c i f i c i t y=l i s t ( )

# indicadora se vamos t r a b a l h a r com e s t i s t i c a# c l a s s i c a − 1 − ou bayes iana − 0 −i f (n<=p) g l oba l=0 else g l oba l=1

for ( l in c ( 1 : 2 0 ) )

X=as .matrix ( read . table ( f i l e=paste ("arvore_inout_" , l , "_" ,p , ".txt" , sep="" ) ) ) [ 1 : n , ]

pva l cor=cor . prob (X)

corX=cor (X)

### CONFERIR O QUE MUDA COM ESSA ALTERACAOalpha=0.01

A=matrix (99 , ncol=p ,nrow=p)for ( i in 1 : ( p−1))for ( j in ( i +1):p )

A[ i , j ]= i f e l s e ( pva l cor [ i , j ]<alpha , corX [ i , j ] , 0 )A[ j , i ]=A[ i , j ]

diag (A)=0


i f ( c l a s s i c o==1)

i f ( g l oba l==1)pva lor=read . table (paste ("pvalor_arvore_inout_" , l ,

"_" ,p , "_" ,n , ".txt" , sep="" ) ) elsepva lor=read . table (paste ("pvalor_arvore_inout_local_

" , l , "_" ,p , "_" ,n , "_2.txt" , sep="" ) ) else

i f ( g l oba l==1)pva lor=read . table (paste ("

fbst_global_arvore_piloto_inout_" , l , "_" ,p ,"_" ,n , ".txt" , sep="" ) ) else

pva lor=read . table (paste (" fbst_local_arvore_inout_" , l , "_" ,p , "_" ,n ,

"_2.txt" , sep="" ) )

i f ( i s .na( pva lor [ 1 , 1 ] ) )

A1<−read . table (paste ("original_inout_" , l , "_" ,p , ".txt" ,sep="" ) )

s e n s i t i v i t y [ [ l ] ]= vector ( )s p e c i f i c i t y [ [ l ] ]= vector ( )

alpha=c ( seq ( 0 , 1 , 0 . 0 5 ) )

dados=matrix ( ,nrow=1,ncol=length ( alpha ) )colnames ( dados)=alpha

for ( k in 1 : length ( alpha ) )

Af=matrix (0 , ncol=p ,nrow=p)

for ( i in 1 : ( p−1))for ( j in ( i +1) :(p ) )i f ( ! i s .na( pva lor [ i , j ] ) & ( pva lor [ i , j ] < alpha [ k ] ) )

Af [ i , j ]=1Af [ j , i ]=1

A2=A∗AfA3=get . adjacency ( graph . adjacency (A2 ,weighted=TRUE) )

###### contando l i g a c o e s f a l s a s e ve rdade i ra s


po s i t i v o . real=sum(A1 [row(A1)<col (A1 ) ] )negat ivo . real=length (A1 [row(A1)<col (A1)])−

po s i t i v o . real

t rue . p o s i t i v e=sum(A1 [row(A1)<col (A1)]+A3 [row(A3)<col (A3)]==2)

true . negat ive=sum(A1 [row(A1)<col (A1)]+A3 [row(A3)<col (A3)]==0)

s e n s i t i v i t y [ [ l ] ] [ k]= true . p o s i t i v e/po s i t i v o . reals p e c i f i c i t y [ [ l ] ] [ k]= true . negat ive/negat ivo . real

## escrevendo a s e n s i b i l i d a d e e a e s p e c i f i c i d a d e de todas#v as arvores separadas

dados . sens=matrix ( ,nrow=1,ncol=length ( alpha ) )dados . spec=matrix ( ,nrow=1,ncol=length ( alpha ) )

dados . sens [1 , ]= s e n s i t i v i t y [ [ l ] ]dados . spec [1 , ]= s p e c i f i c i t y [ [ l ] ]

i f ( l==1)


i f ( g l oba l==1)

write . table ( dados . sens , paste (" sensitivity_inout_" ,p , "_" ,n ,

"_classico_global.txt" , sep="" ) ,col .names=alpha )

write . table ( dados . spec , paste (" specificity_inout_" ,p , "_" ,n ,

"_classico_global.txt" , sep="" ) ,col .names=alpha )

else


"_classico_local.txt" , sep="" ) ,col .names=alpha )



"_classico_local.txt" , sep="" ) ,col .names=alpha )

else

i f ( g l oba l==1)


"_bayesiano_global.txt" , sep="" ) ,col .names=alpha )


"_bayesiano_global.txt" , sep="" ) ,col .names=alpha )

else


"_bayesiano_local.txt" , sep="" ) ,col .names=alpha )


"_bayesiano_local.txt" , sep="" ) ,col .names=alpha )

else


i f ( g l oba l==1)


"_classico_global.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )


"_classico_global.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )

else


write . table ( dados . sens , paste ("

sensitivity_inout_" ,p , "_" ,n ,"_classico_local.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )


"_classico_local.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )

else

i f ( g l oba l==1)


"_bayesiano_global.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )


"_bayesiano_global.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )

else


"_bayesiano_local.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )


"_bayesiano_local.txt" , sep="" ) ,append=TRUE, col .names=FALSE,row .names=l )

mean . s ens=vector ( )mean . spec=vector ( )


median . s ens=vector ( )median . spec=vector ( )

q1 . sens=vector ( )q1 . spec=vector ( )

q3 . sens=vector ( )q3 . spec=vector ( )

for ( k in 1 : length ( alpha ) )

mean . s ens [ k]=mean( sapply ( s e n s i t i v i t y , function ( x )x [ k ] ) )

mean . spec [ k]=mean( sapply ( s p e c i f i c i t y , function ( x )x [ k ] ) )

median . s ens [ k]=median( sapply ( s e n s i t i v i t y , function ( x )x [ k ] ) )

median . spec [ k]=median( sapply ( s p e c i f i c i t y , function ( x )x [ k ] ) )

q1 . sens [ k]=quantile ( sapply ( s e n s i t i v i t y , function ( x )x [ k ] ) , 0 . 2 5 )

q1 . spec [ k]=quantile ( sapply ( s p e c i f i c i t y , function ( x )x [ k ] ) , 0 . 2 5 )

q3 . sens [ k]=quantile ( sapply ( s e n s i t i v i t y , function ( x )x [ k ] ) , 0 . 7 5 )

q3 . spec [ k]=quantile ( sapply ( s p e c i f i c i t y , function ( x )x [ k ] ) , 0 . 7 5 )


i f ( g l oba l==1)write . table (data . frame ( alpha ,mean . sens ,mean . spec ,median . sens ,median . spec , q1 . sens , q1 . spec ,q3 . sens , q3 . spec ) , f i l e=paste ("

media_quantis_sens_spec_inout_" ,p , "_" ,n ,"_classico_global.txt" , sep="" ) ) else

write . table (data . frame ( alpha ,mean . sens ,mean . spec ,median . sens ,median . spec , q1 . sens , q1 . spec ,q3 . sens , q3 . spec ) , f i l e=paste ("

media_quantis_sens_spec_inout_" ,p , "_" ,n ,"_2_classico_local.txt" , sep="" ) )


else

i f ( g l oba l==1)write . table (data . frame ( alpha ,mean . sens ,mean . spec ,median . sens ,median . spec , q1 . sens , q1 . spec ,q3 . sens , q3 . spec ) , f i l e=paste ("

media_quantis_sens_spec_inout_" ,p , "_" ,n ,"_bayesiano_global.txt" , sep="" ) ) else

write . table (data . frame ( alpha ,mean . sens ,mean . spec ,median . sens ,median . spec , q1 . sens , q1 . spec , q3 . sens ,q3 . spec ) , f i l e=paste ("media_quantis_sens_spec_inout_" ,p , "_" ,n , "_bayesiano_local.txt" , sep="" ) )

fim=Sys . time ( )fim− i n i c i o

Referências Bibliográcas

[1] K. Baba, R. Shibata, and M. Sibuya. Partial correlation and conditionalcorrelation as measures of conditional independence. Australian and

New Zealand Journal of Statistics, 46:657664, 2004.

[2] R. Bhatia. Positive Denite Matrices. Princeton University Press, pri-meira edition, 2007.

[3] A.P. Dawid. Conditional independence in statistical theory (with dis-cussion). J. R. Stat. Soc., B 41:131, 1979.

[4] M.H. DeGroot and M.J. Schervish. Probability and Statistics. Addison-Wesley, 3 edition, 2002.

[5] D. Edwards. Introduction to Graphical Modelling. Springer-Verlag NewYork, Inc., 1 edition, 1995.

[6] R.A. Fisher. The distribution of the partial correlation coecien. Me-

tron, 3:329332, 1924.

[7] J.J. Goeman and P. Bühlmann. Analyzing gene expression data in termsof gene sets: methodological issues. Bioinformatics, 23:980987, 2007.

[8] L. Han and J. Zhu. Using matrix of thresholding partial correlationcoecients to infer regulatory network. BioSystems, 91:158165, 2008.

[9] M.G. Kendal and A. Stuart. The Advanced Theory of Statistics, Volume

2. nao sei, 3 edition, 1973.

[10] K. V. Mardia, J. T. Kent, and J. M. Bibby. Multivariate Analysis.Academic Press, London, UK, 10 edition, 1979.

[11] K.L. Mine, N. Shulzhenko, A. Yambartsev, G.F.O. Sanson, S. Varma,Volfovsky N., S.M.F. Brenna, C.R.N. Carvalho, J.C.L. Ribalta, J. Skin-ner, H. Lyng, I.D.C.G. Silva, M. Gerbase-DeLima, and A. Morgun. Re-construction of an integrative gene regulatory meta-network reveals cellcycle and antiviral response as major drivers of cervical cancer. under

review: GENOME Research, 2012.

77

78 REFERÊNCIAS BIBLIOGRÁFICAS

[12] M.E.J Newman. The structure and function of complex networks. SIAMRev., 45:167256, 2003.

[13] J. Peng, W. Pei, N. Zhou, and J. Zhu. Partial correlation estimationby joint sparse regression models. J. Am. Stat. Assoc.x, 104:735746,2009.

[14] C. A. B. Pereira and J. M. Stern. Evidence and credibility: Full bayesiansignicance test for precise hypotheses. Entropy, 1:99110, 1999.

[15] S.J. Press. Applied multivariate analysis: using Bayesian and frequentist

methods of inference. Dover Publications, Inc. Mineola, New York, 2edition, 1982.

[16] H. Raia and R. Schlaifer. Applied Statistical Decision Theory. HarvardPress, Boston, 1 edition, 1961.

[17] J. L. Rodgers and W. A. Nicewander. Thirteen ways to look at corre-lation coecients. American Statisticians, 42:5966, 1988.

[18] T. Shimamura, S. Imoto, R. Yamaguchi, and S. Miyano. Weighted lassoin graphical gaussian modeling for large gene network estimation basedon microarray data. Genome Inform, 19:142153, 2007.

[19] D.Y. Takahashi. Medidas de Fluxo de Informação com Aplicação em

Neurociência. PhD thesis, Programa Interunidades em Bioinformática,Universidade de São Paulo, Brasil, 2008.

[20] H. Toh and K. Horimoto. Inference of a genetic network by a combinedapproach of cluster analysis and graphical gaussian modeling. Bioinfor-matics, 18:287297, 2002.

[21] J. Wang, O. Myklebost, and E. Hovig. Mgraph: graphical models formicroarray data analysis. Bioinformatics, 19:22102211, 2003.

[22] X. Wu and K.R. Subramanian Y. Ye. Interactive analysis of gene inte-ractions using graphical gaussian model. In Proceedings of the 3rd ACM

SIGKDD Workshop on Data Mining in Bioinformatics. Brighton, UK,2003.

Documents

Construção de redes usando estatística clássica e Bayesiana-uma