Bastão de Asclépio & Distribuição Normal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
options(warn=-1)
suppressMessages(library(aplpack, warn.conflicts=FALSE))
suppressMessages(library(calculus, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(cellWise, warn.conflicts=FALSE))
suppressMessages(library(cluster, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(diptest, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(ellipse, warn.conflicts=FALSE))
suppressMessages(library(far, warn.conflicts=FALSE))
suppressMessages(library(fmsb, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(ggfortify, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(HH, warn.conflicts=FALSE))
suppressMessages(library(HSAUR2, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lawstat, warn.conflicts=FALSE))
suppressMessages(library(lfda, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(MatchIt, warn.conflicts=FALSE)) # distancia robusta
suppressMessages(library(matlib, warn.conflicts=FALSE))
suppressMessages(library(matrixcalc, warn.conflicts=FALSE))
suppressMessages(library(MatrixModels, warn.conflicts=FALSE))
suppressMessages(library(matrixStats, warn.conflicts=FALSE))
suppressMessages(library(MomTrunc, warn.conflicts=FALSE))
suppressMessages(library(MVA, warn.conflicts=FALSE)) # An Introduction to Applied Multivariate Analysis with R
suppressMessages(library(MVar.pt, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(mvtnorm, warn.conflicts=FALSE))
suppressMessages(library(plotly, warn.conflicts=FALSE))
suppressMessages(library(pracma, warn.conflicts=FALSE))
suppressMessages(library(rayshader, warn.conflicts=FALSE))
suppressMessages(library(rcompanion, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(rgl, warn.conflicts=FALSE))
suppressMessages(library(scatterplot3d, warn.conflicts=FALSE))
# suppressMessages(library(MVLM, warn.conflicts=FALSE))
# suppressMessages(library(reticulate, warn.conflicts=FALSE))
source("summarySEwithin2.R")
Bastão de Asclépio & Distribuição Normal.
Simple Tools for Examining and Cleaning Dirty Data:
janitor
Trocar \bar
por \overline
Efficient Programming in R - Tsagris - 2025: Rfast
e
Rfast2
Valero-Mora, P., Rodrigo, M. F., Sanchez, M. & SanMartin, J., (2019) “A Plot for the Visualization of Missing Value Patterns in Multivariate Data”, Practical Assessment, Research, and Evaluation 24(1): 9. doi: https://doi.org/10.7275/94ra-1y55
Unwin, A (2015) Graphical data analysis with R. USA: CRC.
Linear Algebra and Its Applications with R - Yoshida - 2021.pdf
library(conics)
library(factoextra)
RPubs
“Não é paradoxo dizer que nos nossos momentos de inspiração mais teórica podemos estar o mais próximo possível de nossas aplicações mais práticas.”
WHITEHEAD, AN apud BOYER, CB (1974) História da Matemática. São Paulo: Blücher/EDUSP, p. 419.
I. GETTING STARTED
Applications of Multivariate Techniques. The Organization of Data. Data Displays and Pictorial Representations. Distance. Final Comments.
Some Basics of Matrix and Vector Algebra. Positive Definite Matrices. A Square-Root Matrix. Random Vectors and Matrices. Mean Vectors and Covariance Matrices. Matrix Inequalities and Maximization. Supplement 2A Vectors and Matrices: Basic Concepts.
The Geometry of the Sample. Random Samples and the Expected Values of the Sample Mean and Covariance Matrix. Generalized Variance. Sample Mean, Covariance, and Correlation as Matrix Operations. Sample Values of Linear Combinations of Variables.
The Multivariate Normal Density and Its Properties. Sampling from a Multivariate Normal Distribution and Maximum Likelihood Estimation. The Sampling Distribution of \(\bar{X}\) and \(S\). Large-Sample Behavior of \(\bar{X}\) and \(S\). Assessing the Assumption of Normality. Detecting Outliners and Data Cleaning. Transformations to Near Normality.
The Plausibility of \(\mu_0\) as a Value for a Normal Population Mean. Hotelling’s \(T^2\) and Likelihood Ratio Tests. Confidence Regions and Simultaneous Comparisons of Component Means. Large Sample Inferences about a Population Mean Vector. Multivariate Quality Control Charts. Inferences about Mean Vectors When Some Observations Are Missing. Difficulties Due To Time Dependence in Multivariate Observations. Supplement 5A Simultaneous Confidence Intervals and Ellipses as Shadows of the \(p\)-Dimensional Ellipsoids.
Paired Comparisons and a Repeated Measures Design. Comparing Mean Vectors from Two Populations. Comparison of Several Multivariate Population Means (One-Way MANOVA). Simultaneous Confidence Intervals for Treatment Effects. Two-Way Multivariate Analysis of Variance. Profile Analysis. Repeated Measures, Designs, and Growth Curves. Perspectives and a Strategy for Analyzing Multivariate Models.
The Classical Linear Regression Model. Least Squares Estimation. Inferences About the Regression Model. Inferences from the Estimated Regression Function. Model Checking and Other Aspects of Regression. Multivariate Multiple Regression. The Concept of Linear Regression. Comparing the Two Formulations of the Regression Model. Multiple Regression Models with Time Dependant Errors. Supplement 7A The Distribution of the Likelihood Ratio for the Multivariate Regression Model.
Population Principal Components. Summarizing Sample Variation by Principal Components. Graphing the Principal Components. Large-Sample Inferences. Monitoring Quality with Principal Components. Supplement 8A The Geometry of the Sample Principal Component Approximation.
The Orthogonal Factor Model. Methods of Estimation. Factor Rotation. Factor Scores. Perspectives and a Strategy for Factor Analysis. Structural Equation Models. Supplement 9A Some Computational Details for Maximum Likelihood Estimation.
Canonical Variates and Canonical Correlations. Interpreting the Population Canonical Variables. The Sample Canonical Variates and Sample Canonical Correlations. Additional Sample Descriptive Measures. Large Sample Inferences.
Separation and Classification for Two Populations. Classifications with Two Multivariate Normal Populations. Evaluating Classification Functions. Fisher’s Discriminant Function. Separation of Populations. Classification with Several Populations. Fisher’s Method for Discriminating among Several Populations. Final Comments.
Similarity Measures. Hierarchical Clustering Methods. Nonhierarchical Clustering Methods. Multidimensional Scaling. Correspondence Analysis. Biplots for Viewing Sample Units and Variables. Procustes Analysis: A Method for Comparing Configurations.
Multivariate statistics: Wikipedia
A análise de dados, embora interessante com uma variável, torna-se verdadeiramente fascinante e desafiadora quando várias variáveis estão envolvidas.
Pesquisadores nas ciências biológicas, físicas e sociais frequentemente coletam medições de diversas variáveis.
Pacotes de computador modernos fornecem prontamente os resultados numéricos para análises estatísticas bastante complexas.
Fornecemos aos estudantes o conhecimento necessário para fazer interpretações adequadas, selecionar técnicas estatísticas multivariadas apropriadas e compreender seus pontos fortes e fracos.
Esperamos que nossas discussões atendam às necessidades de cientistas, em uma ampla variedade de áreas temáticas, como uma introdução legível à análise estatística de observações multivariadas.
Nosso objetivo é apresentar os conceitos e métodos de análise multivariada em um nível que seja facilmente compreensível por leitores que fizeram dois ou mais cursos de estatística.
Enfatizamos as aplicações de métodos multivariados e, conseqüentemente, tentamos tornar a matemática o mais palatável possível.
Evitamos o uso de cálculo diferencial matricial [sic: não deve ser evitado na formação de estatístico aplicado!]. Por outro lado, os conceitos de matriz e de manipulação de matrizes são importantes.
Não assumimos que o leitor esteja familiarizado com a álgebra matricial.
Em vez disso, apresentamos as matrizes conforme elas aparecem naturalmente em nossas discussões e, em seguida, mostramos como elas simplificam a apresentação de modelos e técnicas multivariadas.
A investigação científica é um processo de aprendizagem iterativo. Os objetivos relacionados à explicação de um fenômeno social ou físico devem ser especificados e testados por meio da coleta e análise de dados. A análise dos dados obtidos por meio de experimentação ou observação muitas vezes sugere uma explicação modificada do fenômeno. Ao longo desse processo iterativo, as variáveis são frequentemente adicionadas ou excluídas do estudo. Portanto, investigar a maioria dos fenômenos requer coletar observações sobre várias variáveis. Este livro se concentra em métodos estatísticos projetados para extrair informações desses tipos de conjuntos de dados. Uma vez que os dados envolvem medições simultâneas em múltiplas variáveis, esta coleção de metodologia é conhecida como análise multivariada.
Compreender as relações entre múltiplas variáveis torna a análise multivariada um assunto desafiador. O grande volume de dados muitas vezes sobrecarrega a mente humana. Além disso, desenvolver técnicas estatísticas multivariadas para fazer inferências requer mais matemática em comparação com uma configuração univariada. Neste livro, fornecemos explicações baseadas em conceitos algébricos e evitamos derivações estatísticas que envolvam cálculo diferencial com múltiplas variáveis [sic: não deve ser evitado na formação de estatístico aplicado!]. Nosso objetivo é apresentar várias técnicas multivariadas úteis de maneira clara, contando fortemente com exemplos ilustrativos e minimizando a complexidade matemática. No entanto, um certo nível de proficiência matemática e um desejo de pensamento quantitativo são necessários.
A maior parte do nosso foco será analisar as medições obtidas sem controlar ou manipular ativamente nenhuma das variáveis que estão sendo medidas. Somente nos Capítulos 6 e 7, discutimos alguns delineamentos experimentais (projetos) que envolvem a manipulação ativa de variáveis importantes para gerar dados. Embora o planejamento experimental seja normalmente uma parte crucial de uma investigação científica, muitas vezes é impossível ter controle total sobre a geração de dados apropriados em certas disciplinas, como negócios, economia, ecologia, geologia e sociologia.
Muitos métodos multivariados são baseados em um modelo de probabilidade conhecido como distribuição normal multivariada. Outros métodos são de natureza ad hoc e são justificados por argumentos lógicos ou de senso comum. Independentemente de sua origem, as técnicas multivariadas devem ser implementadas em um computador. Avanços recentes na tecnologia de computadores levaram ao desenvolvimento de linguagens e softwares estatísticos sofisticados, tornando a etapa de implementação mais fácil.
A análise multivariada é um campo diversificado e é difícil estabelecer um esquema de classificação amplamente aceito que indique a adequação das técnicas. Uma classificação distingue técnicas projetadas para estudar relacionamentos interdependentes daquelas projetadas para estudar relacionamentos dependentes. Outra classificação considera o número de populações e conjuntos de variáveis em estudo. Neste texto, os capítulos são divididos em seções que enfocam a inferência sobre meios de tratamento, inferência sobre estrutura de covariância e técnicas de classificação ou agrupamento. No entanto, é importante notar que a escolha dos métodos e tipos de análises empregadas são em grande parte determinadas pelos objetivos da investigação. Na Seção 1.2, fornecemos uma lista de problemas práticos destinados a ilustrar a conexão entre a escolha de um método estatístico e os objetivos do estudo. Esses problemas, juntamente com os exemplos no texto, irão ajudá-lo a avaliar a aplicabilidade de técnicas multivariadas em diferentes campos.
Os objetivos das investigações científicas aos quais os métodos multivariados se prestam mais naturalmente incluem o seguinte:
Concluímos esta breve visão geral da análise multivariada com uma citação de F.H.C. Marriott. Embora a afirmação tenha sido feita no contexto da análise de cluster, acreditamos que ela é aplicável a uma gama mais ampla de métodos. É essencial manter essa perspectiva em mente sempre que conduzir ou ler sobre análise de dados. Permite-nos manter uma perspectiva adequada e não nos deixarmos dominar pela elegância de algumas teorias:
“Se os resultados discordam da opinião informada, não admitem uma interpretação lógica simples e não aparecem claramente em uma apresentação gráfica, provavelmente estão errados. Não há mágica nos métodos numéricos e há muitas maneiras pelas quais eles são uma ajuda valiosa para a interpretação dos dados, não máquinas de salsicha que transformam automaticamente conjuntos de números em pacotes de fatos científicos.”
Ao longo deste texto, estaremos preocupados em analisar medições feitas em várias variáveis ou características.
Essas medições (comumente chamadas de dados) frequentemente precisam ser organizadas e exibidas de várias maneiras. Por exemplo, gráficos e tabelas são auxílios importantes na análise de dados. Números resumidos, que retratam quantitativamente certos aspectos dos dados, também são necessários em qualquer descrição.
Agora introduzimos os conceitos preliminares que fundamentam esses primeiros passos da organização de dados.
Dados multivariados surgem sempre que um pesquisador, buscando entender um fenômeno social ou físico, seleciona um número \(p\ge2\) de variáveis ou características manifestas para registrar. Os valores dessas variáveis são registrados para cada item, indivíduo ou unidade experimental diferente.
Usaremos a notação \(x_{jk}\) para indicar o valor específico da variável \(j\)-ésima observada no \(k\)-ésimo item ou unidade experimental. Ou seja, \(x_{jk}\) = medida da \(j\)-ésima variável no \(k\)-ésimo item.
Consequentemente, \(n\) medições em \(p\) variáveis podem ser exibidas da seguinte forma:
\[ \begin{array}{cccccc} \text{Item} & \text{Variável 1} & \text{Variável 2} & \cdots & \text{Variável}\;k &\cdots & \text{Variável} \;p \\ \hline 1 & x_{11} & x_{12} & \cdots & x_{1k} & \cdots & x_{1p}\\ 2 & x_{21} & x_{22} & \cdots & x_{2k} & \cdots & x_{2p}\\ \vdots & \vdots & \vdots & & \vdots & & \vdots\\ j & x_{j1} & x_{j2} & \cdots & x_{jk} & \cdots & x_{jp}\\ \vdots & \vdots & \vdots & & \vdots & & \vdots\\ n & x_{n1} & x_{n2} & \cdots & x_{nk} & \cdots & x_{np}\\ \end{array} \]
Podemos representar esses dados como uma matriz retangular, denotada como \(\boldsymbol{x}\), com \(n\) linhas e \(p\) colunas:
\[ \boldsymbol{\mathcal{x}} = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1k} & \cdots & x_{1p}\\ x_{21} & x_{22} & \cdots & x_{2k} & \cdots & x_{2p}\\ \vdots & \vdots & & \vdots & & \vdots\\ x_{j1} & x_{j2} & \cdots & x_{jk} & \cdots & x_{jp}\\ \vdots & \vdots & & \vdots & & \vdots\\ x_{n1} & x_{n2} & \cdots & x_{nk} & \cdots & x_{np}\\ \end{bmatrix} \]
A matriz \(\boldsymbol{\mathcal{x}}\) engloba todas as observações sobre todas as variáveis.
Uma seleção de quatro recibos de uma livraria universitária foi obtida para investigar a natureza das vendas de livros. Cada recibo fornecia, entre outras coisas, o número de livros vendidos e o valor total de cada venda. Vamos considerar a primeira variável como vendas totais em dólares (\(x_{1}\)) e a segunda variável como o número de livros vendidos (\(x_{2}\)). Então, podemos considerar os números correspondentes nos recibos como quatro medições em duas variáveis. Suponha que os dados, em forma tabular, sejam os seguintes:
\[ \begin{align} \text{Variável 1 (venda em dólar):} & \quad 42 \quad 52 \quad 48 \quad 58 \\ \text{Variável 2 (número de livros):} & \quad 4 \quad 5 \quad 4 \quad 3 \\ \end{align} \]
Usando a notação acabamos de introduzir, temos:
\[ \begin{align} x_{11} = 42, & \quad x_{21} = 52, & \quad x_{31} = 48, & \quad x_{41} = 58 \\ x_{12} = 4, & \quad x_{22} = 5, & \quad x_{32} = 4, & \quad x_{42} = 3 \\ \end{align} \]
E a matriz de dados \(\boldsymbol{\mathcal{x}}\) é:
\[ \boldsymbol{\mathcal{x}} = \begin{bmatrix} 42 & 4 \\ 52 & 5 \\ 48 & 4 \\ 58 & 3 \\ \end{bmatrix} \]
com quatro linhas e duas colunas.
Considerar os dados na forma de matrizes facilita a exposição do assunto e permite que cálculos numéricos sejam realizados de maneira organizada e eficiente. A eficiência é dupla, pois há ganhos tanto na descrição de cálculos numéricos como operações em matrizes quanto na implementação dos cálculos em computadores, que agora utilizam várias linguagens e pacotes estatísticos para realizar operações em matrizes. Abordaremos a manipulação de matrizes de números no Capítulo 2. Neste ponto, estamos preocupados apenas com o valor delas como dispositivos para exibir dados.
Um conjunto de dados grande é volumoso e sua massa representa um obstáculo sério para qualquer tentativa de extrair informações pertinentes visualmente. Grande parte das informações contidas nos dados pode ser avaliada calculando certos números resumidos, conhecidos como estatísticas descritivas. Por exemplo, a média aritmética, ou média amostral, é uma estatística descritiva que fornece uma medida de localização, ou seja, um “valor central” para um conjunto de números. E a média dos quadrados das distâncias de todos os números em relação à média fornece uma medida da dispersão ou variação dos números.
Vamos nos basear principalmente em estatísticas descritivas que medem localização, variação e associação linear. As definições formais dessas quantidades são as seguintes.
Seja \(x_{1k}, x_{2k}, \ldots, x_{nk}\) as \(n\) medições da \(k\)-ésima variável, \(k=1, 2, \ldots, p\).
Média amostral:
\[ \bar{x}_k=\dfrac{1}{n}\sum_{j=1}^{n}{x_{jk}}\tag{1-1} \]
Variância amostral:
\[ s_{k}^{2}=s_{kk}=\dfrac{1}{n}\sum_{j=1}^{n}{\left(x_{jk} - \bar{x}_k \right)^2}\tag{1-3} \]
Covariância amostral:
\[ s_{ik}=\dfrac{1}{n}\sum_{j=1}^{n}{\left(x_{ji} - \bar{x}_i \right)\left(x_{jk} - \bar{x}_k \right)}\tag{1-4} \]
Coeficiente de correlação de Pearson:
\[ r_{ik}=\dfrac{s_{ik}}{\sqrt{s_{ii}}\sqrt{s_{kk}}}\in[-1,1]\tag{1-5} \]
Valor padronizado:
\[ z_{jk}=\dfrac{x_{jk}-\bar{x}_{k}}{\sqrt{s_{kk}}} \]
Portanto, a covariância de valores padronizados é a correlação dos valores brutos:
\[ s_{ik}=\dfrac{1}{n}\sum_{j=1}^{n}{\left(z_{ji} - \bar{z}_i \right)\left(z_{jk} - \bar{z}_k \right)}=\dfrac{1}{n}\sum_{j=1}^{n}{z_{ji}z_{jk}}=r_{ik} \]
Correlação de Pearson das transformações lineares das variáveis é a correlação de Pearson das variáveis brutas:
\[ y_{ji}=ax_{ji}+b\\ y_{jk}=cx_{jk}+d\\\\ r\left(y_{ji},y_{jk}\right)=r\left(x_{ji},x_{jk}\right)=r_{ik} \]
Matricialmente:
Vetor de médias amostrais:
\[ \bar{\mathbf{x}}= \begin{bmatrix} \overline{x}_1 \\ \overline{x}_2 \\ \vdots \\ \overline{x}_p \end{bmatrix} \]
Matriz de variâncias e covariâncias amostrais:
\[ \mathbf{s}_n= \begin{bmatrix} s_{11} & s_{12} & \cdots & s_{1p} \\ s_{21} & s_{22} & \cdots & s_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ s_{p1} & s_{p2} & \cdots & s_{pp} \end{bmatrix} \]
Matriz de correlações amostrais:
\[ \mathbf{r}= \begin{bmatrix} 1 & r_{12} & \cdots & r_{1p} \\ r_{21} & 1 & \cdots & r_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ r_{p1} & r_{p2} & \cdots & 1 \end{bmatrix} \]
Considere os dados introduzidos no Exemplo 1.1. Cada recibo fornece um par de medições, vendas totais em dólares e número de livros vendidos. Encontre as matrizes \(\bar{x}\), \(s_n\) e \(r\).
Como há quatro recibos, temos um total de quatro medições (observações) em cada variável.
As médias amostrais são:
\[ \begin{align} \bar{x}_1 &= \dfrac{1}{4}\sum_{j=1}^{4}{x_{j1}}=\dfrac{1}{4}(42 + 52 + 48 + 58) = 50 \\ \bar{x}_2 &= \dfrac{1}{4}\sum_{j=1}^{4}{x_{j2}}=\dfrac{1}{4}(4 + 5 + 4 + 3) = 4 \end{align} \]
\[ \bar{\mathbf{x}} = \begin{bmatrix} \bar{x}_1 \\ \bar{x}_2 \\ \end{bmatrix}= \begin{bmatrix} 50 \\ 4 \\ \end{bmatrix} \]
As variâncias amostrais e as covariâncias são:
\[ \begin{align} s_{11} &= \dfrac{1}{4}\sum_{j=1}^{4}{\left(x_{j1} - \bar{x}_1 \right)^2}\\ &=\dfrac{1}{4}\left((42 - 50)^2 + (52 - 50)^2 + (48 - 50)^2 + (58 - 50)^2\right) \\ s_{11}&= 34 \\\\ s_{22} &= \dfrac{1}{4}\sum_{j=1}^{4}{\left(x_{j2} - \bar{x}_2 \right)^2}\\ &=\dfrac{1}{4}\left((4 - 4)^2 + (5 - 4)^2 + (4 - 4)^2 + (3 - 4)^2\right)\\ s_{22} &= 0.5 \\\\ s_{12} &= \dfrac{1}{4}\sum_{j=1}^{4}{\left(x_{j1} - \bar{x}_1 \right)\left(x_{j2} - \bar{x}_2 \right)}\\ &=\dfrac{1}{4}\left((42 - 50)(4 - 4) + (52 - 50)(5 - 4) + (48 - 50)(4 - 4) + (58 - 50)(3 - 4)\right) \\ s_{12}&= -1.5 \\\\ s_{21} &= s_{12} = -1.5 \end{align} \]
Então:
\[ \begin{align} \mathbf{s}_n &= \left[\begin{array}{rr} 34 & -1.5 \\ -1.5 & 0.5 \end{array}\right] \\ \mathbf{r} &= \left[\begin{array}{rr} 1 & -0.36 \\ -0.36 & 1 \end{array}\right] \end{align} \]
Gráficos são auxílios importantes, mas frequentemente negligenciados, na análise de dados.
Embora seja impossível plotar simultaneamente todas as medições feitas em várias variáveis e estudar as configurações, gráficos de variáveis individuais e gráficos de pares de variáveis ainda podem ser muito informativos. Programas de computador sofisticados e equipamentos de exibição permitem o luxo de examinar visualmente os dados em uma, duas ou três dimensões com relativa facilidade. Por outro lado, muitas percepções valiosas podem ser obtidas a partir dos dados construindo gráficos. É uma boa prática estatística plotar pares de variáveis e inspecionar visualmente o padrão de associação.
O papel é fabricado em folhas contínuas com vários pés de largura. Devido à orientação das fibras dentro do papel, ele tem uma resistência diferente quando medido na direção produzida pela máquina em comparação com quando medido em direção cruzada, ou seja, em ângulo reto com a direção da máquina. A Tabela 1.2 mostra os valores medidos de:
Uma apresentação gráfica inovadora desses dados aparece na Figura 1.5. Os gráficos de dispersão são organizados como elementos fora da diagonal principal de uma matriz de covariância e os boxplots como elementos da diagonal principal. Estes últimos estão em uma escala diferente com este software, então usamos apenas a forma geral para fornecer informações sobre simetria e possíveis valores discrepantes para cada característica individual. Os gráficos de dispersão podem ser analisados em busca de padrões e observações incomuns. Na Figura 1.5, há uma observação incomum: a densidade do espécime 25. Alguns dos gráficos de dispersão apresentam padrões que sugerem a presença de dois grupos separados de observações. Esses arranjos de gráficos de dispersão são discutidos em maior detalhe na nossa próxima seção sobre gráficos de software.
Espécime | Densidade | Direção da Máquina | Direção Cruzada |
---|---|---|---|
1 | 0.801 | 121.41 | 70.42 |
2 | 0.824 | 127.70 | 72.47 |
3 | 0.841 | 129.20 | 78.20 |
4 | 0.816 | 131.80 | 74.89 |
5 | 0.840 | 135.10 | 71.21 |
6 | 0.842 | 131.50 | 78.39 |
7 | 0.820 | 126.70 | 69.02 |
8 | 0.802 | 115.10 | 73.10 |
9 | 0.828 | 130.80 | 79.28 |
10 | 0.819 | 124.60 | 76.48 |
11 | 0.826 | 118.31 | 70.25 |
12 | 0.802 | 114.20 | 72.88 |
13 | 0.810 | 120.30 | 68.23 |
14 | 0.802 | 115.70 | 68.12 |
15 | 0.832 | 117.51 | 71.62 |
16 | 0.796 | 109.81 | 53.10 |
17 | 0.759 | 109.10 | 50.85 |
18 | 0.770 | 115.10 | 51.68 |
19 | 0.759 | 118.31 | 50.60 |
20 | 0.772 | 112.60 | 53.51 |
21 | 0.806 | 116.20 | 56.53 |
22 | 0.803 | 118.00 | 70.70 |
23 | 0.845 | 131.00 | 74.35 |
24 | 0.822 | 125.70 | 68.29 |
25 | 0.971 | 126.10 | 72.10 |
26 | 0.816 | 125.80 | 70.64 |
27 | 0.836 | 125.50 | 76.33 |
28 | 0.815 | 127.80 | 76.75 |
29 | 0.822 | 130.50 | 80.33 |
30 | 0.822 | 127.90 | 75.68 |
31 | 0.843 | 123.90 | 78.54 |
32 | 0.824 | 124.10 | 71.91 |
33 | 0.788 | 120.80 | 68.22 |
34 | 0.782 | 107.40 | 54.42 |
35 | 0.795 | 120.70 | 70.41 |
36 | 0.805 | 121.91 | 73.68 |
37 | 0.836 | 122.31 | 74.93 |
38 | 0.788 | 110.60 | 53.52 |
39 | 0.772 | 103.51 | 48.93 |
40 | 0.776 | 110.71 | 53.67 |
41 | 0.758 | 113.80 | 52.42 |
T1.2 <- read.table("JW6Data/T1-2.DAT", quote="\"", comment.char="")
colnames(T1.2) <- c("Densidade",
"Direção_Máquina",
"Direção_Cruzada")
Espécime <- 1:41
df <- cbind(Espécime,
T1.2)
#knitr::kable(head(df), caption="Medições de resistência do papel")
#knitr::kable(tail(df), caption="Medições de resistência do papel")
car::scatterplotMatrix(df[,-1],
diagonal=list(method="boxplot"),
smooth=FALSE,
regLine=FALSE,
upper.panel=NULL,
col="gray65",
bg="black",
cex=1,
pch=1,
main="Figura 1.5: Gráficos de dispersão e boxplots dos dados de \nqualidade do papel da Tabela 1.2.")
GGally::ggpairs(df[, -1],
title = "Gráficos de Dispersão com Densidade",
diag = list(continuous = "densityDiag"))
Um zoólogo obteve medições em \(n = 25\) lagartos conhecidos cientificamente como Cophosaurus texanus. A massa é dada em grama, enquanto o comprimento do focinho até a abertura ventral (SVL) e a envergadura dos membros posteriores (HLS) são dados em milímetro. Os dados são exibidos na Tabela 1.3.
Embora haja três medidas de tamanho, podemos questionar se a maior parte da variação está principalmente restrita a duas dimensões ou até mesmo a uma dimensão. Para ajudar a responder perguntas sobre a dimensionalidade reduzida, construímos o gráfico de dispersão tridimensional na Figura 1.6. Claramente, a maior parte da variação é dispersa em torno de uma linha reta unidimensional. Saber a posição em uma linha ao longo dos principais eixos do conjunto de pontos seria quase tão bom quanto saber as três medidas Massa, SVL e HLS.
No entanto, esse tipo de análise pode ser enganoso se uma variável tiver uma variância muito maior do que as outras. Consequentemente, calculamos primeiro os valores padronizados, \(z_{jk} = \left(x_{jk} - \bar{x}_k\right)/\sqrt{s_{kk}}\), para que as variáveis contribuam igualmente para a variação no gráfico de dispersão. A maior parte da variação pode ser explicada por uma única variável determinada por uma linha através do conjunto de pontos.
Tabela 1.3 - Dados de Tamanho dos Lagartos
Lagarto | Mass | SVL | HLS |
---|---|---|---|
1 | 5.526 | 59.0 | 113.5 |
2 | 10.401 | 75.0 | 142.0 |
3 | 9.213 | 69.0 | 124.0 |
4 | 8.953 | 67.5 | 125.0 |
5 | 7.063 | 62.0 | 129.5 |
6 | 6.610 | 62.0 | 123.0 |
7 | 11.273 | 74.0 | 140.0 |
8 | 2.447 | 47.0 | 97.0 |
9 | 15.493 | 86.5 | 162.0 |
10 | 9.004 | 69.0 | 126.5 |
11 | 8.199 | 70.5 | 136.0 |
12 | 6.601 | 64.5 | 116.0 |
13 | 7.622 | 67.5 | 135.0 |
14 | 10.067 | 73.0 | 136.5 |
15 | 10.091 | 73.0 | 135.5 |
16 | 10.888 | 77.0 | 139.0 |
17 | 7.610 | 61.5 | 118.0 |
18 | 7.733 | 66.5 | 133.5 |
19 | 12.015 | 79.5 | 150.0 |
20 | 10.049 | 74.0 | 137.0 |
21 | 5.149 | 59.5 | 116.0 |
22 | 9.158 | 68.0 | 123.0 |
23 | 12.132 | 75.0 | 141.0 |
24 | 6.978 | 66.5 | 117.0 |
25 | 6.890 | 63.0 | 117.0 |
Variáveis brutas:
T1.3 <- read.table("JW6Data/T1-3.dat",
quote="\"",
comment.char="")
colnames(T1.3) <- c("Mass",
"SVL",
"HLS")
Lagarto <- 1:25
dados <- cbind(Lagarto, T1.3)
print(summary(dados))
Lagarto Mass SVL HLS
Min. : 1 Min. : 2.447 Min. :47.0 Min. : 97.0
1st Qu.: 7 1st Qu.: 6.978 1st Qu.:63.0 1st Qu.:118.0
Median :13 Median : 8.953 Median :68.0 Median :129.5
Mean :13 Mean : 8.687 Mean :68.4 Mean :129.3
3rd Qu.:19 3rd Qu.:10.091 3rd Qu.:74.0 3rd Qu.:137.0
Max. :25 Max. :15.493 Max. :86.5 Max. :162.0
#knitr::kable(head(dados), caption="Tamanho dos Lagartos")
#knitr::kable(tail(dados), caption="Tamanho dos Lagartos")
GGally::ggpairs(dados[, -1],
title = "Gráficos de Dispersão com Densidade",
diag = list(continuous = "densityDiag"))
p <- plotly::plot_ly(dados,
x = ~Mass,
y = ~SVL,
z = ~HLS,
type = "scatter3d",
mode = "markers",
marker = list(size = 5))
p <- plotly::layout(p,
scene = list(
xaxis = list(title = "Mass"),
yaxis = list(title = "SVL"),
zaxis = list(title = "HLS")))
p
Variáveis padronizadas:
Lagarto Mass SVL HLS
Min. : 1 Min. : 2.447 Min. :47.0 Min. : 97.0
1st Qu.: 7 1st Qu.: 6.978 1st Qu.:63.0 1st Qu.:118.0
Median :13 Median : 8.953 Median :68.0 Median :129.5
Mean :13 Mean : 8.687 Mean :68.4 Mean :129.3
3rd Qu.:19 3rd Qu.:10.091 3rd Qu.:74.0 3rd Qu.:137.0
Max. :25 Max. :15.493 Max. :86.5 Max. :162.0
Mass | SVL | HLS |
---|---|---|
-1.18 | -1.18 | -1.16 |
0.64 | 0.83 | 0.93 |
0.20 | 0.08 | -0.39 |
0.10 | -0.11 | -0.32 |
-0.61 | -0.80 | 0.01 |
-0.77 | -0.80 | -0.46 |
Mass | SVL | HLS | |
---|---|---|---|
20 | 0.51 | 0.70 | 0.56 |
21 | -1.32 | -1.11 | -0.98 |
22 | 0.18 | -0.05 | -0.46 |
23 | 1.29 | 0.83 | 0.86 |
24 | -0.64 | -0.24 | -0.90 |
25 | -0.67 | -0.68 | -0.90 |
GGally::ggpairs(dadoz,
title = "Gráficos de Dispersão com Densidade\nVariáveis padronizadas",
diag = list(continuous = "densityDiag"))
p <- plotly::plot_ly(dadoz,
x = ~Mass,
y = ~SVL,
z = ~HLS,
type = "scatter3d",
mode = "markers",
marker = list(size = 5))
p <- plotly::layout(p,
scene = list(
xaxis = list(title = "zMass"),
yaxis = list(title = "zSVL"),
zaxis = list(title = "zHLS")))
p
Variáveis brutas por sexo:
vetor_sexo <- c("fmffmfmfmfmfm", "mmmfmmmffmff")
Sexo <- unlist(strsplit(vetor_sexo, ""))
dados <- cbind(dados, Sexo)
dados$Sexo <- factor(dados$Sexo)
print(summary(dados))
Lagarto Mass SVL HLS Sexo
Min. : 1 Min. : 2.447 Min. :47.0 Min. : 97.0 f:12
1st Qu.: 7 1st Qu.: 6.978 1st Qu.:63.0 1st Qu.:118.0 m:13
Median :13 Median : 8.953 Median :68.0 Median :129.5
Mean :13 Mean : 8.687 Mean :68.4 Mean :129.3
3rd Qu.:19 3rd Qu.:10.091 3rd Qu.:74.0 3rd Qu.:137.0
Max. :25 Max. :15.493 Max. :86.5 Max. :162.0
Lagarto | Mass | SVL | HLS | Sexo |
---|---|---|---|---|
1 | 5.53 | 59.0 | 113.5 | f |
2 | 10.40 | 75.0 | 142.0 | m |
3 | 9.21 | 69.0 | 124.0 | f |
4 | 8.95 | 67.5 | 125.0 | f |
5 | 7.06 | 62.0 | 129.5 | m |
6 | 6.61 | 62.0 | 123.0 | f |
Lagarto | Mass | SVL | HLS | Sexo | |
---|---|---|---|---|---|
20 | 20 | 10.05 | 74.0 | 137 | m |
21 | 21 | 5.15 | 59.5 | 116 | f |
22 | 22 | 9.16 | 68.0 | 123 | f |
23 | 23 | 12.13 | 75.0 | 141 | m |
24 | 24 | 6.98 | 66.5 | 117 | f |
25 | 25 | 6.89 | 63.0 | 117 | f |
GGally::ggpairs(dados[, -1],
mapping = ggplot2::aes(color = Sexo, alpha=0.4),
title = "Gráficos de Dispersão com Densidade",
diag = list(continuous = "densityDiag"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Quando a estatura de uma criança pequena é medida a cada aniversário, os pontos podem ser plotados e conectados por linhas para produzir um gráfico. Isso é um exemplo de uma curva de crescimento. Em geral, medições repetidas da mesma característica na mesma unidade observacional ou participante podem dar origem a uma curva de crescimento se for esperado um padrão crescente, decrescente ou até mesmo crescente seguido de decrescente.
O Departamento de Pesca e Caça do Alasca monitora urso-cinzento ou
pardo (Ursus arctos horribilis) com o objetivo de
manter uma população saudável. Os ursos são anestesiados com um dardo
para induzir o sono e pesados em uma balança suspensa em um tripé.
Medidas de comprimento são feitas com uma fita métrica de aço. A Tabela
1.4 fornece as massas corporais totais (Wt
) em quilograma e
os comprimentos (Length
) em centímetro de sete ursos fêmeas
nas idades de 2, 3, 4 e 5 anos.
Tabela 1.4 Dados de Ursas
Ursa | Wt2 | Wt3 | Wt4 | Wt5 | Lngth2 | Lngth3 | Lngth4 | Lngth5 |
---|---|---|---|---|---|---|---|---|
1 | 48 | 59 | 95 | 82 | 141 | 157 | 168 | 183 |
2 | 59 | 68 | 102 | 102 | 140 | 168 | 174 | 170 |
3 | 61 | 77 | 93 | 107 | 145 | 162 | 172 | 177 |
4 | 54 | 43 | 104 | 104 | 146 | 159 | 176 | 171 |
5 | 45 | 66 | 84 | 112 | 150 | 158 | 168 | 175 |
6 | 68 | 82 | 95 | 118 | 142 | 140 | 178 | 189 |
7 | 68 | 95 | 109 | 111 | 139 | 171 | 176 | 175 |
Primeiro, para cada ursa, plotamos os pesos em relação às idades e, em seguida, conectamos os pesos nos anos sucessivos por linhas retas. Isso dá uma aproximação da curva de crescimento para a massa corporal total.
No campo, os ursos são pesados em uma balança que indica libras.
Por ser difícil inspecionar visualmente as curvas de crescimento individuais em um gráfico combinado, as curvas individuais devem ser plotadas novamente em um arranjo onde semelhanças e diferenças sejam facilmente observadas. Algumas curvas de crescimento parecem lineares, enquanto outras são não lineares.
A ursa 6 parece ter ficado menor dos 2 aos 3 anos de idade, mas o pesquisador sabe que a medição do comprimento com uma fita de aço pode ser afetada pela postura da ursa quando sedada.
Pessoa (2008) usa a estrutura de dados de série temporal (dados equiespaçados) para plotar a curva de crescimento de cada ursa. Este não é método geral, pois os dados da curva de crescimento podem ser não equiespaçados.
Os intervalos de confiança de Mirnam (2014) são incorretos, pois são intervalos para delineamento entre participantes e sem correção de Bonferroni. Os intervalos de confiança têm que ser adequados ao delineamento intraparticipantes (i.e., com medidas repetidas).
# Growth Curve Analysis and Visualization Using R - Mirman - 2014
Year <- rep(2:5, each = 7)
Bear <- rep(1:7, 4)
Wt2 <- c(48, 59, 61, 54, 45, 68, 68)
Wt3 <- c(59, 68, 77, 43, 66, 82, 95)
Wt4 <- c(95, 102, 93, 104, 84, 95, 109)
Wt5 <- c(82, 102, 107, 104, 112, 118, 111)
Wt <- c(Wt2, Wt3, Wt4, Wt5)
Lngth2 <- c(141, 140, 145, 146, 150, 142, 139)
Lngth3 <- c(157, 168, 162, 159, 158, 140, 171)
Lngth4 <- c(168, 174, 172, 176, 168, 178, 176)
Lngth5 <- c(183, 170, 177, 171, 175, 189, 175)
Lngth <- c(Lngth2, Lngth3, Lngth4, Lngth5)
dados_ursas <- data.frame(Year, Bear, Wt, Lngth)
dados_ursas$Bear <- factor(dados_ursas$Bear)
print(str(dados_ursas))
'data.frame': 28 obs. of 4 variables:
$ Year : int 2 2 2 2 2 2 2 3 3 3 ...
$ Bear : Factor w/ 7 levels "1","2","3","4",..: 1 2 3 4 5 6 7 1 2 3 ...
$ Wt : num 48 59 61 54 45 68 68 59 68 77 ...
$ Lngth: num 141 140 145 146 150 142 139 157 168 162 ...
NULL
Year | Bear | Wt | Lngth |
---|---|---|---|
2 | 1 | 48 | 141 |
2 | 2 | 59 | 140 |
2 | 3 | 61 | 145 |
2 | 4 | 54 | 146 |
2 | 5 | 45 | 150 |
2 | 6 | 68 | 142 |
2 | 7 | 68 | 139 |
3 | 1 | 59 | 157 |
3 | 2 | 68 | 168 |
3 | 3 | 77 | 162 |
3 | 4 | 43 | 159 |
3 | 5 | 66 | 158 |
3 | 6 | 82 | 140 |
3 | 7 | 95 | 171 |
4 | 1 | 95 | 168 |
4 | 2 | 102 | 174 |
4 | 3 | 93 | 172 |
4 | 4 | 104 | 176 |
4 | 5 | 84 | 168 |
4 | 6 | 95 | 178 |
4 | 7 | 109 | 176 |
5 | 1 | 82 | 183 |
5 | 2 | 102 | 170 |
5 | 3 | 107 | 177 |
5 | 4 | 104 | 171 |
5 | 5 | 112 | 175 |
5 | 6 | 118 | 189 |
5 | 7 | 111 | 175 |
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Wt, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10) +
ggplot2::facet_wrap(~ Bear)
ggplot2::ggsave("BearWtSep.png", width=3, height=3, dpi=300)
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Lngth, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10) +
ggplot2::facet_wrap(~ Bear)
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Wt, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10)
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Lngth, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10)
alfa <- 0.05
ic <- summarySEwithin2(dados_ursas,
measurevar="Wt",
withinvars="Year",
idvar="Bear",
na.rm=TRUE,
conf.interval=
1-alfa/length(unique(dados_ursas$Year)))
Automatically converting the following non-factors to factors: Year
Anexando pacote: 'data.table'
Os seguintes objetos são mascarados por 'package:reshape2':
dcast, melt
O seguinte objeto é mascarado por 'package:HH':
transpose
Os seguintes objetos são mascarados por 'package:dplyr':
between, first, last
O seguinte objeto é mascarado por 'package:DescTools':
%like%
Year N Wt WtNormed sd se ci
1 2 7 57.6 57.6 4.04 1.53 5.38
2 3 7 70.0 70.0 11.70 4.42 15.57
3 4 7 97.4 97.4 10.35 3.91 13.78
4 5 7 105.1 105.1 9.52 3.60 12.67
grf <- ggplot2::ggplot(ic,
ggplot2::aes(x=Year,
y=Wt)) +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=Wt-ci,
ymax=Wt+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylab("Weight") +
ggplot2::ggtitle("7 Female Bear: Weight\nWithin-subject CI95 Bonferroni") +
ggplot2::theme_bw()
print(grf, digits=3)
ic <- summarySEwithin2(dados_ursas,
measurevar="Lngth",
withinvars="Year",
idvar="Bear",
na.rm=TRUE,
conf.interval=
1-alfa/length(unique(dados_ursas$Year)))
Automatically converting the following non-factors to factors: Year
Year N Lngth LngthNormed sd se ci
1 2 7 143.2857 143.2857 4.993647 1.887421 6.64603
2 3 7 159.2857 159.2857 10.700096 4.044256 14.24073
3 4 7 173.1429 173.1429 4.427547 1.673456 5.89261
4 5 7 177.1429 177.1429 8.369446 3.163353 11.13887
grf <- ggplot2::ggplot(ic,
ggplot2::aes(x=Year,
y=Lngth)) +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=Lngth-ci,
ymax=Lngth+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylab("Length") +
ggplot2::ggtitle("7 Female Bear: Length\nWithin-subject CI95 Bonferroni") +
ggplot2::theme_bw()
print(grf)
Um gráfico de radar, também conhecido como gráfico de teia ou gráfico de estrela, é um tipo de representação gráfica que exibe dados em um formato circular e poligonal. Ele é usado para comparar múltiplas variáveis ou dimensões em um único gráfico, permitindo visualizar padrões, tendências e diferenças entre diferentes conjuntos de dados.
Suponha que cada unidade de dados consista em observações não negativas em \(p\) variáveis. Em duas dimensões, podemos construir círculos de um raio fixo (referência) com \(p\) raios igualmente espaçados emanando do centro do círculo. Os comprimentos dos raios representam os valores das variáveis. As extremidades dos raios podem ser conectadas com linhas retas para formar uma estrela. Cada estrela representa uma observação multivariada, e as estrelas podem ser agrupadas de acordo com suas similaridades (subjetivas).
Frequentemente, ao construir as estrelas, é útil padronizar as observações. Nesse caso, algumas das observações podem ser negativas. As observações podem então ser reexpressas de forma que o centro do círculo represente a menor observação padronizada em todo o conjunto de dados.
O gráfico de radar é útil para destacar diferenças e semelhanças entre múltiplas variáveis. No entanto, é importante lembrar que, se você tiver muitas variáveis, o gráfico pode se tornar confuso e difícil de interpretar. Nesses casos, pode ser necessário considerar outras formas de visualização de dados, como gráficos de barras ou gráficos de dispersão.
As estrelas representando as primeiras 5 das 22 empresas de serviços públicos na Tabela 12.4 são mostradas na Figura 1.16. Existem oito variáveis; consequentemente, as estrelas são octógonos distorcidos.
As observações em todas as variáveis foram padronizadas. Entre as primeiras cinco empresas de serviços públicos, a menor observação padronizada para qualquer variável foi -1.6. Tomando esse valor como zero, as variáveis são plotadas em escalas idênticas ao longo de oito raios equidistantes que se originam no centro do círculo. As variáveis são ordenadas no sentido horário, começando na posição das 12 horas.
À primeira vista, nenhuma dessas empresas de serviços públicos parece ser similar a qualquer outra. No entanto, devido à maneira como as estrelas são construídas, cada variável recebe peso igual na impressão visual. Se nos concentrarmos nas variáveis 6 (vendas em quilowatts-hora [kWh] por ano) e 8 (custos totais de combustível em centavos por kWh), então Boston Edison e Consolidated Edison são similares (pequena variável 6, grande variável 8), e Arizona Public Service, Central Louisiana Electric e Commonwealth Edison são similares (variável 6 moderada, variável 8 moderada).
Figura 1.16: Estrelas para as primeiras cinco empresas de serviços públicos.
# X1: Razão de cobertura de encargos fixos (renda/dívida).
# X2: Taxa de retorno sobre o capital.
# X3: Custo por capacidade KW instalada.
# X4: Fator de carga anual.
# Y5: Crescimento do pico de demanda em kWh de 1974 para 1975.
# X6: Vendas (uso de kWh por ano).
# X7: Percentual nuclear.
# X8: Custos totais de combustível (centavos por kWh).
# Empresa
dados <- read.table("JW6Data/T12-4.DAT", quote="\"", comment.char="")
colnames(dados) <- c("X1", "X2", "X3", "X4",
"X5", "X6", "X7", "X8",
"Empresa")
dados$Empresa <- factor(dados$Empresa)
print(str(dados))
'data.frame': 22 obs. of 9 variables:
$ X1 : num 1.06 0.89 1.43 1.02 1.49 1.32 1.22 1.1 1.34 1.12 ...
$ X2 : num 9.2 10.3 15.4 11.2 8.8 13.5 12.2 9.2 13 12.4 ...
$ X3 : int 151 202 113 168 192 111 175 245 168 197 ...
$ X4 : num 54.4 57.9 53 56 51.2 60 67.6 57 60.4 53 ...
$ X5 : num 1.6 2.2 3.4 0.3 1 -2.2 2.2 3.3 7.2 2.7 ...
$ X6 : int 9077 5088 9212 6423 3300 11127 7642 13082 8406 6455 ...
$ X7 : num 0 25.3 0 34.3 15.6 22.5 0 0 0 39.2 ...
$ X8 : num 0.628 1.555 1.058 0.7 2.044 ...
$ Empresa: Factor w/ 22 levels "Arizona","Boston",..: 1 2 3 4 5 6 7 8 9 10 ...
NULL
X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 | Empresa |
---|---|---|---|---|---|---|---|---|
1.06 | 9.2 | 151 | 54.4 | 1.6 | 9077 | 0.0 | 0.628 | Arizona |
0.89 | 10.3 | 202 | 57.9 | 2.2 | 5088 | 25.3 | 1.555 | Boston |
1.43 | 15.4 | 113 | 53.0 | 3.4 | 9212 | 0.0 | 1.058 | Central |
1.02 | 11.2 | 168 | 56.0 | 0.3 | 6423 | 34.3 | 0.700 | Common |
1.49 | 8.8 | 192 | 51.2 | 1.0 | 3300 | 15.6 | 2.044 | Consolid |
1.32 | 13.5 | 111 | 60.0 | -2.2 | 11127 | 22.5 | 1.241 | Florida |
1.22 | 12.2 | 175 | 67.6 | 2.2 | 7642 | 0.0 | 1.652 | Hawaiian |
1.10 | 9.2 | 245 | 57.0 | 3.3 | 13082 | 0.0 | 0.309 | Idaho |
1.34 | 13.0 | 168 | 60.4 | 7.2 | 8406 | 0.0 | 0.862 | Kentucky |
1.12 | 12.4 | 197 | 53.0 | 2.7 | 6455 | 39.2 | 0.623 | Madison |
0.75 | 7.5 | 173 | 51.5 | 6.5 | 17441 | 0.0 | 0.768 | Nevada |
1.13 | 10.9 | 178 | 62.0 | 3.7 | 6154 | 0.0 | 1.897 | NewEngla |
1.15 | 12.7 | 199 | 53.7 | 6.4 | 7179 | 50.2 | 0.527 | Northern |
1.09 | 12.0 | 96 | 49.8 | 1.4 | 9673 | 0.0 | 0.588 | Oklahoma |
0.96 | 7.6 | 164 | 62.2 | -0.1 | 6468 | 0.9 | 1.400 | Pacific |
1.16 | 9.9 | 252 | 56.0 | 9.2 | 15991 | 0.0 | 0.620 | Puget |
0.76 | 6.4 | 136 | 61.9 | 9.0 | 5714 | 8.3 | 1.920 | SanDiego |
1.05 | 12.6 | 150 | 56.7 | 2.7 | 10140 | 0.0 | 1.108 | Southern |
1.16 | 11.7 | 104 | 54.0 | -2.1 | 13507 | 0.0 | 0.636 | Texas |
1.20 | 11.8 | 148 | 59.9 | 3.5 | 7287 | 41.1 | 0.702 | Wisconsi |
1.04 | 8.6 | 204 | 61.0 | 3.5 | 6650 | 0.0 | 2.116 | United |
1.07 | 9.3 | 174 | 54.3 | 5.9 | 10093 | 26.6 | 1.306 | Virginia |
stars(dados[, 1:8], len = 0.8, key.loc = c(12, 2),
labels = abbreviate(case.names(dados[, 1:8])),
main = "", full = FALSE)
stars(dados[, 1:8], labels = abbreviate(case.names(dados[, 1:8])),
key.loc = c(13, 1.5), main = "", len = 0.8)
stars(dados[, 1:8], len = 0.8, key.loc = c(12, 2),
labels = abbreviate(case.names(dados[, 1:8])),
main = "", draw.segments = TRUE)
# Função para calcular máximo e mínimo de uma coluna
calcular_max_min <- function(x) {
max_val <- max(x)
min_val <- min(x)
return(c(Maximo = max_val, Minimo = min_val))
}
# Aplicar a função calcular_max_min em todas as colunas, exceto a última (Empresa)
resultado <- sapply(dados[, -ncol(dados)], calcular_max_min)
# Converter o resultado em um dataframe
maxmin <- data.frame(t(resultado))
knitr::kable(maxmin, caption="Dados da Empresa")
Maximo | Minimo | |
---|---|---|
X1 | 1.490 | 0.750 |
X2 | 15.400 | 6.400 |
X3 | 252.000 | 96.000 |
X4 | 67.600 | 49.800 |
X5 | 9.200 | -2.200 |
X6 | 17441.000 | 3300.000 |
X7 | 50.200 | 0.000 |
X8 | 2.116 | 0.309 |
dat <- rbind(t(maxmin),dados[,-ncol(dados)])
par(mfrow = c(2,2))
for(i in 3:6) {
fmsb::radarchart(dat[c(1:2, i), ],
title=as.vector(dados[,"Empresa"])[i-2],
pcol = "blue", plwd = 2, plty = 1,
axislabcol = "black", axistype = "radial")
}
As pessoas reagem a faces humanas.
Chernoff, em seu artigo de 1973, propôs o uso de “Chernoff faces” como um método para representar visualmente dados multivariados em um formato bidimensional. Cada observação, representada por um conjunto de variáveis em \(p\) dimensões, é transformada em um “rosto” ao mapear essas variáveis em características faciais específicas. As características faciais, como formato do rosto, curvatura da boca, tamanho dos olhos, etc., são determinadas pelos valores das variáveis.
O método é especialmente útil ao lidar com conjuntos de dados com até 18 variáveis. A atribuição de variáveis às características faciais é subjetiva e geralmente feita pelo pesquisador. Diferentes escolhas na atribuição de variáveis às características podem levar a representações diferentes, portanto, pode ser necessário fazer tentativas até alcançar visualizações satisfatórias.
As Chernoff faces podem servir a diversos propósitos na análise de dados:
Verificação de Grupos: Elas são úteis para verificar agrupamentos iniciais sugeridos pelo conhecimento do assunto ou intuição.
Análise de Clusters: Podem ser usadas para visualizar agrupamentos finais produzidos por algoritmos de clusterização.
Este exemplo ilustra o uso das faces de Chernoff para representar dados de empresas de serviços públicos em 22 empresas públicas. As variáveis são mapeadas em características faciais e rostos similares são agrupados. Ao combinar clusters, é possível obter um número menor de clusters para uma representação mais concisa.
Variável | Característica Facial |
---|---|
X1: Cobertura de Encargos Fixos | Meia altura do rosto |
X2: Taxa de Retorno sobre o Capital | Largura do rosto |
X3: Custo por kW de Capacidade Instalada | Posição do centro da boca |
X4: Fator de Carga Anual | Inclinação dos olhos |
X5: Crescimento de Demanda de kWh Pico desde 1974 | Excentricidade (altura/comprimento) dos olhos |
X6: Vendas (kWh utilizados por ano) | Meio comprimento do olho |
X7: Percentual Nuclear | Curvatura da boca |
X8: Custo Total do Combustível (centavos por kWh) | Comprimento do nariz |
Nesta tabela, as variáveis estão traduzidas para o português, assim como as características faciais utilizadas para criar os rostos de Chernoff para representar os dados multivariados.
Figura 1.17 Faces de Chernoff para 22 empresas de serviços públicos.
Este exemplo demonstra outro uso das Chernoff faces para acompanhar o bem-estar financeiro de uma empresa ao longo do tempo. Cada característica facial representa um único indicador financeiro, facilitando a observação das mudanças longitudinais nos indicadores.
Figura 1.18 Faces de Chernoff ao longo do tempo.
Além disso, as faces de Chernoff podem ser usadas para exibir diferenças em observações multivariadas em duas dimensões. Por exemplo, as localizações geográficas (latitude e longitude) podem ser representadas pelos eixos de coordenadas, enquanto os rostos representam medições de várias cidades.
As faces de Chernoff devem ser construídas com a ajuda de um computador, pois a padronização dos dados e a determinação das localizações, tamanhos e orientações das características faciais fazem parte do processo. Com algum treinamento, esses rostos podem ser usados para comunicar efetivamente similaridades e diferenças nos dados.
Em resumo, as faces de Chernoff oferecem uma forma criativa de visualizar dados multivariados e podem ser benéficas em várias tarefas de análise de dados, especialmente quando se trata de conjuntos de dados relativamente pequenos. À medida que a tecnologia de gráficos de computador avança, o potencial para melhorias e aplicações adicionais das faces de Chernoff aumenta.
# Chernoff, H. (1973) The use of faces to represent statistiscal association,
# JASA, 68, pp 361–368.
means <- lapply(iris[,-5], tapply, iris$Species, mean)
means
$Sepal.Length
setosa versicolor virginica
5.006 5.936 6.588
$Sepal.Width
setosa versicolor virginica
3.428 2.770 2.974
$Petal.Length
setosa versicolor virginica
1.462 4.260 5.552
$Petal.Width
setosa versicolor virginica
0.246 1.326 2.026
m <- t(do.call(rbind, means))
m <- cbind(m, matrix(rep(1, 11*3), nrow=3))
# define the colors, first for all faces the same
col <- replicate(3, c("orchid1", "olivedrab", "goldenrod4",
"peachpuff", "darksalmon", "peachpuff3"))
rownames(col) <- c("nose","eyes","hair","face","lips","ears")
# change haircolor individually for each face
col[3, ] <- c("lightgoldenrod", "coral3", "sienna4")
z <- DescTools::PlotFaces(m, nr=1, nc=3, col=col)
modified.item variable
1 height of face Sepal.Length
2 width of face Sepal.Width
3 structure of face Petal.Length
4 height of mouth Petal.Width
5 width of mouth
6 smiling
7 height of eyes
8 width of eyes
9 height of hair
10 width of hair
11 style of hair
12 height of nose
13 width of nose
14 width of ear
15 height of ear
O arquivo datasets::longley
contém o conjunto de 7
variáveis econômicas, observadas anualmente de 1947 a 1962 (\(n = 16\)).
GNP.deflator
: Deflator de preço implícito do Produto
Nacional Bruto (PNB) (1954 = 100).
GNP
: Produto Nacional Bruto.
Unemployed
: Número de desempregados.
Armed.Forces
: Número de pessoas nas forças
armadas.
Population
: População ‘não institucionalizada’ com
idade ≥ 14 anos.
Year
: Ano (tempo).
Employed
: Número de pessoas empregadas.
modified.item variable
1 height of face GNP.deflator
2 width of face GNP
3 structure of face Unemployed
4 height of mouth Armed.Forces
5 width of mouth Population
6 smiling Year
7 height of eyes Employed
8 width of eyes GNP.deflator
9 height of hair GNP
10 width of hair Unemployed
11 style of hair Armed.Forces
12 height of nose Population
13 width of nose Year
14 width of ear Employed
15 height of ear GNP.deflator
effect of variables:
modified item Var
"height of face " "GNP.deflator"
"width of face " "GNP"
"structure of face" "Unemployed"
"height of mouth " "Armed.Forces"
"width of mouth " "Population"
"smiling " "Year"
"height of eyes " "Employed"
"width of eyes " "GNP.deflator"
"height of hair " "GNP"
"width of hair " "Unemployed"
"style of hair " "Armed.Forces"
"height of nose " "Population"
"width of nose " "Year"
"width of ear " "Employed"
"height of ear " "GNP.deflator"
aplpack::plot.faces(a,
longley[1:16,2],
longley[1:16,3],
width=35,
height=30,
print.info=TRUE,
face.type=0,
cex=1,
na.rm=TRUE)
Embora possam parecer formidáveis à primeira vista, a maioria das técnicas multivariadas são baseadas no conceito simples de distância.
A distância em linha reta, ou distância euclidiana, deve ser familiar. Se considerarmos o ponto \(P = (x_1, x_2)\) no plano, a distância em linha reta de \(P\) à origem \(O = (0, 0)\) é, de acordo com o teorema de Pitágoras:
\[ d(O, P) = \sqrt{x_1^2 + x_2^2} \tag{1-9} \]
A situação é ilustrada na Figura 1.19. Em geral, se o ponto \(P\) tem \(p\) coordenadas, de modo que \(P = (x_1, x_2, \ldots, x_p)\), a distância em linha reta de \(P\) à origem \(O = (0, 0, \ldots, 0)\) é:
\[ d(O, P) = \sqrt{x_1^2 + x_2^2 + \cdots + x_p^2} \tag{1-10} \]
Todos os pontos \((x_1, x_2, \ldots, x_p)\) que estão a uma distância constante ao quadrado, como \(c^2\), da origem satisfazem a equação:
\[ d^2(O, P) = x_1^2 + x_2^2 + \cdots + x_p^2 = c^2 \tag{1-11} \]
Porque esta é a equação de uma hiperesfera (um círculo se \(p = 2\)), os pontos equidistantes da origem estão em uma hiperesfera.
A distância em linha reta entre dois pontos arbitrários \(P\) e \(Q\) com coordenadas \(P = (x_1, x_2, \ldots, x_p)\) e \(Q = (y_1, y_2, \ldots, y_p)\) é dada por:
\[ d(P, Q) = \sqrt{(x_1 - y_1)^2 + (x_2 - y_2)^2 + \cdots + (x_p - y_p)^2} \tag{1-12} \]
A distância em linha reta, ou distância euclidiana, é insatisfatória para a maioria dos propósitos estatísticos. Isso ocorre porque cada coordenada contribui igualmente para o cálculo da distância euclidiana. Quando as coordenadas representam medições sujeitas a flutuações aleatórias de magnitudes diferentes, muitas vezes é desejável ponderar menos intensamente as coordenadas sujeitas a uma grande variabilidade do que aquelas que não são altamente variáveis. Isso sugere uma medida diferente de distância.
Figura 1.19 - Distância dada pelo teorema de Pitágoras.
Nosso objetivo agora é desenvolver uma “distância estatística” que leve em conta as diferenças nas variações e, em devido tempo, a presença de correlação. Como nossa escolha dependerá das variâncias e covariâncias das amostras, neste ponto, usamos o termo distância estatística para distingui-la da distância euclidiana comum. É a distância estatística que é fundamental para a análise multivariada.
Para começar, tomamos como fixo o conjunto de observações plotadas como o gráfico de dispersão em \(p\) dimensões. A partir disso, construiremos uma medida de distância da origem a um ponto \(P = (x_1, x_2, \ldots, x_p)\). Em nossos argumentos, as coordenadas \((x_1, x_2, \ldots, x_p)\) de \(P\) podem variar para produzir diferentes localizações para o ponto. No entanto, os dados que determinam a distância permanecerão fixos.
Para ilustrar, suponha que temos \(n\) pares de medições em duas variáveis, cada uma com média zero. Chame as variáveis de \(x_1\) e \(x_2\) e assuma que as medições de \(x_1\) variam independentemente das medições de \(x_2\). Além disso, suponha que a variabilidade nas medições de \(x_1\) seja maior do que a variabilidade nas medições de \(x_2\). Um gráfico de dispersão dos dados se pareceria com o que está representado na Figura 1.20.
Figura 1.20 Um gráfico de dispersão com maior variabilidade na direção de x1 do que na direção de x2.
Ao olharmos para a Figura 1.20, podemos perceber que valores que estão a uma dada distância da origem na direção de \(x_1\) não são tão “surpreendentes” ou “inusitados” quanto valores equidistantes da origem na direção de \(x_2\). Isso ocorre porque a variabilidade inerente na direção de \(x_1\) é maior do que a variabilidade na direção de \(x_2\). Consequentemente, coordenadas grandes de \(x_1\) (em valor absoluto) não são tão inesperadas quanto coordenadas grandes de \(x_2\). Parece razoável, então, ponderar uma coordenada de \(x_2\) mais fortemente do que uma coordenada de \(x_1\) de mesmo valor ao calcular a “distância” à origem.
Uma maneira de proceder é dividir cada coordenada pelo desvio padrão da amostra. Portanto, após a divisão pelos desvios padrão, temos as coordenadas “padronizadas” \(x^*_1 = \dfrac{x_1}{\sqrt{s_{11}}}\) e \(x_2^* = \dfrac{x_2}{\sqrt{s_{22}}}\). As coordenadas padronizadas estão agora em pé de igualdade entre si. Após levarmos em conta as diferenças na variabilidade, determinamos a distância usando a fórmula euclidiana padrão.
Assim, a distância estatística do ponto \(P = (x_1, x_2)\) da origem \(O = (0, 0)\) pode ser calculada a partir de suas coordenadas padronizadas \(x^* = \dfrac{x_1}{\sqrt{s_{11}}}\) e \(x_2^* = \dfrac{x_2}{\sqrt{s_{22}}}\) como
\[ \begin{align} d(O, P) &= \sqrt{(x_1^*)^2 + (x_2^*)^2} \\ d(O, P) &= \sqrt{\left(\dfrac{x_1}{\sqrt{s_{11}}}\right)^2+ \left(\dfrac{x_2}{\sqrt{s_{22}}}\right)^2}\\ d(O, P) &= \sqrt{\dfrac{x_1^2}{s_{11}}+ \dfrac{x_2^2}{s_{22}}}\end{align} \tag{1-13} \]
Comparando (1-13) com (1-9), percebemos que a diferença entre as duas expressões se deve aos pesos \(k_1 = \dfrac{1}{\sqrt{s_{11}}}\) e \(k_2 = \dfrac{1}{\sqrt{s_{22}}}\) atribuídos a \(x_1\) e \(x_2\) em (1-13). Observe que se as variâncias das amostras são iguais, \(k_1 \approx k_2\), então \(x_1\) e \(x_2\) receberão o mesmo peso. Em casos em que os pesos são iguais, é conveniente ignorar o divisor comum e usar a fórmula usual da distância euclidiana. Em outras palavras, se a variabilidade na direção de \(x_1\) é a mesma que a variabilidade na direção de \(x_2\) e os valores de \(x_1\) variam independentemente dos valores de \(x_2\), a distância euclidiana é apropriada.
Usando (1-13), vemos que todos os pontos que têm coordenadas \((x_1, x_2)\) e estão a uma distância constante ao quadrado \(c^2\) da origem devem satisfazer a equação:
\[ \dfrac{x_1^2}{s_{11}} + \dfrac{x_2^2}{s_{22}} = c^2 \tag{1-14} \]
A Equação (1-14) é a equação de uma elipse centrada na origem, cujos eixos maior e menor coincidem com os eixos de coordenadas. Ou seja, a distância estatística em (1-13) tem uma elipse como a locação de todos os pontos a uma distância constante da origem. Esse caso geral é mostrado na Figura 1.21.
Figura 1.21 - A elipse de distância estatística constant c.
Um conjunto de medidas emparelhadas \((x_1, x_2)\) em duas variáveis resulta em \(\bar{x}_1 = \bar{x}_2 = 0\), \(s_{11} = 4\) e \(s_{22} = 1\). Suponha que as medidas de \(x_1\) não estejam relacionadas às medidas de \(x_2\); ou seja, as medidas dentro de um par variam independentemente uma da outra. Como as variâncias da amostra são desiguais, medimos o quadrado da distância de um ponto arbitrário \(P = (x_1, x_2)\) à origem \(O = (0, 0)\) por:
\[ d^2(O, P) = \dfrac{x_1^2}{4} + \dfrac{x_2^2}{1} \]
Todos os pontos \((x_1, x_2)\) que estão a uma distância constante 1 da origem satisfazem a equação:
\[ \dfrac{x_1^2}{4} + \dfrac{x_2^2}{1} = 1 \]
As coordenadas de alguns pontos a uma distância unitária da origem são apresentadas na seguinte tabela:
Coordenadas \((x_1 , x_2)\) | Distância \((x_1^2/4 + x_2^2/1 = 1)\) |
---|---|
\((0,1)\) | \(0^2/4 + 1^2/1 = 1\) |
\((0,-1)\) | \(0^2/4 + (-1)^2/1 = 1\) |
\((2,0)\) | \(2^2/4 + 0^2/1 = 1\) |
\((1, \sqrt{3}/2)\) | \(1^2/4 + (\sqrt{3}/2)^2/1 = 1\) |
Figura 1.22 - Elipse de distância unitária, x12/4 + x22/1 = 1.
# Definir os parâmetros da elipse
center_x <- 0
center_y <- 0
a <- 2 # semi-eixo maior
b <- 1 # semi-eixo menor
# Criar uma sequência de ângulos
theta <- seq(0, 2*pi, length.out = 100)
# Calcular as coordenadas x e y da elipse
x1 <- center_x + a * cos(theta)
x2 <- center_y + b * sin(theta)
# Plotar a elipse
plot(x1, x2,
type = "l", asp = 1,
main="Elipse",
xlab = "x1", ylab = "x2",
xlim=c(-2,2), ylim=c(-2,2))
abline(h=0,v=0, lty=2)
points(center_x, center_y)
A expressão em (1-13) pode ser generalizada para acomodar o cálculo da distância estatística de um ponto arbitrário \(P = (x_1, x_2)\) para qualquer ponto fixo \(Q = (y_1, y_2)\).
Se assumirmos que as variáveis de coordenadas variam independentemente uma da outra, a distância de \(P\) para \(Q\) é dada por
\[ d(P,Q) = \sqrt{\dfrac{(x_1 - y_1)^2}{s_{11}} + \dfrac{(x_2 - y_2)^2}{s_{22}}} \tag{1-15} \]
A extensão dessa distância estatística para mais de duas dimensões é direta. Deixe os pontos \(P\) e \(Q\) terem \(p\) coordenadas tais que \(P = (x_1, x_2, \ldots, x_p)\) e \(Q = (y_1 , y_2, \ldots, y_p)\).
Suponha que \(Q\) seja um ponto fixo (pode ser a origem \(O = (0, 0,... , 0)\)) e as variáveis de coordenadas variem independentemente umas das outras.
Sejam \(s_{11}, s_{22}, \ldots, s_{pp}\) serem variâncias amostrais construídas a partir de \(n\) medições em \(x_1, x_2, \ldots, x_p\), respectivamente. Então a distância estatística de \(P\) para \(Q\) é
\[ d(P,Q) = \sqrt{\dfrac{(x_1 - y_1)^2}{s_{11}} + \dfrac{(x_2 - y_2)^2}{s_{22}}+ \cdots + \dfrac{(x_p - y_p)^2}{s_{pp}}} \tag{1-16} \]
Todos os pontos \(P\) que estão a uma distância quadrada constante de \(Q\) estão em um hiperelipsoide centrado em \(Q\) cujos eixos maior e menor são paralelos aos eixos de coordenadas.
Notamos o seguinte:
A distância de \(P\) à origem \(O\) é obtida ao definir \(y_1 = y_2 = \cdots = y_p = 0\) em (1-16).
Se \(s_{11} = s_{22} = \cdots = s_{pp}\), a fórmula da distância euclidiana em (1-12) é apropriada.
A distância em (1-16) ainda não inclui a maioria dos casos importantes que encontraremos, por causa da suposição de coordenadas independentes. O gráfico de dispersão na Figura 1-23 mostra uma situação bidimensional na qual as medições \(x_1\) não variam independentemente das medições \(x_2\). Na verdade, as coordenadas dos pares \((x_1, x_2)\) tendem a ser grandes ou pequenas juntas, e o coeficiente de correlação de amostra é positivo. Além disso, a variabilidade na direção \(x_2\) é maior do que a variabilidade na direção \(x_1\).
Qual é uma medida significativa de distância quando a variabilidade na direção \(x_1\) é diferente da variabilidade na direção \(x_2\) e as variáveis \(x_1\) e \(x_2\) estão correlacionadas? Na verdade, podemos usar o que já introduzimos, desde que olhemos para as coisas da maneira correta. Da Figura 1.23, vemos que se rotacionarmos o sistema de coordenadas original pelo ângulo \(\theta\) mantendo a dispersão fixa e rotulando os eixos rotacionados \(x_1^{\prime}\) e \(x_2^{\prime}\), a dispersão em termos dos novos eixos se parece muito com a da Figura 1.20. (Você pode querer girar o livro para colocar os eixos \(x_1\) e \(x_2\) em suas posições habituais.) Isso sugere que calculemos as variâncias amostrais usando as coordenadas \(\tilde{x}_1\) e \(\tilde{x}_2\) e medindo a distância como na Equação (1-13). Ou seja, em referência aos eixos \(\tilde{x}_1\) e \(\tilde{x}_2\), definimos a distância do ponto \(P (\tilde{x}_1, \tilde{x}_2)\) à origem \(O = (0, 0)\):
\[ d(P,Q) = \sqrt{\dfrac{\tilde{x}_1^2}{\tilde{s}_{11}} + \dfrac{\tilde{x}_2^2}{\tilde{s}_{22}}} \tag{1-17} \]
em que \(\tilde{s}_{11}\) e \(\tilde{s}_{22}\) denotam as variâncias amostrais calculadas com as medições \(\tilde{x}_1\) e \(\tilde{x}_2\).
Figura 1.23 - Um gráfico de dispersão para medições positivamente correlacionadas e um sistema de coordenadas rotacionado.
A relação entre as coordenadas originais \(x_1\) e \(x_2\) e as coordenadas rotacionadas \(\tilde{x}_1\) e \(\tilde{x}_2\) é fornecida por:
\[ \begin{align} \tilde{x}_1 &= x_1 \cos(\theta) + x_2 \sin(\theta) \\ \tilde{x}_2 &= -x_1 \sin(\theta) + x_2 \cos(\theta) \end{align} \]
Dadas as relações acima, podemos formalmente substituir \(x_1\) e \(x_2\) em uma expressão anterior e expressar a distância em termos das coordenadas originais.
Após algumas manipulações algébricas simples, a distância de \(P = (x_1, x_2)\) à origem \(O = (0, 0)\) pode ser escrita em termos das coordenadas originais \(x_1\) e \(x_2\) de \(P\) como:
\[ d(O, P) = \sqrt{a_{11} x_1 + 2 a_{12} x_1 x_2 + a_{22} x_2} \tag{1-19} \]
em que os \(a\) são números tais que a distância é não negativa para todos os possíveis valores de \(x_1\) e \(x_2\). Aqui, \(a_{11}\), \(a_{12}\) e \(a_{22}\) são determinados pelo ângulo \(\theta\), e as covariâncias \(s_n\), \(s_{12}\) e \(s_{22}\) calculadas a partir dos dados originais.
A equação (1-19) pode ser comparada com outra expressão. A expressão anterior pode ser considerada um caso especial da expressão (1-19) com \(a_{11} = 1/\sigma_1^2\), \(a_{22} = 1/\sigma_2^2\) e \(a_{12} = 0\).
Em geral, a distância estatística do ponto \(P = (x_1, x_2)\) do ponto fixo \(Q = (y_1, y_2)\) em situações em que as variáveis estão correlacionadas tem a forma geral:
\[ d(P, Q) = \sqrt{a_{11} (x_1 - y_1)^2 + 2 a_{12} (x_1 - y_1)(x_2 - y_2) + a_{22} (x_2 - y_2)^2} \tag{1-20} \]
e sempre pode ser calculada uma vez que \(a_{11}\), \(a_{12}\) e \(a_{22}\) são conhecidos. Além disso, as coordenadas de todos os pontos \(P = (x_1, x_2)\) que estão a uma distância constante \(c^2\) de \(Q\) satisfazem:
\[ a_{11} (x_1 - y_1)^2 + 2 a_{12} (x_1 - y_1)(x_2 - y_2) + a_{22} (x_2 - y_2)^2 = c^2 \tag{1-21} \]
Por definição, esta é a equação de uma elipse centrada em \(Q\).
O gráfico de tal equação é exibido na Figura 1.24. Os eixos principais (longos) e secundários (curtos) são indicados. Eles são paralelos aos eixos \(X_1\) e \(x_2\). Para a escolha de \(B\), \(a_{12}\) e \(a_{22}\) na nota de rodapé 2, os eixos \(X_1\) e \(x_2\) estão em um ângulo \(\theta\) em relação aos eixos \(X_1\) e \(x_2\).
A generalização das fórmulas de distância de (1-19) e (1-20) para \(p\) dimensões é direta. Seja \(P = (x_1, x_2, \ldots, x_p)\) um ponto cujas coordenadas representam variáveis que estão correlacionadas e sujeitas a variabilidade inerente. Seja \(O = (0, 0, \ldots, 0)\) o ponto de origem e \(Q = (y_1, y_2, \ldots, y_p)\) um ponto fixo especificado. Então, as distâncias de \(P\) para \(O\) e de \(P\) para \(Q\) têm as formas gerais:
\[ d(O, P) = \sqrt{a_{11}x_1 + a_{22}x_2 + \cdots + a_{pp}x_p +\\ 2a_{12}x_1x_2 + \\ 2a_{13}x_1x_3 + \cdots + 2a_{p-1,p}x_{p-1}x_p} \tag{1-22} \]
e
\[ d(P, Q) = \sqrt{a_{11}(x_1 - y_1)^2 + a_{22}(x_2 - y_2)^2 + \cdots + a_{pp}(x_p - y_p)^2 + \\ 2a_{12}(x_1 - y_1)(x_2 - y_2) + \\ 2a_{13}(x_1 - y_1)(x_3 - y_3) + \cdots + \\ 2a_{p-1,p}(x_{p-1} - y_{p-1})(x_p - y_p)} \tag{1-23} \]
sendo que os \(a\) são números tais que as distâncias são sempre não negativas.
As expressões algébricas para os quadrados das distâncias em (1-22) e (1-23) são conhecidas como formas quadráticas e, em particular, formas quadráticas definidas positivas. É possível exibir essas formas quadráticas de maneira mais simples usando álgebra de matrizes.
Figura 1.24 - Elipse de pontos a uma distância constante do ponto Q.
A seguir está o gráfico da elipse com centro em (1,2), rotacionada em 30 graus e 50 pontos.
Observamos que as distâncias em \((1-22)\) e \((1-23)\) são completamente determinadas pelos coeficientes (pesos) \(a_{ij}\) para \(i = 1, 2, \ldots, p\) e \(j = 1/2, \ldots, p\). Esses coeficientes podem ser organizados na matriz quadrada:
\[ \left[ \begin{matrix} a_{11} & a_{12} & \cdots & a_{1p} \\ a_{21} & a_{22} & \cdots & a_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ a_{p1} & a_{p2} & \cdots & a_{pp} \end{matrix} \right] \]
em que os elementos \(a_{ij}\) com \(i \neq j\) são exibidos duas vezes, uma vez que são multiplicados por 2 nas fórmulas de distância. Consequentemente, as entradas nesta matriz especificam as funções de distância. Os \(a_{ij}\) não podem ser números arbitrários; eles devem ser tal que a distância calculada seja não negativa para cada par de pontos. (Veja o Exercício \(1.10\).)
As curvas de distâncias constantes calculadas a partir de \((1-22)\) e \((1-23)\) são hiperelipsóides. Um hiperelipsóide se assemelha a uma bola de futebol americano quando \(p = 3\); é difícil visualizar em mais de três dimensões.
Figura 1.25 - Um aglomerado de pontos em relação a um ponto P e a origem.
A necessidade de considerar a distância estatística em vez da distância euclidiana é ilustrada de forma heurística na Figura 1.25. A Figura 1.25 representa um aglomerado de pontos cujo centro de gravidade (média amostral) é indicado pelo ponto \(Q\). Considere as distâncias euclidianas do ponto \(Q\) para o ponto \(P\) e a origem \(O\). A distância euclidiana de \(Q\) para \(P\) é maior do que a distância euclidiana de \(Q\) para \(O\). No entanto, o ponto \(P\) parece ser mais semelhante aos pontos do aglomerado do que a origem. Se levamos em conta a variabilidade dos pontos no aglomerado e medimos a distância pela distância estatística em (1-20), então \(Q\) estará mais próximo de \(P\) do que de \(O\). Esse resultado parece razoável, dada a natureza da dispersão.
Outras medidas de distância podem ser consideradas. (Veja o Exercício 1.12.)
Às vezes, é útil considerar distâncias que não estão relacionadas a circunferência ou elipse. Qualquer medida de distância \(d(P,Q)\) entre dois pontos \(P\) e \(Q\) é válida, desde que satisfaça as seguintes propriedades, em que \(R\) é qualquer outro ponto intermediário:
1.4, 1.5, 1.6, 1.7, 1.8, 1.10, 1.11, 1.12, 1.14, 1.15, 1.16
Vimos no Capítulo 1 que dados multivariados podem ser convenientemente exibidos como uma matriz de números. Em geral, uma matriz retangular de números com, por exemplo, \(n\) linhas e \(p\) colunas é chamada de matriz de dimensão \(n \times p\). O estudo de métodos multivariados é grandemente facilitado pelo uso da álgebra de matrizes.
Os resultados da álgebra de matrizes apresentados neste capítulo nos permitirão declarar de forma concisa modelos estatísticos. Além disso, as relações formais expressas em termos de matrizes podem ser facilmente programadas em computadores para permitir o cálculo rotineiro de importantes quantidades estatísticas.
Começamos introduzindo alguns conceitos muito básicos que são essenciais tanto para nossas interpretações geométricas quanto para as explicações algébricas das técnicas estatísticas subsequentes. Se você ainda não teve contato com os rudimentos da álgebra de matrizes, pode preferir acompanhar o breve resumo na próxima seção antes da revisão mais detalhada fornecida no Suplemento 2A.
Seja \(\mathbf{A}\) uma matriz simétrica de ordem \(k \times k\). Então, \(\mathbf{A}\) possui \(k\) pares de autovalores e autovetores, a saber:
\[ (\lambda_1, \mathbf{e}_1) \quad (\lambda_2, \mathbf{e}_2) \quad \dots \quad (\lambda_k, \mathbf{e}_k) \tag{2-15} \]
Os autovetores podem ser escolhidos de forma que \(\mathbf{e}_1^{\prime}\mathbf{e}_1 = \cdots = \mathbf{e}_k^{\prime}\mathbf{e}_k = 1\) e sejam mutuamente ortogonais. Os autovetores são únicos, exceto quando dois ou mais autovalores forem iguais.
Considere a matriz:
\[ \mathbf{A} = \begin{bmatrix} 1 & -5 \\ -5 & 1 \end{bmatrix} \]
Temos:
\[ \mathbf{A} \begin{bmatrix} \dfrac{1}{\sqrt{2}} \\ -\dfrac{1}{\sqrt{2}} \end{bmatrix} = 6 \begin{bmatrix} \dfrac{1}{\sqrt{2}} \\ -\dfrac{1}{\sqrt{2}} \end{bmatrix} \]
Portanto, \(\lambda_1 = 6\) é um autovalor e
\[ \mathbf{e}_1 = \begin{bmatrix} \dfrac{1}{\sqrt{2}} \\ -\dfrac{1}{\sqrt{2}} \end{bmatrix} \]
é o autovetor normalizado correspondente.
Pode-se verificar que o segundo par autovalor-autovetor é \(\lambda_2 = -4\), com \(\mathbf{e}_2^{\prime} = \left[\dfrac{1}{\sqrt{2}}, \dfrac{1}{\sqrt{2}}\right]\).
# Definindo a matriz A
A <- matrix(c(1, -5,
-5, 1),
nrow = 2,
byrow = TRUE)
# Cálculo dos autovalores e autovetores
eig <- eigen(A)
# Exibindo autovalores
eig$values
[1] 6 -4
[,1] [,2]
[1,] -0.7071068 -0.7071068
[2,] 0.7071068 -0.7071068
[,1] [,2]
[1,] 1 0
[2,] 0 1
O estudo da variação e das inter-relações em dados multivariados é frequentemente baseado em distâncias e na suposição de que os dados seguem distribuição normal multivariada. As distâncias ao quadrado (ver Capítulo 1) e a densidade normal multivariada podem ser expressas em termos de produtos matriciais chamados formas quadráticas (ver Capítulo 4). Consequentemente, não deve surpreender que formas quadráticas desempenhem um papel central na análise multivariada. Nesta seção, consideramos formas quadráticas que são sempre não negativas e as matrizes definidas positivas associadas.
Resultados envolvendo formas quadráticas e matrizes simétricas são, em muitos casos, uma consequência direta de uma expansão de matriz simétrica conhecida como decomposição espectral.
A decomposição espectral de uma matriz simétrica \(\mathbf{A}\) de ordem \(k \times k\) é dada por:
\[ \begin{align} \mathbf{A} &= \lambda_1 \mathbf{e}_1 \mathbf{e}_1^{\prime} + \lambda_2 \mathbf{e}_2 \mathbf{e}_2^{\prime} + \cdots + \lambda_k \mathbf{e}_k \mathbf{e}_k^{\prime}\\ \mathbf{A} &= \sum_{i=1}^{k}{\lambda_i \mathbf{e}_i \mathbf{e}_i^{\prime}} \end{align} \]
onde \(\lambda_1, \lambda_2, \dots, \lambda_k\) são os autovalores de \(\mathbf{A}\) e \(\mathbf{e}_1, \mathbf{e}_2, \dots, \mathbf{e}_k\) são os autovetores normalizados associados (ver também o Resultado 2A.14 no Suplemento 2A). Assim, \(\mathbf{e}_i^{\prime} \mathbf{e}_i = 1\) para \(i = 1, 2, \dots, k\) e \(\mathbf{e}_i^{\prime} \mathbf{e}_j = 0\) para \(i \ne j\).
A decomposição espectral é uma ferramenta analítica importante. Com ela, é possível demonstrar com facilidade certos resultados estatísticos. O primeiro deles é uma explicação matricial da distância, que será apresentada a seguir.
Como a forma \(\mathbf{x}^{\prime}\mathbf{A}\mathbf{x}\) envolve apenas termos quadráticos \(x_i^2\) e produtos \(x_i x_k\), ela é chamada de forma quadrática. Quando uma matriz simétrica \(\mathbf{A}\) de ordem \(k \times k\) satisfaz:
\[ 0 \leq \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} \tag{2-17} \]
para todo vetor \(\mathbf{x}^{\prime} = [x_1, x_2, \dots, x_k]\), tanto a matriz \(\mathbf{A}\) quanto a forma quadrática são ditas definidas não negativas. Se a igualdade na equação (2-17) ocorre apenas para \(\mathbf{x}^{\prime} = [0, 0, \dots, 0]\), então \(\mathbf{A}\) ou a forma quadrática é dita definida positiva. Em outras palavras, \(\mathbf{A}\) é definida positiva se:
\[ 0 < \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} \tag{2-18} \]
para todo vetor \(\mathbf{x} \ne \mathbf{0}\).
Considere a matriz simétrica:
\[ \mathbf{A} = \begin{bmatrix} 1 & 3 & 2 \\ -4 & 13 & -2 \\ 2 & -2 & 10 \end{bmatrix} \]
Os autovalores obtidos da equação característica \(|\mathbf{A} - \lambda \mathbf{I}| = 0\) são: \(\lambda_1 = 9\), \(\lambda_2 = 9\) e \(\lambda_3 = 18\).
Os autovetores correspondentes \(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3\) são soluções normalizadas de \(\mathbf{A} \mathbf{e}_i = \lambda_i \mathbf{e}_i\), para \(i = 1, 2, 3\).
A equação \(\mathbf{A} \mathbf{e}_1 = 9 \mathbf{e}_1\) fornece:
\[ \begin{bmatrix} 13 & -4 & 2 \\ -4 & 13 & -2 \\ 2 & -2 & 10 \end{bmatrix} \begin{bmatrix} e_{11} \\ e_{21} \\ e_{31} \end{bmatrix} = 9 \begin{bmatrix} e_{11} \\ e_{21} \\ e_{31} \end{bmatrix} \]
ou seja:
\[ \begin{align} 13e_{11} - 4e_{21} + 2e_{31} &= 9e_{11} \\ -4e_{11} + 13e_{21} - 2e_{31} &= 9e_{21} \\ 2e_{11} - 2e_{21} + 10e_{31} &= 9e_{31} \end{align} \]
Subtraindo os termos da direita e escolhendo arbitrariamente \(e_{11} = 1\) e \(e_{21} = 1\), obtém-se \(e_{31} = 0\). Assim, o autovetor normalizado é:
\[ \mathbf{e}_1^{\prime} = \left[ \dfrac{1}{\sqrt{2}},\ \dfrac{1}{\sqrt{2}},\ 0 \right] \]
Pode-se verificar que
\[ \mathbf{e}_2^{\prime} = \left[ \dfrac{1}{\sqrt{18}},\ -\dfrac{1}{\sqrt{18}},\ -\dfrac{4}{\sqrt{18}} \right] \]
também é autovetor para \(\lambda = 9\), e que
\[ \mathbf{e}_3^{\prime} = \left[ \dfrac{2}{3},\ -\dfrac{2}{3},\ \dfrac{1}{3} \right] \]
é o autovetor normalizado correspondente a \(\lambda = 18\).
Ademais, \(\mathbf{e}_i^{\prime} \mathbf{e}_j = 0\) para \(i \ne j\).
A decomposição espectral de \(\mathbf{A}\) é dada por:
\[ \mathbf{A} = \lambda_1 \mathbf{e}_1 \mathbf{e}_1^{\prime} + \lambda_2 \mathbf{e}_2 \mathbf{e}_2^{\prime} + \lambda_3 \mathbf{e}_3 \mathbf{e}_3^{\prime} \]
ou seja:
\[ \begin{align} \begin{bmatrix} 13 & -4 & 2 \\ -4 & 13 & -2 \\ 2 & -2 & 10 \end{bmatrix} &= 9 \begin{bmatrix} \dfrac{1}{\sqrt{2}} \\ \dfrac{1}{\sqrt{2}} \\ 0 \end{bmatrix} \begin{bmatrix} \dfrac{1}{\sqrt{2}} & \dfrac{1}{\sqrt{2}} & 0 \end{bmatrix} +\\ &\quad 9 \begin{bmatrix} \dfrac{1}{\sqrt{18}} \\ -\dfrac{1}{\sqrt{18}} \\ -\dfrac{4}{\sqrt{18}} \end{bmatrix} \begin{bmatrix} \dfrac{1}{\sqrt{18}} & -\dfrac{1}{\sqrt{18}} & -\dfrac{4}{\sqrt{18}} \end{bmatrix} +\\ &\quad 18 \begin{bmatrix} \dfrac{2}{3} \\ -\dfrac{2}{3} \\ \dfrac{1}{3} \end{bmatrix} \begin{bmatrix} \dfrac{2}{3} & -\dfrac{2}{3} & \dfrac{1}{3} \end{bmatrix} \end{align} \]
Multiplicando os vetores e somando os produtos:
\[ \begin{align} \begin{bmatrix} 13 & -4 & 2 \\ -4 & 13 & -2 \\ 2 & -2 & 10 \end{bmatrix} &= 9 \begin{bmatrix} \dfrac{1}{2} & \dfrac{1}{2} & 0 \\ \dfrac{1}{2} & \dfrac{1}{2} & 0 \\ 0 & 0 & 0 \end{bmatrix} +\\ &\quad 9 \begin{bmatrix} \dfrac{1}{18} & -\dfrac{1}{18} & -\dfrac{4}{18} \\ -\dfrac{1}{18} & \dfrac{1}{18} & \dfrac{4}{18} \\ -\dfrac{4}{18} & \dfrac{4}{18} & \dfrac{16}{18} \end{bmatrix} +\\ &\quad 18 \begin{bmatrix} \dfrac{4}{9} & -\dfrac{4}{9} & \dfrac{2}{9} \\ -\dfrac{4}{9} & \dfrac{4}{9} & -\dfrac{2}{9} \\ \dfrac{2}{9} & -\dfrac{2}{9} & \dfrac{1}{9} \end{bmatrix} \end{align} \]
Como se pode verificar diretamente, a soma dos três termos reproduz a matriz original \(\mathbf{A}\).
# Definindo a matriz A
A <- matrix(c(13, -4, 2,
-4, 13, -2,
2, -2, 10), nrow = 3, byrow = TRUE)
# Cálculo dos autovalores e autovetores normalizados
eig <- eigen(A)
# Autovalores
eig$values
[1] 18 9 9
[,1] [,2] [,3]
[1,] 0.67 -0.75 0.00
[2,] -0.67 -0.60 0.45
[3,] 0.33 0.30 0.89
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
\[\Diamond\]
Mostrar que a matriz aimétrica associada à seguinte forma quadrática é definida positiva:
\[ 3x_1^2 + 2x_2^2 - 2\sqrt{2}x_1 x_2 \]
Para ilustrar a abordagem geral, primeiro reescrevemos a forma quadrática na notação matricial:
\[ \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} 3 & -\sqrt{2} \\ -\sqrt{2} & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} \]
Pela Definição 2A.30, os autovalores de \(\mathbf{A}\) são soluções da equação \(|\mathbf{A} - \lambda \mathbf{I}| = 0\), ou seja, \((3 - \lambda)(2 - \lambda) - 2 = 0\). As soluções são \(\lambda_1 = 4\), \(\lambda_2 = 1\).
Usando a decomposição espectral (equação 2-16), temos:
\[ \mathbf{A} = \lambda_1 \mathbf{e}_1 \mathbf{e}_1^{\prime} + \lambda_2 \mathbf{e}_2 \mathbf{e}_2^{\prime} = 4 \mathbf{e}_1 \mathbf{e}_1^{\prime} + \mathbf{e}_2 \mathbf{e}_2^{\prime} \]
onde \(\mathbf{e}_1\) e \(\mathbf{e}_2\) são autovetores normalizados e ortogonais associados a \(\lambda_1 = 4\) e \(\lambda_2 = 1\), respectivamente.
Como 4 e 1 são escalares, premultiplicando e postmultiplicando \(\mathbf{A}\) por \(\mathbf{x}^{\prime}\) e \(\mathbf{x}\), com \(\mathbf{x}^{\prime} = [x_1, x_2]\), temos:
\[ \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} = 4 \mathbf{x}^{\prime} \mathbf{e}_1 \mathbf{e}_1^{\prime} \mathbf{x} + \mathbf{x}^{\prime} \mathbf{e}_2 \mathbf{e}_2^{\prime} \mathbf{x} = 4y_1^2 + y_2^2 \geq 0 \]
com:
\[ y_1 = \mathbf{x}^{\prime} \mathbf{e}_1 = \mathbf{e}_1^{\prime} \mathbf{x}, \quad y_2 = \mathbf{x}^{\prime} \mathbf{e}_2 = \mathbf{e}_2^{\prime} \mathbf{x} \]
Mostramos agora que \(y_1\) e \(y_2\) não podem ser ambos nulos, e portanto:
\[ \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} = 4y_1^2 + y_2^2 > 0, \]
ou seja, \(\mathbf{A}\) é definida positiva.
Pelas definições de \(y_1\) e \(y_2\):
\[ \begin{align} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} &= \begin{bmatrix} \mathbf{e}_1^{\prime} \\ \mathbf{e}_2^{\prime} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}\\ \mathbf{y} &= \mathbf{E} \mathbf{x} \end{align} \]
A matriz \(\mathbf{E}\) é ortogonal e, portanto, possui inversa \(\mathbf{E}^{\prime}\). Logo, \(\mathbf{x} = \mathbf{E}^{\prime} \mathbf{y}\). Porém, \(\mathbf{x}\) é vetor não nulo, então \(\mathbf{y} \ne \mathbf{0}\), o que implica \(y_1 \ne 0\) ou \(y_2 \ne 0\).
# Definindo a matriz A da forma quadrática
A <- matrix(c(3, -sqrt(2),
-sqrt(2), 2), nrow = 2, byrow = TRUE)
# Verificando se é definida positiva via autovalores
eig <- eigen(A)
# Autovalores
eig$values
[1] 4 1
[1] TRUE
# Gerando a matriz E de autovetores (ortogonal)
E <- eig$vectors
# Testando a identidade E %*% t(E) = I
E %*% t(E)
[,1] [,2]
[1,] 1 0
[2,] 0 1
# Testando a forma quadrática para vetor não nulo
x <- c(1, 1) # vetor arbitrário não nulo
print(q <- as.numeric(t(x) %*% A %*% x)) # deve ser positivo
[1] 2.171573
[1] TRUE
\[\Diamond\]
Usando a decomposição espectral, podemos facilmente mostrar que uma matriz simétrica \(\mathbf{A}\) de ordem \(k \times k\) é definida positiva se, e somente se, todos os autovalores de \(\mathbf{A}\) forem positivos (ver Exercício 2.17). A matriz \(\mathbf{A}\) é definida não negativa se, e somente se, todos os seus autovalores forem maiores ou iguais a zero.
Admita, por um momento, que os \(p\) elementos \(x_1, x_2, \dots, x_p\) de um vetor \(\mathbf{x}\) são realizações de \(p\) variáveis aleatórias \(X_1, X_2, \dots, X_p\). Como indicado no Capítulo 1, podemos interpretar esses elementos como coordenadas de um ponto em um espaço \(p\)-dimensional. A “distância” do ponto \([x_1, x_2, \dots, x_p]\) até a origem pode, neste caso, ser interpretada em termos de desvios-padrão. Assim, conseguimos incorporar a incerteza inerente (variabilidade) nas observações. Pontos com a mesma “incerteza” são tratados como estando à mesma distância da origem.
Se utilizarmos a fórmula da distância apresentada no Capítulo 1 [ver Equação (1-22)], a distância da origem até os pontos satisfaz a forma geral:
\[ \begin{align} \text{distância}^2 &= a_{11}x_1^2 + a_{22}x_2^2 + \cdots + a_{pp}x_p^2 + \\ &\quad 2(a_{12}x_1x_2 + a_{13}x_1x_3 + \cdots + a_{p-1,p}x_{p-1}x_p) \end{align} \]
Assumindo que \(\text{distância}^2 > 0\) para todo vetor \([x_1, x_2, \dots, x_p] \ne [0, 0, \dots, 0]\) e que \(a_{ij} = a_{ji}\) para \(i \ne j\), temos:
\[ 0 < \text{distância}^2 = \begin{bmatrix} x_1 & x_2 & \cdots & x_p \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1p} \\ a_{21} & a_{22} & \cdots & a_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ a_{p1} & a_{p2} & \cdots & a_{pp} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_p \end{bmatrix} \]
ou seja:
\[ 0 < \text{distância}^2 = \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} \quad \text{para } \mathbf{x} \ne \mathbf{0} \tag{2-19} \]
A partir da equação (2-19), vemos que a matriz simétrica \(\mathbf{A}\) de ordem \(p \times p\) é definida positiva. Em resumo, a distância é determinada por uma forma quadrática definida positiva \(\mathbf{x}^{\prime} \mathbf{A} \mathbf{x}\). Reciprocamente, uma forma quadrática definida positiva pode ser interpretada como uma distância ao quadrado.
Comentário. Seja o quadrado da distância do ponto \(\mathbf{x}^{\prime} = [x_1, x_2, \dots, x_p]\) até a origem dado por \(\mathbf{x}^{\prime} \mathbf{A} \mathbf{x}\), onde \(\mathbf{A}\) é uma matriz simétrica \(p \times p\) definida positiva.
Então, o quadrado da distância de \(\mathbf{x}\) até um ponto fixo arbitrário \(\boldsymbol{\mu}^{\prime} = [\mu_1, \mu_2, \dots, \mu_p]\) é dado pela expressão geral \((\mathbf{x} - \boldsymbol{\mu})^{\prime} \mathbf{A} (\mathbf{x} - \boldsymbol{\mu})\).
Expressar a distância como raiz quadrada de uma forma quadrática definida positiva permite uma interpretação geométrica baseada nos autovalores e autovetores da matriz \(\mathbf{A}\). Por exemplo, suponha \(p = 2\). Então, os pontos \(\mathbf{x}^{\prime} = [x_1, x_2]\) a uma distância constante \(c\) da origem satisfazem:
\[ \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} = a_{11}x_1^2 + a_{22}x_2^2 + 2a_{12}x_1x_2 = c^2 \]
Pela decomposição espectral, como no Exemplo 2.11,
\[ \mathbf{A} = \lambda_1 \mathbf{e}_1 \mathbf{e}_1^{\prime} + \lambda_2 \mathbf{e}_2 \mathbf{e}_2^{\prime} \quad \text{então} \quad \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} = \lambda_1 (\mathbf{x}^{\prime} \mathbf{e}_1)^2 + \lambda_2 (\mathbf{x}^{\prime} \mathbf{e}_2)^2 \]
Assim, \(c^2 = \lambda_1 y_1^2 + \lambda_2 y_2^2\) representa uma elipse em \(y_1 = \mathbf{x}^{\prime} \mathbf{e}_1\) e \(y_2 = \mathbf{x}^{\prime} \mathbf{e}_2\), pois \(\lambda_1, \lambda_2 > 0\) quando \(\mathbf{A}\) é definida positiva (ver Exercício 2.17).
Podemos verificar facilmente que \(\mathbf{x} = c \lambda_1^{-1/2} \mathbf{e}_1\) satisfaz:
\[ \mathbf{x}^{\prime} \mathbf{A} \mathbf{x} = \lambda_1 (c \lambda_1^{-1/2} \mathbf{e}_1^{\prime} \mathbf{e}_1)^2 = c^2 \]
Similarmente, \(\mathbf{x} = c \lambda_2^{-1/2} \mathbf{e}_2\) fornece a distância apropriada na direção de \(\mathbf{e}_2\). Assim, os pontos a uma distância \(c\) da origem estão em uma elipse cujos eixos são dados pelos autovetores de \(\mathbf{A}\), com comprimentos proporcionais aos recíprocos das raízes quadradas dos autovalores. A constante de proporcionalidade é \(c\).
Se \(p > 2\), os pontos \(\mathbf{x}^{\prime} = [x_1, x_2, \dots, x_p]\) a uma distância constante \(c = \sqrt{\mathbf{x}^{\prime} \mathbf{A} \mathbf{x}}\) da origem estão em hiperelipsoides:
\[ c^2 = \lambda_1 (\mathbf{x}^{\prime} \mathbf{e}_1)^2 + \cdots + \lambda_p (\mathbf{x}^{\prime} \mathbf{e}_p)^2 \]
cujos eixos são dados pelos autovetores de \(\mathbf{A}\). O semi-eixo na direção de \(\mathbf{e}_i\) é igual a \(c / \sqrt{\lambda_i}\), para \(i = 1, 2, \dots, p\), onde \(\lambda_1, \lambda_2, \dots, \lambda_p\) são os autovalores de \(\mathbf{A}\).
Lugar geométrico de distância estatística constante c a partir da origem (p = 2, 1 ≤ λ1 < λ2)
A decomposição espectral permite expressar a inversa de uma matriz quadrada em termos de seus autovalores e autovetores, o que leva à definição útil de matriz raiz quadrada.
Seja \(\mathbf{A}\) uma matriz simétrica definida positiva de ordem \(k \times k\), com decomposição espectral dada por:
\[ \mathbf{A} = \sum_{i=1}^{k} \lambda_i \mathbf{e}_i \mathbf{e}_i^{\prime} \]
Sejam os autovetores normalizados as colunas da matriz:
\[ \mathbf{P} = [\mathbf{e}_1, \mathbf{e}_2, \dots, \mathbf{e}_k] \]
Então:
\[ \mathbf{A} = \sum_{i=1}^{k} \lambda_i \mathbf{e}_i \mathbf{e}_i^{\prime} = \mathbf{P} \boldsymbol{\Lambda} \mathbf{P}^{\prime} \tag{2-20} \]
onde \(\mathbf{P} \mathbf{P}^{\prime} = \mathbf{P}^{\prime} \mathbf{P} = \mathbf{I}\) e \(\boldsymbol{\Lambda}\) é a matriz diagonal:
\[ \underset{k\times k}{\boldsymbol{\Lambda}} = \begin{bmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_k \end{bmatrix} \quad \text{com } \lambda_i > 0 \]
Portanto, a inversa de \(\mathbf{A}\) é dada por:
\[ \mathbf{A}^{-1} = \mathbf{P} \boldsymbol{\Lambda}^{-1} \mathbf{P}^{\prime} = \sum_{i=1}^{k} \dfrac{1}{\lambda_i} \mathbf{e}_i \mathbf{e}_i^{\prime} \tag{2-21} \]
Agora, seja \(\boldsymbol{\Lambda}^{1/2}\) a matriz diagonal com \(\sqrt{\lambda_i}\) na diagonal. A matriz
\[ \sum_{i=1}^{k} \sqrt{\lambda_i} \mathbf{e}_i \mathbf{e}_i^{\prime} = \mathbf{P} \boldsymbol{\Lambda}^{1/2} \mathbf{P}^{\prime} \tag{2-22} \]
é chamada de raiz quadrada de \(\mathbf{A}\) e denotada por \(\mathbf{A}^{1/2}\).
A matriz \(\mathbf{A}^{1/2}\) possui as seguintes propriedades:
Um vetor aleatório é um vetor cujos elementos são variáveis aleatórias. Da mesma forma, uma matriz aleatória é uma matriz cujos elementos são variáveis aleatórias. O valor esperado de uma matriz (ou vetor) aleatório é a matriz (ou vetor) composto pelos valores esperados de cada um dos seus elementos. Especificamente, seja \(\boldsymbol{\mathcal{X}}\) uma matriz aleatória de dimensão \(n \times p\).
\[ \boldsymbol{\mathcal{X}} = \begin{bmatrix} X_{11} & X_{12} & \ldots & X_{1p} \\ X_{21} & X_{22} & \ldots & X_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ X_{n1} & X_{n2} & \ldots & X_{np} \\ \end{bmatrix} \]
Então, o valor esperado de \(\boldsymbol{\mathcal{X}}\), denotado por \(\mathbb{E}(\boldsymbol{\mathcal{X}})\), é a matriz \(n \times p\) composta por números (se existirem), conforme a seguir:
\[ \mathbb{E}(\boldsymbol{\mathcal{X}}) = \begin{bmatrix} \mathbb{E}(X_{11}) & \mathbb{E}(X_{12}) & \ldots & \mathbb{E}(X_{1p}) \\ \mathbb{E}(X_{21}) & \mathbb{E}(X_{22}) & \ldots & \mathbb{E}(X_{2p}) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbb{E}(X_{n1}) & \mathbb{E}(X_{n2}) & \ldots & \mathbb{E}(X_{np}) \\ \end{bmatrix} \]
Os resultados envolvendo o valor esperado das somas e produtos de matrizes seguem diretamente da definição do valor esperado de uma matriz aleatória e das propriedades univariadas da esperança, \(\mathbb{E}(X_1 + X_2) = \mathbb{E}(X_1) + \mathbb{E}(X_2)\) e \(\mathbb{E}(cX_1) = c\mathbb{E}(X_1)\), em que \(c\) é uma constante.
Sejam \(\boldsymbol{\mathcal{X}}\) e \(\boldsymbol{\mathcal{Y}}\) matrizes aleatórias de mesma dimensão, e sejam \(\mathbf{A}\) e \(\mathbf{B}\) matrizes de constantes conformáveis. Então, temos:
\[ \begin{align} \mathbb{E}(\boldsymbol{\mathcal{X}} + \boldsymbol{\mathcal{Y}}) &= \mathbb{E}(\boldsymbol{\mathcal{X}}) + \mathbb{E}(\boldsymbol{\mathcal{Y}})\\ \mathbb{E}(\mathbf{A\boldsymbol{\mathcal{X}}B}) &= \mathbf{A}\mathbb{E}(\boldsymbol{\mathcal{X}})\mathbf{B} \end{align} \tag{2-24} \]
Suponha que \(\mathbf{X}_{i}^{\prime} = [X_{i1}, X_{i2}, \ldots, X_{ip}]\) seja um vetor aleatório de dimensão \(p \times 1\) representando a \(i\)-ésima observação multivariada, \(i=1,2, \ldots,n\). Cada elemento de \(\mathbf{X}_i\) é uma variável aleatória com sua própria distribuição de probabilidade marginal. (Veja o Exemplo 2.12.)
As médias marginais e as covariâncias \(\sigma_{ji,ki}\) são definidas como \(\mu_{ki} = \mathbb{E}(X_{ki})\) e \(\sigma_{ji,ki} = \mathbb{E}((X_{ji} - \mu_{ji})(X_{ki} - \mu_{ki}))\) para \(i = 1, 2, \ldots, p\), respectivament, por integral de Lebesgue:
\[ \begin{align} \mu_{ki}&=\int_{\mathbb{R}}{x_{ki}dF(x_{ki})} \\ \sigma_{ki}^{2}&=\int_{\mathbb{R}}{(x_{ki}-\mu_{ki})^{2}dF(x_{ki})} \end{align} \tag{2-25} \]
O comportamento de qualquer par de variáveis aleatórias (\(i\) e \(k\)) de uma unidade observacional (\(j\)), como \(X_{ji}\) e \(X_{jk}\), é descrito por sua função de probabilidade conjunta, e uma medida da associação linear entre elas é fornecida pela covariância:
\[ \begin{align} \mathbb{C}(X_{ji}, X_{jk}) &= \mathbb{E}((X_{ji} - \mu_{ji})(X_{jk} - \mu_{jk}))\\ \mathbb{C}(X_{ji}, X_{jk}) &=\sigma_{ji, jk} \end{align} \]
em que \(\mu_{ji}\) e \(\mu_{jk}\) são as médias marginais, para \(i, k = 1, 2, \ldots, p\). Quando \(i = k\), a covariância se torna a variância marginal.
\[ \begin{align} \mathbb{C}(X_{jk}, X_{jk}) &= \mathbb{E}((X_{jk} - \mu_{jk})(X_{jk} - \mu_{jk}))\\ &= \mathbb{E}((X_{jk} - \mu_{jk})^2)\\ &=\mathbb{V}(X_{jk})\\ \mathbb{C}(X_{jk}, X_{jk})&=\sigma_{jk, jk} \end{align} \]
De maneira mais geral, desconsiderando a unidade observacional, o comportamento coletivo das \(p\) variáveis aleatórias \(X_1, X_2, \ldots, X_p\) ou, equivalentemente, o vetor aleatório \(\mathbf{X}^{\prime} = [X_1\; X_2\; \ldots\; X_p]\), é descrito por uma função de densidade de probabilidade conjunta \(f(x_1, x_2, \ldots, x_p) = f(\mathbf{x})\). Como já mencionado anteriormente neste livro, \(f(\mathbf{x})\) frequentemente será a função de densidade multivariada normal. (Veja o Capítulo 4.)
Se a probabilidade conjunta \(P(X_i = x_i, X_k = x_k)\) pode ser escrita como o produto das probabilidades marginais correspondentes, de modo que:
\[ P(X_i = x_i, X_k = x_k) = P(X_i = x_i)P(X_k = x_k) \tag{2-27} \]
para todos os pares de valores \(x_i, x_k\), as variáveis aleatórias \(X_i\) e \(X_k\) são ditas estatisticamente independentes. Quando \(X_i\) e \(X_k\) são variáveis aleatórias contínuas com densidades conjuntas \(f_{ik}(x_i, x_k)\) e densidades marginais \(f_i(x_i)\) e \(f_k(x_k)\), a condição de independência se torna:
\[ f_{ik}(x_i, x_k) = f_i(x_i)f_k(x_k) \]
para todos os pares \((x_i, x_k)\).
As \(p\) variáveis aleatórias contínuas \(X_1, X_2, \ldots, X_p\) são estatisticamente independentes mutuamente se sua densidade conjunta puder ser fatorada da seguinte forma:
\[ \begin{align} f_{12\cdots p}(x_1, x_2, \ldots, x_p) &= f_1(x_1)f_2(x_2)\cdots f_p(x_p)\\ f_{12\cdots p}(x_1, x_2, \ldots, x_p) &= \prod_{i=1}^{p}{f_i(x_i)} \end{align} \]
para todas as \(p\)-tuplas \((x_1, x_2, \ldots, x_p)\).
A independência estatística tem uma importante implicação para a covariância. A factorização em (2-28) implica que \(\mathbb{C}(X_i, X_k) = 0\). Portanto,
\[ \mathbb{C}(X_i, X_k) = 0 \tag{2-29} \]
se \(X_i\) e \(X_k\) são independentes.
A recíproca de (2-29) não é verdadeira em geral; existem situações em que \(\mathbb{C}(X_i, X_k) = 0\), mas \(X_i\) e \(X_k\) não são independentes. (Ver [5].)
As médias e covariâncias populacionais do vetor aleatório \(\mathbf{X}\) de dimensão \(p \times 1\) podem ser representadas em matrizes.
O valor esperado de cada elemento é contido no vetor de médias \(\mathbb{E}(\mathbf{X}_i)\), \(i=1,2,\ldots,n\), e as \(p\) variâncias \(\sigma_{ii}=\sigma_i^2\) e as \(p(p - 1)/2\) covariâncias distintas \(\sigma_{ik}\) estão contidas na matriz de covariância simétrica \(\mathbf{\Sigma}_i\):
\[ \boldsymbol{\mu}_i = \mathbb{E}(\mathbf{X}_i) = \begin{bmatrix} \mathbb{E}(X_{i1}) \\ \mathbb{E}(X_{i2}) \\ \vdots \\ \mathbb{E}(X_{ip}) \end{bmatrix} \]
Os vetores de média (centróide) populacional podem ser diferentes nas \(n\) unidades observacionais.
\[ \mathbf{\Sigma}_i = \mathbb{C}(\mathbf{X_i})=\mathbb{E}(\mathbf{X_i}-\boldsymbol{\mu}_i)(\mathbf{X_i}-\boldsymbol{\mu}_i)^{\prime}= \begin{bmatrix} \sigma_{i1,i1} & \sigma_{i1,i2} & \ldots & \sigma_{i1,ip} \\ \sigma_{i2,i1} & \sigma_{i2,i2} & \ldots & \sigma_{i2,ip} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{ip,i1} & \sigma_{ip,i2} & \ldots & \sigma_{ip,ip} \end{bmatrix} \]
As covariâncias populacionais são dadas por:
\[ \sigma_{ij,ik} = E[(X_{ij} - \mu_{ij})(X_{ik} - \mu_{ik})]\\\\ j, k = 1, 2, \ldots, p \\ i=1,2,\ldots, n \]
A condição de homocedasticidade multivariada é:
\[ \mathbf{\Sigma}_1=\mathbf{\Sigma}_2=\cdots =\mathbf{\Sigma}_n=\mathbf{\Sigma} \]
A distribuição normal multivariada é completamente especificada quando o vetor médio \(\boldsymbol{\mu}\) e a matriz de covariância \(\mathbf{\Sigma}\) são dados (ver Capítulo 4), portanto, não é surpreendente que essas quantidades desempenhem um papel importante em muitos procedimentos multivariados.
Frequentemente, é informativo separar a informação contida nas variâncias \(\sigma_{ii}\) daquela contida nas medidas de associação, em particular, o coeficiente de correlação populacional \(\rho_{ik}\). O coeficiente de correlação de Pearson, \(\rho_{ik}\), é definido em termos da covariância \(\sigma_{ik}\) e das variâncias \(\sigma_{ii}\) e \(\sigma_{kk}\) como:
\[ \rho_{ik} = \dfrac{\sigma_{ik}}{\sqrt{\sigma_{ii} \sigma_{kk}}} \tag{2-33} \in [-1,1] \]
O coeficiente de correlação mede o grau de associação linear entre as variáveis aleatórias \(X_i\) e \(X_k\). (Ver, por exemplo, [5].)
Seja a matriz de correlação populacional simétrica \(p \times p\):
\[ \boldsymbol{\rho}=\begin{bmatrix} 1 & \rho_{12} & \cdots & \rho_{1p} \\ \rho_{21} & 1 & \cdots & \rho_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{p1} & \rho_{p2} & \cdots & 1 \end{bmatrix} \tag{2-34} \]
e seja a matriz de desvio-padrão \(p \times p\) ser:
\[ \mathbf{V}^{1/2} = \begin{bmatrix} \sqrt{\sigma_{11}} & 0 & \cdots & 0 \\ 0 & \sqrt{\sigma_{22}} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sqrt{\sigma_{pp}} \end{bmatrix} \tag{2-36} \]
Então, é facilmente verificado (ver Exercício 2.23) que:
\[ \mathbf{V}^{1/2} \boldsymbol{\rho} \; \mathbf{V}^{1/2}=\mathbf{\Sigma} \tag{2-36} \]
e
\[ \boldsymbol{\rho} = \left(\mathbf{V}^{1/2}\right)^{-1} \mathbf{\Sigma}\; \left(\mathbf{V}^{1/2}\right)^{-1} \tag{2-37} \]
Ou seja, \(\mathbf{\Sigma}\) pode ser obtida a partir de \(\mathbf{V}^{1/2}\) e \(\boldsymbol{\rho}\), enquanto \(\boldsymbol{\rho}\) pode ser obtida a partir de \(\mathbf{X}_i\). Além disso, expressar essas relações em termos de operações matriciais permite que os cálculos sejam convenientemente implementados em um computador.
Frequentemente, as características medidas em experimentos individuais naturalmente se dividem em dois ou mais grupos. Como exemplos, considere medições de variáveis que representam consumo e renda, ou variáveis que representam traços de personalidade e características físicas. Uma abordagem para lidar com essas situações é permitir que as características que definem os grupos distintos sejam subconjuntos do conjunto total de características. Se o conjunto total é representado por um vetor aleatório de dimensão \(p \times 1\), \(\mathbf{X}_i\), \(i=1,2\ldots,n\), os subconjuntos podem ser considerados como componentes de \(\mathbf{X}_i\) e podem ser organizados por meio da partição de \(\mathbf{X}_i\).
Em geral, podemos particionar as \(p\) características contidas no vetor aleatório \(\mathbf{X}\) de dimensão \(p \times 1\) em, por exemplo, dois grupos de tamanho \(q\) e \(p - q\), respectivamente. Por exemplo, podemos escrever:
Observe que \(\mathbf{\Sigma}_{12} = \mathbf{\Sigma}_{21}^{\prime}\). A matriz de covariância de \(\mathbf{X}^{(1)}\) é \(\mathbf{\Sigma}_{11}\), a de \(\mathbf{X}^{(2)}\) é \(\mathbf{\Sigma}_{22}\), e a dos elementos de \(\mathbf{X}^{(1)}\) e \(\mathbf{X}^{(2)}\) é \(\mathbf{\Sigma}_{12}\) (ou \(\mathbf{\Sigma}_{21}\)).
Às vezes, é conveniente usar a notação
\[ \mathbb{C}\left(\mathbf{X}^{(1)},\mathbf{X}^{(2)}\right)= \mathbf{\Sigma}_{12} \]
Lembramos que se uma única variável aleatória, como \(X_1\), é multiplicada por uma constante \(c\), então
\[ \mathbb{E}(cX_1) = c\mathbb{E}(X_1) = c\mu_1 \]
e
\[ \mathbb{V}(cX_1) = E\left((cX_1 - c\mu_1)^2\right) = c^2\mathbb{V}(X_1) \]
Se \(X_2\) é uma segunda variável aleatória e \(a\) e \(b\) são constantes, então, usando propriedades adicionais da esperança, temos
\[ \begin{align} \mathbb{C}(aX_1, bX_2) &= E\left((aX_1 - a\mu_i)(bX_2 - b\mu_2)\right) \\ &= abE\left(X_1 - \mu_1)(X_2 - \mu_2)\right) \\ &= ab\mathbb{C}(X_1, X_2) \\ \mathbb{C}(aX_1, bX_2) &= ab\sigma_{12} \end{align} \]
Finalmente, para a combinação linear \(aX_1 + bX_2\), temos
\[ \begin{align} \mathbb{E}(aX_1 + bX_2) &= a\mathbb{E}(X_1) + b\mathbb{E}(X_2) \\ \mathbb{E}(aX_1 + bX_2) &= a\mu_1 + b\mu_2 \end{align} \]
\[ \begin{align} \mathbb{V}(aX_1 + bX_2) &= a^2\mathbb{V}(X_1) + 2ab\mathbb{C}(X_1, X_2) + b^2\mathbb{V}(X_2)\\ \mathbb{V}(aX_1 + bX_2) &= a^2\sigma_{11} + 2ab\sigma_{12} + b^2\sigma_{22} \end{align} \tag{2-41} \]
Com \(\mathbf{c}^{\prime} = [a \; b]\), \(aX_1 + bX_2\) pode ser escrito como:
\[ \mathbf{c}^{\prime}\mathbf{X}= \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} X_1 \\ X_2 \end{bmatrix} \]
Da mesma forma, \(\mathbb{E}(aX_1 + bX_2) = a\mu_1 + b\mu_2\) pode ser expresso como:
\[ \mathbf{c}^{\prime}\boldsymbol{\mu}= \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix} \]
Se
\[ \mathbf{\Sigma}= \begin{bmatrix} \sigma_{11} & \sigma_{12}\\ \sigma_{21} & \sigma_{22} \end{bmatrix} \]
é a matriz de covariância de \(\mathbf{X}\), a equação (2-41) torna-se:
\[ \mathbb{V}(aX_1 + bX_2) = \mathbb{V}(\mathbf{c}^{\prime}\mathbf{X}) = \mathbf{c}^{\prime}\mathbf{\Sigma} \mathbf{c} \tag{2-42} \]
Os resultados anteriores podem ser estendidos para uma combinação linear de \(p\) variáveis aleatórias:
A combinação linear \(\mathbf{c}^{\prime}\mathbf{X} = c_1X_1 + \cdots + c_pX_p\) tem
média populacional\(= \mathbb{E}(\mathbf{c}^{\prime}\mathbf{X}) = \mathbf{c}^{\prime}\boldsymbol{\mu}\)
covariância populacional\(= \mathbb{V}(\mathbf{c}^{\prime}\mathbf{X}) = \mathbf{c}^{\prime}\mathbf{\Sigma} \mathbf{c} \tag{2-43}\)
em que \(\boldsymbol{\mu} = \mathbb{E}(\mathbf{X})\) e \(\mathbf{\Sigma} = \mathbb{C}(\mathbf{X})\).
Em geral, considere as \(q\) combinações lineares das \(p\) variáveis aleatórias:
\[ \begin{align} Z_1 &= c_{11}X_1 + c_{12}X_2 + \ldots + c_{1p}X_p\\ Z_2 &= c_{21}X_1 + c_{22}X_2 + \ldots + c_{2p}X_p\\ \vdots\\ Z_q &= c_{q1}X_1 + c_{q2}X_2 + \ldots + c_{qp}X_p \end{align} \]
ou
\[ \mathbf{Z} = \mathbf{CX} \]
em que
\[ \mathbf{C} = \begin{bmatrix} c_{11} & c_{12} & \ldots & c_{1p} \\ c_{21} & c_{22} & \ldots & c_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ c_{q1} & c_{q2} & \ldots & c_{qp} \end{bmatrix} \]
As combinações lineares \(\mathbf{Z} = \mathbf{CX}\) têm:
\[ \boldsymbol{\mu}_{\mathbf{Z}} = \mathbb{E}(\mathbf{Z}) = \mathbb{E}(\mathbf{CX}) = \mathbf{C\boldsymbol{\mu}_{\mathbf{X}}} \]
\[ \mathbf{\Sigma_Z} = \mathbb{C}(\mathbf{Z}) = \mathbb{C}(\mathbf{CX}) = \mathbf{C\Sigma_X C}^{\prime} \tag{2-45} \]
Suponha que \(\boldsymbol{\mu}_{\mathbf{X}}\) e \(\mathbf{\Sigma_X}\) sejam o vetor de média e a matriz covariância de \(\mathbf{X}\), respectivamente, (veja o Exercício 2.28 para o cálculo dos termos fora da diagonal em \(\mathbf{C\Sigma_XC}^{\prime}\)). Iremos depender fortemente do resultado em (2-45) em nossas discussões sobre componentes principais e análise de fatores nos Capítulos 8 e 9.
Seja \(\mathbf{X}^{\prime} = [X_1 \; X_2]\) um vetor aleatório com vetor médio \(\boldsymbol{\mu}_{\mathbf{X}}^{\prime} = [\mu_1 \; \mu_2]\) e matriz de covariância
\[ \mathbf{\Sigma_X} = \begin{bmatrix} \sigma_1^2 & \sigma_{12} \\ \sigma_{12} & \sigma_2^2 \end{bmatrix} \]
Encontre o vetor médio e a matriz de covariância para as combinações lineares
\[ Z_1 = X_1 - X_2\\ Z_2 = X_1 + X_2 \]
ou
\[ \mathbf{Z}=\mathbf{CX} \]
Então:
\[ \mathbf{\boldsymbol{\mu}_{\mathbf{Z}}}=\mathbb{E}(\mathbf{Z})=\mathbf{C\boldsymbol{\mu}_{\mathbf{X}}}= \begin{bmatrix} \mu_1-\mu_2\\ \mu_1+\mu_2 \end{bmatrix} \]
\[ \begin{align} \mathbf{\Sigma_Z} &=\mathbb{C}(\mathbf{Z}) \\ &=\mathbf{C \Sigma_X C}^{\prime} \\ &=\begin{bmatrix} \mathbb{C}(Z_1, Z_1) & \mathbb{C}(Z_1, Z_2) \\ \mathbb{C}(Z_2, Z_1) & \mathbb{C}(Z_2, Z_2) \end{bmatrix} \\ \mathbf{\Sigma_Z} &= \begin{bmatrix} \sigma_1^2 - 2\sigma_{12} + \sigma_2^2 & \sigma_1^2-\sigma_2^2 \\ \sigma_1^2-\sigma_2^2 & \sigma_1^2 + 2\sigma_{12} + \sigma_2^2 \end{bmatrix} \end{align} \]
Note que se \(\sigma_1^2 = \sigma_2^2\), ou seja, se \(X_1\) e \(X_2\) tiverem variâncias populacionais iguais, os termos fora da diagonal se anulam. Isso demonstra o resultado conhecido de que a soma e a diferença de duas variáveis aleatórias com variâncias idênticas (homocedasticidade) são não correlacionadas.
Se \(\sigma_{11}=\sigma_{22}\) e \(\sigma_{12}=\rho\sigma_1\sigma_2\), então:
\[ \begin{align} \mathbf{\Sigma_Z} &= 2\sigma_1^2\begin{bmatrix} 1- \rho & 0 \\ 0 & 1 + \rho \end{bmatrix} \end{align} \]
Princípios de maximização desempenham um papel importante em várias técnicas multivariadas. Por exemplo, a análise discriminante linear está relacionada à alocação de observações a grupos predefinidos. A regra de alocação é frequentemente uma função linear das medições que maximiza a separação entre os grupos em relação à variabilidade dentro de cada grupo. Outro exemplo são as componentes principais, que são combinações lineares de medições com a máxima variabilidade.
Seja \(\mathbf{\Sigma}\) uma matriz \(p \times p\) simétrica e definida positiva, i.e., \(\lambda_1\ge \lambda_2\ge \cdots\ge \lambda_p >0\), e seus respectivos autovetores ortonormalizados \(\mathbf{e}_1, \mathbf{e}_2, \ldots, \mathbf{e}_p\). Então,
\[ \max_{\mathbf{x}} \dfrac{\mathbf{x}^{\prime} \mathbf{\Sigma} \mathbf{x}}{\mathbf{x}^{\prime} \mathbf{x}} = \lambda_1 \quad \text{(alcançado se} \; \mathbf{x}^{*} = \mathbf{e}_1) \tag{2-51} \]
\[ \min_{\mathbf{x}} \dfrac{\mathbf{x}^{\prime} \mathbf{\Sigma} \mathbf{x}}{\mathbf{x}^{\prime} \mathbf{x}} = \lambda_p \quad \text{(alcançado se} \; \mathbf{x}^{*} = \mathbf{e}_p) \tag{2-52} \]
em que \(\mathbf{e}_1\) é o \(1\)-ésimo autovetor associado ao autovalor \(\lambda_1\), e \(\mathbf{e}_p\) é o \(p\)-ésimo autovetor associado ao menor autovalor \(\lambda_p\).
2.1, 2.6, 2.7, 2.9, 2.10, 2.24, 2.25, 2.26, 2.41, 2.42
Com os conceitos vetoriais introduzidos no capítulo anterior, agora podemos explorar mais profundamente as interpretações geométricas das estatísticas descritivas \(\mathbf{\bar{x}}\), \(\mathbf{s}_n\), e \(\mathbf{r}\); fazemos isso na Seção 3.2. Muitas das nossas explicações utilizam a representação das colunas de \(\boldsymbol{\mathcal{x}}\) como \(p\) vetores em \(n\) dimensões. Na Seção 3.3, introduzimos a suposição de que as observações constituem uma amostra aleatória. Em termos simples, amostragem aleatória implica que (1) medições feitas em diferentes itens (ou ensaios) não têm relação entre si e (2) a distribuição conjunta de todas as \(p\) variáveis permanece a mesma para todos os itens. Em última análise, é essa estrutura da amostra aleatória que justifica uma escolha específica de distância e determina a geometria para a representação \(z_i\) dimensional dos dados. Além disso, quando os dados podem ser tratados como uma amostra aleatória, as inferências estatísticas são baseadas em uma base sólida.
Retornando às interpretações geométricas na Seção 3.4, introduzimos um único número, chamado variância generalizada, para descrever a variabilidade. Essa generalização da variância é parte integrante da comparação de médias multivariadas. Em seções posteriores, utilizamos álgebra matricial para fornecer expressões concisas para os produtos e somas de matrizes que nos permitem calcular \(\mathbf{\bar{x}}\) e \(\mathbf{s}_n\) diretamente a partir da matriz de dados \(\boldsymbol{\mathcal{x}}\). A conexão entre \(\mathbf{\bar{x}}\), \(\mathbf{s}_n\), e as médias e covariâncias para combinações lineares de variáveis também é claramente delineada, usando a noção de produtos de matrizes.
Uma única observação multivariada é a coleção de medições em \(p\) diferentes variáveis, obtidas no mesmo item ou ensaio. Como no Capítulo 1, se foram obtidas \(n\) observações, todo o conjunto de dados pode ser colocado em uma matriz \(n \times p\):
\[ \boldsymbol{\mathcal{x}} = \begin{bmatrix} x_{11} & x_{12} & \ldots & x_{1p} \\ x_{21} & x_{22} & \ldots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \ldots & x_{np} \end{bmatrix} \]
Cada linha de \(\boldsymbol{\mathcal{x}}\) representa uma observação multivariada. Como o conjunto inteiro de medições frequentemente é apenas uma realização particular do que poderia ter sido observado, dizemos que os dados são uma amostra de tamanho \(n\) de uma “população” \(p\)-variada. A amostra consiste, portanto, em \(n\) medições, cada uma delas com \(p\) componentes.
Conforme já vimos, os dados podem ser plotados de duas maneiras diferentes. No gráfico de dispersão \(p\)-dimensional, as linhas de \(\boldsymbol{\mathcal{x}}\) representam \(n\) pontos em um espaço \(p\)-dimensional. Podemos escrever:
\[ \underset{n\times p}{\boldsymbol{\mathcal{x}}} = \begin{bmatrix} x_{11} & x_{12} & \ldots & x_{1p} \\ x_{21} & x_{22} & \ldots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \ldots & x_{np} \end{bmatrix} = \begin{bmatrix} \mathbf{x}_{1}^{\prime}\\ \mathbf{x}_{2}^{\prime}\\ \vdots\\ \mathbf{x}_{n}^{\prime}\\ \end{bmatrix} \tag{3-1} \]
O vetor-linha \(\mathbf{x}_i\), representando a \(i\)-ésima observação, contém as coordenadas de um ponto.
O gráfico de dispersão dos \(n\) pontos em um espaço \(p\)-dimensional fornece informações sobre as localizações e variabilidade dos pontos. Se considerarmos os pontos como esferas sólidas, o vetor de média amostral \(\mathbf{\bar{x}}\), dado por (1-8), é o centro de massa. A variabilidade ocorre em mais de uma direção, e ela é quantificada pela matriz de variância-covariância amostral \(\mathbf{s}_n\). Uma única medida numérica de variabilidade é fornecida pelo determinante da matriz de variância-covariância amostral. Quando \(p\) é maior que 3, essa representação de gráfico de dispersão não pode ser efetivamente representada graficamente. No entanto, a consideração dos dados como \(n\) pontos em \(p\) dimensões fornece insights que não estão prontamente disponíveis em expressões algébricas. Além disso, os conceitos ilustrados para \(p = 2\) ou \(p = 3\) permanecem válidos para os outros casos.
Cada linha \(\mathbf{x}_i\), representando a \(i\)-ésima observação, contém as coordenadas de um ponto no espaço \(p\)-dimensional.
Para calcular o vetor médio \(\mathbf{\bar{x}}\) a partir da matriz de dados, dada por
\[ \boldsymbol{\mathcal{x}} = \begin{bmatrix} 4 & 1 \\ -1 & 3 \\ 3 & 5 \end{bmatrix} \]
primeiro somamos as colunas da matriz para obter a soma das coordenadas em cada dimensão:
\[ \text{Soma das coordenadas em dimensão 1} = 4 + (-1) + 3 = 6 \]
\[ \text{Soma das coordenadas em dimensão 2} = 1 + 3 + 5 = 9 \]
Agora, dividimos essas somas pelo número de observações, que neste caso é 3, para obter as coordenadas do vetor médio \(\mathbf{x}\):
\[ \mathbf{\bar{x}} = \begin{bmatrix} \dfrac{6}{3} \\ \dfrac{9}{3} \end{bmatrix} = \begin{bmatrix} 2 \\ 3 \end{bmatrix} \]
Agora, plotamos os 3 pontos de dados \(\mathbf{x}_i\) em um espaço \(p = 2\) e localizamos o vetor médio \(\bar{\mathbf{x}}\) no diagrama resultante. Os pontos de dados são:
\[ \begin{align} \mathbf{x}_1^{\prime} &= [4, 1]\\ \mathbf{x}_2^{\prime} &= [-1, 3]\\ \mathbf{x}_3^{\prime} &= [3, 5] \end{align} \]
O vetor de média amostral \(\mathbf{\bar{x}}^{\prime}\) é [2, 3].
Figura 3.1 Um gráfico da matriz de dados x com n = 3 pontos em espaço bidimensional p = 2.
X <- matrix(c(4, 1, -1, 3, 3, 5), ncol = 2, byrow = TRUE)
x_coords <- X[, 1]
y_coords <- X[, 2]
plot(x_coords, y_coords, pch = 21, col = "black",
xlab = "x1", ylab = "x2", main="Vetor de média amostral")
mean_x <- mean(x_coords)
mean_y <- mean(y_coords)
points(mean_x, mean_y, pch = 16, col = "black")
\[\Diamond\]
A representação geométrica alternativa é construída ao considerarmos os dados como \(p\) vetores em um espaço \(n\)-dimensional. Aqui, tomamos os elementos das colunas da matriz de dados como as coordenadas dos vetores. Seja \(\mathbf{x}\):
\[ \underset{n \times p}{\boldsymbol{\mathcal{x}}} = \begin{bmatrix} x_{11} & x_{12} & \ldots & x_{1p} \\ x_{21} & x_{22} & \ldots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \ldots & x_{np} \end{bmatrix} = \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 & \cdots & \mathbf{y}_p \end{bmatrix} \tag{3-2} \]
Então, as coordenadas do primeiro ponto \(\mathbf{y}_1^{\prime} = [x_{11}\; x_{21}\; \ldots\; x_{n1}]\) são as \(n\) medições na primeira variável. Em geral, o \(j\)-ésimo ponto \(\mathbf{y}_i^{\prime} = [x_{1i}\; x_{2i}\; \ldots\; x_{ni}]\) é determinado pelo n-ple de todas as medições na \(i\)-ésima variável. Nesta representação geométrica, retratamos \(\mathbf{y}_1, \mathbf{y}_2,\ldots, \mathbf{y}_p\) como vetores no espaço \(n\)-dimensional.
Vamos plotar os seguintes dados como \(p = 2\) vetores em espaço \(n = 3\) dimensões:
\[ \boldsymbol{\mathcal{x}} = \begin{bmatrix} 4 & 1 \\ -1 & 3 \\ 3 & 5 \end{bmatrix} \]
Figura 3.2 Um gráfico da matriz de dados x com p = 2 vetores em um espaço n = 3.
Aqui, \(\mathbf{y}_1^{\prime} = [4\; -1 \; 3]\) e \(\mathbf{y}_2^{\prime} = [1\; 3\; 5]\). Esses vetores são mostrados na Figura 3.2.
\[\Diamond\]
Muitas das expressões algébricas que encontraremos na análise multivariada podem ser relacionadas com noções geométricas de comprimento, ângulo e volume. Isso é importante porque representações geométricas normalmente facilitam a compreensão e levam a insights adicionais.
Infelizmente, estamos limitados a visualizar objetos em três dimensões e, consequentemente, a representação \(p\)-dimensional da matriz de dados \(\boldsymbol{\mathcal{x}}\) pode não parecer um dispositivo especialmente útil para \(n > 3\). No entanto, acontece que as relações geométricas e os conceitos estatísticos associados representados por qualquer conjunto de três vetores permanecem válidos, independentemente de sua dimensão. Isso ocorre porque três vetores, mesmo que tenham \(n\) dimensões, podem abranger no máximo um espaço tridimensional, assim como dois vetores com qualquer número de componentes devem estar contidos em um plano. Selecionando uma perspectiva tridimensional apropriada - ou seja, uma porção do espaço \(n\)-dimensional contendo os três vetores de interesse - obtém-se uma visão que preserva tanto os comprimentos quanto os ângulos. Assim, é possível, com a escolha adequada dos eixos, ilustrar certos conceitos estatísticos algébricos em termos de apenas dois ou três vetores de qualquer dimensão \(n\). Uma vez que a escolha específica dos eixos não é relevante para a geometria, sempre rotularemos os eixos de coordenadas como 1, 2 e 3.
É possível dar uma interpretação geométrica ao processo de encontrar uma média amostral. Começamos definindo o vetor \(n \times 1\) como \(\mathbf{1}_{n}^{\prime} = [1\; 1\; \ldots\; 1]\). (Para simplificar a notação, o subscrito \(n\) será omitido quando a dimensão do vetor estiver clara no contexto.) O vetor \(\mathbf{1}_{n}\) forma ângulos iguais com cada um dos eixos de coordenadas, de modo que o vetor \((1/\sqrt{n})\mathbf{1}\) tem comprimento unitário na direção de ângulo igual. Considere o vetor \(\mathbf{y}_i^{\prime} = [x_{1i}, x_{2i}, \ldots, x_{ni}]\). A projeção de \(\mathbf{y}_i\) no vetor unitário \((1/\sqrt{n})\mathbf{1}\) é dada por (2-8):
\[ \text{Projeção de } \mathbf{y}_i \text{ em } \dfrac{1}{\sqrt{n}}\mathbf{1} = \bar{x}_i \mathbf{1} \]
Ou seja, a média amostral
\[ \bar{x}_i = \frac{x_{1i} + x_{2i} + \ldots + x_{ni}}{n} = \frac{\mathbf{y}_i^{\prime} \mathbf{1}}{n} \]
corresponde ao múltiplo de \(\mathbf{1}\) necessário para dar a projeção de \(\mathbf{y}_i\) na reta determinada por \(\mathbf{1}\).
Além disso, para cada \(\mathbf{y}_i\), temos a decomposição:
em que \(\bar{x}_i \mathbf{1}\) é perpendicular a \(\mathbf{y}_i - \bar{x}_i \mathbf{1}\). O vetor de desvio, ou vetor corrigido pela média, é dado por:
\[ \mathbf{d}_i = \mathbf{y}_i - \bar{x}_i \mathbf{1} \tag{3-4} \]
Os elementos de \(\mathbf{d}_i\) são os desvios das medições da \(i\)-ésima variável em relação à sua média amostral. A decomposição dos vetores \(\mathbf{y}_i\) em componentes de média e desvio em relação à média é mostrada na Figura 3.3 para \(p = 3\) e \(n = 3\).
Vamos realizar a decomposição dos vetores \(\mathbf{y}_i\) em suas componentes de média e desvio para os dados fornecidos no Exemplo 3.2:
Dados:
\[ \boldsymbol{\mathcal{x}} = \begin{bmatrix} 4 & 1 \\ -1 & 3 \\ 3 & 5 \end{bmatrix} \]
Calcular a média amostral \(\bar{\mathbf{x}}\) para cada variável.
\[ \bar{\mathbf{x}}^{\prime} = \begin{bmatrix} \dfrac{4+(-1)+3}{3} & \dfrac{1+3+5}{3} \end{bmatrix} = \begin{bmatrix} 2 & 3 \end{bmatrix} \]
\[ \bar{x}_1 \mathbf{1} = 2 \begin{bmatrix} 1\\1\\1 \end{bmatrix} = \begin{bmatrix} 2\\2\\2 \end{bmatrix} \]
\[ \bar{x}_2 \mathbf{1} = 3 \begin{bmatrix} 1\\1\\1 \end{bmatrix} = \begin{bmatrix} 3\\3\\3 \end{bmatrix} \]
\[ \mathbf{d}_1=\mathbf{y}_1-\bar{x}_1 \mathbf{1} = \begin{bmatrix} 4\\-1\\3 \end{bmatrix} - \begin{bmatrix} 2\\2\\2 \end{bmatrix} = \begin{bmatrix} 2\\-3\\1 \end{bmatrix} \]
\[ \mathbf{d}_2=\mathbf{y}_2-\bar{x}_2 \mathbf{1} = \begin{bmatrix} 1\\3\\5 \end{bmatrix} - \begin{bmatrix} 3\\3\\3 \end{bmatrix} = \begin{bmatrix} -2\\0\\2 \end{bmatrix} \]
Observamos que \(\mathbf{x}_{i} \mathbf{1}\) e \(\mathbf{d}_i = \mathbf{y}_i - \bar{x}_1 \mathbf{1}\) são perpendiculares, pois
\[ \left( \mathbf{x}_{1} \mathbf{1} \right)^{\prime} \mathbf{d}_1 = \begin{bmatrix} 2 & 2 & 2\end{bmatrix} \begin{bmatrix} 2 \\ -3 \\ 1 \end{bmatrix} = 4 - 6 +2 = 0 \]
\[ \left( \mathbf{x}_{2} \mathbf{1} \right)^{\prime} \mathbf{d}_2 = 0 \]
Essa propriedade decorre da definição de projeção ortogonal, em que a projeção do vetor \(\mathbf{y}_i\) no subespaço gerado pelo vetor \(\mathbf{x}_{i} \mathbf{1}\) é perpendicular ao vetor de desvio \(\mathbf{d}_i\). O produto escalar entre vetores perpendiculares é igual a zero, o que confirma a perpendicularidade entre \(\mathbf{x}_{i}\) e \(\mathbf{d}_i\).
Traduzimos os vetores de desvio para a origem sem alterar seus comprimentos ou orientações. Agora, consideremos os comprimentos quadrados dos vetores de desvio. Usando as equações (2-5) e (3-4), obtemos
\[ \lVert \mathbf{d}_i \rVert^2 = \mathbf{d}_i^{\prime} \mathbf{d}_i = \sum_{j=1}^{n}{(x_{ji} - \bar{x}_i)^2} \tag{3-5} \]
(Comprimento do vetor de desvio)² = soma dos quadrados dos desvios
A partir da equação (1-3), observamos que o comprimento ao quadrado é proporcional à variância das medições na \(i\)-ésima variável. Equivalentemente, o comprimento é proporcional ao desvio-padrão. Vetores mais longos representam mais variabilidade do que vetores mais curtos.
Vamos denotar por \(\theta_{ik}\) o ângulo formado pelos vetores \(\mathbf{d}_i\) e \(\mathbf{d}_k\). A partir da equação (2-6), obtemos
\[ \cos(\theta_{ik}) = \dfrac{\mathbf{d}_i^{\prime} \mathbf{d}_k}{\lVert \mathbf{d}_i \rVert \lVert \mathbf{d}_k \rVert} \]
ou, usando (3-5) e (3-6), obtemos
\[ \cos(\theta_{ik}) = \dfrac{\sum_{j=1}^{n} (x_{ji} - \bar{x}_i)(x_{jk} - \bar{x}_k)}{\sqrt{\sum_{j=1}^{n} (x_{ji} - \bar{x}_i)^2} \sqrt{\sum_{j=1}^{n} (x_{jk} - \bar{x}_k)^2}} \]
Então, de acordo com (1-5),
\[ r_{ik} = \dfrac{s_{ik}}{\sqrt{s_{ii} s_{kk}}} = \cos(\theta_{ik}) \tag{3-7} \]
O cosseno do ângulo é o coeficiente de correlação de Pearson amostral. Assim, se os dois vetores de desvio tiverem orientações quase iguais, a correlação amostral será próxima de 1. Se os dois vetores forem quase perpendiculares, a correlação amostral será aproximadamente zero. Se os dois vetores estiverem orientados em direções quase opostas, a correlação amostral será próxima de -1.
Dado o vetor de desvio do Exemplo 3.3, vamos calcular a matriz de variância e covariância amostral \(\mathbf{s}_n\) e a matriz de correlação amostral \(\mathbf{r}\) utilizando os conceitos geométricos apresentados anteriormente.
Do Exemplo 3.3, temos:
\[ \begin{align} \mathbf{d}_1 &= \begin{bmatrix} 2 \\ -3 \\ 1 \end{bmatrix}\\ \mathbf{d}_2 &= \begin{bmatrix} -2 \\ 0 \\2\end{bmatrix} \end{align} \]
Passo 1: Calcular a matriz de variância e covariância amostral \(\mathbf{s}_n\).
A variância amostral da primeira variável (variância de \(\mathbf{x}_1\)) é dada pelo comprimento quadrado do vetor de desvio \(\mathbf{d}_1\):
\[ \begin{align} 3s_{11} &= \lVert \mathbf{d}_1 \rVert^2 = 2^2 + (-3)^2 + 1^2= 14 \\ s_{11} &= \dfrac{14}{3} \end{align} \]
A variância amostral da segunda variável (variância de \(x_2\)) é dada pelo comprimento quadrado do vetor de desvio \(\mathbf{d}_2\):
\[ \begin{align} 3s_{22} &= \lVert \mathbf{d}_2 \rVert^2 = (-2)^2 + 0^2 +2^2= 8\\ s_{22}&=\dfrac{8}{3} \end{align} \]
A covariância amostral entre as duas variáveis (covariância entre \(x_1\) e \(x_2\)) é dada pelo produto interno entre os vetores de desvio \(\mathbf{d}_1\) e \(\mathbf{d}_2\):
\[ \begin{align} 3s_{12} &= \mathbf{d}_1^{\prime} \mathbf{d}_2 = 2\times (-2) + (-3)\times 0+2\times 1 = -2\\ s_{12}&=-\dfrac{2}{3} \end{align} \]
Portanto, a matriz de variância e covariância amostral \(\mathbf{s}_n\) é:
\[ \mathbf{s}_n = \begin{bmatrix} s_{11} & s_{12} \\ s_{12} & s_{22} \end{bmatrix} = \begin{bmatrix} \dfrac{14}{3} & -\dfrac{2}{3} \\ -\dfrac{2}{3} & \dfrac{8}{3} \end{bmatrix} \]
Passo 2: Calcular a matriz de correlação amostral \(\mathbf{r}\).
A matriz de correlação amostral \(\mathbf{r}\) é obtida dividindo cada elemento da matriz de variância e covariância \(\mathbf{s}_n\) pelo produto da raiz quadrada das variâncias das variáveis correspondentes:
\[ r_{12} = \frac{s_{12}}{\sqrt{s_{11} \ s_{22}}}= -0.189 \]
Portanto, a matriz de correlação amostral \(\mathbf{R}\) é:
\[ \mathbf{r} = \begin{bmatrix} r_{11} & r_{12} \\ r_{12} & r_{22} \end{bmatrix} = \begin{bmatrix} 1 & -0.189 \\ -0.189 & 1 \end{bmatrix} \]
Essas são as matrizes \(\mathbf{s}_n\) e \(\mathbf{r}\) calculadas a partir dos vetores de desvio \(\mathbf{d}_1\) e \(\mathbf{d}_2\).
Vetores de desvio d1 e d1
\[\Diamond\]
Os conceitos de comprimento, ângulo e projeção nos forneceram uma interpretação geométrica da amostra. Resumimos da seguinte forma:
A projeção de uma coluna \(\mathbf{y}_i\) da matriz de dados \(\boldsymbol{\mathcal{x}}\) no vetor de ângulo igual \(\mathbf{1}\) é o vetor \(\bar{x}_i \mathbf{1}\). O vetor \(\bar{x}_i \mathbf{1}\) tem comprimento \(\sqrt{n}|\bar{x}_i|\). Portanto, a \(i\)-ésima média amostral \(\bar{\mathbf{x}_i}\) está relacionada ao comprimento da projeção de \(\mathbf{y}_i\) no vetor \(\mathbf{1}\).
As informações que compõem \(\mathbf{S}_n\) são obtidas dos vetores de desvio \(\mathbf{d}_i = \mathbf{y}_i -\bar{x}_i \mathbf{1} = [x_{1i} - \bar{x}_i, x_{2i} - \bar{x}_i, \ldots, x_{ni} - \bar{x}_i]\). O quadrado do comprimento de \(\mathbf{d}_i\) é \(ns_{ii}\), e o produto (interno) entre \(\mathbf{d}_i\) e \(\mathbf{d}_k\) é \(ns_{ik}\).
A correlação amostral \(r_{ik}\) é o cosseno do ângulo entre \(\mathbf{d}_i\) e \(\mathbf{d}_k\).
Essa interpretação geométrica nos ajuda a visualizar e entender as características da amostra, como as relações entre as variáveis e a variabilidade dos dados. Os conceitos de comprimento, ângulo e projeção proporcionam uma perspectiva valiosa para análise e interpretação dos dados multivariados.
Para estudar a variabilidade amostral de estatísticas como \(\bar{\mathbf{x}}\) e \(\mathbf{s}_n\) com o objetivo final de fazer inferências, precisamos fazer suposições sobre as variáveis cujos valores observados constituem o conjunto de dados \(\boldsymbol{\mathcal{x}}\).
Suponha, então, que os dados ainda não foram observados, mas pretendemos coletar \(n\) conjuntos de medições em \(p\) variáveis. Antes que as medições sejam feitas, seus valores não podem, em geral, ser previstos exatamente. Consequentemente, tratamos os valores como variáveis aleatórias. Nesse contexto, o \((i, j)\)-ésimo elemento da matriz de dados é a variável aleatória \(X_{ij}\). Cada conjunto de medições \(\mathbf{X}_j\) em \(p\) variáveis é um vetor aleatório, e temos a matriz aleatória
\[ \boldsymbol{\mathcal{X}} = \begin{bmatrix} X_{11} & X_{12} & \dots & X_{1p} \\ X_{21} & X_{22} & \dots & X_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ X_{n1} & X_{n2} & \dots & X_{np} \end{bmatrix} = \begin{bmatrix} \mathbf{X}_{1}^{\prime}\\\mathbf{X}_{2}^{\prime}\\\vdots\\\mathbf{X}_{n}^{\prime}\\ \end{bmatrix} \tag{3-8} \]
Uma amostra aleatória pode ser definida agora.
Se os vetores linha \(\mathbf{X}_1, \mathbf{X}_2, \dots, \mathbf{X}_n\) em (3-8) representam observações independentes de uma distribuição conjunta com função de densidade \(f(\mathbf{x}) = f(x_1, x_2, \dots, x_p)\), então \(\mathbf{X}_1, \mathbf{X}_2, \dots, \mathbf{X}_n\) são considerados uma amostra aleatória de \(f(\mathbf{x})\). Matematicamente, \(\mathbf{X}_1, \mathbf{X}_2, \dots, \mathbf{X}_n\) formam uma amostra aleatória se sua função de densidade conjunta for dada pelo produto \(f(\mathbf{x}_1)f(\mathbf{x}_2) \cdots f(\mathbf{x}_n)\), onde \(f(\mathbf{x}_j) = f(x_{j1}, x_{j2}, \dots, x_{jp})\) é a função de densidade para o \(j\)-ésimo vetor linha.
Dois pontos relacionados à definição de amostra aleatória merecem atenção especial:
As medições das \(p\) variáveis em um único experimento, como \(\mathbf{X}_j = [X_{j1}, X_{j2}, \dots, X_{jp}]\), geralmente estarão correlacionadas. De fato, esperamos que isso seja o caso. As medições de diferentes experimentos, no entanto, devem ser independentes.
A independência das medições de um experimento para outro pode não ser válida quando as variáveis têm probabilidade de variar ao longo do tempo, como no caso de conjuntos de \(p\) preços de ações ou \(p\) indicadores econômicos. Violações da suposição de independência podem ter um impacto sério na qualidade das inferências estatísticas.
Os exemplos a seguir ilustram essas observações.
Como etapa preliminar para projetar um sistema de permissão para utilização de uma área de canoagem em uma região selvagem sem superlotação, um gerente de recursos naturais realizou uma pesquisa com os usuários. A área total da região selvagem foi dividida em sub-regiões, e os entrevistados foram solicitados a fornecer informações sobre as regiões visitadas, duração da estadia e outras variáveis.
O método seguido foi selecionar aleatoriamente pessoas (possivelmente usando uma tabela de números aleatórios) de todos aqueles que entraram na área de canoagem durante uma semana específica. Todas as pessoas tinham igual probabilidade de serem selecionadas para a amostra, portanto, as entradas mais populares foram representadas por uma proporção maior de canoístas.
Nesse caso, esperamos que as observações da amostra se aproximem de perto do critério para uma amostra aleatória da população de usuários ou potenciais usuários. Por outro lado, se um dos entrevistadores esperasse em um acampamento bem no interior da área e entrevistasse apenas os canoístas que chegassem àquele local, as medições sucessivas não seriam independentes. Por exemplo, a duração da estadia na área selvagem para diferentes canoístas desse grupo tenderia a ser grande.
\[\Diamond\]
Devido às preocupações com o futuro do descarte de resíduos sólidos, um estudo em andamento se refere ao peso bruto de resíduos sólidos municipais gerados por ano nos Estados Unidos (Agência de Proteção Ambiental). As quantidades estimadas atribuídas a \(x_1 =\) resíduos de papel e papelão e \(x_2 =\) resíduos plásticos, em milhões de toneladas, são apresentadas para anos selecionados na Tabela 3.1. Essas medições em \(\mathbf{x}^{\prime} = [x_1, x_2]\) devem ser tratadas como uma amostra aleatória de tamanho \(n = 7\)? Não! De fato, exceto por uma ligeira, mas afortunada, queda nos resíduos de papel e papelão em 2003, ambas as variáveis estão aumentando ao longo do tempo.
Tabela 3.1 - Resíduos Sólidos
Ano | Papel | Plástico |
---|---|---|
1960 | 29.2 | 0.4 |
1970 | 44.3 | 2.9 |
1980 | 55.2 | 6.8 |
1990 | 72.7 | 17.1 |
1995 | 81.7 | 18.9 |
2000 | 87.7 | 24.7 |
2003 | 83.1 | 26.7 |
# Dados da tabela
ano <- c(1960, 1970, 1980, 1990, 1995, 2000, 2003)
papel <- c(29.2, 44.3, 55.2, 72.7, 81.7, 87.7, 83.1)
plasticos <- c(0.4, 2.9, 6.8, 17.1, 18.9, 24.7, 26.7)
# Criar a tabela em formato de data frame
dados <- data.frame(Ano = ano, Papel = papel, Plasticos = plasticos)
knitr::kable(dados, caption="Resíduos Sólidos")
Ano | Papel | Plasticos |
---|---|---|
1960 | 29.2 | 0.4 |
1970 | 44.3 | 2.9 |
1980 | 55.2 | 6.8 |
1990 | 72.7 | 17.1 |
1995 | 81.7 | 18.9 |
2000 | 87.7 | 24.7 |
2003 | 83.1 | 26.7 |
\[\Diamond\]
Como argumentamos heuristicamente no Capítulo 1, a noção de independência estatística tem implicações importantes para a medição da distância. A distância euclidiana parece adequada se os componentes de um vetor forem independentes e tiverem a mesma variância. Suponha que consideramos a localização da \(k\)-ésima coluna \(\mathbf{Y}_k^{\prime} = [X_{1k}, X_{2k}, \ldots, X_{nk}]\) de \(\boldsymbol{\mathcal{X}}\), considerada como um ponto em \(n\) dimensões. A localização deste ponto é determinada pela distribuição de probabilidade conjunta \(f(\mathbf{y}_k)\). Quando as medições \(X_{1k}, X_{2k}, \ldots, X_{nk}\) são uma amostra aleatória,
\[ \begin{align} f(\mathbf{y}_k) &= f(x_{1k}, x_{2k}, \ldots, x_{nk}) \\ &= f_k(x_{1k})f_k(x_{2k})\cdots f_k(x_{nk})\\ f(\mathbf{y}_k) &= \prod_{i=1}^{n}{f_k(x_{ik})} \end{align} \]
e, consequentemente, cada coordenada \(x_{jk}\) contribui igualmente para a localização através das distribuições marginais idênticas \(f_k(x_{jk})\).
Se os \(n\) componentes não são independentes ou as distribuições marginais não são idênticas, a influência das medições individuais (coordenadas) na localização é assimétrica. Nesse caso, seríamos levados a considerar uma função de distância na qual as coordenadas são ponderadas de forma desigual, como nas distâncias “estatísticas” ou formas quadráticas introduzidas nos Capítulos 1 e 2.
Determinadas conclusões podem ser alcançadas a respeito das distribuições amostrais de \(\bar{\mathbf{X}}\) e \(\mathbf{S}_n\) sem fazer mais suposições sobre a forma da distribuição conjunta subjacente das variáveis. Em particular, podemos observar como \(\bar{\mathbf{X}}\) e \(\mathbf{S}_n\) se comportam como estimadores pontuais do correspondente vetor de média populacional \(\boldsymbol{\mu}\) e matriz de covariância \(\mathbf{\Sigma}\).
Resultado 3.1. Se \(\mathbf{X} _1, \mathbf{X} _2, \ldots, \mathbf{X} _n\) formarem uma amostra aleatória de uma distribuição conjunta com vetor médio \(\boldsymbol{\mu}\) e matriz de covariância \(\mathbf{\Sigma}\), então \(\bar{\mathbf{X}}\) é um estimador não viesado de \(\boldsymbol{\mu}\) (vetor de média populacional) e \(\dfrac{n}{n-1}\mathbf{S}_n\) é um estimador não viesado de \(\mathbf{\Sigma}\) (matriz de covariância populacional).
Formalmente, a suposição é: As \(n\) observações multivariadas são independentes e identicamente distribuídas (IID).
\[ \left\{\mathbf{X}_i \right\}_{i=1}^{n} \sim \text{IID}(\boldsymbol{\mu}, \mathbf{\Sigma}) \]
Para a matriz de covariância \(\mathbf{S}_n\), \(\mathbb{E}(\mathbf{S}_n) = \dfrac{n-1}{n} \mathbf{\Sigma}\).
Portanto, \(\dfrac{n}{n-1} \mathbf{S}_n\) é um estimador não viesado de \(\mathbf{\Sigma}\), enquanto \(\mathbf{S}_n\) é um estimador viesado com viés igual a \(\mathbb{E}(\mathbf{S}_n) - \mathbf{\Sigma} =-\dfrac{1}{n-1}\mathbf{\Sigma}\).
Prova:
\(\bar{\mathbf{X}} = \dfrac{1}{n}\sum_{i=1}^{n}{\mathbf{X}_i}\). O uso repetido das propriedades da expectativa em (2-24) para dois vetores nos resulta:
\[ \begin{align} \mathbb{E}\left(\bar{\mathbf{X}}\right) &= \mathbb{E}\left(\dfrac{1}{n}\sum_{i=1}^{n}{\mathbf{X}_i}\right) \\ &= \dfrac{1}{n}\sum_{i=1}^{n}{\mathbb{E}\left(\mathbf{X}_i\right)} \\ &= \dfrac{1}{n}\sum_{i=1}^{n}{\boldsymbol{\mu}}\\ \mathbb{E}\left(\bar{\mathbf{X}}\right) &= \boldsymbol{\mu} \end{align} \]
Portanto,
\[ \mathbb{E}\left(\mathbf{S}_n\right) = \dfrac{n-1}{n} \mathbf{\Sigma} \]
Logo, \(\dfrac{n}{n-1} \mathbf{S}_n\) é um estimador não tendencioso de \(\mathbf{\Sigma}\), enquanto \(\mathbf{S}_n\) é um estimador tendencioso com viés igual a \(\mathbb{E}(\mathbf{S}_n) - \mathbf{\Sigma}\).
A consideração do viés motiva uma definição ligeiramente modificada da matriz de covariância amostral.
Pelo Resultado 3.1, tem-se que o estimador \(\mathbf{S}\) é não viesado:
\[ \begin{align} \mathbb{E}(\mathbf{S})&=\mathbf{\Sigma} \\ \mathbf{S} &= \dfrac{n}{n-1} \mathbf{S}_n = \dfrac{1}{n-1} \sum_{j=1}^{n} (\mathbf{X}_i - \bar{\mathbf{X}})(\mathbf{X}_i - \bar{\mathbf{X}})^{\prime} \end{align} \]
Essa definição de matriz de covariância amostral é comumente usada em muitas estatísticas de teste multivariadas. Portanto, ela substitui \(\mathbf{S}_n\) por \(\mathbf{S}\) como a matriz de covariância amostral na maior parte do material ao longo do restante deste livro.
Com uma única variável, a variância amostral é frequentemente usada para descrever a quantidade de variação nas medições dessa variável. Quando \(p\) variáveis são observadas em cada unidade, a variação é descrita pela matriz de covariância amostral:
\[ \underset{p\times p}{\mathbf{s}} = \begin{bmatrix} s_{11} & s_{12} & \cdots & s_{1p} \\ s_{21} & s_{22} & \cdots & s_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ s_{p1} & s_{p2} & \cdots & s_{pp} \end{bmatrix} \]
em que \(s_{ij}\) é a covariância amostral entre as variáveis \(i\) e \(j\), \(1 \leq i,j \leq p\), e \(n\) é o tamanho da amostra.
A matriz de covariância amostral contém \(p\) variâncias e \(p(p - 1)\) covariâncias potencialmente diferentes. Às vezes, é desejável atribuir um único valor numérico para a variação expressa por \(\mathbf{s}\). Uma escolha para esse valor é o determinante de \(\mathbf{s}\), que se reduz à variância amostral usual de uma única característica quando \(p = 1\). Esse determinante é chamado de variância generalizada amostral:
\[ \text{Variância generalizada amostral} = |\mathbf{s}| \tag{3-12} \]
O número de funcionários (\(x_1\)) e lucro por funcionário (\(x_2\)) das 16 maiores empresas dos Estados Unidos estão mostrados na Figura 1.3. A matriz de covariância amostral, obtida a partir dos dados da revista Forbes de 30 de abril de 1990, é:
\[ \mathbf{s} = \begin{bmatrix} 252.04 & -68.43 \\ -68.43 & 123.67 \end{bmatrix} \]
Vamos calcular a variância generalizada.
Nesse caso, podemos calcular \(|\mathbf{s}| = 252.04\times 123.67 - (-68.43\times -68.43) = 26487\).
A variância amostral generalizada fornece uma maneira de resumir as informações de todas as variâncias e covariâncias em um único número. No entanto, quando \(p > 1\), parte da informação sobre a amostra é perdida nesse processo. Uma interpretação geométrica de \(|\mathbf{s}|\) nos ajudará a apreciar seus pontos fortes e fracos como um resumo descritivo.
Considere a área gerada no plano pelos dois vetores de desvio \(\mathbf{d}_1 = \mathbf{y}_1 - \bar{x}_1 \mathbf{1}\) e \(\mathbf{d}_2 = \mathbf{y}_2 - \bar{x}_2 \mathbf{1}\). Seja \(\lVert\mathbf{d}_1\rVert\) o comprimento de \(\mathbf{d}_1\) e \(\lVert\mathbf{d}_2\rVert\) o comprimento de \(\mathbf{d}_2\). Pela geometria elementar, temos o diagrama:
A área do trapézio é \(\lVert \mathbf{d}_1\rVert \sin(\theta) \lVert\mathbf{d}_2\rVert\). Usando \(\cos^2(\theta) + \sin^2(\theta) = 1\), podemos expressar essa área como:
\[ \text{Área} = \lVert\mathbf{d}_1\rVert \lVert\mathbf{d}_2\rVert \sqrt{1 - \cos^2(\theta)} \]
A partir de (3-5) e (3-7), temos:
\[ \lVert \mathbf{d}_1 \rVert^2 = \sum_{j=1}^{n} (x_{j1} - \bar{x}_1)^2 = (n - 1) s_{11} \]
\[ \lVert \mathbf{d}_2 \rVert^2 = \sum_{j=1}^{n} (x_{j2} - \bar{x}_2)^2 = (n - 1) s_{22} \]
E,
\[ \cos(\theta) = r_{12} \]
Portanto, a área pode ser expressa como:
\[ \text{Área} = (n - 1) \sqrt{s_{11} s_{22} (1 - r_{12}^2)} \tag{3-13} \]
Além disso,
\[ \dfrac{s_{12}}{\sqrt{s_{11} s_{22}}} = r_{12} \]
Assim,
\[ |\mathbf{s}| = s_{11} s_{22} (1 - r_{12}^2) \tag{3-14} \]
Se compararmos (3-14) com (3-13), veremos que
\[ |\mathbf{s}| = \left(\dfrac{\text{área}}{n - 1}\right)^2 \]
Assumindo agora que \(|\mathbf{s}| = (n - 1)^{-(p-1)}\text{volume}^2\) é válido para o volume gerado no espaço \(n\) pelas \(p - 1\) vetores de desvio \(\mathbf{d}_1, \mathbf{d}_2, \ldots, \mathbf{d}_{p-1}\), podemos estabelecer o seguinte resultado geral para \(p\) vetores de desvio por indução:
\[ \text{Variância amostral generalizada} = |\mathbf{s}| = (n - 1)^{-p} \;\text{volume}^2 \tag{3-15} \]
Esse resultado nos permite entender a variação geral dos \(p\) vetores de desvio no espaço \(n\)-dimensional.
Figura 3.6 (a) Grande variância generalizada amostral para p = 3. (b) Pequena variância generalizada para p = 3.
A variância generalizada também tem interpretações na representação gráfica de dispersão dos dados em \(p\)-espaços. A interpretação mais intuitiva diz respeito à dispersão em torno do ponto médio da amostra \(\bar{\mathbf{x}}^{\prime} = [\bar{x}_1, \bar{x}_2, \ldots, \bar{x}_2]\). Considere a medida de distância fornecida no comentário abaixo (2-19), com \(\bar{\mathbf{x}}\) desempenhando o papel do ponto fixo \(\boldsymbol{\mu}\) e \(\mathbf{S}^{-1}\) desempenhando o papel de \(A\). Com essas escolhas, as coordenadas \(\mathbf{x}^{\prime} = [x_1, x_2, \ldots, x_p]\) dos pontos a uma distância constante \(c\) de \(\bar{\mathbf{x}}\) satisfazem
\[ (\mathbf{x} - \bar{\mathbf{x}})^{\prime}\mathbf{s}^{-1}(\mathbf{x} - \bar{\mathbf{x}}) = c^2 \tag{3-16} \]
Quando \(p = 1\),
\[ (\mathbf{x} - \bar{\mathbf{x}})^{\prime}\mathbf{s}^{-1}(\mathbf{x} - \bar{\mathbf{x}}) = \left(\dfrac{x_1 - \bar{x}_1}{s_{1}}\right)^2=z_{1}^{2} \]
é a distância ao quadrado de \(x_1\) a \(\bar{x}_1\) em unidades de desvio-padrão.
A Equação (3-16) define um hiperelipsoide (uma elipse se \(p = 2\)) centrado em \(\bar{\mathbf{x}}\). Pode ser demonstrado usando cálculo diferencial que o volume desse hiperelipsoide está relacionado a \(|\mathbf{S}|\), em particular
\[ \text{Volume de }\{(\mathbf{x} - \bar{\mathbf{x}})^{\prime}\mathbf{s}^{-1}(\mathbf{x} - \bar{\mathbf{x}}) \leq c^2\} = k_p \;c^p \sqrt{|\mathbf{s}|} \tag{3-17} \]
ou
\[ (\text{Volume do elipsoide})^2 = (\text{constante}) \times (\text{variância amostral generalizada}) \]
em que a constante \(k_p = \dfrac{2 \pi^{p/2}}{p\Gamma(p/2)}\), sendo que \(\Gamma(z)\) denota a função gama avaliada em \(z\).
Um volume grande corresponde a uma grande variância generalizada.
Embora a variância generalizada tenha algumas interpretações geometricamente agradáveis, ela apresenta uma fraqueza básica como um resumo descritivo da matriz de covariância da amostra \(\mathbf{s}\), como mostra o seguinte exemplo.
Para \(p=2\) e nível de predição de 95%, a área da elipse é:
\[ \text{Área da elipse de 95%} \approx \dfrac{12\pi}{n-1}\sqrt{|\mathbf{s}|} \]
A Figura 3.7 mostra três gráficos de dispersão com padrões de correlação muito diferentes.
Todos os três conjuntos de dados têm \(\bar{\mathbf{x}}^{\prime} = [2,1]\), e as matrizes de covariância são:
\[ \begin{align} \mathbf{s} &= \begin{bmatrix} 5 & 4 \\ 4 & 5 \end{bmatrix} \quad r = 0.8\\ \mathbf{s} &= \begin{bmatrix} 3 & 0 \\ 0 & 3 \end{bmatrix} \quad r = 0\\ \mathbf{s} &= \begin{bmatrix} 5 & -4 \\ -4 & 5 \end{bmatrix} \quad r = -0.8 \end{align} \]
Cada matriz de covariância \(\mathbf{s}\) contém informações sobre a variabilidade das variáveis componentes e também as informações necessárias para calcular o coeficiente de correlação. Nesse sentido, \(\mathbf{s}\) captura a orientação e o tamanho do padrão de dispersão.
Os autovalores e autovetores extraídos de \(\mathbf{s}\) descrevem ainda mais o padrão no gráfico de dispersão.
Figura 3.7 Gráficos de dispersão com três orientações diferentes.
A seguir, para
\[ \mathbf{s} = \begin{bmatrix} 5 & 4 \\ 4 & 5 \end{bmatrix} \]
determinamos os pares de autovalor-autovetor como \(\lambda_1 = 9,\; \mathbf{e}_1^{\prime} = \left[\dfrac{1}{\sqrt{2}}, \dfrac{1}{\sqrt{2}}\right]\) e \(\lambda_2 = 1, \; \mathbf{e}_2^{\prime} = \left[\dfrac{1}{\sqrt{2}}, -\dfrac{1}{\sqrt{2}}\right]\).
A elipse centrada na média, com centro em \(\bar{\mathbf{x}}^{\prime} = [2,1]\) para os três casos, é dada por
\[ (\mathbf{x} - \bar{\mathbf{x}})^{\prime} \mathbf{s}^{-1} (\mathbf{x} - \bar{\mathbf{x}}) \le c^2 \]
Para descrever essa elipse, como na Seção 2.3, com \(\mathbf{A} = \mathbf{S}^{-1}\), notamos que se \((\lambda, \mathbf{e})\) é um par de autovalor-autovetor de \(\mathbf{s}\), então \((\lambda^{-1}, \mathbf{e})\) é um par de autovalor-autovetor para \(\mathbf{s}^{-1}\). Isto é, se \(\mathbf{s}\mathbf{e} = \lambda\mathbf{e}\), então multiplicando à esquerda por \(\mathbf{s}^{-1}\) temos \(\mathbf{s}^{-1}\mathbf{s}\mathbf{e} = \lambda\mathbf{s}^{-1}\mathbf{e}\), ou \(\mathbf{s}^{-1}\mathbf{e} = \lambda^{-1}\mathbf{e}\). Portanto, usando os autovalores de \(\mathbf{s}\), sabemos que a elipse se estende \(c \sqrt{\lambda_i}\) na direção de \(\mathbf{e}_i\) a partir de \(\bar{\mathbf{x}}\).
Em \(p = 2\) dimensões, a escolha \(c^2 = 5.99\) produzirá uma elipse que contém aproximadamente 95% das observações. Os vetores \(3\sqrt{5.99} \mathbf{e_1}\) e \(\sqrt{5.99} \mathbf{e_2}\) são desenhados na Figura 3.8(a). Observe como as direções são os eixos naturais para a elipse, e perceba que os comprimentos desses autovetores escalados são comparáveis ao tamanho do padrão em cada direção.
A seguir, para
\[ \mathbf{s} = \begin{bmatrix} 3 & 0 \\ 0 & 3 \end{bmatrix} \]
os autovalores satisfazem \(0 = (\lambda - 3)^2\) e escolhemos arbitrariamente os autovetores de forma que \(\lambda_1 = 3\), \(\mathbf{e}_1^{\prime} = [1, 0]\) e \(\lambda_2 = 3\), \(\mathbf{e}_2^{\prime} = [0, 1]\). Os vetores \(\sqrt{3} \sqrt{5.99}\mathbf{e}_1\) e \(\sqrt{3}\sqrt{5.99} \mathbf{e}_2\) são desenhados na Figura 3.8(b).
Figura 3.8 Eixos das elipses centradas na média que contêm 95% das observações para os gráficos de dispersão da Figura 3.7.
Por fim, para
\[ \mathbf{s} = \begin{bmatrix} 5 & -4 \\ -4 & 5 \end{bmatrix} \]
os autovalores satisfazem \(0 = (\lambda - 5)^2 - (-4)^2 = (\lambda - 9)(\lambda - 1)\) e determinamos os pares autovalor-autovetor \(\lambda_1 = 9\), \(\mathbf{e}_1^{\prime} = \left[\dfrac{1}{\sqrt{2}}, -\dfrac{1}{\sqrt{2}}\right]\) e \(\lambda_2 = 1\), \(\mathbf{e}_2^{\prime} = \left[\dfrac{1}{\sqrt{2}}, \dfrac{1}{\sqrt{2}}\right]\). Os autovetores escalados \(3\sqrt{5.99} \mathbf{e}_1\) e \(\sqrt{5.99} \mathbf{e}_2\) são mostrados na Figura 3.8(c).
Em duas dimensões, frequentemente podemos traçar os eixos da elipse centrada na média visualmente. No entanto, a abordagem dos autovetores também funciona em dimensões mais altas, onde os dados não podem ser examinados visualmente.
Observe que nossos três padrões de dispersão aparentam cobrir aproximadamente a mesma área. As elipses que resumem a variabilidade \[(\mathbf{x} - \bar{\mathbf{x}})^{\prime} \mathbf{s}^{-1} (\mathbf{x} - \bar{\mathbf{x}}) \le c^2\]
têm exatamente a mesma área [veja (3-17)], já que todas têm \(|\mathbf{s}| = 9\).
Como o Exemplo 3.8 demonstra, estruturas de correlação diferentes não são detectadas por \(|\mathbf{s}|\). A situação para \(p > 2\) pode ser ainda mais obscura.
Consequentemente, muitas vezes é desejável fornecer mais do que o único número \(|\mathbf{s}|\) como um resumo de \(\mathbf{s}\). A partir do Exercício 2.12, \(|\mathbf{s}|\) pode ser expresso como o produto dos autovalores de \(\mathbf{s}\). Além disso, o elipsoide centrado na média com base em \(\mathbf{s}^{-1}\) [veja (3-16)] tem eixos cujos comprimentos são proporcionais às raízes quadradas dos autovalores \(\lambda_i\) (veja a Seção 2.3). Portanto, esses autovalores fornecem informações sobre a variabilidade em todas as direções na representação de espaço \(p\) dos dados. É útil, portanto, relatar seus valores individuais, bem como o seu produto. Abordaremos esse tópico posteriormente quando discutirmos os componentes principais.
A variância amostral generalizada é indevidamente afetada pela variabilidade das medições em uma única variável. Por exemplo, suponha que algum \(s_{ii}\) seja grande ou bastante pequeno. Então, geometricamente, o vetor de desvio correspondente \(\mathbf{d}_i = \mathbf{y}_{i} - \bar{x}_i\mathbf{1}\) será muito longo ou muito curto e, portanto, claramente será um fator importante na determinação do volume. Consequentemente, às vezes é útil escalonar todos os vetores de desvio para que eles tenham o mesmo comprimento.
Escalar os vetores residuais é equivalente a substituir cada observação original \(x_{jk}\) por seu valor padronizado \(\left(x_{jk} - \bar{x}_k\right)/\sqrt{s_{kk}}\). A matriz de covariância amostral das variáveis padronizadas é então \(\mathbf{r}\), a matriz de correlação amostral das variáveis originais. (Veja o Exercício 3.13.)
Então:
\[ \text{Variância amostral generalizada das variáveis padronizadas}=|\mathbf{r}|\tag{3-19} \]
O volume gerado no espaço \(p\) pelas vetores de desvio pode ser relacionado à variância amostral generalizada. Os mesmos passos que levam a (3-15) produzem:
\[ \text{Variância amostral generalizada das variáveis padronizadas}= | \mathbf{r} | = (n - 1)^{-p} \text{volume}^2 \tag{3-20} \]
em que \(n\) é o número de
observações e volume
é o volume no espaço \(p\) gerado pelos vetores de desvio.
Figura 3.10: O volume gerado pelos vetores de desvio de comprimento igual das variáveis padronizadas.
O volume gerado pelos vetores de desvio das variáveis padronizadas é ilustrado na Figura 3.10 para os dois conjuntos de vetores de desvio plotados na Figura 3.6. Uma comparação entre as Figuras 3.10 e 3.6 revela que a influência do vetor \(\mathbf{d}_2\) (alta variabilidade em \(x_2\)) sobre o volume ao quadrado \(|\mathbf{s}|\) é muito maior do que sua influência sobre o volume ao quadrado \(|\mathbf{r}|\).
As quantidades \(|\mathbf{s}|\) e \(|\mathbf{r}|\) estão conectadas pela relação:
\[ |\mathbf{s}| = |\mathbf{r}|\prod_{i=1}^{p}{s_{ii}} \tag{3-21} \]
Portanto,
\[ (n-1)^p |\mathbf{s}| = (n-1)^p |\mathbf{r}|\prod_{i=1}^{p}{s_{ii}} \tag{3-22} \]
[A prova de (3-21) é deixada para o leitor como Exercício 3.12.]
Interpretando (3-22) em termos de volumes, podemos ver a partir de (3-15) e (3-20) que o volume ao quadrado \((n - 1)^p |\mathbf{s}|\) é proporcional ao volume ao quadrado \((n - 1)^p |\mathbf{r}|\). A constante de proporcionalidade é o produto das variâncias, que, por sua vez, é proporcional ao produto dos quadrados dos comprimentos \((n - 1)s_{ii}\) dos \(\mathbf{d}_i\).
A equação (3-21) mostra, de forma algébrica, como uma mudança na escala de medição, por exemplo, afetará a relação entre as variâncias generalizadas. Uma vez que \(|\mathbf{r}|\) é baseado em medições padronizadas, ele não é afetado pela mudança na escala. No entanto, o valor relativo de \(|\mathbf{s}|\) será alterado sempre que o fator multiplicativo \(s_{ii}\) mudar.
Vamos ilustrar a relação em (3-21) para as variâncias generalizadas \(|\mathbf{s}|\) e \(|\mathbf{r}|\) quando \(p = 3\). Suponha:
\[ \mathbf{s} = \begin{bmatrix} 4 & 3 & 1 \\ 3 & 9 & 2 \\ 1 & 2 & 1 \end{bmatrix} \]
Então, temos \(s_{11} = 4\), \(s_{22} = 9\) e \(s_{33} = 1\). Além disso, usando a Definição 2A.24, obtemos:
\[ |\mathbf{s}| = 14\\ |\mathbf{r}|=\dfrac{7}{18} \]
Assim, temos (verificação):
\[ 14=|\mathbf{s}|= s_{11} s_{22} s_{33} |\mathbf{r}| = 4\times 9 \times 1 \times\dfrac{7}{18}= 14 \]
\[\Diamond\]
Concluímos esta discussão mencionando outra generalização de variância. Especificamente, definimos a variância amostral total como a soma dos elementos diagonais da matriz de covariância amostral \(\mathbf{s}\). Assim,
\[ \text{Variância amostral total} =\mathrm{tr}(\mathbf{s})= \sum_{i=1}^{p}{s_{ii}}=\dfrac{1}{n-1}\sum_{i=1}^{p}{\mathbf{d}_i} \tag{3-23} \]
Nessa expressão, \(s_{ii}\) representa os elementos diagonais da matriz de covariância amostral \(\mathbf{s}\) para \(p\) variáveis. Essa medida representa a soma das variâncias amostrais individuais de cada variável (mesmo que elas tenham diferentes unidades de medida!).
Geometricamente, a variância amostral total é a soma dos quadrados dos comprimentos dos \(p\) vetores de desvio \(\mathbf{d}_i = (\mathbf{y}_1 - \bar{x}_1\mathbf{1}), \ldots, \mathbf{d}_p = (\mathbf{y}_p - \bar{x}_p\mathbf{1})\), dividida por \(n - 1\). O critério da variância amostral total não leva em consideração a orientação (estrutura de correlação) dos vetores residuais. Por exemplo, ele atribui os mesmos valores para ambos os conjuntos de vetores residuais (a) e (b) na Figura 3.6.
Em outras palavras, a variância amostral total é uma medida de dispersão que leva em conta apenas os comprimentos dos vetores de desvio, ignorando qualquer relação de correlação entre as variáveis originais. Isso pode ser útil em algumas situações, mas é importante notar que pode perder informações importantes sobre a relação entre as variáveis. Outras medidas, como a matriz de covariância ou a matriz de correlação, são usadas quando se deseja levar em conta a estrutura de correlação das variáveis.
Desenvolvemos representações geométricas da matriz de dados \(\boldsymbol{\mathcal{x}}\) e das estatísticas descritivas derivadas \(\bar{\mathbf{x}}\) e \(\mathbf{s}\). Além disso, é possível vincular algebraicamente o cálculo de \(\bar{\mathbf{x}}\) e \(\mathbf{s}\) diretamente a \(\mathbf{x}\) usando operações matriciais. As expressões resultantes, que descrevem a relação entre \(\bar{\mathbf{x}}\), \(\mathbf{s}\) e o conjunto completo de dados \(\mathbf{s}\) de forma concisa, podem ser facilmente programadas em computadores eletrônicos.
Temos que \(\bar{x}_i = \dfrac{1}{n}\left(1x_{1i} + 1x_{2i} + \cdots + 1x_{ni}\right) = \dfrac{1}{n}\mathbf{y}_{i}^{\prime}\mathbf{1}\). Portanto,
\[ \begin{align} \bar{\mathbf{x}} &= \begin{bmatrix} \bar{x}_1 \\ \bar{x}_2 \\ \vdots \\ \bar{x}_p \end{bmatrix} \\ &= \begin{bmatrix} \dfrac{1}{n}\mathbf{y}_{1}^{\prime}\mathbf{1} \\ \dfrac{1}{n}\mathbf{y}_{i}^{\prime}\mathbf{2} \\ \vdots \\ \dfrac{1}{n}\mathbf{y}_{i}^{\prime}\mathbf{p} \end{bmatrix} \\ &= \dfrac{1}{n} \begin{bmatrix} x_{11} & x_{12} & \ldots & x_{1p} \\ x_{21} & x_{22} & \ldots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \ldots & x_{np} \end{bmatrix}^{\prime} \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} \\ \bar{\mathbf{x}} &=\dfrac{1}{n}\boldsymbol{\mathcal{x}}^{\prime}\mathbf{1}\\\\ \mathbf{s}&=\dfrac{1}{n-1} \boldsymbol{\mathcal{x}}^{\prime} \mathbf{\mathcal{I}}\boldsymbol{\mathcal{x}}\\ \end{align} \]
em que \(\mathbf{\mathcal{I}}=\mathbf{I}-\dfrac{1}{n}\mathbf{1}\mathbf{1}^{\prime}\) sendo que \(\mathbf{I}\) é matriz identidade.
A matriz \(\mathbf{\mathcal{I}}\) é idempotente, i.e., \(\mathbf{\mathcal{I}}^{\prime}\mathbf{\mathcal{I}}=\mathbf{\mathcal{I}}\). Além disso, o traço de \(\mathbf{\mathcal{I}}\) é igual ao número de graus de liberdade associados a \(\mathbf{s}\):
\[ \mathrm{tr}(\mathbf{\mathcal{I}})=n-1 \]
Portanto:
\[ \mathbf{s}=\dfrac{\boldsymbol{\mathcal{x}}^{\prime} \mathbf{\mathcal{I}}\boldsymbol{\mathcal{x}}}{\mathrm{tr}(\mathbf{\mathcal{I}})} \]
A matriz \(\mathbf{s}\) pode ser expressa em função de \(\mathbf{r}\):
\[ \begin{align} \mathbf{s}&=\mathbf{v}^{1/2}\;\mathbf{r}\; \mathbf{v}^{1/2}\\ \mathbf{v}&=\text{diag}[s_1\;s_2\;\cdots\;s_p] \end{align} \]
Resultado 3.5. As combinações lineares
\[ \mathbf{b}^{\prime}\mathbf{x} = b_1x_1 + b_2x_2 + \cdots + b_px_p \\ \mathbf{c}^{\prime}\mathbf{x} = c_1x_1 + c_2x_2 + \cdots + c_px_p \]
têm médias amostrais, variâncias e covariâncias que estão relacionadas com \(\bar{\mathbf{x}}\) e \(\mathbf{s}\) por:
Média amostral de \(\mathbf{b}^{\prime}\mathbf{x} = \mathbf{b}^{\prime} \bar{\mathbf{x}}\)
Média amostral de \(\mathbf{c}^{\prime}\mathbf{x} = \mathbf{c}^{\prime} \bar{\mathbf{x}}\)
Variância amostral de \(\mathbf{b}^{\prime}\mathbf{x} = \mathbf{b}^{\prime}\mathbf{s}\mathbf{b}\)
Variância amostral de \(\mathbf{c}^{\prime}\mathbf{x} = \mathbf{c}^{\prime}\mathbf{s}\mathbf{c}\)
Covariância amostral de \(\mathbf{b}^{\prime}\mathbf{x}\) e \(\mathbf{c}^{\prime}\mathbf{x} = \mathbf{b}^{\prime}\mathbf{s}\mathbf{c}\)
\[\tag{3-36}\]
As relações de média amostral e covariância no Resultado 3.5 são aplicáveis a qualquer número de combinações lineares. Considere as \(q\) combinações lineares:
\[ a_1x_1 + a_2x_2 + \ldots + a_px_p, \;i = 1,2,\ldots,q \tag{3-37} \]
Essas combinações lineares podem ser expressas em notação matricial como:
\[ \begin{align} \begin{bmatrix} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1p}x_p \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2p}x_p \\ \vdots \\ a_{q1}x_1 + a_{q2}x_2 + \cdots + a_{qp}x_p \\ \end{bmatrix} &= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1p} \\ a_{21} & a_{22} & \cdots & a_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ a_{q1} & a_{q2} & \cdots & a_{qp} \\ \end{bmatrix} \begin{bmatrix} X_1 \\ X_2 \\ \vdots \\ X_p \end{bmatrix}\\ &= \mathbf{A}\mathbf{x} \end{align} \tag{3-38} \]
Resultado 3.6. As \(q\) combinações lineares \(\mathbf{Ax}\) em (3-38) têm vetor médio amostral \(\mathbf{A}\bar{\mathbf{x}}\) e matriz de covariância amostral \(\mathbf{AsA}^{\prime}\).
3.2, 3.5, 3.6, 3.7, 3.8, 3.14, 3.16, 3.18, 3.19, 3.20
Uma generalização da familiar densidade normal em formato de sino para várias dimensões desempenha um papel fundamental na análise multivariada. Na verdade, a maioria das técnicas encontradas neste livro é baseada na suposição de que os dados foram gerados a partir de uma distribuição normal multivariada. Embora dados reais nunca sejam exatamente multivariados normais, a densidade normal frequentemente é uma aproximação útil da distribuição “verdadeira” da população.
Uma vantagem da distribuição normal multivariada é que ela é matematicamente tratável e permite a obtenção de resultados “bons”. Isso nem sempre é o caso para outras distribuições geradoras de dados. Claro que a atratividade matemática por si só tem pouco valor para o praticante. No entanto, verifica-se que as distribuições normais são úteis na prática por duas razões: Primeiro, a distribuição normal serve como um modelo de população adequado em alguns casos; segundo, as distribuições amostrais de muitas estatísticas multivariadas são aproximadamente normais, independentemente da forma da população de origem, devido ao efeito do teorema do limite central.
Resumindo, muitos problemas do mundo real naturalmente se encaixam no contexto da teoria normal. A importância da distribuição normal reside em seu papel duplo como modelo de população para certos fenômenos naturais e como uma distribuição amostral aproximada para muitas estatísticas.
Além disso, a distribuição normal multivariada é:
A densidade normal multivariada é uma generalização da densidade normal univariada para \(p\) dimensões. Lembre-se de que a distribuição normal univariada, com média \(\mu\) e variância \(\sigma^2\), i.e., \(X \sim \mathcal{N}(\mu,\sigma^2)\) tem a função densidade de probabilidade:
\[ \begin{align} f(x) &= \dfrac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\dfrac{1}{2}z^2\right)\\ z&=\dfrac{x-\mu}{\sigma}, \quad x \in \mathbb{R} \end{align} \tag{4-1} \]
Figura 4.4: Uma densidade normal com vetor de média μ e variância σ2 e áreas selecionadas sob a curva.
Um gráfico dessa função resulta na conhecida curva em formato de sino mostrada na Figura 4.1. Na figura, também são apresentadas as áreas aproximadas sob a curva dentro de \(\pm1\) desvio-padrão e \(\pm2\) desvios-padrão da média. Essas áreas representam probabilidades e, portanto, para a variável aleatória normal \(X\):
\[ \begin{align} P(\mu - \sigma \leq X \leq \mu + \sigma) &\approx 0.6827\\ P(\mu - 2\sigma \leq X \leq \mu + 2\sigma) &\approx 0.9545\\ P(\mu - 3\sigma \leq X \leq \mu + 3\sigma) &\approx 0.9973 \end{align} \]
# Normal_Uni_Bi.R
# Distribuição Normal Padrão
set.seed(123)
xseq <- seq(-4, 4, 0.01)
densities <- dnorm(xseq, 0, 1)
cumulative <- pnorm(xseq, 0, 1)
randomdeviates <- rnorm(1e3, 0, 1)
par(mfrow=c(2,2))
plot(xseq,
densities,
col="black",
xlab="",
ylab="Density",
type="l",
lwd=2,
cex=2,
main="PDF of Standard Normal",
cex.axis=0.8)
plot(xseq,
cumulative,
col="black",
xlab="",
ylab="Cumulative Probability",
type="l",
lwd=2,
cex=2,
main="CDF of Standard Normal",
cex.axis=0.8)
boxplot(randomdeviates,
horizontal=TRUE,
main="Random draws from Std Normal",
cex.axis=0.8,
xlim=c(-4,4))
plot(density(randomdeviates),
main="Random draws from Std Normal",
cex.axis=0.8,
xlim=c(-4,4))
É conveniente denotar a função densidade normal com média \(\mu\) e variância \(\sigma^2\) por \(\phi(\mu,\sigma^2)\). Portanto, \(\mathcal{N}(10, 4)\) se refere à função em (4-1) com \(\mu = 10\) e \(\sigma^2 = 4\). Essa notação será estendida para o caso multivariado mais tarde.
O termo
\[ \begin{align} z^2 &=\left(\dfrac{x-\mu}{\sigma}\right)^2\\ z^2 &=(x-\mu)\left(\sigma^2\right)^{-1}(x-\mu) \end{align} \tag{4-2} \]
no expoente da função densidade normal univariada mede o quadrado da distância de \(x\) para \(\mu\) em unidades de desvio-padrão: é o escore-z ao quadrado. Isso pode ser generalizado para um vetor \(p \times 1\) \(\mathbf{x}\) de observações em várias variáveis pela distância de Mahalanobis ao quadrado:
\[ d_{M}^{2}=(\mathbf{x} - \boldsymbol{\mu})^{\prime} \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu}) \tag{4-3} \]
O vetor \(\boldsymbol{\mu}\) de tamanho \(p \times 1\) representa o valor esperado do vetor aleatório \(\mathbf{X}\), e a matriz \(\mathbf{\Sigma}\) de tamanho \(p \times p\) é a matriz de covariância de \(\mathbf{X}\). [Ver (2-30) e (2-31).] Assumiremos que a matriz simétrica \(\mathbf{\Sigma}\) é definida positiva, então a expressão em (4-3) é o quadrado da distância generalizada de \(\mathbf{x}\) para \(\boldsymbol{\mu}\).
A densidade normal multivariada é obtida substituindo a distância univariada em (4-2) pela distância generalizada multivariada em (4-3) na função de densidade de (4-1). Quando essa substituição é feita, a constante de normalização univariada \(\dfrac{1}{\sqrt{2\pi\sigma^2}}\) deve ser alterada para uma constante mais geral que torne o volume sob a superfície da função de densidade multivariada igual a um para qualquer valor de \(p\). Isso é necessário porque, no caso multivariado, as probabilidades são representadas por volumes sob a superfície em regiões definidas por intervalos de valores. Pode ser mostrado (ver [1]) que essa constante é \(|\mathbf{\Sigma}|^{1/2}\), e consequentemente, uma densidade normal \(p\)-dimensional para o vetor aleatório \(\mathbf{X} = [X_1\; X_2\; \cdots\; X_p]^{\prime}\) tem a forma:
\[ f(\mathbf{x}) = \dfrac{1}{\sqrt{|2\pi\Sigma|}}\exp\left(-\dfrac{1}{2}d_{M}^{2}\right) \tag{4-4} \]
em que \(\mathbf{x} \in \mathbb{R}^{p}\). Denotaremos esta densidade normal \(p\)-dimensional por \(\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\), que é análoga à densidade normal no caso univariado.
Portanto:
\[ \mathbf{X} = [X_1\; X_2\; \cdots\; X_p]^{\prime}\sim \mathcal{N}_p(\boldsymbol{\mu},\mathbf{\Sigma}) \]
# Normal_Uni_Bi.R
# Criar uma grade de pontos para representar a distribuição bivariada
X1 <- seq(-3, 3, length.out=100)
X2 <- seq(-3, 3, length.out=100)
grid_points <- expand.grid(X1=X1, X2=X2)
# Definir a média (mu) e matriz de covariância (Sigma) da distribuição bivariada
mu <- c(0, 0) # Vetor de médias
sigma <- matrix(c(1, 0.5,
0.5, 1),
nrow=2) # Matriz de covariância
# Calcular a densidade normal bivariada para cada ponto da grade
z <- mvtnorm::dmvnorm(grid_points,
mean=mu,
sigma=sigma)
z_matrix <- matrix(z, nrow = length(X1))
contour(X1, X2, z_matrix,
main="Normal bivariada: curvas de nível",
xlab="X1", ylab="X2")
persp(X1, X2, z_matrix, theta=30, phi=30, col="white",
expand=0.8, ticktype="detailed",
xlab = "X1", ylab = "X2", zlab = "Densidade",
main = "FDP Binormal")
# grid_points$z <- z
# plotly::plot_ly(data = grid_points,
# x = ~X1, y = ~X2, z = ~z,
# type = "contour",
# colorscale = "Viridis",
# showscale = FALSE,
# contours = list(coloring = "lines"))
# fig <- plotly::plot_ly(x = X1, y = X2, z = z_matrix,
# type = "surface",
# colorscale = "Viridis",
# showscale = FALSE)
# plotly::layout(fig,
# scene = list(xaxis_title = "X1",
# yaxis_title = "X2",
# zaxis_title = "Densidade"))
# Definindo a função de densidade da normal bivariada
bivariate_normal <- function(x1, x2, mu1, mu2, sigma1, sigma2, rho) {
z1 <- (x1 - mu1) / sigma1
z2 <- (x2 - mu2) / sigma2
1 / (2 * pi * sigma1 * sigma2 * sqrt(1 - rho^2)) *
exp(-1 / (2 * (1 - rho^2)) * (z1^2 + z2^2 - 2 * rho * z1 * z2))
}
# Parâmetros da normal bivariada
mu1 <- 0
mu2 <- 0
sigma1 <- 1
sigma2 <- 1
rho <- 0.5
# Criando a grade de valores
x1_values <- base::seq(-3, 3, length.out = 100)
x2_values <- base::seq(-3, 3, length.out = 100)
grid <- base::expand.grid(x1 = x1_values, x2 = x2_values)
grid$z <- base::with(grid, bivariate_normal(x1, x2, mu1, mu2, sigma1, sigma2, rho))
# Transformando os dados para o formato longo
grid_melt <- reshape2::melt(grid, id.vars = c("x1", "x2"), measure.vars = "z")
# Gráfico de contorno com ggplot2
ggplot_contour <- ggplot2::ggplot(grid_melt, ggplot2::aes(x = x1, y = x2, z = value)) +
ggplot2::geom_tile(ggplot2::aes(fill = value)) +
ggplot2::geom_contour(color = "black") +
ggplot2::scale_fill_gradientn(colours = grDevices::terrain.colors(10)) +
ggplot2::labs(title = "Gráfico de Contorno", x = "x1", y = "x2", fill = "Densidade") +
ggplot2::theme_minimal()
base::print(ggplot_contour)
# Renderizando a superfície 3D com rayshader
rayshader::plot_gg(ggplot_contour, multicore = TRUE, zoom = 0.8, phi = 40, theta = 30)
rayshader::render_snapshot(clear = FALSE,
background = "white")
# Gráfico de superfície 3D com plotly
fig <- plotly::plot_ly(x = ~x1_values, y = ~x2_values, z = ~base::matrix(grid$z, nrow = 100))
fig <- plotly::add_surface(fig, contours = list(
z = list(
show = TRUE,
usecolormap = TRUE,
highlightcolor = "#ff0000",
project = list(z = TRUE)
)
))
fig <- plotly::layout(fig, scene = list(
xaxis = list(title = "x1"),
yaxis = list(title = "x2"),
zaxis = list(title = "Densidade"),
camera = list(
eye = list(x = 1.87, y = 0.88, z = -0.64)
)
), title = "Superfície 3D da Normal Bivariada")
fig
# Simular dados de uma distribuição normal bivariada
set.seed(123)
n <- 5e3
mu <- c(0, 0)
sigma <- matrix(c(1, 0.5, 0.5, 1), nrow = 2)
simulated_data <- mvtnorm::rmvnorm(n=n, mean=mu, sigma=sigma)
centroidesimul <- colMeans(simulated_data)
# simulated_data <- MASS::mvrnorm(n=n, mu=mu, Sigma=sigma)
# colMeans(simulated_data)
bgp <- DescTools::PlotBag(simulated_data[,1],simulated_data[,2],
main=paste("Bagplot: normal bivariada simulada"),
xlab="X1",
ylab="X2",
na.rm=TRUE,
show.bagpoints=FALSE,
show.looppoints=FALSE,
show.whiskers=FALSE,
show.baghull=TRUE,
col.loophull="white",
col.looppoints="black",
col.baghull="white",
col.bagpoints="black",
asp=1,
add=FALSE,
cex=1)
print(outliers <- as.data.frame(bgp$pxy.outlier))
x y
1 -0.3783736 2.8044520
2 -1.7163324 1.8456164
3 3.4351064 3.4597384
4 2.7148861 3.7671749
5 0.3226809 3.0366593
6 -2.4603639 1.0161382
7 -0.7358310 2.8701042
8 1.4041986 -2.2399320
9 -2.3627500 1.0525603
10 0.8655797 3.1653315
11 3.1261506 0.1364389
12 1.8754778 3.4034411
13 4.0005987 2.4973937
14 -2.7114200 -3.4404977
15 1.2145452 3.2625354
16 -1.3154864 -3.4132021
17 2.0281645 -1.2293051
18 -4.2275025 -2.9105619
19 0.4015734 3.5574154
20 3.2753660 0.8542569
21 -2.5639672 -3.7958796
22 1.9859254 -1.3522451
23 -1.6099207 2.1477640
for (o in 1:nrow(outliers))
{
r <- which(simulated_data[,1]==outliers$x[o] &
simulated_data[,2]==outliers$y[o])
text(outliers$x[o],outliers$y[o], r, pos=1, cex=0.7)
}
# Criar o gráfico da distribuição normal bivariada com as elipses de predição
car::dataEllipse(simulated_data,
plot.points=FALSE,
robust=TRUE,
grid=FALSE,
levels=c(0.50, 0.999),
add=FALSE,
col="black",
xlab = "X1", ylab = "X2",
main = "Distribuição Binormal com Elipses de Predição\nde 50% e 99.9%")
abline(v=0,h=0,lty=2)
car::scatterplot(simulated_data[,1],simulated_data[,2],
regLine=FALSE,
ellipse=list(levels=c(0.5, 0.999),
robust=TRUE, fill=FALSE),
grid=FALSE,
smooth=FALSE,
asp=1,
col=c("black"),
cex=0.1,
pch=1,
xlim = c(-4, 4), ylim = c(-4, 4),
xlab = "X1", ylab = "X2",
main = "Distribuição Binormal com Elipses de Predição\nde 50% e 99.9%")
kd <- MASS::kde2d(simulated_data[,1],simulated_data[,2], n = 50)
fig <- plotly::plot_ly(x = kd$x, y = kd$y, z = kd$z)
plotly::add_contour(fig, color = I("black"))
result <- MVN::mvn(data=data.frame(simulated_data[,1],
simulated_data[,2]),
univariate_test="SW")
print(result$multivariate_normality)
Test Statistic p.value Method MVN
1 Henze-Zirkler 0.646 0.884 asymptotic ✓ Normal
Test Variable Statistic p.value Normality
1 Shapiro-Wilk simulated_data...1. 1 0.682 ✓ Normal
2 Shapiro-Wilk simulated_data...2. 1 0.957 ✓ Normal
# MVN::mvn(data.frame(simulated_data[,1],
# simulated_data[,2]),
# mvn_test = "hz",
# multivariate_plot = "contour")
# MVN::mvn(data.frame(simulated_data[,1],
# simulated_data[,2]),
# mvn_test = "hz",
# multivariate_plot = "persp")
# Distribuição Binormal Padrão (teoria)
bvn <- HH::bivariateNormal(0.6)
print(bvn)
# if (interactive())
# shiny::runApp(file.path(system.file(package="HH"),
# "shiny/bivariateNormal"))
# # Correlação sob binormal (teoria)
# if (interactive())
# shiny::runApp(system.file("shiny/bivariateNormalScatterplot",
# package="HH"))
Vamos avaliar a densidade normal bivariada para \(p = 2\) em termos dos parâmetros individuais: \(\mu_1 = \mathbb{E}(X_1)\), \(\mu_2 = \mathbb{E}(X_2)\), \(\sigma_{11} = \mathbb{V}(X_1)\), \(\sigma_{22} = \mathbb{V}(X_2)\) e \(\rho_{12} = \sigma_{12}/(\sqrt{\sigma_{11}} \sqrt{\sigma_{22}})=\text{Corr}(X_1,X_2)\).
Usando o Resultado 2A.8 (a inversa da matriz de covariância), obtemos que a inversa da matriz de covariância é:
\[ \mathbf{\Sigma}^{-1} = \dfrac{1}{\sigma_{11}\sigma_{22}-\sigma_{12}^{2}} \begin{bmatrix} \sigma_{22} & -\sigma_{12} \\ -\sigma_{12} & \sigma_{11} \end{bmatrix} \]
Introduzindo o coeficiente de correlação \(\rho_{12}\), escrevendo \(\sigma_{12} = \rho_{12}\sqrt{\sigma_{11}}\sqrt{\sigma_{22}}\), temos que \(\sigma_{11}\sigma_{22}-\sigma_{12}^{2} = \sigma_{11} \sigma_{22}\left(1- \rho_{12}^{2}\right)\), e a distância ao quadrado torna-se:
\[ (x - \boldsymbol{\mu})^{\prime} \mathbf{\Sigma}^{-1} (x - \boldsymbol{\mu}) = \dfrac{1}{1-\rho_{12}^{2}}\left( z_{1}^{2}+z_{2}^{2}-2\rho_{12}z_{1}z_{2}\right) \tag{4-5} \]
em que \(z_{i}=\dfrac{x_i-\mu_i}{\sqrt{\sigma_{ii}}}, \; i=1,2\).
Em seguida, dado que \(| \mathbf{\Sigma}| = \sigma_{11} \sigma_{22} - \rho_{12}^2 = \sigma_{11} \sigma_{22}(1 - \rho_{12}^2)\), podemos substituir por \(\mathbf{\Sigma}^{-1}\) e \(| \mathbf{\Sigma}|\) na equação (4-4) para obter a expressão da densidade normal bivariada (\(p = 2\)) envolvendo os parâmetros individuais \(\mu_1, \mu_2, \sigma_{11}, \sigma_{22}\) e \(\rho_{12}\):
\[ \begin{align} f(x_1,x_2) &= \dfrac{1}{2 \pi \sqrt{\sigma_{11} \sigma_{22}(1 - \rho_{12}^2)}} \exp \left( -\dfrac{1}{2(1 - \rho_{12}^2)} \left( z_{1}^{2}+z_{1}^{2}-2\rho_{12}z_{1}z_{1} \right) \right)\\ z_1&=\dfrac{x_1-\mu_1}{\sigma_1}\\ z_2&=\dfrac{x_2-\mu_2}{\sigma_2} \end{align} \tag{4-6} \]
A expressão em (4-6) é um pouco complicada, e a forma geral compacta em (4-4) é mais informativa de várias maneiras. Por outro lado, a expressão em (4-6) é útil para discutir certas propriedades da distribuição normal. Por exemplo, se as variáveis aleatórias \(X_1\) e \(X_2\) são não correlacionadas, de modo que \(\rho_{12} = 0\), a densidade conjunta pode ser escrita como o produto de duas densidades normais univariadas, cada uma no formato de (4-1).
Isso significa que \(f(x_1, x_2) = f(x_1)f(x_2)\) e \(X_1\) e \(X_2\) são independentes. [Ver (2-28).] Esse resultado é válido de forma geral (Ver Resultado 4.5).
Duas distribuições bivariadas com \(\sigma_{11} = \sigma_{22}\) são mostradas na Figura 4.2. Na Figura 4.2(a), \(X_1\) e \(X_2\) são independentes (\(\rho_{12} = 0\)). Na Figura 4.2(b), \(\rho_{12} = 0.75\). Note como a presença de correlação faz com que a probabilidade se concentre ao longo de uma linha.
Figura 4.2 Duas distribuições normais bivariadas.
Figura 4.2b Distribuição normal bivariada, marginais normais e condicional normal.
A partir da expressão em (4-4) para a densidade de uma variável normal p-dimensional, fica claro que as trajetórias dos valores de x que resultam em uma altura constante para a densidade são elipsoides. Ou seja, a densidade normal multivariada é constante em superfícies onde o quadrado da distância \((\mathbf{x} - \boldsymbol{\mu})^{\prime} \mathbf{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\) é constante. Essas trajetórias são chamadas de curvas de nível:
Curva de nível de densidade de probabilidade constante = \(\forall \mathbf{x} \;|\; (\mathbf{x} - \boldsymbol{\mu})^{\prime} \mathbf{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) = c^2\) = superfície de um elipsoide centrado em \(\boldsymbol{\mu}\)
Os eixos de cada elipsoide de densidade constante estão na direção dos autovetores de \(\mathbf{\Sigma}^{-1}\), e seus comprimentos são proporcionais aos recíprocos das raízes quadradas dos autovalores de \(\mathbf{\Sigma}^{-1}\). Felizmente, podemos evitar o cálculo de \(\mathbf{\Sigma}^{-1}\) ao determinar os eixos, pois esses elipsoides também são determinados pelos autovalores e autovetores de \(\mathbf{\Sigma}\). Apresentamos formalmente a correspondência para referência futura.
Resultado 4.1. \(\mathbf{\Sigma}\) é simétrica e definida positiva, i.e., \(\mathbf{\Sigma}>0\), de modo que \(\mathbf{\Sigma}^{-1}\) existe e também é definida positiva, i.e., \(\mathbf{\Sigma}^{-1}>0\), então:
\[ \mathbf{\Sigma}\mathbf{e} = \lambda\mathbf{e}\\\\ \lambda>0\\ \lVert \mathbf{e} \rVert=1 \]
implica
\[ \mathbf{\Sigma}^{-1}\mathbf{e} = \dfrac{1}{\lambda}\mathbf{e} \]
ou seja, \((\lambda, \mathbf{e})\) é um par de autovalor-autovetor para \(\mathbf{\Sigma}\) correspondendo ao par \((\lambda^{-1}, \mathbf{e})\) para \(\mathbf{\Sigma}^{-1}\).
A seguir, resumimos esses conceitos:
As curvas de nível de densidade constante para a distribuição normal \(p\)-dimensional são elipsoides definidos por \(\mathbf{x}\) tal que
\[ (\mathbf{x} - \boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1}(\mathbf{x} - \boldsymbol{\mu}) = c^2 \tag{4-7} \]
Esses elipsoides são centrados em \(\boldsymbol{\mu}\) e têm eixos \(\pm c\sqrt{\lambda_i}\mathbf{e}_i\), em que \(\mathbf{e}_i\) são os autovetores positivos de \(\mathbf{\Sigma}\) e \(\lambda_i\) são os autovalores ortonormais correspondentes para \(i = 1, 2, \ldots, p\).
A extração de autovalores e autovetores de uma matriz de covariância está relacionada ao conceito de transformações lineares e suas invariâncias. Vamos explicar matematicamente esse processo:
Seja \(\mathbf{\Sigma}\) uma matriz de covariância \(p \times p\), simétrica e definida positiva, em que \(p\) é o número de variáveis. As entradas da matriz \(\mathbf{\Sigma}\) representam as covariâncias entre as diferentes variáveis.
Os autovalores \(\lambda_i\) e os autovetores \(\mathbf{v}_i\) da matriz \(\mathbf{\Sigma}\) são encontrados ao resolver a equação característica:
\[ |\mathbf{\Sigma} - \lambda \mathbf{I}| = 0 \]
em que \(\mathbf{I}\) é a matriz identidade \(p \times p\). A equação característica é obtida ao calcular o determinante da matriz \((\mathbf{\Sigma} - \lambda \mathbf{I})\) igual a zero, o que nos dá os autovalores \(\lambda_i\).
Após obter os autovalores, os autovetores correspondentes são encontrados resolvendo o sistema de equações homogêneas:
\[ (\mathbf{\Sigma} - \lambda \mathbf{I})\mathbf{e} = \mathbf{0} \]
As duas equações homogêneas e a equação característica se relacionam da seguinte maneira:
\[ \begin{align}\ \mathbf{\Sigma e} &= \lambda \mathbf{e} \\ \mathbf{\Sigma e} - \lambda \mathbf{e} &= \mathbf{0}\\ (\mathbf{\Sigma} - \lambda \mathbf{I})\mathbf{e} &= \mathbf{0}\\ |\mathbf{\Sigma} - \lambda \mathbf{I}| &= 0 \end{align} \]
Os autovetores \(\mathbf{e}_i\) são os vetores não nulos que satisfazem essa equação, e eles representam as direções (ou eixos) principais do espaço de dados em que a matriz de covariância é diagonal (não tem covariância entre as variáveis).
Esses autovetores são ortonormais, o que significa que eles formam uma base ortonormal do espaço de dados. A magnitude dos autovalores \(\lambda_i\) representa a variância dos dados ao longo dos autovetores \(\mathbf{e}_i\). Os maiores autovalores correspondem às direções de maior variância, enquanto os menores autovalores correspondem às direções de menor variância.
Sim, você está correto! Os autovetores de uma matriz simétrica (como uma matriz de covariância) são ortonormais, o que significa que eles são vetores perpendiculares entre si e têm norma igual a 1.
Matematicamente, se \(\mathbf{\Sigma}\) é uma matriz simétrica e \(\mathbf{e}_i\) é um autovetor associado ao autovalor \(\lambda_i\), então os seguintes dois critérios são satisfeitos:
Ortogonalidade: \(\mathbf{e}_i^{\prime} \mathbf{e}_j = 0\) para \(i \neq j\). Isso significa que os autovetores são perpendiculares entre si, ou seja, o produto interno entre eles é zero.
Normalização: \(\|\mathbf{e}_i\| = 1\). Isso significa que os autovetores têm norma igual a 1, ou seja, eles são vetores unitários.
Essas propriedades tornam os autovetores ortonormais uma base ortonormal do espaço de dados, o que é muito útil em análises multivariadas. Além disso, a matriz de covariância é diagonalizada pelos autovetores, o que simplifica bastante os cálculos em algumas aplicações, como na Análise de Componentes Principais (PCA).
Em resumo, a extração de autovalores e autovetores da matriz de covariância nos permite identificar as direções principais ao longo das quais os dados têm maior e menor variância. Isso é fundamental em análises multivariadas, como a Análise de Componentes Principais (PCA), onde podemos reduzir a dimensionalidade dos dados projetando-os nas direções de maior variância para obter uma representação mais compacta dos dados.
Uma curva de nível de densidade constante para uma distribuição normal bivariada com \(\sigma_{11} = \sigma_{22}\) é obtido no exemplo a seguir.
Vamos obter os eixos das curvas de nível de densidade de probabilidade constante para uma distribuição normal bivariada quando \(\sigma_{11} = \sigma_{22}\). A partir de (4-7), esses eixos são dados pelos autovalores e autovetores de \(\mathbf{\Sigma}\).
\[ \begin{align} |\mathbf{\Sigma} - \lambda \mathbf{I}| &= 0\\ \begin{vmatrix} \sigma_{11}-\lambda &\sigma_{12}\\ \sigma_{12} & \sigma_{11}-\lambda \end{vmatrix}&=0\\ (\lambda-\sigma_{11}-\sigma_{12})(\lambda-\sigma_{11}+\sigma_{12}) &=0 \end{align} \]
Consequentemente, os autovalores são \(\lambda_1 = \sigma_{11} + \sigma_{12}=\sigma_{11}(1+\rho_{12})\) e \(\lambda_2 = \sigma_{11} - \sigma_{12}=\sigma_{11}(1-\rho_{12})\).
O autovetor \(\mathbf{e}_1\) é determinado por:
\[ \begin{align} \begin{bmatrix} \sigma_{11} & \sigma_{12}\\ \sigma_{12} & \sigma_{11}\\ \end{bmatrix} \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix} &= \lambda_1 \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix}\\ \begin{bmatrix} \sigma_{11} & \sigma_{12}\\ \sigma_{12} & \sigma_{11}\\ \end{bmatrix} \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix} &= (\sigma_{11}+\sigma_{12}) \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix} \end{align} \]
\[ \sigma_{11}e_1+\sigma_{12}e_2=(\sigma_{11}+\sigma_{12})e_1\\ \sigma_{12}e_1+\sigma_{11}e_2=(\sigma_{11}+\sigma_{12})e_2 \]
Essas equações implicam que \(e_1 = e_2\), e após normalização, o primeiro par autovalor-autovetor é \(\lambda_1 = \sigma_{11} + \sigma_{12}\) e \(\mathbf{e}_{1}^{\prime}= [1/\sqrt{2}, 1/\sqrt{2}]\).
\[ \begin{align} \begin{bmatrix} \sigma_{11} & \sigma_{12}\\ \sigma_{12} & \sigma_{11}\\ \end{bmatrix} \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix} &= \lambda_2 \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix}\\ \begin{bmatrix} \sigma_{11} & \sigma_{12}\\ \sigma_{12} & \sigma_{11}\\ \end{bmatrix} \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix} &= (\sigma_{11}-\sigma_{12}) \begin{bmatrix} e_1\\ e_2\\ \end{bmatrix} \end{align} \]
\[ \sigma_{11}e_1+\sigma_{12}e_2=(\sigma_{11}-\sigma_{12})e_1\\ \sigma_{12}e_1+\sigma_{11}e_2=(\sigma_{11}-\sigma_{12})e_2 \]
Essas equações implicam que \(e_2 = -e_1\), e após normalização, o primeiro par autovalor-autovetor é \(\lambda_2 = \sigma_{11} - \sigma_{12}\) e \(\mathbf{e}_{2}^{\prime}= [1/\sqrt{2}, -1/\sqrt{2}]\).
Quando a covariância \(\sigma_{12}\) (ou correlação \(\rho_{12}\)) é positiva, \(\lambda_1 = \sigma_{11} + \sigma_{12}\) é o maior autovalor, e seu autovetor associado \(\mathbf{e}_1^{\prime} = [1/\sqrt{2}, 1/\sqrt{2}]\) está ao longo da linha de 45° através do ponto \(\boldsymbol{\mu}' = [\mu_1, \mu_2]\). Isso é verdade para qualquer valor positivo da covariância (correlação).
Como os eixos das elipses de densidade constante são dados por \(\pm c\sqrt{\lambda_1}\mathbf{e}_1\) e \(\pm c\sqrt{\lambda_2}\mathbf{e}_2\) [veja (4-7)], e os autovetores cada um têm comprimento unitário, o eixo maior estará associado ao maior autovalor.
Para variáveis aleatórias normais positivamente correlacionadas, então, o eixo maior das elipses de densidade constante estará ao longo da linha de 45° através de \(\boldsymbol{\mu}^{\prime}\). (Veja a Figura 4.3.)
Figura 4.3 Uma curva de nível de densidade constante para uma distribuição normal bivariada.
Quando a covariância (ou correlação) é negativa, \(\lambda_2 = \sigma_{11} - \sigma_{12}\) será o maior autovalor, e os principais eixos das elipses de densidade constante estarão alinhados com uma linha em ângulo reto em relação à linha de 45° passando por \(\boldsymbol{\mu}\). (Esses resultados são válidos apenas para \(\sigma_{11}= \sigma_{22}\).)
Em resumo, os eixos das elipses de densidade constante para uma distribuição normal bivariada com \(\sigma_{11}= \sigma_{22}\) são determinados por:
\[ \pm c\sqrt{\sigma_{11}(1+\rho_{12})} \begin{bmatrix} \dfrac{1}{\sqrt{2}}\\ \dfrac{1}{\sqrt{2}} \end{bmatrix}\\ \pm c\sqrt{\sigma_{11}(1-\rho_{12})} \begin{bmatrix} \dfrac{1}{\sqrt{2}}\\ \dfrac{-1}{\sqrt{2}} \end{bmatrix} \]
\[\Diamond\]
Mostramos no Resultado 4.7 que a escolha \(c^2 = \chi_{p}^{2}(\gamma)\), em que \(\chi_{p}^{2}(\gamma)\) é o percentil superior de ordem \(100\gamma\) de uma distribuição qui-quadrado com \(p\) graus de liberdade, leva a curva de nível que contêm \(\gamma \times 100\%\) da probabilidade. Especificamente, o seguinte é verdadeiro para uma distribuição normal \(p\)-dimensional:
O elipsoide sólido de valores \(\mathbf{x}\) que satisfazem \((\mathbf{x} - \boldsymbol{\mu})^{\prime} \mathbf{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \leq \chi_{p}^{2}(\gamma)\) (4-8) tem probabilidade \(\gamma\).
As regiões elípticas que contêm 50% e 90% da probabilidade sob as superfícies normais bivariadas na Figura 4.2 são mostradas na Figura 4.4.
Figura 4-4 As regiões elípticas de 50% e 90% para as distribuições normais bivariadas na Figura 4.2.
# Quantil da distribuição qui-quadrado com 2 graus de liberdade para 50% e 90% de significância
q_2_50 <- qchisq(0.5, df = 2)
q_2_90 <- qchisq(0.9, df = 2)
# Quantil da distribuição qui-quadrado com 3 graus de liberdade para 50% e 95% de significância
q_3_50 <- qchisq(0.5, df = 3)
q_3_90 <- qchisq(0.9, df = 3)
# Quantil da distribuição qui-quadrado com 3 graus de liberdade para 50% e 95% de significância
q_10_50 <- qchisq(0.5, df = 10)
q_10_90 <- qchisq(0.9, df = 10)
tabela_quantis <- data.frame(Graus_de_Liberdade = c(2, 3, 10),
Quantil_50 = c(q_2_50, q_3_50, q_10_50),
Quantil_90 = c(q_2_90, q_3_90, q_10_90))
knitr::kable(round(tabela_quantis,2))
Graus_de_Liberdade | Quantil_50 | Quantil_90 |
---|---|---|
2 | 1.39 | 4.61 |
3 | 2.37 | 6.25 |
10 | 9.34 | 15.99 |
A densidade normal \(p\)-variada em (4-4) possui um valor máximo quando a distância ao quadrado em (4-3) é zero - isto é, quando \(\mathbf{x} = \boldsymbol{\mu}\). Assim, \(\boldsymbol{\mu}\) é o ponto de máxima densidade, ou moda, bem como o valor esperado de \(\mathbf{X}\), ou média. O fato de que \(\boldsymbol{\mu}\) é a média da distribuição normal multivariada segue da simetria exibida pelas curvas de nível de densidade constante: essas curvas são centradas, ou equilibradas, em \(\boldsymbol{\mu}\).
Certas propriedades da distribuição normal serão necessárias repetidamente em nossas explicações de modelos e métodos estatísticos. Essas propriedades tornam possível manipular facilmente as distribuições normais e, como sugerimos na Seção 4.1, são em parte responsáveis pela popularidade da distribuição normal. As principais propriedades, que discutiremos em detalhes matemáticos, podem ser declaradas de forma bastante simples.
As seguintes afirmações são verdadeiras para um vetor aleatório \(\mathbf{X}\) com distribuição normal multivariada:
Essas declarações são reproduzidas matematicamente nos resultados que seguem. Muitos desses resultados são ilustrados com exemplos. As provas incluídas devem ajudar a melhorar sua compreensão das manipulações de matrizes e também levá-lo a apreciar a maneira como os resultados se constroem sucessivamente.
O Resultado 4.2 pode ser considerado como uma definição operacional da distribuição normal. Com isso em mãos, as propriedades subsequentes são quase imediatas. Nossa prova parcial do Resultado 4.2 indica como a definição de densidade normal por combinação linear se relaciona com a densidade multivariada em (4-4).
Resultado 4.2. Se \(\mathbf{X}\) tem uma distribuição \(\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\), então qualquer combinação linear das variáveis \(\mathbf{a}^{\prime}\mathbf{X} = a_1X_1 + a_2X_2 + \cdots + a_pX_p\) é distribuída como \(\mathcal{N}(\mathbf{a}^{\prime}\boldsymbol{\mu}, \mathbf{a}^{\prime}\mathbf{\Sigma}\mathbf{a})\). Além disso, se \(\mathbf{a}^{\prime}\mathbf{X}\) tem uma distribuição \(\mathcal{N}(\mathbf{a}^{\prime}\boldsymbol{\mu}, \mathbf{a}^{\prime}\mathbf{\Sigma}\mathbf{a})\) para cada \(\mathbf{a}\), então \(\mathbf{X}\) deve ter distribuição \(\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\).
Resultado 4.3. Se \(\mathbf{X}\sim\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\), então a transformação linear de \(\mathbf{X}\):
\[ \underset{q\times p}{\mathbf{A}}\, \underset{p\times 1}{\mathbf{X}}+ \underset{q\times 1}{\mathbf{b}} \sim \mathcal{N}_q\left(\mathbf{A}\boldsymbol{\mu}+\mathbf{b}, \mathbf{A}\mathbf{\Sigma}\mathbf{A}^{\prime}\right) \]
Em que:
\[ \mathbf{A}\mathbf{X} +\mathbf{b}= \begin{bmatrix} \sum_{i=1}^{p}{a_{1i}X_i + b_1} \\ \sum_{i=1}^{p}{a_{2i}X_i + b_2} \\ \vdots \\ \sum_{i=1}^{p}{a_{qi}X_i + b_q} \end{bmatrix} \]
Resultado 4.5.
Se \(\underset{q_1\times 1}{{\mathbf{X}_{1}}}\) e \(\underset{q_2\times 1}{{\mathbf{X}_{2}}}\) são independentes, então \(\mathbb{C}(\mathbf{X}_1, \mathbf{X}_2) = \mathbf{0}\), ou seja, uma matriz \(q_1 \times q_2\) de zeros.
Se \(\left[\mathbf{X}_1\;\mathbf{X}_2\right]^{\prime}\sim \mathcal{N}_{q_1+q_2}\left(\begin{bmatrix} \boldsymbol{\mu}_1\\\boldsymbol{\mu}_2\end{bmatrix},\begin{bmatrix} \mathbf{\Sigma}_{11} & \mathbf{\Sigma}_{12} \\ \mathbf{\Sigma}_{21} & \mathbf{\Sigma}_{22} \end{bmatrix}\right)\), sendo que \(\mathbf{\Sigma}_{12}=\mathbf{\Sigma}_{21}^{\prime}\), então \(\mathbf{X}_1\) e \(\mathbf{X}_2\) são independentes se e somente se \(\mathbf{\Sigma}_{12} = \mathbf{0}\).
Seja \(\mathbf{X}_{3}\sim\mathcal{N}_3(\boldsymbol{\mu}, \mathbf{\Sigma})\).
\[ \mathbf{\Sigma}= \begin{bmatrix} 4 & 1 & 0 \\ 1 & 3 & 0 \\ 0 & 0 & 2 \end{bmatrix} \]
São \(X_1\) e \(X_2\) independentes? E quanto a \((X_1, X_2)\) e \(X_3\)?
Uma vez que \(X_1\) e \(X_2\) têm covariância \(\sigma_{12} = 1\), eles não são independentes.
Portanto, \((X_1 , X_2)\) e \(X_3\) são independentes de acordo com o Resultado 4.5. Isso implica que \(X_3\) é independente de \(X_1\) e também de \(X_2\).
\[\Diamond\]
Destacamos em nossa discussão da distribuição normal bivariada que \(\rho_{12} = 0\) (correlação nula) implica independência porque a função de densidade conjunta [ver (4-6)] pode então ser escrita como o produto das densidades marginais (normais) de \(X_1\) e \(X_2\). Esse fato, que encorajamos você a verificar diretamente, é simplesmente um caso especial do Resultado 4.5 com \(q_1 = q_2 = 1\).
\[ \begin{align} f(x_1,x_2) &= \dfrac{1}{2 \pi \sqrt{\sigma_{11} \sigma_{22}(1 - \rho_{12}^2)}} \exp \left( -\dfrac{1}{2(1 - \rho_{12}^2)} \left( z_{1}^{2}+z_{1}^{2}-2\rho_{12}z_{1}z_{1} \right) \right) \\ &= \dfrac{1}{2 \pi \sqrt{\sigma_{11} \sigma_{22}}} \exp \left( -\dfrac{1}{2} \left( z_{1}^{2}+z_{1}^{2} \right) \right)\\ &= \dfrac{1}{\sqrt{2 \pi \sigma_{11}}} \exp \left( -\dfrac{1}{2} z_{1}^{2} \right) \dfrac{1}{\sqrt{2 \pi \sigma_{22}}} \exp \left( -\dfrac{1}{2} z_{2}^{2} \right)\\ f(x_1,x_2) &= f(x_1)f(x_2)\\ \end{align} \tag{4-6} \]
Seja \(\mathbf{Z}^{\prime}=[Z_1\;Z_2\;\cdots\;Z_p]\) e \(\{Z_i\}_{i=1}^{p}\sim \mathcal{N}\text{IID}(0, 1)\).
A normal multivariada padrão pode ser definida como:
\[ \mathbf{Z}\sim \mathcal{N}_p(\mathbf{0}, \mathbf{I}) \]
A matriz de covariância populacional pode ser equifatorada:
\[ \mathbf{\Sigma}=\mathbf{\Sigma}^{1/2}\left(\mathbf{\Sigma}^{1/2}\right)^{\prime} \]
Em que \(\mathbf{\Sigma}^{1/2}=\mathbf{P\Lambda}^{1/2}\)
é o espectro da matriz de covariância populacional de tal modo que \(\mathbf{\Sigma}=\mathbf{P\Lambda
P}^{\prime}\) ou \(\mathbf{\Sigma}^{1/2}=\mathbf{C}\) é a
matriz triangular de Cholesky de tal modo que \(\mathbf{\Sigma}=\mathbf{CC}^{\prime}\).
Genz, A. em 1992 publicou o artigo seminal usando matriz triangular de
Cholesky. Ele também é um dos autores do pacote mvtnorm
. A
matriz diagonal de autovalores e a matriz de autovetores ortonormais por
coluna é denotada, respectivamente, por \(\mathbf{\Lambda}\) e \(\mathbf{P}\).
Então:
\[ \mathbf{X}=\mathbf{\Sigma}^{1/2}\mathbf{Z}+\boldsymbol{\mu}\sim \mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma}) \]
Este resultado pode ser usado para simular a distribuição normal multivariada a partir de \(p\) normais padrão univariadas independentes.
[,1] [,2]
[1,] 1.0 0.0000000
[2,] 0.5 0.8660254
[,1] [,2]
[1,] 1.0 0.5
[2,] 0.5 1.0
rMvNorm <- function(mu,Sigma){
r <- length(mu)
C <- t(chol(Sigma))
Z <- rnorm(r)
return(C %*% Z + mu)
}
B <- 5e3
X <- matrix(0, nrow=2, ncol=B)
for(i in 1:B) X[,i] <- rMvNorm(mu, Sigma)
plot(X[1,],X[2,],
main="binormal simulada",
asp=1)
result <- MVN::mvn(data=data.frame(X[1,],X[2,]),
univariate_test="SW")
print(result$multivariate_normality)
Test Statistic p.value Method MVN
1 Henze-Zirkler 0.721 0.752 asymptotic ✓ Normal
Test Variable Statistic p.value Normality
1 Shapiro-Wilk X.1... 1 0.953 ✓ Normal
2 Shapiro-Wilk X.2... 1 0.355 ✓ Normal
Uma aplicação de binormal simulada é apresentada a seguir.
# Simulacao de distribuicao binormal de estatura e MCT
# Estudante masculino do curso de
# Medicina da FMUSP de 2015-18
n <- 5e3
m.estat.m <- 177; m.mct.m <- 73
s.estat.m <- 7; s.mct.m <- 9
r <- 0.6
cov <- r*s.estat.m*s.mct.m
# Parameters for bivariate normal distribution
mv <- c(m.estat.m, m.mct.m) # Mean vector
covar <- matrix(c(s.estat.m^2, cov,
cov, s.mct.m^2), 2) # Covariance matrix
set.seed(123)
estat.mct.m <- MASS::mvrnorm(n, mu = mv, Sigma = covar)
dados <- as.data.frame(estat.mct.m)
names(dados) <- c("Estaturam", "MCTm")
minx <- min(dados$Estaturam)*0.99
maxx <- max(dados$Estaturam)*1.02
miny <- min(dados$MCTm)*0.99
maxy <- max(dados$MCTm)*1.02
plot(dados,
main=paste0("Dados simulados sob binormalidade",
", n(masculino) = ",n,
", rho = ",r),
xlab="Estatura (cm)",
xlim=c(minx,maxx),
ylab="MCT (kg)",
ylim=c(miny,maxy))
# car::dataEllipse(estat.mct.m[,1], estat.mct.m[,2], # suspeito...
# levels=c(.5, .95, .999),
# robust=TRUE, grid=FALSE,
# lwd=2, lty=1,
# col=c("black"),
# cex=0.1,
# main="Elipse de predição de 50, 95 e 99.9% de binormal",
# xlab="Estatura (cm)",
# xlim=c(minx,maxx),
# ylab="MCT (kg)",
# ylim=c(miny,maxy))
car::scatterplot(estat.mct.m[,1], estat.mct.m[,2],
regLine=FALSE,
ellipse=list(levels=c(.5, .95, .999),
robust=TRUE, fill=FALSE),
grid=FALSE,
smooth=FALSE,
col=c("black"),
cex=0.1,
pch=1,
lwd=2, lty=1,
main="Elipse de predição de 50, 95 e 99.9% de binormal",
xlab="Estatura (cm)",
xlim=c(minx,maxx),
ylab="MCT (kg)",
ylim=c(miny,maxy))
# Bag plot
bgp <- DescTools::PlotBag(dados$Estaturam,
dados$MCTm,
main="Estudante Masculino de Medicina - FMUSP (simulado)",
xlab="Estatura (cm)",
xlim=c(minx,maxx),
ylab="MCT (kg)",
ylim=c(miny,maxy),
na.rm=TRUE,
show.bagpoints=FALSE,
show.looppoints=FALSE,
show.whiskers=FALSE,
col.loophull="white",
col.looppoints="black",
col.baghull="white",
col.bagpoints="black",
cex=1)
print(outliers <- as.data.frame(bgp$pxy.outlier),
digits = 2)
x y
1 164 86
2 185 102
3 199 75
4 181 45
5 157 77
6 149 53
7 183 52
8 194 104
for (o in 1:nrow(outliers))
{
r <- which(dados$Estaturam==outliers$x[o] &
dados$MCTm==outliers$y[o])
text(outliers$x[o],outliers$y[o], r, pos=1, cex=0.7)
}
Resultado 4.6. Seja \(\mathbf{X}^{\prime}=\left[\mathbf{X}_1\;\mathbf{X}_2\right]\sim \mathcal{N}_{p}\left(\begin{bmatrix} \boldsymbol{\mu}_1\\\boldsymbol{\mu}_2\end{bmatrix},\begin{bmatrix} \mathbf{\Sigma}_{11} & \mathbf{\Sigma}_{12} \\ \mathbf{\Sigma}_{21} & \mathbf{\Sigma}_{22} \end{bmatrix}\right)\), sendo que \(\mathbf{\Sigma}_{21}=\mathbf{\Sigma}_{12}^{\prime}\), então:
\[ \begin{align} \mathbf{X}_1 |\mathbf{X}_2 &= \mathbf{x}_2 \sim \mathcal{N}_{p}\left(\boldsymbol{\mu}_{1|2} , \mathbf{\Sigma}_{1|2}\right)\\ \mathbb{E}\left(\mathbf{X}_1 |\mathbf{X}_2=\mathbf{x}_2\right) &= \boldsymbol{\mu}_{1|2}=\boldsymbol{\mu}_1+\mathbf{\Sigma}_{12}\mathbf{\Sigma}_{22}^{-1}\left(\mathbf{x}_2-\boldsymbol{\mu}_2\right)\\ \mathbb{C}\left(\mathbf{X}_1 |\mathbf{X}_2=\mathbf{x}_2\right)&=\mathbf{\Sigma}_{1|2}=\mathbf{\Sigma}_{11}-\mathbf{\Sigma}_{12}\mathbf{\Sigma}_{22}^{-1}\mathbf{\Sigma}_{12}^{\prime} \end{align} \]
Sendo que \(\boldsymbol{\beta}_{1,2}=\mathbf{\Sigma}_{12}\mathbf{\Sigma}_{22}^{-1}\).
Então:
\[ \begin{align} \boldsymbol{\mu}_{1|2}&=\boldsymbol{\mu}_1+\boldsymbol{\beta}_{1,2}\left(\mathbf{x}_2-\boldsymbol{\mu}_2\right)\\ \mathbf{\Sigma}_{1|2}&=\mathbf{\Sigma}_{11}-\boldsymbol{\beta}_{1,2}\mathbf{\Sigma}_{12}^{\prime} \end{align} \]
Note que a esperança condicional é uma transformação linear de \(\mathbf{x}_2\) e que a covariância condicional não depende de \(\mathbf{x}_2\), i.e., a distribuição condicional é homocedástica.
Conforme Bickel & Docksum (2001, p. 521), há ganho informacional com o condicionamento, i.e., \(\mathbf{\Sigma}_{1|2}<\mathbf{\Sigma}_{11}\), pois \(\mathbf{\Sigma}_{11}>\mathbf{\Sigma}_{12}\mathbf{\Sigma}_{22}^{-1}\mathbf{\Sigma}_{21}>0\) (desigualdade de Cauchy-Schwarz generalizada).
A densidade condicional de \(X_1\) dado que \(X_2 = x_2\), para qualquer distribuição bivariada, é definida por:
\[ f(x_1|x_2) = \dfrac{f(x_1, x_2)}{f_2(x_2)} \]
em que \(f(x_2)\) é a distribuição marginal de \(X_2\).
Se \(f(x_1, x_2)\) é a densidade binormal
\[ X_1|X_2=x_2\sim \mathcal{N}\left(\mu_1+\dfrac{\sigma_{12}}{\sigma_{22}}(x_2-\mu_2),\;\sigma_{11}-\dfrac{\sigma_{12}^{2}}{\sigma_{22}}\right)\\ X_1|X_2=x_2\sim \mathcal{N}\left(\mu_1+\beta_{1,2}(x_2-\mu_2),\;\sigma_{11}(1-\rho_{12}^{2})\right) \]
\[ \mathbf{X}\sim \mathcal{N}_4(\boldsymbol{\mu}, \mathbf{\Sigma}) \]
Se há duas VD (\(X_1, X_2\)) e duas VE (\(X_3, X_4\)), então \(p=4\) e \(q=2\):
\[ \begin{align} \mathbf{X}^{\prime}&=[\mathbf{X}_1^{\prime}\; \mathbf{X}_2^{\prime}]=[X_1\;X_2\;X_3\;X_4]\\ \mathbf{X}_1^{\prime}&=[X_1\;X_2]\\ \mathbf{X}_2^{\prime}&=[X_3\;X_4] \end{align} \]
Portanto:
\[ \mathbb{E}(\mathbf{X}_1|\mathbf{X}_2=\mathbf{x}_2)= \begin{bmatrix} \mu_1+\beta_{1,3}(x_3-\mu_3)+\beta_{1,4}(x_4-\mu_4)\\ \mu_2+\beta_{2,3}(x_3-\mu_3)+\beta_{2,4}(x_4-\mu_4) \end{bmatrix} \]
Em que \(\beta_{i,j}=\rho_{i,j}\dfrac{\sigma_{i}}{\sigma_{j}}\).
Resultado 4.7. Se \(\mathbf{X}\sim\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\), então:
\((\mathbf{X} - \boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1}(\mathbf{X} - \boldsymbol{\mu})\sim\chi_{p}^{2}\), i.e., segue uma distribuição qui-quadrado com \(p\) graus de liberdade.
A distribuição \(\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\) atribui probabilidade \(1 - \alpha\) ao hiperelipsoide sólido \(\{\mathbf{x} : (\mathbf{x} - \boldsymbol{\mu})'\mathbf{\Sigma}^{-1}(\mathbf{x} - \boldsymbol{\mu}) \leq \chi^2_p(1 - \alpha)\}\), em que \(\chi^2_p(1 - \alpha)\) denota o quantil superior \(100(1 - \alpha)\) da distribuição qui-quadrado com \(p\) graus de liberdade.
\[ P\left(\left(\mathbf{X} - \boldsymbol{\mu}\right)^{\prime}\mathbf{\Sigma}^{-1}\left(\mathbf{X} - \boldsymbol{\mu}\right) \leq \chi^2_p(1 - \alpha)\right) = 1 - \alpha \]
em que \(\chi^2_p(1 - \alpha)\) denota o quantil superior \(1 - \alpha\) da distribuição qui-quadrado com \(p\) graus de liberdade.
O Resultado 4.7 fornece uma interpretação da distância euclidiana estatística (Mahalanobis) ao quadrado. Se \(\mathbf{X}\sim\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\),
\[ (\mathbf{X} - \boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1}(\mathbf{X} - \boldsymbol{\mu}) \]
é a distância estatística ao quadrado de \(\mathbf{X}\) ao vetor de média populacional \(\boldsymbol{\mu}\). Se um componente tiver uma variância muito maior do que outra, ele contribuirá menos para a distância estatística ao quadrado. Além disso, duas variáveis aleatórias altamente correlacionadas contribuirão menos do que duas variáveis que são quase não correlacionadas. Essencialmente, o uso do inverso da matriz de covariância (1) padroniza todas as variáveis e (2) elimina os efeitos da correlação. A partir da prova do Resultado 4.7,
\[ (\mathbf{X} - \boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1}(\mathbf{X} - \boldsymbol{\mu}) = \sum_{i=1}^{p}{Z_i^2} \sim \chi_p^2 \]
em que \(Z_i=\dfrac{\mathbf{e}_{i}^{\prime}}{\sqrt{\lambda_i}}(\mathbf{X}-\boldsymbol{\mu}), \;i=1, 2, \ldots, p\) são variáveis aleatórias independentes e padronizadas (note que os autovalores são ortogonais!), cada uma seguindo uma distribuição qui-quadrado com 1 grau de liberdade, i.e., \(\chi_1^2\). A soma de \(p\) variáveis qui-quadrado independnetes com 1 grau de liberdade tem distribuição qui-quadrado com \(p\) graus de liberdade, i.e., \(\chi_p^2\). Isso significa que a distância estatística ao quadrado pode ser entendida como uma soma ponderada das contribuições das variáveis padronizadas e independentes, e é proporcional à distância quadrática média de cada variável à sua média populacional.
Em termos de \(\mathbf{\Sigma}^{-1/2}\) (veja (2-22)), a variável aleatória \(\mathbf{\Sigma}^{-1/2}(\mathbf{X} - \boldsymbol{\mu})\sim \mathcal{N}_p(\mathbf{0}, \mathbf{I})\), em que \(\mathbf{I}\) é a matriz identidade \(p \times p\). Além disso,
\[ \begin{align} (\mathbf{X} - \boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1}(\mathbf{X} - \boldsymbol{\mu}) &= (\mathbf{X} - \boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1/2}\mathbf{\Sigma}^{-1/2}(\mathbf{X} - \boldsymbol{\mu})\\ &= \mathbf{Z}'\mathbf{Z} \\ (\mathbf{X} - \boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1}(\mathbf{X} - \boldsymbol{\mu})&= \sum_{i=1}^{p}{Z_i^2} \end{align} \]
A distância estatística ao quadrado é calculada como se, inicialmente, o vetor aleatório \(\mathbf{X}\) fosse transformado em \(p\) variáveis aleatórias normais padrão independentes e, em seguida, a distância estatística ao quadrado usual, a soma dos quadrados das variáveis, fosse aplicada. Em outras palavras, a distância estatística ao quadrado é equivalente à soma ponderada dos quadrados das variáveis aleatórias independentes e padronizadas.
A seguir, considere a combinação linear de variáveis aleatórias vetoriais (observações multivariadas independentes):
\[ \sum_{i=1}^{n}{c_i\mathbf{X}_i} = \begin{bmatrix} \mathbf{X}_1 & \mathbf{X}_2 & \cdots & \mathbf{X}_n \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{bmatrix} \tag{4-10} \]
Resultado 4.8. Se \(\mathbf{X}_1, \mathbf{X}_2, \ldots, \mathbf{X}_n\) forem mutuamente independentes, com \(\mathbf{X}_i\sim\mathcal{N}_p(\boldsymbol{\mu}_i, \mathbf{\Sigma})\), em que cada \(\mathbf{X}_i\) tem a mesma matriz de covariância \(\mathbf{\Sigma}\), então:
\[ \mathbf{V_1} = \sum_{i=1}^{n}{c_i\mathbf{X}_i}\sim\mathcal{N}_p(\boldsymbol{\mu}_V, \mathbf{\Sigma}_V) \]
em que \(\boldsymbol{\mu}_{\mathbf{V}_1} = \sum_{i=1}^{n}{c_i\boldsymbol{\mu}_i}\) e \(\mathbf{\Sigma}_{\mathbf{V}_1} = (\mathbf{c}^{\prime}\mathbf{c})\mathbf{\Sigma}\). Além disso, \(\mathbf{V}_1\) e \(\mathbf{V}_2 = \sum_{i=1}^{n}{b_i\mathbf{X}_i}\) são conjuntamente normais multivariadas:
\[ \mathbf{V}= \begin{bmatrix} \mathbf{V}_1 \\\mathbf{V}_2 \end{bmatrix} \sim\mathcal{N}_{2p}\left( \begin{bmatrix} \boldsymbol{\mu}^{\prime}\mathbf{c}\\ \boldsymbol{\mu}^{\prime}\mathbf{b} \end{bmatrix}, \mathbf{\Sigma}_{\mathbf{V}_1,\mathbf{V}_2}\right) \]
em que:
\[ \boldsymbol{\mu}^{\prime} = [\boldsymbol{\mu}_1\;\boldsymbol{\mu}_2\;\cdots\;\boldsymbol{\mu}_n] \]
\[ \mathbf{\Sigma}_{\mathbf{V}_1,\mathbf{V}_2} = \begin{bmatrix} (\mathbf{c}^{\prime}\mathbf{c})\mathbf{\Sigma} & (\mathbf{b}^{\prime}\mathbf{c})\mathbf{\Sigma}\\ (\mathbf{b}^{\prime}\mathbf{c})\mathbf{\Sigma} & (\mathbf{b}^{\prime}\mathbf{b})\mathbf{\Sigma} \end{bmatrix} = \begin{bmatrix} \mathbf{c}^{\prime}\mathbf{c} & \mathbf{b}^{\prime}\mathbf{c}\\ \mathbf{b}^{\prime}\mathbf{c} & \mathbf{b}^{\prime}\mathbf{b} \end{bmatrix}\otimes\mathbf{\Sigma} \]
\[ \mathbf{V}= \begin{bmatrix} \mathbf{V}_1 \\\mathbf{V}_2 \end{bmatrix} \sim\mathcal{N}_{2p}\left( \begin{bmatrix} \boldsymbol{\mu}^{\prime}\mathbf{c}\\ \boldsymbol{\mu}^{\prime}\mathbf{b} \end{bmatrix}, \begin{bmatrix} \mathbf{c}^{\prime}\mathbf{c} & \mathbf{b}^{\prime}\mathbf{c}\\ \mathbf{b}^{\prime}\mathbf{c} & \mathbf{b}^{\prime}\mathbf{b} \end{bmatrix}\otimes\mathbf{\Sigma}\right) \]
Consequentemente, \(\mathbf{V}_1\) e \(\mathbf{V}_2\) são independentes se \(\mathbf{b}^{\prime}\mathbf{c} = 0\).
Se \(c_i=\dfrac{1}{n}\), então \(\mathbf{V}_1=\bar{\mathbf{X}}_1=\dfrac{1}{n}\sum_{i=1}^{n}{\mathbf{X}_i}\sim \mathcal{N}_p\left(\dfrac{1}{n}\sum_{i=1}^{n}{\boldsymbol{\mu}_i},\dfrac{1}{n}\mathbf{\Sigma}\right)\)
O símbolo do produto de Kroneker é \(\otimes\).
O produto de Kronecker, também conhecido como produto tensorial, é uma operação entre duas matrizes que generaliza o produto de matrizes convencional. Se você tem duas matrizes \(A\) e \(B\) de dimensões \(m \times n\) e \(p \times q\), respectivamente, o produto de Kronecker \(A \otimes B\) será uma matriz de dimensões \(mp \times nq\).
A operação é definida de forma que cada elemento \(a_{ij}\) da matriz \(A\) é multiplicado pela matriz completa \(B\), e o resultado é colocado na posição correspondente na matriz resultante. Matematicamente, o produto de Kronecker \(A \otimes B\) é dado por:
\[ A \otimes B = \begin{pmatrix} a_{11} B & a_{12} B & \cdots & a_{1n} B \\ a_{21} B & a_{22} B & \cdots & a_{2n} B \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} B & a_{m2} B & \cdots & a_{mn} B \end{pmatrix} \]
Por exemplo, se \(A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}\) e \(B = \begin{pmatrix} w & x \\ y & z \end{pmatrix}\), então o produto de Kronecker \(A \otimes B\) será:
\[ A \otimes B = \begin{pmatrix} a \times B & b \times B \\ c \times B & d \times B \end{pmatrix} = \begin{pmatrix} aw & ax & bw & bx \\ ay & az & by & bz \\ cw & cx & dw & dx \\ cy & cz & dy & dz \end{pmatrix} \]
Discutimos brevemente amostragem e seleção de amostras aleatórias no Capítulo 3. Nesta seção, estaremos preocupados com amostras de uma população normal multivariada, em particular, com a distribuição amostral de \(\bar{\mathbf{X}}\) e \(\mathbf{S}\).
Vamos assumir que as \(n\) observações multivariadas \(p \times 1\), \(\mathbf{X}_1, \mathbf{X}_2, \ldots, \mathbf{X}_n\), representam uma amostra aleatória de uma população normal multivariada com vetor média \(\boldsymbol{\mu}\) e matriz de covariância \(\mathbf{\Sigma}\), i.e., \(\left\{ \mathbf{X}_i \right\}_{i=1}^{n} \sim \mathcal{N}\text{IID}(\boldsymbol{\mu}, \mathbf{\Sigma})\). Uma vez que eles são mutuamente independentes e cada um tem uma distribuição \(\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma})\), a função de densidade conjunta de todas as observações ou função de verossimilhança é o produtório das densidades marginais normais:
\[ \begin{align} \mathcal{L}(\boldsymbol{\mu}, \mathbf{\Sigma}) &= \prod_{i=1}^{n}{f(\mathbf{x}_i)}\\ \mathcal{L}(\boldsymbol{\mu}, \mathbf{\Sigma})&= |2\pi\mathbf{\Sigma}|^{-n/2}\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathbf{x}_i-\boldsymbol{\mu})^{\prime}\mathbf{\Sigma}^{-1}(\mathbf{x}_i-\boldsymbol{\mu})\right) \end{align} \tag{4-11} \]
Quando os valores numéricos das observações se tornam disponíveis, eles podem ser substituídos pelos \(\mathbf{x}_i\) na Equação (4-11). A expressão resultante, agora considerada como uma função de \(\boldsymbol{\mu}\) e \(\mathbf{\Sigma}\) para o conjunto fixo de observações multivariadas independentes \(\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n\), é chamada de verossimilhança. Muitos bons procedimentos estatísticos empregam valores para os parâmetros populacionais que “melhor” explicam os dados observados. Um significado de “melhor” é selecionar os valores dos parâmetros que maximizam a densidade conjunta avaliada nas observações. Essa técnica é chamada de estimação de máxima verossimilhança (Maximum Likelihood Estimation - MLE), e os valores dos parâmetros que maximizam são chamados de estimativas de máxima verossimilhança.
Neste ponto, consideraremos a estimação de máxima verossimilhança dos parâmetros \(\boldsymbol{\mu}\) e \(\mathbf{\Sigma}\) para uma população normal multivariada. Para fazer isso, tomamos as observações (valores numéricos) multivariadas independentes \(\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n\) como fixas e consideramos a densidade conjunta da Equação (4-11) avaliada nesses valores numéricos. O resultado é a função de verossimilhança. Para simplificar, reescrevemos a função de verossimilhança em outra forma. Precisaremos de algumas propriedades adicionais para o traço de uma matriz quadrada. (O traço de uma matriz é a soma dos elementos diagonais, e as propriedades do traço são discutidas na Definição 2A.28 e no Resultado 2A.12.)
Resultado 4.9. Seja \(\mathbf{\Sigma}\) uma matriz simétrica \(p \times p\) e \(\mathbf{x}\) um vetor \(p \times 1\).
\[ \begin{align} \mathbf{x}^{\prime}\mathbf{\Sigma}\mathbf{x} &= \mathrm{tr}(\mathbf{x}^{\prime}\mathbf{\Sigma}\mathbf{x}) \\ &= \mathrm{tr}(\mathbf{\Sigma}\mathbf{x}\mathbf{x}^{\prime}) \\ \mathbf{x}^{\prime}\mathbf{\Sigma}\mathbf{x} &= \mathrm{tr}(\mathbf{x}\mathbf{x}^{\prime}\mathbf{\Sigma}) \end{align} \]
A igualdade \(\mathbf{\Sigma} = \mathbf{P}^{\prime}\mathbf{\Lambda P}\), em que \(\mathbf{PP}^{\prime} = \mathbf{I}\) (\(\mathbf{P}\) é ortogonal) e \(\mathbf{\Lambda}\) é uma matriz diagonal com entradas \(\lambda_1, \lambda_2, \ldots, \lambda_p\), implica que a soma dos autovalores \(\lambda_i\) é igual à soma dos autovalores da matriz \(\mathbf{P}^{\prime}\mathbf{\Lambda}\mathbf{P}\). Vamos verificar isso passo a passo:
\[ \begin{align} \mathrm{tr}(\mathbf{\mathbf{\Sigma}}) &= \mathrm{tr}(\mathbf{P}^{\prime}\mathbf{\Lambda P}) \\ &= \mathrm{tr}(\mathbf{\Lambda PP}^{\prime}) \\ &= \mathrm{tr}(\mathbf{\Lambda}) \\ \mathrm{tr}(\mathbf{\mathbf{\Sigma}})&= \sum_{i=1}^{p}{\lambda_i} \end{align} \]
Portanto, a soma dos autovalores \(\lambda_i\) é igual à soma dos autovalores da matriz \(\mathbf{P}^{\prime}\mathbf{\Lambda P}\).
\[ \begin{align} \mathbb{E}(\mathbf{X}^{\prime}\mathbf{\Sigma}^{-1}\mathbf{X}) &=\mathbb{E}(\mathrm{tr}(\mathbf{X}^{\prime}\mathbf{\Sigma}^{-1}\mathbf{X}))\\ &= \mathrm{tr}(\mathbf{\Sigma}^{-1}\mathbf{\Sigma}) + \boldsymbol{\mu}^{\prime}\mathbf{\Sigma}^{-1}\boldsymbol{\mu}\\ &=\mathrm{tr}(\mathbf{I}) + \boldsymbol{\mu}^{\prime}\mathbf{\Sigma}^{-1}\boldsymbol{\mu}\\ \mathbb{E}(\mathbf{X}^{\prime}\mathbf{\Sigma}^{-1}\mathbf{X})&= p + \boldsymbol{\mu}^{\prime}\mathbf{\Sigma}^{-1}\boldsymbol{\mu} \end{align} \]
\[\Diamond\]
A nova função de verossimilhança é dada por:
\[ \mathcal{L}(\boldsymbol{\mu}, \mathbf{\Sigma}) = |2\pi\mathbf{\Sigma}|^{-n/2}\exp\left(-\frac{1}{2}\left(\;\mathrm{tr}\left(\left(\dfrac{\mathbf{\Sigma}}{n}\right)^{-1}\mathbf{s}_n\right) +\\(\bar{\mathbf{x}}-\boldsymbol{\mu})^{\prime}\left(\dfrac{\mathbf{\Sigma}}{n}\right)^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu})\right)\right) \tag{4-15} \]
A obtenção das estimativas de máxima verossimilhança (MLE) da densidade conjunta de uma amostra aleatória de uma população normal multivariada envolve derivar os parâmetros que maximizam a função de verossimilhança. A demonstração completa é complexa e requer conhecimento avançado em cálculo diferencial multivariado e otimização.
O objetivo é encontrar os valores de \(\boldsymbol{\mu}\) e \(\mathbf{\Sigma}\) que maximizam a função de verossimilhança \(\mathcal{L}(\boldsymbol{\mu}, \mathbf{\Sigma})\).
As estimativas de máxima verossimilhança resultantes são chamados de estimativas de máxima verossimilhança da média (\(\hat{\boldsymbol{\mu}}\)) e da matriz de covariância (\(\hat{\mathbf{\Sigma}}\)) da população normal multivariada.
Condição necessária de primeira ordem
Isso envolve derivar a função log-verossimilhança e resolver as equações de derivadas parciais para \(\boldsymbol{\mu}\) e \(\mathbf{\Sigma}\).
A função log-verossimilhança é dada por:
\[ \ln \mathcal{L}(\boldsymbol{\mu}, \mathbf{\Sigma}) = -\dfrac{np}{2} \ln(2\pi) - \dfrac{n}{2} \ln |\mathbf{\Sigma}| - \dfrac{1}{2} \mathrm{tr}\left(\left(\dfrac{\mathbf{\Sigma}}{n}\right)^{-1}\mathbf{s}_n\right) + \\(\bar{\mathbf{x}}-\boldsymbol{\mu})^{\prime}\left(\dfrac{\mathbf{\Sigma}}{n}\right)^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}) \]
As derivadas parciais em relação a \(\boldsymbol{\mu}\) e \(\mathbf{\Sigma}\) são então calculadas e igualadas a zero para encontrar o ponto crítico, conforme Reis (2001, p. 105-8):
\[ \begin{align} \dfrac{\partial \ln \mathcal{L}}{\partial \boldsymbol{\mu}} &= \mathbf{\Sigma}^{-1} (\bar{\mathbf{x}} - \boldsymbol{\mu})\\ \frac{\partial\ln\mathcal{L}}{\partial \mathbf{\Sigma}} &= \dfrac{n}{2}(2\mathbf{M}-\text{diag}(\mathbf{M}))\\ \mathbf{M}&=\mathbf{\Sigma}-\mathbf{s}_n-(\bar{\mathbf{x}} - \boldsymbol{\mu})(\bar{\mathbf{x}} - \boldsymbol{\mu})^{\prime} \end{align} \]
Igualando as duas derivadas parciais a \(\mathbf{0}\) para determinação do ponto crítico, tem-se:
\[ \begin{align} \hat{\mathbf{\Sigma}}^{-1} (\bar{\mathbf{x}} - \hat{\boldsymbol{\mu}})&=\mathbf{0}\\ \dfrac{n}{2}(2\hat{\mathbf{M}}-\text{diag}(\hat{\mathbf{M}}))&=\mathbf{0}\\ \hat{\mathbf{M}}&=\hat{\mathbf{\Sigma}}-\mathbf{s}_n-(\bar{\mathbf{x}} - \hat{\boldsymbol{\mu}})(\bar{\mathbf{x}} - \hat{\boldsymbol{\mu}})^{\prime} \end{align} \]
Resolvendo \(\hat{\mathbf{\Sigma}}^{-1} (\bar{\mathbf{x}} - \hat{\boldsymbol{\mu}})=\mathbf{0}\), obtemos a estimativa candidata à de máxima verossimilhança para \(\boldsymbol{\mu}\):
\[ \begin{align} \hat{\boldsymbol{\mu}} &= \dfrac{1}{n} \sum_{i=1}^{n} \mathbf{x}_i\\ \hat{\boldsymbol{\mu}} &=\bar{\mathbf{x}} \end{align} \]
Como a última parcela de \(\hat{\mathbf{M}}\) é nula, então \(\hat{\mathbf{M}}=\hat{\mathbf{\Sigma}}-\mathbf{s}_n\).
Para que \(\dfrac{n}{2}(2\hat{\mathbf{M}}-\text{diag}(\hat{\mathbf{M}}))=\mathbf{0}\) é necessário que \(\hat{\mathbf{M}}=\hat{\mathbf{\Sigma}}-\mathbf{s}_n=\mathbf{0}\).
Portanto, a estimativa candidata à de máxima verossimilhança de matriz de covariância populacional é:
\[ \hat{\mathbf{\Sigma}}=\mathbf{s}_n \]
Condição suficiente de segunda ordem
A demonstração foi realizada por Lucas Ayres B. de C. Barros no
arquivo intitulado MaxVerossMultivariada.pdf
.
Temos que o valor máximo da função de verossimilhança é:
\[ \begin{align} \mathcal{L}(\hat{\boldsymbol{\mu}},\hat{\mathbf{\Sigma}})&= |2\pi \hat{\mathbf{\Sigma}}|^{-n/2}\exp\left(-\dfrac{np}{2}\right) \\ &= (2\pi)^{-np/2}\exp\left(-\dfrac{np}{2}\right)|\hat{\mathbf{\Sigma}}|^{-n/2}\\ &= k|\hat{\mathbf{\Sigma}}|^{-n/2}\\ &= k\left(\dfrac{n-1}{n}\right)^{p}|\mathbf{s}|^{-n/2}\\ \mathcal{L}(\hat{\boldsymbol{\mu}},\hat{\mathbf{\Sigma}})&= c|\mathbf{s}|^{-n/2}\\ \mathcal{L}(\hat{\boldsymbol{\mu}},\hat{\mathbf{\Sigma}})&\propto \text{variância generalizada}^{-n/2} \end{align} \tag{4-18} \]
A variância generalizada determina o “pico” da função de verossimilhança e, consequentemente, é uma medida natural de variabilidade quando a população de origem é multivariada normal.
Resultado 4.11. Se \(\mathbf{X} = [X_1\; X_2\; \cdots\; X_p]^{\prime}\sim \mathcal{N}_p(\boldsymbol{\mu},\mathbf{\Sigma})\), então: (i) os estimadores (variáveis aleatórias) de máxima verossimilhança de \(\boldsymbol{\mu}\) e \(\mathbf{\Sigma}\) são, respectivamente:
\[ \begin{align} \hat{\boldsymbol{\mu}} &= \bar{\mathbf{X}}\\ \hat{\mathbf{\Sigma}} &=\mathbf{S}_n=\dfrac{n-1}{n}\mathbf{S} \end{align} \]
\[ \begin{align} \hat{\boldsymbol{\mu}} &= \bar{\mathbf{x}}\\ \hat{\mathbf{\Sigma}} &=\mathbf{s}_n=\dfrac{n-1}{n}\mathbf{s} \end{align} \]
Os estimadores de máxima verossimilhança possuem uma propriedade de invariância. Seja \(\hat{\boldsymbol{\theta}}\) o estimador de máxima verossimilhança para \(\boldsymbol{\theta}\), e considere a a função \(h(\boldsymbol{\theta})\), que é uma função de \(\boldsymbol{\theta}\). Então, a estimativa de máxima verossimilhança para \(h(\boldsymbol{\theta})\) é dada por \(h(\hat{\boldsymbol{\theta}})\), em que \(\hat{\boldsymbol{\theta}}\) é o estimador de máxima verossimilhança para \(\boldsymbol{\theta}\).
Exemplos com \(\boldsymbol{\theta}=(\boldsymbol{\mu}, \mathbf{\Sigma})\):
O estimador de máxima verossimilhança para \(\boldsymbol{\mu}^{\prime}\mathbf{\Sigma}^{-1}\boldsymbol{\mu}\) é \(\hat{\boldsymbol{\mu}}^{\prime}\hat{\mathbf{\Sigma}}^{-1}\hat{{\boldsymbol{\mu}}}\).
O estimador de máxima verossimilhança para \(\sqrt{\sigma}_{ii}\) é \(\sqrt{\hat\sigma}_{ii}\).
A estatística é suficiente em relação a um modelo estatístico e seu parâmetro desconhecido associado se “nenhuma outra estatística que possa ser calculada a partir da mesma amostra fornece qualquer informação adicional quanto ao valor do parâmetro”, conforme Fisher (1922).
A partir da expressão (4-15), a densidade conjunta depende do conjunto completo de observações \(\mathbf{X}_1, \mathbf{X}_2, \ldots, \mathbf{X}_n\) somente através da média amostral \(\mathbf{\bar{X}}\) e da matriz de soma de quadrados e produtos cruzados \((\mathbf{X}_i - \mathbf{\bar{X}})(\mathbf{X}_i - \mathbf{\bar{X}})^{\prime}\) ponderada por \(n - 1\), também conhecida como \(\mathbf{S}\). Expressamos esse fato dizendo que \(\mathbf{\bar{X}}\) e \(\mathbf{S}\) (ou \((n - 1)\mathbf{S}\)) são estatísticas suficientes:
Se \(\mathbf{X}_1, \mathbf{X}_2, \ldots, \mathbf{X}_n\) é uma amostra aleatória de uma população normal multivariada com média \(\boldsymbol{\mu}\) e covariância \(\mathbf{\Sigma}\), então
\[ \mathbf{\bar{X}} \; \text{e} \; \mathbf{S} \; \text{são estatísticas suficientes} \tag{4-21} \]
A importância das estatísticas suficientes para populações normais é que toda a informação sobre \(\boldsymbol{\mu}\) e \(\mathbf{\Sigma}\) na matriz de dados \(\boldsymbol{\mathcal{X}}\) está contida em \(\mathbf{\bar{X}}\) e \(\mathbf{S}\), independentemente do tamanho da amostra \(n\). Isso geralmente não é verdade para populações não normais. Como muitas técnicas multivariadas começam com médias e covariâncias amostrais, é prudente verificar a adequação da suposição de normalidade multivariada. Se os dados não puderem ser considerados normalmente distribuídos multivariados, as técnicas que dependem apenas de \(\mathbf{\bar{X}}\) e \(\mathbf{S}\) podem estar ignorando outras informações úteis da amostra (assimetria, curtose etc).
A suposição preliminar de que \(\mathbf{X}_1, \mathbf{X}_2, \ldots, \mathbf{X}_n\) constituem uma amostra aleatória de uma população normal com média \(\boldsymbol{\mu}\) e covariância \(\mathbf{\Sigma}\) determina completamente as distribuições amostrais de \(\bar{\mathbf{X}}\) e \(\mathbf{S}\). Apresentamos aqui os resultados sobre as distribuições amostrais de \(\bar{\mathbf{X}}\) e \(\mathbf{S}\) traçando um paralelo com as conclusões familiares univariadas.
No caso univariado (\(p = 1\)), sabemos que
\[ \bar{X}\sim\mathcal{N}\left(\mu,\dfrac{\sigma^2}{n}\right) \]
O resultado para o caso multivariado (\(p \ge 2\)) é análogo, pois
\[ \bar{\mathbf{X}} \sim \mathcal{N}_p\left(\boldsymbol{\mu},\dfrac{\mathbf{\Sigma}}{n}\right) \]
Para a variância amostral, lembre-se de que \((n - 1)S^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2\) é distribuída como \(\sigma^2\) vezes uma variável qui-quadrado com \(n - 1\) graus de liberdade. Por sua vez, essa qui-quadrado é a distribuição de uma soma de quadrados de variáveis aleatórias normais padrão independentes. Ou seja, \((n - 1)S^2\) é distribuída como \(\sigma^2\sum_{i=1}^{n}{Z_i^2} = \sum_{i=1}^{n}{(\sigma Z_i)^2}\). Os termos individuais \(\sigma Z_i\) são distribuídos de forma independente como \(\mathcal{N}(0, \sigma^2)\). É esta última forma que é convenientemente generalizada para a distribuição amostral básica para a matriz de covariância amostral.
Resumindo para \(p=1\):
\[ \begin{align} \bar{X}&\sim\mathcal{N}\left(\mu,\dfrac{\sigma^2}{n}\right) \\ (n - 1)\dfrac{S^2}{\sigma^2} &\sim \chi^2_{n-1} \end{align} \]
Para \(p\ge2\):
\[ \begin{align} \bar{\mathbf{X}}&\sim\mathcal{N}_p\left(\boldsymbol{\mu},\dfrac{\mathbf{\Sigma}}{n}\right) \\ (n - 1)\mathbf{S} &\sim \mathcal{W}_{p}(n-1,\mathbf{\Sigma}) \end{align}\tag{4-23} \]
\(\mathcal{W}_{p}(n-1,\mathbf{\Sigma})\) é a distribuição de Wishart com \(n-1\) graus de liberdade. Se \(p=1\), então:
\[ (n - 1)S^2 \sim \sigma^2\chi^2_{n-1}=\mathcal{W}_{1}(n-1,\sigma^2) \]
\[ \sum_{i=1}^{k}{\mathbf{W}_i}\sim \mathcal{W}_p\left(\mathbf{\Sigma}, \sum_{i=1}^{k}{n_i}\right) \]
2.Se \(\mathbf{W}\sim \mathcal{W}_p\left(\mathbf{\Sigma}, n-1\right)\) e \(\underset{p\times 1}{\mathbf{a}}\) é um vetor de constantes, então:
\[ \dfrac{\mathbf{a}^{\prime}\mathbf{W}\mathbf{a}}{\mathbf{a}^{\prime}\mathbf{\Sigma}\mathbf{a}} \sim \chi^2_{n-1} \]
\[ \mathbf{C}\mathbf{W}\mathbf{C}^T \sim \mathcal{W}_q\left(\mathbf{C}\mathbf{\Sigma}\mathbf{C}^T, m\right) \]
Aqui está o significado dos termos envolvidos:
Esta propriedade é frequentemente usada em estatísticas multivariadas para simplificar expressões complexas envolvendo matrizes com distribuição Wishart e é fundamental no estudo de estruturas de covariância.
, , 1
[,1] [,2]
[1,] 7.162881 2.480879
[2,] 2.480879 1.386713
, , 2
[,1] [,2]
[1,] 0.3579103 0.4934989
[2,] 0.4934989 0.7119320
, , 3
[,1] [,2]
[1,] 2.499046 -3.262622
[2,] -3.262622 4.710663
, , 4
[,1] [,2]
[1,] 0.1598054 -0.6526172
[2,] -0.6526172 2.7345820
Além disso, \(\bar{\mathbf{X}}\) e \(\mathbf{S}\) são independentes, i.e, \(\bar{\mathbf{X}}\perp \mathbf{S}\).
A prova desta independência está na fatoração da função de verossimilhança com variáveis aleatórias em dois termos separam \(\bar{\mathbf{X}}\) e \(\mathbf{S}\) (adaptada de 4-15):
\[ \begin{align} \text{L}(\boldsymbol{\mu}, \mathbf{\Sigma}) &= |2\pi\mathbf{\Sigma}|^{-n/2}\exp\left(-\dfrac{1}{2}\mathrm{tr}\left(\left(\dfrac{\mathbf{\Sigma}}{n}\right)^{-1}\dfrac{n-1}{n}\mathbf{S}\right)\right)\times\\ &\times\exp\left(-\dfrac{1}{2}(\bar{\mathbf{X}}-\boldsymbol{\mu})^{\prime}\left(\dfrac{\mathbf{\Sigma}}{n}\right)^{-1}(\bar{\mathbf{X}}-\boldsymbol{\mu})\right)\\ \text{L}(\boldsymbol{\mu}, \mathbf{\Sigma}) &\propto g(\mathbf{S})\times h(\bar{\mathbf{X}}) \end{align} \]
Porque \(\mathbf{\Sigma}\) é desconhecida, a distribuição de \(\bar{\mathbf{X}}\) não pode ser usada diretamente para fazer inferências sobre \(\boldsymbol{\mu}\). No entanto, \(\mathbf{S}\) fornece informações independentes sobre \(\mathbf{\Sigma}\), e a distribuição de \(\mathbf{S}\) não depende de \(\boldsymbol{\mu}\). Isso nos permite construir uma estatística para fazer inferências sobre \(\boldsymbol{\mu}\), como veremos no Capítulo 5.
Suponha que a quantidade \(X\) seja determinada por um grande número de causas independentes \(V_1, V_2,\ldots,V_n\), onde as variáveis aleatórias que representam as causas têm aproximadamente a mesma variabilidade. Se \(X\) for a soma dessas variáveis, \(X=\sum_{i=1}^{n}{V_i}\), então o teorema do limite central (TLC) se aplica, e concluímos que \(X\) possui uma distribuição que é quase normal. Isso é verdade para praticamente quaisquer distribuições das variáveis \(V_i\), desde que \(n\) seja suficientemente grande (\(n>30\)).
O TLC univariado também nos diz que a distribuição amostral da média amostral, \(\bar{X}\), para um tamanho de amostra grande (\(n>30\)) é quase normal, independentemente da forma da distribuição subjacente da população. Um resultado semelhante vale para muitas outras estatísticas importantes univariadas (e.g., variância amostral).
Acontece que certas estatísticas multivariadas, como \(\bar{\mathbf{X}}\) e \(\mathbf{S}\), possuem propriedades em amostras grandes análogas às suas contrapartes univariadas. À medida que o tamanho da amostra aumenta indefinidamente, certas regularidades governam a variação amostral em \(\bar{\mathbf{X}}\) e \(\mathbf{S}\), independentemente da forma da população original. Portanto, as conclusões apresentadas nesta seção não exigem populações multivariadas normais. Os únicos requisitos são que a população original, seja qual for a sua forma, tenha média \(\boldsymbol{\mu}\) e covariância finita (não-singular) \(\mathbf{\Sigma}\).
Resultado 4.13 (Teorema do limite central multivariado). Sejam \(\mathbf{X}_1, \mathbf{X}_2,\ldots,\mathbf{X}_n\) observações independentes de quaisquer populações (distribuições multivariadas) com média \(\boldsymbol{\mu}\) e covariância não-singular \(\mathbf{\Sigma}\). Então:
\[ \sqrt{n}(\bar{\mathbf{X}} - \boldsymbol{\mu}) \underset{a}{\sim} \mathcal{N}_p(\mathbf{0}, \mathbf{\Sigma}) \]
para tamanho de amostra grande, i.e., \(n-p>30\).
A aproximação fornecida pelo TLC se aplica tanto a populações multivariadas discretas quanto contínuas. Matematicamente, o limite é exato, e a tendência para a normalidade muitas vezes é bastante rápida. Além disso, com base nos resultados na Seção 4.4, sabemos que \(\bar{\mathbf{X}}\) é exatamente distribuído normalmente quando a população subjacente é normal. Portanto, podemos esperar que a aproximação do TLC seja bastante precisa para tamanhos de amostra moderados quando a população é quase normal.
Como vimos, quando \(n\) é grande, \(\mathbf{S}\) se aproxima de \(\mathbf{\Sigma}\) com alta probabilidade. Consequentemente, substituir \(\mathbf{\Sigma}\) por \(\mathbf{S}\) na distribuição normal aproximada para \(\bar{\mathbf{X}}\) terá um efeito desprezível nos cálculos de probabilidade subsequentes.
O Resultado 4.7 pode ser usado para mostrar que \(n(\bar{\mathbf{X}} - \boldsymbol{\mu})^{\prime} \mathbf{\Sigma}^{-1} (\bar{\mathbf{X}} - \boldsymbol{\mu})\) segue a distribuição qui-quadrado com \(p\) graus de liberdade quando \(\bar{\mathbf{X}}\) é distribuído como \(\mathcal{N}_p(\boldsymbol{\mu}, \mathbf{\Sigma}/n)\) ou equivalente, quando \(\sqrt{n} (\bar{\mathbf{X}} - \boldsymbol{\mu})\) segue uma distribuição \(\mathcal{N}_p(\mathbf{0},\mathbf{\Sigma})\). A distribuição qui-quadrado com \(p\) graus de liberdade é uma aproximação da distribuição amostral de \(n(\bar{\mathbf{X}} - \boldsymbol{\mu})^{\prime} \mathbf{\Sigma}^{-1} (\bar{\mathbf{X}} - \boldsymbol{\mu})\) quando \(\bar{\mathbf{X}}\) está aproximadamente distribuído normalmente. Substituir \(\mathbf{\Sigma}^{-1}\) por \(\mathbf{S}^{-1}\) não afeta seriamente essa aproximação para \(n-p>30\).
Resumimos as principais conclusões desta seção da seguinte forma:
Sejam \(\mathbf{X}_1, \mathbf{X}_2,\ldots,\mathbf{X}_n\) observações independentes de uma população com média \(\boldsymbol{\mu}\) e covariância não singular \(\mathbf{\Sigma}\). Então,
\[ \sqrt{n}(\bar{\mathbf{X}} - \boldsymbol{\mu}) \underset{a}{\sim} \mathcal{N}_p(\mathbf{0}, \mathbf{\Sigma}) \]
e
\[ n(\bar{\mathbf{X}} - \boldsymbol{\mu})^{\prime} \mathbf{\Sigma}^{-1} (\bar{\mathbf{X}} - \boldsymbol{\mu})\underset{a}{\sim}\chi_{p}^{2} \]
para \(n - p>30\).
Nas próximas três seções, consideramos maneiras de verificar a
suposição de normalidade e métodos para transformação não-linear das
observações não normais em observações que são aproximadamente
simétricas, e.g., transformação potência de Tukey
(rcompanion::transformTukey
), Box-Cox
(MASS::boxcox
, cellWise::transfo
e
MVN::mvn
) ou Yeo–Johnson (cellWise::transfo
e
e MVN::mvn
).
Se \(n-p>30\), não é necessário aplicar transformação de simetrização em observações não normais, pois TCL é aplicável.
Se \(n-p\le30\), a determinação da transformação de simetrização ótima é prejudicada pelo tamanho da amostra pequeno.
Como já apontamos, a maioria das técnicas estatísticas discutidas nos capítulos subsequentes assume que cada observação vetorial \(\mathbf{X}_i\) provém de uma distribuição normal multivariada. Por outro lado, em situações em que o tamanho da amostra é grande e as técnicas dependem exclusivamente do comportamento de \(\bar{\mathbf{X}}\), ou distâncias envolvendo \(\bar{\mathbf{X}}\) na forma de \(n(\bar{\mathbf{X}} - \boldsymbol{\mu})^{\prime}\mathbf{S}^{-1}(\bar{\mathbf{X}} - \boldsymbol{\mu})\), a suposição de normalidade das observações individuais é menos crucial. No entanto, até certo ponto, a qualidade das inferências feitas por esses métodos depende de quão de perto a verdadeira população de origem se assemelha à forma normal multivariada. É imperativo, então, que existam procedimentos para detectar casos em que os dados apresentam desvios moderados a extremos do que é esperado sob a normalidade multivariada.
Queremos responder a esta pergunta: As observações \(\mathbf{X}_i\) parecem violar a suposição de que elas provêm de uma população normal? Com base nas propriedades das distribuições normais, sabemos que todas as combinações lineares de variáveis normais são normais e que as curvas de densidade normal multivariada são elipsóides. Portanto, abordamos as seguintes questões:
Ficará claro que nossas investigações sobre a normalidade se concentrarão no comportamento das observações em uma ou duas dimensões (por exemplo, distribuições marginais e gráficos de dispersão). Como seria de se esperar, tem sido difícil construir um “bom” teste geral de normalidade conjunta em mais de duas dimensões devido ao grande número de problemas que podem surgir. Até certo ponto, devemos pagar um preço por nos concentrarmos em exames univariados e bivariados da normalidade: nunca podemos ter certeza de que não deixamos passar alguma característica que só é revelada em dimensões superiores. (É possível, por exemplo, construir uma distribuição bivariada não normal com marginais normais. [Ver Exercício 4.8.]) No entanto, muitos tipos de não normalidade são frequentemente refletidos nas distribuições marginais e nos gráficos de dispersão. Além disso, para a maioria dos trabalhos práticos, investigações unidimensionais e bidimensionais são normalmente suficientes. Felizmente, conjuntos de dados patológicos que são normais em representações de dimensão inferior, mas não normais em dimensões superiores, não são frequentemente encontrados na prática.
Histograma não é adequado para avaliar a normalidade de variável intervalar, conforme Silveira e Siqueira (2022).
O departamento de controle de qualidade de um fabricante de fornos de micro-ondas é obrigado pelo governo federal a monitorar a quantidade de radiação emitida quando as portas dos fornos estão fechadas. Observações da radiação emitida através das portas fechadas de \(n = 42\) fornos selecionados aleatoriamente foram feitas. Os dados estão listados na Tabela 4.1.
T4.1 <- read.table("JW6Data/T4-1.DAT",
quote="\"", comment.char="")
colnames(T4.1) <- c("Radiation")
print(psych::describe(T4.1$Radiation))
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 42 0.13 0.1 0.1 0.12 0.07 0.01 0.4 0.39 1.19 0.67 0.02
Hartigans' dip test for unimodality / multimodality
data: T4.1$Radiation
D = 0.059524, p-value = 0.2921
alternative hypothesis: non-unimodal, i.e., at least bimodal
m-out-of-n bootstrap symmetry test by Miao, Gel, and Gastwirth (2006)
data: T4.1$Radiation
Test statistic = 2.8182, p-value = 0.01
alternative hypothesis: the distribution is asymmetric.
sample estimates:
bootstrap optimal m
26
Shapiro-Wilk normality test
data: T4.1$Radiation
W = 0.85793, p-value = 9.902e-05
Distância euclidiana.
Distância de Mahalanabis.
As fórmulas dadas envolvem cálculos relacionados à distância de Mahalanobis ao quadrado \(d_M^2\) para um par de variáveis \(\mathbf{x}=(x_1,x_2)\) em relação ao centróide \(\bar{\mathbf{x}}=(\bar{x}_1,\bar{x}_2)\), em que:
\[ d_M^2 (\mathbf{x}) = [z_1 \quad z_2] \begin{bmatrix} 1 & r \\ r & 1 \end{bmatrix}^{-1} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} \]
Esta fórmula pode ser simplificada como:
\[ d_M^2 (\mathbf{x}) = \frac{z_1^2 + z_2^2 - 2r z_1 z_2}{1 - r^2} \]
Em que: - \(z_1\) e \(z_2\) são os valores normalizados das variáveis \(x_1\) e \(x_2\), respectivamente, utilizando as médias \(\bar{x}_1\) e \(\bar{x}_2\) e os desvios-padrão \(s_1\) e \(s_2\) como normalizadores. Ou seja, \(z_1 = (x_1 - \bar{x}_1) / s_1\) e \(z_2 = (x_2 - \bar{x}_2) / s_2\), em que \(r\) é o coeficiente de correlação de Pearson entre as variáveis \(x_1\) e \(x_2\), e \(d_M^2\) é a medida de distância de Mahalanobis ao quadrado.
Essas fórmulas são usadas para calcular a distância de Mahalanobis entre duas variáveis considerando seus desvios-padrão, médias e correlação amostrais. Essa medida é útil para avaliar a distância das observações em relação ao centróide em um espaço multivariado, levando em consideração a relação entre as variáveis.
O conjunto de dados “iris” é um dos conjuntos de dados incorporados ao pacote “datasets” no ambiente de programação estatística R. Este conjunto de dados é frequentemente utilizado como exemplo em análises estatísticas e aprendizado de máquina. Ele contém informações sobre várias medidas de flores de três diferentes espécies de íris: setosa, versicolor e virginica. Cada espécie é representada por 50 observações.
As colunas deste conjunto de dados incluem:
Este conjunto de dados é frequentemente utilizado para ilustrar técnicas de visualização de dados, classificação de padrões e outros métodos estatísticos. É uma escolha popular para iniciantes em análise de dados e aprendizado de máquina devido à sua simplicidade e natureza bem estruturada.
Dados <- as.matrix(subset(x=datasets::iris,
subset=Species=="setosa",
select=-Species))
DM2 <- mahalanobis(Dados,
colMeans(Dados),
cov(Dados))
print(psych::describe(DM2))
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 50 3.92 3.21 2.97 3.45 2.51 0.34 12.33 11.98 1.12 0.31 0.45
# https://cran.r-project.org/web/packages/MVN/vignettes/MVN.pdf
data(iris)
Dados <- iris
head(Dados)
# A tibble: 6 × 5
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
<dbl> <dbl> <dbl> <dbl> <fct>
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
# A tibble: 6 × 5
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
<dbl> <dbl> <dbl> <dbl> <fct>
1 6.7 3.3 5.7 2.5 virginica
2 6.7 3 5.2 2.3 virginica
3 6.3 2.5 5 1.9 virginica
4 6.5 3 5.2 2 virginica
5 6.2 3.4 5.4 2.3 virginica
6 5.9 3 5.1 1.8 virginica
# https://cran.r-project.org/web/packages/MVN/vignettes/MVN.pdf
# setosa <- subset(x=iris,
# subset=Species=="setosa",
# select=-Species)
# result <- MVN::mvn(data = setosa,
# mvn_test = "hz")
# result$multivariate_normality
# result <- try(MVN::mvn(data=iris,
# subset="Species",
# mvn_test="hz",
# univariate_test="SW"))
# result$multivariate_normality
# result$univariate_normality
# setosa <- subset(x=iris,
# subset=Species=="setosa",
# select=c(Sepal.Length:Sepal.Width))
# # perspective plot
# result <- MVN::mvn(setosa,
# mvnTest = "hz",
# multivariatePlot = "persp")
# # contour plot
# result <- MVN::mvn(setosa,
# mvnTest = "hz",
# multivariatePlot = "contour")
# bagplot
# bgp <- DescTools::PlotBag(setosa[,1],setosa[,2],
# main=paste("Bagplot: setosa"),
# xlab="Sepal.Length",
# ylab="Sepal.Width",
# na.rm=TRUE,
# show.bagpoints=FALSE,
# show.looppoints=FALSE,
# show.whiskers=FALSE,
# show.baghull=TRUE,
# col.loophull="white",
# col.looppoints="black",
# col.baghull="white",
# col.bagpoints="black",
# asp=1,
# add=FALSE,
# cex=1)
# print(outliers <- as.data.frame(bgp$pxy.outlier))
# for (o in 1:nrow(outliers))
# {
# r <- which(setosa[,1]==outliers$x[o] &
# setosa[,2]==outliers$y[o])
# text(outliers$x[o],outliers$y[o], r, pos=1, cex=0.7)
# }
# # Mahalanobis distance
# result <- MVN::mvn(data = setosa,
# mvn_test = "hz",
# multivariate_outlier_method = "quan")
# # Adjusted Mahalanobis distance
# result <- MVN::mvn(data = setosa,
# mvn_test = "hz",
# multivariate_outlier_method = "adj")
result <- MVN::mvn(data = iris[-4],
subset = "Species",
power_family="none")
summary(result, select = "mvn")
── Multivariate Normality Test Results ─────────────────────────────────────────
Group Test Statistic p.value MVN
1 setosa Henze-Zirkler 0.524 0.831 ✓ Normal
2 versicolor Henze-Zirkler 0.714 0.326 ✓ Normal
3 virginica Henze-Zirkler 0.726 0.299 ✓ Normal
result <- MVN::mvn(data = iris[-4],
subset = "Species",
power_family="bcPower")
summary(result, select = "mvn")
── Multivariate Normality Test Results ─────────────────────────────────────────
Group Test Statistic p.value MVN
1 setosa Henze-Zirkler 0.559 0.747 ✓ Normal
2 versicolor Henze-Zirkler 0.550 0.771 ✓ Normal
3 virginica Henze-Zirkler 0.703 0.351 ✓ Normal
result <- MVN::mvn(data = iris[-4],
subset = "Species",
power_family="yjPower")
summary(result, select = "mvn")
── Multivariate Normality Test Results ─────────────────────────────────────────
Group Test Statistic p.value MVN
1 setosa Henze-Zirkler 0.577 0.700 ✓ Normal
2 versicolor Henze-Zirkler 0.526 0.827 ✓ Normal
3 virginica Henze-Zirkler 0.709 0.337 ✓ Normal
4.2, 4.3, 4.8, 4.9, 4.17, 4.18, 4.19, 4.21, 4.22, 4.31, 4.39
Soch, J et al. (2020) The Book of Statistical Proofs: https://statproofbook.github.io/
Batschelet, E (1978) Introdução à matemática para biocientistas. Tradução da 2ª ed. São Paulo: EDUSP e Rio de Janeiro: Interciência.
Batschelet, E (1979) Introduction to mathematics for life scientists. 3rd ed. NY: Springer.
Bickel, PJ & Doksum, KA (2001) Mathematical Statistics: Basic Ideas and Selected Topics, Volume I. USA: CRC.
Bilodeau, M & Brenner, D (1999) Theory of Multivariate Statistics. USA: Springer. T^2 de Hotelling
Denis, DJ (2021) Applied Univariate, Bivariate, and Multivariate Statistics Using Python: A Beginners Guide to Advanced Data Analysis. NJ: Wiley.
Genz, A (1992) Numerical Computation of Multivariate Normal Probabilities. Journal of Computational and Graphical Statistics 1(2): 141–149. https://doi.org/10.2307/1390838
Hardle, W & Hlavka, Z (2007) Multivariate Statistics - Exercises and Solutions. USA: Springer.
KASSAMBARA, A (2017) Practical guide to cluster analysis in R: Unsupervised machine learning. STHDA: http://www.sthda.com.
Kollo. T & von Rosen, D (2005) Advanced Multivariate Statistics with Matrices. USA: Springer.
Kutner, MH et al. (2005) Applied linear statistical model. 5th ed. NY: McGraw-Hill/Irwin.
MAIR, P (2018) Modern psychometrics with R. USA: Springer.
Manly, BFJ & Alberto, JAN (2017) Multivariate Statistical
Methods: A Primer using R. 4th ed. USA: CRC.
Matloff, N (2020) Probability and Statistics for Data Science: Math + R Data_. USA: CRC.
Mirman, D (2014) Growth Curve Analysis and Visualization Using R. USA: CRC.
Pessoa, DGC (2008) Exemplos do livro de Johnson e Wichern.
Rawlings, JO et al. (1998) Applied Regression Analysis: A Research Tool. USA: Springer. T^2 de Hotelling
Schumacker, RE (2016) Using R With Multivariate Statistics. USA: SAGE.
Siqueira, JO (2012) Fundamentos de métodos quantitativos: aplicados em Administração, Economia, Contabilidade e Atuária usando WolframAlpha e SCILAB. São Paulo: Saraiva. Soluções dos exercícios em https://www.researchgate.net/publication/326533772_Fundamentos_de_metodos_quantitativos_-_Solucoes.
Tay, A (2018) OLS using Matrix Algebra. http://www.mysmu.edu/faculty/anthonytay/MFE/OLS_using_Matrix_Algebra.pdf
Teichroew, D (1964) An introduction to management science: deterministic models. NJ: Wiley.
Timm, NH (2002) Applied Multivariate Analysis. USA: Springer. T^2 de Hotelling
Zelterman, D (2015) Applied Multivariate Statistics with R. USA: Springer.
Nguyen, M (2020) A Guide on Data Analysis. Bookdown. https://bookdown.org/mike/data_analysis/
Chaves Neto, A. CE-704 ANÁLISE MULTIVARIADA APLICADA À PESQUISA. https://docs.ufpr.br/~soniaisoldi/ce090/CE076AM_2010.pdf
ME 731 - Métodos em Análise Multivariada em R. https://www.ime.unicamp.br/~cnaber/Material_AM_2S_2020.htm
Nogueira, FE (2007) Modelos de regressão multivariada. Dissertação (Mestrado). https://teses.usp.br/teses/disponiveis/45/45133/tde-25062007-163150/publico/dissertacao_4.pdf
Análise Multivariada. https://www.professores.uff.br/samuelcampos/analise-multivariada/
Severn, K (2023) Multivariate Statistics. https://rich-d-wilkinson.github.io/MATH3030/
Powell, J (2019) Multivariate statistics in R. https://www.hiercourse.com/multivariate
Plotting PCA/clustering results using ggplot2 and ggfortify: https://rpubs.com/sinhrks/plot_pca
ggfortify : Extension to ggplot2 to handle some popular packages - R software and data visualization: http://www.sthda.com/english/wiki/ggfortify-extension-to-ggplot2-to-handle-some-popular-packages-r-software-and-data-visualization
A Little Book of Python for Multivariate Analysis: https://python-for-multivariate-analysis.readthedocs.io/index.html