invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))

suppressMessages(library(readxl))
suppressMessages(library(DescTools))
suppressMessages(library(effectsize))
suppressMessages(library(psych))
suppressMessages(library(car))
suppressMessages(library(bootES))
suppressMessages(library(HH))
suppressMessages(library(shiny))
suppressMessages(library(zeroEQpart))
suppressMessages(library(RVAideMemoire))
suppressMessages(library(MVN))
suppressMessages(library(MASS))
suppressMessages(library(RVAideMemoire))
suppressMessages(library(cocor))
suppressMessages(library(dplyr))

Material

HTML de R Markdown em RPubs
Arquivos em GitHub

Introdução

Esta prática é relativa aos capítulos 6 do livro adotado.

Análise de Correlação

Dados <- readxl::read_excel(path="Biometria_FMUSP.xlsx",
                            sheet="dados",
                            na=c("NA","na","Na","nA"))

# exclusao de dados com erros de digitacao
Dados$MCT[Dados$MCT==658] <- NA
Dados$Estatura[Dados$Estatura==120] <- NA

Dados <- dplyr::mutate_if(Dados, is.character, as.factor)
Dados$AtivFisica <- factor(Dados$AtivFisica,
                           levels=c("sempre_inativo",
                                    "atualmente_inativo",
                                    "baixa_intensidade",
                                    "media_intensidade",
                                    "alta_intensidade"),
                           labels=c("sempre_inativo",
                                    "atualmente_inativo",
                                    "baixa_intensidade",
                                    "media_intensidade",
                                    "alta_intensidade"))
Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")
Dados.M$Estaturaz <- scale(Dados.M$Estatura)[,1]
Dados.M$MCTz <- scale(Dados.M$MCT)[,1]
Dados.F$Estaturaz <- scale(Dados.F$Estatura)[,1]
Dados.F$MCTz <- scale(Dados.F$MCT)[,1]
Dados <- rbind(Dados.M,Dados.F)
DadosSM <- Dados[!is.na(Dados$Estatura),c("Estatura","MCT","Sexo","Sedentarismo")]
DadosSM <- DadosSM[!is.na(DadosSM$MCT),]
alfa <- 0.05

# Delineamento entre participante
print(psych::describeBy(x=Dados[c("Estatura","MCT")],
                        group=Dados$Sexo,
                        mat=1,
                        digits=2))

          item group1 vars   n   mean    sd median trimmed   mad min max range
Estatura1    1      F    1 229 164.34  6.38    163  163.94  5.93 153 186    33
Estatura2    2      M    1 313 175.77  7.26    175  175.88  7.41 155 195    40
MCT1         3      F    2 230  57.63  9.05     56   56.75  7.41  41 120    79
MCT2         4      M    2 312  71.59 12.09     72   70.88 11.12  43 125    82
           skew kurtosis   se
Estatura1  0.61     0.02 0.42
Estatura2 -0.13     0.21 0.41
MCT1       2.43    12.04 0.60
MCT2       0.90     2.37 0.68

car::densityPlot(Estatura~Sexo, data=Dados)

boxplot(Estatura~Sexo, data=Dados)

car::densityPlot(MCT~Sexo, data=Dados)

boxplot(MCT~Sexo, data=Dados)

cat("\nEstimativas pontuais das correlações entre Estatura e MCT por sexo:")


Estimativas pontuais das correlações entre Estatura e MCT por sexo:

corM <- cor(Dados.M$Estatura, Dados.M$MCT, use="complete")
cat("\n\tCorrelação (Masculino): ",corM,sep="")


    Correlação (Masculino): 0.5343547

corF <- cor(Dados.F$Estatura, Dados.F$MCT, use="complete")
cat("\n\tCorrelação (Feminino): ",corF,sep="")


    Correlação (Feminino): 0.4206021

cat("\n\nMatrizes:")



Matrizes:

corM <- cor(Dados.M[c("Estatura","MCT")], use="complete")
cat("\n\tMatriz de correlação (Masculino):\n")


    Matriz de correlação (Masculino):

print(corM)

          Estatura       MCT
Estatura 1.0000000 0.5343547
MCT      0.5343547 1.0000000

varM <- var(Dados.M[c("Estatura","MCT")], use="complete")
cat("\n\tMatriz de covariância (Masculino):\n")


    Matriz de covariância (Masculino):

