Análise de Correlação em R
MSP4112: Estatística Aplicada à Medicina
14 outubro 2025 20:51h
Bastão de Asclépio & Distribuição Normal
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(psych, warn.conflicts=FALSE))
suppressMessages(library(psychTools, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(XICOR, 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")Código R
Correlacao_Original_e_Padronizada.RCorrelacao_r_de_Pearson.RCorrelacaoPearsonDuasCorrelacoesDepTest.RCorrelacaoPearsonDuasCorrelacoesIndep.RCorrelacaoPearsonDuasCorrelacoesIndepSemDadosBrutos.RCorrelacaoPearsonDuasCorrelacoesIndepTest.RCorrelacaoPearsonDuasCorrelacoesIndepZ.RCorrelacaoPearsonParcialGenero.RCorrelacaoPearsonParcialGeneroIdade.RCorrelacaoPearsonParcialSemDadosBrutos.R- CorrelacaoPearsonUmaCorrelacaoRho0.R
- `CorrelacaoPearsonUmaCorrelacaoRhoQualquer.R``
demo_Bagplot.Rdemo_covEstMCT.Rdemo_EstMCT.Rdemo_padronizaMCT.Rdemo_r_de_Spearman.Rdemo_r_de_Spearman_boot.Rdemo_r_de_Spearman_boot_monotonico.Rdemo_r_de_Spearman_boot_monotonicoD.Rdemo_r_de_Spearman_boot_quasemonotonico.RGestantes_corrBoot.RGestantes_descritiva.RGestantes_matrizcorrelacoes.RGestantes_rPearson.RGestantes_rPearson_boot_z.R
Material
- HTML de Rmarkdown 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, Spearman e \(\xi_n\).
- Testar uma correlação, uma correlação parcial, duas correlações independentes e duas correlações dependentes.
História de regressão & correlação
Galton (1988) introduz o termo co-relation e descreve, com base em medidas antropométricas, o conceito de associação entre variáveis intervalares, especialmente entre pais-filhos humanos.
Galton (1986) introduz o conceito de inclinação de regressão linear simples chamando-o de “ratio of regression” ou “rate of regression”. Ele observou empiricamente e erroneamente a “lei da regressão à mediocridade” (regression toward mediocrity), mas não tinha uma métrica precisa. Sua análise era essencialmente geométrica e descritiva.
Dados_height <- psychTools::galton # Mid parent child height Galton (1886)
sunflowerplot(data=Dados_height, child~parent)
lines(lowess(Dados_height$parent, Dados_height$child))
parent child
parent 1.0000000 0.4587624
child 0.4587624 1.0000000
Call:
lm(formula = child ~ parent, data = Dados_height)
Residuals:
Min 1Q Median 3Q Max
-7.8050 -1.3661 0.0487 1.6339 5.9264
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.94153 2.81088 8.517 <2e-16 ***
parent 0.64629 0.04114 15.711 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.239 on 926 degrees of freedom
Multiple R-squared: 0.2105, Adjusted R-squared: 0.2096
F-statistic: 246.8 on 1 and 926 DF, p-value: < 2.2e-16
Call:
lm(formula = child ~ parent, data = data.frame(scale(Dados_height)))
Residuals:
Min 1Q Median 3Q Max
-3.09976 -0.54256 0.01934 0.64889 2.35368
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.983e-15 2.918e-02 0.00 1
parent 4.588e-01 2.920e-02 15.71 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.889 on 926 degrees of freedom
Multiple R-squared: 0.2105, Adjusted R-squared: 0.2096
F-statistic: 246.8 on 1 and 926 DF, p-value: < 2.2e-16Dados_peas <- psychTools::peas # parent and child generation of 700 sweet peas Galton (1886)
sunflowerplot(data=Dados_peas, child~parent)
lines(lowess(Dados_peas$parent, Dados_peas$child))
parent child
parent 1.0000000 0.3463319
child 0.3463319 1.0000000
Call:
lm(formula = child ~ parent, data = Dados_peas)
Residuals:
Min 1Q Median 3Q Max
-2.6443 -1.5050 -0.3114 1.5129 5.6886
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.11429 0.63663 15.887 <2e-16 ***
parent 0.34286 0.03515 9.754 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.86 on 698 degrees of freedom
Multiple R-squared: 0.1199, Adjusted R-squared: 0.1187
F-statistic: 95.13 on 1 and 698 DF, p-value: < 2.2e-16
Call:
lm(formula = child ~ parent, data = data.frame(scale(Dados_peas)))
Residuals:
Min 1Q Median 3Q Max
-1.3346 -0.7596 -0.1572 0.7635 2.8711
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.795e-16 3.548e-02 0.000 1
parent 3.463e-01 3.551e-02 9.754 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9388 on 698 degrees of freedom
Multiple R-squared: 0.1199, Adjusted R-squared: 0.1187
F-statistic: 95.13 on 1 and 698 DF, p-value: < 2.2e-16Dados_cubit <- psychTools::heights # height and cubit Galton (1888)
sunflowerplot(data=Dados_cubit, height~cubit)
lines(lowess(Dados_cubit$cubit, Dados_cubit$height))
height cubit
height 1.0000000 0.7539688
cubit 0.7539688 1.0000000
Call:
lm(formula = height ~ cubit, data = Dados_cubit)
Residuals:
Min 1Q Median 3Q Max
-4.4259 -1.1521 -0.1521 1.3002 4.9849
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.9278 1.9291 13.44 <2e-16 ***
cubit 2.2739 0.1065 21.35 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.55 on 346 degrees of freedom
Multiple R-squared: 0.5685, Adjusted R-squared: 0.5672
F-statistic: 455.8 on 1 and 346 DF, p-value: < 2.2e-16
Call:
lm(formula = height ~ cubit, data = data.frame(scale(Dados_cubit)))
Residuals:
Min 1Q Median 3Q Max
-1.87860 -0.48899 -0.06454 0.55188 2.11585
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.267e-15 3.526e-02 0.00 1
cubit 7.540e-01 3.532e-02 21.35 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6579 on 346 degrees of freedom
Multiple R-squared: 0.5685, Adjusted R-squared: 0.5672
F-statistic: 455.8 on 1 and 346 DF, p-value: < 2.2e-16Pearson (1896) introduz formalmente o coeficiente de correlação, que mede a associação linear entre duas variáveis intervalares. Ele formalizou o coeficiente de correlação como a covariância padronizada, obtida pela soma dos desvios em relação à média aritmética. Ele resolveu um problema fundamental deixado por Galton: quantificar rigorosamente a relação linear entre duas variáveis biológicas associadas, como esetatura de pais e filhos, massa corporal total e comprimentos das partes do corpo humano etc. Galton
Pearson introduziu o coeficiente \(r\) para resolver isso: um número adimensional entre –1 e +1 que mede a força e direção da associação linear entre duas variáveis intervalares. Assim, ele:
- Formalizou a ideia de regressão linear como \(\hat{y} = a + b\,x\), onde \(b = r \dfrac{s_y}{s_x}\);
- Definiu \(r\) como a covariância normalizada, permitindo comparar relações entre variáveis de diferentes escalas;
- Eliminou o problema conceitual da assimetria galtoniana (em que a “regressão dos filhos nos pais” não era igual à “dos pais nos filhos”), mostrando que \(r_{xy} = r_{yx}\).
Pearson resolveu o problema de como medir, de forma objetiva e simétrica, o grau de dependência linear entre duas características biológicas. Este foi um passo decisivo na transição da biometria descritiva para a estatística inferencial.
Resumindo:
– Coeficiente de correlação é o coeficiente angular de regressão padronizado; – Coeficiente angular de regressão é a correlação escalada pela razão das variâncias.
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 (RLS) estima e testa a relação entre as variáveis intervalares dependente (VD) e independente (VI).
Conforme Bhagat (2018) e Mohammed (2018), as diferenças entre as análises de correlação e de regressão são apresentadas na tabela a seguir.
| Comparação por: | Correlação de Pearson | RLS |
|---|---|---|
|
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 |
|
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 |
|
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 |
|
Adimensional. Indica em qual medida as duas variáveis variam conjuntamente | Dimensional. Avalia o impacto da mudança de uma unidade de medida da variável independente na variável dependente |
|
Encontrar um valor numérico para expressar o grau de associação entre duas variáveis | Estimar o valor da variável dependente com base no valor da variável independente |
O coeficiente de correlação de Pearson é uma medida da associação linear 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 intervalar pode ser predito ou explicado pelo(s) valor(es) de outra(s) variável(is) intervalar(es) pela equação de regressão. Tem os seguintes tipos:
- Regressão linear simples (RLS): 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 (RLM): a variável dependente \(Y\) depende de duas ou mais variáveis independentes.
- Regressão não linear (RNL): se a relação entre a variável dependente \(Y\) e a(s) variável(is) independente(s) é não linear como acontece em relações exponenciais.
Exemplo
Apresentamos os dados utilizados numa dissertação de mestrado de Shinohara (1989) 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. O ácido fólico é a forma sintética e oxidada da vitamina B9 (folato). As deficiências de ácido fólico e B12 estão relacionadas com anemia megaloblástica e defeitos do tubo neural. Anemia megaloblástica (défict de folato/B12) é uma anemia macrocítica. No pré-natal, anemia se HB <11 g/dL. 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 = negra; 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 |
Há algumas variáveis quantitativas nestes dados, obtidos de 40 gestantes e 25 saudáveis. Usaremos apenas algumas destas variáveis do início do pré-natal:
- HT … Hematócrito (%)
- HB … Hemoglobina (g/dl)l
- HEM … Contagem de Hemácias (milhões/mm³)
- LEUC … Contagem de Leucócitos (milhares/mm³)
Há associação entre porcentagem do hematócrito e concentração de hemoglobina?
Há associação entre porcentagem do hematócrito e contagem de hemácias (eritrócitos)?
Há associação entre porcentagem do hematócrito e contagem de leucócitos?
A coleta da medida de porcentam de hematócrito é mais simples e fácil de ser obtida, pois seu o processamento requer uma centrífuga e uma régua.
A coleta de sangue venoso é mais trabalhosa, pois requer equipamento e um computador (Blood Testing at Warp Speed). Existem equipamentos mais sofisticados para automatizar o processo.
Panorama
Este diagrama pode ser útil para definir os testes adequados para cada situação. Para correlação e RLS, 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 <- data.frame(readxl::read_excel("Gestantes.xlsx"))
Dados$Grupo <- factor(Dados$Grupo)
Dados$COR <- factor(Dados$COR)
Dados$Primigesta <- factor(Dados$Primigesta)
Dados$TABAGISMO <- factor(Dados$TABAGISMO)
# Anemia no início do pré-natal (HB < 11 g/dL)
Dados$AnemiaIPN <- ifelse(is.na(Dados$HB), NA_integer_, as.integer(Dados$HB < 11))
Dados$AnemiaIPN <- factor(Dados$AnemiaIPN, levels=c(0,1), labels=c("Não","Sim"))
print(summary(Dados[,c(1,5:8,25)]), digits=1) Grupo HB HT HEM LEUC
controle:25 Min. :11.20 Min. :34.00 Min. :3.530 Min. : 3.600
Gestante:40 1st Qu.:12.50 1st Qu.:36.00 1st Qu.:4.240 1st Qu.: 6.400
Median :13.30 Median :38.00 Median :4.440 Median : 7.600
Mean :13.43 Mean :38.58 Mean :4.498 Mean : 7.902
3rd Qu.:14.20 3rd Qu.:41.00 3rd Qu.:4.650 3rd Qu.: 9.000
Max. :16.40 Max. :45.00 Max. :5.830 Max. :15.400
AnemiaIPN
Não:65
Sim: 0
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", "LEUC")
cat("Número de gestantes = ", nrow(Dados), "\n\n", sep="")
print(summary(Dados[Vars]), digits=2)
cat("\nMatriz de correlação:\n")
print(cor(Dados[Vars]), digits=2)
sunflowerplot(HB~HT,
pch=21,
xlab="Hematócrito (%)", ylab="Hemoglobina (mg/dl)",
data=Dados)
lines(lowess(Dados$HT, Dados$HB), lty=2)
sunflowerplot(HEM~HT,
pch=21,
xlab="Hematócrito (%)", ylab="Hemácias (milhão/mm³)",
data=Dados)
lines(lowess(Dados$HT, Dados$HEM), lty=2)
sunflowerplot(LEUC~HT,
pch=21,
xlab="Hematócrito (%)", ylab="Leucócitos (milhar/mm³)",
data=Dados)
lines(lowess(Dados$HT, Dados$LEUC), lty=2)Número de gestantes = 40
HT HB HEM LEUC
Min. :34.0 Min. :11.20 Min. :3.530 Min. : 4.20
1st Qu.:36.0 1st Qu.:12.40 1st Qu.:4.112 1st Qu.: 7.60
Median :37.0 Median :12.75 Median :4.280 Median : 8.40
Mean :37.4 Mean :13.04 Mean :4.318 Mean : 8.76
3rd Qu.:39.0 3rd Qu.:13.60 3rd Qu.:4.490 3rd Qu.: 9.70
Max. :45.0 Max. :16.40 Max. :5.380 Max. :15.40
Matriz de correlação:
HT HB HEM LEUC
HT 1.00 0.865 0.69 0.114
HB 0.86 1.000 0.59 0.035
HEM 0.69 0.587 1.00 0.168
LEUC 0.11 0.035 0.17 1.000
Observando os gráficos, cada ponto é um par de medidas de uma gestante. 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.
É 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 associadas. 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 programas estatísticos os três coeficientes de correlação: Coeficiente de correlação de Pearson
Karl Pearson (1857-1936) foi 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 (1863-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 (1907-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. |
Correlação de Pearson
O coeficiente de correlação 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
A correlação de Pearson avalia o grau de linearidade de duas variáveis intervalares sob a a suposição 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 da correlação de Pearson são 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) necessitam da independência entre os pares de observações (condição necessária), mas não necessitam da suposição de binormalidade (condição suficiente).
Observação: coeficiente de correlação de Pearson é adequado apenas para associações com formato de reta.
Princípio da correlação de Pearson
Admita duas variáveis intervalares como, e.g., estatura e massa
corpórea total (MCT), observadas (hipoteticamente) em nove 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)
Dados <- data.frame(estatura, massa)
print(Dados)
print(summary(Dados))
sunflowerplot(massa~estatura,
xlab = "Estatura (cm)",
ylab="MCT (kg)",
pch=21,
col="black",
data=Dados)
lines(lowess(Dados$estatura,
Dados$massa), lty=2)obtendo-se:
estatura massa
1 165 61
2 169 73
3 171 68
4 172 74
5 173 65
6 174 82
7 176 69
8 178 81
9 180 83
estatura massa
Min. :165.0 Min. :61.00
1st Qu.:171.0 1st Qu.:68.00
Median :173.0 Median :73.00
Mean :173.1 Mean :72.89
3rd Qu.:176.0 3rd Qu.:81.00
Max. :180.0 Max. :83.00 
Aparentemente, há uma relação crescente, aproximadamente linear, da MCT 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 nula de ausência de associação linear na população, i.e., \(\rho=0\) (do qual r é o estimador) para concluirmos se, populacionalmente, a associação linear existe:
\[ \begin{cases} H_0: \rho = 0 \\ H_1: \rho \ne 0 \end{cases}\\ \alpha=0.05 \]
Calculamos \(\hat{\rho}=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:
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
Dados <- data.frame(estatura, massa)
correlacao <- cor.test(~massa+estatura,
data=Dados)
print(correlacao)
Pearson's product-moment correlation
data: massa and estatura
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
Significado de correlação de Pearson
A correlação de Pearson é adimensional e adirecional, servindo para medir a intensidade da associação linear entre duas variáveis intervalares.
Apesar de utilizarmos cor.test com os valores brutos, a
correlação de Pearson é um coeficiente padronizado: é a intensidade da
associação das duas variáveis intervalares 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} = \dfrac{x - \bar{x}}{s_x} \]
e
\[ z_{y} = \dfrac{y - \bar{y}}{s_y} \]
Os respectivos escores \(z\) são adimensionais. Os formatos das respectivas distribuições não são afetadas.
|
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 imaginar que a padronização transforma a distribuição numa normal. 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 dos dados. A distribuição normal padronizada origina-se da 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 ( obtendo-se:
|
Podemos verificar que padronização não altera o valor da 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)
Dados <- data.frame(estatura, massa)
print(Dados)
print(summary(Dados))
cat("\nCorrelação com variáveis originais:\n")
correlacao <- cor.test(~massa+estatura,
data=Dados)
print(correlacao)
Dadoz <- data.frame(scale(Dados))
cat("\nVariáveis padronizadas:\n")
print(Dadoz, digits=1)
print(summary(Dadoz, digits=1))
cat("\nCorrelação com variáveis padronizadas:\n")
correlazao <- cor.test(~massa+estatura,
data=Dados)
print(correlazao) estatura massa
1 165 61
2 169 73
3 171 68
4 172 74
5 173 65
6 174 82
7 176 69
8 178 81
9 180 83
estatura massa
Min. :165.0 Min. :61.00
1st Qu.:171.0 1st Qu.:68.00
Median :173.0 Median :73.00
Mean :173.1 Mean :72.89
3rd Qu.:176.0 3rd Qu.:81.00
Max. :180.0 Max. :83.00
Correlação com variáveis originais:
Pearson's product-moment correlation
data: massa and estatura
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
Variáveis padronizadas:
estatura massa
1 -1.77 -1.51
2 -0.89 0.01
3 -0.46 -0.62
4 -0.24 0.14
5 -0.02 -1.00
6 0.19 1.16
7 0.63 -0.49
8 1.06 1.03
9 1.50 1.29
estatura massa
Min. :-1.77 Min. :-1.51
1st Qu.:-0.46 1st Qu.:-0.62
Median :-0.02 Median : 0.01
Mean : 0.00 Mean : 0.00
3rd Qu.: 0.63 3rd Qu.: 1.03
Max. : 1.50 Max. : 1.29
Correlação com variáveis padronizadas:
Pearson's product-moment correlation
data: massa and estatura
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) \] scale, que faz o mesmo de forma mais
simples.
|
Valores de correlação 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:
e escolha o jogo r in R. |
Correlação 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")
source("eiras.correg.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="Hematócrito (%)", 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 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
-----------
Pearson's product-moment correlation
data: Hematócrito (%) 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 da correlação de Pearson.
Observa-se, portanto, que a concentração de hemoglobina e a contagem de hemácias têm asssociação linear 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 de 95%, quando este não inclui o zero.
Correlação de Pearson com bootstrapping
Elaboramos Gestantes_rPearson_boot_z.R. Utilizamos
a mesma função (correg) adicionando o parâmetro
B reamostragens. As diferenças são: (i) a exibição de
várias elipses, cada uma vinda de uma reamostragem, e (ii) 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.6920725 0.9414919
sample estimates:
r (bootstrap)
0.866198 
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.3901000 0.8487659
sample estimates:
r (bootstrap)
0.6857914 
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.1318038 0.3956893
sample estimates:
r (bootstrap)
0.1195828 Na versão com bootstrapping, a função correg
mostra uma sombra das elipses variantes, evidenciando 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.
O gráfico a seguir mostra a evolução do intervalo de confiança da correlação de Pearson à medida que o tamanho da amostra aumenta.
Menor n com semi-amplitude <= 0,05 ao redor de ρ=0,5 (IC95%, Fisher): n = 809Tamanho de efeito
O tamanho do efeito é a magnitude da diferença entre condições ou o grau de relacionamento entre variáveis.
As propriedades da estimativa de tamanho de efeito são:
- Independe do tamanho da amostra (obrigatória)
- Em tabela de contigência, independe da dimensão (obrigatória)
- Adimensional
- Tem limites inferior (0) e superior (1)
Conforme Jones et al. (2025):
“Tamanhos de efeito padronizados são índices sem unidade usados para descrever a magnitude de uma associação.
Ao contrário dos valores de p, que são frequentemente usados para avaliar significância estatística, os tamanhos de efeito não dependem do tamanho da amostra (Betensky, 2019). Uma crítica bem conhecida aos valores de p e aos testes de significância é que, para tamanhos de amostra grandes, efeitos muito pequenos serão considerados significativos, mesmo que esses efeitos possam ser negligenciáveis na aplicação prática (Wasserstein & Lazar, 2016).
Em contraste, tamanhos de efeito comunicam a força do efeito, e não a mera existência de um efeito de tamanho arbitrário, o que pode ser mais significativo na prática (Sullivan & Feinn, 2012). Embora o aumento do tamanho da amostra ajude a melhorar a precisão da estimativa de um tamanho de efeito, o tamanho de efeito é um parâmetro que não depende do tamanho da amostra (Kang et al., 2023). Revistas e diretrizes estatísticas estão cada vez mais incentivando os autores a reportar tamanhos de efeito e seus intervalos de confiança juntamente com, ou no lugar de, valores de p (Wasserstein & Lazar, 2016; Wilkinson, 1999; Association, 1994; American Psychological Association, 2001, 2010; Althouse et al., 2021). No entanto, eles ainda não são comumente reportados (Fritz et al., 2012; Amaral & Line, 2021) e, quando reportados, frequentemente não incluem intervalos de confiança (Fritz et al., 2012).
Há quatro desafios para relatar tamanhos de efeito que limitam seu uso generalizado.
Primeiro, existem muitas medidas de tamanho de efeito disponíveis (Cohen, 1988; Hedges & Olkin, 1985; Rosenthal, 1994; Zhang & Schoeps, 1997; Serdar et al., 2021), mas elas são tipicamente definidas no contexto de um parâmetro populacional específico, o que dificulta a comparação de efeitos em uma ampla variedade de modelos (Vandekar et al., 2020).
Segundo, muitas medidas de tamanho de efeito disponíveis não permitem parâmetros de controle (nuisance parameters) ou covariáveis (Vandekar et al., 2020).
Terceiro, muitas medidas de tamanho de efeito não possuem procedimentos acurados para intervalos de confiança, o que impede a quantificação da incerteza em torno da estimativa do tamanho de efeito (Kang et al., 2023).
Finalmente, muitas funções padrão de resumo de modelos disponíveis em softwares estatísticos produzem automaticamente valores de p, mas poucas também reportam tamanhos de efeito.
Correlação de Pearson é uma medida natural de 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 entre VD e VE.
Ambos são medidas de tamanho de efeito, pois:
- São adimensionais;
- Variam entre 0 e 1;
- Não dependem do tamanho da amostra.
set.seed(123)
## Gerador não-normal: x ~ U(-1,1), y = 0.7*x + ruído uniforme
gen_sample <- function(n){
x <- runif(n, -1, 1)
y <- 0.7*x + runif(n, -0.5, 0.5) # relação linear + ruído NÃO-normal
list(x=x, y=y)
}
## Duplicar os dados não muda r
df <- data.frame(gen_sample(20))
sunflowerplot(scale(y)~scale(x), pch=21, data=df)
lines(lowess(scale(df$x),scale(df$y)))
r <- cor(df$x, df$y)
df2 <- rbind(df,df)
sunflowerplot(scale(y)~scale(x), pch=21, data=df2)
lines(lowess(scale(df2$x),scale(df2$y)))
r_dup <- cor(df2$x, df2$y) # duplica pares => mesma informação
cat(sprintf("r(n=20)=%.4f | r(n=40 duplicado)=%.4f\n", r, r_dup))r(n=20)=0.8846 | r(n=40 duplicado)=0.8846Classificação de tamanho de efeito de correlação de Pearson para planejamento
Ellis, 2010
Coeficiente de correlação de Spearman
Para distingui-lo da correlação 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".
\(s\) de Spearman também é aplicável se pelo menos uma das variáveis é considerada ordinal. Neste caso, a representação gráfica por meio do gráfico de dispersão não é uma representação válida.
Da mesma forma que a correlação de Pearson, \(s\) varia de \(-1\) a \(+1\).
A correlação de Spearman é a correlação de Pearson dos postos das duas variáveis. Desta forma, a padronização não é necessária.
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
Dados <- data.frame(estatura, massa)
# Spearman "oficial"
rho_s <- cor(Dados$estatura, Dados$massa,
method = "spearman")
# Pearson dos postos
rx <- rank(Dados$estatura)
ry <- rank(Dados$massa)
rho_p_ranks <- cor(rx, ry, method = "pearson")
# Opcional: mostrar que padronizar é redundante
rho_p_ranks_z <- as.numeric(cor(scale(rx), scale(ry),
method = "pearson"))
print(list(
spearman = rho_s,
pearson_dos_postos = rho_p_ranks,
pearson_dos_postos_padronizados = rho_p_ranks_z
))$spearman
[1] 0.7
$pearson_dos_postos
[1] 0.7
$pearson_dos_postos_padronizados
[1] 0.7A 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} \]
alfa <- 0.05
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
Dados <- data.frame(estatura, massa)
sunflowerplot(massa~estatura,
xlab = "Estatura (cm)",
ylab="MCT (kg)",
pch=21,
col="black",
data=Dados)
lines(lowess(Dados$estatura,
Dados$massa), lty=2)
Pearson's product-moment correlation
data: massa and estatura
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
Spearman's rank correlation rho
data: massa and estatura
S = 36, p-value = 0.04325
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.7 rho lwr.ci upr.ci
0.70000 0.06705 0.93121 O coeficiente de correlação de Spearman é uma medida natural de tamanho de efeito.
set.seed(123)
## Gerador não-normal: x ~ U(-1,1), y = 0.7*x + ruído uniforme
gen_sample <- function(n){
x <- runif(n, -1, 1)
y <- 0.7*x + runif(n, -0.5, 0.5) # relação linear + ruído NÃO-normal
list(x=x, y=y)
}
## Duplicar os dados não muda r
df <- data.frame(gen_sample(20))
sunflowerplot(scale(y)~scale(x), pch=21, data=df)
lines(lowess(scale(df$x),scale(df$y)))
r <- cor(df$x, df$y, method="spearman")
df2 <- rbind(df,df)
sunflowerplot(scale(y)~scale(x), pch=21, data=df2)
lines(lowess(scale(df2$x),scale(df2$y)))
r_dup <- cor(df2$x, df2$y, method="spearman") # duplica pares => mesma informação
cat(sprintf("s(n=20)=%.4f | s(n=40 duplicado)=%.4f\n", r, r_dup))s(n=20)=0.8707 | s(n=40 duplicado)=0.8707O 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
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
Dados <- data.frame(estatura, massa)
col <- friendlycolor(31) # preto
bg <- friendlycolor(24) # amarelo
pch <- 21 # circulo
lst <- correg(Dados$estatura, Dados$massa,
method="spearman",
main="",
xlab="Estatura (cm)",
ylab="MCT (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)
}
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 MCT (kg) 74 7.8652 61 73 68 65 82 9 0
lambda = 0.855927523022535
Assuming:
- Explanatory variable: Estatura (cm);
- Dependent variable: MCT (kg).
-----------
Correlation
-----------
Spearman's rank correlation rho
data: Estatura (cm) and MCT (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 motivo de \(|s| < 1\): à medida que a estatura (\(x\)) aumenta, os valores de massa corporal (\(y\)) ora crescem, ora decrescem.
|
Implementamos, também, uma versão em bootstrapping em
O teorema central do limite está razoavelmente atendido. |
O coeficiente de correlação de Spearman é máximo se o crescimento é
monotônico (\(s=1\)) ou se o
decrescimento é monotônico (\(s=-1\)).
E.g., executando: demo_r_de_Spearman_boot_monotonico.R
obtemos:
# demo_r_de_Spearman_boot_monotonico.R
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
# 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,2),"\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")x: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
y: 1.97 2.32 10.18 15.46 21.02 35.49 39.05 39.21 50.33 69.53 144.08 266.42 351.39 496.11 693.31 894.6 1178.87 1208.26 1422.7 1441.13 
Variable Mean Sd Min. 1st Qu. Median 3rd Qu.
1 x 4 5.9161 1 2 3 5
2 y 15.4564 522.7931 1.9718 2.3159 10.182 21.0211078875745
Max. n NA's
1 6 20 0
2 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, modificamos um valor do exemplo anterior com:
Obtendo:
# demo_r_de_Spearman_boot_quasemonotonico.R
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
# 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 reduzidox: 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 
Variable Mean Sd Min. 1st Qu. Median 3rd Qu.
1 x 4 5.9161 1 2 3 5
2 y 15.4564 523.758 1.9718 2.3159 10.182 21.0211078875745
Max. n NA's
1 6 20 0
2 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, |
O valor é \(r=-1\) para o decrescimento monotônico:
# demo_r_de_Spearman_boot_monotonicoD.R
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
# 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")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 
Variable Mean Sd Min. 1st Qu. Median 3rd Qu.
1 x 4 5.9161 1 2 3 5
2 y 1681.6555 391.1633 1697.175 1693.5109 1692.3845 1662.21456341912
Max. n NA's
1 6 20 0
2 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 As correlações de Pearson e Spearman são nulas no caso de pontos em formato perfeitamente parabólico e simétrico em torna da origem.
x <- -10:10
y <- x^2
cat("x:",x,"\n")
cat("y:",round(y,3),"\n")
plot(x,y)
print(cor.test(x,y))
print(cor.test(x,y, method="spearman"))
print(DescTools::SpearmanRho(x,y, conf.level=0.95), digits=4)
set.seed(123)
xin <- XICOR::xicor(x,y,pvalue=TRUE)
cat("\nksi_n = ", round(xin$xi,2), " p = ", signif(xin$pval,2), sep="")x: -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10
y: 100 81 64 49 36 25 16 9 4 1 0 1 4 9 16 25 36 49 64 81 100 
Pearson's product-moment correlation
data: x and y
t = 0, df = 19, p-value = 1
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4316868 0.4316868
sample estimates:
cor
0
Spearman's rank correlation rho
data: x and y
S = 1540, p-value = 1
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0
rho lwr.ci upr.ci
0.0000 -0.4317 0.4317
ksi_n = 0.73 p = 1e-07Matriz de correlação de Pearson no exemplo de Gestante
Exemplificamos em Gestantes_matrizcorrelacoes.R a
montagem de uma matriz de correlações de Pearson e de Spearman. Usamos a
função GGally::ggcorr. 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
Gestantes <- readRDS("Gestante.rds")
Dados <- data.frame(Gestantes$IDADE,
Gestantes$HT,
Gestantes$HB,
Gestantes$HEM,
Gestantes$LEUC,
Gestantes$FOLICO,
Gestantes$B12)
names(Dados) <- c("IDADE","HT","HB","HEM","LEUC","FOLICO","B12")
cat("\nMatriz de correlações:\n")
print(cor(Dados), digits=1) # matriz de correlacoes
Dadoz <- scale(Dados)
car::scatterplotMatrix(~IDADE+HT+HB+HEM+LEUC+FOLICO+B12,
col="black",
smooth=FALSE,
ellipse=list(levels=c(.68),
robust=TRUE,
fill=FALSE),
regLine=FALSE,
cex=0.8,
plot.points=TRUE,
data=Dadoz)
print(GGally::ggcorr(Dadoz,
name="Pearson",
nbreaks = 6,
digits = 3,
label = TRUE,
label_size = 4,
color = "#888888"))
print(cor(Dados, method="spearman"), digits=1) # matriz de correlacoes
print(GGally::ggcorr(Dadoz,
name="Spearman",
method=c("complete.obs", "spearman"),
nbreaks = 6,
digits = 3,
label = TRUE,
label_size = 4,
color = "#888888"))
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
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
IDADE HT HB HEM LEUC FOLICO B12
IDADE 1.00 -0.02 -0.06 0.06 0.07 0.05 -0.05
HT -0.02 1.00 0.78 0.61 0.19 0.27 -0.03
HB -0.06 0.78 1.00 0.52 0.11 -0.02 -0.16
HEM 0.06 0.61 0.52 1.00 0.19 0.24 0.36
LEUC 0.07 0.19 0.11 0.19 1.00 0.03 -0.21
FOLICO 0.05 0.27 -0.02 0.24 0.03 1.00 0.32
B12 -0.05 -0.03 -0.16 0.36 -0.21 0.32 1.00
Variância compartilhada
O coeficiente de correlação 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, correg, 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, 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 entre VD e VE.”
Qual é o significado disto?
É fácil ver por meio de um diagrama de Venn. Como as variáveis são compartilhadas (i.e., todas têm média nula e variância unitária), todas são representadas por círculos de áreas iguais e unitárias. 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\)). Num exemplo hipotético de hematócrito contra hemoglobina, hemácias e leucócitos, os valores de r e \(r^2\) são:
Níveis de mensuração e outros tipos de correlação
Existem alguns outros 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 testa a hipótese nula de
que não há associação linear populacional entre duas variáveis
intervalares (\(H_0: \rho = 0\)).
Em publicações científicas, 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. E.g., 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 \]
O resultado é:
alfa <- 0.05
n <- 200
r <- 0.5
# H₀: ρ = 0 vs. H₁: ρ ≠ 0
# https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
cat("Teste bilateral:\n")
cat("\tH₀: ρ = 0 vs. H₁: ρ ≠ 0\n")
t <- r*sqrt((n-2)/(1-r^2))
p <- 2*pt(-abs(t), df=n-2)
cat("t(",n-2,") = ",round(t,2)," , p = ",round(p,4),"\n",sep="")
r <- DescTools::CorCI(r, n=n, conf.level = 1-alfa,
alternative = "two.sided")
print(r, digits=2)
cat("\n")
cat("Teste unilateral à direita:\n")
cat("\tH₀: ρ = 0 vs. H₁: ρ > 0\n")
p <- 1-pt(t, df=n-2)
cat("t(",n-2,") = ",round(t,2)," , p = ",round(p,4),"\n",sep="")
r <- DescTools::CorCI(r, n=n, conf.level = 1-alfa,
alternative = "greater")
print(r, digits=2)
cat("\n")
cat("Teste unilateral à esquerda:\n")
cat("\tH₀: ρ = 0 vs. H₁: ρ < 0\n")
p <- pt(t, df=n-2)
cat("t(",n-2,") = ",round(t,2)," , p = ",round(p,4),"\n",sep="")
r <- DescTools::CorCI(r, n=n, conf.level = 1-alfa,
alternative = "less")
print(r, digits=2)Teste bilateral:
H₀: ρ = 0 vs. H₁: ρ ≠ 0
t(198) = 8.12 , p = 0
cor lwr.ci upr.ci
0.50 0.39 0.60
Teste unilateral à direita:
H₀: ρ = 0 vs. H₁: ρ > 0
t(198) = 8.12 , p = 0
cor.cor cor.lwr.ci cor.upr.ci lwr.ci.cor upr.ci.lwr.ci
0.50 0.39 0.60 0.41 0.28
Teste unilateral à esquerda:
H₀: ρ = 0 vs. H₁: ρ < 0
t(198) = 8.12 , p = 1
cor.cor.cor cor.cor.lwr.ci cor.cor.upr.ci cor.lwr.ci.cor
0.50 0.39 0.60 0.41
cor.upr.ci.lwr.ci lwr.ci upr.ci.cor.cor
0.28 -1.00 0.58 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 é correlação positiva).
A função DescTools::CorCI produz apenas o intervalo de
confiança da correlação de Pearson. O valor p é obtido por
fórmula (CorrelacaoPearsonUmaCorrelacaoRho0.R):
alfa <- 0.05
n <- 200
r <- 0.5
cat("\n")
cat("Teste de correlação:\n")
cat("\tH₀: ρ = 0 vs. H₁: ρ ≠ 0\n",sep="")
test <- DescTools::CorCI(rho=r,
n=n,
conf.level=1-alfa,
alternative="two.sided")
print(test, digits=3)
z <- sqrt((n-3)/2)*log(1/((1-r)/(1+r)))
p <- 2*pnorm(-abs(z))
cat("z = ",round(z,2),", p = ",round(p,4), 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,") = ",round(t,2)," , p = ",round(p,4),"\n",sep="")
Teste de correlação:
H₀: ρ = 0 vs. H₁: ρ ≠ 0
cor lwr.ci upr.ci
0.500 0.388 0.597
z = 10.9, p = 0
Alternativamente:
t(198) = 8.12 , p = 0No 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("\tH₀: ρ = ",ro," vs. H₁: ρ ≠ ",ro,"\n",sep="")
z <- sqrt((n-3)/2)*log(((1-ro)/(1+ro))/((1-r)/(1+r)))
p <- 2*pnorm(-abs(z))
cat("z = ",round(z,2),", p = ",round(p,4), sep="")
Teste de correlação:
H₀: ρ = 0.6 vs. H₁: ρ ≠ 0.6
z = -2.86, p = 0.0043Quando temos os dados brutos, há diversas funções em R que podem fornecer o intervalo de confiança 95%. E.g., no caso da correlação de Pearson 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("\tH₀: ρ = ",ro," vs. H₁: ρ ≠ ",ro,"\n",sep="")
z <- sqrt((n-3)/2)*log(((1-ro)/(1+ro))/((1-r)/(1+r)))
p <- 2*pnorm(-abs(z))
cat("z = ",round(z,2),", p = ",round(p,4), sep="")Hematócrito e Hemoglobina:
95.00% bca Confidence Interval, 10000 replicates
Stat CI (Low) CI (High) bias SE
0.865 0.698 0.941 -0.009 0.063
Hematócrito e Hemácias:
95.00% bca Confidence Interval, 10000 replicates
Stat CI (Low) CI (High) bias SE
0.686 0.424 0.860 -0.019 0.117
Hematócrito e Leucócitos:
95.00% bca Confidence Interval, 10000 replicates
Stat CI (Low) CI (High) bias SE
0.114 -0.170 0.369 0.010 0.136 Comparação entre duas correlações de Pearson independentes
Exemplo 1: Animal humano
Padronizar as variáveis e proceder à regressão linear pode ajudar na
exploração das relações entre VD e VE. 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.
Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosF <- subset(Dados, Sexo=="Feminino")
DadosM <- subset(Dados, Sexo=="Masculino")
n <- nrow(Dados)
nF <- nrow(DadosF)
nM <- nrow(DadosM)
cat("MCT:\n")
print(FSA::Summarize(MCT~Sexo,
digits=2,
data=Dados))
cat("Estatura:\n")
print(FSA::Summarize(Estatura~Sexo,
digits=2,
data=Dados))
## ----- IC a partir de cor.test (Pearson) e gráfico lado a lado -----
# Testes por grupo (cor.test já devolve IC95% em $conf.int)
alfa <- 0.05
g <- nlevels(Dados$Sexo)
alfaBonf <- alfa/g
ctF <- cor.test(~DadosF$Estatura+DadosF$MCT,
method = "pearson",
conf.level = 1 - alfaBonf,
use = "complete.obs",
data=Dados)
ctM <- cor.test(~DadosM$Estatura+DadosM$MCT,
method = "pearson",
conf.level = 1 - alfaBonf,
use = "complete.obs",
data=Dados)
rF <- unname(ctF$estimate)[1]
ciF <- unname(ctF$conf.int) # vetor c(lo, hi)
nF_eff <- sum(complete.cases(DadosF$Estatura, DadosF$MCT))
rM <- unname(ctM$estimate)[1]
ciM <- unname(ctM$conf.int)
nM_eff <- sum(complete.cases(DadosM$Estatura, DadosM$MCT))
TAB.ic <- data.frame(
Sexo = c("Feminino","Masculino"),
n = c(nF_eff, nM_eff),
r = c(rF, rM),
lo = c(ciF[1], ciM[1]),
hi = c(ciF[2], ciM[2])
)
cat("\nIC95% Bonferroni de correlação de Pearson\n")
print(TAB.ic, digits = 3, row.names = FALSE)
## Gráfico dos ICs lado a lado (base R)
pal <- c("blue","black")
ypos <- 1:2
plot(NA, xlim = c(0,1), ylim = c(0.5,2.5),
xlab = "r",
ylab = "", yaxt = "n",
main = "IC95% Bonferroni de correlação")
abline(v = 0, lty = 2, col = "grey70")
segments(TAB.ic$lo, ypos, TAB.ic$hi, ypos, lwd = 3, col = pal)
points(TAB.ic$r, ypos, pch = 19, cex = 1.2, col = pal)
axis(2, at = ypos, labels = paste0(" n=", TAB.ic$n), las = 2)
text(TAB.ic$hi, ypos + 0.08, labels = sprintf("r=%.3f", TAB.ic$r), cex = 0.9, pos = 2)
# Elipse de predicao 95%
matriz <- as.matrix(Dados[, 3:4])
car::dataEllipse(matriz[,1], matriz[,2],
groups=Dados$Sexo,
group.labels=c("F", "M"),
col=c("blue","black"),
grid=FALSE,
levels=c(.68),
robust=TRUE,
main=paste("Elipse de predição de 68%\n",
"n = ",n," (Fem=",nF,", Masc=",nM,")",
sep=""),
xlab="Estatura (m)",
ylab="MCT (kg)",
xlim=c(1.4, 2.0),
ylim=c(30, 110),
lwd=2, lty=1)MCT:Registered S3 methods overwritten by 'FSA':
method from
confint.boot car
hist.boot car Sexo n mean sd min Q1 median Q3 max
1 Feminino 38 54.63 6.19 43 50.25 53.5 59.75 70
2 Masculino 51 74.47 12.10 48 66.00 73.0 82.00 105
Estatura:
Sexo n mean sd min Q1 median Q3 max
1 Feminino 38 1.64 0.07 1.50 1.60 1.64 1.70 1.78
2 Masculino 51 1.76 0.08 1.56 1.72 1.76 1.82 1.93
IC95% Bonferroni de correlação de Pearson
Sexo n r lo hi
Feminino 38 0.637 0.358 0.812
Masculino 51 0.646 0.418 0.798
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
|
A correlação de Peaarson 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 <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosF <- subset(Dados, Sexo=="Feminino")
DadosM <- subset(Dados, Sexo=="Masculino")
n <- nrow(Dados)
nF <- nrow(DadosF)
nM <- nrow(DadosM)
# 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)
df <- data.frame(Dados$Estatura.z,Dados$MCT.z)
matriz <- as.matrix(df)
car::dataEllipse(matriz[,1], matriz[,2],
grid=FALSE,
groups=factor(Dados$Sexo),
group.labels=c("F", "M"),
col=c("darkorange3","royalblue3"),
levels=0.68,
robust=TRUE,
main=paste("Elipse de predição de 68%\n",
"n = ",n," (Fem=",nF,", Masc=",nM,")",sep=""),
xlab="Estatura (z)",
ylab="MCT (z)",
xlim=c(-3,3),
ylim=c(-3,3),
lwd=2, lty=1)
Visualmente, a correlação de Pearson parece igual para indivíduos dos
dois sexos (feminino e masculino). Esta observação é acompanhada pelas
elipses de predição 68% com as variáveis padronizadas que estão bem
sobrepostas. Podemos testar estatisticamente com a mesma função
bootES::bootES utilizada para uma correlação:
\[ \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 <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
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.304 0.233 -0.002 0.137 Como o zero está contido no intervalo dado pelo teste, não rejeitamos \(H_0\).
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, obtemos:
nF <- 38
corrF <- 0.637148
nM <- 51
corrM <- 0.646062
out <- psych::r.test(n=nF, n2=nM, corrF, corrM)
cat("z = ",round(out$z,2),", p = ",round(out$p,4),"\n",sep="")z = 0.07, p = 0.9457A conclusão sobre \(H_0\) é a mesma: não rejeitar.
Exemplo 2: Flor iris
A seguir é apresentado o exemplo sobre a correlação entre largura e
comprimento de sépala para três espécies de iris.
Sepal.Length Sepal.Width Species
Min. :4.300 Min. :2.000 setosa :50
1st Qu.:5.100 1st Qu.:2.800 versicolor:50
Median :5.800 Median :3.000 virginica :50
Mean :5.843 Mean :3.057
3rd Qu.:6.400 3rd Qu.:3.300
Max. :7.900 Max. :4.400
item group1 vars n mean sd median trimmed mad min max
Sepal.Length1 1 setosa 1 50 5.01 0.35 5.0 5.00 0.30 4.3 5.8
Sepal.Length2 2 versicolor 1 50 5.94 0.52 5.9 5.94 0.52 4.9 7.0
Sepal.Length3 3 virginica 1 50 6.59 0.64 6.5 6.57 0.59 4.9 7.9
Sepal.Width1 4 setosa 2 50 3.43 0.38 3.4 3.42 0.37 2.3 4.4
Sepal.Width2 5 versicolor 2 50 2.77 0.31 2.8 2.78 0.30 2.0 3.4
Sepal.Width3 6 virginica 2 50 2.97 0.32 3.0 2.96 0.30 2.2 3.8
range skew kurtosis se
Sepal.Length1 1.5 0.11 -0.45 0.05
Sepal.Length2 2.1 0.10 -0.69 0.07
Sepal.Length3 3.0 0.11 -0.20 0.09
Sepal.Width1 2.1 0.04 0.60 0.05
Sepal.Width2 1.4 -0.34 -0.55 0.04
Sepal.Width3 1.6 0.34 0.38 0.05
IC95% Bonferroni de correlação de Pearson
Especie n r lo hi
setosa 50 0.743 0.542 0.863
versicolor 50 0.526 0.231 0.732
virginica 50 0.457 0.144 0.687
Call:
lm(formula = Sepal.Length ~ scale(Sepal.Width, scale = FALSE),
data = I.S)
Residuals:
Min 1Q Median 3Q Max
-0.52476 -0.16286 0.02166 0.13833 0.44428
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.00600 0.03373 148.392 < 2e-16 ***
scale(Sepal.Width, scale = FALSE) 0.69049 0.08990 7.681 6.71e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2385 on 48 degrees of freedom
Multiple R-squared: 0.5514, Adjusted R-squared: 0.542
F-statistic: 58.99 on 1 and 48 DF, p-value: 6.71e-10
Call:
lm(formula = Sepal.Length ~ scale(Sepal.Width, scale = FALSE),
data = I.Vc)
Residuals:
Min 1Q Median 3Q Max
-0.73497 -0.28556 -0.07544 0.43666 0.83805
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.93600 0.06273 94.627 < 2e-16 ***
scale(Sepal.Width, scale = FALSE) 0.86508 0.20194 4.284 8.77e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4436 on 48 degrees of freedom
Multiple R-squared: 0.2766, Adjusted R-squared: 0.2615
F-statistic: 18.35 on 1 and 48 DF, p-value: 8.772e-05
Call:
lm(formula = Sepal.Length ~ scale(Sepal.Width, scale = FALSE),
data = I.Vg)
Residuals:
Min 1Q Median 3Q Max
-1.26067 -0.36921 -0.03606 0.19841 1.44917
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.58800 0.08081 81.529 < 2e-16 ***
scale(Sepal.Width, scale = FALSE) 0.90153 0.25311 3.562 0.000843 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5714 on 48 degrees of freedom
Multiple R-squared: 0.2091, Adjusted R-squared: 0.1926
F-statistic: 12.69 on 1 and 48 DF, p-value: 0.0008435
# Teste de paralelismo (igualdade de inclinações entre 3 grupos)
Analysis of Variance Table
Model 1: Sepal.Length ~ Sepal.Width + Species
Model 2: Sepal.Length ~ Sepal.Width * Species
Res.Df RSS Df Sum of Sq F Pr(>F)
1 146 28.0
2 144 27.9 2 0.157 0.41 0.67
# Interceptos iguais supondo paralelismo (3 grupos)
Analysis of Variance Table
Model 1: Sepal.Length ~ Sepal.Width
Model 2: Sepal.Length ~ Sepal.Width + Species
Res.Df RSS Df Sum of Sq F Pr(>F)
1 148 101
2 146 28 2 72.8 190 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Correlações por grupo e testes pareados (independentes)
cor(setosa) = 0.743
cor(versicolor)= 0.526
cor(virginica) = 0.457
Diferenças de correlação (teste de duas correlações independentes):
Correlation tests
Call:psych::r.test(n = nS, r12 = rS, r34 = rVc, n2 = nVc)
Test of difference between two independent correlations
z value 1.8017 with probability 0.0716
Correlation tests
Call:psych::r.test(n = nS, r12 = rS, r34 = rVg, n2 = nVg)
Test of difference between two independent correlations
z value 2.2412 with probability 0.025
Correlation tests
Call:psych::r.test(n = nVc, r12 = rVc, r34 = rVg, n2 = nVg)
Test of difference between two independent correlations
z value 0.4396 with probability 0.6603Comparação entre duas correlações de Pearson 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 = ",round(out$t,2),", p = ",round(out$p,4),"\n",sep="")
t = 1.42, p = 0.1599Conclusão:
A correlação entre inteligência e habilidades aritméticas (0.53) não é estatisticamente diferente da correlação entre inteligência e compreensão de leitura (0.41).
Em termos práticos: ambas as relações são semelhantes em magnitude; não há diferença significante entre o grau de associação de “inteligência” com “aritmética” e com “leitura”.
Correlação de Pearson parcial
Retomando os dados de Adm2008.xlsx, podemos verificar a
correlação de Pearson 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 <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
# codifica Feminino==1, Masculino==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(round(rXY,2)," (r² = ",round(rXY^2,2),")\n",sep="")
cat("Correlação da Estatura com Sexo: ")
rXZ <- cor(Dados$Estatura,Dados$SexoNum)
cat(round(rXZ,2)," (r² = ",round(rXZ^2,2),")\n",sep="")
cat("Correlação da MCT com Sexo: ")
rYZ <- cor(Dados$MCT,Dados$SexoNum)
cat(round(rYZ,2)," (r² = ",round(rYZ^2,2),")\n",sep="")
cat("\nTeste 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")
cat("\nCorrelação Parcial de MCT e Estatura controlando por Sexo:\n")
print(rXY.Z, digits=2)
cat("\n(rXY.Z)² = ",round(rXY.Z$estimate,2), "^2 = ", round(as.numeric(rXY.Z$estimate)^2,2),"\n",sep="")Correlação da Estatura com MCT: 0.79 (r² = 0.62)
Correlação da Estatura com Sexo: 0.63 (r² = 0.4)
Correlação da MCT com Sexo: 0.7 (r² = 0.5)
Teste de correlação parcial
(controlando pelo efeito linear de Sexo):
Correlação Parcial de MCT e Estatura controlando por Sexo:
estimate p.value statistic n gp Method
1 0.63 5.1e-11 7.5 89 1 pearson
(rXY.Z)² = 0.63^2 = 0.4Já 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, diminui a correlação entre estatura e MCT, mas continua significante. 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(round(rXY,2)," (r²=",round(rXY^2,2),")\n",sep="")
cat("Correlação da Estatura com Sexo: ")
rXZ <- 0.63
cat(round(rXZ,2)," (r²=",round(rXZ^2,2),")\n",sep="")
cat("Correlação da MCT com Sexo: ")
rYZ <- 0.70
cat(round(rYZ,2)," (r²=",round(rYZ^2,2),")\n",sep="")
cat("\nCorrelação Parcial de MCT e Estatura controlando por Sexo:\n")
r.Parcial.Pearson <- (rXY - rXZ*rYZ)/(sqrt(1 - rXZ^2)*sqrt(1 - rYZ^2))
cat("rXY.Z = ",round(r.Parcial.Pearson,2)," ((rXY.Z)² = ",round(r.Parcial.Pearson^2,2),")\n",sep="") Correlação da Estatura com MCT: 0.79 (r²=0.62)
Correlação da Estatura com Sexo: 0.63 (r²=0.4)
Correlação da MCT com Sexo: 0.7 (r²=0.49)
Correlação Parcial de MCT e Estatura controlando por Sexo:
rXY.Z = 0.63 ((rXY.Z)² = 0.4)Note que a correlação parcial de estatura e MCT controlando por sexo de 0.63 é aproximadamente igual às correlações para as condições de sexo masculino e feminino de 0.65 e 0.64, respectivamente.
Quando temos os dados brutos podemos fazer mais. Por exemplo, controlando para mais de uma variável. 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 <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
# codifica mulher==1, homem==2
Dados$SexoNum <- NA
Sexos <- factor(levels(Dados$Sexo))
for (g in 1:nlevels(Dados$Sexo))
{
Dados$SexoNum[Dados$Sexo==Sexos[g]] <- g
}
cat("Correlação da Estatura com MCT: ")
rXY <- cor(Dados$Estatura,Dados$MCT)
cat(round(rXY,2)," (r² = ",round(rXY^2,2),")\n",sep="")
cat("Correlação da Estatura com Sexo: ")
rXZ1 <- cor(Dados$Estatura,Dados$SexoNum)
cat(round(rXZ1,2)," (r² = ",round(rXZ1^2,2),")\n",sep="")
cat("Correlação da Estatura com Idade: ")
rXZ2 <- cor(Dados$Estatura,Dados$Idade)
cat(round(rXZ2,2)," (r² = ",round(rXZ2^2,2),")\n",sep="")
cat("Correlação da MCT com Sexo: ",sep="")
rYZ1 <- cor(Dados$MCT,Dados$SexoNum)
cat(round(rYZ1,2)," (r² = ",round(rYZ1^2,2),")\n",sep="")
cat("Correlação da MCT com Idade: ",sep="")
rYZ2 <- cor(Dados$MCT,Dados$Idade)
cat(round(rYZ2,2)," (r² = ",round(rYZ2^2,2),")\n",sep="")
cat("Correlação do Sexo com Idade: ")
rZ1Z2 <- cor(Dados$SexoNum,Dados$Idade)
cat(round(rZ1Z2,2)," (r²=",round(rZ1Z2^2,2),")\n",sep="")
cat("\nTeste 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, digits=2)
cat("\n(rXY.Z1Z2)² = ",round(rXY.Z1Z2$estimate,2), "^2 = ", round(as.numeric(rXY.Z1Z2$estimate)^2,2),"\n",sep="")Correlação da Estatura com MCT: 0.79 (r² = 0.62)
Correlação da Estatura com Sexo: 0.63 (r² = 0.4)
Correlação da Estatura com Idade: 0.11 (r² = 0.01)
Correlação da MCT com Sexo: 0.7 (r² = 0.5)
Correlação da MCT com Idade: 0.34 (r² = 0.12)
Correlação do Sexo com Idade: 0.26 (r²=0.07)
Teste de Correlação parcial
(controlando pelo efeito linear de Sexo e Idade):
estimate p.value statistic n gp Method
1 0.79 4.4e-39 17 178 1 pearson
(rXY.Z1Z2)² = 0.79^2 = 0.62É difícil entender o porquê da correlação parcial 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. # Correlação de Pearson sem dados brutos
- Teste de correlação de Pearson para uma condição
alfa <- 0.05
n <- 100
r <- 0.5
# H₀: ρ = 0 vs. H₁: ρ ≠ 0
test <- DescTools::CorCI(rho=r,
n=n,
conf.level=1-alfa,
alternative="two.sided")
print(test, digits=3) cor lwr.ci upr.ci
0.500 0.337 0.634 # Teste t
t <- r*sqrt((n-2)/(1-r^2))
p <- 2*pt(-abs(t), df=n-2)
cat("\nt(",n-2,") = ",round(t,2),", p = ",p,"\n\n",sep="")
t(98) = 5.72, p = 1.180492e-07# Teste z
z <- sqrt((n-3)/2)*log(1/((1-r)/(1+r)))
p <- 2*pnorm(-abs(z))
cat("z = ",round(z,2),", p = ",p, sep="")z = 7.65, p = 1.995011e-14# H₀: ρ = ρ₀ vs. H₁: ρ ≠ ρ₀
test_r_rho0 <- function(r, n, rho0 = 0, alternative = c("two.sided","less","greater"),
conf.level = 0.95){
alternative <- match.arg(alternative)
z <- atanh(r); z0 <- atanh(rho0); se <- 1/sqrt(n - 3)
zstat <- (z - z0)/se
p <- switch(alternative,
"two.sided" = 2*pnorm(-abs(zstat)),
"greater" = 1 - pnorm(zstat),
"less" = pnorm(zstat))
zalpha <- qnorm(0.5 + conf.level/2)
ci <- tanh(c(z - zalpha*se, z + zalpha*se))
list(z = zstat, p = p, conf.int = ci, conf.level = conf.level)
}
# exemplo:
test_r_rho0(r = 0.5, n = 100, rho0 = 0.10)$z
[1] 4.42185
$p
[1] 9.78596e-06
$conf.int
[1] 0.3366433 0.6341398
$conf.level
[1] 0.95- Correlação de Pearson parcial
rXY <- 0.40
rXZ <- 0.48
rYZ <- 0.85
# Metodo 1
r_XY.Z <- (rXY - rXZ*rYZ) / sqrt((1 - rXZ^2)*(1 - rYZ^2))
print(r_XY.Z, digits=3)[1] -0.0173# Metodo 2
R <- matrix(c(1, rXY, rXZ,
rXY, 1, rYZ,
rXZ, rYZ, 1),
nrow = 3, byrow = TRUE)
P <- solve(R) # matriz de precisao
r_XY.Z <- -P[1,2] / sqrt(P[1,1]*P[2,2]) # parcial de X–Y|Z
print(r_XY.Z, digits=3)[1] -0.0173- Correlações de Pearson independentes
nF <- 38
corrF <- 0.637
nM <- 51
corrM <- 0.646
out <- psych::r.test(n=nF, n2=nM, corrF, corrM)
cat("z = ",round(out$z,2),", p = ",round(out$p,4),"\n",sep="")z = 0.07, p = 0.9451- Correlações de Pearson dependentes
t = 0.81, p = 0.4227Um novo coeficiente de correlação
Conforme Chatterjee (2021),
RESUMO: É possível definir um coeficiente de correlação que seja (a) tão simples quanto os coeficientes clássicos, como o de Pearson ou o de Spearman, e ainda (b) estime de forma consistente alguma medida simples e interpretável do grau de dependência entre as variáveis — que seja igual a 0 se, e somente se, as variáveis forem independentes, e igual a 1 se, e somente se, uma for uma função mensurável da outra — e (c) possua uma teoria assintótica simples sob a hipótese de independência, como os coeficientes clássicos? Este artigo responde afirmativamente a essa questão, apresentando tal coeficiente. Nenhuma suposição é necessária sobre as distribuições das variáveis. Há diversos coeficientes na literatura que convergem para 0 se, e somente se, as variáveis forem independentes, mas nenhum que satisfaça as demais propriedades mencionadas acima. Materiais suplementares para este artigo estão disponíveis online.
Esta nova proposta de coeficiente de correlação está implementada no
pacote R denominado XICOR.
## ------------------------------------------------------------
## Mosaico dos 3 casos + tabela resumo (r, s, xi_n, p)
## ------------------------------------------------------------
set.seed(123)
# Utilitário: ξ_n e p (robusto a nomes do p-valor)
safe_xi <- function(x, y){
out <- try(XICOR::xicor(x, y, pvalue = TRUE), silent = TRUE)
if (inherits(out,"try-error") || is.null(out)) return(list(xi=NA_real_, p=NA_real_))
xi <- if (is.list(out)) suppressWarnings(as.numeric(out$xi))[1] else suppressWarnings(as.numeric(out))[1]
p <- NA_real_
if (is.list(out)) {
if (!is.null(out$p.value)) p <- suppressWarnings(as.numeric(out$p.value))[1]
if ((is.na(p) || is.null(p)) && !is.null(out$pval)) p <- suppressWarnings(as.numeric(out$pval))[1]
if ((is.na(p) || is.null(p)) && !is.null(out$p)) p <- suppressWarnings(as.numeric(out$p))[1]
}
list(xi=xi, p=p)
}
# Dados base
x <- -10:10
y_parab <- x^2
y_up <- x
y_down <- -x
# Painel 1×3
oldpar <- par(no.readonly=TRUE); on.exit(par(oldpar))
par(mfrow=c(1,3), mar=c(3,3,3,1))
plot(x, y_parab, pch=1, main="(a) y = x^2", xlab="x", ylab="y")
r <- cor(x, y_parab); s <- cor(x, y_parab, method="spearman"); xi <- safe_xi(x, y_parab)
legend("topleft", bty="n",
legend=sprintf("r=%.2f s=%.2f ξ_n=%.2f p=%s",
r, s, xi$xi, ifelse(is.na(xi$p),"NA", formatC(xi$p, format="g", digits=3))))
plot(x, y_up, pch=1, main="(b) y = x", xlab="x", ylab="y")
r <- cor(x, y_up); s <- cor(x, y_up, method="spearman"); xi <- safe_xi(x, y_up)
legend("topleft", bty="n",
legend=sprintf("r=%.2f s=%.2f ξ_n=%.2f p=%s",
r, s, xi$xi, ifelse(is.na(xi$p),"NA", formatC(xi$p, format="g", digits=3))))
plot(x, y_down, pch=1, main="(c) y = -x", xlab="x", ylab="y")
r <- cor(x, y_down); s <- cor(x, y_down, method="spearman"); xi <- safe_xi(x, y_down)
legend("topleft", bty="n",
legend=sprintf("r=%.2f s=%.2f ξ_n=%.2f p=%s",
r, s, xi$xi, ifelse(is.na(xi$p),"NA", formatC(xi$p, format="g", digits=3))))
# Tabela resumo
summ_row <- function(lbl, x, y){
r <- suppressWarnings(cor(x, y, method="pearson"))
s <- suppressWarnings(cor(x, y, method="spearman"))
xi <- safe_xi(x, y)
data.frame(caso=lbl, n=length(x), r=r, s=s, xi_n=xi$xi, p_xi=xi$p)
}
TAB <- rbind(
summ_row("(a) y = x^2", x, y_parab),
summ_row("(b) y = x", x, y_up),
summ_row("(c) y = -x", x, y_down)
)
TAB$p_xi <- signif(TAB$p_xi, 3)
print(TAB, row.names = FALSE, digits = 3) caso n r s xi_n p_xi
(a) y = x^2 21 0 0 0.727 1.00e-07
(b) y = x 21 1 1 0.864 2.39e-10
(c) y = -x 21 -1 -1 0.864 2.39e-10
set.seed(123)
## Gerador não-normal: x ~ U(-1,1), y = 0.7*x + ruído uniforme
gen_sample <- function(n){
x <- runif(n, -1, 1)
y <- 0.7*x + runif(n, -0.5, 0.5) # relação linear + ruído NÃO-normal
list(x=x, y=y)
}
## Duplicar os dados não muda r
df <- data.frame(gen_sample(20))
sunflowerplot(scale(y)~scale(x), pch=21, data=df)
lines(lowess(scale(df$x),scale(df$y)))
$xi
[1] 0.4586466
$sd
[1] 0.1422166
$pval
[1] 0.0006298935$xi
[1] 0.4586466
$sd
[1] 0.1422166
$pval
[1] 0.0006298935O coeficiente de correlação \(\xi_n\) não uma estatística de tamanho de efeito.
set.seed(123)
## Gerador não-normal: x ~ U(-1,1), y = 0.7*x + ruído uniforme
gen_sample <- function(n){
x <- runif(n, -1, 1)
y <- 0.7*x + runif(n, -0.5, 0.5) # relação linear + ruído NÃO-normal
list(x=x, y=y)
}
## Duplicar os dados não muda r
df <- data.frame(gen_sample(20))
sunflowerplot(scale(y)~scale(x), pch=21, data=df)
lines(lowess(scale(df$x),scale(df$y)))
xin <- XICOR::xicor(df$x, df$y, pvalue = TRUE)
df2 <- rbind(df,df)
sunflowerplot(scale(y)~scale(x), pch=21, data=df2)
lines(lowess(scale(df2$x),scale(df2$y)))
xin_dup <- XICOR::xicor(df2$x, df2$y, pvalue = TRUE) # duplica pares => mesma informação
cat(sprintf("ksi_n(n=20)=%.4f | ksi_n(n=40 duplicado)=%.4f\n", xin$xi, xin_dup$xi))ksi_n(n=20)=0.4586 | ksi_n(n=40 duplicado)=0.7293Exemplo: As ervilhas de Galton revisitado
Sir Francis Galton coletou, em 1875, um dos conjuntos de dados mais
antigos e famosos da estatística. O banco contém 700 observações dos
diâmetros médios das ervilhas-de-cheiro nas plantas-mãe e nas
plantas-filha. O processo exato de coleta não foi devidamente
registrado; sabe‑se apenas que Galton enviou pacotes de sementes a
amigos, que plantaram as sementes, cultivaram as plantas e enviaram as
sementes das novas plantas de volta a Galton (ver Stigler 1986, p. 296,
para mais detalhes). O conjunto de dados está disponível em
psychTools::peas.
Seja \(X\) o diâmetro médio das ervilhas na planta‑mãe e \(Y\) o diâmetro médio na planta‑filha. Como já observado por Pearson há muito tempo, a correlação entre \(X\) e \(Y\) é cerca de 0.35. Os \(X_i\) apresentam muitos empates, o que torna \(\xi_n(X,Y)\) aleatório devido ao desempate randômico. A média de 10000 simulações forneceu um valor próximo de 0.11 para \(\xi_n(X,Y)\). O valor‑p para o teste de independência, usando os Teoremas 2.2 e 2.3, foi menor que 0.0001; portanto, \(\xi_n(X,Y)\) teve sucesso em detectar dependência entre \(X\) e \(Y\).
Até aqui, nada surpreendente. A verdadeira surpresa foi que o valor
de \(\xi_n(Y,X)\) (em vez de \(\xi_n(X,Y)\)) ficou em aproximadamente 0.92
e pareceu ser independente do procedimento de desempate. Pelo Teorema
1.1, isso significa que \(X\) está
próximo de ser uma função sem ruído de \(Y\). Pelo diagrama de dispersão dos dados
(gráfico de dispersão abaixo por sunflowerplot), isso não é
evidente. O mistério se resolve examinando a tabela de contingência dos
dados (Tabela abaixo por table). Cada linha da tabela
corresponde a um valor de \(Y\) e cada
coluna corresponde a um valor de \(X\).
Observa‑se que cada coluna tem várias células com contagens não nulas, o
que significa que, para cada valor de \(X\), há muitos valores diferentes de \(Y\) nos dados. Por outro lado, cada linha
contém exatamente uma célula com contagem não nula (muitas vezes
bastante grande). Ou seja, para qualquer valor de \(Y\), o valor de \(X\) nos dados é sempre o mesmo.
Por exemplo, entre as plantas‑mãe com diâmetro médio 15, houve 46 casos em que a planta‑filha tinha diâmetro 13.77, 14 casos com 14.77, 11 casos com 16.77, 14 casos com 17.77 e 4 casos com 18.77. Por outro lado, para todas as 46 plantas‑filha com diâmetro 13.77, as plantas‑mãe tinham diâmetro 15. De modo semelhante, para todas as 34 plantas‑filha com diâmetro 14.28, as plantas‑mãe tinham diâmetro 16.
O bom senso sugere que a razão por trás desse fenômeno estranho é algum artifício do método de coleta ou registro dos dados, e não um fato biológico profundo. Provavelmente não é apenas um efeito de arredondamento simples; por exemplo, em todos os 46 casos em que \(Y=13.77\), temos \(X=15\), mas em todos os 37 casos em que \(Y=13.92\) (apenas um pouco diferente de 13.77), temos \(X=17\). Entretanto, se admitirmos que os valores registrados nos dados são exatamente os valores medidos e que as observações são IID (nenhuma das suposições é estritamente verdadeira), então, olhando a tabela , não há como escapar da conclusão de que o diâmetro médio das ervilhas na planta‑mãe pode ser previsto com considerável certeza a partir do diâmetro médio das ervilhas na planta‑filha (mas não o contrário). O coeficiente \(\xi_n(Y,X)\) revela numericamente esse fato ao assumir valor próximo de 1. É provável que essa característica dos dados de Galton já tenha sido notada antes; se foi, é certamente difícil de encontrar. Não localizei referência que a mencione, apesar de muita procura.
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
Dados <- psychTools::peas
sunflowerplot(scale(Dados$parent), scale(Dados$child),
main = "Galton: diâmetros médios (mãe vs filha)",
xlab = "X = parent", ylab = "Y = child") # Figura 1
lines(lowess(scale(Dados$parent), scale(Dados$child)), lwd=2, col="black")
child
parent 13.77 13.92 14.07 14.28 14.35 14.66 14.67 14.77 14.92 15.07 15.28 15.35
15 46 0 0 0 0 0 0 14 0 0 0 0
16 0 0 0 34 0 0 0 0 0 0 15 0
17 0 37 0 0 0 0 0 0 16 0 0 0
18 0 0 0 0 34 0 0 0 0 0 0 12
19 0 0 35 0 0 0 0 0 0 16 0 0
20 0 0 0 0 0 23 0 0 0 0 0 0
21 0 0 0 0 0 0 22 0 0 0 0 0
child
parent 15.66 15.67 15.77 15.92 16.07 16.28 16.35 16.66 16.67 16.77 16.92 17.07
15 0 0 9 0 0 0 0 0 0 11 0 0
16 0 0 0 0 0 18 0 0 0 0 0 0
17 0 0 0 13 0 0 0 0 0 0 16 0
18 0 0 0 0 0 0 13 0 0 0 0 0
19 0 0 0 0 12 0 0 0 0 0 0 13
20 10 0 0 0 0 0 0 12 0 0 0 0
21 0 8 0 0 0 0 0 0 10 0 0 0
child
parent 17.28 17.35 17.66 17.67 17.77 17.92 18.07 18.28 18.35 18.66 18.67 18.77
15 0 0 0 0 14 0 0 0 0 0 0 4
16 16 0 0 0 0 0 0 13 0 0 0 0
17 0 0 0 0 0 13 0 0 0 0 0 0
18 0 17 0 0 0 0 0 0 16 0 0 0
19 0 0 0 0 0 0 11 0 0 0 0 0
20 0 0 17 0 0 0 0 0 0 20 0 0
21 0 0 0 18 0 0 0 0 0 0 21 0
child
parent 18.92 19.07 19.28 19.35 19.66 19.67 19.77 19.92 20.07 20.28 20.35 20.66
15 0 0 0 0 0 0 2 0 0 0 0 0
16 0 0 3 0 0 0 0 0 0 1 0 0
17 4 0 0 0 0 0 0 1 0 0 0 0
18 0 0 0 6 0 0 0 0 0 0 2 0
19 0 10 0 0 0 0 0 0 2 0 0 0
20 0 0 0 0 13 0 0 0 0 0 0 3
21 0 0 0 0 0 13 0 0 0 0 0 0
child
parent 20.67 22.07 22.66 22.67
15 0 0 0 0
16 0 0 0 0
17 0 0 0 0
18 0 0 0 0
19 0 1 0 0
20 0 0 2 0
21 6 0 0 2
Pearson's product-moment correlation
data: Dados$parent and Dados$child
t = 9.7536, df = 698, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2793996 0.4099146
sample estimates:
cor
0.3463319
Spearman's rank correlation rho
data: Dados$parent and Dados$child
S = 36495384, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.3615955 rho lwr.ci upr.ci
0.3615955 0.2954082 0.4243283 $xi
[1] 0.1153208
$sd
[1] 0.02408662
$pval
[1] 8.432983e-07 parent child
parent 0.99625 0.1056318
child 0.92250 0.9958986
Referências
- Archer, KJ et al. (2008) A disattenuated correlation estimate when variables are measured with error: Illustration estimating cross-platform correlations. Statistics in Medicine 27(7): 1026–39. https://doi.org/10.1002/sim.2984
- Betensky, RA (2019) The p-Value Requires Context, Not a Threshold. The American Statistician, 73(sup1): 115-117. DOI:10.1080/00031305.2018.1529624
- Bhagat, V (2018) Aspiring Algo Trader Tech Enthusiast at Self-Employment. Traduzido e modificado de https://www.quora.com/What-is-the-difference-between-correlation-analysis-and-regression-analysis
- Chatterjee, S (2021) A new coefficient of correlation. Journal of the American Statistical Association 116(536): 2009-22. https://doi.org/10.1080/01621459.2020.1758115
- Dancey, CP & Reidy, J (2019) Estatística sem Matemática para Psicologia 7ª ed. Porto Alegre: Penso.
- Ellis, PD (2010) The essential guide to effect sizes. UK: Cambridge.
- Erdem, O et al. (2014) A new correlation coefficient for bivariate time-series data. Physica A 414: 274–84.
- Galton, F (1886) Regression towards mediocrity in hereditary stature. The Journal of the Anthropological Institute of Great Britain and Ireland 15: 246–63. https://doi.org/10.2307/2841583
- Galton, F (1888) Co-relations and their measurement, chiefly from anthropometric data. Proceedings of the Royal Society of London 45: 135–45.
- Jones, M et al. (2025) RESI: An R package for robust effect sizes. Journal of Statistical Software 112(3). doi: 10.18637/jss.v112.i03.
- Judkins, DR & Porter, KE (2015) Robustness of ordinary least squares in randomized clinical trials. Statistics in Medicine 35(11): 1763-73. doi: 10.1002/sim.6839.
- Pearson, K (1896) Mathematical contributions to the theory of evolution. III. Regression, heredity and panmixia. Philosophical Transactions of the Royal Society of London. Series A 187: 253–318.
- Mohammed, N (2018) Traduzido e modificado de https://www.researchgate.net/post/What_is_the_key_differences_between_correlation_and_regression. University of Anbar.
- Revelle, W (2014) Personality Project. http://personality-project.org/revelle.html.
- Shinohara, EMG (1989) Prevalência de anemia em gestantes de primeira consulta em centros de saúde do estado no Subdistrito de Paz do Butantã, Município de São Paulo [dissertação]. SP: USP, Faculdade de Ciências Farmacêuticas. doi:10.11606/D.9.1989.tde-27032008-142216. Os dados utilizados nos exemplos deste capítulo são dados parciais relativos a este trabalho, gentilmente fornecidos pelo Prof. Raymundo Soares de Azevedo Neto, Departamento de Patologia, Faculdade de Medicina da USP.