Análise de Regressão Linear Simples e
Regressão de Deming em R

MSP4112: Estatística Aplicada à Medicina

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

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

suppressMessages(library(bootES, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts = FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts = FALSE))
suppressMessages(library(eirasagree, warn.conflicts = FALSE))
suppressMessages(library(estimatr, warn.conflicts = FALSE))
suppressMessages(library(ggplot2, warn.conflicts = FALSE))
suppressMessages(library(HH, warn.conflicts = FALSE)) # HH::ci.plot
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts = FALSE))
suppressMessages(library(mcr, warn.conflicts = FALSE))
suppressMessages(library(MVN, warn.conflicts = FALSE))
suppressMessages(library(ppcor, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts = FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(RColorBrewer, warn.conflicts=FALSE))
suppressMessages(library(sp, warn.conflicts = FALSE))
suppressMessages(library(WebPower, warn.conflicts = FALSE))

source("BiNormal.R")
source("eiras.bartitle.R")
source("eiras.col2rgbstring.R")
source("eiras.ConfidenceBandBoot.R")
source("eiras.correg.R")
source("eiras.createobj.htest.R")  
source("eiras.deming.regression.R")
source("eiras.density_and_normal.R")
source("eiras.ellipseaxis.R")
source("eiras.friendlycolor.R")
source("eiras.jitter.R")
source("eiras.LambdaEstimate.R")
source("eiras.tukey.try.R")

Adicionar

• bruceR
  • bruceR::cor_diff
  • bruceR::Corr
  • bruceR::Describe
  • bruceR::RECODE
  • bruceR::model_summary
  • bruceR::regress
  • bruceR::EMMEANS
  • bruceR::GLM_summary
  • bruceR::HLM_summary
  • bruceR::TTEST
  • bruceR::MANOVA: https://psychbruce.github.io/bruceR/reference/MANOVA.html
• apaTables::apa.aov.table

apaTables::apa.aov.table(fit,
                         conf.level=1-alfa,
                         type=3)

Pacote eirasagree

Instalar antes os pacotes eirasdata e eiras

• Pacote eirasagree

Funções

• AllStructuralTests Structural Tests, execute all, in sequence
• ConfidenceBandBoot Bootstrapping Confidence Band
• DemingPrepare Computing true values for Deming regression
• DemingRegression Deming Regression and Statistical Test
• DemingRsquared R-squared for Deming Regression
• DescriptiveStat Show data frame content and summary
• EllipticalRegion Analytic Elliptical Confidence Region
• Hello Hello from eiras
• LambdaEstimate Estimation of lambda
• plotBlandAltman Bland Altman plot method
• plotRawData Plot raw data
• RobustConfidenceBand Robust confidence band
• RobustEllipticalRegion Bootstrap Elliptical Confidence Region
• SymbolsAgreement Symbols for Agreement
• testAccuracy Accuracy test
• testDeming Deming Statistical Test
• testPrecision Precision test
• testReliability Reliability test

Material

• HTML de Rmarkdown em RPubs

Objetivos

• Distinguir os conceitos de regressão linear simples e correlação de Pearson.
• Calcular e interpretar os coeficientes de correlação de Pearson e de determinação.
• Definir e construir diagrama de dispersão.
• Calcular e interpretar a reta de regressão (inclinação, intercepto, equação de regressão).
• Avaliar o relacionamento entre a variável dependente/desfecho (VD) intervalar e uma variável explicativa/previsora (VE) intervalar.
• Prever o escore de uma unidade experimental na VD a partir dos valores conhecidos da VE.
• Representar graficamente a relação entre as variáveis de um estudo e a reta de regressão a partir da equação estimada.
• Testar a significância dos coeficientes obtidos na regressão linear simples.
• Determinar o tamanho de efeito da VE.
• Comparar retas de regressão entre condições independentes.

História de regressão & correlação

Galton (1988) introduz o termo co-relation e descreve, com base em medidas antropométricas, o conceito de associação entre variáveis intervalares, especialmente entre pais-filhos humanos.

Galton (1986) introduz o conceito de inclinação de regressão linear simples chamando-o de “ratio of regression” ou “rate of regression”. Ele observou empiricamente e erroneamente a “lei da regressão à mediocridade” (regression toward mediocrity), mas não tinha uma métrica precisa. Sua análise era essencialmente geométrica e descritiva.

• galtondata.R

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)

print(summary(lm(data=Dados_height, child~parent)), digits=3)

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

print(summary(lm(data=data.frame(scale(Dados_height)), child~parent)), digits=4)

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

print(cor(Dados_height), digits=4)

       parent  child
parent 1.0000 0.4588
child  0.4588 1.0000

Dados_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)

print(summary(lm(data=Dados_peas, child~parent)), digits=3)

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

print(summary(lm(data=data.frame(scale(Dados_peas)), child~parent)), digits=4)

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

print(cor(Dados_peas), digits=4)

       parent  child
parent 1.0000 0.3463
child  0.3463 1.0000

Dados_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))

print(summary(lm(data=Dados_cubit, height~cubit)), digits=3)

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

print(summary(lm(data=data.frame(scale(Dados_cubit)), height~cubit)), digits=4)

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

print(cor(Dados_cubit), digits=4)

       height cubit
height  1.000 0.754
cubit   0.754 1.000

Pearson (1896) introduz formalmente o coeficiente de correlação, que mede a associação linear entre duas variáveis intervalares. Ele formalizou o coeficiente de correlação como a covariância padronizada, obtida pela soma dos desvios em relação à média aritmética. Ele resolveu um problema fundamental deixado por Galton: quantificar rigorosamente a relação linear entre duas variáveis biológicas associadas, como esetatura de pais e filhos, massa corporal total e comprimentos das partes do corpo humano etc. Galton

Pearson introduziu o coeficiente \(r\) para resolver isso: um número adimensional entre –1 e +1 que mede a força e direção da associação linear entre duas variáveis intervalares. Assim, ele:

1. Formalizou a ideia de regressão linear como \(\hat{y} = a + b\,x\), onde \(b = r \dfrac{s_y}{s_x}\);
2. Definiu \(r\) como a covariância normalizada, permitindo comparar relações entre variáveis de diferentes escalas;
3. Eliminou o problema conceitual da assimetria galtoniana (em que a “regressão dos filhos nos pais” não era igual à “dos pais nos filhos”), mostrando que \(r_{xy} = r_{yx}\).

Pearson resolveu o problema de como medir, de forma objetiva e simétrica, o grau de dependência linear entre duas características biológicas. Este foi um passo decisivo na transição da biometria descritiva para a estatística inferencial.

Resumindo:

• Coeficiente de correlação é inclinação de regressão padronizado;

• Inclinação de regressão é a correlação escalada pela razão das variâncias.

Regressão linear simples

A análise de regressão linear simples (RLS) é uma extensão da análise de correlação de Pearson.

Ao contrário da correlação, a RLS é direcional e dimensional.

Queremos estabelecer uma equação (uma reta, no caso da RLS) que sirva para, a partir de uma variável explicativa (VE) conhecida, estimar o valor médio de uma variável de desfecho (VD). A equação de regressão preserva as respectivas unidades de medida da VD e da VE.

A pergunta é direcional:

Quanto varia a VD, em média, se a VE variar?

Consideramos, portanto, que RLS pode ser utilizada para estimar a massa corpórea a partir da estatura de animal humano adulto. Assim, ajustamos uma reta (demo_regEstMCT.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(summary(Dados), digits=1)
print(cor.test(~ massa + estatura,
               data=Dados), 
      digits=2)
sunflowerplot(massa ~ estatura,
     xlim = c(163,182), ylim = c(55,90),
     xlab="Estatura (cm)", ylab="Massa (kg)",
     pch=21, col="black", bg="blue",
     data=Dados)
lines(lowess(estatura,massa), lty=2)
modelo <- lm(massa ~ estatura)
sumario <- summary(modelo)
print(sumario, digits=3)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
massa.medio <- intercepto + inclinacao*estatura
lines(estatura, massa.medio, col="black", lwd=3, lty=1)
title(main = sprintf("Massa média = %.2f + %.2f × Estatura", intercepto, inclinacao),
      sub  = bquote("Estatura" %in% "[" * .(min(estatura)) * "," * .(max(estatura)) * "]"))

    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  

    Pearson's product-moment correlation

data:  massa and estatura
t = 3, df = 7, p-value = 0.02
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.14 0.94
sample estimates:
 cor 
0.73 

Call:
lm(formula = massa ~ estatura)

Residuals:
   Min     1Q Median     3Q    Max 
 -7.75  -2.24   1.46   2.51   8.00 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -144.53      76.13   -1.90    0.099 .
estatura        1.26       0.44    2.86    0.024 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.71 on 7 degrees of freedom
Multiple R-squared:  0.538, Adjusted R-squared:  0.472 
F-statistic: 8.16 on 1 and 7 DF,  p-value: 0.0245

Neste gráfico:

• A variável em \(x\) é a VE (estatura (cm)); variável explicativa, preditora, independente.
• A variável em \(y\) é a VD (massa corporal (kg)); variável de desfecho, dependente.

A reta foi estimada a partir da função lm, a qual é parte do pacote nativo stats. Computa o modelo linear da massa (eixo \(y\), VD) em função da estatura (eixo \(x\), VE). Para traçar a reta, foram estimados os coeficientes (intercepto e inclinação) que utilizamos para calcular os valores médios da VD, usando a equação da reta. Esta equação corresponde, no código, à linha

massa.medio <- intercepto + inclinacao*estatura

A variável estatura é um vetor numérico com os valores observados de estatura. Os estimadores de intercepto (\(\hat{\beta}_0\)) e inclinacao (\(\hat{\beta}_1\)) são ótimos e únicos. Consequentemente, esta linha de código computa massa.medio (\(\hat{y}\)) como um vetor numérico de mesmo comprimento de estatura. Com estes valores, foi traçada a reta de regressão.

O intercepto fica em modelo$coefficients[1]=-144.53 (\(\hat{\beta}_0\)) e a inclinação em modelo$coefficients[2]=1.26 (\(\hat{\beta}_1\)). A equação de uma reta é

\[ \hat{y} = \hat{\beta}_0 + \hat{\beta}_1 \,x \]

portanto a linha sólida do gráfico é :

\(\hat{y}=\) -144.53 \(+\) 1.26 \(x\)

\(\min\{x_i\} \le x \le \max\{x_i\}\)

ou, particularizando este exemplo:

\(\widehat{massa}=\) -144.53 \(+\) 1.26 \(\text{estatura}\)

\(165 \le \text{estatura} \le 180\)

Note o circunflexo sobre o valor \(\hat{y}\) (ou \(\widehat{massa}\)), indica que este é o valor médio estimado da VD condicionado ao valor da VE (\(x\) ou \(\text{estatura}\)).

Sinonímia Dependendo do contexto da análise encontramos: \(x\) é a variável explicativa (VE) quantitativa, sem erro de mensuração, também chamada de variável independente (VI), variável preditiva, variável regressora, variável previsora, variável de exposição, covariável ou variável auxiliar. \(y\) é a variável de desfecho (VD) quantitativa, com erro de mensuração, também conhecida como variável dependente, variável-alvo (target), critério ou variável de resposta.

A função que desenvolvemos, internamente, automatiza este procedimento e ainda adiciona a banda de confiança de 95% (que é larga porque a amostra, neste exemplo, é pequena) (demo_corregEstMCT.R):

source("eiras.friendlycolor.R")
source("eiras.correg.R")
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
lst <- correg(estatura, massa, method="lm", 
              xlab="Estatura (cm)", ylab="Massa (kg)",
              pch=21, bg="white")
centro_estatura <- mean(estatura)
centro_massa <- mean(massa)
massa.H0 <- rep(centro_massa,length(estatura))
lines(estatura, massa.H0, col="blue", 
      lwd=2, lty=1)
points(centro_estatura, centro_massa, 
       pch=21, col="black", bg="black")
text(centro_estatura-6, centro_massa + 1.5, expression(H[0]), col="blue", cex=1.2, pos=2)

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

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = Massa (kg) ~ Estatura (cm))"

Residuals:
   Min     1Q Median     3Q    Max 
-7.749 -2.237  1.459  2.507  7.995 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)  
(Intercept)   -144.5250    76.1288  -1.898   0.0994 .
Estatura (cm)    1.2559     0.4396   2.857   0.0245 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.713 on 7 degrees of freedom
Multiple R-squared:  0.5383,    Adjusted R-squared:  0.4723 
F-statistic: 8.161 on 1 and 7 DF,  p-value: 0.02445


Equation:
  mean[Massa (kg)]=-144.525 + 1.25592105 . Estatura (cm)

Significado da reta de regressão

A reta de regressão estimada é o conjunto das médias de VD condicionadas a todos os valores reais da VE no seu intervalo observado. A reta de regressão passa sempre pelo centróide (médias amostrais da VD e VE). A reta de regressão estimada é válida apenas dentro do intervalo dos valores mínimo e máximo observados da VE.

\[ \widehat{\text{VD}} = \hat{\alpha}+\hat{\beta} \; \text{VE} \]

sendo que:

\[ \text{VE} \in \left[\min\{\text{VE}_{i}\}, \max\{\text{VE}_{i}\}\right]\\ i=1,2,\ldots,n \]

Portanto, a reta de regressão é uma interpolação média da VD no domínio observado da VE que passa sempre pelo centróide.

Testes estatísticos

• As suposições principais para lm são:
  • independência entre os pares de observações.
  • ausência de erro de mensuração da VE.
  • distribuição normal da VD.
  • homocedasticidade da VD nos níveis da VE.
  • linearidade da relação entre VD e VE.

Repare que existem estatísticas t e seus valores p calculados nas linhas dos coeficientes. O valor encontrado em sumario$coefficients[2,4]=0.0244505 está associado com a inclinação da reta. O intercepto não será interpretado por enquanto. O teste foi sobre a reta estar suficientemente inclinada para ser utilizada para predizer a massa corporal a partir da estatura:

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

A hipótese nula indica que, caso não houvesse relação entre estatura e massa corpórea, esperaríamos uma horizontal: a massa corpórea oscilaria ao redor de seu valor médio, insensível à variação do valor da estatura. Adicionamos, então, uma linha horizontal na altura da média das massas corporais observadas para representar \(H_0\) (demo_regressao.R):

# demo_regressao.R

source("eiras.friendlycolor.R")

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
sunflowerplot(estatura, massa,
     xlim = c(163,182), ylim = c(55,95),
     xlab="Estatura (cm)", ylab="Massa (kg)",
     pch=21, col="black", bg="white")
modelo <- lm(massa ~ estatura)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
massa.medio <- intercepto + inclinacao*estatura
lines(estatura, massa.medio, col=friendlycolor(8), 
      lwd=3, lty=1)
title(main = sprintf("Massa média = %.2f + %.2f × Estatura", intercepto, inclinacao),
      sub  = bquote("Estatura" %in% "[" * .(min(estatura)) * "," * .(max(estatura)) * "]"))


# media de estatura e media de massa
centro_estatura <- mean(estatura)
centro_massa <- mean(massa)
# media de massa, traca linha horizontal
massa.H0 <- rep(centro_massa,length(estatura))
lines(estatura, massa.H0, col=friendlycolor(30), 
      lwd=3, lty=1)
# centroide, plota o ponto
points(centro_estatura, centro_massa, pch=21, col="black", bg="black")
text(centro_estatura-6, centro_massa + 1.5, expression(H[0]), col=friendlycolor(30), cex=1.2, pos=2)

# parte explicada pelo modelo
for(i in 1:length(estatura))
{
  lines(c(estatura[i]-0.05,estatura[i]-0.05), c(massa.medio[i],massa.H0[i]), 
        lwd=2, lty=3, col=friendlycolor(8))
}

# resíduo
for(i in 1:length(estatura))
{
  lines(c(estatura[i]+0.05,estatura[i]+0.05), c(massa[i],massa.medio[i]), 
        lwd=2, lty=4, col=friendlycolor(25))
}

legend("topleft",
        c("ponto observado",
          "média de massa",
          "reta de regressão",
          "parte explicada",
          "resíduo"
        ), 
        lty=c(NA, 1, 1, 3, 4), 
        lwd=c(NA, 3, 3, 2, 2), 
        pch=c(21, NA, NA, NA, NA),
        col=c("black",friendlycolor(30), friendlycolor(8), 
              friendlycolor(8), friendlycolor(25)), 
        pt.bg=c(friendlycolor(24),friendlycolor(30), friendlycolor(8), 
                friendlycolor(8), friendlycolor(25)), 
        cex=0.9, bty="n")

sumario <- summary(modelo)
print(sumario, digits=4)

Call:
lm(formula = massa ~ estatura)

Residuals:
   Min     1Q Median     3Q    Max 
-7.749 -2.237  1.459  2.507  7.995 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -144.5250    76.1288  -1.898   0.0994 .
estatura       1.2559     0.4396   2.857   0.0245 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.713 on 7 degrees of freedom
Multiple R-squared:  0.5383,    Adjusted R-squared:  0.4723 
F-statistic: 8.161 on 1 and 7 DF,  p-value: 0.02445

Como vimos acima, a hipótese nula é rejeitada para \(\alpha=0.05\), pois \(p=0.0245\) e, portanto, admitimos que a reta tem inclinação diferente de zero. Em outras palavras, existe modelo linear que nos permite, com a equação estimada pela RLS, calcular o valor médio da massa corporal a partir da estatura.

Em RLS, testar \(H_0:\ \beta_1=0\) pode ser feito pela estatística \[ F=(n-2)\frac{R^2}{1-R^2}\ \sim\ F_{1,n-2} \] sendo que \(R^2\) é o coeficiente de determinação e \(n\) o tamanho amostral. Esse teste é equivalentemente expresso pelo teste \(t\) do coeficiente angular, com a identidade exata \(F=t^2\).

# Exemplo com dados de Galton (cúbito e estatura em cm)
Dados <- psychTools::heights
Dados$estatura <- Dados$height * 2.54
Dados$cubito   <- Dados$cubit  * 2.54

fit <- lm(estatura ~ cubito, data = Dados)

sm  <- summary(fit)
n   <- nobs(fit)
R2  <- sm$r.squared
t1  <- sm$coefficients["cubito","t value"]
F_a <- anova(fit)[1,"F value"]
F_r <- (n - 2) * R2 / (1 - R2)

cat(sprintf("t^2 = %.6f | F_ANOVA = %.6f | F_R2 = %.6f\n", t1^2, F_a, F_r))

t^2 = 455.796184 | F_ANOVA = 455.796184 | F_R2 = 455.796184

Os testes estatísticos para correlação e RLS produzem resultados coincidentes. Veja o valor p computado para a correlação com cor.test associado à uma estatística t com \(n-2\) graus de liberdade, e o valor p computado por lm para a inclinação da reta, igualmente associado à esta estatística t (na RLS, também coincidente com o valor p da estatística \(F\) com \(n-2\) graus de liberdade no denominador).

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

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

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

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 

Isto é dizer que, para RLS, é equivalente testar a correlação (\(H_0: \rho=0\)) ou a inclinação da reta de regressão (\(H_0: \beta_1=0\)): variáveis não correlacionadas (\(\rho=0\)) produzem reta de regressão horizontal (\(\beta_1=0\)). No caso da RLS (em que existe somente uma VE), a existência do modelo está diretamente ligada à não nulidade da inclinação da reta populacional, i.e., se a inclinação desta reta é nula, a VD é insensível à variação da VE.

Incluímos uma banda de confiança de 95% para a regressão. Esta banda de confiança foi estimada por bootstrapping. (demo_regressao_BC_boot.R):

# demo_regressao_BC_boot.R

source("eiras.correg.R")
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
lst <- correg(estatura, massa, 
              method="lm", lowess=FALSE,
              B=1e4,
              xlab="Estatura (cm)", ylab="Massa (kg)",
              pch=21, bg="white")
centro_estatura <- mean(estatura)
centro_massa <- mean(massa)
massa.H0 <- rep(centro_massa,length(estatura))
lines(estatura, massa.H0, col="blue", 
      lwd=1, lty=1)
points(centro_estatura, centro_massa, 
       pch=21, col="black", bg="black")
text(centro_estatura-6, centro_massa + 1.5, expression(H[0]), col="blue", cex=1.2, pos=2)

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

------------------------
Simple linear regression
------------------------

    SLR: 10000 resamples

data:  Massa (kg) ~ Estatura (cm)

alternative hypothesis: true slope is not equal to 0

               median         min       1st.Q       3rd.Q        max
intercept -145.544715 -669.704918 -174.852882 -108.722786 370.828283
slope        1.263636   -1.757576    1.049891    1.432402   4.311475
                IC95.1    IC95.2
intercept -239.8415466 15.611511
slope        0.3277106  1.810661

Equation:
  mean[Massa (kg)] = -145.54471545 [-239.84154656, 15.61151148] + 
        1.26363636 [0.32771064, 1.81066078] . Estatura (cm)

Dentro desta banda, com confiança (probabilidade) de 95%, está a reta de regressão populacional. Coerente com os testes estatísticos anteriores, uma outra forma de testar a hopótese nula de que não há uma reta horizontal que possa ser traçada dentro banda de confiança no dominío observado da VE. Portanto, a reta populacional de onde esta amostra é oriunda, tem inclinação diferente de zero.

Coeficiente de determinação

Quando imprimimos o resultado de lm, observamos uma nova estatística (R-squared), $R^2=$0.54, conhecida como coeficiente de determinação (desconsidere o adjetivo Multiple que aparece aqui, pois lm serve, também, para modelar regressão linear múltipla, usada quando há mais de uma VE para predizer a VD). Observe que, na RLS, seu valor é igual ao coeficiente r elevado ao quadrado e armazenado na variável correlacao$estimate=0.73, calculado anteriormente:

0.54 \(= R^2=r^2=\) correlacao$estimate2 \(=\) 0.73² \(=\) 0.54

O valor de \(R^2\) varia de 0 a 1. Neste exemplo, o coeficiente de determinação indica que aproximadamente 53.83% da variância amostral da VD é explicada pela VE. Os restantes 46.17% da variancia amostral da VD são explicados pelo termo de erro, i.e., por outras variáveis ou por outro formato de função (e.g., uma regressão não linear que se ajuste melhor aos dados). Assim como o \(r^2\), \(R^2\) também é medida de tamanho de efeito, equivalente ao \(\eta^2\) omnibus, pois RLS é um caso particular de modelo linear geral (GLM).

Observe, também, o valor p associado ao R-squared por uma estatística F omnibus. Seu valor é igual ao valor p computado para a inclinação da reta e para a correlação de Pearson.

Ellis, 2010

R2 ajustado e ω2 Não confunda o coeficiente de determinação, \(R^2\), com o \(R^2_{\text{ajustado}}\). \(R^2_{\text{ajustado}}\) leva em conta o número de preditores e o tamanho da amostra, penalizando modelos complexos (mais preditores). Portanto, \(R^2_{\text{ajustado}}\) não é uma estatística de tamanho de efeito. Ele é útil na comparação de modelos ajustados aos mesmos dados, pois pode diminuir quando variáveis irrelevantes são incluídas. Entre modelos comparáveis, quanto maior o \(R^2_{\text{ajustado}}\), melhor o equilíbrio entre ajuste e parcimônia. O coeficiente \(\omega^2\) não é uma medida de tamanho de efeito, pois depende do tamanho da amostra. Assim como o \(R^2_{\text{ajustado}}\), o \(\omega^2\) atua como critério de seleção de modelo, e não como estimador da proporção de variância explicada na população. Ambos penalizam a complexidade do modelo (mais preditores), favorecendo modelos mais parcimoniosos, e por isso podem diminuir quando variáveis irrelevantes são incluídas.

Ajuste de reta

O centróide \((\bar{x},\bar{y})\) é o ponto localizado no cruzamento da reta ajustada com a reta horizontal \(H_0\).