print(varM)

         Estatura       MCT
Estatura 52.88720  46.97123
MCT      46.97123 146.10124

cat("\nConversão de covariância para correlação:")


Conversão de covariância para correlação:

rM <- varM[1,2]/sqrt(varM[1,1]*varM[2,2])
cat("\n\tCorrelação (Masculino, usando fórmula): ",rM,sep="")


    Correlação (Masculino, usando fórmula): 0.5343547

cat("\n\tMatriz de correlação (Masculino, usando cov2cor):\n")


    Matriz de correlação (Masculino, usando cov2cor):

print(cov2cor(varM))

          Estatura       MCT
Estatura 1.0000000 0.5343547
MCT      0.5343547 1.0000000

corF <- cor(Dados.F[c("Estatura","MCT")], use="complete")
cat("\n\tMatriz de correlação (Feminino):\n")


    Matriz de correlação (Feminino):

print(corF)

          Estatura       MCT
Estatura 1.0000000 0.4206021
MCT      0.4206021 1.0000000

varF <- var(Dados.F[c("Estatura","MCT")], use="complete")
cat("\n\tMatriz de covariância (Feminino):\n")


    Matriz de covariância (Feminino):

print(varF)

         Estatura      MCT
Estatura 40.68172 21.65820
MCT      21.65820 65.17823

cat("\nConversão de covariância para correlação:")


Conversão de covariância para correlação:

rF <- varF[1,2]/sqrt(varF[1,1]*varF[2,2])
cat("\n\tCorrelação (Feminino, usando fórmula): ",rF,sep="")


    Correlação (Feminino, usando fórmula): 0.4206021

cat("\n\tMatriz de correlação (Feminino, usando cov2cor):\n")


    Matriz de correlação (Feminino, usando cov2cor):

print(cov2cor(varM))

          Estatura       MCT
Estatura 1.0000000 0.5343547
MCT      0.5343547 1.0000000

# Bag plot
bgp <- DescTools::PlotBag(Dados.M$Estatura, 
                          Dados.M$MCT,
                          main="Estudante Masculino de Medicina - FMUSP",
                          xlab="Estatura (cm)",
                          ylab="MCT (kg)",
                          na.rm=TRUE,
                          show.bagpoints=FALSE,
                          show.looppoints=FALSE,
                          show.whiskers=FALSE,
                          col.loophull="white",
                          col.looppoints="black", 
                          col.baghull="white",
                          col.bagpoints="black",
                          cex=1)

[1] "Warning: NA elements have been removed!!"

print(outliers <- as.data.frame(bgp$pxy.outlier))

for (o in 1:nrow(outliers))
{
  r <- which(Dados.M$Estatura==outliers$x[o] & 
               Dados.M$MCT==outliers$y[o])
  text(outliers$x[o],outliers$y[o], r, pos=1, cex=0.7)
}

bgp <- DescTools::PlotBag(Dados.F$Estatura, 
                          Dados.F$MCT,
                          main="Estudante Feminino de Medicina - FMUSP",
                          xlab="Estatura (cm)",
                          ylab="MCT (kg)",
                          na.rm=TRUE,
                          show.bagpoints=FALSE,
                          show.looppoints=FALSE,
                          show.whiskers=FALSE,
                          col.loophull="white",
                          col.looppoints="black", 
                          col.baghull="white",
                          col.bagpoints="black",
                          cex=1)

[1] "Warning: NA elements have been removed!!"

print(outliers <- as.data.frame(bgp$pxy.outlier))

for (o in 1:nrow(outliers))
{
  r <- which(Dados.F$Estatura==outliers$x[o] & 
               Dados.F$MCT==outliers$y[o])
  text(outliers$x[o],outliers$y[o], r, pos=1, cex=0.7)
}

sunflowerplot(MCT~Estatura,
              data=Dados.F,
              main="Estudante de Medicina - FMUSP\n(sunflowerplot)",
              xlab="Estatura (cm)",
              ylab="MCT (kg)",
              rotate=TRUE, 
              size=.1,
              col="red", 
              seg.col="red", 
              seg.lwd=.8,
              xlim=c(150,200),
              ylim=c(40,140),
              axes=FALSE)
