Análise de Correlação em R
MSP4112: Estatística Aplicada à Medicina
27 outubro 2026 12:46h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(bootES, warn.conflicts=FALSE))
suppressMessages(library(bruceR, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(CovCorTest, 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(metafor, warn.conflicts=FALSE))
suppressMessages(library(ppcor, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(psychTools, warn.conflicts=FALSE))
suppressMessages(library(qgraph, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(TOSTER, 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.RCorrelacaoPearsonUmaCorrelacao.RCorrelacaoPearsonUmaCorrelacaoRho0.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 (1888) 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 (1886) 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, asp=1, main="Human height")
lines(lowess(Dados_height$parent, Dados_height$child))
curve(x*1, add=TRUE, lty=2)
Call:
lm(formula = child ~ parent, data = Dados_height)
Residuals:
Min 1Q Median 3Q Max
-7.805 -1.366 0.049 1.634 5.926
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.9415 2.8109 8.52 <2e-16 ***
parent 0.6463 0.0411 15.71 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.24 on 926 degrees of freedom
Multiple R-squared: 0.21, Adjusted R-squared: 0.21
F-statistic: 247 on 1 and 926 DF, p-value: <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-16 parent child
parent 1.0000 0.4588
child 0.4588 1.0000Dados_peas <- psychTools::peas # parent and child generation of 700 sweet peas Galton (1886)
sunflowerplot(data=Dados_peas, child~parent, asp=1, main="Sweet peas")
lines(lowess(Dados_peas$parent, Dados_peas$child))
curve(x*1, add=TRUE, lty=2)
Call:
lm(formula = child ~ parent, data = Dados_peas)
Residuals:
Min 1Q Median 3Q Max
-2.644 -1.505 -0.311 1.513 5.689
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.1143 0.6366 15.89 <2e-16 ***
parent 0.3429 0.0352 9.75 <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.12, Adjusted R-squared: 0.119
F-statistic: 95.1 on 1 and 698 DF, p-value: <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-16 parent child
parent 1.0000 0.3463
child 0.3463 1.0000Dados_cubit <- psychTools::heights # height and cubit Galton (1888)
sunflowerplot(data=Dados_cubit, height~cubit, main="Human height and cubit")
lines(lowess(Dados_cubit$cubit, Dados_cubit$height))
Call:
lm(formula = height ~ cubit, data = Dados_cubit)
Residuals:
Min 1Q Median 3Q Max
-4.426 -1.152 -0.152 1.300 4.985
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.928 1.929 13.4 <2e-16 ***
cubit 2.274 0.107 21.4 <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.568, Adjusted R-squared: 0.567
F-statistic: 456 on 1 and 346 DF, p-value: <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-16 height cubit
height 1.000 0.754
cubit 0.754 1.000Pearson (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 estatura de pais e filhos, massa corporal total e comprimentos das partes do corpo humano etc.
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:
O coeficiente de correlação é a inclinação da reta de regressão padronizado;
A inclinação da reta 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: Gestante do Centro de Saúde Geraldo de Paula Souza
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. O grupo selecionado para análise o de gestantes. Usaremos apenas algumas destas variáveis do início do pré-natal:
HT… Hematócrito (%)HB… Hemoglobina (g/dl)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.
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
individualizados, 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 <- 1e3
# 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, 1000 resamples
data: Hematócrito (%) and Hemoglobina (mg/dl)
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6795834 0.9414803
sample estimates:
r (bootstrap)
0.8682372 
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, 1000 resamples
data: Hematócrito (%) and Hemácias (milhão/mm3)
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4265115 0.8510835
sample estimates:
r (bootstrap)
0.6879649 
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, 1000 resamples
data: Hematócrito (%) and Leucócitos (milhar/mm3)
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1418519 0.3959828
sample estimates:
r (bootstrap)
0.1203115 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
- Em tabela de contigência, independe da dimensão
- 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 variável dependente (VD) e variável explicativa (VE).
Ambos são medidas de tamanho de efeito, pois:
- Não dependem do tamanho da amostra;
- São adimensionais;
- Variam entre 0 e 1.
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"
)
# MATRIZ DE DISPERSÃO
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::ggpairs(Dadoz))
# CORRELAÇÕES DE PEARSON
cat("\nMatriz de correlações de Pearson:\n")
print(
R.pearson <- cor(Dados),
digits = 3
)
# GGally: Pearson
print(GGally::ggcorr(
Dadoz,
name = "Pearson",
nbreaks = 6,
digits = 3,
label = TRUE,
label_size = 4,
color = "#888888"
))
# bruceR: Pearson com Holm
bruceR::Corr(
Dados,
method = "pearson",
digits = 3,
p.adjust = "holm",
plot = FALSE
)
# CORRELAÇÕES DE SPEARMAN
cat("\nMatriz de correlações de Spearman:\n")
print(
R.spearman <- cor(Dados, method = "spearman"),
digits = 3
)
# GGally: Spearman
print(GGally::ggcorr(
Dadoz,
name = "Spearman",
method = c("complete.obs", "spearman"),
nbreaks = 6,
digits = 3,
label = TRUE,
label_size = 4,
color = "#888888"
))
# bruceR: Spearman com Holm
bruceR::Corr(
Dados,
method = "spearman",
digits = 3,
p.adjust = "holm",
plot = FALSE
)
# QGRAPH: COMPARAÇÃO PEARSON E SPEARMAN
# Layout comum calculado pela média
# das magnitudes das duas matrizes
R.media <- (
abs(R.pearson) + abs(R.spearman)
) / 2
set.seed(123)
G.layout <- qgraph::qgraph(
R.media,
layout = "spring",
DoNotPlot = TRUE
)$layout
# Gráficos lado a lado
op <- par(mfrow = c(1, 2))
qgraph::qgraph(
R.pearson,
layout = G.layout,
labels = colnames(Dados),
weighted = TRUE,
maximum = 1,
minimum = 0,
cut = 0,
edge.labels = TRUE,
edge.label.cex = 0.8,
theme = "colorblind",
title = "Pearson"
)
qgraph::qgraph(
R.spearman,
layout = G.layout,
labels = colnames(Dados),
weighted = TRUE,
maximum = 1,
minimum = 0,
cut = 0,
edge.labels = TRUE,
edge.label.cex = 0.8,
theme = "colorblind",
title = "Spearman"
)
par(op)
Matriz de correlações de Pearson:
IDADE HT HB HEM LEUC FOLICO B12
IDADE 1.00000 0.00723 -0.0492 0.0144 0.2125 0.0464 -0.0293
HT 0.00723 1.00000 0.8645 0.6862 0.1141 0.2526 -0.0537
HB -0.04923 0.86453 1.0000 0.5870 0.0351 0.0831 -0.1100
HEM 0.01442 0.68621 0.5870 1.0000 0.1681 0.2594 0.2619
LEUC 0.21249 0.11415 0.0351 0.1681 1.0000 -0.0491 -0.1234
FOLICO 0.04636 0.25258 0.0831 0.2594 -0.0491 1.0000 0.2986
B12 -0.02934 -0.05373 -0.1100 0.2619 -0.1234 0.2986 1.0000
Pearson's r and 95% confidence intervals:
p values and 95% CIs are adjusted using the "holm" method.
─────────────────────────────────────────────────
r [95% CI] p N
─────────────────────────────────────────────────
IDADE-HT 0.007 [-0.305, 0.318] 1.000 40
IDADE-HB -0.049 [-0.455, 0.374] 1.000 40
IDADE-HEM 0.014 [-0.340, 0.365] 1.000 40
IDADE-LEUC 0.212 [-0.257, 0.601] 1.000 40
IDADE-FOLICO 0.046 [-0.360, 0.438] 1.000 40
IDADE-B12 -0.029 [-0.399, 0.349] 1.000 40
HT-HB 0.865 [ 0.670, 0.948] <.001 *** 40
HT-HEM 0.686 [ 0.331, 0.871] <.001 *** 40
HT-LEUC 0.114 [-0.338, 0.523] 1.000 40
HT-FOLICO 0.253 [-0.221, 0.630] 1.000 40
HT-B12 -0.054 [-0.465, 0.376] 1.000 40
HB-HEM 0.587 [ 0.177, 0.823] .001 ** 40
HB-LEUC 0.035 [-0.359, 0.418] 1.000 40
HB-FOLICO 0.083 [-0.356, 0.492] 1.000 40
HB-B12 -0.110 [-0.517, 0.337] 1.000 40
HEM-LEUC 0.168 [-0.296, 0.568] 1.000 40
HEM-FOLICO 0.259 [-0.217, 0.636] 1.000 40
HEM-B12 0.262 [-0.217, 0.639] 1.000 40
LEUC-FOLICO -0.049 [-0.449, 0.367] 1.000 40
LEUC-B12 -0.123 [-0.534, 0.334] 1.000 40
FOLICO-B12 0.299 [-0.182, 0.664] 1.000 40
─────────────────────────────────────────────────
p values ABOVE the diagonal are adjusted using the "holm" method.
Matriz de correlações de Spearman:
IDADE HT HB HEM LEUC FOLICO B12
IDADE 1.0000 -0.0249 -0.0582 0.0596 0.0748 0.0504 -0.0504
HT -0.0249 1.0000 0.7803 0.6113 0.1906 0.2668 -0.0250
HB -0.0582 0.7803 1.0000 0.5165 0.1101 -0.0209 -0.1562
HEM 0.0596 0.6113 0.5165 1.0000 0.1868 0.2446 0.3605
LEUC 0.0748 0.1906 0.1101 0.1868 1.0000 0.0267 -0.2074
FOLICO 0.0504 0.2668 -0.0209 0.2446 0.0267 1.0000 0.3174
B12 -0.0504 -0.0250 -0.1562 0.3605 -0.2074 0.3174 1.0000
Spearman's r and 95% confidence intervals:
p values and 95% CIs are adjusted using the "holm" method.
─────────────────────────────────────────────────
r [95% CI] p N
─────────────────────────────────────────────────
IDADE-HT -0.025 [-0.374, 0.331] 1.000 40
IDADE-HB -0.058 [-0.463, 0.366] 1.000 40
IDADE-HEM 0.060 [-0.371, 0.469] 1.000 40
IDADE-LEUC 0.075 [-0.364, 0.486] 1.000 40
IDADE-FOLICO 0.050 [-0.357, 0.441] 1.000 40
IDADE-B12 -0.050 [-0.450, 0.366] 1.000 40
HT-HB 0.780 [ 0.498, 0.913] <.001 *** 40
HT-HEM 0.611 [ 0.211, 0.836] <.001 *** 40
HT-LEUC 0.191 [-0.275, 0.584] 1.000 40
HT-FOLICO 0.267 [-0.209, 0.641] 1.000 40
HT-B12 -0.025 [-0.396, 0.353] 1.000 40
HB-HEM 0.517 [ 0.077, 0.788] .012 * 40
HB-LEUC 0.110 [-0.337, 0.517] 1.000 40
HB-FOLICO -0.021 [-0.330, 0.293] 1.000 40
HB-B12 -0.156 [-0.554, 0.300] 1.000 40
HEM-LEUC 0.187 [-0.275, 0.578] 1.000 40
HEM-FOLICO 0.245 [-0.229, 0.624] 1.000 40
HEM-B12 0.360 [-0.114, 0.701] .402 40
LEUC-FOLICO 0.027 [-0.366, 0.411] 1.000 40
LEUC-B12 -0.207 [-0.598, 0.262] 1.000 40
FOLICO-B12 0.317 [-0.159, 0.674] .781 40
─────────────────────────────────────────────────
p values ABOVE the diagonal are adjusted using the "holm" method.
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=100\) e \(r=0.2\), 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 <- 100
r <- 0.2
# 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="")
rd <- DescTools::CorCI(r, n=n, conf.level = 1-alfa,
alternative = "two.sided")
print(rd, 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="")
rd <- DescTools::CorCI(r, n=n, conf.level = 1-alfa,
alternative = "greater")
print(rd, 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="")
rd <- DescTools::CorCI(r, n=n, conf.level = 1-alfa,
alternative = "less")
print(rd, digits=2)Teste bilateral:
H₀: ρ = 0 vs. H₁: ρ ≠ 0
t(98) = 2.02 , p = 0.046
cor lwr.ci upr.ci
0.2000 0.0037 0.3814
Teste unilateral à direita:
H₀: ρ = 0 vs. H₁: ρ > 0
t(98) = 2.02 , p = 0.023
cor lwr.ci upr.ci
0.200 0.036 1.000
Teste unilateral à esquerda:
H₀: ρ = 0 vs. H₁: ρ < 0
t(98) = 2.02 , p = 0.977
cor lwr.ci upr.ci
0.20 -1.00 0.35 Observe que, mesmo sem os dados individualizados, é 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.2\) é significantemente diferente de zero no teste bilateral (zero não está contido no intervalo de confiança 95%) e o intervalo de confiança é assimétrico.
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.2 é 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 <- 100
r <- 0.2
cat("\n")
cat("Teste de correlação:\n")
cat("\tH0: rho = 0 vs. H1: rho != 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) * atanh(r)
p <- 2 * (1 - 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_t <- 2 * (1 - pt(abs(t), df = n - 2))
cat("t(", n - 2, ") = ", round(t, 2), ", p = ",
round(p_t, 4), "\n", sep = "")
Teste de correlação:
H0: rho = 0 vs. H1: rho != 0
cor lwr.ci upr.ci
0.20000 0.00373 0.38143
z = 2, p = 0.0459
Alternativamente:
t(98) = 2.02, p = 0.046O teste da hipótese nula bilateral pode realizado para qualquer valor de correlação de Person populacional, 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 <- 100
r <- 0.7
ro <- 0.6
cat("r = ", r, "\n", sep = "")
cat("Teste de correlação:\n")
cat("\tH₀: ρ = ", ro, " vs. H₁: ρ ≠ ", ro, "\n", sep = "")
z <- sqrt(n - 3) * (atanh(r) - atanh(ro))
p <- 2 * (1 - pnorm(abs(z)))
cat("z = ", round(z, 2), ", p = ", round(p, 4), sep = "")
zr <- atanh(r)
se_z <- 1 / sqrt(n - 3)
z_crit <- qnorm(1 - alfa/2)
li_z <- zr - z_crit * se_z
ls_z <- zr + z_crit * se_z
li_r <- tanh(li_z)
ls_r <- tanh(ls_z)
cat("\nIC95% bilateral para ρ: [",
round(li_r, 4), ", ", round(ls_r, 4), "]\n", sep = "")r = 0.7
Teste de correlação:
H₀: ρ = 0.6 vs. H₁: ρ ≠ 0.6
z = 1.72, p = 0.0863
IC95% bilateral para ρ: [0.5839, 0.7881]O teste da hipótese nula unilateral pode realizado para qualquer valor de correlação de Person populacional, e.g., \(\rho = 0.6\):
\[ \begin{cases} H_0: \rho \le 0.6\\ H_1: \rho > 0.6 \end{cases}\\ \alpha = 0.05 \]
Podemos usar o teste implementado em CorrelacaoPearsonUmaCorrelacaoRhoQualquerUnilateral.R:
Pearson's product-moment correlation
data: x and y
z = 1.7152, N = 100, p-value = 0.04315
alternative hypothesis: true correlation is greater than 0.6
90 percent confidence interval:
0.6045524 0.7756312
sample estimates:
cor
0.7 cor lwr.ci upr.ci
0.7000 0.6046 1.0000 TESTE UNILATERAL DE CORRELAÇÃOn = 100; r = 0.700H0: ρ <= 0.600 vs. H1: ρ > 0.600z = 1.7152; valor-p = 0.04315IC95% unilateral para ρ: [0.6046, 1]Decisão: Rejeitar H0Ao nível de significância de 5%, há evidência estatisticamente significante de que a correlação linear populacional é superior a 0.600. Entretanto, a evidência é limítrofe, pois o valor-p está próximo de 0.05 e o limite inferior do intervalo de confiança supera 0.600 por apenas 0.0046.
Quando temos os dados individualizados, 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:
Dados <- readRDS("Gestante.rds")
B <- 1e3
cat("Hematócrito e Hemoglobina:")
ic <- bootES::bootES(Dados[c("HT","HB")], R=B)
print(ic)
cat("Hematócrito e Hemácias:")
ic <- bootES::bootES(Dados[c("HT","HEM")], R=B)
print(ic)
cat("Hematócrito e Leucócitos:")
ic <- bootES::bootES(Dados[c("HT","LEUC")], R=B)
print(ic)Hematócrito e Hemoglobina:
95.00% bca Confidence Interval, 1000 replicates
Stat CI (Low) CI (High) bias SE
0.865 0.721 0.942 -0.009 0.061
Hematócrito e Hemácias:
95.00% bca Confidence Interval, 1000 replicates
Stat CI (Low) CI (High) bias SE
0.686 0.406 0.856 -0.011 0.113
Hematócrito e Leucócitos:
95.00% bca Confidence Interval, 1000 replicates
Stat CI (Low) CI (High) bias SE
0.114 -0.176 0.377 0.015 0.141 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 <- 1e3
rIC95 <- bootES::bootES(Dados[c("Estatura","MCT","Sexo")],
group.col="Sexo",
R=B)
print(rIC95)
95.00% bca Confidence Interval, 1000 replicates
Stat CI (Low) CI (High) bias SE
-0.009 -0.326 0.240 0.003 0.138 Como o zero está contido no intervalo dado pelo teste, não rejeitamos \(H_0\).
Existe alternativa para os casos em que os dados individualizados 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.
Comparação de três ou mais correlações independentes
Quando uma investigação envolve três ou mais grupos independentes, pode ser necessário verificar se a intensidade da associação linear entre duas variáveis é homogênea entre os grupos.
A comparação individual das correlações por testes pareados aumenta o número de testes e exige procedimentos de controle da multiplicidade. Uma alternativa é realizar inicialmente um teste global (omnibus) da igualdade das correlações populacionais.
Teste Q de Cochran para homogeneidade de correlações
O teste Q de Cochran, aplicado aos coeficientes de correlação de Pearson transformados por Fisher, permite avaliar simultaneamente a homogeneidade de \(k\) correlações independentes.
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.
alfa <- 0.05
# dados: três espécies
I3 <- iris
I3$Species <- droplevels(I3$Species)
I.S <- subset(I3, Species=="setosa")
I.Vc <- subset(I3, Species=="versicolor")
I.Vg <- subset(I3, Species=="virginica")
print(summary(I3[, c("Sepal.Length","Sepal.Width","Species")])) 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## ===== IC das correlações de Pearson por espécie (cor.test) =====
alfa <- 0.05
g <- nlevels(I3$Species) # 3 espécies
alfaBonf <- alfa/g # opcional (Bonferroni)
conf_bonf <- 1 - alfaBonf
# cor.test por espécie
ctS <- with(I.S , cor.test(Sepal.Length, Sepal.Width, method="pearson", conf.level=conf_bonf))
ctVc <- with(I.Vc, cor.test(Sepal.Length, Sepal.Width, method="pearson", conf.level=conf_bonf))
ctVg <- with(I.Vg, cor.test(Sepal.Length, Sepal.Width, method="pearson", conf.level=conf_bonf))
rS <- unname(ctS$estimate)[1]; ciS <- unname(ctS$conf.int); nS <- sum(complete.cases(I.S$Sepal.Length, I.S$Sepal.Width))
rVc <- unname(ctVc$estimate)[1]; ciVc <- unname(ctVc$conf.int); nVc <- sum(complete.cases(I.Vc$Sepal.Length, I.Vc$Sepal.Width))
rVg <- unname(ctVg$estimate)[1]; ciVg <- unname(ctVg$conf.int); nVg <- sum(complete.cases(I.Vg$Sepal.Length, I.Vg$Sepal.Width))
TAB.ic <- data.frame(
Especie = c("setosa","versicolor","virginica"),
n = c(nS, nVc, nVg),
r = c(rS, rVc, rVg),
lo = c(ciS[1], ciVc[1], ciVg[1]),
hi = c(ciS[2], ciVc[2], ciVg[2])
)
cat("\nIC95% Bonferroni de correlação de Pearson\n")
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.687xlim_min <- max(-1, min(TAB.ic$lo) - 0.02)
xlim_max <- min( 1, max(TAB.ic$hi) + 0.02)
ypos <- 1:3
plot(NA, xlim=c(-.1, xlim_max), ylim=c(0.5, 3.5),
xlab="r",
ylab="", yaxt="n",
main="IC95% Bonferroni: SepalL×SepalW")
abline(v=0, lty=2, col="grey75")
# paleta segura para daltônicos
pal <- RColorBrewer::brewer.pal(3, "Dark2")
cols <- pal # já definido antes (Dark2)
segments(TAB.ic$lo, ypos, TAB.ic$hi, ypos, lwd=3, col=cols, lend="round")
points(TAB.ic$r, ypos, pch=19, cex=1.2, col=cols)
axis(2, at=ypos,
labels=paste0("n=", TAB.ic$n),
las=1, tick=FALSE)
# rótulo numérico opcional
text(TAB.ic$hi, ypos + 0.1, labels=sprintf("r=%.3f", TAB.ic$r), cex=0.9, pos=2)
car::dataEllipse(I3$Sepal.Length, I3$Sepal.Width,
grid=FALSE,
groups=factor(I3$Species),
group.labels=c("setosa", "versicolor", "virginica"),
col=c("black","blue", "red"),
levels=0.68,
robust=TRUE,
main="Elipses de predição de 68%",
xlab="Sepal Length (cm)",
ylab="Sepal Width (cm)",
lwd=2, lty=1)
I3$Sepal.Length.z <- unsplit(lapply(split(I3$Sepal.Length, I3$Species), scale),
I3$Species)
I3$Sepal.Width.z <- unsplit(lapply(split(I3$Sepal.Width, I3$Species), scale),
I3$Species)
df <- data.frame(I3$Sepal.Length.z,I3$Sepal.Width.z)
matriz <- as.matrix(df)
car::dataEllipse(matriz[,1], matriz[,2],
grid=FALSE,
groups=factor(I3$Species),
group.labels=c("setosa", "versicolor", "virginica"),
col=c("black","blue", "red"),
levels=0.68,
robust=TRUE,
main="Elipses de predição de 68%",
xlab="Sepal Length (z)",
ylab="Sepal Width (z)",
xlim=c(-3,3),
ylim=c(-3,3),
lwd=2, lty=1)
# Correlações por grupo e testes pareados (independentes)rS <- cor(I.S$Sepal.Width, I.S$Sepal.Length); nS <- nrow(I.S)
rVc <- cor(I.Vc$Sepal.Width, I.Vc$Sepal.Length); nVc <- nrow(I.Vc)
rVg <- cor(I.Vg$Sepal.Width, I.Vg$Sepal.Length); nVg <- nrow(I.Vg)
cat("cor(setosa) = ", round(rS ,3), "\n", sep="")cor(setosa) = 0.743cor(versicolor)= 0.526cor(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.80167 with probability 0.0716Correlation tests
Call:psych::r.test(n = nS, r12 = rS, r34 = rVg, n2 = nVg)
Test of difference between two independent correlations
z value 2.24124 with probability 0.02501Correlation tests
Call:psych::r.test(n = nVc, r12 = rVc, r34 = rVg, n2 = nVg)
Test of difference between two independent correlations
z value 0.43956 with probability 0.66025Comparação das correlações entre as espécies de Iris
Foram estimados os coeficientes de correlação de Pearson entre o comprimento e a largura da sépala, separadamente para as três espécies:
- Setosa: \(r=0.743\)
- Versicolor: \(r=0.526\)
- Virginica: \(r=0.457\)
As três correlações são positivas, sendo a associação linear mais intensa em setosa e menos intensa em virginica.
Testes de igualdade de correlações independentes
Foram realizados três testes baseados na transformação \(z\) de Fisher, com as hipóteses:
\[ H_0:\rho_i=\rho_j \qquad H_1:\rho_i\neq\rho_j \]
Os resultados foram:
| Comparação | Estatística \(z\) | Valor-p |
|---|---|---|
| Setosa × Versicolor | 1.80167 | 0.07160 |
| Setosa × Virginica | 2.24124 | 0.02501 |
| Versicolor × Virginica | 0.43956 | 0.66025 |
Adotando-se \(\alpha=0.05\), sem correção para múltiplas comparações, somente a diferença entre as correlações de setosa e virginica foi estatisticamente significante (\(p=0.02501\)).
Correção de Holm para múltiplas comparações
Considerando os três testes simultaneamente, aplicou-se a correção de Holm para controlar a probabilidade familiar de erro do tipo I (family-wise error rate, FWER) em 5%.
p <- c(
"Setosa vs Versicolor" = 0.07160,
"Setosa vs Virginica" = 0.02501,
"Versicolor vs Virginica" = 0.66025
)
p.holm <- stats::p.adjust(p, method = "holm")
print(round(p.holm, 5)) Setosa vs Versicolor Setosa vs Virginica Versicolor vs Virginica
0.14320 0.07503 0.66025 Os valores-p ajustados foram:
| Comparação | Valor-p ajustado (Holm) |
|---|---|
| Setosa × Versicolor | 0.14320 |
| Setosa × Virginica | 0.07503 |
| Versicolor × Virginica | 0.66025 |
Após a correção de Holm, nenhuma das três comparações apresentou diferença estatisticamente significante ao nível de 5%.
Conclusão
As correlações amostrais apresentaram a ordenação:
\[ r_{\mathrm{Setosa}}> r_{\mathrm{Versicolor}}> r_{\mathrm{Virginica}} \]
Embora a correlação observada em setosa tenha sido maior que em virginica, essa diferença não permaneceu estatisticamente significante após a correção de Holm.
Conclui-se que não foram encontradas evidências estatísticas suficientes para rejeitar as hipóteses de igualdade das correlações populacionais em nenhuma das três comparações pareadas, considerando o controle do erro familiar em 5%.
Entretanto, a ausência de diferenças estatisticamente significantes não demonstra equivalência ou igualdade das correlações populacionais.
Teste Q de Cochran para homogeneidade de correlações independentes
Hipóteses estatísticas
O teste Q de Cochran, implementado por metafor::rma(),
avalia a homogeneidade das correlações populacionais entre o comprimento
e a largura da sépala nas três espécies de Iris.
\[ H_0:\rho_{\mathrm{Setosa}}= \rho_{\mathrm{Versicolor}}= \rho_{\mathrm{Virginica}} \]
\[ H_1:\text{Pelo menos uma correlação populacional difere das demais.} \]
# Teste omnibus de homogeneidade das correlações
teste <- metafor::rma(
measure = "ZCOR",
ri = c(rS, rVc, rVg),
ni = c(nS, nVc, nVg),
method = "EE"
)
print(teste, digits = 5)
Equal-Effects Model (k = 3)
I^2 (total heterogeneity / total variability): 64.55%
H^2 (total variability / sampling variability): 2.82
Test for Heterogeneity:
Q(df = 2) = 5.64158, p-val = 0.05956
Model Results:
estimate se zval pval ci.lb ci.ub
0.67814 0.08422 8.05242 <.00001 0.51308 0.84319 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Q = 5.641584 GL = 2 Valor-p = 0.05955876 Teste global de homogeneidade
O teste apresentou os seguintes resultados:
- Estatística de Cochran: \(Q=5.64158\).
- Graus de liberdade: \(gl=2\).
- Valor-p: \(p=0.05956\).
Adotando-se o nível de significância de 5%, não se rejeita a hipótese nula de homogeneidade das três correlações populacionais.
Entretanto, o valor-p relativamente próximo de 0.05 sugere cautela na interpretação, particularmente porque o teste envolve apenas três grupos independentes.
Medidas de heterogeneidade
Foram obtidas as seguintes medidas:
- \(I^2=64.55\%\)
- \(H^2=2.82\)
O índice \(I^2\) estima a proporção da variabilidade observada das correlações transformadas por Fisher que excede a variabilidade esperada exclusivamente pelo erro amostral.
O valor de \(I^2=64.55\%\) sugere heterogeneidade de magnitude relevante entre as correlações das três espécies.
O índice \(H^2=2.82\) indica que a variabilidade total observada é aproximadamente 2.82 vezes a variabilidade esperada por erro amostral.
Embora essas medidas sugiram heterogeneidade, o teste Q não forneceu evidência estatisticamente significante ao nível de 5%.
Não existe contradição entre esses resultados: o teste Q avalia a evidência estatística contra a homogeneidade, enquanto \(I^2\) e \(H^2\) descrevem a magnitude relativa da heterogeneidade observada.
Além disso, com apenas três grupos, a estimativa de \(I^2\) pode apresentar considerável instabilidade, e o teste Q pode ter baixo poder estatístico.
Estimativa combinada da correlação
Sob o modelo de efeito comum (equal-effects model), a estimativa combinada na escala da transformação \(z\) de Fisher foi:
\[ \hat z_F=0.67814 \]
com erro-padrão de 0.08422 e intervalo de confiança de 95%:
\[ \text{IC}_{95\%}(z_F)=[0.51308;\ 0.84319] \]
Aplicando-se a transformação inversa de Fisher:
\[ r=\tanh(z_F) \]
obtém-se:
\[ \boxed{\hat\rho=0.5903} \]
com intervalo de confiança de 95%:
\[ \text{IC}_{95\%}(\rho)=[0.4723;\ 0.6875] \]
O teste da hipótese de correlação comum igual a zero apresentou:
\[ z=8.05242,\qquad p<0.00001 \]
Portanto, a estimativa combinada da correlação foi positiva e estatisticamente significante.
Esse teste não avalia a homogeneidade das correlações entre as espécies.
Ele verifica se a correlação populacional comum, pressuposta pelo modelo, é igual a zero.
A interpretação dessa estimativa como uma correlação comum depende da adequação da hipótese de homogeneidade, que não foi demonstrada pelo teste Q.
Conclusão
Ao nível de significância de 5%, não foram encontradas evidências estatísticas suficientes para rejeitar a hipótese global de igualdade das correlações populacionais entre as três espécies de Iris (\(Q(2)=5.64158\); \(p=0.05956\)).
Entretanto, o índice \(I^2=64.55\%\) sugere heterogeneidade de magnitude relevante, apesar da ausência de significância estatística no teste global. Esse resultado recomenda cautela, considerando o pequeno número de grupos e a limitada precisão da avaliação da heterogeneidade.
Sob a hipótese de correlação populacional comum, a estimativa combinada foi de 0.5903, com IC95% de [0.4723; 0.6875], indicando associação linear positiva entre comprimento e largura da sépala.
Conclui-se que a hipótese de homogeneidade não foi rejeitada, mas os resultados não permitem afirmar que as três correlações populacionais sejam iguais. A estimativa combinada deve ser interpretada condicionalmente à adequação do modelo de efeito comum.
Homogeneidade das matrizes de correlação entre espécies de Iris
Serão comparadas as matrizes de correlação de Pearson das quatro
medidas morfológicas do conjunto iris:
- Sepal.Length
- Sepal.Width
- Petal.Length
- Petal.Width
Cada matriz possui seis correlações distintas, correspondentes às combinações das quatro variáveis tomadas duas a duas.
Preparação dos dados
Dados <- datasets::iris
vars <- names(Dados)[1:4]
# Separação das três espécies
Grupos <- split(Dados[, vars], Dados$Species)
# Tamanhos amostrais
sapply(Grupos, nrow) setosa versicolor virginica
50 50 50 $setosa
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.000 0.743 0.267 0.278
Sepal.Width 0.743 1.000 0.178 0.233
Petal.Length 0.267 0.178 1.000 0.332
Petal.Width 0.278 0.233 0.332 1.000
$versicolor
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.000 0.526 0.754 0.546
Sepal.Width 0.526 1.000 0.561 0.664
Petal.Length 0.754 0.561 1.000 0.787
Petal.Width 0.546 0.664 0.787 1.000
$virginica
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.000 0.457 0.864 0.281
Sepal.Width 0.457 1.000 0.401 0.538
Petal.Length 0.864 0.401 1.000 0.322
Petal.Width 0.281 0.538 0.322 1.000Comparação entre duas espécies
Serão comparadas as matrizes de correlação das espécies setosa e versicolor.
As hipóteses são:
\[ H_0:\mathbf R_{\mathrm{Setosa}}= \mathbf R_{\mathrm{Versicolor}} \]
\[ H_1:\mathbf R_{\mathrm{Setosa}}\neq \mathbf R_{\mathrm{Versicolor}} \]
Como existem quatro variáveis, a hipótese nula envolve seis restrições independentes.
Teste de Jennrich
O teste de Jennrich utiliza uma estatística com distribuição assintótica qui-quadrado com seis graus de liberdade, sob a hipótese de igualdade das matrizes.
teste.jennrich <- psych::cortest.jennrich(
R1 = R$setosa,
R2 = R$versicolor,
n1 = nrow(Grupos$setosa),
n2 = nrow(Grupos$versicolor)
)
print(teste.jennrich)$chi2
[1] 30.305
$prob
[1] 3.439247e-05O teste de Jennrich identificou diferença estatisticamente significante entre as matrizes de correlação das quatro medidas morfológicas das espécies setosa e versicolor (\(\chi^2(6)=30.305\); \(p<0.001\)).
Portanto, rejeita-se a hipótese de homogeneidade das matrizes de correlação, ao nível de significância de 5%.
Teste ATS com bootstrap
Os testes de Jennrich e ATS avaliam a mesma hipótese nula, mas empregam estatísticas e procedimentos de inferência distintos.
set.seed(123)
teste.2 <- CovCorTest::test_correlation(
X = Grupos[c("setosa", "versicolor")],
hypothesis = "equal-correlated",
method = "BT",
repetitions = 5000
)
print(teste.2)
Correlation Test
2 groups
Hypothesis: equal-correlated
Teststatistic value: 7.7065
p-value: p < 2e-04
p-value computed using Bootstrap with B=5000 repetitions ATS = 7.706507 Valor-p = 0 O teste ATS com bootstrap (\(B=5000\)) identificou diferença estatisticamente significante entre as matrizes de correlação das quatro medidas morfológicas das espécies setosa e versicolor (ATS = 7.7065; \(p<0.0002\)).
Portanto, rejeita-se a hipótese de homogeneidade das matrizes de correlação, ao nível de significância de 5%, corroborando o resultado do teste de Jennrich.
Comparação entre as três espécies
Neste caso, pretende-se avaliar a igualdade global das matrizes de correlação das três espécies.
As hipóteses são:
\[ H_0: \mathbf R_{\mathrm{Setosa}}= \mathbf R_{\mathrm{Versicolor}}= \mathbf R_{\mathrm{Virginica}} \]
\[ H_1: \text{Pelo menos uma matriz de correlação difere das demais.} \]
O número de restrições independentes é:
\[ gl_{\mathrm{restrições}}= (k-1)\frac{p(p-1)}{2} \]
Para \(k=3\) grupos e \(p=4\) variáveis:
\[ gl_{\mathrm{restrições}}= (3-1)\frac{4(4-1)}{2}=12 \]
Teste global ATS com bootstrap
O procedimento testa simultaneamente a igualdade das seis correlações populacionais entre as três espécies, considerando as dependências entre as correlações estimadas dentro de cada grupo.
set.seed(123)
teste.3 <- CovCorTest::test_correlation(
X = Grupos,
hypothesis = "equal-correlated",
method = "BT",
repetitions = 5000
)
print(teste.3)
Correlation Test
3 groups
Hypothesis: equal-correlated
Teststatistic value: 4.9635
p-value: p < 2e-04
p-value computed using Bootstrap with B=5000 repetitions ATS = 4.963533 Valor-p = 0 O teste ATS com bootstrap (\(B=5000\)) identificou diferença estatisticamente significante entre as matrizes de correlação das quatro medidas morfológicas das três espécies de Iris (ATS = 4.9635; \(p<0.0002\)).
Portanto, rejeita-se a hipótese global de homogeneidade das matrizes de correlação, ao nível de significância de 5%, indicando que pelo menos uma das três matrizes populacionais difere das demais.
Comparaçã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 individualizados 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 individualizados (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 individualizados 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 individualizados
Teste de correlação de Pearson para uma condição
alfa <- 0.05
n <- 100
r <- 0.2
# 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.20000 0.00373 0.38143 # Teste t
t <- r * sqrt((n - 2) / (1 - r^2))
p_t <- 2 * pt(-abs(t), df = n - 2)
cat("\nt(", n - 2, ") = ", round(t, 2), ", p = ", p_t, "\n\n", sep = "")
t(98) = 2.02, p = 0.04603629# Teste z de Fisher (aproximação assintótica)
z <- sqrt(n - 3) * atanh(r)
p_z <- 2 * pnorm(-abs(z))
cat("z = ", round(z, 2), ", p = ", p_z, sep = "")z = 2, p = 0.04585951# 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)
zr <- atanh(r)
z0 <- atanh(rho0)
se <- 1 / sqrt(n - 3)
zstat <- (zr - z0) / se
p <- switch(alternative,
"two.sided" = 2 * pnorm(-abs(zstat)),
"greater" = 1 - pnorm(zstat),
"less" = pnorm(zstat))
if (alternative == "two.sided") {
zcrit <- qnorm(1 - (1 - conf.level) / 2)
ci <- tanh(c(zr - zcrit * se, zr + zcrit * se))
} else if (alternative == "greater") {
zcrit <- qnorm(conf.level)
ci <- c(tanh(zr - zcrit * se), 1)
} else {
zcrit <- qnorm(conf.level)
ci <- c(-1, tanh(zr + zcrit * se))
}
list(z = zstat, p = p, conf.int = ci, conf.level = conf.level)
}
# exemplo
test_r_rho0(r = r, n = n, rho0 = 0)$z
[1] 1.996684
$p
[1] 0.04585951
$conf.int
[1] 0.003728345 0.381434008
$conf.level
[1] 0.95Correlaçã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.0173Correlaçõ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.9451Correlaçõ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), main="n=20", 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), main="n=20", 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), main="n=40", 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 arquivo de dados 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.