Infinitas retas passam pelo centróide. Construindo quadrados cujos lados sejam a distância entre os valores observados \(y\) e valores os estimados \(\hat{y}\) por uma reta, a que melhor se ajusta é aquela cuja somatória das áreas dos quadrados é mínima. Esta é a reta de regressão ajustada pelo método de mínimos quadrados ordinários (MQO ou OLS) representada nos gráficos. Toda reta de regressão ajustada por MQO passa obrigatoriamente pelo centróide.

Execute demo_minimos_quadrados.R para uma demonstração:

source("demo_minimos_quadrados.R")

Resíduo

Também adicionamos ao último gráfico de estatura e massa corporal algumas linhas verticais.

As linhas tracejadas estão:

• em azul e com tracejado mais fino, entre a reta horizontal na altura do valor médio da massa corporal (\(\bar{y}\)) e os valores estimados pela reta de regressão (\(\hat{y}_i\));
• em preto e com tracejado mais grosseiro, entre os valores estimados pela reta de regressão (\(\hat{y}_i\)) e os valores observados (\(y_i\)).

Como a reta de regressão é inclinada, diferencia da horizontal \(H_0\) e, portanto, com a equação encontrada, dada a VE, podemos calcular a VD. No entanto, observando os valores da amostra (pontos no gráfico), vemos que nem tudo foi capturado pela RLS.

Denominamos, aqui, de “parte explicada” pelo modelo, as distâncias verticais entre o valor médio calculado pela reta de regressão e a horizontal (\(\hat{y}_i - \bar{y}\)). Considere que o tanto de inclinação da reta de regressão em relação à horizontal corresponde a quanto da associação dos pontos observados foram capturados pela RLS. Portanto, as partes “não explicadas” pelo modelo são as que sobram, os resíduos, correspondendo às distâncias verticais entre os pontos observados e os estimados pela reta de regressão (\(y_i - \hat{y}_i\)).

Intercepto com significado

Até aqui desprezamos o intercepto igual a -144.53 quando fizemos a regressão linear. Este valor parece não fazer sentido, pelo motivo de que, neste contexto, não faz sentido mesmo.

No entanto, a mesma regressão linear pode ser feita centrando a VE (i.e., subtraindo a média da estatura de todos os valores de estatura observados):

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
media <- mean(estatura)
estatura_centrada <- estatura - media
cat(estatura_centrada)

-8.111111 -4.111111 -2.111111 -1.111111 -0.1111111 0.8888889 2.888889 4.888889 6.888889

mean(estatura_centrada)

[1] -3.157871e-15

A regressão com a estatura centrada está implementada em demo_regressao_centrada.R, resultando em:

# demo_regressao.R

source("eiras.friendlycolor.R")

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
estatura <- estatura - mean(estatura)
sunflowerplot(estatura, massa,
     xlim = c(-10,10), ylim = c(55,95),
     xlab="Estatura (cm)", ylab="Massa (kg)",
     pch=21, col="black", bg="white")
modelo <- lm(massa ~ estatura)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
massa.medio <- intercepto + inclinacao*estatura
lines(estatura, massa.medio, col=friendlycolor(8), 
      lwd=3, lty=1)
title(main = sprintf("Massa média = %.2f + %.2f × EstaturaC", intercepto, inclinacao),
      sub  = bquote("EstaturaC" %in% "[" * .(min(estatura)) * "," * .(max(estatura)) * "]"))


# media de estatura e media de massa
centro_estatura <- mean(estatura)
centro_massa <- mean(massa)
# media de massa, traca linha horizontal
massa.H0 <- rep(centro_massa,length(estatura))
lines(estatura, massa.H0, col=friendlycolor(30), 
      lwd=3, lty=1)
# centroide, plota o ponto
points(centro_estatura, centro_massa, pch=21, col="black", bg="black")
text(centro_estatura-6, centro_massa + 1.5, expression(H[0]), col=friendlycolor(30), cex=1.2, pos=2)

# parte explicada pelo modelo
for(i in 1:length(estatura))
{
  lines(c(estatura[i]-0.05,estatura[i]-0.05), c(massa.medio[i],massa.H0[i]), 
        lwd=2, lty=3, col=friendlycolor(8))
}

# resíduo
for(i in 1:length(estatura))
{
  lines(c(estatura[i]+0.05,estatura[i]+0.05), c(massa[i],massa.medio[i]), 
        lwd=2, lty=4, col=friendlycolor(25))
}

legend("topleft",
        c("ponto observado",
          "média de massa",
          "reta de regressão",
          "parte explicada",
          "resíduo"
        ), 
        lty=c(NA, 1, 1, 3, 4), 
        lwd=c(NA, 3, 3, 2, 2), 
        pch=c(21, NA, NA, NA, NA),
        col=c("black",friendlycolor(30), friendlycolor(8), 
              friendlycolor(8), friendlycolor(25)), 
        pt.bg=c(friendlycolor(24),friendlycolor(30), friendlycolor(8), 
                friendlycolor(8), friendlycolor(25)), 
        cex=0.9, bty="n")

sumario <- summary(modelo)
print(sumario, digits=4)

Call:
lm(formula = massa ~ estatura)

Residuals:
   Min     1Q Median     3Q    Max 
-7.749 -2.237  1.459  2.507  7.995 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  72.8889     1.9044  38.273 2.16e-09 ***
estatura      1.2559     0.4396   2.857   0.0245 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.713 on 7 degrees of freedom
Multiple R-squared:  0.5383,    Adjusted R-squared:  0.4723 
F-statistic: 8.161 on 1 and 7 DF,  p-value: 0.02445

Observe que os valores da regressão foram preservados (coeficiente angular, resíduos, \(R^2\) e estatísticas), exceto pelo intercepto que, agora, corresponde ao valor médio da massa corporal.

Padronização de variável

Outra observação interessante decorre de padronizarmos VD e VE variáveis (i.e., subtrair a média e dividir pelo desvio-padrão). A regressão com as variáveis padronizadas está implementada em demo_regressao_z.R, resultando em:

# demo_regressao.R

source("eiras.friendlycolor.R")

estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
estatura <- scale(estatura)
massa <- scale(massa)
cat("z_estatura:",round(estatura,2),"\n")
cat("z_massa:",round(massa,2),"\n")
sunflowerplot(estatura, massa,
     xlim = c(-2,2), ylim = c(-2,2),
     xlab="zEstatura", ylab="zMassa",
     pch=21, col="black", bg="white")
modelo <- lm(massa ~ estatura)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
massa.medio <- intercepto + inclinacao*estatura
lines(estatura, massa.medio, col=friendlycolor(8), 
      lwd=3, lty=1)
title(main = sprintf("zMassa média = %.2f + %.2f × zEstatura", intercepto, inclinacao),
      sub  = bquote("zEstatura" %in% "[" * .(min(estatura)) * "," * .(max(estatura)) * "]"))

# media de estatura e media de massa
centro_estatura <- mean(estatura)
centro_massa <- mean(massa)
# media de massa, traca linha horizontal
massa.H0 <- rep(centro_massa,length(estatura))
lines(estatura, massa.H0, col=friendlycolor(30), 
      lwd=3, lty=1)
# centroide, plota o ponto
points(centro_estatura, centro_massa, pch=21, col="black", bg="black")
text(centro_estatura-1.2, centro_massa + 0.3, expression(H[0]), col=friendlycolor(30), cex=1.2, pos=2)

# parte explicada pelo modelo
for(i in 1:length(estatura))
{
  lines(c(estatura[i]-0.05,estatura[i]-0.05), c(massa.medio[i],massa.H0[i]), 
        lwd=2, lty=3, col=friendlycolor(8))
}

# resíduo
for(i in 1:length(estatura))
{
  lines(c(estatura[i]+0.05,estatura[i]+0.05), c(massa[i],massa.medio[i]), 
        lwd=2, lty=4, col=friendlycolor(25))
}

legend("topleft",
        c("zponto observado",
          "média de zmassa",
          "reta de regressão padronizada",
          "parte explicada",
          "zresíduo"
        ), 
        lty=c(NA, 1, 1, 3, 4), 
        lwd=c(NA, 3, 3, 2, 2), 
        pch=c(21, NA, NA, NA, NA),
        col=c("black",friendlycolor(30), friendlycolor(8), 
              friendlycolor(8), friendlycolor(25)), 
        pt.bg=c(friendlycolor(24),friendlycolor(30), friendlycolor(8), 
                friendlycolor(8), friendlycolor(25)), 
        cex=0.9, bty="n")

sumario <- summary(modelo)
print(sumario, digits=4)

z_estatura: -1.77 -0.89 -0.46 -0.24 -0.02 0.19 0.63 1.06 1.5 
z_massa: -1.51 0.01 -0.62 0.14 -1 1.16 -0.49 1.03 1.29 

Call:
lm(formula = massa ~ estatura)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.9853 -0.2845  0.1855  0.3187  1.0165 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.061e-15  2.421e-01   0.000   1.0000  
estatura    7.337e-01  2.568e-01   2.857   0.0245 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7264 on 7 degrees of freedom
Multiple R-squared:  0.5383,    Adjusted R-squared:  0.4723 
F-statistic: 8.161 on 1 and 7 DF,  p-value: 0.02445

Com a padronização, não é surpresa que o centróide seja \((0,0)\) (as médias de \(x\) e \(y\) são nulas). Em RLS, a inclinação da reta é estimada por:

\[ \hat{\beta}_1 = r \;\dfrac{s_y}{s_x} \]

e o intercepto é estimado por:

\[ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \;\bar{x} \]

Com a padronização, os desvios-padrão são unitários e as médias são nulas. Observe os valores obtidos pela regressão:

\[\hat{\beta}_1 = r\]

A inclinação da reta é a correlação de Pearson.

\[\hat{\beta}_0 = 0\] O intercepto é nulo.

O valor de \(R^2\) continua inalterado. Concluímos que, em RLS, o coeficiente de determinação é a variância compartilhada das variáveis padronizadas, uma medida de tamanho de efeito que indica quanto uma reta ajustada explica da associação entre duas variáveis intervalares.

Planejamento amostral

O planejamento amostral pode ser feito com a função WebPower::wp.regression. E.g., adotando tamanho de efeito médio de acordo com os dados compilados por Ellis (2010), poder de 80% e nível de significância de 5%, obtemos

R2 <- 0.13
f2 <- R2/(1-R2) # = 0.15 
cat("f^2 = ",round(f2,2),"\n",sep="") 
wp <- WebPower::wp.regression(power=0.8, 
                              alpha=0.05, 
                              p1=1, 
                              f2=f2)
print(wp)

f^2 = 0.15

Power for multiple regression

           n p1 p2        f2 alpha power
    54.51598  1  0 0.1494253  0.05   0.8

URL: http://psychstat.org/regression

Para poder de 90%, obtemos

R2 <- 0.13
f2 <- R2/(1-R2) # = 0.15 
cat("f^2 = ",round(f2,2),"\n",sep="") 
wp <- WebPower::wp.regression(power=0.9, 
                              alpha=0.05, 
                              p1=1, 
                              f2=f2)
print(wp)

f^2 = 0.15

Power for multiple regression

           n p1 p2        f2 alpha power
    72.29472  1  0 0.1494253  0.05   0.9

URL: http://psychstat.org/regression

Esta função serve para planejamento de diversos modelos de regressão linear múltipla. O parâmetro p1 é o número de VE, que no caso da RLS é de apenas 1. O parâmetro f2 é \(f^2\) de Cohen, calculado a partir do \(R^2\). A implementação desta função está de acordo com Cohen (1992).

Cuidados com RLS

Modelo linear é adequado?

RLS supõe que o modelo linear é adequado. Portanto, o primeiro passo a fazer é a estatística descritiva, tanto numérica quanto gráfica. O segundo é verificar se as duas variáveis estão correlacionadas.

Seus dados parecem seguir uma reta?

Um exemplo que ilustra a necessidade da estatística descritiva é o quarteto de Anscombe.

Quarteto de Anscombe é o nome dado a quatro conjuntos de dados que têm estatísticas descritivas quase idênticas (como a média e a variância), mas que têm distribuições muito diferentes e aparências muito distintas quando exibidos graficamente. Cada conjunto de dados consiste de onze pontos (x,y). Eles foram construídos em 1973 pelo estatístico Francis Anscombe, com o objetivo de demonstrar tanto a importância de se visualizar os dados antes de analisá-los, quanto o efeito dos outliers e outras observações influentes nas propriedades estatísticas. Ele descreveu o artigo como tendo a finalidade de combater a impressão entre os estatísticos de que “cálculos numéricos são exatos, mas gráficos são aproximados e grosseiros.”

Os quatro conjuntos (implementados em demo_Anscombe.R) resultam em:

   x1 x2 x3 x4    y1   y2    y3    y4
1  10 10 10  8  8.04 9.14  7.46  6.58
2   8  8  8  8  6.95 8.14  6.77  5.76
3  13 13 13  8  7.58 8.74 12.74  7.71
4   9  9  9  8  8.81 8.77  7.11  8.84
5  11 11 11  8  8.33 9.26  7.81  8.47
6  14 14 14  8  9.96 8.10  8.84  7.04
7   6  6  6  8  7.24 6.13  6.08  5.25
8   4  4  4 19  4.26 3.10  5.39 12.50
9  12 12 12  8 10.84 9.13  8.15  5.56
10  7  7  7  8  4.82 7.26  6.42  7.91
11  5  5  5  8  5.68 4.74  5.73  6.89

  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    9 3.3166   10       8     13      11   14 11    0
2        y 8.81 2.0316 8.04    6.95   7.58    8.33 9.96 11    0
lambda =  1.80086371549847 
Assuming:
    - Explanatory variable: x;
    - Dependent variable: y.

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = y ~ x)"

Residuals:
     Min       1Q   Median       3Q      Max 
-1.92127 -0.45577 -0.04136  0.70941  1.83882 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0001     1.1247   2.667  0.02573 * 
x             0.5001     0.1179   4.241  0.00217 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared:  0.6665,    Adjusted R-squared:  0.6295 
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217


Equation:
  mean[y]=3.00009091 + 0.50009091 . x

  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    9 3.3166   10       8     13      11   14 11    0
2        y 8.77 2.0317 9.14    8.14   8.74    9.26  8.1 11    0
lambda =  1.79950642693625 
Assuming:
    - Explanatory variable: x;
    - Dependent variable: y.

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = y ~ x)"

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9009 -0.7609  0.1291  0.9491  1.2691 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)    3.001      1.125   2.667  0.02576 * 
x              0.500      0.118   4.239  0.00218 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared:  0.6662,    Adjusted R-squared:  0.6292 
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179


Equation:
  mean[y]=3.00090909 + 0.5 . x

  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    9 3.3166   10       8     13      11   14 11    0
2        y 7.11 2.0304 7.46    6.77  12.74    7.81 8.84 11    0
lambda =  1.80347027989639 
Assuming:
    - Explanatory variable: x;
    - Dependent variable: y.

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = y ~ x)"

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1586 -0.6146 -0.2303  0.1540  3.2411 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0025     1.1245   2.670  0.02562 * 
x             0.4997     0.1179   4.239  0.00218 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared:  0.6663,    Adjusted R-squared:  0.6292 
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176


Equation:
  mean[y]=3.00245455 + 0.49972727 . x

  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    8 3.3166    8       8      8       8    8 11    0
2        y 8.84 2.0306 6.58    5.76   7.71    8.47 7.04 11    0
lambda =  1.80441610212107 
Assuming:
    - Explanatory variable: x;
    - Dependent variable: y.

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = y ~ x)"

Residuals:
   Min     1Q Median     3Q    Max 
-1.751 -0.831  0.000  0.809  1.839 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   3.0017     1.1239   2.671  0.02559 * 
x             0.4999     0.1178   4.243  0.00216 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared:  0.6667,    Adjusted R-squared:  0.6297 
F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165


Equation:
  mean[y]=3.00172727 + 0.49990909 . x

Em todos os gráficos, computamos as RLS com as respectivas bandas de confiança de 95%.

RLS é adequada apenas para o primeiro conjunto de dados. Porém, é somente a saída gráfica que revela este fato. Média, mediana, quartis, os valores de r e \(R^2\), os valores de intercepto e inclinação e os valores p associados são todos similares (demo_Anscombe2.R):

####--- média, d.p. e quartis

Conjunto  1 :
  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    9 3.3166   10       8     13      11   14 11    0
2        y 8.81 2.0316 8.04    6.95   7.58    8.33 9.96 11    0


Conjunto  2 :
  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    9 3.3166   10       8     13      11   14 11    0
2        y 8.77 2.0317 9.14    8.14   8.74    9.26  8.1 11    0


Conjunto  3 :
  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    9 3.3166   10       8     13      11   14 11    0
2        y 7.11 2.0304 7.46    6.77  12.74    7.81 8.84 11    0


Conjunto  4 :
  Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1        x    8 3.3166    8       8      8       8    8 11    0
2        y 8.84 2.0306 6.58    5.76   7.71    8.47 7.04 11    0


####--- reta de regressão (intercepto, inclinação e estatísticas)

Conjunto  1 :
            Estimate Std. Error t value Pr(>|t|)
(Intercept)      3.0      1.125    2.67  0.02573
x                0.5      0.118    4.24  0.00217


Conjunto  2 :
            Estimate Std. Error t value Pr(>|t|)
(Intercept)      3.0      1.125    2.67  0.02576
x                0.5      0.118    4.24  0.00218


Conjunto  3 :
            Estimate Std. Error t value Pr(>|t|)
(Intercept)      3.0      1.124    2.67  0.02562
x                0.5      0.118    4.24  0.00218


Conjunto  4 :
            Estimate Std. Error t value Pr(>|t|)
(Intercept)      3.0      1.124    2.67  0.02559
x                0.5      0.118    4.24  0.00216

RLS é válida apenas no domínio observado da VE?

RLS é válida apenas no domínio observado da variável explicativa.
Fora desse intervalo, o modelo realiza extrapolação, cuja validade depende de pressupostos não verificáveis nos dados.
Portanto, interpretações e previsões devem restringir-se ao domínio amostral da VE.

Repare, em todos os gráficos anteriores, que as retas de regressão e respectivos intervalos de confiança não vão além dos valores de \(x_i\) observados.

O motivo é que não sabemos se o fenômeno manterá o comportamento linear além deste ponto.

E.g., vamos utilizar os dados MASS::Boston, contendo dados de 506 residências nos subúrbios na década de 70 de Boston, MA, EUA. Duas variáveis nos interessam aqui:

• medv … o valor mediano das casas (em milhares de dólares americanos)
• lstat … a porcentagem da população que tem baixo status social (BSS)

Suponha que queiramos usar BSS como preditor do valor das residências com valor igual ou menor que 20 mil dólares (demo_Boston_ate20000.R):

source("eiras.correg.R")
par(mfrow=c(1,1))
Boston <- MASS::Boston
x <- Boston$lstat[Boston$medv<=20]
y <- Boston$medv[Boston$medv<=20]
lst <- correg(x, y,
              B=0, 
              method="lm_robust",
              xlab="BSS (%)", 
              ylab="MedianaCasa (x1000 USD)", 
              bg="transparent", col="black", pch=1)

                 Variable  Mean     Sd  Min. 1st Qu. Median 3rd Qu. Max.   n
1                 BSS (%) 13.27 6.2885 29.93    17.1  20.45   10.26 8.47 215
2 MedianaCasa (x1000 USD)  18.9 3.7873  16.5    18.9     15    18.2 19.9 215
  NA's
1    0
2    0
lambda =  1.87741178290278 
Assuming:
    - Explanatory variable: BSS (%);
    - Dependent variable: MedianaCasa (x1000 USD).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = MedianaCasa (x1000 USD) ~ BSS (%))"

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper  DF
(Intercept)  22.2352    0.66252  33.561 5.395e-87  20.9292  23.5411 213
BSS (%)      -0.3811    0.03919  -9.725 9.903e-19  -0.4583  -0.3038 213

Multiple R-squared:  0.4004 ,   Adjusted R-squared:  0.3976 
F-statistic: 94.57 on 1 and 213 DF,  p-value: < 2.2e-16

Equation:
  mean[MedianaCasa (x1000 USD)]=22.23517138 - 0.38108589 . BSS (%)

Encontramos a reta de regressão:

\[ \widehat{\text{medv}} = 22.24 - 0.38 \; \text{lstat} \]

Como era de se esperar, a relação é negativa: quando maior a porcentagem de BSS, menor o valor da casa do bairro. Podemos, com esta equação, predizer, por exemplo, que em um bairro com 14% de habitantes com baixo status social, o valor mediano da casa é

\[ \widehat{\text{medv}} = 22.24 - 0.38 \times 14 \approx 16.9 \]

I.e., \(\widehat{\text{medv}}(14)\) é quase 17 mil dólares.

Qual o valor mediano de uma casa em um baixo mais rico, como por exemplo, onde 1% dos habitantes são BSS?

A equação prediz:

\[ \widehat{\text{medv}} = 22.24 - 0.38 \times 1 \approx 21.85 \]

I.e., \(\widehat{\text{medv}}(1)\) é quase 22 mil dólares.

No entanto, a consequência de termos estudado as casas com valor de 20 mil dólares ou menos é que só encontramos uma regra válida para esta faixa de valores. Ao estudarmos as casas de menor valor, sem perceber separamos os bairros mais pobres com 7.79% ou mais de pessoas BSS:

cat(min(Boston$lstat[Boston$medv<=20]),"\n")

7.79 

Portanto, os bairros com menos do que 7.79% da população com BSS estão fora da faixa que avaliamos e, portanto, a equação da regressão pode não ser boa preditora para um bairro com 1% de BSS. Caso consideremos todos os dados disponíveis (demo_Boston_todas.R), temos:

source("eiras.correg.R")
Boston <- MASS::Boston
lst <- correg(Boston$lstat, Boston$medv,
              B=0, 
              method="lm_robust",
              xlab="BSS (%)", 
              ylab="MedianaCasa (x1000 USD)", 
              bg="transparent", col="black", pch=1)

                 Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.   n NA's
1                 BSS (%) 2.94 7.1411 4.98    9.14   4.03    5.33 5.21 506    0
2 MedianaCasa (x1000 USD) 33.4 9.1971   24    21.6   34.7    36.2 28.7 506    0
lambda =  0.748982103853726 
Assuming:
    - Explanatory variable: BSS (%);
    - Dependent variable: MedianaCasa (x1000 USD).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = MedianaCasa (x1000 USD) ~ BSS (%))"

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value   Pr(>|t|) CI Lower CI Upper  DF
(Intercept)    34.55    0.75562   45.73 1.742e-181   33.069  36.0384 504
BSS (%)        -0.95    0.04979  -19.08  1.719e-61   -1.048  -0.8522 504

Multiple R-squared:  0.5441 ,   Adjusted R-squared:  0.5432 
F-statistic: 364.1 on 1 and 504 DF,  p-value: < 2.2e-16

Equation:
  mean[MedianaCasa (x1000 USD)]=34.55384088 - 0.95004935 . BSS (%)

A primeira equação não pode ser usada porque subestima o valor da casa ao utilizar lstat fora da faixa de valores que foi utilizada para estimar a reta de regressão. No entanto, esta segunda reta de regressão também não pode ser usada porque, claramente, a relação entre as duas variáveis não é linear.

Há pontos de alavancagem?

O terceiro conjunto de dados do Quarteto de Anscombe mostra o efeito de ponto de alavancagem, com um único ponto visivelmente deslocando a inclinação da reta de regressão. Veja o efeito, em dois exemplos (modificações de demo_regressao.R) nas quais apenas o ponto \((174, 68)\) foi deslocado:

• em demo_regressao_3.R para \((172, 78)\):

# demo_regressao_3.R
# alavancagem produzida por mudança de coordenadas
# do ponto (x3,y3)

source("eiras.friendlycolor.R")

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

