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

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(bootES, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(ppcor, warn.conflicts=FALSE))
suppressMessages(library(ppcor, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
source("eiras.bartitle.R")
source("eiras.col2rgbstring.R")
source("eiras.cor.test.boot.R")
source("eiras.correg.R")
source("eiras.createobj.htest.R")
source("eiras.ellipseaxis.R")
source("eiras.findcommonchars.R")
source("eiras.friendlycolor.R")
source("eiras.jitter.R")
source("eiras.text.leading.R")
alfa <- 0.05

Material

  • HTML de R Markdown em RPubs

Objetivos

  • Avaliar a associação entre duas variáveis intervalares.
  • Definir e construir diagrama de dispersão.
  • Distinguir os conceitos de correlação de Pearson e de regressão linear simples.
  • Calcular e interpretar o coeficiente de correlação de Pearson.
  • Distinguir calcular e testar os coeficientes de correlação de Pearson e Spearman.
  • Testar uma correlação, uma correlação parcial, duas independentes e duas dependentes.

Correlação e Regressão: semelhanças e diferenças

A análise de correlação estima e testa a associação entre duas variáveis intervalares.

A análise de regressão linear simples estima e testa a relação entre as variáveis intervalares de desfecho e independente.

Bhagat, 2018

A diferença entre as análises de correlação e de regressão é uma das perguntas mais frequentes em entrevistas.

Muitas pessoas confundem os dois conceitos.

Comparação por: Correlação Regressão
  • Significado
Medida estatística que determina a associação entre duas variáveis Descreve como uma variável independente é numericamente relacionada com a variável dependente
  • Uso
Representação linear da relação entre duas variáveis Ajuste da melhor função para estimar uma variável dependente com base na independente
  • Variável independente (explicativa, previsora) e dependente (desfecho)
Adirecional. Não importa definir qual variável é independente ou dependente Direcional, da variável independente sem erro de mensuração para a dependente com erro de mensuração
  • Indicação
Adimensional. Indica em qual medida as duas variáveis variam conjuntamente Dimensional. Indica o impacto da mudança de uma unidade de medida da variável conhecida (independente) na variável estimada (dependente)
  • Objetivo
Encontrar um valor numérico para expressar a relação entre duas variáveis Estimar o valor de uma variável aleatória com base nos valores de uma variável fixada

Mohammed, 2018

Uma correlação é uma medida da relação entre duas variáveis intervalares. Cada variável representa um fenômeno específico; quando um dos fenômenos muda em uma direção específica, a segunda variável muda na direção do primeiro ou na direção oposta do primeiro.

A mudança dos dois fenômenos na mesma direção no sentido do aumento no primeiro deslocamento por um aumento no segundo ou vice-versa na redução do primeiro acompanhado por uma diminuição no segundo relacionamento é positiva ou crescente (associação positiva). Ao contrário, quando o aumento no primeiro deslocamento é associado com uma redução no segundo ou vice-versa, uma redução no primeiro fenômeno seja acompanhada por um aumento no segundo, dizemos que a ligação é inversa ou decrescente (associação negativa).

A regressão é um método pelo qual o valor de uma variável pode ser estimado pelo valor da outra variável pela equação de regressão. Tem os seguintes tipos:

  • Regressão linear simples: a palavra “simples” significa que a variável dependente \(Y\) depende de uma variável independente \(X\) e a palavra “linear” significa que a relação entre as variáveis \(Y\) e \(X\) é linear.
  • Regressão linear múltipla: se a variável \(Y\) depende de duas ou mais variáveis independentes.
  • Regressão não linear: se a relação entre a variável \(Y\) e a(s) variável(is) independente(s) é não linear como acontece em relações exponenciais.

Exemplo

Apresentamos dados de um estudo realizado no Centro de Saúde Geraldo de Paula Souza, da Faculdade de Saúde Pública da USP, no qual gestantes foram avaliadas quanto à deficiência de ácido fólico e vitamina B12 no início do pré-natal e no 7º mês de gestação. Um grupo controle de mulheres saudáveis foi constituído para servir de referência para as medidas das vitaminas acima mencionadas. Durante o parto, foram obtidas amostras de sangue do recém-nascido para verificar a concentração de anticorpos contra a rubéola.

A planilha Gestantes.xlsx contém os seguintes campos (colunas):

Variável Descrição
Grupo Gestantes ou Controle
NOME Iniciais do nome das pacientes
IDADE Idade das pacientes em anos
COR Etnia das pacientes: br = branca; pd = mulata; am = amarela
HB Hemoglobina (g/dL) no início do pré-natal
HT Hematócrito (%) no início do pré-natal
HEM Contagem de Hemácias (milhões/mm³) no início do pré-natal
LEUC Contagem de Leucócitos (milhares/mm³) no início do pré-natal
RET Reticulócitos (%) no início do pré-natal
PLAQUET Plaquetas (milhares/mm³) no início do pré-natal
FOLICO Concentração de Ácido Fólico (ng/mL) no início do pré-natal (valores abaixo de 3 ng/mL indicam deficiência de Ácido Fólico)
B12 Concentração de Vitamina B12 (pg/mL) no início do pré-natal (valores abaixo de 400 pg/mL indicam deficiência de Vitamina B12)
FOL_7m Concentração de Ácido Fólico (ng/mL) no sétimo mês de gestação
B12_7m Concentração de Vitamina B12 (pg/mL) no sétimo mês de gestação
Hb_7m Hemoglobina no sétimo mês de gestação
Ht_7m Hematócrito no sétimo mês de gestação
Hm_7m Contagem de Hemácias no sétimo mês de gestação
Leu_7m Contagem de Leucócitos no sétimo mês de gestação
Ret_7m Reticulócitos no sétimo mês de gestação
Plq_7m Contagem de Plaquetas no sétimo mês de gestação
AC_Rub_MAE Anticorpos IgG contra o vírus da Rubéola da gestante (UI/mL)
AC_Rub_RN Anticorpos IgG contra o vírus da Rubéola do recém-nascido (UI/mL)
Primigesta Primeira gestação? sim ou não
TABAGISMO Gestante fuma? sim ou não

Shinohara, 1989

Há diversas variáveis quantitativas nestes dados, obtidos de 40 gestantes e 25 controles. Usaremos apenas algumas das variáveis:

  • HT … Hematócrito (%) no início do pré-natal
  • HB … Hemoglobina (g/dl) no início do pré-natal
  • HEM … Contagem de Hemácias (milhões/mm³) no início do pré-natal
  • LEUC … Contagem de Leucócitos (milhares/mm³) no início do pré-natal

para responder à pergunta:

Podemos estimar a concentração de hemoglobina e as contagens de eritrócitos e de leucócitos no sangue pela medida do hematócrito?

hematócrito

A coleta é mais simples; o processamento requer uma centrífuga e uma régua:

https://www.youtube.com/watch?v=gyJLaXO4RcU

hemograma

A coleta do sangue é mais trabalhosa:

https://www.youtube.com/watch?v=T5oqEUOVYhk

Mais trabalhosa, esta medida requer equipamento e um computador:

https://www.youtube.com/watch?v=7XRIXm4Nl4s

Existem equipamentos mais sofisticados para automatizar o processo:

https://www.youtube.com/watch?v=2ORzUhhBe9Q

Panorama

Este diagrama pode ser útil para definir os testes adequados para cada situação. Para correlação e regressão linear simples, seguimos o caminho no qual temos duas variáveis para as quais precisamos verificar associação.

Modificado de Dancey & Reidy (2019)

Arquivo de dados

Dados <- readxl::read_excel("Gestantes.xlsx")
Dados <- subset(Dados,
                Grupo=="Gestante")
saveRDS(Dados, "Gestante.rds")

Avaliação gráfica

Usaremos o arquivo de dados Gestante.rds, para realizar análise descritiva, textual e gráfica, executando Gestantes_descritiva.R.

Dados <- readRDS("Gestante.rds")
Vars <- c("HT", "HB", "HEM")
print(summary(Dados[Vars]))
print(cor(Dados[Vars]), digits=2)
sunflowerplot(HB~HT, 
              xlab="Hematócrito (%)", ylab="Hemoglobina (mg/dl)",
              data=Dados)
lines(lowess(Dados$HT~Dados$HB), lty=2)
sunflowerplot(HEM~HT, 
              xlab="Hematócrito (%)", ylab="Hemácias (milhão/mm³)",
              data=Dados)
lines(lowess(Dados$HT~Dados$HEM), lty=2)
sunflowerplot(LEUC~HT, 
              xlab="Hematócrito (%)", ylab="Leucócitos (milhar/mm³)",
              data=Dados)
lines(lowess(Dados$HT~Dados$LEUC), lty=2)
       HT             HB             HEM       
 Min.   :34.0   Min.   :11.20   Min.   :3.530  
 1st Qu.:36.0   1st Qu.:12.40   1st Qu.:4.112  
 Median :37.0   Median :12.75   Median :4.280  
 Mean   :37.4   Mean   :13.04   Mean   :4.318  
 3rd Qu.:39.0   3rd Qu.:13.60   3rd Qu.:4.490  
 Max.   :45.0   Max.   :16.40   Max.   :5.380  
      HT   HB  HEM
HT  1.00 0.86 0.69
HB  0.86 1.00 0.59
HEM 0.69 0.59 1.00

Observando os gráficos, cada ponto é um par de medidas de um indivíduo. Para referência, adicionamos uma linha pontilhada obtida com a função lowess, capaz de encontrar uma curva que passa ponderadamente entre os pontos de um gráfico para nos ajudar a julgar o tipo de relação entre as variáveis.

O hematócrito parece ter relação positiva e linear com hemoglobina e com contagem de hemácias, sendo mais bem associado com hemoglobina do que com o número de hemácias, pois os pontos estão mais agrupados ao redor da linha pontilhada. Para os leucócitos a dispersão é grande.

Parece possível assumir que existe uma relação aproximadamente linear entre hematócrito (HT) e hemoglobina (HB) e entre hematócrito e hemácia (HEM) e, portanto, a relação pode ser adequadamente analisada por correlação. Entre hematócrito e leucócito a relação linear é duvidosa.

Onde estamos?

Ainda que o interesse seja prever o valor médio da variável dependente a partir de uma variável independente conhecida através de uma regressão linear simples (a veremos, adiante), não há sentido neste procedimento para variáveis que não estejam sequer relacionadas. O primeiro passo, portanto, é verificar se, e com qual intensidade, as variáveis são associadas linearmente. O caminho percorrido no diagrama, portanto, é:

Três coeficientes de correlação

É comum encontrarmos nos softwares estatísticos 3 coeficientes de correlação:

Coeficiente de correlação de Pearson

Karl Pearson (27 de março de 1857 a 27 de abril de 1936) era um matemático e bioestatístico inglês. A ele se atribui o estabelecimento da disciplina de estatística matemática e a fundação do primeiro departamento de estatística do mundo, no University College, London em 1911.

Coeficiente de correlação de Spearman

Charles Edward Spearman, (10 de setembro de 1863 - 17 de setembro de 1945) foi um psicólogo inglês conhecido pelo trabalho em estatística, como pioneiro na análise fatorial e pelo coeficiente de correlação de Spearman. Ele também fez um trabalho seminal em modelos de inteligência humana, incluindo sua teoria de que resultados díspares de testes cognitivos refletem um único fator de inteligência geral, cunhando o termo fator g.

\(\tau\) de Kendall

Sir Maurice George Kendall (6 de setembro de 1907 - 29 de março de 1983) foi um estatístico britânico, amplamente conhecido por sua contribuição para a estatística. Desenvolveu o coeficiente de correlação \(\tau\) (letra grega, tau) de Kendall em 1938 para mensurar a proporção de diferença entre pares concordantes e discordantes com empates. É uma medida de desordem (desarray) e não será discutida neste contexto.

r de Pearson

  • Também chamado de:
    • Coeficiente de correlação produto-momento de Pearson (CCPMP)
    • Coeficiente de correlação de Pearson (CCP)
    • r de Pearson
    • r
  • Avalia o grau de linearidade de duas variáveis intervalares sob as suposições de:
    • Independência entre os pares de observações (condição necessária)
  • O teste da hipótese nula e o intervalo de confiança computados sob as suposições de:
    • Independência entre os pares de observações (condição necessária)
    • Binormalidade (condição suficiente)

O teste da hipótese nula e o intervalo de confiança obtidos por reamostragem (bootstrapping) não necessitam da suposição de binormalidade.

Observe: coeficiente de correlação (r) de Pearson é adequado apenas para relações com formato de reta.

Princípio do r de Pearson

Admita duas variáveis intervalares como, por exemplo, estatura e massa corpórea, observadas em 9 pessoas (implementado em demo_EstMCT.R):

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
cat("estaturas:",estatura,"\n")
cat("massas corporais:",massa,"\n")
sunflowerplot(estatura, massa, 
              xlab = "Estatura (cm)", ylab="Massa corporal (kg)",
              pch=21, col="black")
lines(lowess(estatura, massa), lty=2)

obtendo-se:

estaturas: 165 169 171 172 173 174 176 178 180 
massas corporais: 61 73 68 74 65 82 69 81 83

Aparentemente há uma relação crescente, aproximadamente linear da massa corporal em função da estatura (i.e., pessoas mais altas são, também, mais pesadas e, vice-versa, esperamos que as mais leves sejam as mais baixas).

Para verificar se existe correlação testamos a hipótese da nulidade de \(\rho\) (do qual r é o estimador) para concluirmos se, populacionalmente, a correlação existe:

\[H_0: \rho = 0\] \[H_1: \rho \ne 0\] \[\alpha=0.05\]

Calculamos r e testamos \(H_0\) com a função cor.test (esta função assume \(\alpha=0.05\) por default), como no exemplo fornecido em Correlacao_r_de_Pearson.R:

# Correlacao_r_de_Pearson.R

# valores 
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
correlacao <- cor.test(estatura, massa) 
print(correlacao)

    Pearson's product-moment correlation

data:  estatura and massa
t = 2.8568, df = 7, p-value = 0.02445
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1356665 0.9398558
sample estimates:
     cor 
0.733684

O valor de r é dado pela covariância entre as variáveis em relação ao produto de seus desvios-padrão:

\[ r = \dfrac{\text{cov}(x,y)}{\text{sd}(x) \;\text{sd}(y)} \]

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
covariancia <- cov(estatura,massa)
sd_estatura <- sd(estatura)
sd_massa <- sd(massa)
r <- covariancia / (sd_estatura * sd_massa)
cat("r de Pearson = ", round(covariancia,3), "/(",
    round(sd_estatura,3), "*", round(sd_massa,3), ") = ", r, "\n", sep="")
r de Pearson = 26.514/(4.595*7.865) = 0.733684
print(cor.test(estatura, massa)$estimate)
     cor 
0.733684
Neste exemplo obtivemos a mesma estimativa pontual que cor.test fornece. A função cor.test, porém, além desta estimativa pontual de r=0.7337, calcula o intervalo de confiança 95%, \(\text{IC}95\%\)=[0.1357, 0.9399], largo por conta do número reduzido de pares, mas mesmo assim com o valor 0 fora do intervalo, a estatística de teste \(t\) e seu valor-\(p\) correspondente. Como \(p\)=0.0245 é menor que \(\alpha\), rejeitamos \(H_0\) e assumimos que as variáveis \(x\) e \(y\) estão correlacionadas.

Significado de r de Pearson

A correlação é adimensional e adirecional, servindo para medir a intensidade da associação entre duas variáveis.

Apesar de utilizarmos cor.test com os valores brutos, o r de Pearson é um coeficiente padronizado: é a intensidade da associação das duas variáveis padronizadas que é computada para r. Padronizar é computar o escore z, subtraindo-se a média e dividindo-se todos os valores pelo desvio-padrão. Considerando a estatura no eixo \(x\) e a massa corporal no eixo \(y\), temos:

\[ z_{x,i} = \dfrac{x_i - \bar{x}}{s_x} \]

e

\[ z_{y,i} = \dfrac{y_i - \bar{y}}{s_y} \]

Os respectivos escores \(z\) são adimensionais e a forma das respectivas distribuições não é afetada.

Padronização é um procedimento de transformação linear que substitui valores com unidades de medida pelos seus respectivos escores \(z\). Como a distribuição normal padronizada usa escores \(z\), é um erro comum as pessoas pensarem que a padronização ‘normaliza’ uma distribuição. Não é assim: o que é expresso em valores \(z\) tem média igual a 0 e desvio-padrão igual a 1, qualquer que seja a distribuição. A distribuição normal padronizada vem de uma distribuição normal não padronizada, apenas retendo sua forma original.

Usando o R como laboratório, podemos verificar que a padronização não altera o formato da distribuição (demo_padronizaMCT.R):

# variavel nao padronizada
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
d.massa <- density(massa)
xfit<-seq(50,100,length=100)
yfit<-dnorm(xfit,mean=mean(massa),sd=sd(massa))
ymax <- max(c(d.massa$y,yfit),na.rm=TRUE)

# padronizacao usando a função scale()
massa.z <- scale(massa)
d.massa.z <- density(massa.z)
xfit.z<-seq(-4,4,length=100)
yfit.z<-dnorm(xfit.z,mean=0,sd=1)
ymax.z <- max(c(d.massa.z$y,yfit.z),na.rm=TRUE)

o.par <- par(no.readonly=TRUE)
m <- matrix(1:2, nrow=1, ncol=2, byrow=TRUE)
layout(m)
plot(d.massa, main="Não padronizada", 
     xlab="Massa (kg)",
     ylab="Densidade",
     ylim=c(0,ymax))
lines(xfit,yfit,lty=2)
plot(d.massa.z, main="Padronizada", 
     xlab="Massa (z)",
     ylab="Densidade",
     ylim=c(0,ymax.z))
lines(xfit.z,yfit.z,lty=2)
par(o.par)

obtendo-se:

Observe a mudança da escala no eixo das abscissas.


Podemos verificar que padronização não muda a correlação. Com o código em Correlacao_Original_e_Padronizada.R executamos a correlação com os valores originais e com os valores padronizados.

Compare:

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)

cat("\nValores originais:\n")
cat("estaturas:",estatura,"\n")
cat("massas corporais:",massa,"\n")
cat("\nCorrelação com as variáveis originais:\n")
correlacao <- cor.test(estatura, massa)
print(correlacao)

# padronizando x e y
z_estatura <- scale(estatura)
z_massa <- scale(massa)
cat("\nValores padronizados:\n")
cat("estatura padronizada:",round(z_estatura,3),"\n")
cat("massa corporal padronizada:",round(z_massa,3),"\n")
cat("\nCorrelação com as variáveis padronizadas:\n")
correlacao_padrao <- cor.test(z_estatura, z_massa)
print(correlacao_padrao)

Valores originais:
estaturas: 165 169 171 172 173 174 176 178 180 
massas corporais: 61 73 68 74 65 82 69 81 83 

Correlação com as variáveis originais:

    Pearson's product-moment correlation

data:  estatura and massa
t = 2.8568, df = 7, p-value = 0.02445
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1356665 0.9398558
sample estimates:
     cor 
0.733684 


Valores padronizados:
estatura padronizada: -1.765 -0.895 -0.459 -0.242 -0.024 0.193 0.629 1.064 1.499 
massa corporal padronizada: -1.512 0.014 -0.622 0.141 -1.003 1.158 -0.494 1.031 1.286 

Correlação com as variáveis padronizadas:

    Pearson's product-moment correlation

data:  z_estatura and z_massa
t = 2.8568, df = 7, p-value = 0.02445
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1356665 0.9398558
sample estimates:
     cor 
0.733684

Os valores de r e de seu intervalo de confiança 95% são os mesmos quando computamos com os valores originais ou com seus respectivos escores \(z\).

O valor de r é dado por:

\[ r = \dfrac{\text{cov}(x,y)}{\text{sd}(x) \;\text{sd}(y)} \]

No entanto, quando as variáveis são padronizadas, os desvios-padrão são iguais a 1. Consequentemente, percebemos que r também é a covariância entre as variáveis padronizadas:

\[ r = \text{cov}(z_x,z_y) \]

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
mean_estatura <- round(mean(estatura),3)
sd_estatura <- round(sd(estatura),3)
mean_massa <- round(mean(massa),3)
sd_massa <- round(sd(massa),3)
cat("Variáveis não padronizadas\n", sep="")
cat("\tmédia (MCT) = ",mean_massa,"\n", sep="")
cat("\tst.dev. (MCT) = ",sd_massa,"\n", sep="")
cat("\tmédia (EST) = ",mean_estatura,"\n", sep="")
cat("\tst.dev. (EST) = ",sd_estatura,"\n", sep="")
covariancia <- cov(estatura,massa)
cat("cov(EST,MCT) = ",covariancia,"\n", sep="")
r <- cov(estatura,massa)/(sd_massa*sd_estatura)
cat("r de Pearson = ",r,"\n", sep="")

cat("\n")
cat("Variáveis padronizadas\n", sep="")
z_estatura <- (estatura-mean_estatura)/sd_estatura
z_massa <- (massa-mean_massa)/sd_massa
mean_z_estatura <- round(mean(z_estatura),3)
sd_z_estatura <- round(sd(z_estatura),3)
mean_z_massa <- round(mean(z_massa),3)
sd_z_massa <- round(sd(z_massa),3)
cat("\tmédia (MCT) = ",mean_z_massa,"\n", sep="")
cat("\tst.dev. (MCT) = ",sd_z_massa,"\n", sep="")
cat("\tmédia (EST) = ",mean_z_estatura,"\n", sep="")
cat("\tst.dev. (EST) = ",sd_z_estatura,"\n", sep="")
z_r <- cov(z_estatura,z_massa)
cat("cov(EST,MCT) = r de Pearson = ",z_r,"\n", sep="")
Variáveis não padronizadas
    média (MCT) = 72.889
    st.dev. (MCT) = 7.865
    média (EST) = 173.111
    st.dev. (EST) = 4.595
cov(EST,MCT) = 26.51389
r de Pearson = 0.7336505

Variáveis padronizadas
    média (MCT) = 0
    st.dev. (MCT) = 1
    média (EST) = 0
    st.dev. (EST) = 1
cov(EST,MCT) = r de Pearson = 0.7336505
Note que, neste código, fizemos a padronização explicitamente, subtraindo a média e dividindo pelo desvio-padrão. No código anterior usamos a função scale, que faz o mesmo de forma mais simples.

Valores de r de Pearson

O coeficiente r varia de -1 (correlação negativa perfeita) a +1 (correlação positiva perfeita), passando pelo 0 (ausência de correlação).

Graficamente, quanto mais bem alinhados estiverem os pontos, mais \(|r|\) se aproxima de 1. Valores negativos indicam que o valor de uma das variáveis decresce quando a outra cresce; positivos quando ambos os valores decrescem ou crescem concordantemente:

\(r = -1\) \(r = -0.8\) \(r = -0.4\) \(r = -0.2\)
\(r = 0\)

\(r = 0.2\) \(r = 0.4\) \(r = 0.8\) \(r = 1\)

É possível treinar-se para acertar o valor de r observando um gráfico desde que as variáveis sejam padronizadas antes; caso contrário, a distorção causada pela unidades de medida das variáveis \(x\) e \(y\) atrapalharão.

Uma implementação do jogo está no APEx. A partir do menu principal acesse:

Estatística: jogos

e escolha o jogo r in R.

r de Pearson no exemplo de Gestante

Executando Gestantes_rPearson.R, revisitamos os gráficos de hemoglobina, hemácias e leucócitos em função do hematócrito (agora padronizados), verificando suas respectivas correlações:

source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")

col_HB <- friendlycolor(30) # tijolo
pch_HB <- 22
col_HEM <- friendlycolor(28) # bordo
pch_HEM <- 23
col_LEUC <- friendlycolor(9) # azul
pch_LEUC <- 24

Gestantes <- readRDS("Gestante.rds")

# HT x HB
with(Gestantes,correg(HT, HB, method="pearson",
                      xlab="Hematcrito (%)", ylab="Hemoglobina (mg/dl)",
                      col=col_HB, bg="transparent", pch=pch_HB))

# HT x HEM
with(Gestantes,correg(HT, HEM, method="pearson",
                      xlab="Hematócrito (%)", ylab="Hemácias (milhão/mm3)",
                      col=col_HEM, bg="transparent", pch=pch_HEM))

# HT x LEUC
with(Gestantes,correg(HT, LEUC, method="pearson",
                      xlab="Hematócrito (%)", ylab="Leucócitos (milhar/mm3)",
                      col=col_LEUC, bg="transparent", pch=pch_LEUC))

             Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1      Hematcrito (%)   34 2.5899   40      36     34      37   39 40    0
2 Hemoglobina (mg/dl) 11.9 1.0446   14    11.2   12.2    13.2 13.6 40    0
lambda =  3.86866005069142 
Assuming:
    - Explanatory variable: Hematcrito (%);
    - Dependent variable: Hemoglobina (mg/dl).

-----------
Correlation
-----------

    Pearson's product-moment correlation

data:  Hematcrito (%) and Hemoglobina (mg/dl)
t = 10.604, df = 38, p-value = 6.523e-13
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7568524 0.9265212
sample estimates:
      cor 
0.8645336 

               Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1       Hematócrito (%)   34 2.5899   40      36     34      37   39 40    0
2 Hemácias (milhão/mm3) 3.65 0.3831 4.53    4.89   4.15    4.32 4.57 40    0
lambda =  13.8068135732994 
Assuming:
    - Explanatory variable: Hematócrito (%);
    - Dependent variable: Hemácias (milhão/mm3).

-----------
Correlation
-----------

    Pearson's product-moment correlation

data:  Hematócrito (%) and Hemácias (milhão/mm3)
t = 5.8153, df = 38, p-value = 1.02e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4765742 0.8220069
sample estimates:
      cor 
0.6862105 

                 Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1         Hematócrito (%)   34 2.5899   40      36     34      37   39 40    0
2 Leucócitos (milhar/mm3)  8.4 2.0243  8.3      11    8.7      10  7.1 40    0
lambda =  2.60777879175739 
Assuming:
    - Explanatory variable: Hematócrito (%);
    - Dependent variable: Leucócitos (milhar/mm3).

-----------
Correlation
-----------

    Pearson's product-moment correlation

data:  Hematócrito (%) and Leucócitos (milhar/mm3)
t = 0.70829, df = 38, p-value = 0.4831
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2046370  0.4110423
sample estimates:
      cor 
0.1141489

Com este código, os valores são padronizados e é adicionada uma elipse de predição de 95% com seu eixo principal (não é uma reta de regressão!). Compare esta elipse com as nuvens de pontos que foram vistas na sessão valores do r de Pearson.

Observa-se, portanto, que a concentração de hemoglobina e a contagem de hemácias têm correlação significante com o hematócrito para \(\alpha=0.05\), mas a contagem de leucócitos não tem.

Há mais de uma forma para se chegar a esta decisão estatística. Recorde que testa-se \(H_0: \rho=0\). Podemos rejeitar a hipótese nula quando o valor \(p < \alpha = 5\%\) ou pelo intervalo de confiança para r, quando este não inclui o zero.

r de Pearson com bootstrapping

Elaboramos Gestantes_rPearson_boot_z.R. Utilizamos a mesma função (correg) adicionando o parâmetro B. As diferenças são a exibição de várias elipses, cada uma vinda de uma reamostragem, e a substituição das saídas textuais de cor.test pelos valores obtidos por bootstrapping. Para \(B=1e4\) reamostragens, obtém-se:

source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")
source("eiras.correg.R")

Gestantes <- readRDS("Gestante.rds")

col_HB <- friendlycolor(30) # tijolo
pch_HB <- 22
col_HEM <- friendlycolor(28) # bordo
pch_HEM <- 23
col_LEUC <- friendlycolor(9) # azul
pch_LEUC <- 24
B <- 1e4

# HT x HB
correg(Gestantes$HT, Gestantes$HB, method="pearson", B=B,
              col=col_HB, bg=paste(col_HB,"50",sep=""), pch=pch_HB,
              xlab="Hematócrito (%)", ylab="Hemoglobina (mg/dl)")

# HT x HEM
correg(Gestantes$HT, Gestantes$HEM, method="pearson", B=B,
       col=col_HEM, bg=paste(col_HEM,"50",sep=""), pch=pch_HEM,
       xlab="Hematócrito (%)", ylab="Hemácias (milhão/mm3)")

# HT x LEUC
correg(Gestantes$HT, Gestantes$LEUC, method="pearson", B=B,
       col=col_LEUC, bg=paste(col_LEUC,"50",sep=""), pch=pch_LEUC,
       xlab="Hematócrito (%)", ylab="Leucócitos (milhar/mm3)")

             Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1     Hematócrito (%)   34 2.5899   40      36     34      37   39 40    0
2 Hemoglobina (mg/dl) 11.9 1.0446   14    11.2   12.2    13.2 13.6 40    0
lambda =  3.86866005069142 
Assuming:
    - Explanatory variable: Hematócrito (%);
    - Dependent variable: Hemoglobina (mg/dl).

-----------
Correlation
-----------

    Bootstrapp r, 10000 resamples

data:  Hematócrito (%) and Hemoglobina (mg/dl)

alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6938923 0.9424315
sample estimates:
r (bootstrap) 
     0.866473 

               Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1       Hematócrito (%)   34 2.5899   40      36     34      37   39 40    0
2 Hemácias (milhão/mm3) 3.65 0.3831 4.53    4.89   4.15    4.32 4.57 40    0
lambda =  13.8068135732994 
Assuming:
    - Explanatory variable: Hematócrito (%);
    - Dependent variable: Hemácias (milhão/mm3).

-----------
Correlation
-----------

    Bootstrapp r, 10000 resamples

data:  Hematócrito (%) and Hemácias (milhão/mm3)

alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4020068 0.8459592
sample estimates:
r (bootstrap) 
    0.6862866 

                 Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1         Hematócrito (%)   34 2.5899   40      36     34      37   39 40    0
2 Leucócitos (milhar/mm3)  8.4 2.0243  8.3      11    8.7      10  7.1 40    0
lambda =  2.60777879175739 
Assuming:
    - Explanatory variable: Hematócrito (%);
    - Dependent variable: Leucócitos (milhar/mm3).

-----------
Correlation
-----------

    Bootstrapp r, 10000 resamples

data:  Hematócrito (%) and Leucócitos (milhar/mm3)

alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1334425  0.3939304
sample estimates:
r (bootstrap) 
    0.1182038

Na versão com bootstrapping, a função correg mostra uma sombra das elipses variantes, sugerindo sua dispersão. A decisão sobre a correlação é feita pelo intervalo de confiança, pois não existe valor \(p\) calculado com este método. O teorema central do limite é aplicável neste caso, como se verifica pela distribuição dos valores de r apresentados.

Ainda assim, compare com os resultados convencionais e observe que a conclusão é a mesma: somente a concentração de hemoglobina e a contagem de hemácias têm correlação com o hematócrito.

Tamanho de efeito

  • \(|r|\), o valor absoluto de r, é a intensidade da linearidade.
  • \(r^2\), o quadrado de r, é a proporção da variância compartilhada.

Ambos são medidas de tamanho de efeito:

  • São adimensionais
  • Variam entre 0 e 1
  • Não dependem do tamanho da amostra

Ellis, 2010

\(s\) de Spearman

Para distingui-lo do r de Pearson, denotaremos o coeficiente de correlação de Spearman como \(s\). A função em R é a mesma, cor.test, computando \(s\) quando o parâmetro é method="spearman".

O \(s\) de Spearman também serve para variáveis ordinais. Neste caso, a representação gráfica por meio do gráfico de dispersão não terá sentido.

Da mesma forma que o r de Pearson, \(s\) varia de \(-1\) a \(+1\).

A informação provida por este coeficiente é sobre quanto o (de)crescimento dos valores é monotônico (i.e., o valor de uma variável cresce sempre ou decresce sempre quando o valor da outra aumenta). Assim:

\[ \begin{cases} H_0:&\text{as duas variáveis têm relação monotônica}\\ H_1:&\text{as duas variáveis não têm relação monotônica} \end{cases} \]

O exemplo hipotético com estatura e massa corpórea implementado em demo_r_de_Spearman.R, com o método de Spearman, mostra:

# demo_r_de_Spearman_boot.R

# disable warnings
options(warn=-1)

alfa <- 0.05
B <- 0
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")

# cores
col <-  friendlycolor(31) # preto
bg <- friendlycolor(24) # amarelo
pch <- 21 # circulo

# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
cat("estatura:",estatura,"\n")
cat("massa corporal:",massa,"\n")
lst <- correg(estatura, massa, 
              method = "spearman",
              main="",
              xlab="Estatura (cm)", 
              ylab="Massa corporal (kg)", 
              jitter=0,
              col=col, bg=bg, pch=pch,
              alpha=alfa,
              B=B)
# arrows
for (i in 1:(length(estatura)-1))
{
  arrows(estatura[i],massa[i],
         estatura[i+1],massa[i+1], length=0.15, angle=20)
}
# enable warnings
options(warn=0)
estatura: 165 169 171 172 173 174 176 178 180 
massa corporal: 61 73 68 74 65 82 69 81 83 
165 169 171 172 173 174 176 178 180 
 61 73 68 74 65 82 69 81 83 

             Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1       Estatura (cm)  172 4.5947  165     169    171     173  174 9    0
2 Massa corporal (kg)   74 7.8652   61      73     68      65   82 9    0
lambda =  0.855927523022535 
Assuming:
    - Explanatory variable: Estatura (cm);
    - Dependent variable: Massa corporal (kg).

-----------
Correlation
-----------

    Spearman's rank correlation rho

data:  Estatura (cm) and Massa corporal (kg)
S = 36, p-value = 0.04325
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho 
0.7

As setas foram adicionadas para evidenciar o que torna \(|s| < 1\): com o crescimento da estatura (\(x\)), os valores de massa corporal (\(y\)) ora crescem, ora descrescem.

Implementamos, também, uma versão em bootstrapping em demo_r_de_Spearman_boot.R, que tem a vantagem de permitir a decisão através do intervalo de confiança (ao preço de não ter valor \(p\)) executado com.

# demo_r_de_Spearman_boot.R

# disable warnings
options(warn=-1)

alfa <- 0.05
B <- 1e4

source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")

# cores
col <-  friendlycolor(31) # preto
bg <- friendlycolor(24) # amarelo
pch <- 21 # circulo

# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
cat("estatura:",estatura,"\n")
cat("massa corporal:",massa,"\n")
lst <- correg(estatura, massa, 
              method = "spearman",
              main="",
              xlab="Estatura (cm)", 
              ylab="Massa corporal (kg)", 
              jitter=0,
              col=col, bg=bg, pch=pch,
              alpha=alfa,
              B=B)

# enable warnings
options(warn=0)
estatura: 165 169 171 172 173 174 176 178 180 
massa corporal: 61 73 68 74 65 82 69 81 83 

             Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1       Estatura (cm)  172 4.5947  165     169    171     173  174 9    0
2 Massa corporal (kg)   74 7.8652   61      73     68      65   82 9    0
lambda =  0.855927523022535 
Assuming:
    - Explanatory variable: Estatura (cm);
    - Dependent variable: Massa corporal (kg).

-----------
Correlation
-----------

    Bootstrapp s, 10000 resamples

data:  Estatura (cm) and Massa corporal (kg)

alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.01785714 1.00000000
sample estimates:
s (bootstrap) 
    0.7142857

O teorema central do limite parece razoavelmente atendido (poderia estar melhor).

O coeficiente de Spearman é máximo quando o crescimento é monotônico (\(s=1\)) ou quando o decrescimento é monotônico (\(s=-1\)). Por exemplo, executando: demo_r_de_Spearman_boot_monotonico.R obtemos:

# demo_r_de_Spearman_boot_monotonico.R

# disable warnings
options(warn=-1)

alfa <- 0.05
B <- 0

source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.correg.R")

# cores
col <- friendlycolor(31) # preto
bg <- friendlycolor(24) # amarelo
pch <- 21 # circulo

# valores
set.seed(27)
x <- 1:20
y <- c(runif(1,1,2))
for (i in 2:length(x))
{
  y <- c(y,runif(1,y[i-1]+0.01,y[i-1]+i^2))
}
cat("x:",x,"\n")
cat("y:",round(y,3),"\n")
correg(x, y,
       alpha=alfa, B=B, 
       method="spearman",
       xlab="x", ylab="y", 
       col=col, bg=bg, pch=pch)
lines(x,y,type="b")

# enable warnings
options(warn=0)
x: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
y: 1.972 2.316 10.182 15.456 21.021 35.486 39.048 39.215 50.328 69.527 144.079 266.425 351.394 496.108 693.31 894.604 1178.866 1208.258 1422.695 1441.131 
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
 1.97175 2.315943 10.18203 15.45644 21.02111 35.48643 39.04813 39.21492 50.32816 69.52715 144.0785 266.4249 351.3942 496.108 693.31 894.6037 1178.866 1208.258 1422.695 1441.131 

  Variable    Mean       Sd   Min. 1st Qu. Median          3rd Qu.             Max.  n NA's
1        x       4   5.9161      1       2      3                5                6 20    0
2        y 15.4564 522.7931 1.9718  2.3159 10.182 21.0211078875745 35.4864272550144 20    0
lambda =  0.0101364986583543 
Assuming:
    - Explanatory variable: x;
    - Dependent variable: y.

-----------
Correlation
-----------

    Spearman's rank correlation rho

data:  x and y
S = 0, p-value = 5.976e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho 
  1

Basta que algum valor não seja maior que o anterior para \(r < 1\). Para isto, modifico um valor do exemplo anterior com:

y[14] <- y[13]-100 # criando valor menor que o anterior

Obtendo:

# demo_r_de_Spearman_boot_quasemonotonico.R

# disable warnings
options(warn=-1)

alfa <- 0.05
B <- 0

source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.correg.R")

# cores
col <-  friendlycolor(31) # preto
bg <- friendlycolor(24) # amarelo
pch <- 21 # circulo

# valores
set.seed(27)
x <- 1:20
y <- c(runif(1,1,2))
for (i in 2:length(x))
{
  y <- c(y,runif(1,y[i-1]+0.01,y[i-1]+i^2))
}
y[14] <- y[13]-100 # criando valor menor que o anterior
cat("x:",x,"\n")
cat("y:",round(y,3),"\n")
correg(x, y,
       alpha=alfa, B=B, 
       method="spearman",
       jitter=0,
       xlab="x", ylab="y", 
       col=col, bg=bg, pch=pch)
lines(x,y,type="b")
points(x[14],y[14],
       col="black",bg=friendlycolor(30),pch=pch,
       cex=1.4) # exibindo o valor reduzido

# enable warnings
options(warn=0)
x: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
y: 1.972 2.316 10.182 15.456 21.021 35.486 39.048 39.215 50.328 69.527 144.079 266.425 351.394 251.394 693.31 894.604 1178.866 1208.258 1422.695 1441.131 
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
 1.97175 2.315943 10.18203 15.45644 21.02111 35.48643 39.04813 39.21492 50.32816 69.52715 144.0785 266.4249 351.3942 251.3942 693.31 894.6037 1178.866 1208.258 1422.695 1441.131 

  Variable    Mean      Sd   Min. 1st Qu. Median          3rd Qu.             Max.  n NA's
1        x       4  5.9161      1       2      3                5                6 20    0
2        y 15.4564 523.758 1.9718  2.3159 10.182 21.0211078875745 35.4864272550144 20    0
lambda =  0.00993023255371815 
Assuming:
    - Explanatory variable: x;
    - Dependent variable: y.

-----------
Correlation
-----------

    Spearman's rank correlation rho

data:  x and y
S = 6, p-value = 6.135e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.9954887

Neste código e nos dois seguintes deixamos B=0, significando que foram executados sem bootstrapping. Caso deseje, altere este parâmetro (geralmente recomenda-se \(1e4\) ou mais reamostragens para um resultado confiável) para comparar com esta versão tradicional.

O valor será de \(r=-1\) para o decrescimento monotônico:

# demo_r_de_Spearman_boot_monotonicoD.R

# disable warnings
options(warn=-1)

alfa <- 0.05
B <- 0

source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.correg.R")

# cores
col <- friendlycolor(31) # preto
bg <- friendlycolor(4) # violeta
pch <- 21 # circulo

# valores
set.seed(27)
x <- 1:20
y <- c(runif(1,1600,1700))
for (i in 2:length(x))
{
  y <- c(y,runif(1,y[i-1]-i^2,y[i-1]+0.01))
}
cat("x:",x,"\n")
cat("y:",round(y,3),"\n")
correg(x, y,
       alpha=alfa, B=B, 
       method="spearman",
       jitter=0,
       xlab="x", ylab="y", 
       col=col, bg=bg, pch=pch)
lines(x,y,type="b")

# enable warnings
options(warn=0)
x: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
y: 1697.175 1693.511 1692.384 1681.655 1662.215 1640.678 1595.231 1531.388 1461.494 1380.687 1334.24 1312.594 1228.563 1177.282 1149.491 1094.791 1090.062 795.446 648.886 267.313 
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
 1697.175 1693.511 1692.384 1681.655 1662.215 1640.678 1595.231 1531.388 1461.494 1380.687 1334.24 1312.594 1228.563 1177.282 1149.491 1094.791 1090.062 795.446 648.8857 267.3127 

  Variable      Mean       Sd     Min.   1st Qu.    Median          3rd Qu.             Max.  n NA's
1        x         4   5.9161        1         2         3                5                6 20    0
2        y 1681.6555 391.1633 1697.175 1693.5109 1692.3845 1662.21456341912 1640.67791575092 20    0
lambda =  0.0140842894040421 
Assuming:
    - Explanatory variable: x;
    - Dependent variable: y.

-----------
Correlation
-----------

    Spearman's rank correlation rho

data:  x and y
S = 2660, p-value = 5.976e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho 
 -1

Exemplo: Gestante

Exemplificamos em Gestantes_matrizcorrelacoes.R a montagem de uma matriz de correlações e sugerimos a função GGally::ggcorr que auxilia com informação visual. Além de hematócrito, hemoglobina, hemácias e leucócitos, adicionamos mais variáveis (idade, ácido fólico e vitamina B12). Resulta em:

# Gestantes_matrizcorrelacoes.R

options(warn=-1)

Gestantes <- readRDS("Gestante.rds")
mc <- data.frame(Gestantes$IDADE, 
                 Gestantes$HT, 
                 Gestantes$HB, 
                 Gestantes$HEM, 
                 Gestantes$LEUC, 
                 Gestantes$FOLICO, 
                 Gestantes$B12)
names(mc) <- c("Idade","HT","HB","HEM","LEUC","FOLICO","B12")
print(head(mc))
cat("\n...\n")
print(tail(mc, addrownums = FALSE, n=2L))
cat("\nMatriz de correlações:\n")
print(cor(mc), digits=1) # matriz de correlacoes
# grafico da matriz
print(GGally::ggcorr(mc,
                     nbreaks = 6,
                     digits = 3,
                     label = TRUE,
                     label_size = 4,
                     color = "#888888"))


options(warn=0)
  Idade HT   HB  HEM LEUC FOLICO B12
1    24 40 14.0 4.53  8.3   11.6 324
2    34 36 11.2 4.89 11.0   12.8 740
3    27 34 12.2 4.15  8.7   13.9 360
4    37 34 11.9 3.65  8.4   13.3 291
5    26 37 13.2 4.32 10.0    6.6 303
6    28 39 13.6 4.57  7.1    8.0 292

...
   Idade HT   HB  HEM LEUC FOLICO B12
39    19 37 13.8 4.25  9.6    2.9 235
40    34 39 13.9 4.20 10.6    6.9 208

Matriz de correlações:
        Idade     HT    HB  HEM  LEUC FOLICO   B12
Idade   1.000  0.007 -0.05 0.01  0.21   0.05 -0.03
HT      0.007  1.000  0.86 0.69  0.11   0.25 -0.05
HB     -0.049  0.865  1.00 0.59  0.04   0.08 -0.11
HEM     0.014  0.686  0.59 1.00  0.17   0.26  0.26
LEUC    0.212  0.114  0.04 0.17  1.00  -0.05 -0.12
FOLICO  0.046  0.253  0.08 0.26 -0.05   1.00  0.30
B12    -0.029 -0.054 -0.11 0.26 -0.12   0.30  1.00

Variância compartilhada

A correlação computa o r de Pearson quantifica o grau de associação linear entre duas variáveis intervalares. No entanto, não traçamos retas de regressão linear entre as variáveis; o máximo que fizemos foi utilizar a função que construímos, ellipseaxis, que exibe uma elipse e uma linha em seu eixo maior servindo-nos para auxiliar o julgamento sobre associação linear.

Também comentamos que a correlação é adirecional (não falamos em variável explicativa, VE, nem dependente ou de desfecho, VD), pois os resultados seriam os mesmos se invertêssemos os eixos \(x\) e \(y\). Também é adimensional, pois o procedimento não necessita das unidades de medida, e seu resultado é o mesmo quando as duas variáveis são padronizadas.

O que, então, a correlação mede?

Dissemos, acima: “o valor absoluto do r é a intensidade da linearidade; seu quadrado é a proporção da variância compartilhada.”

Qual o significado disto?

É fácil ver através de um diagrama de Venn. Como as variáveis são compartilhadas (i.e., todas têm média igual a zero e variância igual a um), todas são representadas por círculos de dimensões iguais. O \(r^2\) é a proporção da intersecção entre as duas variáveis (dois círculos completamente sobrepostos equivale a \(r^2 = 1\); dois círculos isolados representam \(r^2 = 0\)). No exemplo do hematócrito contra hemoglobina, hemácias e leucócitos, os valores de r e \(r^2\) correspondem aproximadamente a:

Como vimos acima, apenas a correlação entre hematócrito e a contagem de leucócitos não foi significante.

Níveis de mensuração e outros tipos de correlação

Existem vários tipos de correlação além das duas que exploramos aqui. A correlação de Pearson é aplicável para variáveis intervalares. A correlação de Spearman usa a mesma fórmula, adaptada para variáveis ordinais. Das outras, algumas delas são também correlações de Pearson, como a ponto-bisserial e phi, quando há variáveis dicotômicas autênticas envolvidas. As demais são para situações mais complexas, fora do escopo deste capítulo.

Revelle, 2014

Testes para uma correlação

Por default, cor.test usa a hipótese nula de que não há associação populacional entre duas variáveis (\(\rho = 0\)).

Como fizemos em outros assuntos, em publicações muitas vezes só temos as medidas-resumo. O código R implementado em CorrelacaoPearsonUmaCorrelacao.R permite que se verifique, com dado tamanho da amostra, se um determinado valor de r é significantemente diferente de zero. Por exemplo, para \(n=200\) e \(r=0.5\), as hipóteses são:

\[ \begin{cases} H_0:& \rho = 0\\ H_1:& \rho \ne 0 \end{cases}\\ \alpha = 0.05 \]

A saída é:

alfa <- 0.05
n <- 200
r <- 0.5  
# H0: ro = 0 vs. H1: ro != 0
# https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
cat("Teste bilateral:\n")
cat("\tH0: ro = 0 vs. H1: ro <> 0\n")
t <- r*sqrt((n-2)/(1-r^2))
p <- 2*pt(-abs(t), df=n-2)
cat("t(",n-2,") = ",t," , p = ",p,"\n",sep="")
correlacao <- DescTools::CorCI(r, n=n, conf.level = 1-alfa, 
                               alternative = "two.sided")
print(correlacao)
cat("\n")
cat("Teste unilateral a direita:\n")
cat("\tH0: ro = 0 vs. H1: ro > 0\n")
p <- 1-pt(t, df=n-2)
cat("t(",n-2,") = ",t," , p = ",p,"\n",sep="")
correlacao <- DescTools::CorCI(r, n=n, conf.level = 1-alfa, 
                               alternative = "greater")
print(correlacao)
cat("\n")
cat("Teste unilateral a esquerda:\n")
cat("\tH0: ro = 0 vs. H1: ro < 0\n")
p <- pt(t, df=n-2)
cat("t(",n-2,") = ",t," , p = ",p,"\n",sep="")
correlacao <- DescTools::CorCI(r, n=n, conf.level = 1-alfa, 
                               alternative = "less")
print(correlacao)
Teste bilateral:
    H0: ro = 0 vs. H1: ro <> 0
t(198) = 8.124038 , p = 4.773812e-14
      cor    lwr.ci    upr.ci 
0.5000000 0.3881878 0.5973056 

Teste unilateral a direita:
    H0: ro = 0 vs. H1: ro > 0
t(198) = 8.124038 , p = 2.38698e-14
      cor    lwr.ci    upr.ci 
0.5000000 0.4070875 1.0000000 

Teste unilateral a esquerda:
    H0: ro = 0 vs. H1: ro < 0
t(198) = 8.124038 , p = 1
      cor    lwr.ci    upr.ci 
 0.500000 -1.000000  0.582671

Observe que, mesmo sem os dados brutos, é possível computar o valor \(p\) e os respectivos intervalos de confiança para os testes bilateral, unilateral à esquerda e unilateral à direita. Neste exemplo, \(r=0.5\) é significantemente diferente de zero no teste bilateral (zero não está contido no intervalo de confiança 95%) e o intervalo não está centrado na estimativa pontual.

No teste unilateral à esquerda não rejeitamos a hipótese nula (a alternativa era dizer que o valor é negativo). No teste unilateral à direita rejeitamos a hipótese nula (então 0.5 é positivo).

Como a função DescTools::CorCI não fornece o valor \(p\), por fórmula podemos obtê-lo (CorrelacaoPearsonUmaCorrelacaoRho0.R):

alfa <- 0.05
n <- 200 
r <- 0.5
cat("\n")
cat("Teste de correlação:\n")
cat("\tH0: r = 0 vs. H1: r <> 0\n",sep="")
z <- sqrt((n-3)/2)*log(1/((1-r)/(1+r)))
p <- 2*pnorm(-abs(z))
cat("z = ",z,", p = ",p, sep="")

cat("\n\nAlternativamente:\n")
t <- r*sqrt((n-2)/(1-r^2))
p <- 2*pt(-abs(t), df=n-2)
cat("t(",n-2,") = ",t," , p = ",p,"\n",sep="")

Teste de correlação:
    H0: r = 0 vs. H1: r <> 0
z = 10.90342, p = 1.110102e-27

Alternativamente:
t(198) = 8.124038 , p = 4.773812e-14

No entanto, outras condições podem ser verificadas. Caso quiséssemos saber se o valor observado deste exemplo, \(r=0.5\), difere de um outro valor qualquer (e.g., \(\rho = 0.6\)):

\[ \begin{cases} H_0:& \rho = 0.6\\ H_1:& \rho \ne 0.6 \end{cases}\\ \alpha = 0.05 \]

Podemos usar o teste implementado em CorrelacaoPearsonUmaCorrelacaoRhoQualquer.R:

alfa <- 0.05
n <- 200
r <- 0.5 
ro <- 0.6
cat("\n")
cat("Teste de correlação:\n")
cat("\tH0: ro = ",ro," vs. H1: ro <> ",ro,"\n",sep="")
z <- sqrt((n-3)/2)*log(((1-ro)/(1+ro))/((1-r)/(1+r)))
p <- 2*pnorm(-abs(z))
cat("z = ",z,", p = ",p, sep="")

Teste de correlação:
    H0: ro = 0.6 vs. H1: ro <> 0.6
z = -2.855163, p = 0.004301474

Quando temos os dados brutos, há diversas funções em R que podem lhe fornecer o intervalo de confiança 95%. Por exemplo, no caso da correlação do Hematócrito com Hemoglobina, Hemácias ou Leucócitos, testando-se:

\[ \begin{cases} H_0:& \rho = 0\\ H_1:& \rho \ne 0 \end{cases}\\ \alpha = 0.05 \]

A função bootES::bootES() encontra o intervalo por bootstrapping, como implementamos em Gestantes_corrBoot.R:

alfa <- 0.05
n <- 200
r <- 0.5 
ro <- 0.6
cat("\n")
cat("Teste de correlação:\n")
cat("\tH0: ro = ",ro," vs. H1: ro <> ",ro,"\n",sep="")
z <- sqrt((n-3)/2)*log(((1-ro)/(1+ro))/((1-r)/(1+r)))
p <- 2*pnorm(-abs(z))
cat("z = ",z,", p = ",p, sep="")
source("Gestantes_corrBoot.R")
Hematócrito e Hemoglobina:
95.00% bca Confidence Interval, 10000 replicates
Stat        CI (Low)    CI (High)   bias        SE          
0.865       0.700       0.943       -0.012      0.065       

Hematócrito e Hemácias:
95.00% bca Confidence Interval, 10000 replicates
Stat        CI (Low)    CI (High)   bias        SE          
0.686       0.420       0.858       -0.017      0.115       

Hematócrito e Leucócitos:
95.00% bca Confidence Interval, 10000 replicates
Stat        CI (Low)    CI (High)   bias        SE          
0.114       -0.158      0.367       0.010       0.134

Comparação entre duas correlações independentes

Padronizar as variáveis e proceder à regressão linear não é apenas uma curiosidade. Pode ajudar em explorar melhor as relações. Observe o que acontece com os dados de Adm2008.xlsx quando exibimos separadamente as medidas de estatura e massa corpórea total (MCT) de homens e mulheres obtidos entre estudantes da Faculdade de Economia, Administração e Contabilidade da Universidade de São Paulo.

Implementamos em CorrelacaoPearsonDuasCorrelacoesIndep.R:

Dados <- readxl::read_excel("Adm2008.xlsx")
DadosF <- Dados[Dados$Sexo=="Feminino",]
DadosM <- Dados[Dados$Sexo=="Masculino",]
nF <- nrow(DadosF)
nM <- nrow(DadosM)

# saida textual 
source("eiras.bartitle.R")

cat(bartitle("Estatística descritiva"))
cat(bartitle("- sexo feminino"))
print(summary(Dados[Dados$Sexo=="Feminino",3:4]))
cat("n: ",nF,"\n",sep="")
cat("Correlacao: ",cor(Dados[Dados$Sexo=="Feminino",3:4])[1,2],"\n",sep="")
cat(bartitle("- sexo masculino"))
print(summary(Dados[Dados$Sexo=="Masculino",3:4]))
cat("n: ",nM,"\n",sep="")
cat("Correlacao: ",cor(Dados[Dados$Sexo=="Masculino",3:4])[1,2],"\n",sep="")

# Elipse de predicao 95%
matriz <- as.matrix(Dados[, 3:4])
n <- nrow(matriz)
car::dataEllipse(matriz[,1], matriz[,2],
                 groups=factor(Dados$Sexo),
                 group.labels=c("F", "M"),
                 col=c("darkorange3","royalblue3"),
                 levels=c(.95),
                 robust=TRUE,
                 main=paste("Elipse de predição de 95%\n",
                            "n = ",n," (Fem=",nF,", Masc=",nM,")",sep=""),
                 xlab="Estatura (m)",
                 ylab="MCT (kg)",
                 xlim=c(1.3, 2.1),
                 ylim=c(25, 120),
                 lwd=1.5, lty=2)

----------------------
Estatística descritiva
----------------------

---------------
- sexo feminino
---------------
    Estatura          MCT       
 Min.   :1.500   Min.   :43.00  
 1st Qu.:1.600   1st Qu.:50.25  
 Median :1.640   Median :53.50  
 Mean   :1.641   Mean   :54.63  
 3rd Qu.:1.698   3rd Qu.:59.75  
 Max.   :1.780   Max.   :70.00  
n: 38
Correlacao: 0.637148

----------------
- sexo masculino
----------------
    Estatura          MCT        
 Min.   :1.560   Min.   : 48.00  
 1st Qu.:1.720   1st Qu.: 66.00  
 Median :1.760   Median : 73.00  
 Mean   :1.764   Mean   : 74.47  
 3rd Qu.:1.820   3rd Qu.: 82.00  
 Max.   :1.930   Max.   :105.00  
n: 51
Correlacao: 0.646062

Pela localização e inclinação das elipses, parece que os homens são mais altos e mais pesados que as mulheres, e também são maiores suas variabilidades em estatura e em peso. Para avaliar visualmente as correlações, no entanto, devemos verificar com as variáveis padronizadas.

Para detectar outlier bivariado, John Tukey criou o bagplot. Implementamos em demo_Bagplot.R um bagplot que relata quem são os outliers; neste exemplo aparece um outlier bidimensional entre as mulheres:

Dados <- readxl::read_excel("Adm2008.xlsx")
Dados.F <- subset(Dados, Sexo=="Feminino")
Dados.M <- subset(Dados, Sexo=="Masculino")
n.F <- min(sum(!is.na(Dados.F$Estatura)), 
           sum(!is.na(Dados.F$MCT)))
n.M <- min(sum(!is.na(Dados.M$Estatura)), 
           sum(!is.na(Dados.M$MCT)))
bgp.F <- DescTools::PlotBag(Dados.F$Estatura, 
                            Dados.F$MCT, 
                            main=paste0("Feminino"," n=",n.F),
                            xlab="Estatura (m)",
                            ylab="Massa (kg)",
                            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.F <- as.data.frame(bgp.F$pxy.outlier))
for (o in 1:nrow(outliers.F))
{
  r.F <- which(Dados.F$Estatura==outliers.F$x[o] & 
                 Dados.F$MCT==outliers.F$y[o])
  text(outliers.F$x[o],outliers.F$y[o], r.F, pos=1, cex=0.7)
}
bgp.M <- DescTools::PlotBag(Dados.M$Estatura, 
                            Dados.M$MCT,
                            main=paste0("Masculino"," n=",n.M),
                            xlab="Estatura (m)",
                            ylab="Massa (kg)",
                            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.M <- as.data.frame(bgp.M$pxy.outlier))
if(nrow(outliers.M)>0){
for (o in 1:nrow(outliers.M)){
  r.M <- which(Dados.M$Estatura==outliers.M$x[o] & 
                 Dados.M$MCT==outliers.M$y[o])
  text(outliers.M$x[o], outliers.M$y[o], r.M, pos=1, cex=0.7)
} }
if(nrow(outliers.M)<=0){cat("Sem outlier bidimensional")}

     x  y
1 1.78 51

data frame com 0 coluna e 0 linha
Sem outlier bidimensional

Será que a correlação entre estatura e massa corpórea é diferente para os dois sexos?

Considerando que são membros da mesma espécie, não deveríamos observar a mesma correlação?

Implementamos em CorrelacaoPearsonDuasCorrelacoesIndepZ.R:

Dados <- readxl::read_excel("Adm2008.xlsx")
# padronizacao
Dados$Estatura.z <- unsplit(lapply(split(Dados$Estatura, Dados$Sexo), scale), 
                            Dados$Sexo)
Dados$MCT.z <- unsplit(lapply(split(Dados$MCT, Dados$Sexo), scale), 
                            Dados$Sexo)
nF <- sum(Dados$Sexo=="Feminino")
nM <- sum(Dados$Sexo=="Masculino")
# saida textual
source("eiras.bartitle.R")
 
# Elipse de predicao 95%
df <- data.frame(Dados$Estatura.z,Dados$MCT.z)

matriz <- as.matrix(df)
n <- nrow(matriz)
car::dataEllipse(matriz[,1], matriz[,2],
                 groups=factor(Dados$Sexo),
                 group.labels=c("F", "M"),
                 col=c("darkorange3","royalblue3"),
                 levels=0.95,
                 robust=TRUE,
                 main=paste("Elipse de predição de 95%\n",
                            "n = ",n," (Fem=",nF,", Masc=",nM,")",sep=""),
                 xlab="Estatura (z)",
                 ylab="MCT (z)",
                 xlim=c(-4,4),
                 ylim=c(-4,4),
                 lwd=1.5, lty=2)

Visualmente, a correlação parece igual para indivíduos dos dois sexos (feminino e masculino). Esta observação é acompanhada pelas elipses de predição 95% bastante coincidentes com as variáveis padronizadas sobrepostas. Podemos testar estatisticamente com a mesma função bootES::bootES utilizada para uma correlação, testando-se:

\[ \begin{cases} H_0:& \rho_F - \rho_M = 0\\ H_1:& \rho_F - \rho_M \ne 0 \end{cases}\\ \alpha = 0.05 \]

Este procedimento está implementado em CorrelacaoPearsonDuasCorrelacoesIndepTest.R:

Dados <- readxl::read_excel("Adm2008.xlsx")
B <- 1e4
rIC95 <- bootES::bootES(Dados[c("Estatura","MCT","Sexo")],
                        group.col="Sexo",
                        R=B)
print(rIC95) 

95.00% bca Confidence Interval, 10000 replicates
Stat        CI (Low)    CI (High)   bias        SE          
-0.009      -0.319      0.235       -0.004      0.139

Como o zero está contido no intervalo dado pelo teste, não rejeitamos \(H_0\). Concluímos que estatura e massa corpórea para mulheres e homens são, portanto, igualmente correlacionadas.

Existe alternativa para os casos em que os dados brutos não estão disponíveis, utilizando psych::r.test, como implementado em CorrelacaoPearsonDuasCorrelacoesIndepSemDadosBrutos.R. Supondo o exemplo que acabamos de ver caso somente soubéssemos os tamanhos dos grupos e suas correlações, obteríamos:

nF <- 38
corrF <- 0.637148
nM <- 51
corrM <- 0.646062 

out <- psych::r.test(n=nF, n2=nM, corrF, corrM)
cat("z = ",out$z,", p = ",out$p,"\n",sep="")
z = 0.06816731, p = 0.9456524

A conclusão sobre \(H_0\) é a mesma.

Comparação entre duas correlações dependentes

Suponha que 85 crianças do terceiro ano foram testadas com testes de inteligência (1), habilidades aritméticas (2) e compreensão de leitura (3). A correlação entre inteligência e habilidades aritméticas equivale a \(r_{1,2}\) = 0.53, inteligência e leitura se correlacionam com \(r_{1,3}\) = 0.41 e aritmética e leitura com \(r_{2,3}\) = 0.59. A correlação entre inteligência e habilidade aritmética é diferente da correlação entre inteligência e compreensão de leitura?

\[ \begin{cases} H_0:& \rho_{1,2} - \rho_{1,3} = 0\\ H_1:& \rho_{1,2} - \rho_{1,3} \ne 0 \end{cases}\\ \alpha = 0.05 \]

Esta situação pode ser verificada com a função psych::r.test, como implementada em CorrelacaoPearsonDuasCorrelacoesDepTest.R:

out <- psych::r.test(n=85, r12=0.53, r13=0.41, r23=0.59); cat("\n")
cat("t = ",out$t,", p = ",out$p,"\n",sep="") 

t = 1.418182, p = 0.159927

Pelo valor \(p\) verificamos que as correlações não são diferentes.

Correlação parcial

Retomando os dados de Adm2008.xlsx podemos verificar a correlação entre estatura e massa corpórea total removendo o efeito do sexo.

\[ \begin{cases} H_0:& {\rho_{(\text{estatura},\text{MCT}) \cdot \text{sexo}}} = 0\\ H_1:& {\rho_{(\text{estatura},\text{MCT}) \cdot \text{sexo}}} \ne 0 \end{cases}\\ \alpha = 0.05 \]

Está implementado em CorrelacaoPearsonParcialGenero.R:

Dados <- readxl::read_excel("Adm2008.xlsx")

# codifica mulher==1, homem==2
Dados$SexoNum <- NA
Sexos <- sort(unique(Dados$Sexo))
for (g in 1:length(Sexos)) 
{
  Dados$SexoNum[Dados$Sexo==Sexos[g]] <- g
}

cat("Correlação da Estatura com MCT: ")
rXY <- cor(Dados$Estatura,Dados$MCT)
cat(rXY," (r²=",rXY^2,")\n",sep="")
cat("Correlação da Estatura com Sexo: ")
rXZ <- cor(Dados$Estatura,Dados$SexoNum)
cat(rXZ," (r²=",rXZ^2,")\n",sep="")
cat("Correlação da MCT com Sexo: ")
rYZ <- cor(Dados$MCT,Dados$SexoNum)
cat(rYZ," (r²=",rYZ^2,")\n",sep="")

cat("Teste de correlação parcial\n",
    "(controlando pelo efeito linear de Sexo):\n")
rXY.Z <- ppcor::pcor.test(x=Dados$Estatura,
                         y=Dados$MCT,
                         z=Dados$SexoNum, 
                         method = "pearson")
print(rXY.Z)
cat("... rXY.Z²=",as.numeric(rXY.Z[1])^2,"\n",sep="")
Correlação da Estatura com MCT: 0.7902592 (r²=0.6245097)
Correlação da Estatura com Sexo: 0.628854 (r²=0.3954574)
Correlação da MCT com Sexo: 0.7037259 (r²=0.4952301)
Teste de correlação parcial
 (controlando pelo efeito linear de Sexo):
  estimate      p.value statistic  n gp  Method
1 0.629459 5.070604e-11  7.512368 89  1 pearson
... rXY.Z²=0.3962187

Já tinha sido verificado que as correlações entre Estatura e MCT eram altas para homens e mulheres. No entanto, a correlação de estatura com o sexo e de MCT com o sexo também são altas. Ao remover o efeito do sexo, sobre pouco entre estatura e MCT. Isto pode ser visto com um diagrama de Venn:

Quando os dados brutos não estiverem disponíveis, há uma solução em CorrelacaoPearsonParcialSemDadosBrutos.R. Simulamos a situação utilizando os valores aproximados das correlações computadas no exemplo com os dados brutos (respectivamente \(rXY=0.79\), \(rXZ=0.63\) e \(rYZ=0.70\)):

cat("Correlação da Estatura com MCT: ")
rXY <- 0.79 
cat(rXY," (r²=",rXY^2,")\n",sep="")
cat("Correlação da Estatura com Sexo: ")
rXZ <- 0.63
cat(rXZ," (r²=",rXZ^2,")\n",sep="")
cat("Correlação da MCT com Sexo: ")
rYZ <- 0.70
cat(rYZ," (r²=",rYZ^2,")\n",sep="")

cat("\nCorrelação Parcial de MCT e Estatura controlando por Sexo: ")
r.Parcial.Pearson <- (rXY - rXZ*rYZ)/(sqrt(1 - rXZ^2)*sqrt(1 - rYZ^2))
cat(r.Parcial.Pearson," (r²=",r.Parcial.Pearson^2,")\n",sep="") 
Correlação da Estatura com MCT: 0.79 (r²=0.6241)
Correlação da Estatura com Sexo: 0.63 (r²=0.3969)
Correlação da MCT com Sexo: 0.7 (r²=0.49)

Correlação Parcial de MCT e Estatura controlando por Sexo: 0.6292825 (r²=0.3959965)

Quando temos os dados brutos podemos fazer mais, por exemplo controlando para mais de uma variável. Por exemplo, podemos controlar simultaneamente para sexo e Idade.

\[ \begin{cases} H_0:& {\rho_{(\text{estatura},\text{MCT}) \cdot (\text{sexo},\text{idade})}} = 0\\ H_1:& {\rho_{(\text{estatura},\text{MCT}) \cdot (\text{sexo},\text{idade})}} \ne 0 \end{cases}\\ \alpha = 0.05 \]

Está implementado em CorrelacaoPearsonParcialGeneroIdade.R:

Dados <- readxl::read_excel("Adm2008.xlsx")
 
# codifica mulher==1, homem==2
Dados$SexoNum <- NA
Sexos <- sort(unique(Dados$Sexo))
for (g in 1:length(Sexos))
{
  Dados$SexoNum[Dados$Sexo==Sexos[g]] <- g
}

cat("Correlação da Estatura com MCT: ")
rXY <- cor(Dados$Estatura,Dados$MCT)
cat(rXY," (r²=",rXY^2,")\n",sep="")

cat("Correlação da Estatura com Sexo: ")
rXZ1 <- cor(Dados$Estatura,Dados$SexoNum)
cat(rXZ1," (r²=",rXZ1^2,")\n",sep="")
cat("Correlação da Estatura com Idade: ")
rXZ2 <- cor(Dados$Estatura,Dados$Idade)
cat(rXZ2," (r²=",rXZ2^2,")\n",sep="")
cat("Correlação da MCT com Sexo: ",sep="")
rYZ1 <- cor(Dados$MCT,Dados$SexoNum)
cat(rYZ1," (r²=",rYZ1^2,")\n",sep="")
cat("Correlação da MCT com Idade: ",sep="")
rYZ2 <- cor(Dados$MCT,Dados$Idade)
cat(rYZ2," (r²=",rYZ2^2,")\n",sep="")

cat("Correlação do Sexo com Idade: ")
rZ1Z2 <- cor(Dados$SexoNum,Dados$Idade)
cat(rZ1Z2," (r²=",rZ1Z2^2,")\n",sep="")

cat("Teste de Correlação parcial\n",
    "(controlando pelo efeito linear de Sexo e Idade):\n")
rXY.Z1Z2 <- ppcor::pcor.test(x=Dados$Estatura,
                          y=Dados$MCT,
                          z=c(Dados$SexoNum,Dados$Idade), 
                          method = "pearson")
print(rXY.Z1Z2)
cat("... rXY.Z1Z2²=",as.numeric(rXY.Z1Z2[1])^2,"\n",sep="")
Correlação da Estatura com MCT: 0.7902592 (r²=0.6245097)
Correlação da Estatura com Sexo: 0.628854 (r²=0.3954574)
Correlação da Estatura com Idade: 0.1065655 (r²=0.01135621)
Correlação da MCT com Sexo: 0.7037259 (r²=0.4952301)
Correlação da MCT com Idade: 0.3417262 (r²=0.1167768)
Correlação do Sexo com Idade: 0.2575778 (r²=0.06634632)
Teste de Correlação parcial
 (controlando pelo efeito linear de Sexo e Idade):
   estimate      p.value statistic   n gp  Method
1 0.7903242 4.441521e-39  17.06414 178  1 pearson
... rXY.Z1Z2²=0.6246123

É difícil entender o porquê da correlação parcial pode aumentar quando um novo controle é adicionado, pois uma área maior da intersecção entre Estatura e MCT tem que ser subtraída. A explicação está em que área fora da intersecção (originalmente pertencente a Estatura ou MCT) também foi retirada. Um diagrama de Venn aproximado é este:

O resultado final é que a área sombreada em cinza é cerca de 62% da área total remanescente.

Encontramos um vídeo sobre correlação parcial em Partial and semipartial correlation: YouTube.

FIM

Referências