axis(1)
axis(2)
sunflowerplot(MCT~Estatura,
              data=Dados.M,
              rotate=TRUE, 
              size=.1,
              col="blue", 
              seg.col="blue", 
              seg.lwd=.8,
              add=TRUE)

sunflowerplot(Estaturaz~MCTz,
              data=Dados.F,
              main="Estudante de Medicina - FMUSP\n(sunflowerplot)",
              xlab="Estatura (z)",
              ylab="MCT (z)",
              rotate=TRUE, 
              size=.1,
              col="red", 
              seg.col="red", 
              seg.lwd=.8,
              ylim=c(-2.5,5),
              xlim=c(-3,3.5),
              axes=FALSE)
axis(1, pos=0)
axis(2, pos=0)
sunflowerplot(Estaturaz~MCTz,
              data=Dados.M,
              rotate=TRUE, 
              size=.1,
              col="blue", 
              seg.col="blue", 
              seg.lwd=.8,
              add=TRUE)

car::scatterplot(MCT~Estatura, 
                 main="Estudante de Medicina - FMUSP\nElipses de predição 68%\n(car::scatterplot)",
                 xlab="Estatura (cm)",
                 ylab="MCT (kg)",
                 groups=Dados$Sexo,
                 regLine=FALSE, 
                 smooth=FALSE, 
                 boxplots=TRUE, 
                 ellipse=list(levels=c(0.68), 
                              robust=TRUE, 
                              fill=FALSE, 
                              fill.alpha=0.2),
                 grid=FALSE,
                 col=c("red","blue"),
                 xlim=c(150,200),
                 ylim=c(40,140),
                 data=Dados,
                 axes=FALSE)
axis(1)
axis(2)

car::scatterplot(MCTz~Estaturaz, 
                 main="Estudante de Medicina - FMUSP\nElipses de predição 68%\n(car::scatterplot)",
                 xlab="Estatura (z)",
                 ylab="MCT (z)",
                 groups=Dados$Sexo,
                 regLine=FALSE, 
                 smooth=FALSE, 
                 boxplots=TRUE, 
                 ellipse=list(levels=c(0.68), 
                              robust=TRUE, 
                              fill=FALSE, 
                              fill.alpha=0.2),
                 grid=FALSE,
                 col=c("red","blue"),
                 ylim=c(-2.5,5),
                 xlim=c(-3,3.5),
                 data=Dados,
                 axes=FALSE)
axis(1, pos=0)
axis(2, pos=0)

car::dataEllipse(Dados$Estatura, Dados$MCT, 
                 groups=Dados$Sexo, 
                 group.labels=c("Feminino", "Masculino"),
                 levels=0.68,
                 robust=TRUE,
                 col=c("red","blue"),
                 main="Estudante de Medicina - FMUSP\nElipse de predicao de 68%\n(car::dataEllipse)",
                 xlab="Estatura (cm)", 
                 ylab="MCT (kg)",
                 xlim=c(150, 200),
                 ylim=c(40, 125))

car::dataEllipse(Dados$Estaturaz, Dados$MCTz, 
                 groups=Dados$Sexo, 
                 group.labels=c("Feminino", "Masculino"),
                 levels=0.68,
                 robust=TRUE,
                 col=c("red","blue"),
                 main="Estudante de Medicina - FMUSP\nElipse de predicao de 68%\n(car::dataEllipse)",
                 xlab="Estatura (z)", 
                 ylab="MCT (z)",
                 ylim=c(-2.5,5),
                 xlim=c(-3,3.5))

# Distribuição Binormal Padrão (teoria)
bvn <- HH::bivariateNormal(0.6) 
print(bvn)

update(bvn[3], layout=c(1,1))

# if (interactive())
#   shiny::runApp(file.path(system.file(package="HH"), 
#                           "shiny/bivariateNormal")) 
# # Correlação sob binormal (teoria)
# if (interactive())
#   shiny::runApp(system.file("shiny/bivariateNormalScatterplot", 
#                             package="HH")) 

# Testes de normalidade
result <- MVN::mvn(data=Dados[,c("Estatura","MCT","Sexo")], 
                   subset="Sexo", 
                   mvnTest="hz", 
                   univariateTest="SW")
cat("\nBinormalidade:\n")


Binormalidade:

print(result$multivariateNormality)

$F
           Test       HZ      p value MVN
