Análise de Regressão Linear Simples & Regressão de Deming
em R
13 agosto 2025 21:16h
Bastão de Asclépio & Distribuição Normal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))options(warn=-1)
suppressMessages(library(bootES, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts = FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts = FALSE))
suppressMessages(library(estimatr, warn.conflicts = FALSE))
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(sp, warn.conflicts = FALSE))
suppressMessages(library(WebPower, warn.conflicts = FALSE))
source("BiNormal.R")
source("eiras.bartitle.R")
source("eiras.col2rgbstring.R")
source("eiras.ConfidenceBand.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")Pacote eirasagree
Instalar antes os pacotes eirasdata e eiras
- Pacote
eirasagree
Funções
AllStructuralTestsStructural Tests, execute all, in sequenceConfidenceBandAnalytic Confidence BandDemingPrepareComputing true values for Deming regressionDemingRegressionDeming Regression and Statistical TestDemingRsquaredR-squared for Deming RegressionDescriptiveStatShow data frame content and summaryEllipticalRegionAnalytic Elliptical Confidence RegionHelloHello from eirasLambdaEstimateEstimation of lambda, \(\lambda=\dfrac{V[\delta]}{V[\varepsilon]}\)plotBlandAltmanBland Altman plot methodplotRawDataPlot raw dataRobustConfidenceBandRobust confidence bandRobustEllipticalRegionBootstrap Elliptical Confidence RegionSymbolsAgreementSymbols for AgreementtestAccuracyAccuracy testtestDemingDeming Statistical TesttestPrecisionPrecision testtestReliabilityReliability test
Material
- HTML de R Markdown em
RPubs
Incluir
HH::ci.plot: Plot confidence and prediction intervals for simple linear regression
tmp <- data.frame(x=rnorm(20), y=rnorm(20)) tmp.lm <- lm(y ~ x, data=tmp) HH::ci.plot(tmp.lm)
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 ou mais variáveis explicativas/previsoras (VE) intervalares.
- Prever o escore de uma unidade observacional na VD a partir dos valores conhecidos de uma VE.
- Representar graficamente a relação entre as variáveis de um estudo e a reta de regressão a partir da equação de regressão obtida.
- Testar a significância dos coeficientes obtidos em um estudo de regressão linear.
- Determinar o tamanho de efeito da VE.
Regressão linear simples
A análise de regressão linear simples (RLS) é uma extensão da análise de correlação.
Ao contrário da correlação, a RLS é direcional e dimensional. Queremos estabelecer uma equação (uma reta, no caso da regressão linear simples) 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.
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):
source("eiras.friendlycolor.R")
# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
# scatterplot
plot(estatura, massa,
xlim = c(163,182), ylim = c(55,90),
xlab="Estatura (cm)", ylab="Massa (kg)",
pch=21, col="black", bg=friendlycolor(24)
)
lines(lowess(estatura,massa), lty=2)
# modelo linear, plota a reta
modelo <- lm(massa ~ estatura)
print(modelo)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
massa.medio <- intercepto + inclinacao*estatura
lines (estatura, massa.medio, col=friendlycolor(8),
lwd=3, lty=1)
sumario <- summary(modelo)
print(sumario)
Call:
lm(formula = massa ~ estatura)
Coefficients:
(Intercept) estatura
-144.525 1.256
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.02445Neste gráfico:
- A variável em \(x\) é a VE (estatura); variável explicativa, independente.
- A variável em \(y\) é a VD (massa corporal); variável de desfecho, dependente.
A reta foi estimada a partir da função lm, a qual é
parte do pacote 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 parâmetros intercepto
(corresponde a \(\hat{\beta_0}\)) e
inclinacao (corresponde a \(\hat{\beta_1}\)) são valores ótimos e
únicos. Consequentemente, esta linha de código computa
massa.medio (corresponde a \(\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{mc}\) \(=\) -144.53 \(+\) 1.26 \(\text{estatura}\)
\(172 \le \text{estatura} \le 180\)
Note o circunflexo sobre o valor \(\hat{y}\) (ou \(\widehat{mc}\)), indicando que este é o valor médio estimado da VD dado o valor da VE (\(x\) ou \(\text{estatura}\)).
|
Sinonímia Dependendo do contexto da análise encontramos:
|
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, é muito pequena) (demo_corregEstMCT.R):
source("eiras.friendlycolor.R")
source("eiras.correg.R")
# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
# scatterplot e regressao linear simples
lst <- correg(estatura, massa, method="lm",
xlab="Estatura (cm)", ylab="Massa (kg)",
pch=21, col="#000000", bg=friendlycolor(24))
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 é o conjunto das médias de VD condicionadas a todos os valores contínuos da VE no seu intervalo observado.
\[ \widehat{\text{VD}} = \hat{\alpha}+\hat{\beta} \times \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.
Testes estatísticos
- As suposições principais para
lmsã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 tracejada horizontal na altura
da média das massas corporais observadas para representar \(H_0\); como há muitos elementos no gráfico,
também incluímos uma legenda (demo_regressao.R):
# demo_regressao.R
source("eiras.friendlycolor.R")
source("eiras.ConfidenceBand.R")
# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
# scatterplot
plot(estatura, massa,
xlim = c(163,182), ylim = c(55,95),
xlab="Estatura (cm)", ylab="Massa (kg)",
pch=21, col="black", bg=friendlycolor(24)
)
# modelo linear, traca a reta
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)
# 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=2)
# centroide, plota o ponto
points(centro_estatura, centro_massa, pch=21, col="black", bg="black")
# 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))
}
# legenda
legend("topleft",
c("valor observado","média de massa",
paste("reta: massa_media = ",round(intercepto,3)," + ",
round(inclinacao,3)," x estatura",sep=""),
"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)
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.02445Como vimos acima, a hipótese nula é rejeitada para \(\alpha=0.05\), pois \(p=0.0328\) 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.
Os testes estatísticos para correlação e RLS dão 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 da nulidade 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:
# Correlacao_r_de_Pearson.R
# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
correlacao <- cor.test(estatura, massa)
print(correlacao)
Pearson's product-moment correlation
data: estatura and massa
t = 2.8568, df = 7, p-value = 0.02445
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.1356665 0.9398558
sample estimates:
cor
0.733684 Isto é dizer que, para modelos lineares simples, é 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 (\(r=0\)) produzem retas de regressão horizontais (\(\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 for 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.friendlycolor.R")
source("eiras.ConfidenceBand.R")
# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
# scatterplot
plot(estatura, massa,
xlim = c(163,182), ylim = c(55,90),
xlab="Estatura (cm)", ylab="Massa (kg)",
pch=21, col="black", bg=friendlycolor(24)
)
# modelo linear, traca a reta
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)
# 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=2)
# centroide, plota o ponto
points(centro_estatura, centro_massa, pch=21, col="black", bg="black")
# banda de confianca
alfa <- 0.05
B <- 1e3
lst <- ConfidenceBand(x=estatura, y=massa, alpha=alfa, B=B)
band <- lst[[2]]
lines(band$X, band$LB, col=friendlycolor(8), lwd=2, lty=2)
lines(band$X, band$UB, col=friendlycolor(8), lwd=2, lty=2)
Dentro desta banda, com confiança de 95%, é onde está a reta de regressão populacional. Coerente com os testes estatísticos anteriores, uma outra forma de é verificar que não há uma reta horizontal que possa ser traçada para ficar inteiramente contida nesta banda. Portanto a reta populacional de onde esta amostra foi retirada, seja ela qual for, deve ter inclinação diferente de zero.
Coeficiente de determinação
Quando imprimimos o resultado de lm observamos uma nova
medida (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
estimar 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.732 =
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 foi 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\) global pois a RLS é um caso particular do modelo linear geral.
Observe, também, o valor p associado ao R-squared por uma estatística F. Seu valor é igual aos valores p computados para a inclinação da reta e para a correlação.
Ellis, 2010
|
Não confunda o coeficiente de determinação, \(R^2\), com o \(R^2\) ajustado ( |
Ajuste de reta
Note, no gráfico acima, que adicionamos um centróide em \((\bar{x},\bar{y})\) localizado no cruzamento da reta ajustada com a reta horizontal pontilhada (um ponto preto). Seu valor é dado pela média dos pares de valores de estatura e de massa corpórea observados.
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 for mínima. Esta reta ajustada pelo método de mínimos quadrados ordinários é a que foi representada no gráfico (toda reta ajustada por este método passa obrigatoriamente pelo centróide).
Execute demo_minimos_quadrados.R para uma
demonstração:
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 sólida é inclinada, diferencia da horizontal \(H_0\) e, portanto, com a equação encontrada, dada a VE, poderemos calcular a VD. No entanto, observando os valores da amostra (círculos amarelos), 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}\)). É como considerar que o tanto de inclinação da reta em relação à horizontal corresponde a quanto da associação dos pontos observados foram capturados pelo modelo linear. 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(172, 173, 174, 175, 176, 177, 178, 179, 180)
media <- mean(estatura)
estatura <- estatura - media
cat(estatura)-4 -3 -2 -1 0 1 2 3 4A regressão com a estatura centrada está implementada em demo_regressao_centrada.R, resultando
em:
# demo_regressao.R
source("eiras.friendlycolor.R")
source("eiras.ConfidenceBand.R")
# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
# centrando x
estatura <- estatura-mean(estatura)
cat("estatura(centrada): ",estatura,"\n")
cat("massa(original): ",massa,"\n")
# scatterplot
plot(estatura, massa,
xlim = c(min(estatura),max(estatura)),
ylim = c(55,95),
xlab="Estatura centrada (cm)", ylab="Massa (kg)",
pch=21, col="black", bg=friendlycolor(24),
axes=FALSE
)
axis(1)
axis(2, pos=0)
# modelo linear, traca a reta
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)
# 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=2)
# centroide, plota o ponto
points(centro_estatura, centro_massa, pch=21, col="black", bg="black")
# banda de confianca
alfa <- 0.05
B <- 1e3
lst <- ConfidenceBand(x=estatura, y=massa, alpha=alfa, B=B)
band <- lst[[2]]
lines(band$X, band$LB, col=friendlycolor(8), lwd=2, lty=2)
lines(band$X, band$UB, col=friendlycolor(8), lwd=2, lty=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))
}
# legenda
legend("topleft",
c("valor observado","massa média",
paste("reta: massa_media = ",round(intercepto,3)," + ",
round(inclinacao,3)," x estatura centrada",sep=""),
"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)estatura(centrada): -8.111111 -4.111111 -2.111111 -1.111111 -0.1111111 0.8888889 2.888889 4.888889 6.888889
massa(original): 61 73 68 74 65 82 69 81 83 
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.02445Observe 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 as duas
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")
source("eiras.ConfidenceBand.R")
# valores
estatura <- c(165, 169, 171, 172, 173, 174, 176, 178, 180)
massa <- c(61, 73, 68, 74, 65, 82, 69, 81, 83)
# padronizando x e y
estatura <- scale(estatura)
massa <- scale(massa)
cat("z_estatura: ",round(estatura,2),"\n")
cat("z_massa: ",round(massa,2),"\n")
# scatterplot
plot(estatura, massa,
xlim = c(-2,2),
ylim = c(-2,2),
xlab="Estatura (z)", ylab="Massa (z)",
pch=21, col="black", bg=friendlycolor(24),
axes=0
)
axis(1, pos=0)
axis(2, pos=0)
# modelo linear, traca a reta
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)
# 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=2)
# centroide, plota o ponto
points(centro_estatura, centro_massa, pch=21, col="black", bg="black")
# banda de confianca
alfa <- 0.05
B <- 1e3
lst <- ConfidenceBand(x=estatura, y=massa, alpha=alfa, B=B)
band <- lst[[2]]
lines(band$X, band$LB, col=friendlycolor(8), lwd=2, lty=2)
lines(band$X, band$UB, col=friendlycolor(8), lwd=2, lty=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))
}
# legenda
legend("topleft",
c("valor observado","massa centrada média",
paste("reta: massa_media = ",round(intercepto,3)," + ",
round(inclinacao,3)," x estatura centrada",sep=""),
"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)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.02445Com 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) e
\(~~~~~~~~\hat{\beta}_0 = 0\) (intercepto nulo).
O valor de \(R^2\) continua inalterado e assim concluímos que, em uma regressão linear qualquer, 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.
Planejamento amostral
O planejamento amostral pode ser feito com a função
WebPower::wp.regression. Por exemplo, 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 = ",f2,"\n",sep="")
wp <- WebPower::wp.regression(power=0.8,
alpha=0.05,
p1=1,
f2=f2)
print(wp)f^2 = 0.1494253Power for multiple regression
n p1 p2 f2 alpha power
54.51598 1 0 0.1494253 0.05 0.8
URL: http://psychstat.org/regressionPara poder de 90%, obtemos
R2 <- 0.13
f2 <- R2/(1-R2) # = 0.15
cat("f^2 = ",f2,"\n",sep="")
wp <- WebPower::wp.regression(power=0.9,
alpha=0.05,
p1=1,
f2=f2)
print(wp)f^2 = 0.1494253Power for multiple regression
n p1 p2 f2 alpha power
72.29472 1 0 0.1494253 0.05 0.9
URL: http://psychstat.org/regressionEsta função serve para planejamento de diversos modelos de regressão.
O parâmetro p1 é o número de VEs, que no caso da RLS é de
apenas 1. O parâmetro f2 é o \(f^2\) de Cohen, calculado a partir do \(R^2\) que foi fornecido. A implementação
desta função está de acordo com Cohen (1992).
Cuidados com RLS
Verificar se um 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:
# Anscombe.R
source("eiras.friendlycolor.R")
source("eiras.correg.R")
col <- friendlycolor(c(1,7,13,19))
pch <- 21:24
Anscombe <- datasets::anscombe
min.x <- 4 # min(c(Anscombe$x1,Anscombe$x2,Anscombe$x3,Anscombe$x4),na.rm=TRUE)
max.x <- 19 # max(c(Anscombe$x1,Anscombe$x2,Anscombe$x3,Anscombe$x4),na.rm=TRUE)
min.y <- 3 # min(c(Anscombe$y1,Anscombe$y2,Anscombe$y3,Anscombe$y4),na.rm=TRUE)
max.y <- 17 # max(c(Anscombe$y1,Anscombe$y2,Anscombe$y3,Anscombe$y4),na.rm=TRUE)
res <- list()
for (i in 1:4)
{
if (i==1)
{
x <- Anscombe$x1
y <- Anscombe$y1
}
if (i==2)
{
x <- Anscombe$x2
y <- Anscombe$y2
}
if (i==3)
{
x <- Anscombe$x3
y <- Anscombe$y3
}
if (i==4)
{
x <- Anscombe$x4
y <- Anscombe$y4
}
res <- c(res, correg (x, y,
alpha=0.05, B=0,
method = "lm",
jitter=0,
main=paste("conjunto ",i,sep=""),
xlab="x", ylab="y",
xlim=c(min.x,max.x),
ylim=c(min.y,max.y),
bg="transparent", col=col[i], pch=pch[i]))
}
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 . xEm todos os gráficos computamos as regressões lineares simples com as respectivas bandas de confiança.
Um modelo linear só é adequado para o primeiro conjunto, mas é
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 coeficiente angular e os valores p
associados são todos similares (demo_Anscombe2.R):
# Anscombe.R
cat("\n\n####--- media, d.p. e quartis")
for (r.idx in 0:3)
{
cat("\n\nConjunto ",r.idx+1,":\n")
print(res[[1+6*r.idx]])
}
cat("\n\n####--- retas de regressão (intercepto, coef. angular e estatistica)")
for (r.idx in 0:3)
{
cat("\n\nConjunto ",r.idx+1,":\n")
print(res[[5+6*r.idx]]$coefficients)
}
cat("\n\n####--- lambda (razao entre variancias)")
for (r.idx in 0:3)
{
cat("\n\nConjunto ",r.idx+1,":\n")
print(res[[3+6*r.idx]])
}
####--- media, 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
####--- retas de regressão (intercepto, coef. angular e estatistica)
Conjunto 1 :
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.0000909 1.1247468 2.667348 0.025734051
x 0.5000909 0.1179055 4.241455 0.002169629
Conjunto 2 :
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.000909 1.1253024 2.666758 0.025758941
x 0.500000 0.1179637 4.238590 0.002178816
Conjunto 3 :
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.0024545 1.1244812 2.670080 0.025619109
x 0.4997273 0.1178777 4.239372 0.002176305
Conjunto 4 :
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.0017273 1.1239211 2.670763 0.025590425
x 0.4999091 0.1178189 4.243028 0.002164602
####--- lambda (razao entre variancias)
Conjunto 1 :
[,1]
lambda "1.80086371549847"
status "Assuming:\n\t- Explanatory variable: x;\n\t- Dependent variable: y.\n"
Conjunto 2 :
[,1]
lambda "1.79950642693625"
status "Assuming:\n\t- Explanatory variable: x;\n\t- Dependent variable: y.\n"
Conjunto 3 :
[,1]
lambda "1.80347027989639"
status "Assuming:\n\t- Explanatory variable: x;\n\t- Dependent variable: y.\n"
Conjunto 4 :
[,1]
lambda "1.80441610212107"
status "Assuming:\n\t- Explanatory variable: x;\n\t- Dependent variable: y.\n"Atenção aos limites
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.
Por exemplo, vamos utilizar os dados disponíveis no pacote
MASS, contendo o valor de 506 residências nos subúrbios de
Boston, MA, Estados Unidos. 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 (Rscript em demo_Boston_ate20000.R):
source("eiras.friendlycolor.R")
source("eiras.correg.R")
# Load the data
data("Boston", package = "MASS")
x <- Boston$lstat[Boston$medv<=20]
y <- Boston$medv[Boston$medv<=20]
alfa <- 0.05
B <- 0
lst <- correg(x, y,
alpha=alfa, B=B,
method="lm_robust",
xlab="BSS (%)",
ylab="MedianaCasa (x1000 USD)",
bg="transparent", col=friendlycolor(14), pch=21)
Variable Mean Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1 BSS (%) 13.27 6.2885 29.93 17.1 20.45 10.26 8.47 215 0
2 MedianaCasa (x1000 USD) 18.9 3.7873 16.5 18.9 15 18.2 19.9 215 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 \times \text{lstat} \]
Como era de se esperar, a relação é negativa: quando maior a porcentagem de BSS, menor o valor das casas 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 das casas é
\[ \widehat{\text{medv}} = 22.24 - 0.38 \cdot 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 de classe social mais baixa?
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 mais baratas, sem perceber separamos os bairros mais pobres com 7.79% ou mais de pessoas de classe social baixa:
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 estimadora para um bairro com 1% de BSS. Caso
considerássemos todos os dados disponíveis (demo_Boston_todas.R), teríamos:
source("eiras.correg.R")
# Load the data
data("Boston", package = "MASS")
alfa <- 0.05
B <- 0
lst <- correg(Boston$lstat, Boston$medv,
alpha=alfa, B=B,
method="lm_robust",
xlab="BSS (%)",
ylab="MedianaCasa (x1000 USD)",
bg="transparent", col=friendlycolor(14), pch=21)
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. No entanto, esta segunda reta também não pode ser usada porque, claramente, a relação entre as duas variáveis não é linear.
Alavancagem
O terceiro conjunto de dados do Quarteto de Anscombe mostra este
efeito, 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.Rpara \((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.Rpara \((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")
Métodos robustos
Para os métodos mais robustos:
estimatr::lm_robustnão precisa da suposição de homocedasticidade.métodos por bootstrapping só requerem independência entre os pares de observações.
Sempre que possível, utilize método robusto.
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 (dos indivíduos nos 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. Recomenda adotá-la sempre no lugar de
lm.
Judkins, 2015, https://doi.org/10.1002/sim.6839
Para situações ainda mais difíceis, lancaremos mão de bootstrapping, que tolera todas as violações, com exceção de dependência entre os pares de observações.
Quanto à normalidade, verificaremos a aparência de seus density plots (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). É possível testá-la, assim como a binormalidade (algo como um chapéu mexicano em 3 dimensões), que é ainda mais difícil de abordar.
- hematócrito (
demo_HTnormal.R):
Gestantes <- readRDS("Gestante.rds")
densidade <- density(Gestantes$HT)
media <- mean(Gestantes$HT)
desvpad <- sd(Gestantes$HT)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Hematocrito (%)",
ylab="densidade",
type="l")
lines(x,y,lty=2)
- concentração de hemoglobina (
demo_HBnormal.R):
Gestantes <- readRDS("Gestante.rds")
densidade <- density(Gestantes$HB)
media <- mean(Gestantes$HB)
desvpad <- sd(Gestantes$HB)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Hemoglobina (mg/dl)",
ylab="densidade",
type="l")
lines(x,y,lty=2)
- contagem de hemácias (
demo_HEMnormal.R):
Gestantes <- readRDS("Gestante.rds")
densidade <- density(Gestantes$HEM)
media <- mean(Gestantes$HEM)
desvpad <- sd(Gestantes$HEM)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Hemacias (milhoes/mm^3)",
ylab="densidade",
type="l")
lines(x,y,lty=2)
- contagem de leucócitos (
demo_LEUCnormal.R):
Gestantes <- readRDS("Gestante.rds")
densidade <- density(Gestantes$LEUC)
media <- mean(Gestantes$LEUC)
desvpad <- sd(Gestantes$LEUC)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Leucocitos (milhares/mm^3)",
ylab="densidade",
type="l")
lines(x,y,lty=2)
|
Um teste satisfatório de binormalidade (normalidade bivariada) e de
normalidade são os de Henze-Zirkler e de Shapiro-Wilk, respectivamente.
Estão implementados em A saída é: Observe que a não rejeição da binormalidade não implica na não-rejeição das uninormalidades (e.g., HT x LEUC). Observe os valores p. As hipóteses nulas de binormalidade e a uninormalidade são rejeitadas em alguns casos. No entanto, como \(n\)=40, seguiremos com a regressão linear simples com base no teorema central do limite (o número de pares independentes é maior que 30). Pode acontecer o reverso, quando a binormalidade for rejeitada, mas para cada variável isoladamente não houver rejeição. Neste caso temos um situação que poderíamos chamar de binormalidade fraca. A binormalidade é condição suficiente, mas não necessária para a validade do teste da regressão. Para amostras de tamanho pequeno, pode executar a regressão com cautela, mas não há garantias. Uma alternativa é utilizar bootstrapping, que prescinde da suposição de binormalidade. |
Neste caso, vamos inicialmente assumir que as regressões lineares serão válidas, mesmo porque verificamos que, visualmente, as quatro variáveis têm distribuição, se não normal, razoavelmente simétricas.
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) ou
bootstrapping, 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 regressao populacional esteja contida.
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, há modelo estatisticamente válido 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 regressão linear simples.
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. A função mais robusta leva em conta quanto os valores originais divergem de uma distribuição normal.
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 <- 1e3
# 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: 1000 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 IC95.2
intercept 0.005796432 -3.0577195 -0.6748573 0.7786868 5.4884615 -1.8222610 2.703369
slope 0.348530741 0.1961538 0.3279310 0.3669662 0.4230878 0.2749798 0.395818
Equation:
mean[Hemoglobina (mg/dl)] = 0.00579643 [-1.82226098, 2.70336926] +
0.34853074 [0.27497976, 0.39581798] . 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: 1000 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 IC95.2
intercept 0.5001719 -1.1764859 0.08294684 1.0243454 3.4327640 -0.63836337 2.2435998
slope 0.1021272 0.0236311 0.08806780 0.1129421 0.1483878 0.05548873 0.1320752
Equation:
mean[Hemácia (milhão/mm³)] = 0.50017193 [-0.63836337, 2.24359982] +
0.10212717 [0.05548873, 0.13207515] . 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: 1000 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 IC95.2
intercept 5.48929940 -11.6814256 2.82248859 8.0821697 19.7148707 -2.8373180 13.7492447
slope 0.08678702 -0.2674569 0.01713941 0.1590913 0.5630927 -0.1214843 0.3141092
Equation:
mean[Leucócito (milhar/mm³)] = 5.4892994 [-2.837318, 13.74924475] +
0.08678702 [-0.12148429, 0.31410916] . Hematócrito (%)Esta modalidade robusta não reporta valores p. A decisão é tomada pelo intervalo de confiança: 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, aparece o valor 0 no intervalo.
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. A banda é a região dentro da qual espera, com confiança 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 que fique inteiramente contida na banda, 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 mais tempo 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
| Variáveis | Método | Coeficiente angular | Equação obtida |
|---|---|---|---|
| HT x 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 (%) | |
| Bootstrapping | 0.349 [0.275,0.396] | mean[Hemoglobina (mg/dl)] = 0.00579643 [-1.82226098, 2.70336926] + 0.34853074 [0.27497976, 0.39581798] . Hematócrito (%) | |
| HT x 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 (%) | |
| Bootstrapping | 0.102 [0.055,0.132] | mean[Hemácia (milhão/mm³)] = 0.50017193 [-0.63836337, 2.24359982] + 0.10212717 [0.05548873, 0.13207515] . Hematócrito (%) | |
| HT x 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 (%) | |
| Bootstrapping | 0.087 [-0.121,0.314] | mean[Leucócito (milhar/mm³)] = 5.4892994 [-2.837318, 13.74924475] + 0.08678702 [-0.12148429, 0.31410916] . Hematócrito (%) |
- A função
lmnão fornece o intervalo de confiança, mas pudemos obtê-los comconfint(lm(y~x)).
Segundo estudo: corpo e cérebro
Podemos usar o tamanho do corpo para prever o tamanho do cérebro em mamíferos a partir de alguns animais conhecidos?
Vamos utilizar os dados colecionados em CorpoCerebroVonBonin.xlsx (von Bonin,
1937).
Esta planilha tem a massa, em gramas, de corpos e cérebros de 119 animais. Adicionamos um Homo sapiens, que não fazia parte da tabela original.
A regressão linear dada por lm_robust mostra o
seguinte:
CorpoCerebro <- readxl::read_excel("CorpoCerebroVonBonin.xlsx")
modelo <- estimatr::lm_robust(CorpoCerebro$Cerebro ~ CorpoCerebro$Corpo)
sumario <- summary(modelo)
print(sumario)
Call:
estimatr::lm_robust(formula = CorpoCerebro$Cerebro ~ CorpoCerebro$Corpo)
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
(Intercept) 2.51e+02 6.437e+01 3.899 0.0001607 1.235e+02 3.785e+02
CorpoCerebro$Corpo 7.31e-05 2.030e-05 3.601 0.0004657 3.289e-05 1.133e-04
DF
(Intercept) 118
CorpoCerebro$Corpo 118
Multiple R-squared: 0.4637 , Adjusted R-squared: 0.4592
F-statistic: 12.96 on 1 and 118 DF, p-value: 0.0004657A saída textual mostra, aproximadamente:
- \(R^2\) = 46.4% e
- p = 4.66e-04
Portanto, para \(\alpha=0.05\), rejeitamos \(H_0\) e assumimos que existe modelo. 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 gramas, aproximadamente, com a reta de regressão:
\[ \widehat{\text{cérebro}} = 251 + 7.31 \cdot 10^{-5} \times \text{corpo} \]
Pode parecer que está tudo bem. 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 <- 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)
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)Os dados estão muito concentrados no lado inferior-esquerdo do
gráfico (difícil aferir linearidade). Há elefantes (com cérebros maiores
do que 3.5 kg) e baleias (com corpos maiores que 40 toneladas) afastados
dos demais; os valores parecem tornar mais dispersos com os corpos
maiores (sinal de heterocedasticidade). A banda de confiança 95% é uma
versão analítica e podemos ver boa parte dos dados longe dela, outro
sinal de má escolha do modelo linear 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 <- 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 torna 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 parece uma boa reta para predizer seus valores (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 <- 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)
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 IC95.1 IC95.2
intercept 2.386879e+02 4.756744e+01 1.934587e+02 2.849918e+02 4.673195e+02 1.061900e+02 3.795772e+02
slope 7.434957e-05 3.712535e-05 6.175337e-05 9.160721e-05 2.743132e-03 3.984389e-05 1.112845e-03
Equation:
mean[Cérebro (g)] = 238.68790591 [106.18996699, 379.57723026] +
7.435e-05 [3.984e-05, 0.00111284] . Corpo (g)Além disto, as distribuições dos pesos dos corpos e dos cérebros não
são aproximadamente normais. 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 <- readxl::read_excel("CorpoCerebroVonBonin.xlsx")
# 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")
Este exemplo, portanto, não atende às premissas da RLS clássica usada no exemplo anterior.
Transformação potência de Tukey
A transformação potência de Tukey de variável intervalar é uma transformação não linear que pode alterar 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 minimiza a 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 o peso do corpo (mostrando todas as transformações):
CorpoCerebro <- readxl::read_excel("CorpoCerebroVonBonin.xlsx")
source("eiras.tukey.try.R")
potencia <- tukey.try(CorpoCerebro$Corpo, xlab="Corpo (g)", show="all")
-------------------
Tukey's lambda = -3
-------------------
mean: -3.194846e-06
st.dev.: 1.776165e-05
median: -2.264176e-11
mode: 7.126102e-10
iqr: 9.2136e-10
Asymmetry coefficient: 3468.306
p < 1e-18
---------------------
Tukey's lambda = -2.5
---------------------
mean: -1.529461e-05
st.dev.: 8.134489e-05
median: -1.34606e-09
mode: 2.283494e-08
iqr: 2.952413e-08
Asymmetry coefficient: 518.8109
p < 1e-18
-------------------
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.5
---------------------
mean: -0.0004028596
st.dev.: 0.001748723
median: -4.757796e-06
mode: -6.292558e-07
iqr: 3.004904e-05
Asymmetry coefficient: 13.3858
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.5
---------------------
mean: -0.0297666
st.dev.: 0.0419449
median: -0.01681864
mode: -0.01104701
iqr: 0.02446507
Asymmetry coefficient: 0.7651561
p = 1.47e-11
------------------
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 = 0.5
--------------------
mean: 349.2362
st.dev.: 1225.994
median: 59.45932
mode: 34.44313
iqr: 116.5643
Asymmetry coefficient: -2.700596
p < 1e-18
------------------
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 = 1.5
--------------------
mean: 11879142737
st.dev.: 77007968634
median: 210229
mode: -2513543
iqr: 3249959
Asymmetry coefficient: -3655.941
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
--------------------
Tukey's lambda = 2.5
--------------------
mean: 7.453301e+17
st.dev.: 5.204421e+18
median: 743377129
mode: -56060540177
iqr: 72482735225
Asymmetry coefficient: -10282865
p < 1e-18
------------------
Tukey's lambda = 3
------------------
mean: 6.021923e+21
st.dev.: 4.336581e+22
median: 44206269206
mode: -8.335041e+12
iqr: 1.077668e+13
Asymmetry coefficient: -558791971
p < 1e-18
-------------------
Best transformation
-------------------
Tukey's lambda = 0
(logarithm transformation)
Asymmetry coefficient = -0.17784156613391
p = 0.0180
equation: ln[Corpo (g)]
deve ser transformado por logaritmo.
- para o peso do cérebro (mostrando só a melhor escolha ao final):
CorpoCerebro <- 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)]
também deve ser transformado por logaritmo.
Então executamos a regressão implementada em CorpoCerebro_RLS_log.R, que cria duas
colunas adicionais com a transformação sugerida:
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.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 <- readxl::read_excel("CorpoCerebroVonBonin.xlsx")
# transformacao log
CorpoCerebro$Cerebro_log <- log(CorpoCerebro$Cerebro)
CorpoCerebro$Corpo_log <- log(CorpoCerebro$Corpo)
lst <- correg(CorpoCerebro$Corpo_log, CorpoCerebro$Cerebro_log,
alpha=alfa, B=B,
method="lm_robust",
xlab="Ln[Corpo (g)]", ylab="Ln[Cérebro (g)]",
bg="transparent", col=col_CxC, pch=pch_CxC)
Variable Mean Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1 Ln[Corpo (g)] 9.158 2.9119 11.0572 10.0008 11.4134 8.63550941823428 8.60410456340553 120 0
2 Ln[Cérebro (g)] 4.8675 2.0633 5.9915 5.8435 6.0521 4.57779898919196 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)]Observe que os dados estão muito mais bem arranjados no gráfico, mais uniformemente distribuídos ao longo da reta. Para \(\alpha=0.05\) o modelo é estatisticamente válido (\(p<2.2 \cdot 10^{-16}\)). O coeficiente de determinação é \(R^2 \approx 85.7\%\) (grande). A equação é (ambos, cérebro e corpo em gramas):
\[ \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} \]
Interessantemente, no trabalho original de 1936/1937, Gerhardt von Bonin utilizando apenas papel e paciência chegou à expressão:
\[ \widehat{\text{Cérebro}} = 0.18 \times {\text{Corpo}}^{0.655} \]
Este resultado é um feito notável!
|
Outro feito notável é
|
Outro aspecto interessante é que, usando uma regressão linear, 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 <- 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:
# A tibble: 3 × 4
Grupo Animal Cerebro Corpo
<chr> <chr> <dbl> <dbl>
1 Cetacea Megaptera boops 3531 42372000
2 Cetacea Balaenoptera rostrata 2490 62250000
3 Cetacea Balaenoptera rostrata 7000 73997000Sã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,])# A tibble: 3 × 4
Grupo Animal Cerebro Corpo
<chr> <chr> <dbl> <dbl>
1 Ungulata Elephas africanus 4370 1642000
2 Ungulata Elephas indicus 5045 2547000
3 Cetacea Balaenoptera musculus 3636 5090400São os elefantes africano e indiano, e a baleia azul:
E qual é este com corpo acima de uma tolenada, mas cérebro abaixo de 1 kg?
# A tibble: 1 × 4
Grupo Animal Cerebro Corpo
<chr> <chr> <dbl> <dbl>
1 Ungulata Hippopotamus amphibiu IS 582 582 1749730Interessantemente, 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,])# A tibble: 2 × 4
Grupo Animal Cerebro Corpo
<chr> <chr> <dbl> <dbl>
1 Cetacea Lagenorrhynchus albirostris 1126 67560
2 Primates Homo sapiens 1320 62000Dois animais: um cetáceo e o ser humano. Com uma rápida pesquisa no Google encontramos um golfinho:
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,])# A tibble: 1 × 4
Grupo Animal Cerebro Corpo
<chr> <chr> <dbl> <dbl>
1 Cetacea Tursiops tursio 1886 278000Outro golfinho:
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,])# A tibble: 6 × 4
Grupo Animal Cerebro Corpo
<chr> <chr> <dbl> <dbl>
1 Primates Anthropopithecus troblodytes 345 22045
2 Primates Papio doguerra 182 4990
3 Primates Papio anubis 213. 6810
4 Primates Papio sphinx 174. 4763
5 Primates Papio bahuin 180. 10546
6 Primates Papio spec. 213 22220Encontramos, 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.
Revisitando a correlação
Diferença entre correlação e regressão
A correlação é sempre confundida com a regressão linear simples. 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 a retas de regressão
horizontais, i.e., \(\hat{\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: a correlação mede o grau de linearidade entre duas variáveis, sem direção e sem dimensão, enquanto a regressão linear simples pretende obter uma equação para, a partir do valor de uma variável explicativa (ou independente), estimar o valor médio de uma variável de desfecho (ou dependente), 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 é computada com as variáveis padronizadas. Como o procedimento de padronização elimina as unidades de medida, dizer que a correlação utiliza as variáveis padronizadas implica que a correlação é adimensional. A regressão linear simples, ao contrário, busca 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 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 círculo quando \(r=0\). Todas as elipses têm eixo principal com inclinação próxima a 45o 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 45o.
No entanto, na seção padronização das variáveis encontramos para a regressão linear simples com as variáveis padronizadas os valores \(\hat{\beta_1} = r\) (a inclinação da reta é o coeficiente angular) e \(\hat{\beta_0} = 0\) (o intercepto é igual a zero). 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 45o. Que retas e que ângulos são estes? O entendimento destes comportamentos da correlação e da regressão linear simples é 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 traz a estatura e a massa corpórea relatadas por 89 estudantes
de uma classe de graduação de uma faculdade de administração em 2008.
Uma relação linear pode ser considerada? 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 <- readxl::read_excel("Adm2008.xlsx")
lst <- correg(Dados$Estatura, Dados$MCT,
method="preview",
alpha=alfa, B=B,
xlab="Estatura (cm)", ylab="Massa (kg)",
col=col, bg="transparent", pch=pch)
Variable Mean Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1 Estatura (cm) 1.57 0.0967 1.61 1.56 1.72 1.68 1.69 89 0
2 Massa (kg) 44 14.0235 53 50 60 57 60 89 0
lambda = 0.00542906452245859
Assuming:
- Explanatory variable: Estatura (cm);
- Dependent variable: Massa (kg).Admitindo que sim,
- a correlação entre estas duas variáveis (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 <- readxl::read_excel("Adm2008.xlsx")
lst <- correg(Dados$Estatura, Dados$MCT,
method="pearson",
alpha=alfa, B=B,
xlab="Estatura (cm)", ylab="Massa (kg)",
col=col, bg="transparent", pch=pch)
Variable Mean Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1 Estatura (cm) 1.57 0.0967 1.61 1.56 1.72 1.68 1.69 89 0
2 Massa (kg) 44 14.0235 53 50 60 57 60 89 0
lambda = 0.00542906452245859
Assuming:
- Explanatory variable: Estatura (cm);
- Dependent variable: Massa (kg).
-----------
Correlation
-----------
Pearson's product-moment correlation
data: Estatura (cm) and Massa (kg)
t = 12.029, df = 87, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6966559 0.8574067
sample estimates:
cor
0.7902592 para o qual computamos o coeficiente de correlação igual a r=0.79.
- a regressão linear (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 <- readxl::read_excel("Adm2008.xlsx")
lst <- correg(Dados$Estatura, Dados$MCT,
method="lm_robust",
alpha=alfa, B=B,
xlab="Estatura (cm)", ylab="Massa (kg)",
col=col, bg="transparent", pch=pch)
Variable Mean Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1 Estatura (cm) 1.57 0.0967 1.61 1.56 1.72 1.68 1.69 89 0
2 Massa (kg) 44 14.0235 53 50 60 57 60 89 0
lambda = 0.00542906452245859
Assuming:
- Explanatory variable: Estatura (cm);
- Dependent variable: Massa (kg).
------------------------
Simple linear regression
------------------------
Call:
"lm_robust(formula = Massa (kg) ~ Estatura (cm))"
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) -130.2 11.751 -11.08 2.635e-18 -153.5 -106.8 87
Estatura (cm) 114.6 7.111 16.12 7.334e-28 100.5 128.8 87
Multiple R-squared: 0.6245 , Adjusted R-squared: 0.6202
F-statistic: 259.8 on 1 and 87 DF, p-value: < 2.2e-16
Equation:
mean[Massa (kg)]=-130.17999098 + 114.63475279 . Estatura (cm)\(\widehat{\text{massa}}\, \text{(kg)} =\) -130.18 + 114.63 \(~ \times ~ \text{estatura (cm)}\)
- e a regressão linear com as variáveis padronizadas (implementada em
Adm2008_RLS_z.R) é:
# Biometria_rPearson.R
source("eiras.friendlycolor.R")
source("eiras.correg.R")
alfa <- 0.05
B <- 0
col <- friendlycolor(7) # azul
pch <- 24
Dados <- readxl::read_excel("Adm2008.xlsx")
lst <- correg(Dados$Estatura, Dados$MCT,
method="lm_robust",
standardize=TRUE,
alpha=alfa, B=B,
xlab="Estatura (cm)", ylab="Massa (kg)",
col=col, bg="transparent", pch=pch)
Variable Mean Sd Min. 1st Qu. Median 3rd Qu. Max. n NA's
1 Estatura (cm) -1.4621 1 -1.0483 -1.5656 0.0895 -0.32426785231236 -0.220827569675084 89 0
2 Massa (kg) -1.5688 1 -0.927 -1.1409 -0.4279 -0.641778990570569 -0.42785266038038 89 0
lambda = 2.26579971435054
Assuming:
- Explanatory variable: Estatura (cm);
- Dependent variable: Massa (kg).
------------------------
Simple linear regression
------------------------
Call:
"lm_robust(formula = Massa (kg) ~ Estatura (cm))"
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 5.462e-16 0.06516 8.383e-15 1.000e+00 -0.1295 0.1295 87
Estatura (cm) 7.903e-01 0.04902 1.612e+01 7.334e-28 0.6928 0.8877 87
Multiple R-squared: 0.6245 , Adjusted R-squared: 0.6202
F-statistic: 259.8 on 1 and 87 DF, p-value: < 2.2e-16
Equation:
mean[Massa (kg)]=0 + 0.79025925 . Estatura (cm)\(\widehat{\text{massa}}(z) =\) 5.462e-16 + 0.79 \(~ \times ~ \text{estatura(z)}\)
Exceto por arredondamento, a reta calculada comprova o que esperávamos da teoria: é praticamente igual a
\[ \widehat{\text{massa}}(z) = 0 + r \times \text{estatura(z)} \]
É difícil visualizar a elipse neste exemplo por causa da distorção da escala. Então, vamos exibir os dados padronizados com os eixos utilizando escalas idênticas:
col <- friendlycolor(2) # violeta
pch <- 23
Dados <- readxl::read_excel("Adm2008.xlsx")
# padroniza os valores
Dados$Estatura <- scale(Dados$Estatura)
Dados$MCT <- scale(Dados$MCT)
plot(add.jitter(Dados$Estatura),add.jitter(Dados$MCT),
xlim=c(-4,4), ylim=c(-4,4),
xlab="Estatura (z)", ylab="Massa (z)",
pch=pch,col=col,bg=col,cex=0.6)
# bissetriz
lines(c(-4,4),c(-4,4),lty=2)
Note que a dispersão dos pontos é próxima a \(r=0.8\) 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 (coeficiente angular igual a \(1\)), como nos gráficos que observamos na seção de padronização das variáveis. No entanto, o valor de r \(\approx\) 0.79 desvia da bissetriz, como teorizamos.
Para verificar como a reta de regressão se relaciona com a correlação, adicionamos ao código anterior a regressão linear simples dos dados gerados e a plotamos:
plot(add.jitter(Dados$Estatura),add.jitter(Dados$MCT),
xlim=c(-4,4), ylim=c(-4,4),
xlab="Estatura (z)", ylab="Massa (z)",
pch=pch,col=col,bg=col,cex=0.6)
# bissetriz
lines(c(-4,4),c(-4,4),lty=2)
# regressao linear y~x dos valores gerados
rls <- lm_robust(Dados$MCT~Dados$Estatura)
xhat <- seq(min(Dados$Estatura,na.rm=TRUE),
max(Dados$Estatura,na.rm=TRUE), length.out = 10)
yhat <- rls$coefficients[1] + rls$coefficients[2]*xhat
lines(xhat,yhat,lwd=4) 
Para melhorar a visualização, podemos adicionar a elipse de confiança 90% (poderia ser qualquer outra) com seu eixo principal e duas verticais pontilhadas em seus limites, à esquerda e à direita:
plot(add.jitter(Dados$Estatura),add.jitter(Dados$MCT),
xlim=c(-4,4), ylim=c(-4,4),
xlab="Estatura (z)", ylab="Massa (z)",
pch=pch,col=col,bg=col,cex=0.6)
# bissetriz
lines(c(-4,4),c(-4,4),lty=2)
# regressao linear y~x dos valores gerados
lines(xhat,yhat,lwd=4)
ellipse <- ellipseaxis(x=Dados$Estatura,
y=Dados$MCT,
draw=FALSE, col=col, level=0.9)
lines(ellipse, col=col, lwd=2,lty=2)
abline(v=min(ellipse[,1]),col=col, lwd=2,lty=2)
abline(v=max(ellipse[,1]),col=col, lwd=2,lty=2)
É 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 é adirecional, poderíamos
computar uma segunda reta de regressão, invertendo as variáveis, 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 regressão linear simples, mostrando as duas retas de
regressão y~x e x~y tocando os pontos extremos
da elipse de confiança 95%.
|
Exemplo de aplicação
Padronizar as variáveis e proceder à regressão linear não é apenas
uma curiosidade. Pode ajudar em explorar melhor as relações. Observe o
que acontece com os dados de Adm2008.xlsx quando exibimos
separadamente as medidas de estatura e massa corpórea total (MCT) de
homens e mulheres obtidos entre estudantes do curso de graduação de
Administração da Faculdade de Economia e Administração da Universidade
de São Paulo.
Dados <- readxl::read_excel("Adm2008.xlsx")
DadosF <- Dados[Dados$Sexo=="Feminino",]
DadosM <- Dados[Dados$Sexo=="Masculino",]
nF <- min(sum(!is.na(DadosF$Estatura)),
sum(!is.na(DadosF$MCT)))
nM <- min(sum(!is.na(DadosM$Estatura)),
sum(!is.na(DadosM$MCT)))
# Elipse de predicao 95%
matriz <- as.matrix(Dados[, 3:4])
n <- nrow(matriz)
car::dataEllipse(matriz[,1], matriz[,2],
groups=factor(Dados$Sexo),
group.labels=c("F", "M"),
levels=c(.95,.999),
robust=TRUE,
main=paste("Elipses de predicao de 95% e 99.9%\n",
"n = ",n," (Fem.=",nF,", Masc.=",nM,")",sep=""),
xlab="Estatura (m)",
ylab="MCT (kg)",
xlim=c(1.3, 2.1),
ylim=c(25, 120),
lwd=0.8, lty=2)
----------------------
Estatística descritiva
----------------------
---------------
- sexo feminino
--------------- Estatura MCT
Min. :1.500 Min. :43.00
1st Qu.:1.600 1st Qu.:50.25
Median :1.640 Median :53.50
Mean :1.641 Mean :54.63
3rd Qu.:1.698 3rd Qu.:59.75
Max. :1.780 Max. :70.00
----------------
- sexo masculino
---------------- Estatura MCT
Min. :1.560 Min. : 48.00
1st Qu.:1.720 1st Qu.: 66.00
Median :1.760 Median : 73.00
Mean :1.764 Mean : 74.47
3rd Qu.:1.820 3rd Qu.: 82.00
Max. :1.930 Max. :105.00 As elipses externas (confiança de 99.9%) mostram que não há outliers se considerarmos que as distribuições dos grupos são binormais.
A elipse tinha uma aparência estranha porque há dois grupos misturados. Pelas localização e inclinação das elipses parece que os homens são mais altos e mais pesados que as mulheres, e também são maiores suas variabilidades em estatura e em peso.
Será que a relação entre estatura (m) e peso (kg) é a mesma em mulheres e homens?
Há duas formas de observar o comportamento das variáveis.
Regressões lineares simples (para cada grupo)
# split F & M
DadosF <- Dados[Dados$Sexo=="Feminino",]
DadosM <- Dados[Dados$Sexo=="Masculino",]
# regressao
rlsF <- estimatr::lm_robust(DadosF$MCT ~ DadosF$Estatura)
xF <- seq(min(DadosF$Estatura,na.rm=TRUE),
max(DadosF$Estatura,na.rm=TRUE),length.out=10)
yF <- rlsF$coefficients[1] + rlsF$coefficients[2]*xF
rlsM <- estimatr::lm_robust(DadosM$MCT ~ DadosM$Estatura)
xM <- seq(min(DadosM$Estatura,na.rm=TRUE),
max(DadosM$Estatura,na.rm=TRUE),length.out=10)
yM <- rlsM$coefficients[1] + rlsM$coefficients[2]*xM
plot(NA,
main="Regressão Linear Simples",
xlab="Estatura (m)",
ylab="MCT (kg)",
xlim=c(1.3, 2.1),
ylim=c(25, 120))
text(min(xF),min(yF),col="black",pos=2,"F")
points(DadosF$Estatura,DadosF$MCT,col="black",pch=1)
lines(xF,yF,col="black",lwd=2)
text(max(xM),max(yM),col="blue",pos=3,"M")
points(DadosM$Estatura,DadosM$MCT,col="blue",pch=2)
lines(xM,yM,col="blue",lwd=2)
As retas de regressão linear acompanham o perfil das elipses, com inclinação aparentemente maior para os homens. Então, para os homens, a variação em estatura corresponde a uma variação maior na massa corpórea do que a observada entre as mulheres.
Correlações (para cada grupo)
Será que a correlação entre estatura e massa corpórea é diferente para os dois sexos?
Considerando que são membros da mesma espécie, não deveríamos observar a mesma correlação?
Podemos, para responder a esta questão, comparar as correlações.
Para a parte gráfica, pelo que vimos acima, é mais coerente expressar graficamente com as variáveis padronizadas (lembre: não é apropriado traçar retas para a correlação; retas são para RLS). Assim, observamos:
# padronizacao
Dados$Estatura_z <- NA
Dados$MCT_z <- NA
DadosF$Estatura_z <- scale(DadosF$Estatura)
Dados$Estatura_z[Dados$Sexo=="Feminino"] <- scale(Dados$Estatura[Dados$Sexo=="Feminino"])
DadosF$MCT_z <- scale(DadosF$MCT)
Dados$MCT_z[Dados$Sexo=="Feminino"] <- scale(Dados$MCT[Dados$Sexo=="Feminino"])
DadosM$Estatura_z <- scale(DadosM$Estatura)
Dados$Estatura_z[Dados$Sexo=="Masculino"] <- scale(Dados$Estatura[Dados$Sexo=="Masculino"])
DadosM$MCT_z <- scale(DadosM$MCT)
Dados$MCT_z[Dados$Sexo=="Masculino"] <- scale(Dados$MCT[Dados$Sexo=="Masculino"])car::dataEllipse(Dados$Estatura_z, Dados$MCT_z,
groups=factor(Dados$Sexo),
group.labels=c("F", "M"),
levels=c(.95),
robust=TRUE,
main=paste("Elipses predição 95% e RLS padronizadas\n",
"n = ",n," (Fem=",nF,", Masc=",nM,")",sep=""),
xlab="Estatura (z)",
ylab="MCT (z)",
xlim=c(-4,4),
ylim=c(-4,4))
# RLS padronizada
rlsF <- estimatr::lm_robust(DadosF$MCT_z ~ DadosF$Estatura_z)
xF <- seq(min(DadosF$Estatura_z,na.rm=TRUE),
max(DadosF$Estatura_z,na.rm=TRUE),length.out=10)
yF <- rlsF$coefficients[1] + rlsF$coefficients[2]*xF
lines(xF,yF,col="black",lwd=2)
rlsM <- estimatr::lm_robust(DadosM$MCT_z ~ DadosM$Estatura_z)
xM <- seq(min(DadosM$Estatura_z,na.rm=TRUE),
max(DadosM$Estatura_z,na.rm=TRUE),length.out=10)
yM <- rlsM$coefficients[1] + rlsM$coefficients[2]*xM
lines(xM,yM,col="blue",lwd=2)
Visualmente, a regressão linear simples da massa corporal em função da estatura parece diferente e a correlação parece igual para mulheres e homens. Esta observação é acompanhada pelas elipses de predição 95% bastante coincidentes e retas de regressão com as variáveis padronizadas sobrepostas.
Implementamos, com testes estatísticos apropriados em demo_estatmct.R,
resultando:
# demo_estatmct.R
source("eiras.bartitle.R")
source("eiras.friendlycolor.R")
source("eiras.ConfidenceBand.R")
alfa <- 0.05
B <- 1e3
colM <- "#3a5fcd" # royalblue3
colH <- "#cd6600" # darkorange2
df_adm <- readxl::read_excel("Adm2008.xlsx")
df_admF <- df_adm[df_adm$Sexo=="Feminino",]
df_admM <- df_adm[df_adm$Sexo=="Masculino",]
nF <- nrow(df_admF)
nM <- nrow(df_admM)
# Elipse de predicao 95%
matriz <- as.matrix(df_adm[, 3:4])
n <- nrow(matriz)
car::dataEllipse(matriz[,1], matriz[,2],
groups=factor(df_adm$Sexo),
group.labels=c("Feminino", "Masculino"),
col=c("darkorange3","royalblue3"),
levels=c(0.95),
robust=TRUE,
main=paste("Elipse de predição de 95%\n",
"n = ",n," (Fem=",nF,", Masc=",nM,")",sep=""),
xlab="Estatura (m)",
ylab="MCT (kg)",
xlim=c(1.3, 2.1),
ylim=c(25, 120),
lwd=0.8, lty=2
)
# Correlacao, teste
correlF <- cor.test(df_admF$Estatura,df_admF$MCT)
correlM <- cor.test(df_admM$Estatura,df_admM$MCT)
crl <- psych::r.test(n=nF, n2=nM, correlF$estimate, correlM$estimate)
# RLS
# regressao
rlsF <- estimatr::lm_robust(df_admF$MCT ~ df_admF$Estatura)
xF <- seq(min(df_admF$Estatura,na.rm=TRUE),
max(df_admF$Estatura,na.rm=TRUE),length.out=10)
yF <- rlsF$coefficients[1] + rlsF$coefficients[2]*xF
rlsM <- estimatr::lm_robust(df_admM$MCT ~ df_admM$Estatura)
xM <- seq(min(df_admM$Estatura,na.rm=TRUE),
max(df_admM$Estatura,na.rm=TRUE),length.out=10)
yM <- rlsM$coefficients[1] + rlsM$coefficients[2]*xM
plot(NA,
main="Regressão Linear Simples",
xlab="Estatura (m)",
ylab="MCT (kg)",
xlim=c(1.3, 2.1),
ylim=c(25, 120))
text(min(xF),min(yF),col=colH,pos=2,"Feminino")
points(df_admF$Estatura,df_admF$MCT,col=paste(colH,"60",sep=""),pch=1)
lines(xF,yF,col=colH,lwd=2)
text(max(xM),max(yM),col=colM,pos=3,"Masculino")
points(df_admM$Estatura,df_admM$MCT,col=paste(colM,"60",sep=""),pch=2)
lines(xM,yM,col=colM,lwd=2)
# band
lst <- ConfidenceBand(x=df_admF$Estatura,y=df_admF$MCT, alpha=alfa, B=B)
bandF <- lst[[2]]
lst <- ConfidenceBand(x=df_admM$Estatura,y=df_admM$MCT, alpha=alfa, B=B)
bandM <- lst[[2]]
lines(bandF$X, bandF$LB, col=colH, lty=2, lwd=1.3)
lines(bandF$X, bandF$UB, col=colH, lty=2, lwd=1.3)
lines(bandM$X, bandM$LB, col=colM, lty=2, lwd=1.3)
lines(bandM$X, bandM$UB, col=colM, lty=2, lwd=1.3)
# RLS, test (p.value da interacao, paralelismo)
df_adm$SexoNum <- NA
df_adm$SexoNum[df_adm$Sexo=="Feminino"] <- 0
df_adm$SexoNum[df_adm$Sexo=="Masculino"] <- 1
rls_paralelismo <- estimatr::lm_robust(df_adm$MCT ~
df_adm$Estatura+df_adm$SexoNum+
df_adm$Estatura:df_adm$SexoNum)
# RLS, test (p.value do fator, coincidencia)
rls_coincidencia <- estimatr::lm_robust(df_adm$MCT ~
df_adm$Estatura+df_adm$SexoNum)
# padronizacao
df_adm$Estatura_z <- NA
df_adm$MCT_z <- NA
m <- mean(df_admF$Estatura,na.rm=TRUE)
s <- sd(df_admF$Estatura,na.rm=TRUE)
df_admF$Estatura_z <- (df_admF$Estatura-m)/s
df_adm$Estatura_z[df_adm$Sexo=="Feminino"] <- (df_adm$Estatura[df_adm$Sexo=="Feminino"]-m)/s
m <- mean(df_admF$MCT,na.rm=TRUE)
s <- sd(df_admF$MCT,na.rm=TRUE)
df_admF$MCT_z <- (df_admF$MCT-m)/s
df_adm$MCT_z[df_adm$Sexo=="Feminino"] <- (df_adm$MCT[df_adm$Sexo=="Feminino"]-m)/s
m <- mean(df_admM$Estatura,na.rm=TRUE)
s <- sd(df_admM$Estatura,na.rm=TRUE)
df_admM$Estatura_z <- (df_admM$Estatura-m)/s
df_adm$Estatura_z[df_adm$Sexo=="Masculino"] <- (df_adm$Estatura[df_adm$Sexo=="Masculino"]-m)/s
m <- mean(df_admM$MCT,na.rm=TRUE)
s <- sd(df_admM$MCT,na.rm=TRUE)
df_admM$MCT_z <- (df_admM$MCT-m)/s
df_adm$MCT_z[df_adm$Sexo=="Masculino"] <- (df_adm$MCT[df_adm$Sexo=="Masculino"]-m)/s
car::dataEllipse(df_adm$Estatura_z, df_adm$MCT_z,
groups=factor(df_adm$Sexo),
group.labels=c("Feminino", "Masculino"),
col=c(colH,colM),
levels=c(1-alfa),
robust=TRUE,
main=paste("Elipses predição 95% e RLS padronizadas\n",
"n = ",n," (Fem=",nF,", Masc=",nM,")",sep=""),
xlab="Estatura (z)",
ylab="MCT (z)",
xlim=c(-4,4),
ylim=c(-4,4))
# RLS padronizada
rlsF <- estimatr::lm_robust(df_admF$MCT_z ~ df_admF$Estatura_z)
xF <- seq(min(df_admF$Estatura_z,na.rm=TRUE),max(df_admF$Estatura_z,na.rm=TRUE),length.out=10)
yF <- rlsF$coefficients[1] + rlsF$coefficients[2]*xF
lines(xF,yF,col=colH,lwd=2)
rlsM <- estimatr::lm_robust(df_admM$MCT_z ~ df_admM$Estatura_z)
xM <- seq(min(df_admM$Estatura_z,na.rm=TRUE),max(df_admM$Estatura_z,na.rm=TRUE),length.out=10)
yM <- rlsM$coefficients[1] + rlsM$coefficients[2]*xM
lines(xM,yM,col=colM,lwd=2)
# saida numericaS
cat(bartitle("Estatística descritiva"))
cat(bartitle("- sexo feminino"))
print(summary(df_adm[df_adm$Sexo=="Feminino",3:4]))
cat(bartitle("- sexo masculino"))
print(summary(df_adm[df_adm$Sexo=="Masculino",3:4]))
cat(bartitle("RLS"))
cat(bartitle("- sexo feminino"))
print(summary(rlsF))
cat(paste("\nEquation: media[MCT(kg)] = ",round(rlsF$coefficients[1],4)," + ",
round(rlsF$coefficients[2],4)," * Estatura(m)\n",sep=""))
cat(bartitle("- sexo masculino"))
print(summary(rlsM))
cat(paste("\nEquation: media[MCT(kg)] = ",round(rlsM$coefficients[1],4)," + ",
round(rlsM$coefficients[2],4)," * Estatura(m)\n",sep=""))
cat(bartitle("Correlação"))
cat(bartitle("- sexo feminino"))
print(correlF)
cat(bartitle("- sexo masculino"))
print(correlM)
cat(bartitle("Paralelismo das RLS"))
# print(rls_paralelismo)
cat("H0: retas populacionais com inclinações iguais\n",sep="")
cat("\tp = ",rls_paralelismo$p.value[4],sep="")
cat(bartitle("Coincidência das RLS"))
# print(rls_coincidencia)
cat("H0: retas populacionais com interceptos iguais\n",sep="")
cat("\tp = ",rls_coincidencia$p.value[3],sep="")
cat(bartitle("Correlação"))
# print(rls_coincidencia)
cat("H0: correlações de Pearson populacionais iguais\n",sep="")
cat("\tp = ",crl$p,sep="")
----------------------
Estatística descritiva
----------------------
---------------
- sexo feminino
---------------
Estatura MCT
Min. :1.500 Min. :43.00
1st Qu.:1.600 1st Qu.:50.25
Median :1.640 Median :53.50
Mean :1.641 Mean :54.63
3rd Qu.:1.698 3rd Qu.:59.75
Max. :1.780 Max. :70.00
----------------
- sexo masculino
----------------
Estatura MCT
Min. :1.560 Min. : 48.00
1st Qu.:1.720 1st Qu.: 66.00
Median :1.760 Median : 73.00
Mean :1.764 Mean : 74.47
3rd Qu.:1.820 3rd Qu.: 82.00
Max. :1.930 Max. :105.00
---
RLS
---
---------------
- sexo feminino
---------------
Call:
estimatr::lm_robust(formula = df_admF$MCT_z ~ df_admF$Estatura_z)
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
(Intercept) 1.349e-15 0.1276 1.057e-14 1.0000000 -0.2588 0.2588
df_admF$Estatura_z 6.371e-01 0.1570 4.059e+00 0.0002536 0.3188 0.9555
DF
(Intercept) 36
df_admF$Estatura_z 36
Multiple R-squared: 0.406 , Adjusted R-squared: 0.3895
F-statistic: 16.48 on 1 and 36 DF, p-value: 0.0002536
Equation: media[MCT(kg)] = 0 + 0.6371 * Estatura(m)
----------------
- sexo masculino
----------------
Call:
estimatr::lm_robust(formula = df_admM$MCT_z ~ df_admM$Estatura_z)
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
(Intercept) -6.033e-16 0.10739 -5.618e-15 1.000e+00 -0.2158 0.2158
df_admM$Estatura_z 6.461e-01 0.07428 8.698e+00 1.669e-11 0.4968 0.7953
DF
(Intercept) 49
df_admM$Estatura_z 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: media[MCT(kg)] = 0 + 0.6461 * Estatura(m)
----------
Correlação
----------
---------------
- sexo 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
95 percent confidence interval:
0.3986675 0.7949180
sample estimates:
cor
0.637148
----------------
- sexo 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
95 percent confidence interval:
0.4507278 0.7823524
sample estimates:
cor
0.646062
-------------------
Paralelismo das RLS
-------------------
H0: retas populacionais com inclinações iguais
p = 0.03505788
--------------------
Coincidência das RLS
--------------------
H0: retas populacionais com interceptos iguais
p = 1.327539e-05
----------
Correlação
----------
H0: correlações de Pearson populacionais iguais
p = 0.9456523Podemos observar os dois últimos gráficos. O teste de duas correlações independentes foi feito no capítulo de Análise de Correlação. Também é possível testar se as duas regressões são iguais, i.e. paralelas e coincidentes. Os testes aparecem no final da saída textual.
Concluímos, estatisticamente, que estatura e massa corpórea para mulheres e homens são, portanto, igualmente correlacionadas, mas as respectivas regressões lineares simples diferem. São dois teste hierárquicos: o primeiro verifica se as retas são paralelas e o segundo, caso sejam paralelas, verifica se os interceptos são 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 no domínio comum aos dois grupos:
As bandas de confiança 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 homens e mulheres 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) foi rejeitada.
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.84Em 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 scaleTodos 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.3356De 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): rejectedO 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
- 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.
- Cohen, J (1992) A power primer. Quantitative methods in Psychology 112(1): 155-159.
- 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
- 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
- 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.
- 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.
- Zhang, Z & Yuan, K-H (2018). Practical statistical power analysis using WebPower and R (Eds). Granger, IN: ISDSA Press.