# para nao ter um ponto na coordenada do centroide
# (melhora o exemplo)
centro_x <- mean(x)
x[x==centro_x] <- x[x==centro_x]-0.5
# cat("x:",x,"\n")
# cat("y:",y,"\n")

# centroide (media de x, media de y)
centro_x <- mean(x)
centro_y <- mean(y)
# para nao ter um ponto na coordenada do centroide
x[x==centro_x] <- x[x==centro_x]-0.5
centro_x <- mean(x)

# scatterplot
plot(NA, 
     xlab="Estatura (cm)", ylab="Massa (kg)",
     xlim = c(min(x),max(x)), ylim = c(min(y),max(y)),
     pch=21, col="black", bg=friendlycolor(24)
     )

# media de y, traca linha horizontal
y.H0 <- rep(centro_y,length(x))
lines(x, y.H0, col=friendlycolor(30), 
      lwd=3, lty=2)

# modelo linear, traca a reta
modelo <- lm(y ~ x)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
y.medio <- intercepto + inclinacao*x
lines(x, y.medio, 
      col=paste(friendlycolor(8),"44",sep=""), 
      lwd=3, lty=1
)
points(x, y,
       pch=21, col="darkgray", 
       bg=paste(friendlycolor(24),"44",sep="")
)

# centroide
points(centro_x, centro_y, pch=21, col="black", bg="black")

# modelo linear, traca a reta
original.x <- x[3]
original.y <- y[3]
x[3] <- 167
y[3] <- 78
modelo <- lm(y ~ x)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
y.medio <- intercepto + inclinacao*x
lines(x, y.medio, 
      col=friendlycolor(8), 
      lwd=3, lty=1
)
points(x, y,
       pch=21, col="black", 
       bg=friendlycolor(24)
)
arrows(original.x, original.y, x[3], y[3], length = 0.2, lty=2)

# legenda
legend("topleft",
       c("y observado","média de y",
         paste("reta: y_médio = ",round(intercepto,3)," + ",
               round(inclinacao,3),"x",sep="")
       ), 
       lwd=c(NA, 3, 3), 
       pch=c(21, NA, NA),
       col=c("black",friendlycolor(30), friendlycolor(8)), 
       pt.bg=c(friendlycolor(24),friendlycolor(30), friendlycolor(8)), 
       cex=0.9, bty="n")

• em demo_regressao_3b.R para \((180, 63)\):

# demo_regressao_3.R
# alavancagem produzida por mudança de coordenadas
# do ponto (x3,y3)

source("eiras.friendlycolor.R")

# valores
x <- estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
y <- massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
# para nao ter um ponto na coordenada do centroide
# (melhora o exemplo)
centro_x <- mean(x)
x[x==centro_x] <- x[x==centro_x]-0.5
# cat("x:",x,"\n")
# cat("y:",y,"\n")

# centroide (media de x, media de y)
centro_x <- mean(x)
centro_y <- mean(y)
# para nao ter um ponto na coordenada do centroide
x[x==centro_x] <- x[x==centro_x]-0.5
centro_x <- mean(x)

# scatterplot
plot(NA, 
     xlab="Estatura (cm)", ylab="Massa (kg)",
     xlim = c(min(x),max(x)), ylim = c(min(y),max(y)),
     pch=21, col="black", bg=friendlycolor(24)
     )

# media de y, traca linha horizontal
y.H0 <- rep(centro_y,length(x))
lines (x, y.H0, col=friendlycolor(30), 
       lwd=3, lty=2)

# modelo linear, traca a reta
modelo <- lm(y ~ x)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
y.medio <- intercepto + inclinacao*x
lines (x, y.medio, 
       col=paste(friendlycolor(8),"44",sep=""), 
       lwd=3, lty=1
       )
points(x, y,
     pch=21, col="darkgray", 
     bg=paste(friendlycolor(24),"44",sep="")
     )

# centroide
points(centro_x, centro_y, pch=21, col="black", bg="black")

# modelo linear, traca a reta
original.x <- x[3]
original.y <- y[3]
x[3] <- 180
y[3] <- 63
modelo <- lm(y ~ x)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
y.medio <- intercepto + inclinacao*x
lines (x, y.medio, 
       col=friendlycolor(8), 
       lwd=3, lty=1
)
points(x, y,
       pch=21, col="black", 
       bg=friendlycolor(24)
)
arrows(original.x, original.y, x[3], y[3], length = 0.2, lty=2)

# legenda
legend ("topleft",
        c("y observado","média de y",
          paste("reta: y_médio = ",round(intercepto,3)," + ",
                round(inclinacao,3),"x",sep="")
        ), 
        lwd=c(NA, 3, 3), 
        pch=c(21, NA, NA),
        col=c("black",friendlycolor(30), friendlycolor(8)), 
        pt.bg=c(friendlycolor(24),friendlycolor(30), friendlycolor(8)), 
        cex=0.9, bty="n")

• RLS_alavancagem.R

Call:
lm(formula = massa ~ estatura)

Residuals:
   Min     1Q Median     3Q    Max 
 -7.75  -2.24   1.46   2.51   8.00 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -144.53      76.13   -1.90    0.099 .
estatura        1.26       0.44    2.86    0.024 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.71 on 7 degrees of freedom
Multiple R-squared:  0.538, Adjusted R-squared:  0.472 
F-statistic: 8.16 on 1 and 7 DF,  p-value: 0.0245

Call:
lm(formula = massa ~ estatura)

Residuals:
   Min     1Q Median     3Q    Max 
-11.13  -6.72  -2.26   6.01  12.78 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   99.800    118.082    0.85     0.43
estatura      -0.137      0.682   -0.20     0.85

Residual standard error: 8.86 on 7 degrees of freedom
Multiple R-squared:  0.00572,   Adjusted R-squared:  -0.136 
F-statistic: 0.0403 on 1 and 7 DF,  p-value: 0.847

RLS robusta

RLS robusta estimatr::lm_robust não precisa da suposição de homocedasticidade.

RLS robusta por bootstrapping necessita apenas da suposição de independência entre os pares de observações.

Sempre que possível, utilize RLS robusta.

Primeiro estudo: hematócrito

• Arquivo de dados: Gestante.rds

Já observamos nos gráficos de dispersão do hematócrito contra hemoglobina, hemácias e leucócitos porque os pares de medida (das gestantes nas quais as quatro variáveis foram medidas) podem ser considerados independentes e seus gráficos de dispersão pareciam acomodar bem o ajuste de uma reta.

No entanto, não verificamos as outras duas suposições necessárias para os testes estatísticos:

• uninormalidade da VD (regressão) e binormalidade (correlação)
• homocedasticidade da VD nos níveis da VE (regressão)

A função utilizada até aqui foi lm do pacote stats. Existe uma alternativa, mais robusta, estimatr::lm_robust que dispensa a suposição de homocedasticidade. Recomendamos adotá-la sempre que possível no lugar de lm.

Judkins, 2015, https://doi.org/10.1002/sim.6839

Para situações ainda mais difíceis, usaremos o método bootstrapping, que tolera todas as violações, com exceção de dependência entre os pares de observações.

O teste de normalidade na RLS é univariado, aplicado aos resíduos ou diretamente aos valores da VD.

Quanto à normalidade da VD, verificaremos a aparência de gráficos de densidade (linhas sólidas) sobre os quais apresentaremos a distribuição normal que tem a mesma média e desvio-padrão dos dados (linhas pontilhadas).

A não rejeição de hipótese nula não prova que a distruição da variável é normal.

• demo_LEUCnormal.R

Gestantes <- readRDS("Gestante.rds")

plot_normcmp <- function(v, xlab) {
        v <- na.omit(v)
        d  <- density(v)
        m  <- mean(v); s <- sd(v)
        x  <- seq(m-3*s, m+3*s, length.out=200)
        y  <- dnorm(x, mean=m, sd=s)
        plot(d,
             xlim=c(min(x), max(x)),
             ylim=c(0, max(y, d$y, na.rm=TRUE)),
             main="",
             xlab=xlab, ylab="densidade", type="l")
        lines(x, y, lty=2)
        rug(v)
        legend("topright", lty=c(1,2), bty="n",
               legend=c("estimada", "normal"))
}

op <- par(no.readonly=TRUE)
par(mfrow=c(2,2), mar=c(4,4,2,1))

plot_normcmp(Gestantes$HT,   "Hematócrito (%)")
plot_normcmp(Gestantes$HB,   "Hemoglobina (mg/dL)")
plot_normcmp(Gestantes$HEM,  "Hemácias (milhões/mm³)")
plot_normcmp(Gestantes$LEUC, "Leucócitos (milhares/mm³)")

par(op)  # volta ao 1x1 padrão

O teste de uninormalidade aplicado é o teste de Shapiro-Wilk que está implementado em MVN::mvn. demo_normalidade.R A saída é: Gestantes <- readRDS("Gestante.rds") cat("\n\nTotal de Gestantes:",nrow(Gestantes)) cat("\n\nHT x HB:") result <- MVN::mvn(data=Gestantes[,c("HT","HB","HEM","LEUC")], univariate_test="SW") cat("\nTeste de uninormalidade:\n") print(result$univariate_normality, digits=4) Total de Gestantes: 40 HT x HB: Teste de uninormalidade: Test Variable Statistic p.value Normality 1 Shapiro-Wilk HT 0.886 <0.001 ✗ Not normal 2 Shapiro-Wilk HB 0.928 0.014 ✗ Not normal 3 Shapiro-Wilk HEM 0.951 0.08 ✓ Normal 4 Shapiro-Wilk LEUC 0.948 0.067 ✓ Normal Observe os valores p. A hipótese nula de uninormalidade são rejeitadas em alguns casos: HT e HB. A não rejeição da uninormalidade da VD não prova que a distribuição da variável é normal. No entanto, como $n=$40, RLS é válida com base no teorema central do limite, pois o número de pares independentes, 40, é maior que 30. Conforme Wooldridge (2016, p. 155), “Sabemos que a normalidade não desempenha nenhum papel na não-viesamento do OLS, nem afeta a conclusão de que o OLS é o melhor estimador linear não-viesado sob as suposições de Gauss-Markov. Mas a inferência exata baseada nas estatísticas t e F requer normalidade de y. Isso significa que, em nossa análise anterior no Exemplo 4.6, devemos abandonar a estatística t para determinar quais variáveis são estatisticamente significantes? Felizmente, a resposta a esta pergunta é não. Mesmo que y não seja de uma distribuição normal, podemos usar o teorema do limite central do Apêndice C para concluir que os estimadores OLS satisfazem a normalidade assintótica, o que significa que eles são distribuídos aproximadamente normalmente em tamanhos de amostra suficientemente grandes. Isso sugere que a normalidade dos erros é realmente uma condição suficiente, mas não uma condição necessária. Suficiente para quê? Suficiente para garantir a normalidade (aproximada) das distribuições amostrais dos estimadores dos parâmetros do modelo (ou seja, os coeficientes).” Uma VE binária, e.g. Sexo, pode ser usada na análise de RLS porque ela é intervalar. Uma justificativa teórica e exemplo de aplicação com massa corporal total e sexo encontram-se em Howell (2013, Chapter 10: Alternative Correlational Techniques, p. 305-7). Para amostras de tamanho pequeno, RLS pode ser válida, pois a suposição de normalidade é condição suficiente, não é condição necessária para sua validade. Uma alternativa é utilizar o método bootstrapping, que prescinde da suposição de uninormalidade da VD e pode ser usada para qualquer tamanho de amostra.

A função disponível em eiras.correg.R pode ser chamada para correlação e regressão, usando lm (método clássico), estimatr::lm_robust (método robusto à heterocedasticidade) ou bootstrapping (robusto à heterocedasticidade, não normalidade e amostra pequena), bastando alterar seus parâmetros.

Vejamos a diferença de resultados nas três situações. Gestantes_RLS.R implementa a chamada com lm (consulte o código para verificar como isto foi feito):

# Gestantes_RLS.R
# (versao nao padronizada, para regressão)

alfa <- 0.05

source("eiras.jitter.R")
source("eiras.friendlycolor.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")
B <- 0
# HT x HB
reslm_hb <- correg(Gestantes$HT, Gestantes$HB,
                   alpha=alfa, B=B, 
                   method="lm",
                   xlab="Hematócrito (%)", ylab="Hemoglobina (mg/dl)", 
                   bg="transparent", col=col_HB, pch=pch_HB)
# HT x HEM
reslm_hem <- correg(Gestantes$HT, Gestantes$HEM,
                    alpha=alfa, B=B, 
                    method="lm",
                    xlab="Hematócrito (%)", ylab="Hemácia (milhão/mm³)", 
                    bg="transparent", col=col_HEM, pch=pch_HEM)
# HT x LEUC
reslm_leuc <- correg(Gestantes$HT, Gestantes$LEUC,
                     alpha=alfa, B=B, 
                     method="lm",
                     xlab="Hematócrito (%)", ylab="Leucócito (milhar/mm³)", 
                     bg="transparent", col=col_LEUC, 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).

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = Hemoglobina (mg/dl) ~ Hematócrito (%))"

Residuals:
     Min       1Q   Median       3Q      Max 
-1.40302 -0.16455  0.04438  0.31045  0.89698 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.001109   1.232719   0.001    0.999    
Hematócrito (%) 0.348700   0.032884  10.604 6.52e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5319 on 38 degrees of freedom
Multiple R-squared:  0.7474,    Adjusted R-squared:  0.7408 
F-statistic: 112.4 on 1 and 38 DF,  p-value: 6.523e-13


Equation:
  mean[Hemoglobina (mg/dl)]=0.00110856 + 0.34870031 . Hematócrito (%)

              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ácia (milhão/mm³) 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ácia (milhão/mm³).

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = Hemácia (milhão/mm³) ~ Hematócrito (%))"

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6459 -0.1843 -0.0059  0.1362  0.7141 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      0.52196    0.65429   0.798     0.43    
Hematócrito (%)  0.10150    0.01745   5.815 1.02e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2823 on 38 degrees of freedom
Multiple R-squared:  0.4709,    Adjusted R-squared:  0.457 
F-statistic: 33.82 on 1 and 38 DF,  p-value: 1.02e-06


Equation:
  mean[Hemácia (milhão/mm³)]=0.52195719 + 0.10149847 . Hematócrito (%)

                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ócito (milhar/mm³)  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ócito (milhar/mm³).

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = Leucócito (milhar/mm³) ~ Hematócrito (%))"

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3459 -1.0601 -0.3189  0.8865  6.7649 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)      5.42317    4.72211   1.148    0.258
Hematócrito (%)  0.08922    0.12597   0.708    0.483

Residual standard error: 2.037 on 38 degrees of freedom
Multiple R-squared:  0.01303,   Adjusted R-squared:  -0.01294 
F-statistic: 0.5017 on 1 and 38 DF,  p-value: 0.4831


Equation:
  mean[Leucócito (milhar/mm³)]=5.42316514 + 0.08922018 . Hematócrito (%)

Observamos nas saídas as retas de regressão, a equação da reta, o coeficiente de determinação (\(R^2\)), e a banda de confiança de 95% (linhas tracejadas), onde esperamos que a reta de regressão populacional esteja contida com alta probabilidade.

O valor p é referente ao teste da hipótese nula:

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

Portanto, RLS é estatisticamente plausível para estimar a concentração de hemoglobina e a contagem de hemácias a partir do hematócrito:

• para hemoglobina: \(\widehat{\text{HB}} = 0.00 + 0.35 \times \text{HT}, \; R^2=75\%\)
• para hemácia: \(\widehat{\text{HEM}} = 0.52 + 0.10 \times \text{HT}, \; R^2=47\%\)

A impressão visual de que a predição para hemoglobina devia ser melhor que a para hemácias tem apoio na análise de RLS.

Compare os resultados com o uso de estimatr::lm_robust implementada em Gestantes_RLS_robust.R:

# Gestantes_RLS.R
# (versao nao padronizada, para regressão)

alfa <- 0.05

source("eiras.jitter.R")
source("eiras.friendlycolor.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
B <- 0
reslmr_hb <- correg(Gestantes$HT, Gestantes$HB,
                    alpha=alfa, B=B, 
                    method="lm_robust",
                    xlab="Hematócrito (%)", ylab="Hemoglobina (mg/dl)", 
                    bg="transparent", col=col_HB, pch=pch_HB)
# HT x HEM
reslmr_hem <- correg(Gestantes$HT, Gestantes$HEM,
                     alpha=alfa, B=B, 
                     method="lm_robust",
                     xlab="Hematócrito (%)", ylab="Hemácia (milhão/mm³)", 
                     bg="transparent", col=col_HEM, pch=pch_HEM)
# HT x LEUC
reslmr_leuc <- correg(Gestantes$HT, Gestantes$LEUC,
                      alpha=alfa, B=B, 
                      method="lm_robust",
                      xlab="Hematócrito (%)", ylab="Leucócitos (milhar/mm³)", 
                      bg="transparent", col=col_LEUC, 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).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Hemoglobina (mg/dl) ~ Hematócrito (%))"

Standard error type:  HC2 

Coefficients:
                Estimate Std. Error   t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)     0.001109    1.10896 9.996e-04 9.992e-01  -2.2439   2.2461 38
Hematócrito (%) 0.348700    0.02956 1.180e+01 2.847e-14   0.2889   0.4085 38

Multiple R-squared:  0.7474 ,   Adjusted R-squared:  0.7408 
F-statistic: 139.2 on 1 and 38 DF,  p-value: 2.847e-14

Equation:
  mean[Hemoglobina (mg/dl)]=0.00110856 + 0.34870031 . Hematócrito (%)

              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ácia (milhão/mm³) 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ácia (milhão/mm³).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Hemácia (milhão/mm³) ~ Hematócrito (%))"

Standard error type:  HC2 

Coefficients:
                Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)       0.5220    0.72153  0.7234 4.739e-01 -0.93870   1.9826 38
Hematócrito (%)   0.1015    0.01933  5.2512 6.042e-06  0.06237   0.1406 38

Multiple R-squared:  0.4709 ,   Adjusted R-squared:  0.457 
F-statistic: 27.58 on 1 and 38 DF,  p-value: 6.042e-06

Equation:
  mean[Hemácia (milhão/mm³)]=0.52195719 + 0.10149847 . Hematócrito (%)

                 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/mm³)  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/mm³).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Leucócitos (milhar/mm³) ~ Hematócrito (%))"

Standard error type:  HC2 

Coefficients:
                Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)      5.42317     4.0947  1.3244   0.1933  -2.8662  13.7125 38
Hematócrito (%)  0.08922     0.1066  0.8368   0.4080  -0.1266   0.3051 38

Multiple R-squared:  0.01303 ,  Adjusted R-squared:  -0.01294 
F-statistic: 0.7002 on 1 and 38 DF,  p-value: 0.408

Equation:
  mean[Leucócitos (milhar/mm³)]=5.42316514 + 0.08922018 . Hematócrito (%)

Note que, neste caso, há pequenas diferenças na estimativa dos valores p.

Em Gestantes_RLS_bootstrapping.R chamamos as mesmas funcões, mas agora usamos bootstrapping (veja o código para ver como as chamadas foram modificadas):

# Gestantes_RLS.R
# (versao nao padronizada, para regressão)

alfa <- 0.05

source("eiras.jitter.R")
source("eiras.friendlycolor.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")
B <- 1e4
# HT x HB
reslmb_hb <- correg(Gestantes$HT, Gestantes$HB,
                   alpha=alfa, B=B, 
                   method="lm",
                   xlab="Hematócrito (%)", ylab="Hemoglobina (mg/dl)", 
                   bg="transparent", col=col_HB, pch=pch_HB)
# HT x HEM
reslmb_hem <- correg(Gestantes$HT, Gestantes$HEM,
                    alpha=alfa, B=B, 
                    method="lm",
                    xlab="Hematócrito (%)", ylab="Hemácia (milhão/mm³)", 
                    bg="transparent", col=col_HEM, pch=pch_HEM)
# HT x LEUC
reslmb_leuc <- correg(Gestantes$HT, Gestantes$LEUC,
                     alpha=alfa, B=B, 
                     method="lm",
                     xlab="Hematócrito (%)", ylab="Leucócito (milhar/mm³)", 
                     bg="transparent", col=col_LEUC, 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).

------------------------
Simple linear regression
------------------------

    SLR: 10000 resamples

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

alternative hypothesis: true slope is not equal to 0

             median        min      1st.Q     3rd.Q       max     IC95.1
intercept 0.1011784 -4.3394297 -0.5867508 0.8382916 6.5450053 -1.7869133
slope     0.3461634  0.1683491  0.3263087 0.3645140 0.4534232  0.2780912
             IC95.2
intercept 2.5955660
slope     0.3947574

Equation:
  mean[Hemoglobina (mg/dl)] = 0.10117845 [-1.78691327, 2.59556597] + 
        0.34616339 [0.27809123, 0.39475743] . Hematócrito (%)

              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ácia (milhão/mm³) 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ácia (milhão/mm³).

------------------------
Simple linear regression
------------------------

    SLR: 10000 resamples

data:  Hemácia (milhão/mm³) ~ Hematócrito (%)

alternative hypothesis: true slope is not equal to 0

             median         min      1st.Q     3rd.Q       max      IC95.1
intercept 0.5241211 -1.30399051 0.10481483 1.0294513 3.7109177 -0.57802789
slope     0.1013580  0.01573971 0.08793501 0.1126387 0.1492341  0.05780478
             IC95.2
intercept 2.1475746
slope     0.1306509

Equation:
  mean[Hemácia (milhão/mm³)] = 0.52412108 [-0.57802789, 2.14757463] + 
        0.10135803 [0.05780478, 0.13065088] . Hematócrito (%)

                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ócito (milhar/mm³)  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ócito (milhar/mm³).

------------------------
Simple linear regression
------------------------

    SLR: 10000 resamples

data:  Leucócito (milhar/mm³) ~ Hematócrito (%)

alternative hypothesis: true slope is not equal to 0

              median         min      1st.Q     3rd.Q       max     IC95.1
intercept 5.27812964 -12.2029412 2.52334868 7.9545219 24.830250 -3.7385011
slope     0.09191454  -0.4079691 0.02326986 0.1646922  0.575463 -0.1134645
              IC95.2
intercept 13.2846824
slope      0.3313723

Equation:
  mean[Leucócito (milhar/mm³)] = 5.27812964 [-3.73850106, 13.28468236] + 
        0.09191454 [-0.11346449, 0.33137228] . Hematócrito (%)

Esta modalidade robusta não reporta valores p. A decisão é tomada pelo intervalo de confiança de 95%: veja que o intervalo de \(R^2\) não inclui o zero para hemoglobina ou contagem de hemácias em função do hematócrito. Portanto, podemos assumir que suas inclinações são diferentes de zero. No caso de hematócrito vs. contagem de leucócitos, o valor 0 pertence ao intervalo de confiança de 95%.

Esta versão de nossa função com bootstrapping cria uma sombra feita por várias retas e a banda de confiança 95% baseada nestas retas de regressão. A banda é a região dentro da qual espera, com confiança (probablidade) de 95%, encontrar o segmento de reta populacional.

Portanto, a outra forma de se ver a significância destes três casos é verificar se é possível traçar uma reta horizontal (\(H_0\))que fique inteiramente contida na banda de confiança, o que apenas acontece quando não rejeitamos \(H_0\). Verifique que não é possível acomodar uma reta horizontal no caso da hemoblogina ou da contagem de hemácias, mas podemos fazê-lo no caso da contagem de leucócitos.

O preço a pagar em bootstrapping é tempo de processamento. Em geral, precisamos de 10.000 (dez mil, 1e4) a 1.000.000 (um milhão, 1e6) reamostragens, o que demanda muito mais tempo de processamento do que as versões analíticas de lm e lm_robust.

Neste exemplo, chegamos às mesmas conclusões com os três métodos:

