7.1 Correlação

Correlação é uma medida do grau de associação linear entre duas variáveis aleatórias. No senso comum possui uma ampla gama de significados, mas o termo usualmente se refere ao coeficiente de correlação de Pearson, detalhado na Seção 7.1.2.

Exemplo 7.1 A medida chamada frequência fundamental ou altura (pitch) pode indicar a frequência em Hz de vibração das cordas vocais. É largamente reportada na literatura especializada, da qual consideram-se os algoritmos de (Schäfer-Vincent 1983) e (Jackson 1995). Utilizando o pacote voice (Zabala 2022) é possível extrair estes dois traços fonéticos, comparados a seguir.

wavFiles <- dir(system.file('extdata', package = 'wrassp'),
                pattern <- glob2rx('*.wav'), full.names = TRUE)
M <- voice::extract_features(wavFiles, features = c('f0','f0_mhs'))
M <- M[M$f0_mhs < 125 | is.na(M$f0_mhs),] # removendo pontos anômalos
plot(M$f0, M$f0_mhs, las = 1)
txt <- paste('r =', round(cor(M$f0, M$f0_mhs, use = 'pair'), 3))
legend(55, 110, txt, bty = 'n')

Exercício 7.1

1. Analise a documentação do objeto anscombe.
   
2. Acesse os links datasauRus e https://marcusnunes.me/posts/correlacoes-e-graficos-de-dispersao/. \(\\\)

Exercício 7.2 Acesse o link sobre Bivariate Correlation Simulation and Bivariate Scatter Plotting e faça simulações variando os parâmetros, observando os gráficos resultantes.

7.1.1 \(\rho\) \(\cdot\) Correlação universal

(Samuel Kotz et al. 2005, 1375) definem a correlação universal de duas v.a. \(X\) e \(Y\) por \[\begin{equation} \rho = Cor(X,Y) = \dfrac{Cov(X,Y)}{D(X) D(Y)}, \tag{7.1} \end{equation}\] onde \[\begin{equation} Cov(X,Y) = E \left[ (X-E(X)) (Y-E(Y)) \right] = \frac{\sum_{i=1}^N (x_{i} - \mu_{X})(y_{i} - \mu_{Y})}{N} \tag{7.2} \end{equation}\]
é a covariância entre \(X\) e \(Y\), \(D(X)\) e \(D(Y)\) são respectivamente os desvios padrão universais de \(X\) e \(Y\) conforme Eq. (2.34) e \[\begin{equation} -1 \le \rho \le +1 \end{equation}\] \[ \Updownarrow \] \[\begin{equation} 0 \le |\rho| \le +1. \end{equation}\]

Se \(|\rho| = +1\), então existe uma relação linear da forma \(Y = \beta_0 + \beta_1 X\). Se \(\rho = +1\), \(\beta_1 >0\); se \(\rho = -1\), \(\beta_1 <0\). Se \(X\) é independente de \(Y\), então \(\rho = 0\), mas o contrário não é necessariamente verdadeiro.

7.1.2 \(r\) \(\cdot\) Coeficiente de Correlação de Pearson

O coeficiente de correlação (amostral) de Pearson foi desenvolvido por Karl Pearson e William Fleetwood Sheppard “antes do final de 1897”, conforme relatado por (Pearson 1920, 44). Denotado por \(r\), pode ser obtido por qualquer das equações a seguir.

\[\begin{equation} r = \frac{1}{n-1} \sum_{i=1}^{n} \left( \dfrac{x_{i} - \bar{x}}{s_{x}} \right) \left( \dfrac{y_{i} - \bar{y}}{s_{y}} \right) \tag{7.3} \end{equation}\]

\[\begin{equation} r = \frac{\sum (x - \bar{x})(y - \bar{y})}{\sqrt{\sum (x - \bar{x})^2 \sum (y - \bar{y})^2}} \tag{7.4} \end{equation}\]

\[\begin{equation} r = \frac{n \sum{x y} - \sum{x} \sum{y}}{\sqrt{[n \sum{x^2} - (\sum{x})^2] [n \sum{y^2} - (\sum{y})^2]}} \tag{7.5} \end{equation}\]

onde \[ \bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_{i}, \;\; s_{x}^2 = \frac{1}{n-1} \sum_{i=1}^{n} (x_{i} - \bar{x}) ^2 \] \[ \bar{y} = \frac{1}{n} \sum_{i=1}^{n} y_{i}, \;\; s_{y}^2 = \frac{1}{n-1} \sum_{i=1}^{n} (y_{i} - \bar{y}) ^2 \] \[ -1 \le r \le +1 \]