1 Henze-Zirkler 3.506795 2.798703e-08  NO

$M
           Test       HZ     p value MVN
1 Henze-Zirkler 1.365331 0.007788244  NO

cat("\nUninormalidade:\n")


Uninormalidade:

print(result$univariateNormality)

$F
          Test  Variable Statistic   p value Normality
1 Shapiro-Wilk Estatura     0.9626  <0.001      NO    
2 Shapiro-Wilk    MCT       0.9055  <0.001      NO    

$M
          Test  Variable Statistic   p value Normality
1 Shapiro-Wilk Estatura     0.9907  0.0446      NO    
2 Shapiro-Wilk    MCT       0.9541  <0.001      NO

# Simulacao de distribuicao binormal de estatura e MCT 
# Estudante masculino do curso de 
# Medicina da FMUSP de 2015-18
n <- 5e3
m.estat.m <- 177; m.mct.m <- 73
s.estat.m <- 7; s.mct.m <- 9
r <- 0.6 
cov <- r*s.estat.m*s.mct.m
# Parameters for bivariate normal distribution
mv <- c(m.estat.m, m.mct.m) # Mean vector
covar <- matrix(c(s.estat.m^2, cov, 
                  cov, s.mct.m^2), 2) # Covariance matrix
set.seed(123)
estat.mct.m <- MASS::mvrnorm(n, mu = mv, Sigma = covar)
dados <- as.data.frame(estat.mct.m)
names(dados) <- c("Estaturam", "MCTm")
minx <- min(dados$Estaturam)*0.99
maxx <- max(dados$Estaturam)*1.02 
miny <- min(dados$MCTm)*0.99
maxy <- max(dados$MCTm)*1.02
plot(dados,
     main=paste0("Dados simulados sob binormalidade",
                 ", n(masculino) = ",n,
                 ", rho = ",r),
     xlab="Estatura (cm)", 
     xlim=c(minx,maxx),
     ylab="MCT (kg)",
     ylim=c(miny,maxy))

# car::dataEllipse(estat.mct.m[,1], estat.mct.m[,2], # suspeito...
#                  levels=c(.5, .95, .999),
#                  robust=TRUE, grid=FALSE, 
#                  lwd=2, lty=1, 
#                  col=c("black"),
#                  cex=0.1,
#                  main="Elipse de predição de 50, 95 e 99.9% de binormal",
#                  xlab="Estatura (cm)", 
#                  xlim=c(minx,maxx),
#                  ylab="MCT (kg)",
#                  ylim=c(miny,maxy))
car::scatterplot(estat.mct.m[,1], estat.mct.m[,2],
                 regLine=FALSE,
                 ellipse=list(levels=c(.5, .95, .999), 
                              robust=TRUE, fill=FALSE), 
                 grid=FALSE,
                 smooth=FALSE,
                 col=c("black"),
                 cex=0.1,
                 pch=1,
                 lwd=2, lty=1, 
                 main="Elipse de predição de 50, 95 e 99.9% de binormal",
                 xlab="Estatura (cm)", 
                 xlim=c(minx,maxx),
                 ylab="MCT (kg)",
                 ylim=c(miny,maxy))

# Bag plot
bgp <- DescTools::PlotBag(dados$Estaturam, 
                          dados$MCTm,
                          main="Estudante Masculino de Medicina - FMUSP (simulado)",
                          xlab="Estatura (cm)", 
                          xlim=c(minx,maxx),
                          ylab="MCT (kg)",
                          ylim=c(miny,maxy),
                          na.rm=TRUE,
                          show.bagpoints=FALSE,
                          show.looppoints=FALSE,
                          show.whiskers=FALSE,
                          col.loophull="white",
                          col.looppoints="black", 
                          col.baghull="white",
                          col.bagpoints="black",
                          cex=1)
print(outliers <- as.data.frame(bgp$pxy.outlier))

         x         y
1 163.6222  85.73380
2 184.9955 102.49039
3 198.6472  75.42137
4 181.1799  45.23651
5 157.0260  76.94900
6 149.4119  52.62527
7 183.3122  51.80230
8 194.1720 104.19144

for (o in 1:nrow(outliers))
{
  r <- which(dados$Estaturam==outliers$x[o] & 
               dados$MCTm==outliers$y[o])
  text(outliers$x[o],outliers$y[o], r, pos=1, cex=0.7)
}