VE × VD	Método	Inclinação	Reta de regressão
HT × HB	Convencional (lm)	0.349 [0.282,0.415] p=6.52e-13	mean[Hemoglobina (mg/dl)]=0.00110856 + 0.34870031 . Hematócrito (%)
	Robusto (lm_robust)	0.349 [0.289,0.409] p=2.85e-14	mean[Hemoglobina (mg/dl)]=0.00110856 + 0.34870031 . Hematócrito (%)
	Robusto (Bootstrapping)	0.346 [0.278,0.395]	mean[Hemoglobina (mg/dl)] = 0.10117845 [-1.78691327, 2.59556597] + 0.34616339 [0.27809123, 0.39475743] . Hematócrito (%)
HT × HEM	Convencional (lm)	0.101 [0.066,0.137] p=1.02e-06	mean[Hemácia (milhão/mm³)]=0.52195719 + 0.10149847 . Hematócrito (%)
	Robusto (lm_robust)	0.101 [0.062,0.141] p=6.04e-06)	mean[Hemácia (milhão/mm³)]=0.52195719 + 0.10149847 . Hematócrito (%)
	Robusto (Bootstrapping)	0.101 [0.058,0.131]	mean[Hemácia (milhão/mm³)] = 0.52412108 [-0.57802789, 2.14757463] + 0.10135803 [0.05780478, 0.13065088] . Hematócrito (%)
HT × LEUC	Convencional (lm)	0.089 [-0.166,0.344] p=4.83e-01)	mean[Leucócito (milhar/mm³)]=5.42316514 + 0.08922018 . Hematócrito (%)
	Robusto (lm_robust)	0.089 [-0.127,0.305] p=4.08e-01	mean[Leucócitos (milhar/mm³)]=5.42316514 + 0.08922018 . Hematócrito (%)
	Robusto (Bootstrapping)	0.092 [-0.113,0.331]	mean[Leucócito (milhar/mm³)] = 5.27812964 [-3.73850106, 13.28468236] + 0.09191454 [-0.11346449, 0.33137228] . Hematócrito (%)

• A função lm não fornece o intervalo de confiança, mas pudemos obtê-los com confint(lm(y~x, …)).

Segundo estudo: relação entre massas corporal e cerebral em mamífero

A motivação central do artigo de Gerhardt von Bonin (1937) é esclarecer a relação quantitativa entre massa do cérebro (\(E\)) e massa corporal (\(S\)) em mamíferos, buscando uma lei empírica que descreva essa relação de forma estatisticamente válida, sem recorrer a hipóteses sobre “inteligência” ou funções “psíquicas”. Ele critica os autores Snell (1891) e Dubois (1898, 1914) por assumirem, sem base empírica suficiente, que a inteligência correlaciona-se diretamente com o tamanho do cérebro.

A conjectura do expoente da função potência ser aproximadamente \(2/3=6/9\) vem de Snell, que argumentou teoricamente que, se o corpo cresce em três dimensões, mas a massa encefálica cresce com limitações geométricas e funcionais (como de superfície cortical e conexões neurais), então o cérebro aumentaria com uma taxa sublinear em relação ao corpo, próxima a

\[ E = k S^{2/3} \] Dubois depois estimou empiricamente o expoente como \(5/9 \approx 0.56\), enquanto von Bonin, usando dados mais amplos e análise de regressão log-log, encontrou \(E = 0.18\, S^{0.655}\), portanto muito próximo de \(2/3\). A função potência estimada é inelástica com elasticidade constante de \(2/3\). A deficiência do análise estatística de von Bonin é que ela é apenas descritiva, pois faltou determinar o intervalo de confiança de 95% do expoente.

Em resumo:

• Motivação: descrever empiricamente a relação entre massa corporal e cerebral sem pressupostos psicológicos.
• Conjectura \(2/3\): derivada de argumentos geométricos e fisiológicos sobre escalonamento alométrico; o cérebro cresce menos que proporcionalmente ao corpo, seguindo uma lei de potência com expoente aproximadamente \(2/3\).

Von Bonin (1937) utilizou espécies adultas de mamíferos com dados confiáveis de massa corporal e massa encefálica, escolhidas para representar ampla diversidade taxonômica.

Ele compilou dados principalmente de Warncke (1908) e Hrdlička (1925), abrangendo cerca de 90 espécies, agrupadas em ordens:

Incluindo humanos e cetáceos, a amostra de mamíferos utilizada ou referida por Gerhardt von Bonin (1937) compreende as seguintes ordens e exemplos representativos:

• Primates: Homo sapiens (excluído da regressão principal, mas citado), Pan troglodytes, Gorilla gorilla, Hylobates spp., Macacus spp., Cebus capucinus, Ateles spp., Lemur spp., Nycticebus tardigradus, Galago galago, entre outros.

• Carnivora: Ursus arctos, Canis lupus, Felis domestica, Panthera leo, Panthera tigris, Vulpes spp., Hyaena spp.

• Rodentia: Mus musculus, Rattus norvegicus, Sciurus spp., Castor canadensis, Cavia porcellus, Oryctolagus cuniculus.

• Chiroptera: Pteropus edulis, Rhinolophus spp., Vespertilio serotinus.

• Insectivora: Erinaceus europaeus, Talpa europaea.

– Ungulata: Equus caballus, Bos taurus, Hippopotamus amphibius, Cervus canadensis, Elephas africanus, Elephas indicus.

• Cetacea (mencionados e excluídos da regressão por razões metodológicas): Phocaena communis, Lagenorrhynchus albirostris, Tursiops tursio, Globicephalus melas, Balaenoptera musculus, Balaenoptera rostrata, Megaptera boops.

• Edentata: Bradypus tridactylus, Myrmecophaga jubata, Dasypus sexcinctus.

• Marsupialia: Didelphis marsupialis, Trichosurus vulpecula, Onychogale frenata.

ou popularmente,

• Primatas: ser humano, chimpanzé, gorila, gibão, macaco-rhesus, macaco-prego, bugio, macaco-aranha, sagui, lêmure, loris, galago.

• Carnívoros: urso-pardo, lobo, raposa, hiena, gato-doméstico, leão, tigre, leopardo.

• Roedores: rato, camundongo, esquilo, castor, cobaia, coelho.

• Morcegos (quirópteros): raposa-voadora, morcego-ferradura, morcego-serotino.

• Insetívoros: ouriço-cacheiro, toupeira-europeia, musaranho-de-Java.

• Ungulados (herbívoros de casco): cavalo, zebra, boi, veado, hipopótamo, elefante-africano, elefante-indiano.

• Cetáceos: boto-comum, golfinho-nariz-de-garrafa, golfinho-branco, baleia-negra, baleia-minke, baleia-azul, jubarte.

• Edentados: preguiça-de-três-dedos, tamanduá-bandeira, tatu-peba.

• Marsupiais: gambá, canguru, cuíca.

Critérios de seleção:

1. Só indivíduos adultos, com massa do cérebro e do corpo medidas.
2. Dados de fontes comparáveis (cérebro fresco, mesma unidade).
3. Exclusão de humano e cetáceos — o primeiro por especialização cognitiva extrema; o segundo por massa corporal inflada por gordura subcutânea.
4. Uso de médias por espécie quando possível, evitando valores isolados.

Essas espécies foram usadas não por interesse zoológico específico, mas para testar empiricamente a lei alométrica geral entre cérebro e corpo em mamíferos, cobrindo desde pequenos roedores até grandes primatas, ungulados e cetáceos.

Incluímos humano e cetáceos na análise estatística e o resultado de van Bonin se manteve.

Podemos usar o tamanho do corpo para explicar e prever o tamanho do cérebro em mamíferos a partir de alguns animais conhecidos?

A resposta pode ser obtida de duas maneiras: teórica e empírica.

A maneira teórica é partir da relação dimensional alométrica em mamífero:

\[ M_E \propto M_S^b \] sendo \(M_E\) a massa cerebral, \(M_S\) a massa corporal total e \(0<b<1\) o parâmetro alométrico. Dessa forma, a relação entre as massas é inelástica constante igual a \(b\). Alguns pesquisadores conjecturam \(b=2/3\).

A relação alométrica produz a função potência:

\[ E=kS^b \] A função potência é não linear, mas linearizável:

\[ \ln(E)=\ln(k)+b\ln(S) \] Dessa forma, a equação de regressão linear simples para estimar o intervalo de confiança de 95% de \(b\) é:

\[ \ln(E_i)=a+b\ln(S_i)+\varepsilon_i\\ i=1,2,\ldots,n \] A segunda maneira é simetrizar as variáveis \(E\) e \(S\) por meio de transformação potência de Tukey. A transformação que maximiza a simetria das duas variáveis é muito provavelmente a logarítmica. Esta transformação torna a relação entre as duas variáveis mais linear, homocedástica e com menos pontos de alavancagem e discreptantes.

Vamos utilizar os dados colecionados em CorpoCerebroVonBonin.xlsx (von Bonin, 1937).

Esta planilha tem a massa (g) de corpos e cérebros de 119 animais.

RLS por lm_robust e lm mostram os seguintes resultados:

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))
CorpoCerebro$Grupo <- factor(CorpoCerebro$Grupo)
CorpoCerebro$Animal <- factor(CorpoCerebro$Animal)
print(summary(CorpoCerebro))

       Grupo                             Animal       Cerebro       
 Primates :41   Balaenoptera rostrata       :  2   Min.   :   0.31  
 Carnivora:18   Alactaga saliens            :  1   1st Qu.:  11.96  
 Rodentia :17   Alouata nigra               :  1   Median :  64.83  
 Ungulata :12   Anthropopithecus troblodytes:  1   Mean   : 368.88  
 Prosimiae: 9   Asinus asinus               :  1   3rd Qu.: 168.25  
 Cetacea  : 8   Ateles ater                 :  1   Max.   :7000.00  
 (Other)  :15   (Other)                     :113                    
     Corpo         
 Min.   :      19  
 1st Qu.:    1028  
 Median :    3536  
 Mean   : 1612502  
 3rd Qu.:   22089  
 Max.   :73997000  
                   

modelo <- estimatr::lm_robust(CorpoCerebro$Cerebro ~ CorpoCerebro$Corpo)
sumario <- summary(modelo)
print(sumario, digits=5)

Call:
estimatr::lm_robust(formula = CorpoCerebro$Cerebro ~ CorpoCerebro$Corpo)

Standard error type:  HC2 

Coefficients:
                     Estimate Std. Error t value   Pr(>|t|)   CI Lower
(Intercept)        2.5101e+02 6.4370e+01  3.8995 0.00016065 1.2354e+02
CorpoCerebro$Corpo 7.3096e-05 2.0302e-05  3.6005 0.00046568 3.2894e-05
                     CI Upper  DF
(Intercept)        3.7848e+02 118
CorpoCerebro$Corpo 1.1330e-04 118

Multiple R-squared:  0.46371 ,  Adjusted R-squared:  0.45916 
F-statistic: 12.964 on 1 and 118 DF,  p-value: 0.00046568

modelm <- lm(CorpoCerebro$Cerebro ~ CorpoCerebro$Corpo)
sumariolm <- summary(modelm)
print(sumariolm, digits=5)

Call:
lm(formula = CorpoCerebro$Cerebro ~ CorpoCerebro$Corpo)

Residuals:
     Min       1Q   Median       3Q      Max 
-2311.25  -239.73  -191.20  -101.52  4607.81 

Coefficients:
                     Estimate Std. Error t value  Pr(>|t|)    
(Intercept)        2.5101e+02 6.9867e+01  3.5927 0.0004785 ***
CorpoCerebro$Corpo 7.3096e-05 7.2366e-06 10.1009 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 754.6 on 118 degrees of freedom
Multiple R-squared:  0.46371,   Adjusted R-squared:  0.45916 
F-statistic: 102.03 on 1 and 118 DF,  p-value: < 2.22e-16

A saída textual mostra, aproximadamente:

• \(R^2 =\) 46.4%
• \(p =\) 4.66e-04

Portanto, para \(\alpha=0.05\), rejeitamos \(H_0\) e RLS é estatisticamente plausível. Podemos, então, utilizar os valores de intercepto e inclinação da reta para prevermos a massa cerebral (mais complicado para medir, especialmente com o animal vivo) a partir da massa corporal (talvez mais fácil, dependendo do animal), ambas dadas em grama, aproximadamente, com a reta de regressão:

\[ \widehat{\text{cérebro}} = 251 + 7.31 \cdot 10^{-5} \times \text{corpo} \]

RLS parece ser adequada. No entanto, não fizemos uma análise descritiva e gráfica, a qual é sempre recomendada. Observe, verificando com a RLS com lm_robust e as funções gráficas desenvolvidas (CorpoCerebro_RLS.R):

# CorpoCerebro_RLS.R
# (versao nao padronizada, para regressão)

alfa <- 0.05
B <- 0
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.correg.R")

col_Corpo <- friendlycolor(14) # verde
pch_Corpo <- 22
col_Cerebro <- friendlycolor(39) # cinza
pch_Cerebro <- 23
col_CxC <- friendlycolor(39) # cinza
pch_CxC <- 21

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))

correg(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
       main="estimatr::lm_robust",
       alpha=alfa, B=B, 
       method="lm_robust",
       xlab="Corpo (g)", ylab="Cérebro (g)", 
       bg="transparent", col="black", pch=pch_CxC)

correg(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
       main="lm",
       alpha=alfa, B=B, 
       method="lm",
       xlab="Corpo (g)", ylab="Cérebro (g)", 
       bg="transparent", col="black", pch=pch_CxC)

fit <- lm(data=CorpoCerebro,
          Cerebro~Corpo)
shapiro.test(resid(fit))
lmtest::harvtest(fit)
lmtest::bptest(fit)
plot(fit, ask=FALSE, which = 5)
DescTools::PlotBag(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
                   xlab="Corpo (g)", ylab="Cérebro (g)", cex=0.6)
print(REAT::curvefit(CorpoCerebro$Corpo, CorpoCerebro$Cerebro, print.results=FALSE)$models_comp, digits=3)

source("CorpoCerebro_RLS.R")

     Variable Mean           Sd  Min. 1st Qu. Median 3rd Qu. Max.   n NA's
1   Corpo (g) 9490 9558998.2198 63400   22045  90525    5628 5454 120    0
2 Cérebro (g)  130    1026.0897   400     345    425    97.3  103 120    0
lambda =  14228.7891413801 
Assuming:
    - Explanatory variable: Corpo (g);
    - Dependent variable: Cérebro (g).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Cérebro (g) ~ Corpo (g))"

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value  Pr(>|t|)  CI Lower  CI Upper  DF
(Intercept) 2.51e+02  6.437e+01   3.899 0.0001607 1.235e+02 3.785e+02 118
Corpo (g)   7.31e-05  2.030e-05   3.601 0.0004657 3.289e-05 1.133e-04 118

Multiple R-squared:  0.4637 ,   Adjusted R-squared:  0.4592 
F-statistic: 12.96 on 1 and 118 DF,  p-value: 0.0004657

Equation:
  mean[Cérebro (g)]=251.01075929 + 7.31e-05 . Corpo (g)

     Variable Mean           Sd  Min. 1st Qu. Median 3rd Qu. Max.   n NA's
1   Corpo (g) 9490 9558998.2198 63400   22045  90525    5628 5454 120    0
2 Cérebro (g)  130    1026.0897   400     345    425    97.3  103 120    0
lambda =  14228.7891413801 
Assuming:
    - Explanatory variable: Corpo (g);
    - Dependent variable: Cérebro (g).

------------------------
Simple linear regression
------------------------

Call:
"lm(formula = Cérebro (g) ~ Corpo (g))"

Residuals:
    Min      1Q  Median      3Q     Max 
-2311.2  -239.7  -191.2  -101.5  4607.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.510e+02  6.987e+01   3.593 0.000478 ***
Corpo (g)   7.310e-05  7.237e-06  10.101  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 754.6 on 118 degrees of freedom
Multiple R-squared:  0.4637,    Adjusted R-squared:  0.4592 
F-statistic:   102 on 1 and 118 DF,  p-value: < 2.2e-16


Equation:
  mean[Cérebro (g)]=251.01075929 + 7.31e-05 . Corpo (g)

                 a         b Std. Error a Std. Error b t value a t value b
Linear      251.01  7.31e-05      69.8670     7.24e-06      3.59     10.10
Power         0.18  6.56e-01       0.0974     2.47e-02     -7.65     26.54
Exponential  44.60  7.75e-08       0.1791     1.85e-08     21.21      4.18
Logistic      5.02 -1.29e-07       0.1896     1.96e-08     26.50     -6.55
            Pr(>|t|) a Pr(>|t|) b R squared Adj. R squared F value   Pr(>F)
Linear        4.78e-04   1.16e-17     0.464          0.459   102.0 1.16e-17
Power         5.88e-12   1.44e-51     0.856          0.855   704.1 1.44e-51
Exponential   4.55e-42   5.71e-05     0.129          0.121    17.4 5.71e-05
Logistic      1.62e-51   1.60e-09     0.266          0.260    42.9 1.60e-09

Os dados estão muito concentrados no lado inferior-esquerdo do gráfico (difícil assumir linearidade). Há elefantes (com cérebros maiores do que 3.5 kg) e baleias (com corpos maiores que 40 toneladas) muito distantes dos demais; os pontos tornam-se mais dispersos com os corpos maiores (evidência de heterocedasticidade). Neste gráfico, a banda de confiança 95% é analítica e podemos ver boa parte dos dados está fora dela, sendo este um outro sinal de má escolha de RLS para ajustar estes dados (CorpoCerebro_RLS_problema.R).

# CorpoCerebro_RLS_problema.R
# (versao nao padronizada, para regressão)

alfa <- 0.05
B <- 0

source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.correg.R")
source("eiras.density_and_normal.R")
col_CxC <- friendlycolor(39) # cinza
pch_CxC <- 21

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))

correg(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
       alpha=alfa, B=B, 
       method="lm_robust",
       xlab="Corpo (g)", ylab="Cérebro (g)", 
       bg="transparent", col=col_CxC, pch=pch_CxC,
       suppress.text=TRUE)
# marca os problemas
points(200000,500,cex=5,lwd=4)
lines(c(-2000000,3000000),c(1500,11000),lwd=4,lty=5)
lines(c(300000,70000000),c(-50,2000),lwd=4,lty=5)

Usar versão mais robusta, com bootstrapping, não resolve bem o problema. A banda de confiança de 95% torna-se muito mais larga e inclui melhor os dados, mas a reta de regressão muda para passar onde não há dados. Portanto, não é uma reta de regressão adequada para predizer e explicar seus pontos (CorpoCerebro_RLS_bootstrapping.R):

# CorpoCerebro_RLS_bootstrapping.R

alfa <- 0.05
B <- 1e3

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

col_Corpo <- friendlycolor(14) # verde
pch_Corpo <- 22
col_Cerebro <- friendlycolor(39) # cinza
pch_Cerebro <- 23
col_CxC <- friendlycolor(39) # cinza
pch_CxC <- 21

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))

correg (CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
                     alpha=alfa, B=B, 
                     method="lm_robust",
                     xlab="Corpo (g)", ylab="Cérebro (g)", 
                     ylim=c(0,9000),
                     bg="transparent", col=col_CxC, pch=pch_CxC)

source("CorpoCerebro_RLS_bootstrapping.R")

     Variable Mean           Sd  Min. 1st Qu. Median 3rd Qu. Max.   n NA's
1   Corpo (g) 9490 9558998.2198 63400   22045  90525    5628 5454 120    0
2 Cérebro (g)  130    1026.0897   400     345    425    97.3  103 120    0
lambda =  14228.7891413801 
Assuming:
    - Explanatory variable: Corpo (g);
    - Dependent variable: Cérebro (g).

------------------------
Simple linear regression
------------------------

    SLR: 1000 resamples

data:  Cérebro (g) ~ Corpo (g)

alternative hypothesis: true slope is not equal to 0

                median          min        1st.Q        3rd.Q          max
intercept 2.426792e+02 4.217911e+01 1.992007e+02 2.895633e+02 4.701943e+02
slope     7.461084e-05 3.710690e-05 6.129055e-05 9.133222e-05 3.576278e-03
                IC95.1       IC95.2
intercept 1.192914e+02 3.964797e+02
slope     4.050237e-05 9.210436e-04

Equation:
  mean[Cérebro (g)] = 242.67920517 [119.29143779, 396.47974915] + 
        7.461e-05 [4.05e-05, 0.00092104] . Corpo (g)

Além disto, as distribuições das massas dos corpos e dos cérebros não são uninormais. Observe os gráficos (CorpoCerebro_distribuicoes.R):

# CorpoCerebro_RLS_problema.R
# (versao nao padronizada, para regressão)

source("eiras.friendlycolor.R")
source("eiras.density_and_normal.R")

col_Corpo <- friendlycolor(14) # verde
col_Cerebro <- friendlycolor(39) # cinza

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))

boxplot(CorpoCerebro$Corpo,
        xlab="Corpo (g)", ylab="densidade")
boxplot(CorpoCerebro$Cerebro,
        xlab="Cérebro (g)", ylab="densidade")

# density plots
density_and_normal(CorpoCerebro$Corpo, 
                   col=col_Corpo,
                   xlab="Corpo (g)", ylab="densidade")
density_and_normal(CorpoCerebro$Cerebro, 
                   col=col_Cerebro,
                   xlab="Cérebro (g)", ylab="densidade")

Neste exemplo, portanto, viola as suposições de uninormalidade (qualquer que seja a VD), homocedasticidade e linearidade da RLS.

Transformação potência de Tukey

A transformação potência de Tukey de variável intervalar é uma transformação não linear que altera o formato da distribuição dos dados. Sugerimos utilizar o gráfico de densidade para visualizar o formato da distribuição de variável intervalar.

Desenvolvemos uma função que seleciona a transformação potência de Tukey da variável intervalar que maximiza sua simetria: eiras.tukey.try.R.

Este método de simetrização é indicado se a amostra tem pelo menos 30 observações independentes.

Com esta função verificamos:

• para massa do corpo (mostrando todas as transformações potência):

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))

source("eiras.tukey.try.R")
potencia <- tukey.try(CorpoCerebro$Corpo,
                      xlab="Corpo (g)", 
                      show="all")

-------------------
Tukey's lambda = -2
-------------------
mean: -7.56089e-05 
st.dev.: 0.0003752577 
median: -8.002574e-08 
mode: 7.307591e-07 
iqr: 9.448253e-07 
Asymmetry coefficient: 80.79764 
    p < 1e-18

-------------------
Tukey's lambda = -1
-------------------
mean: -0.002630764 
st.dev.: 0.008322571 
median: -0.0002828739 
mode: -0.0001323697 
iqr: 0.0009277772 
Asymmetry coefficient: 2.692882 
    p < 1e-18

------------------
Tukey's lambda = 0
------------------
mean: 8.599484 
st.dev.: 2.91189 
median: 8.17056 
mode: 8.053916 
iqr: 3.067716 
Asymmetry coefficient: -0.1778416 
    p = 0.0180

------------------
Tukey's lambda = 1
------------------
mean: 1612502 
st.dev.: 9558998 
median: 3535.5 
mode: -16270.35 
iqr: 21061 
Asymmetry coefficient: -77.33594 
    p < 1e-18

------------------
Tukey's lambda = 2
------------------
mean: 9.321316e+13 
st.dev.: 6.296546e+14 
median: 12501021 
mode: -376555023 
iqr: 486862295 
Asymmetry coefficient: -191457.7 
    p < 1e-18

-------------------
Best transformation
-------------------
Tukey's lambda = 0
    (logarithm transformation)
Asymmetry coefficient = -0.17784156613391
    p = 0.0180
equation: ln[Corpo (g)]

Massa do corpo deve ser transformada por logaritmo (potência tendendo a zero).

• para massa do cérebro (mostrando só a melhor transformação potência):

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))
potencia <- tukey.try(CorpoCerebro$Cerebro, 
                      xlab="Cérebro (g)", 
                      show="final")