Note pela Equação (7.3) que \(r\) é uma média dos produtos dos pares ordenados \((x_{i}, y_{i})\) padronizados, com \(i \in \lbrace 1, 2, \ldots, n \rbrace\). Se os pares de produto positivo predominarem, \(r\) será positivo; se houver predominância dos pares de produto negativo, \(r\) será negativo. Esta estrutura é chamada de momento-produto. A Equação (7.4) remete à definição (7.1), enquanto a Equação (7.5) é útil para a realização dos cálculos.

Exemplo 7.2 Considere a idéia de estimar o número de garrafas de bebida a serem geladas dependendo da temperatura máxima do dia. Seja \(X\): ‘temperatura máxima do dia em \(^{\circ}{\rm C}\)’ e \(Y\): ‘número de garrafas de bebida consumidas’, observadas conforme tabela a seguir.

\(i\)	\(x_{i}\)	\(y_{i}\)	\(i\)	\(x_{i}\)	\(y_{i}\)	\(i\)	\(x_{i}\)	\(y_{i}\)
1	29.5	146	11	28.5	183	21	40.9	233
2	31.3	170	12	28.0	158	22	28.6	169
3	34.7	167	13	36.7	181	23	36.1	192
4	40.4	244	14	31.5	123	24	27.1	106
5	28.4	159	15	38.1	223	25	29.5	170
6	40.3	195	16	33.5	176	26	31.6	167
7	41.1	225	17	37.2	196	27	25.2	133
8	36.2	206	18	41.9	238	28	31.5	138
9	35.7	200	19	31.5	184	29	39.8	199
10	26.1	134	20	38.2	213	30	30.8	172

Um gráfico de dispersão pode ajudar a explorar o comportamento da temperatura e garrafas consumidas.

dr <- read.table('https://filipezabala.com/data/drinks.txt', header = T)
dim(dr)  # dimensão do banco de dados

## [1] 30  2

head(dr) # início do banco de dados

##   temp gar
## 1 29.5 146
## 2 31.3 170
## 3 34.7 167
## 4 40.4 244
## 5 28.4 159
## 6 40.3 195

plot(dr$temp, dr$gar)

O grau de alinhamento das variáveis pode ser estimado pelo coeficiente de correlação de Pearson, bastando calcular \[ \sum{x} = 1010, \; \sum{x}^2 = 34730, \] \[ \sum{y} = 5400, \; \sum{y}^2 = 1006954, \] \[ \sum{x y} = 186117, \; n=30 \]

e substituir na Equação (7.5), resultando em

\[\begin{equation} r = \frac{30 \times 186117 - 1010 \times 5400}{\sqrt{[30 \times 34730 - 1010^2] [30 \times 1006954 - 5400^2]}} \nonumber \\ r = \frac{130056}{\sqrt{21988.49 \times 1048620}} \nonumber \\ r \approx 0.8565. \nonumber \end{equation}\]

cor(dr$temp, dr$gar)

## [1] 0.8564926

É possível obter o valor de \(r\) a partir da função stats::lm. Ela calcula o coeficiente de determinação \(r^2\) conforme Seção 7.2.3.1, que no modelo RLS (Seção 7.2) é o quadrado de \(r\). Logo \(r=\sqrt{r^2}\).

f <- lm(dr$gar ~ dr$temp)                   # Modelo RLS
r2 <- summary(f)$r.squared                  # Extraindo r2
sqrt(r2)                                    # r, a correlação de Pearson

## [1] 0.8564926

all.equal(cor(dr$temp, dr$gar), sqrt(r2))   # Comparando r via stats::cor e stats::lm

## [1] TRUE

Exercício 7.3 Calcule o valor de \(r\) com os dados do Exemplo 7.2 a partir das Equações (7.3) e (7.4). Compare com os valores obtidos via Eq. (7.5) e pela função stats::cor.

Teste para \(\rho\)

O teste básico é comparar \(\rho\) com zero, que indica ausência completa de alinhamento entre as variáveis. Assim, testa-se \(H_0: \rho = 0\) (não há correlação) vs \(H_1: \rho \ne 0\) (há correlação), denotado por