# Teste de correlacao de Pearson
cat("\nTeste de correlação de Pearson para uma condição (Masculino):\n")


Teste de correlação de Pearson para uma condição (Masculino):

cat("\n\t- utilizando cor.test() -\n")


    - utilizando cor.test() -

print(cor.test(Dados.M$Estatura, Dados.M$MCT, na.rm=TRUE))


    Pearson's product-moment correlation

data:  Dados.M$Estatura and Dados.M$MCT
t = 11.131, df = 310, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4500174 0.6092445
sample estimates:
      cor 
0.5343547

cat("\n\t- utilizando bootES::bootES() -\n")


    - utilizando bootES::bootES() -

print(bootES::bootES(na.omit(Dados.M[c("Estatura","MCT")]), R=1e4))


95.00% bca Confidence Interval, 10000 replicates
Stat        CI (Low)    CI (High)   bias        SE          
0.534       0.438       0.614       0.000       0.045

cat("\nTeste de correlação para H0: rho = #qualquer valor\n")


Teste de correlação para H0: rho = #qualquer valor

teste_correl_bilateral <- function(r,ro,n)
{
  z <- sqrt((n-3)/2)*log(((1-ro)/(1+ro))/((1-r)/(1+r)))
  p <- 2*pnorm(-abs(z))
  return(c(z,p))
}
ro0 <- 0.6
cat("\n\tCorrelação (hipotetizada): ",ro0,sep="")


    Correlação (hipotetizada): 0.6

corM <- cor(Dados.M$Estatura, Dados.M$MCT, use="complete")
cat("\n\tCorrelação (observada): ",corM,sep="")


    Correlação (observada): 0.5343547

fit <- teste_correl_bilateral(r=corM,
                              ro=ro0,
                             n=min(sum(!is.na(Dados.M$Estatura)),
                                   sum(!is.na(Dados.M$MCT))))
cat("\nz =",fit[1], "p =",fit[2])


z = -2.409559 p = 0.01597182

# Teste de igualdade de correlacoes em condicoes independentes
cat("\n\nTeste de correlação de Pearson para duas condições independentes (Masculino e Feminino):\n")



Teste de correlação de Pearson para duas condições independentes (Masculino e Feminino):

cat("\n\t- utilizando bootES::bootES() -\n")


    - utilizando bootES::bootES() -

print(bootES::bootES(na.omit(Dados[c("Estatura","MCT","Sexo")]),
                     group.col="Sexo", R=1e4))


95.00% bca Confidence Interval, 10000 replicates
Stat        CI (Low)    CI (High)   bias        SE          
-0.114      -0.260      0.032       0.004       0.073

cat("\n\t- utilizando cocor::cocor() -\n")


    - utilizando cocor::cocor() -

# https://f-santos.gitlab.io/2020-04-01-comparing-correlation-coefficients.html
print(cocor::cocor(~ Estatura + MCT | Estatura + MCT, 
                   data=list(data.frame(Dados.M),
                       data.frame(Dados.F))))


  Results of a comparison of two correlations based on independent groups

Comparison between r1.jk (Estatura, MCT) = 0.5344 and r2.hm (Estatura, MCT) = 0.4206
Difference: r1.jk - r2.hm = 0.1138
Data: list(data.frame(Dados.M), data.frame(Dados.F)): j = Estatura, k = MCT, h = Estatura, m = MCT
Group sizes: n1 = 312, n2 = 229
Null hypothesis: r1.jk is equal to r2.hm
Alternative hypothesis: r1.jk is not equal to r2.hm (two-sided)
Alpha: 0.05

fisher1925: Fisher's z (1925)
  z = 1.6886, p-value = 0.0913
  Null hypothesis retained

zou2007: Zou's (2007) confidence interval
  95% confidence interval for r1.jk - r2.hm: -0.0180 0.2492
  Null hypothesis retained (Interval includes 0)

cat("\n\t- implementando função fisher.z() -\n")


    - implementando função fisher.z() -