-------------------
Best transformation
-------------------
Tukey's lambda = 0
    (logarithm transformation)
Asymmetry coefficient = 0.225156260573833
    p = 0.0781
equation: ln[Cérebro (g)]

Massa do cérebro também deve ser transformada por logaritmo.

Regressão log-log

Então, executamos a regressão implementada em CorpoCerebro_RLS_log.R, que cria duas colunas adicionais com a transformação logarítmica:

CorpoCerebro$Corpo_log <- log(CorpoCerebro$Corpo)
CorpoCerebro$Cerebro_log <- log(CorpoCerebro$Cerebro)

obtendo:

# CorpoCerebro_RLS_log.R

alfa <- 0.05
B <- 0
source("eiras.jitter.R")
source("eiras.correg.R")

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))
# transformacao log
CorpoCerebro$Cerebro_log <- log(CorpoCerebro$Cerebro)
CorpoCerebro$Corpo_log <- log(CorpoCerebro$Corpo)

correg(CorpoCerebro$Corpo_log, CorpoCerebro$Cerebro_log,
       main="estimatr::lm_robust",
       alpha=alfa, B=B, 
       method="lm_robust",
       xlab="Ln[Corpo (g)]", ylab="Ln[Cérebro (g)]", 
       bg="transparent", col="black", pch=21)

MVN::mvn(CorpoCerebro[3:4],
         univariate_test = "SW")$univariate_normality
nml <- MVN::mvn(CorpoCerebro[3:4],
                univariate_test = "SW",
                power_family="bcPower",
                power_transform_type="rounded")
nml$univariate_normality
nml$power_transform_lambda

nml_log <- MVN::mvn(CorpoCerebro[5:6],
                    univariate_test = "SW")
nml_log$univariate_normality

fit <- lm(data=CorpoCerebro,
          Cerebro~Corpo)
shapiro.test(resid(fit))
lmtest::harvtest(fit)
lmtest::bptest(fit)
plot(fit, ask=FALSE, which = 5)
DescTools::PlotBag(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
                   xlab="Corpo (g)", ylab="Cérebro (g)", cex=0.6)
print(REAT::curvefit(CorpoCerebro$Corpo, CorpoCerebro$Cerebro, print.results=FALSE)$models_comp, digits=3)

fitlog <- lm(data=CorpoCerebro,
             Cerebro_log~Corpo_log)
summary(fitlog)
shapiro.test(resid(fitlog))
lmtest::harvtest(fitlog)
lmtest::bptest(fitlog)
plot(fitlog, ask=FALSE, which = 5)
DescTools::PlotBag(log(CorpoCerebro$Corpo_log), log(CorpoCerebro$Cerebro_log),
                   xlab="ln[Corpo (g)]", ylab="ln[Cérebro (g)]", cex=0.6)
print(REAT::curvefit(log(CorpoCerebro$Corpo+1), log(CorpoCerebro$Cerebro+1), print.results=FALSE)$models_comp, digits=3)

source("CorpoCerebro_RLS_log.R")

         Variable   Mean     Sd    Min. 1st Qu.  Median          3rd Qu.
1   Ln[Corpo (g)]  9.158 2.9119 11.0572 10.0008 11.4134 8.63550941823428
2 Ln[Cérebro (g)] 4.8675 2.0633  5.9915  5.8435  6.0521 4.57779898919196
              Max.   n NA's
1 8.60410456340553 120    0
2 4.63472898822964 120    0
lambda =  2.37998013164249 
Assuming:
    - Explanatory variable: Ln[Corpo (g)];
    - Dependent variable: Ln[Cérebro (g)].

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Ln[Cérebro (g)] ~ Ln[Corpo (g)])"

Standard error type:  HC2 

Coefficients:
              Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper  DF
(Intercept)    -1.7165    0.30900  -5.555 1.743e-07  -2.3284  -1.1046 118
Ln[Corpo (g)]   0.6558    0.03659  17.922 1.802e-35   0.5833   0.7282 118

Multiple R-squared:  0.8565 ,   Adjusted R-squared:  0.8553 
F-statistic: 321.2 on 1 and 118 DF,  p-value: < 2.2e-16

Equation:
  mean[Ln[Cérebro (g)]]=-1.71652078 + 0.65575759 . Ln[Corpo (g)]

                 a         b Std. Error a Std. Error b t value a t value b
Linear      251.01  7.31e-05      69.8670     7.24e-06      3.59     10.10
Power         0.18  6.56e-01       0.0974     2.47e-02     -7.65     26.54
Exponential  44.60  7.75e-08       0.1791     1.85e-08     21.21      4.18
Logistic      5.02 -1.29e-07       0.1896     1.96e-08     26.50     -6.55
            Pr(>|t|) a Pr(>|t|) b R squared Adj. R squared F value   Pr(>F)
Linear        4.78e-04   1.16e-17     0.464          0.459   102.0 1.16e-17
Power         5.88e-12   1.44e-51     0.856          0.855   704.1 1.44e-51
Exponential   4.55e-42   5.71e-05     0.129          0.121    17.4 5.71e-05
Logistic      1.62e-51   1.60e-09     0.266          0.260    42.9 1.60e-09

[1] "Warning: NA elements have been exchanged by median values!!"

                  a      b Std. Error a Std. Error b t value a t value b
Linear      -1.2485  0.612       0.2025       0.0223     -6.17      27.5
Power        0.0889  1.744       0.0631       0.0684    -16.66      25.5
Exponential  0.7045  0.184       0.1060       0.0117     -3.30      15.8
Logistic     3.6929 -0.408       0.1544       0.0170     23.92     -24.0
            Pr(>|t|) a Pr(>|t|) b R squared Adj. R squared F value   Pr(>F)
Linear        1.01e-08   4.45e-53     0.865          0.864     754 4.45e-53
Power         8.98e-33   8.34e-50     0.846          0.845     650 8.34e-50
Exponential   1.27e-03   7.78e-31     0.678          0.676     249 7.78e-31
Logistic      4.38e-47   3.44e-47     0.830          0.828     575 3.44e-47

Observe, no gráfico, que os dados transformados atendem melhor às suposições de normalidade, homocedasticidade e linearidade. Para \(\alpha=0.05\), RLS é estatisticamente significante (\(p<2.2 \cdot 10^{-16}\)). O tamanho de efeito é grande, \(R^2 \approx 85.7\%\). A equação da reta de regressão log-log é:

\[ \widehat{\ln(\text{cérebro})} = -1.7165208 + 0.6557576 \times \ln(\text{corpo}) \]

ou

\[ \widehat{\text{cérebro}} = \exp(-1.7165208 + 0.6557576 \times \ln(\text{corpo})) \]

ou, ainda, utilizando https://www.wolframalpha.com/, podemos encontrar que esta fórmula é equivalente a

\[ \widehat{\text{cérebro}} = 0.17969 \times {\text{corpo}}^{0.655758} \]

• ln(y) = -1.7165208 + 0.6557576 ln(x)

Interessantemente, no artigo científico publicado em 1936/1937, Gerhardt von Bonin utilizando apenas papel e paciência chegou à expressão:

\[ \begin{align} \widehat{\text{cérebro}} &= 0.18 \times {\text{corpo}}^{0.655}\\ \widehat{\text{cérebro}} &= 0.18 \times {\text{corpo}}^{2/3} \end{align} \]

Este resultado é um feito notável!

• Bonin.R

                    Modelo Inclinação IC_inferior IC_superior
1                OLS_todos     0.6558      0.6068      0.7047
2            Robusto_todos     0.6558      0.5833      0.7282
3     OLS_semHumanoCetaceo     0.6558      0.6068      0.7047
4 Robusto_semHumanoCetaceo     0.6558      0.5833      0.7282

O parâmetro de inclinação representa o expoente r no modelo E = k S^r.
Valores próximos de 0.65 indicam crescimento inelástico com elasticidade constante < 1.
A inclusão de humanos tende a elevar r (por encefalização), e a de cetáceos tende a reduzi-lo (por massa corporal inflada).
Comparar os intervalos: se o IC de r sem humanos/cetáceos inclui 0.65, o resultado reproduz von Bonin (1937).

Usando regressão linear simples robusta, o intervalo de confiança de 95% do parâmetro alométrico \(b\) é \([0.583, 0.728]\). Note que este IC95% contém o valor \(2/3=6/9=0.6\bar{6}\). Mostrando que a conjectura de Snell (1891) é plausível estatisticamente, mesmo com a inclusão de humano e cetáceos. A estimativa pontual com e sem humanos/cetáceos do parâmetro alométrico é \(b=0.656\). O IC95% de \(b\) com e sem humanos/cetáceos não inclui \(5/9 = 0.5\bar5\). Além disso, IC95% de \(b\) não contém o valor unitário, sendo inelástica a relação alométrica entre as massas.

Venditti et al. (2024) realizaram comparativo amplo entre mamíferos que mostrou que a relação cérebro-corpo não segue uma simples lei de potência universal, mas apresenta dependência curvilinear da massa corporal. Este resultado desafia a suposição de que a relação é estritamente linear no log-log.

Outro feito notável é WolframAlpha. Há muitos recursos nele. Neste exemplo, para obtermos a simplificação da equação \(\widehat{\text{cérebro}} = \exp(-1.7165208 + 0.6557576 \times \ln(\text{corpo}))\) para \(\widehat{\text{cérebro}} = 0.17969 \times {\text{corpo}}^{0.655758}\), tudo que fizemos foi incluir a equação na mesma notação que usamos em R, assim: Observe que o sistema tenta “adivinhar” nossas possíveis intenções, oferecendo a versão simplificada que transcrevemos, elaborando gráficos e encontrado propriedades sobre a equação solicitada.

Outro aspecto interessante é que, usando RLS, por causa das transformações não lineares, encontramos um ajuste de uma função não linear para a predição da massa do cérebro a partir da massa corpórea. Podemos verificar que a solução é adequada, gerando o gráfico:

CorpoCerebro <- data.frame(readxl::read_excel("CorpoCerebroVonBonin.xlsx"))

# equacao obtida
corpo_range <- seq(min(CorpoCerebro$Corpo,na.rm=TRUE), 
                   max(CorpoCerebro$Corpo,na.rm=TRUE),
                   length.out = 10000)
cerebro_estimado <- 0.183502*(corpo_range^0.651679)

# grafico
plot(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
     xlab="Corpo (g)", ylab="Cérebro (g)", 
     pch=21, col="#000088cc", bg="transparent")
lines(corpo_range,cerebro_estimado,
      col="black",lwd=2)

A predição parece funcionar bem para a maior parte dos animais, excetuando os três que têm corpos maiores (acima de \(4 \cdot 10^7 = 40000000~g = 40000~kg~~\text{ou}~~40~\text{toneladas}\)). Para saber quais são, usamos:

print(CorpoCerebro[CorpoCerebro$Corpo>=4e7,])

      Grupo                Animal Cerebro    Corpo
108 Cetacea       Megaptera boops    3531 42372000
109 Cetacea Balaenoptera rostrata    2490 62250000
110 Cetacea Balaenoptera rostrata    7000 73997000

São baleias:

• uma Megaptera boops e
• duas Balaenoptera rostrata

Para melhor observarmos o acerto da solução, podemos alterar a escala dos eixos, excluindo estes animais muito grandes que estão à direita e parecem outliers.

plot(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
     xlab="Corpo (g)", ylab="Cérebro (g)", 
     xlim=c(0,6e6),
     pch=21, col="#000088cc", bg="transparent")
lines(corpo_range,cerebro_estimado,
      col="black",lwd=2)

Quem são os três animais com cérebro acima de 3 kg?

print(CorpoCerebro[CorpoCerebro$Cerebro>3000 & 
                   CorpoCerebro$Cerebro<6000 & 
                   CorpoCerebro$Corpo<=6e6,])

       Grupo                Animal Cerebro   Corpo
93  Ungulata     Elephas africanus    4370 1642000
94  Ungulata       Elephas indicus    5045 2547000
111  Cetacea Balaenoptera musculus    3636 5090400

São os elefantes africano e indiano, e a baleia azul:

• Elephas africanus
• Elephas indicus
• Balaenoptera musculus

E qual é este com corpo acima de uma tolenada, mas cérebro abaixo de 1 kg?

print(CorpoCerebro[CorpoCerebro$Cerebro<1000 & 
                   CorpoCerebro$Corpo>=1e6,])

       Grupo                       Animal Cerebro   Corpo
100 Ungulata Hippopotamus amphibiu IS 582     582 1749730

Interessantemente, outro ungulado: o hipopótamo.

Ainda, para observar melhor o canto inferior-esquerdo onde está a maioria dos animais considerados (corpos de até 300 kg e cérebros de até 2 kg), alteramos novamente a escala do gráfico:

plot(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
     xlab="Corpo (g)", ylab="Cérebro (g)", 
     xlim=c(0,300000), ylim=c(0,2000),
     pch=21, col="#000088cc", bg="transparent")
lines(corpo_range,cerebro_estimado,
      col="black",lwd=2)

Quais são os dois animais de 70 kg com cérebro entre 1 e 2 kg? Podemos encontrá-los facilmente:

print(CorpoCerebro[CorpoCerebro$Corpo>=50000 & CorpoCerebro$Corpo<=100000 &
                   CorpoCerebro$Cerebro>=1000 & CorpoCerebro$Cerebro<=1500,])

       Grupo                      Animal Cerebro Corpo
105  Cetacea Lagenorrhynchus albirostris    1126 67560
120 Primates                Homo sapiens    1320 62000

Dois animais: um cetáceo e o ser humano. Com uma rápida pesquisa no Google, encontramos um golfinho:

• Lagenorrhynchus albirostris

Qual é o outro, com cérebro de quase 2 kg e maior peso que aparece no canto superior-direito do gráfico?

print(CorpoCerebro[CorpoCerebro$Corpo>=250000 & CorpoCerebro$Corpo<=300000 &                                 CorpoCerebro$Cerebro>=1800 & CorpoCerebro$Cerebro<=2000,])

      Grupo          Animal Cerebro  Corpo
106 Cetacea Tursiops tursio    1886 278000

Outro golfinho:

• Tursiops tursio

E podemos prosseguir mais, explorando o canto inferior-esquerdo do gráfico:

plot(CorpoCerebro$Corpo, CorpoCerebro$Cerebro,
     xlab="Corpo (g)", ylab="Cerebro (g)", 
     xlim=c(0,50000), ylim=c(0,400),
     pch=21, col="#000088cc", bg="transparent")
lines(corpo_range,cerebro_estimado,
      col="black",lwd=2)

Procurando os maiores cérebros deste trecho:

print(CorpoCerebro[CorpoCerebro$Corpo<=50000 &
                   CorpoCerebro$Cerebro>170 & CorpoCerebro$Cerebro<400,])

      Grupo                       Animal Cerebro Corpo
2  Primates Anthropopithecus troblodytes   345.0 22045
35 Primates               Papio doguerra   182.0  4990
36 Primates                 Papio anubis   212.7  6810
37 Primates                 Papio sphinx   174.5  4763
39 Primates                 Papio bahuin   179.5 10546
40 Primates                  Papio spec.   213.0 22220

Encontramos, agora, primatas:

• chimpanzé: Anthropopithecus troblodytes é o nome arcaico, hoje renomeado como Pan troglodytes

• Babuíno: Papio sp.

Não encontramos somente uma função média que relaciona massa de corpo e cérebro, mas também podemos utilizá-la para localizar animais atípicos. Cetáceos, elefantes e primatas próximos aos animais humanos são atípicos. Golfinhos têm uma relação tamanho corpo/cérebro parecida com a de nossa espécie. Em uma planilha com 120 animais, encontrar uma função geral e linear é muito interessante preditivamente e explicativamente.

Terceiro estudo: Relação entre massa encefálica e massa corporal em mamíferos por Venditti et al (2024)

O estudo de Venditti et al (2024) utilizou um conjunto de 1.504 espécies de mamíferos com medidas emparelhadas de massa encefálica e massa corporal, abrangendo toda a diversidade filogenética dos mamíferos placentários e marsupiais.

Os dados estão disponíveis em 41559_2024_2451_MOESM3_ESM.xlsx.

• RLS_correl_comp_BBM.R

Carregando pacotes exigidos: ggplot2

Carregando pacotes exigidos: RColorBrewer

Carregando pacotes exigidos: eirasagree

Anexando pacote: 'eirasagree'

O seguinte objeto é mascarado _por_ '.GlobalEnv':

    LambdaEstimate

Intervalos de confiança e predição

O IC95% mostra onde deve estar a média verdadeira da estatura para pessoas com um determinado valor de cúbito. Ele mede a incerteza sobre a média prevista pelo modelo. É estreito porque considera apenas o erro da estimativa da média.

O IP95% mostra onde deve cair uma nova observação individual com aquele cúbito. Ele é mais largo porque inclui também a variabilidade natural entre indivíduos.

Quando se usa predict com interval="confidence", obtém-se o IC95%; com interval="prediction", o IP95%. Ambos usam a distribuição t de Student e valem para cada ponto separadamente, pois são intervalos pontuais. Os intervalos pontuais com t são exatamente os produzidos por predict, HH::ci.plot e investr::plotFit.

As bandas de Working–Hotelling (WH) são diferentes: em vez de considerar um ponto de cada vez, elas garantem que toda a reta de regressão esteja dentro da faixa com 95% de confiança. São intervalos simultâneos, mais conservadores e, portanto, mais largos. Não há pacote R para bandas de WH.

Resumindo:

• predict(..., "confidence"): confiança da média (pontual)
• predict(..., "prediction"): predição individual (pontual)
• Working–Hotelling: confiança ou predição simultânea ao longo de todo o eixo x

Exemplo: Cúbito de Galton

No artigo de Galton (1888), o cúbito é usado para predizer estatura de estudante masculino inglês de 21 anos.

• ICIP95_cubito.R)

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
# Desenhar IC95% e IP95% em x0 = 46 cm, com reta e bandas (sem pontos)
# Base R; auto-contido.

# Dados e conversão
Dados_cubit <- psychTools::heights
Dados_cubit$cubit_cm  <- Dados_cubit$cubit  * 2.54
Dados_cubit$height_cm <- Dados_cubit$height * 2.54

# Regressao
fit <- stats::lm(height_cm ~ cubit_cm, data = Dados_cubit)
print(smm <- summary(fit), digits=3)
cat("\nErro-padrão da estimativa (SEE) = ", round(smm$sigma, 2), "\n", sep="")

# Bandas ao longo do eixo x
alpha  <- 0.05
x_seq  <- seq(min(Dados_cubit$cubit_cm), max(Dados_cubit$cubit_cm), length.out = 400)
nd     <- data.frame(cubit_cm = x_seq)
pred_c <- stats::predict(fit, newdata = nd, interval = "confidence",  level = 1 - alpha)
pred_p <- stats::predict(fit, newdata = nd, interval = "prediction", level = 1 - alpha)

# Grafico base: reta e bandas (sem pontos)
yl <- range(pred_p)
graphics::plot(NA, xlim = range(x_seq), ylim = yl,
               xlab = "Cúbito (cm)", ylab = "Estatura (cm)",
               main = "Reta e bandas 95%\ncúbito = 46 cm")
graphics::abline(fit, lwd = 2)
graphics::lines(x_seq, pred_c[, "lwr"], lty = 2, lwd = 2)  # IC 95% (media)
graphics::lines(x_seq, pred_c[, "upr"], lty = 2, lwd = 2)
graphics::lines(x_seq, pred_p[, "lwr"], lty = 3, lwd = 2)  # IP 95% (nova obs.)
graphics::lines(x_seq, pred_p[, "upr"], lty = 3, lwd = 2)
graphics::legend("topleft",
                 legend = c("Reta de regressão", "IC 95% (média)", "IP 95% (nova obs.)"),
                 lwd = 2, lty = c(1, 2, 3), bty = "n")

# x0 e previsoes pontuais
x0   <- 46
novo <- data.frame(cubit_cm = x0)
IC95 <- stats::predict(fit, newdata = novo, interval = "confidence",  level = 0.95)
IP95 <- stats::predict(fit, newdata = novo, interval = "prediction", level = 0.95)
yhat <- as.numeric(IC95[1]); IC_l <- as.numeric(IC95[2]); IC_u <- as.numeric(IC95[3])
IP_l <- as.numeric(IP95[2]); IP_u <- as.numeric(IP95[3])

# Linha vertical em x0
graphics::abline(v = x0, lty = 4)

# Ticks horizontais para IC e IP em torno de x0
dx <- diff(range(x_seq)) * 0.01
graphics::segments(x0 - dx, IC_l, x0 + dx, IC_l, lwd = 2, lty = 2) # IC lower
graphics::segments(x0 - dx, IC_u, x0 + dx, IC_u, lwd = 2, lty = 2) # IC upper
graphics::segments(x0 - dx, IP_l, x0 + dx, IP_l, lwd = 2, lty = 3) # IP lower
graphics::segments(x0 - dx, IP_u, x0 + dx, IP_u, lwd = 2, lty = 3) # IP upper

# Marcador da media prevista em x0
graphics::points(x0, yhat, pch = 16, cex = 1.1)

# Rotulos compactos (opcional)
graphics::text(x0 + 1.2*dx, yhat,  labels = sprintf("yhat=%.2f", yhat), pos = 4, cex = 0.9)
#graphics::text(x0 + 1.2*dx, IC_u, labels = sprintf("IC_up %.2f", IC_u), pos = 4, cex = 0.8)
#graphics::text(x0 + 1.2*dx, IC_l, labels = sprintf("IC_lo %.2f", IC_l), pos = 4, cex = 0.8)
graphics::text(x0 + 1.2*dx, IP_u, labels = sprintf("IP_up %.2f", IP_u), pos = 4, cex = 0.8)
graphics::text(x0 + 1.2*dx, IP_l, labels = sprintf("IP_lo %.2f", IP_l), pos = 4, cex = 0.8)

# Saida numerica
cat(sprintf("x0 = %.2f cm | yhat = %.2f cm | IC95 = [%.2f, %.2f] cm | IP95 = [%.2f, %.2f] cm\n",
            x0, yhat, IC_l, IC_u, IP_l, IP_u))

source("ICIP95_cubito.R")

Call:
stats::lm(formula = height_cm ~ cubit_cm, data = Dados_cubit)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.242  -2.926  -0.386   3.303  12.662 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   65.857      4.900    13.4   <2e-16 ***
cubit_cm       2.274      0.107    21.4   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.94 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


Erro-padrão da estimativa (SEE) = 3.94

x0 = 46.00 cm | yhat = 170.45 cm | IC95 = [170.04, 170.87] cm | IP95 = [162.70, 178.21] cm

• ICIP95_WH_cubito.R)

x0=46.00 cm | yhat=170.45 cm
IC95% t (pontual):   [170.04, 170.87] cm
IC95% WH (simult.):  [169.94, 170.97] cm
IP95% t (pontual):   [162.70, 178.21] cm
IP95% WH (simult.):  [160.76, 180.15] cm

Diferença entre correlação e regressão

A correlação de Pearson é sempre confundida com RLS. A confusão é natural, por causa de suas similaridades: ambas tratam de relações entre duas variáveis intervalares (quantitativas), relação esta que deve ser linear; calculam, respectivamente, \(r\) e \(R^2\), sendo \(r^2=R^2\); são testáveis estatisticamente para \(H_0: \rho=0\) e \(H_0: \beta_1=0\) chegando à mesma conclusão (variáveis não correlacionadas, i.e., \(\rho=0\), correspondem à reta de regressão horizontal, i.e., \(\beta_1=0\)), com o mesmo valor p quando computamos com os métodos clássicos implementados em cor.test e lm.

É necessário, então, ter clareza sobre as diferenças, além da finalidade propriamente dita: correlação de Pearson mede o grau de linearidade entre duas variáveis, sem direção e sem dimensão, enquanto RLS pretende obter uma equação para, a partir do valor de uma variável explicativa (VE), explicar e estimar o valor médio de uma variável de desfecho (VD), nesta direção e com dimensão.