\[ \left\{ \begin{array}{l} H_0: \rho = 0 \\ H_1: \rho \ne 0\\ \end{array} \right. . \]

Sob \(H_0\) considerando o modelo RLS (Seção 7.2), \[\begin{equation} T = \frac{r}{\sqrt{\frac{1-r^2}{n-2}}} \sim t_{n-2}. \tag{7.6} \end{equation}\]

Exercício 7.4 A partir da Eq. (7.6), mostre que \(T = r \sqrt{\frac{n-2}{1-r^2}}\).

Exemplo 7.3 (Verificando o alinhamento no modelo RLS) Considere novamente as informações apresentadas no Exemplo 7.2. Considerando a estimativa do modelo \(y = -19.11 + 5.91 x\), pode-se testar \(H_0: \rho = 0\) vs \(H_1: \rho \ne 0\). Neste caso \(gl=30-2=28\), \(T \sim t_{28}\) e sob \(H_0\) \[ T = 0.8565 \sqrt{\frac{30-2}{1-0.8565^2}} \approx 8.78. \]

n <- 30
gl <- n-2
(r <- cor(dr$temp, dr$gar))

## [1] 0.8564926

(Tt <- r*sqrt(gl/(1-r^2)))

## [1] 8.780493

(p_value <- 2*pt(-abs(Tt), gl))

## [1] 1.565265e-09

cor.test(dr$temp, dr$gar)         # Função que realiza o teste de hipótese

## 
##  Pearson's product-moment correlation
## 
## data:  dr$temp and dr$gar
## t = 8.7805, df = 28, p-value = 1.565e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7176748 0.9298423
## sample estimates:
##       cor 
## 0.8564926

all.equal(r, sqrt(r2))            # Comparado r (calculado) e r2 (via stats::lm)

## [1] TRUE

IC para \(\rho\)

Passo 1: Obtenha a transformação z de Fisher (Ronald A. Fisher 1921, 210).

\[\begin{equation} z_r = \frac{1}{2} \ln \left( \frac{1+r}{1-r} \right) = \tanh^{-1}(r) \tag{7.7} \end{equation}\]

Passo 2: Obtenha os limitantes inferior e superior transformados com a confiança desejada \(1-\alpha\).

\[\begin{equation} L = z_r - \frac{z_{1-\alpha/2}}{\sqrt{n-3}} \tag{7.8} \end{equation}\]

\[\begin{equation} U = z_r + \frac{z_{1-\alpha/2}}{\sqrt{n-3}} \tag{7.9} \end{equation}\]

Passo 3: Obtenha o intervalo de confiança \(1-\alpha\) para \(\rho\).

\[\begin{equation} IC[\rho, 1-\alpha] = \left[ \frac{e^{2L}-1}{e^{2L}+1}, \frac{e^{2U}-1}{e^{2U}+1} \right] = [ \tanh L, \tanh U ] \tag{7.10} \end{equation}\]

Exemplo 7.4 Considerando os dados do Exemplo 7.3 pode-se obter o IC 95% para a correlação.

# Passo 1
(zr <- atanh(r))

## [1] 1.280029

# Passo 2
(ep <- 1.96/sqrt(n-3))

## [1] 0.3772022

(L <- zr-ep); (U <- zr+ep)

## [1] 0.9028267

## [1] 1.657231

# Passo 3
tanh(L); tanh(U)

## [1] 0.7176715

## [1] 0.9298433

Exercício 7.5

1. Verifique a identidade das Equações (7.7) e (7.10). Aplique as equações alternativas aos dados do Exemplo 7.4 e compare com a saída da função stats::cor.test.
   
2. Obtenha \(IC[ \rho, 80\% ]\).
   
3. Escreva uma função que retorne o IC para \(\rho\) em função de dois vetores e de uma confiança arbitrária.

7.1.3 \(\rho_{RPO}\) e \(r_{RPO}\) \(\cdot\) Correlação na RPO

Existe um caso especial de cálculo de correlação chamado Regressão Pela Origem (RPO), baseado no modelo conforme Eq. (7.17). Neste caso calcula-se \(r_{RPO}\) através da expressão

\[\begin{equation} r_{RPO} = \sqrt{\dfrac{\sum \hat{y}_{i}^2}{\sum y_{i}^2}}, \tag{7.11} \end{equation}\]

tal que \(\hat{y}_{i} = \hat{\beta}_1^{RPO}x_i\).

Exemplo 7.5 (Correlação na RPO) Considere as informações do Exemplo 7.2. Pode-se calcular \[ r_{RPO} = \sqrt{\dfrac{997410}{1006954}} = \sqrt{0.9905} \approx 0.9952. \]

f <- lm(gar ~ temp-1, data = dr)
yhat <- sum(f$fitted.values^2)
y2 <- sum(dr$gar^2)
yhat/y2

## [1] 0.9905222

É possível obter o valor de \(r_{RPO}\) a partir da função stats::lm. Ela calcula o coeficiente de determinação \(r_{RPO}^2\), que no modelo RPO é o quadrado de \(r_{RPO}\). Logo \(r_{RPO}=\sqrt{r_{RPO}^2}\).

f_rpo <- lm(dr$gar ~ dr$temp - 1)         # Modelo RPO
r2_rpo <- summary(f_rpo)$r.squared        # Função para obter r2_RPO
sqrt(r2_rpo)                              # r_RPO

## [1] 0.9952498

all.equal(sqrt(yhat/y2), sqrt(r2_rpo))    # Comparado r_rpo (calculado) e r2_rpo (via stats::lm)

## [1] TRUE

Exercício 7.6 Considere novamente as informações apresentadas no Exemplo 7.2.
a. Do Exemplo 7.12 sabe-se que \(\hat{y}_{i} = 5.36 x_i\), tal que \(x_i \in \{29.5, 31.3, \ldots, 30.8 \}\). Mostre que \(\sum \hat{y}_{i}^2 = 997410\).
b. Inspecione o objeto f$fitted.values do Exemplo 7.5.

Teste para \(\rho_{RPO}\)

Se considerarmos o modelo RPO conforme Eq. (7.17), pode-se testar as hipóteses

\[ \left\{ \begin{array}{l} H_0: \rho_{RPO} = 0 \\ H_1: \rho_{RPO} \ne 0\\ \end{array} \right. . \]

Sob \(H_0\), a estatística do teste é \[\begin{equation} T_{RPO} = \frac{r_{RPO}}{\sqrt{\frac{1-r_{RPO}^2}{n-1}}} \sim t_{n-1}. \tag{7.12} \end{equation}\]

Exemplo 7.6 (Verificando o alinhamento no modelo RPO) Considere novamente as informações apresentadas no Exemplo 7.2. Se a estimativa do modelo RPO é \(y = 5.36 x\), pode-se testar

\[ \left\{ \begin{array}{l} H_0: \rho^{RPO} = 0 \\ H_1: \rho^{RPO} \ne 0\\ \end{array} \right. \]

Neste caso \(gl=30-1=29\), \(T \sim t_{29}\) e sob \(H_0\) \[ T = 0.9952 \sqrt{\frac{30-1}{1-0.9952^2}} \approx 55.05. \]

n <- 30
gl_rpo <- n-1
(r_rpo <- sqrt(yhat/y2))

## [1] 0.9952498

(Tt_rpo <- r_rpo*sqrt(gl_rpo/(1-r_rpo^2)))

## [1] 55.05267

(p_value <- 2*pt(-abs(Tt_rpo), gl_rpo))

## [1] 6.778938e-31

IC para \(\rho_{RPO}\)

7.1.4 Coeficiente de correlação (ordinal) de Spearman

7.1.5 Tau de Kendall

Referências

Fisher, Ronald A. 1921. “014: On the "Probable Error" of a Coefficient of Correlation Deduced from a Small Sample.” Metron. https://digital.library.adelaide.edu.au/dspace/bitstream/2440/15169/1/14.pdf.

Jackson, Jeffrey C. 1995. “The Harmonic Sieve: A Novel Application of Fourier Analysis to Machine Learning Theory and Practice.” Carnegie-Mellon University Pittsburgh PA Schoo of Computer Science. https://apps.dtic.mil/sti/pdfs/ADA303368.pdf.

Kotz, Samuel, Narayanaswamy Balakrishnan, Campbell B Read, and Brani Vidakovic. 2005. Encyclopedia of Statistical Sciences. John Wiley & Sons.

———. 1920. “Notes on the History of Correlation.” Biometrika 13 (1): 25–45. https://www.jstor.org/stable/2331722.

Schäfer-Vincent, Kurt. 1983. “Pitch Period Detection and Chaining: Method and Evaluation.” Phonetica 40 (3): 177–202. https://doi.org/10.1159/000261691.

———. 2022. voice: Tools for Voice Analysis, Speaker Recognition and Mood Inference. https://CRAN.R-project.org/package=voice.