fisher.z<- function (r1,r2,n1,n2)
{
  return (((0.5*log((1+r1)/(1-r1)))-(0.5*log((1+r2)/(1-r2))))/((1/(n1-3))+(1/(n2-3)))^0.5)
}
cor1 <- cor(Dados.M$Estaturaz, Dados.M$MCTz, use="complete.obs")
cor2 <- cor(Dados.F$Estaturaz, Dados.F$MCTz, use="complete.obs")
z.score <- fisher.z(cor1,cor2,nrow(Dados.M),nrow(Dados.F))
p.value <- 2*(1-pnorm(abs(z.score)))
cat("\n\nComparação das correlações (Fisher's z, versão analítica):")



Comparação das correlações (Fisher's z, versão analítica):

cat("\n\t- Masculino: ", cor1,sep="")


    - Masculino: 0.5343547

cat("\n\t- Feminino: ", cor2,sep="")


    - Feminino: 0.4206021

cat("\nDiferença: ", cor1-cor2,
    ", z = ",z.score,", p = ",p.value,sep="")


Diferença: 0.1137526, z = 1.700783, p = 0.08898375

cat("\n\n\n\t- implementando por bootstrapping percentílico -\n")




    - implementando por bootstrapping percentílico -

# diferenca de correlacoes por bootstrap
diff_corr <- function(data1, data2) 
{
  r.data <- sample(1:nrow(data1),size=nrow(data1),replace=TRUE)  
  data <- data1[r.data,c("Estaturaz","MCTz")]
  cor1 <- cor(data$Estaturaz, data$MCTz, use="complete.obs")
  r.data <- sample(1:nrow(data2),size=nrow(data2),replace=TRUE)  
  data <- data2[r.data,c("Estaturaz","MCTz")]
  cor2 <- cor(data$Estaturaz, data$MCTz, use="complete.obs")
  return(cor1 - cor2)
}
B <- 1e5
res_boot <- replicate(n=B, diff_corr(data1=Dados.M, data2=Dados.F))
d <- density(res_boot)
intconf <- HDInterval::hdi(res_boot,credMass = 1-alfa)
lwr <- as.numeric(intconf[1])
upr <- as.numeric(intconf[2])
cy <- max(d$y)/20
dy <- max(d$y)/40
densidade <- plot(d, 
                  main=paste0("Distribuição das Diferenças (",B," reamostragens)"),
                  xlab="Correlação M - aCorrelação F", ylab="Densidade")
lines(c(lwr,lwr,lwr,upr,upr,upr),
      c(cy+dy,cy-dy,cy,cy,cy+dy,cy-dy))
abline(v=0,lty=2)

cat("\n\nComparação das correlações (",sprintf("%d",B),
    " reamostragens):",sep="")



Comparação das correlações (100000 reamostragens):

cat("\nDiferença: ", mean(res_boot)," IC",(1-alfa)*100,
    "(diferenca pop)=[",lwr,",",upr,"]",sep="")


Diferença: 0.1100015 IC95(diferenca pop)=[-0.03941457,0.2512093]

# Tests for partial association/correlation between paired samples
cat("\n\nCorrelação entre Estatura e MCT (ambos os sexos)\n")



Correlação entre Estatura e MCT (ambos os sexos)

print(cor.test(DadosSM$Estatura,DadosSM$MCT))


    Pearson's product-moment correlation

data:  DadosSM$Estatura and DadosSM$MCT
t = 20.997, df = 539, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6216231 0.7146609
sample estimates:
      cor 
0.6707724

cat("\nCorrelação entre Estatura e MCT (controlando para sexo)\n")


Correlação entre Estatura e MCT (controlando para sexo)

print(RVAideMemoire::pcor.test(Dados$Estatura, 
                               Dados$MCT, 
                               as.numeric(Dados$Sexo),
                               nrep=1e4))


    Pearson's product-moment partial correlation

data:  Dados$Estatura and Dados$MCT, controlling for Dados.Sexo
t = 13.258, df = 538, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4298763 0.5572959
sample estimates:
      cor 
0.4962538

cat("\nTeste de efeito do sexo na correlação\n\n")


Teste de efeito do sexo na correlação

test <- zeroEQpart::pzcor(DadosSM$Estatura, 
                          DadosSM$MCT, 
                          as.numeric(DadosSM$Sexo), 
                          k=1e4)
print(test)

pzcor.default(x = DadosSM$Estatura, y = DadosSM$MCT, z = as.numeric(DadosSM$Sexo), 
    k = 10000)

 * * * * * * * * Test Summary * * * * * * * *                                             
Hypothesis                p.xy - p.xy.z  =  0
Semi_Partial                            FALSE
Correlation                           pearson
p.xy                                  0.67077
p.xy.z                                0.49630
difference                            0.17447
p.value  less than                    0.00031
Bootstrap_Size                          10000


 * * * * * * * * * Distribution Summary * * * * * * * * * * 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1152  0.1605  0.1736  0.1744  0.1873  0.2631

cat("\n\nCorrelação entre Estatura e MCT (controlando para sexo e sedentarismo)\n")



Correlação entre Estatura e MCT (controlando para sexo e sedentarismo)

print(RVAideMemoire::pcor.test(Dados$Estatura, 
                               Dados$MCT, 
                               list(as.numeric(Dados$Sexo),
                                    as.numeric(Dados$Sedentarismo)),
                               nrep=1e4))


    Pearson's product-moment partial correlation

data:  Dados$Estatura and Dados$MCT, controlling for Dados.Sexo, Dados.Sedentarismo
t = 13.003, df = 537, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4223607 0.5510423
sample estimates:
      cor 
0.4893607

cat("\nTeste de efeito do sexo e sedentarismo na correlação\n\n")


Teste de efeito do sexo e sedentarismo na correlação

test <- zeroEQpart::pzcor(DadosSM$Estatura, 
                          DadosSM$MCT, 
                          as.matrix(cbind(
                            as.numeric(DadosSM$Sexo), 
                            as.numeric(DadosSM$Sedentarismo)
                            )),
                          k=1e4)
print(test)

pzcor.default(x = DadosSM$Estatura, y = DadosSM$MCT, z = as.matrix(cbind(as.numeric(DadosSM$Sexo), 
    as.numeric(DadosSM$Sedentarismo))), k = 10000)

 * * * * * * * * Test Summary * * * * * * * *                                             
Hypothesis                p.xy - p.xy.z  =  0
Semi_Partial                            FALSE
Correlation                           pearson
p.xy                                  0.67077
p.xy.z                                0.48942
difference                            0.18135
p.value  less than                    0.00028
Bootstrap_Size                          10000


 * * * * * * * * * Distribution Summary * * * * * * * * * * 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1055  0.1673  0.1802  0.1812  0.1941  0.2653

cat("\n\nCorrelação entre Estatura e MCT (controlando para sedentarismo)\n")



Correlação entre Estatura e MCT (controlando para sedentarismo)

print(RVAideMemoire::pcor.test(Dados$Estatura, 
                               Dados$MCT, 
                               as.numeric(Dados$Sedentarismo),
                               nrep=1e4))


    Pearson's product-moment partial correlation

data:  Dados$Estatura and Dados$MCT, controlling for Dados.Sedentarismo
t = 20.832, df = 538, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6186985 0.7124053
sample estimates:
      cor 
0.6681935

cat("\nTeste de efeito do sedentarismo na correlação\n\n")


Teste de efeito do sedentarismo na correlação

test <- zeroEQpart::pzcor(DadosSM$Estatura, 
                          DadosSM$MCT, 
                          as.numeric(DadosSM$Sedentarismo),
                          k=1e4)
print(test)

pzcor.default(x = DadosSM$Estatura, y = DadosSM$MCT, z = as.numeric(DadosSM$Sedentarismo), 
    k = 10000)

 * * * * * * * * Test Summary * * * * * * * *                                             
Hypothesis                p.xy - p.xy.z  =  0
Semi_Partial                            FALSE
Correlation                           pearson
p.xy                                  0.67077
p.xy.z                                0.66820
difference                            0.00257
p.value  equal to                     0.18910
Bootstrap_Size                          10000


 * * * * * * * * * Distribution Summary * * * * * * * * * * 
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.0152217  0.0006196  0.0021811  0.0025574  0.0041128  0.0228292

Prática R: Análise de Correlação

José Siqueira
Paulo SP Silveira (silveira@usp.br)

09 agosto 2024 20:35h

Material

Introdução

Análise de Correlação

Prática R: Análise de Correlação

José Siqueira Paulo SP Silveira (silveira@usp.br)

09 agosto 2024 20:35h

Material

Introdução

Análise de Correlação

José Siqueira
Paulo SP Silveira (silveira@usp.br)