Ser adimensional ou dimensional neste contexto corresponde, respectivamente, a tratar com números puros ou com unidades de medida. Como mostramos acima, a correlação de Pearson é computada com as variáveis padronizadas. Como o procedimento de padronização elimina as unidades de medida, dizer que a correlação de Pearson utiliza as variáveis padronizadas implica que ela é adimensional. RLS, ao contrário, busca explicação e previsão e, portanto, precisa manter as unidades de medida das variáveis originais, i.e., é dimensional.

Volte a observar os gráficos da seção valores do r de Pearson. Mostramos exemplos de nuvens de pontos com \(-1 \le r \le 1\). Repare que não procuramos traçar retas de regressão entre os pontos; as nuvens esboçam elipses, mais estreitas quanto mais próximos os valores de r se aproximam de \(-1\) ou \(1\) e com aspecto de circunferência quando \(r=0\). Todas as elipses têm eixo principal com inclinação próxima a 45^o em relação ao eixo das abscissas, pois este é o efeito causado pela padronização. Caso alguém tentasse imaginar uma reta que passasse entre os pontos destas nuvens, provavelmente traçaria a reta data por \(\hat{y}=x\), i.e., com \(\hat{\beta}_1=1, \hat{\beta}_0=0\), correspondendo à bissetriz do primeiro quadrante do plano cartesiano, inclinada em 45^o.

No entanto, na seção padronização das variáveis encontramos para a RLS com as variáveis padronizadas os valores \(\hat{\beta}_1 = r\) (a inclinação da reta é a correlação de Pearson) e \(\hat{\beta}_0 = 0\) (o intercepto é nula). Não parece estranho? As retas de regressão calculadas para as nuvens, com exceção de \(r=-1\) e \(r=1\), não têm inclinação correspondente ao eixo principal das nuvens de pontos que observamos, sempre com ângulos de 45^o. Que retas e que ângulos são estes? O entendimento destes comportamentos da correlação de Pearson e da RLS é o que mostra mais claramente suas diferenças.

Para entender o que ocorre, vamos usar os dados da planilha Adm2008.xlsx. Entre outras variáveis, esta planilha tem a estatura e a massa corpórea relatadas por 89 estudantes masculinos e femininos de uma turma de graduação do curso de administração da FEA-USP em 2008.

Uma relação linear pode ser considerada para estudantes do sexo masculino?

Executando Adm2008_descritiva.R:

# Biometria_rPearson.R
source("eiras.friendlycolor.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
col <- friendlycolor(7) # azul
pch <- 24

Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosM <- subset(Dados, Sexo=="Masculino")
lst <- correg(DadosM$Estatura, DadosM$MCT,
              method="preview",
              alpha=alfa, B=B, 
              xlab="EstaturaM (m)", ylab="MassaM (kg)", 
              col=col, bg="transparent", pch=pch)

source("Adm2008_descritiva.R")

       Variable Mean      Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 EstaturaM (m) 1.72  0.0808 1.88    1.91   1.72    1.82 1.88 51    0
2   MassaM (kg)   82 12.0985   82      83     58      77   82 51    0
lambda =  0.00428473520865842 
Assuming:
    - Explanatory variable: EstaturaM (m);
    - Dependent variable: MassaM (kg).

Admitindo relação linear,

• Correlação de Pearson (implementada em Adm2008_rPearson.R) é:

# Biometria_rPearson.R
source("eiras.friendlycolor.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
col <- friendlycolor(7) # azul
pch <- 24

Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosM <- subset(Dados, Sexo=="Masculino")
lst <- correg(DadosM$Estatura, DadosM$MCT,
              method="pearson",
              alpha=alfa, B=B, 
              xlab="EstaturaM (m)", ylab="MassaM (kg)", 
              col=col, bg="transparent", pch=pch)

source("Adm2008_rPearson.R")

       Variable Mean      Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 EstaturaM (m) 1.72  0.0808 1.88    1.91   1.72    1.82 1.88 51    0
2   MassaM (kg)   82 12.0985   82      83     58      77   82 51    0
lambda =  0.00428473520865842 
Assuming:
    - Explanatory variable: EstaturaM (m);
    - Dependent variable: MassaM (kg).

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

    Pearson's product-moment correlation

data:  EstaturaM (m) and MassaM (kg)
t = 5.925, df = 49, p-value = 3.054e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4507278 0.7823524
sample estimates:
     cor 
0.646062 

para o qual computamos \(r=\) 0.65.

• RLS (implementada em Adm2008_RLS.R) é:

# Biometria_rPearson.R
source("eiras.friendlycolor.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
col <- friendlycolor(7) # azul
pch <- 24

Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosM <- subset(Dados, Sexo=="Masculino")
lst <- correg(DadosM$Estatura, DadosM$MCT,
              method="lm_robust",
              alpha=alfa, B=B, 
              xlab="EstaturaM (cm)", ylab="MassaM (kg)", 
              col=col, bg="transparent", pch=pch)

source("Adm2008_RLS.R")

        Variable Mean      Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 EstaturaM (cm) 1.72  0.0808 1.88    1.91   1.72    1.82 1.88 51    0
2    MassaM (kg)   82 12.0985   82      83     58      77   82 51    0
lambda =  0.00428473520865842 
Assuming:
    - Explanatory variable: EstaturaM (cm);
    - Dependent variable: MassaM (kg).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = MassaM (kg) ~ EstaturaM (cm))"

Standard error type:  HC2 

Coefficients:
               Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)      -96.20      19.61  -4.906 1.068e-05  -135.61   -56.79 49
EstaturaM (cm)    96.78      11.13   8.698 1.669e-11    74.42   119.14 49

Multiple R-squared:  0.4174 ,   Adjusted R-squared:  0.4055 
F-statistic: 75.66 on 1 and 49 DF,  p-value: 1.669e-11

Equation:
  mean[MassaM (kg)]=-96.19878084 + 96.77716058 . EstaturaM (cm)

Estatura e massa corpórea têm unidades de medida, respectivamente cm e kg.

RLS encontrou a equação:

\(\widehat{\text{massa}}\, \text{(kg)} =\) -96.2 \(+\) 96.78 \(~ \times ~ \text{estatura (cm)}\)

• e RLS com as variáveis padronizadas (implementada em Adm2008_RLS_z.R) é:

# Adm2008_RLS_z.R
source("eiras.friendlycolor.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
col <- friendlycolor(7) # azul
pch <- 24

Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosM <- subset(Dados, Sexo=="Masculino")
lst <- correg(DadosM$Estatura, DadosM$MCT,
              method="lm_robust",
              standardize=TRUE,
              alpha=alfa, B=B, 
              xlab="EstaturaM (m)", ylab="MassaM (kg)", 
              col=col, bg="transparent", pch=pch)

source("Adm2008_RLS_z.R")

       Variable   Mean Sd   Min. 1st Qu.  Median           3rd Qu.
1 EstaturaM (m) -0.539  1 1.4421  1.8135  -0.539 0.699179733978704
2   MassaM (kg) 0.6223  1 0.6223   0.705 -1.3614  0.20906791549181
               Max.  n NA's
1  1.44205820133107 51    0
2 0.622341701929108 51    0
lambda =  2.24303391346275 
Assuming:
    - Explanatory variable: EstaturaM (m);
    - Dependent variable: MassaM (kg).

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = MassaM (kg) ~ EstaturaM (m))"

Standard error type:  HC2 

Coefficients:
                Estimate Std. Error    t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)   -6.033e-16    0.10739 -5.618e-15 1.000e+00  -0.2158   0.2158 49
EstaturaM (m)  6.461e-01    0.07428  8.698e+00 1.669e-11   0.4968   0.7953 49

Multiple R-squared:  0.4174 ,   Adjusted R-squared:  0.4055 
F-statistic: 75.66 on 1 and 49 DF,  p-value: 1.669e-11

Equation:
  mean[MassaM (kg)]=0 + 0.64606202 . EstaturaM (m)

encontrando

\(\widehat{\text{zmassa}} =\) -6.033e-16 \(+\) 0.65 \(~ \times ~ \text{zestatura}\)

Exceto por arredondamento, a reta calculada comprova o que esperávamos da teoria: é praticamente igual a

\[ \widehat{\text{zmassa}} = 0 + r \times \text{zestatura} \]

Vamos exibir os dados padronizados com os eixos utilizando escalas idênticas:

Note que a correlação dos pontos é próxima a \(r=\) 0.65, vista na seção valores do r de Pearson. A orientação de uma reta imaginária que passa pelo meio da nuvem de pontos poderia ser próxima à bissetriz (inclinação igual a \(1\)), como nos gráficos que observamos na seção de padronização das variáveis. No entanto, com \(r=\) 0.65 a reta de regressão desvia da bissetriz, como teorizamos.

Para verificar como a reta de regressão se relaciona com a correlação de Pearson, adicionamos ao código anterior a RLS dos dados gerados e a plotamos:

Para melhorar a visualização, podemos adicionar a elipse de predição de 68% (poderia ser qualquer outra) com seu eixo principal e duas verticais pontilhadas em seus limites, à esquerda e à direita:

É pelos pontos em que estas verticais tocam a elipse que passa a reta de regressão computada pelo modelo lm_robust(y~x).

Para ser mais completo, como a correlação de Pearson é adirecional, poderíamos computar uma segunda reta de regressão, invertendo as variáveis, com lm_robust(x~y). Caso queira, pode experimentar com demo_r.R. Este programa permite gerar distribuições artificiais com correlação próxima ao valor de r desejado e o número de pontos \(n\), para observar da RLS, mostrando as duas retas de regressão y~x e x~y tocando os pontos extremos da elipse de predição de 68%.

Comparação das relações de duas variáveis em duas condições independentes

Observe o que ocorre com os dados de Adm2008.xlsx quando exibimos separadamente as observações de estatura e massa corpórea total (MCT) de homens e mulheres obtidos entre estudantes do curso de graduação de Administração da FEA-USP.

As elipses externas de predição de 99.9% mostram que não há outliers bidimensionais se considerarmos que as distribuições dos grupos são binormais.

Pelas localização e inclinação das elipses, parece que homem é mais alto e mais pesado que mulher, e também são maiores suas variabilidades em estatura e em MCT.

A regressão e a correlaçao entre estatura e MCT são as mesmas em mulher e homem?

Comparação de correlações de Pearson

A correlação entre estatura e MCT é diferente para os dois sexos?

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

Podemos, para responder a esta questão, comparar as correlações de Pearson. Observar que os intervalos de confiança de 95% com correção de Bonferroni das correlações de Pearson estão bastante sobrepostos.

Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosF <- subset(Dados, Sexo=="Feminino")
DadosM <- subset(Dados, Sexo=="Masculino")
alpha <- 0.05
g <- nlevels(Dados$Sexo)
alphaBonf <- alpha/g
print(rF <- cor.test(DadosF$Estatura,DadosF$MCT,
                     conf.level=1-alphaBonf), 
      digits=2)

    Pearson's product-moment correlation

data:  DadosF$Estatura and DadosF$MCT
t = 5, df = 36, p-value = 2e-05
alternative hypothesis: true correlation is not equal to 0
97.5 percent confidence interval:
 0.36 0.81
sample estimates:
 cor 
0.64 

print(rM <- cor.test(DadosM$Estatura,DadosM$MCT,
                     conf.level=1-alphaBonf), 
      digits=2)

    Pearson's product-moment correlation

data:  DadosM$Estatura and DadosM$MCT
t = 6, df = 49, p-value = 3e-07
alternative hypothesis: true correlation is not equal to 0
97.5 percent confidence interval:
 0.42 0.80
sample estimates:
 cor 
0.65 

Graficamente, é mais coerente expressar com as variáveis padronizadas (lembre: não é apropriado traçar retas de regressão para a correlação). Assim, observamos:

As elipses de predição de 68% e retas de regressão com as variáveis padronizadas são, ambas, quase coincidentes.

O teste de igualdade de correlações de Pearson é apresentado a seguir:

Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)
DadosF <- subset(Dados, Sexo=="Feminino")
DadosM <- subset(Dados, Sexo=="Masculino")
print(nF <- nrow(DadosF))

[1] 38

print(nM <- nrow(DadosM))

[1] 51

print(rF <- cor(DadosF$Estatura,DadosF$MCT), digits=3)

[1] 0.637

print(rM <- cor(DadosM$Estatura,DadosM$MCT), digits=3)

[1] 0.646

tr <- psych::r.test(n=nF, r12=rF, n2=nM, r34=rM)
print(tr, digits=4)

Correlation tests 
Call:psych::r.test(n = nF, r12 = rF, r34 = rM, n2 = nM)
Test of difference between two independent correlations 
 z value 0.0682    with probability  0.9457

Comparação das retas de regressão

Também é possível testar se as duas regressões populacionais são iguais, i.e., paralelas e coincidentes.

Primeiramente, concluímos que a diferença das correlações de Pearson entre estatura e MCT para mulher e homem é não significante.

Para a comparação das duas retas de regressão, há dois testes hierárquicos. O primeiro testa se as retas são paralelas (inclinações iguais) e o segundo, caso sejam paralelas, verifica se são coincidentes (interceptos iguais). Como são dois testes, adotamos \(\alpha=0.025\) (correção de Bonferroni). Assim, não rejeitamos a hipótese nula de que são paralelas, mas os interceptos diferem e, portanto, as retas não são coincidentes.

Graficamente, podemos acomodar retas paralelas dentro das respectivas bandas de confiança de 95% com correção de Bonferroni no domínio comum aos dois sexos:

As bandas de confiança de 95% com correção de Bonferroni delimitam a região onde a reta de regressão populacional pode existir. Portanto, a possibilidade de haver pelo menos este par de retas significa que a hipótese nula (sempre populacional) de paralelismo da regressão para homem e mulher não é rejeitada. No entanto, nenhum par de retas paralelas com interceptos iguais podem ser acomodada neste caso, e a igualdade populacional dos interceptos (consequentemente a coincidência completa das retas populacionais) é rejeitada.

Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Sexo <- factor(Dados$Sexo)

## Método 1
cat("\nTeste de paralelismo (igualdade de inclinações): valor-p do efeito de interaçao Estatura:Sexo > 0.025\n")

Teste de paralelismo (igualdade de inclinações): valor-p do efeito de interaçao Estatura:Sexo > 0.025

rls_paralel <- lm(MCT ~ Estatura*Sexo,
                  data=Dados)
print(anova(rls_paralel), digits=3)

Analysis of Variance Table

Response: MCT
              Df Sum Sq Mean Sq F value Pr(>F)    
Estatura       1  10808   10808  179.93 <2e-16 ***
Sexo           1   1224    1224   20.38  2e-05 ***
Estatura:Sexo  1    169     169    2.81  0.097 .  
Residuals     85   5106      60                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nTeste de coincidência (igualdade de interceptos): valor-p do fator Sexo > 0.025\n")

Teste de coincidência (igualdade de interceptos): valor-p do fator Sexo > 0.025

rls_coincid <- lm(MCT ~ Estatura+Sexo,
                  data=Dados)
print(anova(rls_coincid), digits=3)

Analysis of Variance Table

Response: MCT
          Df Sum Sq Mean Sq F value  Pr(>F)    
Estatura   1  10808   10808     176 < 2e-16 ***
Sexo       1   1224    1224      20 2.4e-05 ***
Residuals 86   5274      61                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Método 2
cat("\nTeste de paralelismo (igualdade de inclinações):\nvalor-p do efeito de interaçao Estatura:Sexo > 0.025\n")

Teste de paralelismo (igualdade de inclinações):
valor-p do efeito de interaçao Estatura:Sexo > 0.025

m_full <- lm(MCT ~ Estatura*Sexo,
             data=Dados)
m_part  <- lm(MCT ~ Estatura+Sexo,
              data=Dados)   # sem interação
print(anova(m_part, m_full), digits=4)

Analysis of Variance Table

Model 1: MCT ~ Estatura + Sexo
Model 2: MCT ~ Estatura * Sexo
  Res.Df  RSS Df Sum of Sq    F Pr(>F)  
1     86 5274                           
2     85 5106  1     168.8 2.81 0.0973 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nTeste de coincidência (igualdade de interceptos):\nvalor-p do fator Sexo > 0.025\n")

Teste de coincidência (igualdade de interceptos):
valor-p do fator Sexo > 0.025

m_full <- lm(MCT ~ Estatura+Sexo,
             data=Dados)  # interceptos por grupo
m_part   <- lm(MCT ~ Estatura,
               data=Dados)            # intercepto comum
print(anova(m_part, m_full), digits=4)

Analysis of Variance Table

Model 1: MCT ~ Estatura
Model 2: MCT ~ Estatura + Sexo
  Res.Df  RSS Df Sum of Sq     F  Pr(>F)    
1     87 6498                               
2     86 5274  1      1224 19.96 2.4e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note que os testes de paralelismo e coincidência implementam os modelos completos (full) com três e duas variáveis independentes. No teste de paralelismo, as fórmulas do modelo linear geral (GLM) são equivalentes:

• MCT ~ Estatura*Sexo
• MCT ~ Estatura + Sexo + Estatura:Sexo

Portanto, este GLM combina os efeitos de fator (variável nominal) e covariável (variável intervalar). Há dois efeitos principais (covariável Estatura e fator Sexo) e um efeito de interação entre covariável e fator.

Implementamos, com testes estatísticos apropriados em demo_estatmct.R, resultando:

source("demo_estatmct.R")

-------------------
Paralelismo das RLS
-------------------
H0: retas populacionais com inclinações iguais
    p = 0.0973
--------------------
Coincidência das RLS
--------------------
H0: retas populacionais com interceptos iguais
    p = 0
----------
Correlação
----------
H0: correlações de Pearson populacionais iguais
    p = 0.9457
---
RLS
---

----------
- feminino
----------

Call:
lm(formula = df_admF$MCT_z ~ df_admF$Estatura_z)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.88651 -0.47436  0.02501  0.39027  1.46502 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        1.407e-15  1.268e-01    0.00        1    
df_admF$Estatura_z 6.371e-01  1.285e-01    4.96  1.7e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7814 on 36 degrees of freedom
Multiple R-squared:  0.406, Adjusted R-squared:  0.3895 
F-statistic:  24.6 on 1 and 36 DF,  p-value: 1.698e-05


Equation: media[MCT(kg)] = 0 + 0.6371 * Estatura(m)

-----------
- masculino
-----------

Call:
lm(formula = df_admM$MCT_z ~ df_admM$Estatura_z)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1519 -0.4866 -0.2106  0.2866  2.0717 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -6.079e-16  1.080e-01   0.000        1    
df_admM$Estatura_z  6.461e-01  1.090e-01   5.925 3.05e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.771 on 49 degrees of freedom
Multiple R-squared:  0.4174,    Adjusted R-squared:  0.4055 
F-statistic: 35.11 on 1 and 49 DF,  p-value: 3.054e-07


Equation: media[MCT(kg)] = 0 + 0.6461 * Estatura(m)

----------
Correlação
----------

----------
- feminino
----------

    Pearson's product-moment correlation

data:  df_admF$Estatura and df_admF$MCT
t = 4.96, df = 36, p-value = 1.698e-05
alternative hypothesis: true correlation is not equal to 0
97.5 percent confidence interval:
 0.3579143 0.8117791
sample estimates:
     cor 
0.637148 


-----------
- masculino
-----------

    Pearson's product-moment correlation

data:  df_admM$Estatura and df_admM$MCT
t = 5.925, df = 49, p-value = 3.054e-07
alternative hypothesis: true correlation is not equal to 0
97.5 percent confidence interval:
 0.4177730 0.7976172
sample estimates:
     cor 
0.646062 

Comparação de três retas de regressão e correlações com iris

Comparar, no arquivo iris, as três retas de regressão e correlações de comprimento da sépala em função da largura da sépala para as espécies setosa, versicolor e virginica. Verificar primeiramente o paralelismo das retas por meio da comparação entre um modelo com interação entre largura da sépala e espécie e um modelo sem interação. Se o paralelismo não for rejeitado, testar a igualdade dos interceptos sob o modelo com retas paralelas. Apresentar bandas de confiança simultâneas para cada reta com correção de Bonferroni considerando três grupos, adotar nível de significância de 5% e restringir a interpretação ao domínio observado da largura da sépala.

No código a seguir é mostrado que as hipóteses nulas de igualdade de duas correlações de Pearson entre largura e comprimento de sépala para as três combinações de espécies de íris não são rejeitadas para \(alpha=5\%\) com correção de Bonferroni.

A hipótese nula omnibus de igualdade das três inclinações (paralelismo) não é rejeita para \(\alpha=5\%\).

Supondo paralelismo, a hipótese nula omnibus de igualdade dos três interceptos (coincidência) é rejeita para \(\alpha=5\%\).

Supondo paralelismo, as hipóteses nulas post hoc de igualdade de dois interceptos (coincidência) nas três combinações de espécies de íris são rejeitas para \(\alpha=5\%\) com correção de Bonferroni.

Implementamos, com testes estatísticos apropriados em RLS_correl_comp_iris.R, resultando:

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")]))
print(psych::describeBy(I3[, c("Sepal.Length","Sepal.Width")],
                        group=I3$Species, mat=1, digits=2))

## ===== 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")
print(TAB.ic, digits=3, row.names=FALSE)

xlim_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)

# bandas de confiança (Bonferroni) por grupo no ggplot
g <- nlevels(I3$Species)                # = 3
alfa.Bonf <- alfa/g

p <- ggplot2::ggplot(I3, ggplot2::aes(x=Sepal.Width, y=Sepal.Length,
                    group=Species, linetype=Species,
                    fill=Species, shape=Species, color=Species)) +
  ggplot2::geom_point() +
  ggplot2::scale_color_manual(values=pal) +
  ggplot2::scale_fill_manual(values=pal) +
  #ggplot2::geom_smooth(method = stats::lm,
  #ggplot2::geom_smooth(method = MASS::rlm,
  ggplot2::geom_smooth(method = estimatr::lm_robust,
                       formula = y ~ x, 
                       na.rm = TRUE,
                       color = "black", 
                       ggplot2::aes(fill = Species),
                       level = 1 - alfa.Bonf) +
  ggplot2::labs(title = "Bandas de confiança Bonferroni (95%)",
       x = "Largura da sépala (cm)",
       y = "Comprimento da sépala (cm)") +
  ggplot2::theme_bw() +
  ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
        panel.grid.minor = ggplot2::element_blank(),
        panel.background  = ggplot2::element_blank(),
        plot.background   = ggplot2::element_rect(fill="white", color=NA),
        axis.line         = ggplot2::element_line(color="black"))
print(p)

invisible(eirasagree::ConfidenceBand(IV=I.S$Sepal.Width, 
                                     DV=I.S$Sepal.Length,
                                     method="lm",
                                     alpha = alfa.Bonf,
                                     plot = TRUE,
                                     add = FALSE,
                                     xlim = c(2,4.5),
                                     ylim = c(4,8),
                                     xlab = "Largura da sépala (cm)",
                                     ylab = "Comprimento da sépala (cm)"))
invisible(eirasagree::ConfidenceBand(IV=I.Vc$Sepal.Width, 
                                     DV=I.Vc$Sepal.Length,
                                     method="lm",
                                     alpha = alfa.Bonf,
                                     plot = TRUE,
                                     add = TRUE,
                                     xlab = "Largura da sépala (cm)",
                                     ylab = "Comprimento da sépala (cm)"))
invisible(eirasagree::ConfidenceBand(IV=I.Vg$Sepal.Width, 
                                     DV=I.Vg$Sepal.Length,
                                     method="lm",
                                     alpha = alfa.Bonf,
                                     add = TRUE,
                                     xlab = "Largura da sépala (cm)",
                                     ylab = "Comprimento da sépala (cm)"))

# regressões centradas por grupo (relatório)
regS  <- lm(Sepal.Length ~ scale(Sepal.Width, scale=FALSE), data=I.S)
regVc <- lm(Sepal.Length ~ scale(Sepal.Width, scale=FALSE), data=I.Vc)
regVg <- lm(Sepal.Length ~ scale(Sepal.Width, scale=FALSE), data=I.Vg)
print(summary(regS))
print(summary(regVc))
print(summary(regVg))

# elipses 95% dos coeficientes (intercepto × inclinação) por grupo
car::confidenceEllipse(lm(Sepal.Length ~ Sepal.Width, data=I.S),
                       main="Elipses de Confiança de 95% (Bonferroni)\nsetosa(verde) versicolor(laranja) virginica(roxo)",
                       xlim=c(0,7), ylim=c(0,2),
                       col=pal[1], levels=1-alfa.Bonf,
                       Scheffe=FALSE, center.cex=1, grid=FALSE,
                       fill=TRUE, fill.alpha=0.05)
car::confidenceEllipse(lm(Sepal.Length ~ Sepal.Width, data=I.Vc),
                       col=pal[2], levels=1-alfa.Bonf,
                       Scheffe=FALSE, center.cex=1, grid=FALSE,
                       fill=TRUE, fill.alpha=0.05, add=TRUE)
car::confidenceEllipse(lm(Sepal.Length ~ Sepal.Width, data=I.Vg),
                       col=pal[3], levels=1-alfa.Bonf,
                       Scheffe=FALSE, center.cex=1, grid=FALSE,
                       fill=TRUE, fill.alpha=0.05, add=TRUE)
abline(h=0, v=0, lty=2)

cat("\n# Teste de paralelismo (igualdade de inclinações entre 3 grupos)\n")
m_full <- lm(Sepal.Length ~ Sepal.Width * Species, data=I3)
m_par  <- lm(Sepal.Length ~ Sepal.Width + Species, data=I3)   # sem interação
#print(summary(m_full), digits=3)
print(anova(m_par, m_full), digits=3)

cat("\n# Interceptos iguais supondo paralelismo (3 grupos)\n")
m_int_diff <- lm(Sepal.Length ~ Sepal.Width + Species, data=I3)  
m_int_eq   <- lm(Sepal.Length ~ Sepal.Width, data=I3)            
print(anova(m_int_eq, m_int_diff), digits=3)

I.VcVg <- droplevels(subset(I3, Species!="setosa"))
I.VgS <- droplevels(subset(I3, Species!="versicolor"))
I.VcS <- droplevels(subset(I3, Species!="virginica"))

m_int_diff <- lm(Sepal.Length ~ Sepal.Width + Species, data=I.VcVg)  
m_int_eq   <- lm(Sepal.Length ~ Sepal.Width, data=I.VcVg)            
print(anova(m_int_eq, m_int_diff), digits=3)

m_int_diff <- lm(Sepal.Length ~ Sepal.Width + Species, data=I.VgS)  
m_int_eq   <- lm(Sepal.Length ~ Sepal.Width, data=I.VgS)            
print(anova(m_int_eq, m_int_diff), digits=3)

m_int_diff <- lm(Sepal.Length ~ Sepal.Width + Species, data=I.VcS)  
m_int_eq   <- lm(Sepal.Length ~ Sepal.Width, data=I.VcS)            
print(anova(m_int_eq, m_int_diff), digits=3)

cat("\n# Correlações por grupo e testes pareados (independentes)\n")
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="")
cat("cor(versicolor)= ", round(rVc,3), "\n", sep="")
cat("cor(virginica) = ", round(rVg,3), "\n", sep="")
cat("\nDiferenças de correlação (teste de duas correlações independentes):\n")
print(psych::r.test(n=nS , r12=rS , n2=nVc, r34=rVc), digits=5)  # setosa vs versicolor
cat("\n\n")
print(psych::r.test(n=nS , r12=rS , n2=nVg, r34=rVg), digits=5)  # setosa vs virginica
cat("\n\n")
print(psych::r.test(n=nVc, r12=rVc, n2=nVg, r34=rVg), digits=5)  # versicolor vs virginica

source("RLS_correl_comp_iris.R")

  Sepal.Length    Sepal.Width          Species  
 Min.   :4.300   Min.   :2.000   setosa    :50  
 1st Qu.:5.100   1st Qu.:2.800   versicolor:50  
 Median :5.800   Median :3.000   virginica :50  
 Mean   :5.843   Mean   :3.057                  
 3rd Qu.:6.400   3rd Qu.:3.300                  
 Max.   :7.900   Max.   :4.400                  
              item     group1 vars  n mean   sd median trimmed  mad min max
Sepal.Length1    1     setosa    1 50 5.01 0.35    5.0    5.00 0.30 4.3 5.8
Sepal.Length2    2 versicolor    1 50 5.94 0.52    5.9    5.94 0.52 4.9 7.0
Sepal.Length3    3  virginica    1 50 6.59 0.64    6.5    6.57 0.59 4.9 7.9
Sepal.Width1     4     setosa    2 50 3.43 0.38    3.4    3.42 0.37 2.3 4.4
Sepal.Width2     5 versicolor    2 50 2.77 0.31    2.8    2.78 0.30 2.0 3.4
Sepal.Width3     6  virginica    2 50 2.97 0.32    3.0    2.96 0.30 2.2 3.8
              range  skew kurtosis   se
Sepal.Length1   1.5  0.11    -0.45 0.05
Sepal.Length2   2.1  0.10    -0.69 0.07
Sepal.Length3   3.0  0.11    -0.20 0.09
Sepal.Width1    2.1  0.04     0.60 0.05
Sepal.Width2    1.4 -0.34    -0.55 0.04
Sepal.Width3    1.6  0.34     0.38 0.05

IC95% Bonferroni de correlação de Pearson
    Especie  n     r    lo    hi
     setosa 50 0.743 0.542 0.863
 versicolor 50 0.526 0.231 0.732
  virginica 50 0.457 0.144 0.687

Call:
lm(formula = Sepal.Length ~ scale(Sepal.Width, scale = FALSE), 
    data = I.S)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.52476 -0.16286  0.02166  0.13833  0.44428 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        5.00600    0.03373 148.392  < 2e-16 ***
scale(Sepal.Width, scale = FALSE)  0.69049    0.08990   7.681 6.71e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2385 on 48 degrees of freedom
Multiple R-squared:  0.5514,    Adjusted R-squared:  0.542 
F-statistic: 58.99 on 1 and 48 DF,  p-value: 6.71e-10


Call:
lm(formula = Sepal.Length ~ scale(Sepal.Width, scale = FALSE), 
    data = I.Vc)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.73497 -0.28556 -0.07544  0.43666  0.83805 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        5.93600    0.06273  94.627  < 2e-16 ***
scale(Sepal.Width, scale = FALSE)  0.86508    0.20194   4.284 8.77e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4436 on 48 degrees of freedom
Multiple R-squared:  0.2766,    Adjusted R-squared:  0.2615 
F-statistic: 18.35 on 1 and 48 DF,  p-value: 8.772e-05


Call:
lm(formula = Sepal.Length ~ scale(Sepal.Width, scale = FALSE), 
    data = I.Vg)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.26067 -0.36921 -0.03606  0.19841  1.44917 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        6.58800    0.08081  81.529  < 2e-16 ***
scale(Sepal.Width, scale = FALSE)  0.90153    0.25311   3.562 0.000843 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5714 on 48 degrees of freedom
Multiple R-squared:  0.2091,    Adjusted R-squared:  0.1926 
F-statistic: 12.69 on 1 and 48 DF,  p-value: 0.0008435

# Teste de paralelismo (igualdade de inclinações entre 3 grupos)
Analysis of Variance Table

Model 1: Sepal.Length ~ Sepal.Width + Species
Model 2: Sepal.Length ~ Sepal.Width * Species
  Res.Df  RSS Df Sum of Sq    F Pr(>F)
1    146 28.0                         
2    144 27.9  2     0.157 0.41   0.67

# Interceptos iguais supondo paralelismo (3 grupos)
Analysis of Variance Table

Model 1: Sepal.Length ~ Sepal.Width
Model 2: Sepal.Length ~ Sepal.Width + Species
  Res.Df RSS Df Sum of Sq   F Pr(>F)    
1    148 101                            
2    146  28  2      72.8 190 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table

Model 1: Sepal.Length ~ Sepal.Width
Model 2: Sepal.Length ~ Sepal.Width + Species
  Res.Df  RSS Df Sum of Sq    F  Pr(>F)    
1     98 30.1                              
2     97 25.1  1      5.03 19.4 2.7e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table

Model 1: Sepal.Length ~ Sepal.Width
Model 2: Sepal.Length ~ Sepal.Width + Species
  Res.Df  RSS Df Sum of Sq   F Pr(>F)    
1     98 84.3                            
2     97 18.5  1      65.8 344 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table

Model 1: Sepal.Length ~ Sepal.Width
Model 2: Sepal.Length ~ Sepal.Width + Species
  Res.Df  RSS Df Sum of Sq   F Pr(>F)    
1     98 39.0                            
2     97 12.3  1      26.8 212 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Correlações por grupo e testes pareados (independentes)
cor(setosa)    = 0.743
cor(versicolor)= 0.526
cor(virginica) = 0.457

Diferenças de correlação (teste de duas correlações independentes):
Correlation tests 
Call:psych::r.test(n = nS, r12 = rS, r34 = rVc, n2 = nVc)
Test of difference between two independent correlations 
 z value 1.80167    with probability  0.0716

Correlation tests 
Call:psych::r.test(n = nS, r12 = rS, r34 = rVg, n2 = nVg)
Test of difference between two independent correlations 
 z value 2.24124    with probability  0.02501

Correlation 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.66025

Elasticidade da regressão

Elasticidade, \(\eta\), é definida como a variação percentual média da VD decorrente da variação de 1% da VE. No caso da reta de regressão, \(\hat{y} = a + b\,x\), corresponde a:

\[ \eta(x) = \dfrac{b\,x}{a + b\,x} \] Se \(\eta(x)=1\) a função é isoelástica em \(x\).

Se \(\eta(x)<1\) a função é inelástica em \(x\).

Se \(\eta(x)>1\) a função é elástica em \(x\).

• elasticidade.R

Modelo: MCT ~ Estatura
Coeficientes: a = -96.198781, b = 96.777161
x_min = 1.560, x_max = 1.930
x̄ = 1.763529
η(x̄) = 2.291769
η média = 2.306439
η máxima = 2.756299 (x = 1.560000)
η mínima = 2.062018 (x = 1.930000)

O modelo ajustado foi:

\[ \widehat{MCT} = -96.2 + 96.8\,\text{Estatura} \]

A elasticidade (não linear) é dada por:

\[ \eta(x) = \dfrac{96.8\,x}{-96.2 + 96.8\,x} \]

Como \(a < 0\) e \(b > 0\), a derivada \(\dfrac{d\eta}{dx} = \dfrac{a\,b}{(a + b\,x)^2}\) é negativa. Logo, a elasticidade é decrescente ao longo de \(x\).

No domínio observado (\(1.56 \le x \le 1.93\)):

\[ \begin{aligned} \eta_{\max} &= 2.76 \text{ em } x = 1.56 \\ \eta_{\min} &= 2.06 \text{ em } x = 1.93 \\ \eta(\bar{x}) &= 2.29 \end{aligned} \]

Interpretação:

A relação é elástica em todo o domínio, pois \(\eta(x) > 1\).
Isso significa que um aumento de 1% na Estatura está associado a um aumento médio de cerca de 2.3% na MCT. MCT cresce proporcionalmente mais do que a Estatura.

A elasticidade cai conforme a Estatura aumenta, indicando diminuição relativa de sensibilidade em indivíduos mais altos.

Note que a reta de regressão é linear e crescente, mas a elasticidade é não linear e decrescente.

A seguinte relação é inelástica com elasticidade constante \(\eta(\text{Corpo})=\dfrac{2}{3}\):

\[ \widehat{\text{cérebro}} = 0.18 \times {\text{corpo}}^{2/3} \]

Portanto, a função é potência crescente (não linear) com elasticidade linear e constante igual a \(2/3\).

• elasticidade_bodybrain.R

Modelo: Cérebro = 0.18 * Corpo^{0.667}
Elasticidade constante: eta = 0.667
x_min = 100.00, x_max = 80000.00

Regressão de Deming

• Silveira PSP, Vieira JE, Siqueira JO (2024) Is the Bland-Altman plot method useful without inferences for accuracy, precision, and agreement? Rev Saude Publica 58:01. DOI 10.11606/s1518-8787.2024058005430: PDF

Há situações em que queremos estabelecer a equivalência de duas técnicas de medição, especialmente quando queremos substituir um procedimento tradicional por outro que seja mais barato, mais rápido ou menos invasivo.

Por exemplo, suponha que queremos medir a massa corporal total de indivíduos usando uma balança profissional e pedindo-lhes que relatem seus pesos. A segunda técnica de medição é menos cara e possivelmente mais fácil de aplicar.

Outras técnicas candidatas também podem ser levadas em consideração. Avaliadores treinados podem estimar o peso dos indivíduos por sua aparência ou, ainda mais barato, esses árbitros podem ser recrutados sem serem treinados anteriormente. Outra opção é uma balança doméstica não profissional como instrumento de medida.

Em um experimento hipotético, voluntários foram mensurados com as quatro técnicas disponíveis (MCT4tecnicas.xlsx). A balança profissional é a técnica de referência.

Uma simples inspeção desta tabela (MCT4tecnicas.xlsx) mostra:

MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))
print(MCT)
print(psych::describe(MCT[2:6]))

   Subject Professional scale Self-report Trained referees Untrained referees
1        1               70.5          69               71                 66
2        2               70.5          69               70                 66
3        3               70.7          73               72                 72
4        4               71.7          72               71                 73
5        5               71.8          72               72                 73
6        6               72.2          71               70                 68
7        7               72.2          73               72                 71
8        8               72.5          73               73                 71
9        9               73.8          75               73                 72
10      10               74.1          73               72                 71
11      11               74.3          73               73                 72
12      12               74.3          76               74                 71
13      13               74.4          74               76                 72
14      14               74.4          74               75                 73
15      15               76.4          76               75                 73
16      16               77.1          78               77                 75
17      17               77.3          78               76                 75
18      18               79.0          80               81                 78
19      19               79.0          77               80                 80
20      20               79.3          77               81                 81
21      21               79.5          76               79                 77
22      22               80.5          81               84                 82
23      23               80.6          81               85                 83
24      24               80.8          80               81                 80
25      25               81.9          82               78                 79
26      26               82.0          82               84                 83
27      27               82.2          82               81                 80
28      28               83.1          84               85                 82
29      29               83.3          84               84                 81
30      30               85.6          85               87                 87
   Domestic scale
1              73
2              73
3              72
4              72
5              72
6              73
7              72
8              72
9              75
10             73
11             73
12             75
13             78
14             78
15             80
16             82
17             82
18             81
19             81
20             81
21             81
22             81
23             81
24             81
25             83
26             81
27             83
28             82
29             82
30             88
                   vars  n  mean   sd median trimmed  mad  min  max range  skew
Professional scale    1 30 76.83 4.46  76.75   76.72 5.86 70.5 85.6  15.1  0.17
Self-report           2 30 76.67 4.57  76.00   76.58 5.19 69.0 85.0  16.0  0.17
Trained referees      3 30 77.07 5.22  76.00   76.83 5.93 70.0 87.0  17.0  0.31
Untrained referees    4 30 75.57 5.46  74.00   75.58 5.19 66.0 87.0  21.0  0.14
Domestic scale        5 30 78.03 4.60  80.50   77.96 3.71 72.0 88.0  16.0 -0.03
                   kurtosis   se
Professional scale    -1.35 0.82
Self-report           -1.17 0.83
Trained referees      -1.34 0.95
Untrained referees    -1.03 1.00
Domestic scale        -1.31 0.84

Em relação à balança profissional, o autorrelato forneceu média e desvio padrão semelhantes. Árbitros treinados conseguiram aproximar os pesos individuais em média, mas mostraram maior variância. Árbitros não treinados subestimaram os pesos e também mostraram maior variância. A balança doméstica apresentou variância ligeiramente maior, mas parece enviesada para superestimar os pesos.

Como avaliar qual técnica de medição candidata pode substituir o uso de uma balança profissional?

Intuitivamente, avaliar a correlação ou executar regressão linear simples poderia resolver o problema. Iniciamos pelos gráficos de correlação e regressão para verificar se um modelo linear parece adequado (demo_MCT4correg.R):

# Dados
MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))

cor <- list()
reg <- list()
for (c.aux in 3:6)
{
  # scatterplot e regressao linear simples
  cor <- c(cor,correg(MCT$`Professional scale`, as.numeric(unlist(MCT[,c.aux])), 
                method="pearson",
                xlab="Professional scale", ylab=names(MCT)[c.aux],
                pch=21, col="#000000", bg=friendlycolor(24)))
  
  # scatterplot e regressao linear simples
  reg <- c(reg,correg(MCT$`Professional scale`, as.numeric(unlist(MCT[,c.aux])), 
                method="lm_robust",
                xlab="Professional scale", ylab=names(MCT)[c.aux],
                pch=21, col="#000000", bg=friendlycolor(24)))
}

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2        Self-report   72 4.5662   69      69     73      72   71 30    0
lambda =  1.62534363730857 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Self-report.

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

    Pearson's product-moment correlation

data:  Professional scale and Self-report
t = 18.535, df = 28, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9200416 0.9817465
sample estimates:
      cor 
0.9615822 

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2        Self-report   72 4.5662   69      69     73      72   71 30    0
lambda =  1.62534363730857 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Self-report.

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Self-report ~ Professional scale)"

Standard error type:  HC2 

Coefficients:
                   Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)          1.1083    3.80319  0.2914 7.729e-01   -6.682    8.899 28
Professional scale   0.9834    0.04902 20.0608 3.746e-18    0.883    1.084 28

Multiple R-squared:  0.9246 ,   Adjusted R-squared:  0.9219 
F-statistic: 402.4 on 1 and 28 DF,  p-value: < 2.2e-16

Equation:
  mean[Self-report]=1.1082923 + 0.98340617 . Professional scale

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2   Trained referees   71 5.2189   71      70     72      72   70 30    0
lambda =  0.466272486651781 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Trained referees.

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

    Pearson's product-moment correlation

data:  Professional scale and Trained referees
t = 15.839, df = 28, p-value = 1.66e-15
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8935247 0.9754312
sample estimates:
      cor 
0.9484719 

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2   Trained referees   71 5.2189   71      70     72      72   70 30    0
lambda =  0.466272486651781 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Trained referees.

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Trained referees ~ Professional scale)"

Standard error type:  HC2 

Coefficients:
                   Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)          -8.114    5.05682  -1.605 1.198e-01 -18.4721    2.245 28
Professional scale    1.109    0.06779  16.354 7.376e-16   0.9698    1.247 28

Multiple R-squared:  0.8996 ,   Adjusted R-squared:  0.896 
F-statistic: 267.5 on 1 and 28 DF,  p-value: 7.376e-16

Equation:
  mean[Trained referees]=-8.11370519 + 1.10863825 . Professional scale

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2 Untrained referees   73 5.4563   66      66     72      73   68 30    0
lambda =  0.363013075021059 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Untrained referees.

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

    Pearson's product-moment correlation

data:  Professional scale and Untrained referees
t = 14.063, df = 28, p-value = 3.239e-14
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8685208 0.9693503
sample estimates:
      cor 
0.9359348 

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2 Untrained referees   73 5.4563   66      66     72      73   68 30    0
lambda =  0.363013075021059 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Untrained referees.

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Untrained referees ~ Professional scale)"

Standard error type:  HC2 

Coefficients:
                   Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)         -12.312    7.04865  -1.747 9.166e-02  -26.750    2.127 28
Professional scale    1.144    0.09068  12.613 4.554e-13    0.958    1.329 28

Multiple R-squared:  0.876 ,    Adjusted R-squared:  0.8715 
F-statistic: 159.1 on 1 and 28 DF,  p-value: 4.554e-13

Equation:
  mean[Untrained referees]=-12.3115857 + 1.14375166 . Professional scale

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2     Domestic scale   72 4.5975   73      73     72      72   73 30    0
lambda =  2.42236873463687 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Domestic scale.

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

    Pearson's product-moment correlation

data:  Professional scale and Domestic scale
t = 13.118, df = 28, p-value = 1.772e-13
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8516701 0.9651813
sample estimates:
      cor 
0.9273883 

            Variable Mean     Sd Min. 1st Qu. Median 3rd Qu. Max.  n NA's
1 Professional scale 71.7 4.4649 70.5    70.5   70.7    71.8 72.2 30    0
2     Domestic scale   72 4.5975   73      73     72      72   73 30    0
lambda =  2.42236873463687 
Assuming:
    - Explanatory variable: Professional scale;
    - Dependent variable: Domestic scale.

------------------------
Simple linear regression
------------------------

Call:
"lm_robust(formula = Domestic scale ~ Professional scale)"

Standard error type:  HC2 

Coefficients:
                   Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)          4.6634    4.72843  0.9862 3.325e-01   -5.022   14.349 28
Professional scale   0.9549    0.06148 15.5332 2.716e-15    0.829    1.081 28

Multiple R-squared:   0.86 ,    Adjusted R-squared:  0.8551 
F-statistic: 241.3 on 1 and 28 DF,  p-value: 2.716e-15

Equation:
  mean[Domestic scale]=4.66338749 + 0.95492337 . Professional scale

Todos os quatro métodos estão fortemente correlacionados com as medidas da balança profissional de acordo com os tamanhos de efeito de Cohen. Os coeficientes de correlação de Pearson [com intervalos de confiança de 95%] são (A) autorrelato: 0.962 [0.92, 0.982]; (B) árbitros treinados: 0.948 [0.894, 0.975]; (C) árbitros não treinados: 0.936 [0.869, 0.969]; e (D) balança doméstica: 0.927 [0.852, 0.965].

Uma segunda abordagem possível é a regressão linear simples. Se escolhida, não pode ser aplicada de forma grosseira. Para avaliar a concordância, espera uma linha de regressão próxima à bissetriz, indicando que as técnicas de referência e candidata fornecem a mesma medida.

Para a reta de regressão, dada por \[y = \alpha + \beta x + \varepsilon\] onde, para a bissetriz, \(\alpha=0\) e \(\beta=1\), não adianta verificar a existência de modelo, como é habitual, testando \(H_0: \beta=0\), que é uma reta horizontal. O artifício para testar \(H_0: \beta=1\) é subtrair \(x\) dos dois lados da equação, obtendo:

\[ \begin{align} y-x &= \alpha + \beta x -x + \varepsilon \\ y-x &= \alpha + (\beta-1)x + \varepsilon \\ y-x &= \alpha + \beta'x + \varepsilon \end{align} \]

e testamos:

\[ H_0: \beta^{\prime}=0 \]

de forma que, agora, o teste da inclinação da regressão linear equivale a \(H_0: \beta=1\). Isto é testável com funções implementadas em demo_MCT4linreg.R:

# Dados
MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))

for (c.aux in 3:6)
{
  out_lr_rm <- mcr::mcreg(MCT$`Professional scale`, 
                          as.numeric(unlist(MCT[,c.aux])),
                          method.reg="LinReg", 
                          method.ci="bootstrap", nsamples = 1e4, rng.seed = 123, 
                          mref.name = "Professional scale", 
                          mtest.name = names(MCT)[c.aux], na.rm=TRUE)
  mcr::printSummary(out_lr_rm)
  mcr::plot(out_lr_rm)
  mcr::compareFit(out_lr_rm)
}

------------------------------------------

Reference method: Professional scale
Test method:     Self-report
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%


------------------------------------------

LINEAR REGRESSION FIT:

                EST SE        LCI      UCI
Intercept 1.1082923 NA -5.7870758 9.108592
Slope     0.9834062 NA  0.8798697 1.072006


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept    1.10829        1.18291  0.07462      3.79323
Slope        0.98341        0.98237 -0.00104      0.04909

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

------------------------------------------

Reference method: Professional scale
Test method:     Trained referees
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%


------------------------------------------

LINEAR REGRESSION FIT:

                EST SE         LCI      UCI
Intercept -8.113705 NA -18.1207564 1.888061
Slope      1.108638 NA   0.9730358 1.241214


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept   -8.11371       -8.29447 -0.18076      5.08957
Slope        1.10864        1.11099  0.00235      0.06813

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

------------------------------------------

Reference method: Professional scale
Test method:     Untrained referees
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%


------------------------------------------

LINEAR REGRESSION FIT:

                 EST SE         LCI       UCI
Intercept -12.311586 NA -25.7163746 0.8810338
Slope       1.143752 NA   0.9737208 1.3162437


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept  -12.31159      -12.34667 -0.03509      6.84761
Slope        1.14375        1.14419  0.00044      0.08820

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

------------------------------------------

Reference method: Professional scale
Test method:     Domestic scale
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%


------------------------------------------

LINEAR REGRESSION FIT:

                EST SE        LCI       UCI
Intercept 4.6633875 NA -4.5666159 13.977038
Slope     0.9549234 NA  0.8355433  1.075789


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept    4.66339        4.58732 -0.07607      4.73994
Slope        0.95492        0.95601  0.00109      0.06152

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

Observe que todos os testes da regressão não rejeitam a hipótese nula, testando \(H_{0,1}: \alpha=0\) e \(H_{0,2}: \beta=1\). Portanto, poderíamos considerar os quatro métodos como vicários para a balança profissional. No entanto, em dois casos a bissetriz não está totalmente contida na banda de confiança 95%.

Os autores citados afirmam que a correlação de Pearson (Bland & Altman, 1986) e a regressão linear simples (Bland & Altman, 2003) não são adequadas para resolver este tipo de problema. Por isso, propuseram seu famoso método (implementado em demo_MCT4baclassic.R):

# Dados
MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))

for (c.aux in 3:6)
{
  out_lr_rm <- mcr::mcreg(MCT$`Professional scale`, as.numeric(unlist(MCT[,c.aux])),
                    method.reg="LinReg", 
                    method.ci="bootstrap", nsamples = 1e4, rng.seed = 123, 
                    mref.name = "Professional scale", 
                    mtest.name = names(MCT)[c.aux], na.rm=TRUE)
  mcr::plotDifference(out_lr_rm)
}

O método gráfico de Bland-Altman não tem teste estatístico e sua interpretação é subjetiva. Em suma, os gráficos de Bland-Altman avaliam os limites de concordância de 95% (LoA) dados pela diferença média das medições individuais fornecidas por duas técnicas mais ou menos 1.96 vezes seu desvio padrão plotado como linhas horizontais. Quando os LoA são menores do que uma dada tolerância (diferenças pequenas o suficiente para serem consideradas clinicamente sem importância), as duas técnicas são consideradas como tendo equivalência aceitável. Mais recentemente, intervalos de confiança foram adicionados em torno da LoA superior e inferior. Embora forneça um pouco mais de espaço para a tolerância, as LoAs podem ser consideradas um teste estatístico para os limites da banda, mas não fornecem qualquer decisão para a equivalência das técnicas.

Nos gráficos de Bland-Altman a distância entre a linha sólida em vermelho e o valor zero em preto é o viés (bias, em inglês). Precisamos definir se este viés é nulo. Este teste pode ser feito com funções implementadas em demo_MCT4bias.R:

# Dados
MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))

for (c.aux in 3:6)
{
  out_lr_rm <- mcr::mcreg(MCT$`Professional scale`, as.numeric(unlist(MCT[,c.aux])),
                    method.reg="LinReg", 
                    method.ci="bootstrap", nsamples = 1e4, rng.seed = 123, 
                    mref.name = "Professional scale", 
                    mtest.name = names(MCT)[c.aux], na.rm=TRUE)
  mcr::plotBias(out_lr_rm, 
                xlab="Professional scale",
                ylab=paste0(names(MCT)[c.aux]," - Professional scale"))
  
}

Note que nestes gráficos o eixo das ordenadas é o mesmo da regressão linear simples modificada, \(y-x\). De acordo com esta análise, os árbitros não treinados (subestimando) e a balança doméstica (superestimando) não obtiveram os mesmos valores que a balança profissional em média porque a linha horizontal pontilhada na altura da diferença nula (i.e., ausência de viés) não está totalmente contida na banda de confiança de 95%.

Até aqui as análises não consideraram que toda técnica de mensuração tem erro e, portanto, a balança profissional que adotamos como referência não é exceção. A regressão linear simples não é suficiente porque assume que a VE (\(x\)) não tem erro de medida. A regressão que considera erros de medida em \(x\) e \(y\) e permite testar a equivalência de técnicas é a Regressão de Deming.

De acordo com Shukla (1973), ao calcular a regressão linear como \(y-x = \alpha + \beta(y+x) + \theta\), a inclinação será diferente de zero quando a variância de duas técnicas de medição for diferente, sendo a variância um sucedâneo das precisões da técnica de medição.

Aqui observamos que os eixos propostos por Shukla são os mesmos usados pelo conceito original de Bland e Altman em 1986: a diferença entre as medidas representadas no eixo \(y\) e a soma (ou a média, que não altera a regressão) na eixo \(x\). Este arranjo, no entanto, é projetado para comparar a variância entre as técnicas de medição e não uma equivalência completa como os gráficos de Bland-Altman sugeririam.

A precisão de uma técnica de medida é a variância de seus erros de medida. Supõe que os erros de medida são independentes intra e entre técnicas (aqui, cada técnica é aplica somente uma vez em cada sujeito). Então, para avaliar igualdade das precisões das duas técnicas de medida estimaremos o intervalo de confiança 95% de \(\lambda\). Assuma lambda como a razão das precisões das duas técnicas e, portanto, serão iguais quando \(\lambda=1\) estiver contido ao IC95. O cálculo de lambda está implementado em demo_MCT4shukla.R:

# Estimacao do IC95% do lambda sem medidas repetidas (medida única)
# A ordem de x e y importa
# Shukla (1973), p. 376.

# Dados
MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))

for (c.aux in 3:6)
{
  cat("\n\nProfessional scale x ",names(MCT)[c.aux],sep="")

  dt_xy <- data.frame(MCT$`Professional scale`,as.numeric(unlist(MCT[,c.aux])))
  names(dt_xy) <- c("x","y")
  dt_xy <- na.omit(dt_xy)
  n <- nrow(dt_xy)
  a <- var(dt_xy$x) - cov(dt_xy$x,dt_xy$y)
  b <- var(dt_xy$y) - cov(dt_xy$x,dt_xy$y)
  gl <- n - 2
  alpha <- 0.05
  t <- qt(1 - alpha/2, gl)
  P <- (t^2)*(var(dt_xy$x)*var(dt_xy$y) - 
                cov(dt_xy$x,dt_xy$y)^2)/(n - 2)
  LI <- abs((b - sqrt(P))/(a + sqrt(P)))
  LS <- abs((b + sqrt(P))/(a - sqrt(P)))
  cat("\nIC95%(lambda) = [", round(LI,4), ", ", round(LS,4),"]",sep="")
  lambda <- (LI + LS)/2
  cat("\nEstimativa pontual de lambda = ",round(lambda,4),sep="") 
}

Professional scale x Self-report
IC95%(lambda) = [0.3686, 1.8589]
Estimativa pontual de lambda = 1.1138

Professional scale x Trained referees
IC95%(lambda) = [3.2888, 1.5911]
Estimativa pontual de lambda = 2.44

Professional scale x Untrained referees
IC95%(lambda) = [8.0101, 1.6634]
Estimativa pontual de lambda = 4.8368

Professional scale x Domestic scale
IC95%(lambda) = [0.2254, 2.4457]
Estimativa pontual de lambda = 1.3356

De acordo com esta análise, os árbitros treinados e não treinados têm precisão menor do que a balança profissional (intervalo de confiança de \(\lambda\) está acima de 1), enquanto o autorrelato e a balança doméstica tem precisão estatisticamente não diferente da balança profissional (\(\lambda=1\) está contido no intervalo de confiança de 95%).

Esta relação entre as precisões é calculada porque a regressão de Deming a leva em conta. As regressões de Deming estão implementadas em demo_MCT4deming.R:

# Dados
MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))

for (c.aux in 3:6)
{
  cat("\n\nProfessional scale x ",names(MCT)[c.aux],sep="")

  dt_xy <- data.frame(MCT$`Professional scale`,as.numeric(unlist(MCT[,c.aux])))
  names(dt_xy) <- c("x","y")
  dt_xy <- na.omit(dt_xy)
  n <- nrow(dt_xy)
  a <- var(dt_xy$x) - cov(dt_xy$x,dt_xy$y)
  b <- var(dt_xy$y) - cov(dt_xy$x,dt_xy$y)
  gl <- n - 2
  alpha <- 0.05
  t <- qt(1 - alpha/2, gl)
  P <- (t^2)*(var(dt_xy$x)*var(dt_xy$y) - 
                cov(dt_xy$x,dt_xy$y)^2)/(n - 2)
  LI <- abs((b - sqrt(P))/(a + sqrt(P)))
  LS <- abs((b + sqrt(P))/(a - sqrt(P)))
  lambda <- (LI + LS)/2
  # regressao de Deming
  out_boot_rm <- mcr::mcreg(dt_xy$x,dt_xy$y,
                            error.ratio=lambda,
                            method.reg="Deming", 
                            # method.ci="analytical",  
                            method.ci="bootstrap", nsamples = 1e5, rng.seed = 123,                             
                            mref.name = "Professional scale", 
                            mtest.name = names(MCT)[c.aux], 
                            na.rm=TRUE)
  mcr::printSummary(out_boot_rm)
  mcr::plot(out_boot_rm)
}

Professional scale x Self-report

------------------------------------------

Reference method: Professional scale
Test method:     Self-report
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%
Error ratio: 1.113771

------------------------------------------

DEMING REGRESSION FIT:

                EST SE        LCI      UCI
Intercept -2.149930 NA -8.4225592 4.907569
Slope      1.025813 NA  0.9351107 1.106124


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept   -2.14993       -2.08661  0.06332      3.41496
Slope        1.02581        1.02494 -0.00087      0.04372

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

Professional scale x Trained referees

------------------------------------------

Reference method: Professional scale
Test method:     Trained referees
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%
Error ratio: 2.439952

------------------------------------------

DEMING REGRESSION FIT:

                EST SE        LCI       UCI
Intercept -15.38611 NA -26.740344 -5.548465
Slope       1.20329 NA   1.071791  1.352757


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept  -15.38611      -15.54506 -0.15894      5.38265
Slope        1.20329        1.20525  0.00196      0.07132

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

Professional scale x Untrained referees

------------------------------------------

Reference method: Professional scale
Test method:     Untrained referees
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%
Error ratio: 4.836759

------------------------------------------

DEMING REGRESSION FIT:

                EST SE       LCI        UCI
Intercept -23.22027 NA -35.07557 -10.240593
Slope       1.28573 NA   1.11768   1.440764


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept  -23.22027      -23.11475  0.10552      6.28738
Slope        1.28573        1.28441 -0.00132      0.08165

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

Professional scale x Domestic scale

------------------------------------------

Reference method: Professional scale
Test method:     Domestic scale
Number of data points: 30

------------------------------------------

The confidence intervals are calculated with
 bootstrap  ( quantile ) method.
Confidence level: 95%
Error ratio: 1.335583

------------------------------------------

DEMING REGRESSION FIT:

                EST SE         LCI      UCI
Intercept -2.153940 NA -13.1508293 7.519870
Slope      1.043652 NA   0.9169883 1.189889


------------------------------------------

BOOTSTRAP SUMMARY

          global.est bootstrap.mean     bias bootstrap.se
Intercept   -2.15394       -2.29466 -0.14072      5.25148
Slope        1.04365        1.04570  0.00204      0.06932

Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123

Observando as bissetrizes e as bandas de confiança, concluímos que o autorrelato é a única técnica equivalente ao uso da balança profissional.

Tomamos as decisões por inspeção visual. Ainda podemos formalizar com um teste estatístico.

Note que a saída textual das regressões de Deming exibiram intervalos de confiança 95% para o intercepto (não rejeitando a hipótese nula quando o intervalo contém \(0\)) e coeficiente angular (não rejeitando a hipótese nula quando o intervalo contém o valor \(1\)). Não é suficiente:

• No caso do autorrelato, não se rejeita \(H_0\) nos dois casos e a bissetriz está graficamente dentro da banda, coerentemente.
• Nos casos dos árbitros treinados ou não treinados, ambas \(H_0\) são rejeitadas e a bissetriz não está inteiramente dentro das bandas de confiança 95%, dando a falsa impressão de que se pode confiar nesta saída textual.
• No entanto, no caso da balança doméstica, não se rejeita \(H_0\) pelos intervalos de confiança do intercepto e do coeficiente angular, mas a bissetriz está quase inteiramente fora da banda de confiança 95%.

Para testar conjuntamente o intercepto e o coeficiente angular, é necessário verificar se a hipótese nula \(H_0: (\alpha,\beta) = (0,1)\) encontra dentro da região elíptica de confiança 95% (Francq & Govaerts B, 2014 e 2016). Está implementado em demo_MCT4elipse.R:

source("eiras.deming.regression.R")
# Dados
MCT <- as.data.frame(readxl::read_excel("MCT4tecnicas.xlsx"))

for (c.aux in 3:6)
{
  dt_xy <- data.frame(MCT$`Professional scale`,as.numeric(unlist(MCT[,c.aux])))
  xlab <- "Professional scale"
  ylab <- names(MCT)[c.aux]
  names(dt_xy) <- c("x","y")
  dt_xy <- na.omit(dt_xy)
  deming.regression(dt_xy,xlab=names(MCT)[2],ylab=names(MCT)[c.aux])
}

Computing regressions: Professional scale x Self-report

---------------------------------------------

Simple linear regression:
___ modified to test y-x ~ x ___

Call:
estimatr::lm_robust(formula = diff ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)   1.1083    3.80319  0.2914   0.7729   -6.682    8.899 28
x             0.9834    0.04902 -0.3385   0.7375    0.883    1.084 28

Multiple R-squared:  0.003481 , Adjusted R-squared:  -0.03211 
F-statistic: 0.1146 on 1 and 28 DF,  p-value: 0.7375
                 2.5 %   97.5 %
(Intercept) -6.6821900 8.898775
x            0.8829908 1.083822

---------------------------------------------

Deming regression (analytical):

lambda = 1.113771
intercept = -1.811883
slope = 1.021413

Call:
estimatr::lm_robust(formula = y ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)   -1.812    3.50371 -0.5171   0.6091  -8.9889    5.365 28
x              1.021    0.04473  0.4787   0.6358   0.9298    1.113 28

Multiple R-squared:  0.005233 , Adjusted R-squared:  -0.03029 
F-statistic: 0.2292 on 1 and 28 DF,  p-value: 0.6358
                 2.5 %   97.5 %
(Intercept) -8.9889065 5.365141
true.x       0.9297908 1.113035

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): not rejected


---------------------------------------------

Deming regression (bootstrapping):
    n=30, B=10000, alpha=0.05

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): not rejected
H0: regression line is bisector
    (not rejected when (0,1) is inside 95 elliptical confidence region)
    Decision (alpha=0.05): not rejected

Computing regressions: Professional scale x Trained referees

---------------------------------------------

Simple linear regression:
___ modified to test y-x ~ x ___

Call:
estimatr::lm_robust(formula = diff ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)   -8.114    5.05682  -1.605   0.1198 -18.4721    2.245 28
x              1.109    0.06779   1.603   0.1202   0.9698    1.247 28

Multiple R-squared:  0.07922 ,  Adjusted R-squared:  0.04634 
F-statistic: 2.568 on 1 and 28 DF,  p-value: 0.1202
                  2.5 %   97.5 %
(Intercept) -18.4721349 2.244725
x             0.9697786 1.247498

---------------------------------------------

Deming regression (analytical):

lambda = 2.439952
intercept = -11.37846
slope = 1.15113

Call:
estimatr::lm_robust(formula = y ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)  -11.378    6.42247  -1.772  0.08733 -24.5343    1.777 28
x              1.151    0.08584   1.761  0.08924   0.9753    1.327 28

Multiple R-squared:  0.08806 ,  Adjusted R-squared:  0.05549 
F-statistic:   3.1 on 1 and 28 DF,  p-value: 0.08924
                  2.5 %   97.5 %
(Intercept) -24.5343008 1.777372
true.x        0.9752894 1.326970

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): not rejected


---------------------------------------------

Deming regression (bootstrapping):
    n=30, B=10000, alpha=0.05

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): not rejected
H0: regression line is bisector
    (not rejected when (0,1) is inside 95 elliptical confidence region)
    Decision (alpha=0.05): not rejected

Computing regressions: Professional scale x Untrained referees

---------------------------------------------

Simple linear regression:
___ modified to test y-x ~ x ___

Call:
estimatr::lm_robust(formula = diff ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)  -12.312    7.04865  -1.747  0.09166  -26.750    2.127 28
x              1.144    0.09068   1.585  0.12413    0.958    1.329 28

Multiple R-squared:  0.1004 ,   Adjusted R-squared:  0.06824 
F-statistic: 2.513 on 1 and 28 DF,  p-value: 0.1241
                 2.5 %   97.5 %
(Intercept) -26.750087 2.126916
x             0.958007 1.329496

---------------------------------------------

Deming regression (analytical):

lambda = 4.836759
intercept = -15.02433
slope = 1.179059

Call:
estimatr::lm_robust(formula = y ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)  -15.024    10.6695  -1.408   0.1701 -36.8799    6.831 28
x              1.179     0.1374   1.303   0.2032   0.8976    1.461 28

Multiple R-squared:  0.06609 ,  Adjusted R-squared:  0.03274 
F-statistic: 1.698 on 1 and 28 DF,  p-value: 0.2032
                  2.5 %   97.5 %
(Intercept) -36.8798737 6.831215
true.x        0.8975859 1.460531

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): rejected


---------------------------------------------

Deming regression (bootstrapping):
    n=30, B=10000, alpha=0.05

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): rejected
H0: regression line is bisector
    (not rejected when (0,1) is inside 95 elliptical confidence region)
    Decision (alpha=0.05): rejected

Computing regressions: Professional scale x Domestic scale

---------------------------------------------

Simple linear regression:
___ modified to test y-x ~ x ___

Call:
estimatr::lm_robust(formula = diff ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)   4.6634    4.72843  0.9862   0.3325   -5.022   14.349 28
x             0.9549    0.06148 -0.7332   0.4695    0.829    1.081 28

Multiple R-squared:  0.01351 ,  Adjusted R-squared:  -0.02172 
F-statistic: 0.5376 on 1 and 28 DF,  p-value: 0.4695
                 2.5 %    97.5 %
(Intercept) -5.0223688 14.349144
x            0.8289952  1.080852

---------------------------------------------

Deming regression (analytical):

lambda = 1.335583
intercept = -0.3729917
slope = 1.020473

Call:
estimatr::lm_robust(formula = y ~ x)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept)   -0.373    5.37490 -0.0694   0.9452 -11.3830   10.637 28
x              1.020    0.07062  0.2899   0.7740   0.8758    1.165 28

Multiple R-squared:  0.002108 , Adjusted R-squared:  -0.03353 
F-statistic: 0.08403 on 1 and 28 DF,  p-value: 0.774
                  2.5 %    97.5 %
(Intercept) -11.3829681 10.636985
true.x        0.8758042  1.165141

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): rejected


---------------------------------------------

Deming regression (bootstrapping):
    n=30, B=10000, alpha=0.05

H0: regression line is bisector
    (not rejected when bisector is inside 95 confidence band)
    Decision (alpha=0.05): rejected
H0: regression line is bisector
    (not rejected when (0,1) is inside 95 elliptical confidence region)
    Decision (alpha=0.05): rejected

O gráfico da regressão de Deming mostra a bissetriz em linha pontilhada e a reta de regressão em linha sólida preta. As bandas de confiança 95% foram estimadas analiticamente e por bootstrapping. As elipses de confiança foram estimadas a partir da regressão de Deming analiticamente e por bootstrapping; para comparação, a elipse obtida pela regressão linear simples tradicional também é exibida.

Por esta última análise, as duas técnicas mais elegíveis para reproduzir as medições da balança profissional são o autorrelato e o emprego de árbitros treinados.

Os árbitros não treinados estão no limite, mas \(H_0\) é rejeitada.

A balança doméstica mostra \(H_0\) claramente fora da elipse de confiança 95%, possivelmente devido a serem estes os dois métodos que têm viés das medidas.

Neste exemplo é necessário notar que em todos os casos a regressão linear simples também conduz às mesmas conclusões. Isto nem sempre ocorre.

FIM

Referências

• Anscombe, FJ (1973) Graphs in Statistical Analysis. The American Statistician 27(1): 17-21
• Bland JM, Altman DG (2003) Applying the right statistics: analyses of measurement studies. Ultrasound Obst Gyn. 22(1):85-93.
• Bland JM, Altman DG (1986) Statistical methods for assessing agreement between two methods of clinical measurement. Lancet 1(8476):307‑10.
• von Bonin, G (1937) Brain-Weight and BodyWeight of Mammals. The Journal of General Psychology 16(2): 379-89. DOI: 10.1080/00221309.1937.9917959
• Cohen, J (1992) A power primer. Quantitative methods in Psychology 112(1): 155-9.
• 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. 1st ed. Cambridge University Press.
• Francq BG, Govaerts B (2014) Hyperbolic confidence bands of errors-in-variables regression lines applied to method comparison studies. Journal de la Société Française de Statistique 155(1). http://www.numdam.org/item/JSFS_2014__155_1_23_0.pdf
• Francq BG, Govaerts B (2016) How to regress and predict in a Bland-Altman plot? Review and contribution based on tolerance intervals and correlated-errors-in-variables models. Statist. Med., 35: 2328-58. https://doi.org/10.1002/sim.6872
• 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.
• Howell, D (2013) Statistical Methods in Psychology. 8th ed. USA: CENGAGE.
• 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.
• Nihad Mohammed, University of Anbar. Traduzido e modificado de https://www.researchgate.net/post/What_is_the_key_differences_between_correlation_and_regression, 08Aug2018
• 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.
• Revelle W (2014) Personality Project. http://personality-project.org/revelle.html.
• Shinohara, Elvira Maria Guerra. 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]. São Paulo: Universidade de São Paulo, Faculdade de Ciências Farmacêuticas; 1989 [citado 2019-04-23]. 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.
• Silveira PSP, Vieira JE & Siqueira JO (2024) Is the Bland-Altman plot method useful without inferences for accuracy, precision, and agreement? Rev Saude Publica 58:01. DOI 10.11606/s1518-8787.2024058005430: PDF
• Shukla GK (1973) Some exact tests of hypotheses about Grubbs’s estimators. Biometrics 29(2):373-7.
• Venditti, C et al. (2024) Co-evolutionary dynamics of mammalian brain and body size. Nat Ecol Evol 8, 1534–42. https://doi.org/10.1038/s41559-024-02451-3
• Viraj Bhagat, 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, 13Apr2018.
• von Bonin, G (1937). Brain-weight and body-weight of mammals. The Journal of General Psychology 16(2): 379-89. https://doi.org/10.1080/00221309.1937.9917959.
• Wooldridge, JM (2016) Introductory econometrics: A modern approach. 6ª ed. Boston: Cengage Learning.
• Zhang, Z & Yuan, K-H (2018). Practical statistical power analysis using WebPower and R (Eds). IN: ISDSA.