Fundamentos de Inferência Estatística em R
MSP4112: Estatística Aplicada à Medicina
21 agosto 2026 20:52h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(conf, warn.conflicts=FALSE))
suppressMessages(library(confintr, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(EnvStats, warn.conflicts=FALSE))
suppressMessages(library(epiR, warn.conflicts=TRUE))
suppressMessages(library(fmsb, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(gplots, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(labelled, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(openxlsx, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(Rmisc, warn.conflicts=FALSE))
suppressMessages(library(RVAideMemoire, warn.conflicts=FALSE))
source("eiras.friendlycolor.R")
source("eiras.create.population.R")
source("eiras.plot.density.withmeansd.R")
source("eiras.sampling.R")
source("summarySEwithin2.R")Código R
demo_binomial.Rdemo_binomialChocolate.Rdemo_binomialCondicoes.Rdemo_binomialDiff.Rdemo_binomialMethods.Reiras.create.population.Reiras.exit.Reiras.friendlycolor.Reiras.plot.barmeansd.Reiras.plot.density.empty.Reiras.plot.density.withmeansd.Reiras.sampling.Reiras.text.leading.Rfi_Amostra.Rfi_boot02.Rfi_boot02b.Rfi_bootstrapping.Rfi_bootstrapping2.Rfi_bootstrapping3.Rfi_criapopsample.Rfi_IntervaloConfiancaAnalitico.R`- fi_IntervaloPredicao.R`
fi_IntervalosDadosEntre.Rfi_IntervalosDadosIntra.Rfi_Normal.Rfi_Normal_e_Amostra.Rfi_ToDoMargemErro.RsimulaIC95pivotal.RsimulaIC95pivotal_Dieta.RsimulaIC95pivotal_MF.RsimulaIC95pivotal2.R
Arquivo de dados
Material
Conteúdo
- Amostragem
- Intervalo de predição
- Intervalo de confiança
- Distribuição normal
- Teorema central do limite
- Reamostragem (Bootstrapping)
- Intervalo de confiança de proporção
- Intervalo de credibilidade bayesiano
Amostragem e inferência estatística
Intervalo de predição de variável
Descrição da amostra
O intervalo simétrico centrado na média com amplitude de 4 desvios-padrão contém pelo menos 75% das observações da amostra.
Desigualdade de Chebyshev
\[ \text{IP}^{75\%}(X) = \left[\bar{x} \pm 2s\right] \]
- \(X\): variável intervalar
- \(\bar{x}\): média amostral
- \(s\): desvio-padrão amostral
Indício de distribuição normal
A distribuição normal contém aproximadamente 95% dos valores da amostra no intervalo simétrico centrado na média com amplitude de 4 desvios-padrão.
Ler Biometria_FMUSP.rds
Os dados do arquivo Biometria_FMUSP.rds foram coletados
pelos docentes dos estudantes do segundo ano de uma mesma disciplina do
curso de Medicina da FMUSP em três anos consecutivos.
As variáveis do arquivo são:
- ID: idenficador do(a) estudante
- Ano da coleta dos dados: 1, 2, 3
- Turma: A, B
- Sexo: Feminino, Masculino
- Mao: Destro, Canhoto, Ambidestro
- TipoSang: A+, A-, …
- ABO: A, B, AB, O
- Rh: +, -
- AtivFisica: nível de atividade física
- Sedentarismo: Não, Sim
- MCT: massa corporal total (kg)
- Estatura: cm
- IMC
'data.frame': 549 obs. of 13 variables:
$ ID : Factor w/ 549 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Ano : num 2 2 1 1 1 3 3 2 3 3 ...
$ Turma : Factor w/ 2 levels "A","B": 2 2 2 2 2 1 1 2 2 2 ...
$ Sexo : Factor w/ 2 levels "F","M": 2 2 2 1 1 2 2 2 2 2 ...
$ Mao : Factor w/ 3 levels "A","C","D": 3 3 3 3 3 3 3 3 3 3 ...
$ TipoSang : Factor w/ 8 levels "A+","O+","B+",..: 1 1 1 1 6 1 1 3 2 2 ...
$ ABO : Factor w/ 4 levels "A","O","B","AB": 1 1 1 1 2 1 1 3 2 2 ...
$ AtivFisica : Factor w/ 5 levels "sempre_inativo",..: 3 2 4 4 2 5 5 4 5 5 ...
$ Sedentarismo: Factor w/ 2 levels "N","S": 1 2 1 1 2 1 1 1 1 1 ...
$ MCT : num 68 82 79 49 52 73 73 60 75 75 ...
$ Estatura : num 173 175 172 160 164 175 175 173 182 182 ...
$ Rh : Factor w/ 2 levels "-","+": 2 2 2 2 1 2 2 2 2 2 ...
$ IMC : num 22.7 26.8 26.7 19.1 19.3 ...
NULL[1] 232[1] 317tabela <- knitr::kable(head(Dados), format="html", caption="Biometria FMUSP")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela| ID | Ano | Turma | Sexo | Mao | TipoSang | ABO | AtivFisica | Sedentarismo | MCT | Estatura | Rh | IMC |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2 | B | M | D | A+ | A | baixa_intensidade | N | 68 | 173 |
|
22.72044 |
| 2 | 2 | B | M | D | A+ | A | atualmente_inativo | S | 82 | 175 |
|
26.77551 |
| 3 | 1 | B | M | D | A+ | A | media_intensidade | N | 79 | 172 |
|
26.70362 |
| 4 | 1 | B | F | D | A+ | A | media_intensidade | N | 49 | 160 |
|
19.14062 |
| 5 | 1 | B | F | D | O- | O | atualmente_inativo | S | 52 | 164 |
|
19.33373 |
| 6 | 3 | A | M | D | A+ | A | alta_intensidade | N | 73 | 175 |
|
23.83673 |
tabela <- knitr::kable(tail(Dados), format="html", caption="Biometria FMUSP")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela| ID | Ano | Turma | Sexo | Mao | TipoSang | ABO | AtivFisica | Sedentarismo | MCT | Estatura | Rh | IMC | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 544 | 544 | 3 | B | F | C | NA | NA | alta_intensidade | N | 54 | 157 | NA | 21.90758 |
| 545 | 545 | 3 | B | F | D | B+ | B | alta_intensidade | N | 47 | 157 |
|
19.06771 |
| 546 | 546 | 3 | B | F | D | B+ | B | sempre_inativo | S | 48 | 162 |
|
18.28989 |
| 547 | 547 | 3 | B | M | A | A+ | A | alta_intensidade | N | 70 | 175 |
|
22.85714 |
| 548 | 548 | 3 | B | M | D | O- | O | atualmente_inativo | S | 67 | 171 |
|
22.91303 |
| 549 | 549 | 3 | B | F | D | A+ | A | media_intensidade | N | 45 | 160 |
|
17.57812 |
Ano Turma Sexo Mao TipoSang ABO
Min. :1.000 A:285 F:232 A : 28 A+ :171 A :220
1st Qu.:1.000 B:264 M:317 C : 37 O+ :132 O :177
Median :2.000 D :442 B+ : 65 B : 99
Mean :2.095 NAs: 42 A- : 49 AB : 42
3rd Qu.:3.000 O- : 45 NAs: 11
Max. :3.000 (Other): 76
NAs : 11
AtivFisica Sedentarismo MCT Estatura
sempre_inativo : 71 N:316 Min. : 41.00 Min. :153.0
atualmente_inativo:152 S:233 1st Qu.: 56.00 1st Qu.:164.0
baixa_intensidade : 59 Median : 64.00 Median :172.0
media_intensidade : 99 Mean : 65.67 Mean :170.9
alta_intensidade :168 3rd Qu.: 73.00 3rd Qu.:177.0
Max. :125.00 Max. :195.0
NAs :7 NAs :7
Rh IMC
- :130 Min. :14.53
+ :408 1st Qu.:20.03
NAs: 11 Median :21.79
Mean :22.32
3rd Qu.:24.07
Max. :40.35
NAs :8 Usando os dados de Biometria_FMUSP.rds, executamos fi_IntervaloPredicao.R:
Referências:
- White, C (1953) Sampling in medical research. BMJ 12: 1284-8.
- Roach, KE (2001) A clinician’s guide to specification and sampling. Journal of Orthopaedic & Sports Physical Therapy 31(12): 753-8.
- Ranstam, J (2009) Sampling uncertainty in medical research. Osteoarthritis and Cartilage 17:1416-9.
Intervalo de confiança de média populacional
Intervalo de confiança da média de uma característica populacional é uma estimativa por intervalo que tem grande probabilidade de conter seu valor médio populacional, usando apenas os dados de uma amostra com observações independentes.
Por exemplo, se uma pesquisa estima que a média de pressão arterial em uma amostra com 30 homens adultos é 120 mmHg, com intervalo de confiança de 95% (IC95%) entre 117 e 123 mmHg, isso quer dizer que há 95% de confiança (probabilidade) de que a média da população está entre 117 e 123 mmHg.
Ou seja, o intervalo de confiança:
- não garante que a média populacional está contida no IC95%, mas
- se repetíssemos o estudo muitas vezes, em 95% dos casos o intervalo conteria o valor da média populacional.
Intervalo de confiança é uma maneira de expressar a incerteza da estimativa, levando em conta a variabilidade dos dados e o tamanho da amostra. Ceteris paribus, quanto maior o tamanho da amostra, menor a incerteza e mais estreito o intervalo.
- Margem de erro: Wikipédia
- Intervalo de confiança: Wikipédia
- Vídeo: Confidence Interval Creation: YouTube
- Confidence Interval Creation: Demo
Referências:
- Banjanovic, ES & Osborne, JW (2016) Confidence intervals for effect sizes: Applying bootstrap resampling. Practical Assessment, Research & Evaluation 21(5).
- Cumming, G & Finch, S (2005) Inference by eye: Confidence intervals and how to read pictures of data. American Psychologist 60(2): 170-80.
- Krzywinski, M & Altman, N (2013) Importance of being uncertain. Nature Methods 10(9): 809-10.
- Puth, M et al. (2015) On the variety of methods for calculating confidence intervals by bootstrapping. Journal of Animal Ecology 84: 892-7.
- Weissgerber TL, Milic NM, Winham SJ, Garovic VD (2015) Beyond bar and line graphs: Time for a new data presentation paradigm. PLoS Biol 13(4): e1002128. doi:10.1371/journal.
- Wonnacott, TH & Wonnacott, RJ (1969) Introductory
Statistics. NJ: Wiley.
- Wonnacott, TH & Wonnacott, RJ (1981) Estatística aplicada à Economia e à Administração. RJ: LTC.
- Wonnacott, TH & Wonnacott, RJ (1990) Introductory Statistics for Business and Economics. NJ: Wiley.
Fórmula do intervalo de confiança da média
\[ \text{IC}^{95\%}(\mu) = \left[\bar{x} \pm t^{0.975}_{n-1}\dfrac{s}{\sqrt{n}}\right] \]
- \(\mu\): média populacional
- \(n\): tamanho da amostra
- \(\bar{x}\): média amostral
- \(s\): desvio-padrão amostral
- \(t^{0.975}_{n-1} =\)
qt(p=.975, df=n-1): quantil da distribuição t - \(n = 2\):
qt(p=.975, df=1)= 12.71 - \(n = 3\):
qt(p=.975, df=2)= 4.3 - \(n = 12\):
qt(p=.975, df=11)= 2.2 - \(n = 20\):
qt(p=.975, df=19)= 2.09 - \(n = 30\):
qt(p=.975, df=29)= 2.05 - \(n = 500\):
qt(p=.975, df=499)= 1.96
Estimação de intervalo de confiança de MCT
- sd: desvio-padrão amostral
- se: erro-padrão (desvio-padrão da média amostral)
| group | n | mean | sd | se |
|---|---|---|---|---|
| F | 230 | 57.6 | 9.1 | 0.6 |
| M | 312 | 71.6 | 12.1 | 0.7 |
Para comparar os intervalos de confiança dos grupos masculino e feminino de MCT, é aconselhável o uso da correção de Bonferroni no nível de significância: \(\alpha_{\text{Bonferroni}}=\alpha/2\).
Executando fi_ToDoMargemErro_Bonf.R, resulta:
---- Feminino:
Media de MCT Feminino = 57.6
n Feminino = 230
Erro-padrao da media de MCT Feminino = 0.6
Margem de erro da media de MCT Feminino = 1.3
IC95Bonf(media populacional de MCT Feminino) = [ 56.3 59 ]
---- Masculino:
Media de MCT Masculino = 71.6
n Masculino = 312
Erro-padrao da media de MCT Masculino = 0.7
Margem de erro da media de MCT Masculino = 1.5
IC95(media populacional de MCT Masculino) = [ 70 73.1 ]Com dados individualizados
alpha <- alfa <- 0.05
k <- nlevels(Dados$Sexo)
alfa.Bonf <- alpha/k # Correção de Bonferroni
# Com dados individualizados
CI_Fmean_t <- confintr::ci_mean(x = Dados.F$MCT,
probs = c(alfa.Bonf/k, 1-alfa.Bonf/k),
type = "t")
print(CI_Fmean_t$interval, digits=4)[1] 56.28 58.98CI_Mmean_t <- confintr::ci_mean(x = Dados.M$MCT,
probs = c(alfa.Bonf/k, 1-alfa.Bonf/k),
type = "t")
print(CI_Mmean_t$interval, digits=4)[1] 70.05 73.13 mean lwr.ci upr.ci
57.63 56.28 58.98 mean lwr.ci upr.ci
71.59 70.05 73.13 gplots::plotmeans(MCT ~ Sexo,
main="Intervalo de confiança 95% Bonferroni: MCT",
barwidth=2,
p=1-alfa.Bonf,
cex=1,
col="black",
barcol="black",
connect=FALSE,
data=Dados)
Sem dados individualizados
tmp <- split(Dados$MCT, Dados$Sexo)
means <- sapply(tmp, mean, na.rm=TRUE)
stdev <- sqrt(sapply(tmp, var, na.rm=TRUE))
n <- sapply(tmp, length)
ciw <- qt(1-alfa.Bonf/k, n-1) * stdev / sqrt(n)
g <- length(unique(Dados$Sexo))
gplots::plotCI(x=means,
uiw=ciw,
main="Intervalo de confiança Bonferroni 95%: MCT",
ylab="MCT",
xlab="Sexo",
p=1-alfa.Bonf,
cex=1,
col="black",
barcol="black",
labels=round(means,-3),
xaxt="n")
axis(side=1, at=1:g, labels=names(tmp), cex=0.7)
Intervalo de confiança da média: condições de validade
O intervalo de confiança tem validade se:
- Delineamento entre participantes: observações independentes
- Variável intervalar
- Distribuição normal da variável intervalar
A seguir, são realizadas as análises descritivas e o teste de normalidade de \(MCT\) feminino e masculino:
item group1 vars n mean sd median trimmed mad min max range skew
X1*1 1 F 1 230 15.73 6.81 15 15.46 7.41 1 34 33 0.32
X1*2 2 M 1 312 28.18 10.82 29 28.21 10.38 1 53 52 -0.06
kurtosis se
X1*1 -0.51 0.45
X1*2 -0.28 0.61
m.F <- mean(Dados.F$MCT, na.rm=TRUE)
s.F <- sd(Dados.F$MCT, na.rm=TRUE)
plot(density(Dados.F$MCT, na.rm=TRUE),
ylim=c(0, 0.06),
main="Feminino",
xlab="MCT (kg)",
ylab="Densidade")
rug(jitter(Dados.F$MCT))
x.F <- seq(m.F-4*s.F, m.F+4*s.F, length.out=1e3)
y.F <- dnorm(x.F, m.F, s.F)
lines(x.F, y.F,lty=2)
legend("topright", legend=c("estimada", "normal"),
lty=c(1,2), bty="n")
Results of Hypothesis Test
--------------------------
Alternative Hypothesis:
Test Name: Shapiro-Wilk normality test
Data: Dados.F$MCT
Test Statistic: W = 0.832
P-value: 4.94e-15m.M <- mean(Dados.M$MCT, na.rm=TRUE)
s.M <- sd(Dados.M$MCT, na.rm=TRUE)
plot(density(Dados.M$MCT, na.rm=TRUE),
ylim=c(0, 0.04),
main="Masculino",
xlab="MCT (kg)",
ylab="Densidade")
rug(jitter(Dados.M$MCT))
x.M <- seq(m.M-4*s.M, m.M+4*s.M, length.out=1e3)
y.M <- dnorm(x.M, m.M, s.M)
lines(x.M, y.M,lty=2)
legend("topright", legend=c("estimada", "normal"),
lty=c(1,2), bty="n")
Results of Hypothesis Test
--------------------------
Alternative Hypothesis:
Test Name: Shapiro-Wilk normality test
Data: Dados.M$MCT
Test Statistic: W = 0.954
P-value: 2.61e-08Região de confiança da média e desvio-padrão populacionais
A região de confiança para média e desvio-padrão populacionais é o análogo bidimensional do intervalo de confiança: em vez de dar apenas um intervalo para a média ou para o desvio-padrão separadamente, ela mostra todas as combinações plausíveis (μ, σ) que são consistentes com os dados a um certo nível de confiança (e.g.: 95%), assumindo distribuição normal.
Geometricamente, a curva no gráfico é a fronteira da região de confiança: qualquer ponto dentro representa um par (média e desvio-padrão populacionais) compatível com os dados, segundo a estatística de máxima verossimilhança para a normal.
## plot the 95% confidence region for normal parameters
MCT.M <- Dados.M$MCT
print(DescTools::MeanCI(MCT.M, na.rm=TRUE), digits=4) mean lwr.ci upr.ci
71.59 70.24 72.94 var lwr.ci upr.ci
12.09 11.21 13.12 invisible(conf::crplot(dataset=sort(MCT.M),
distn="norm",
alpha=alfa,
pts=FALSE,
animate=FALSE,
info=TRUE,
silent=TRUE,
# animate=TRUE,
main="MCT Masculino"))
mean lwr.ci upr.ci
57.63 56.45 58.81 var lwr.ci upr.ci
9.053 8.294 9.965 invisible(conf::crplot(dataset=sort(MCT.F),
distn="norm",
alpha=alfa,
pts=FALSE,
animate=FALSE,
info=TRUE,
silent=TRUE,
# animate=TRUE,
main="MCT Feminino"))
invisible(outM <- conf::crplot(sort(Dados.M$MCT),
distn="norm",
alpha=alfa.Bonf,
pts=FALSE,
animate=FALSE,
info=TRUE,
silent=TRUE,
origin=TRUE))
invisible(outF <- conf::crplot(sort(Dados.F$MCT),
distn="norm",
alpha=alfa.Bonf,
pts=FALSE,
animate=FALSE,
info=TRUE,
silent=TRUE,
origin=TRUE))
plot(outM$mu,
outM$sigma,
type="l",
col="blue",
xlab=expression(mu), ylab=expression(sigma),
xlim=c(55,75),
ylim=c(7,14),
main="Regiões de confiança 95% Bonf para μ e σ de MCT")
lines(outF$mu,
outF$sigma,
col="red")
points(outM$muhat, outM$sigmahat, pch=19, col="blue", cex=1.2)
points(outF$muhat, outF$sigmahat, pch=19, col="red", cex=1.2)
legend("topleft",
legend=c("Masculino", "Feminino"),
col=c("blue", "red"),
lty=1,
bty="n")
Suposição de normalidade: condição suficiente
Uma das suposições dos seguintes intervalos de confiança da média populacional (por fórmula) é a distribuição normal de MCT. A hipótese nula de distribuição normal de MCT foi rejeitada com \(\alpha=0.05\) para estudantes das condições Masculino e Feminino.
A suposição de normalidade é uma condição suficiente para a validade do intervalo de confiança (por fórmula). Ela não é uma condição necessária para sua validade.
O intervalo de confiança da média populacional por fórmula é válido, mesmo com distribuição de MCT não normal, se o tamanho da amostra por condição de sexo é maior que 30 (Teorema Central do Limite) ou usando o método de estimação chamado de bootstrapping. TCL e bootstrapping não necessitam da suposição de normalidade para a validade do intervalo de confiança. O método bootstrapping produz intervalo de confiança válida para tamanho de amostra a partir de duas observações independentes (\(n\ge 2\)).
# Com dados individualizados
R <- 1e4
seed <- 123
CI_Fmean_t <- confintr::ci_mean(x = Dados.F$MCT,
probs = c(alfa.Bonf/k, 1-alfa.Bonf/k),
type = "t")
print(CI_Fmean_t$interval, digits=4)[1] 56.28 58.98CI_Fmean_B <- confintr::ci_mean(x = Dados.F$MCT,
probs = c(alfa.Bonf/k, 1-alfa.Bonf/k),
type = "bootstrap",
R = R,
seed = seed)
print(CI_Fmean_B$interval, digits=4)[1] 56.43 59.20CI_Mmean_t <- confintr::ci_mean(x = Dados.M$MCT,
probs = c(alfa.Bonf/k, 1-alfa.Bonf/k),
type = "t")
print(CI_Mmean_t$interval, digits=4)[1] 70.05 73.13CI_Mmean_B <- confintr::ci_mean(x = Dados.M$MCT,
probs = c(alfa.Bonf/k, 1-alfa.Bonf/k),
type = "bootstrap",
R = R,
seed = seed)
print(CI_Mmean_B$interval, digits=4)[1] 70.10 73.22print(DescTools::MeanCI(x=Dados.F$MCT,
conf.level=1-alfa.Bonf,
method="classic",
na.rm=TRUE), digits=4) mean lwr.ci upr.ci
57.63 56.28 58.98 set.seed(123)
print(DescTools::MeanCI(x=Dados.F$MCT,
conf.level=1-alfa.Bonf,
method="boot",
R=1e4,
sides="two.sided",
na.rm=TRUE), digits=4) mean lwr.ci upr.ci
57.63 56.25 58.90 print(DescTools::MeanCI(x=Dados.M$MCT,
conf.level=1-alfa.Bonf,
method="classic",
na.rm=TRUE), digits=4) mean lwr.ci upr.ci
71.59 70.05 73.13 set.seed(123)
print(DescTools::MeanCI(x=Dados.M$MCT,
conf.level=1-alfa.Bonf,
method="boot",
R=1e4,
sides="two.sided",
na.rm=TRUE), digits=4) mean lwr.ci upr.ci
71.59 70.04 73.12 Intervalo de predição vs. Intervalo de confiança
Pacote confintr
O pacote confintr calcula intervalos de confiança
para:
ci_chisq_ncp: parâmetro de não centralidade da distribuição qui-quadradoci_cor: coeficiente de correlaçãoci_cramersv: V de Cramerci_f_ncp: parâmetro de não centralidade da distribuição Fci_IQR: intervalo interquartilci_kurtosis: curtoseci_mad: desvio absoluto mediano (MAD)ci_mean: médiaci_mean_diff: diferença de médias populacionaisci_median: medianaci_median_diff: diferença de medianas populacionais entre duas amostrasci_oddsratio: odds ratioci_proportion: proporçãoci_quantile: quantilci_quantile_diff— diferença de quantis populacionais entre duas amostrasci_rsquared: R²ci_sd: desvio-padrãoci_skewness: assimetriaci_var: variância
Pacote DescTools
MeanCI: Confidence Intervals for the MeanMeanCIn:Sample Size for a Given Width of a Confidence Interval for a MeanMeanDiffCI: Confidence Interval For Difference of MeansMedianCI: Confidence Interval for the MedianVarCI: Confidence Intervals for the VarianceBootCI: Simple Bootstrap Confidence Interval
set.seed(123)
print(DescTools::BootCI(Dados.F$MCT,
FUN=mean,
conf.level=1-alpha,
sides="two.sided",
na.rm=TRUE,
R=R), digits=4) mean lwr.ci upr.ci
57.63 56.45 58.80 set.seed(123)
print(DescTools::BootCI(Dados.F$MCT,
FUN=DescTools::Var,
conf.level=1-alpha,
sides="two.sided",
na.rm=TRUE,
R=R), digits=4)DescTools::Var lwr.ci upr.ci
81.95 42.13 122.56 spearman <- function(x,y) cor(x, y, method="spearman", use="p")
set.seed(123)
print(DescTools::BootCI(Dados.F$MCT,
Dados.F$Estatura,
FUN=spearman,
conf.level=1-alpha,
sides="two.sided",
R=R), digits=4)spearman lwr.ci upr.ci
0.4918 0.3866 0.6003 Teorema Central do Limite
O Teorema Central do Limite (TCL) assume:
- observações independentes,
- variável intervalar,
- amostra com pelo menos 30 observações.
Conforme Dancey & Reidy, 2019, p. 108:
“[…] o Teorema Central do Limite […] declara que à medida que o tamanho das amostras que selecionamos aumenta, mais próxima da população [sic, da média populacional] será a média destas médias amostrais e mais próxima de uma normal estará a distribuição das médias amostrais.”
Segundo Wikipedia:
“Seja uma amostra aleatória simples \(X_1, X_2, \ldots ,X_n\) de tamanho \(n\) dada a partir de uma população com média \(\mu\) e variância \(\sigma ^{2}\) finitas. À medida que \(n\) cresce, a distribuição amostral da média \(\dfrac{\sum_{i=1}^{n}{X_i}}{n} = \bar{X}\) aproxima-se de uma distribuição normal com média \(\mu\) e variância \(\dfrac{\sigma^2}{n}\).”
Em nossas palavras:
Independentemente do formato da distribuição da variável na população, a distribuição das médias amostrais tem distribuição assintoticamente normal com média igual à média populacional e erro-padrão da média igual ao desvio-padrão populacional dividido pela raiz quadrada do tamanho da amostra.
Erro-padrão
O erro-padrão é fundamental na definição do TCL.
Conforme Dancey & Reidy, 2019:
“… erro-padrão [é igual ao] desvio-padrão populacional dividido pela raiz quadrada do tamanho da amostra.”
Este erro-padrão é estimado por:
\[\text{EP}=\dfrac{s}{\sqrt{n}}\]
sendo que \(s\) é o desvio-padrão amostral e \(n\) é o tamanho da amostra.
Note que o desvio-padrão populacional \(\sigma\) não aparece neste estimador porque é desconhecido.
Estas definições podem não ser muito esclarecedoras a quem nunca as encontrou. São muitas afirmações entrelaçadas:
Qual a relação entre a média populacional \(\mu\) e a média amostral \(\bar{X}\)?
Qual a relação entre o desvio-padrão populacional \(\sigma\) e o desvio-padrão amostral \(s\)?
Qual é a origem da expressão \({s}\big/{\sqrt{n}}\) e qual é seu significado?
A quais médias amostrais (no plural) ou a qual distribuição amostral da média o enunciado se refere, uma vez que há apenas uma amostra?
Qual deve ser o formato da distribuição da variável na população?
De qual variável é a distribuição normal mencionada neste enunciado?
Wonnacott & Wonnacott (1969, 1990) descrevem o TCL textualmente e graficamente:
Esta figura pode ser muito esclarecedora, mas somente depois que alguém é apresentado ao processo do TCL. De início parece aumentar o enigma. Observe a figura. Há duas distribuições desenhadas: população à esquerda e amostra à direita. A distribuição da população é bimodal, aparecendo as letras gregas \(\mu\) (denotando média) em um eixo \(x\) e \(\sigma\) (denotando desvio-padrão populacional). Há vários \(\bar{x}\) em cinza e um em preto (amostra observada). Dois (um cinza, um preto) aparecem entre a população e a amostra. Vários em cinza e um em preto estão agrupados abaixo da distribuição amostral. A distribuição amostral, denotada como gaussiana (sinônimo de distribuição normal) tem novamente a letra \(\mu\) (que é a média populacional, apesar de ser a figura da amostra) e o eixo da distribuição amostral é \(\bar{x}\). A expressão \({\sigma}\big/{\sqrt{n}}\) reaparece no lugar de um desvio-padrão. Voltaremos a esta figura adiante.
Simulação
A proposta, aqui, é utilizar uma simulação pode mostrar o que o erro-padrão representa e como o TCL funciona.
Suponha que a distribuição do nível de hemoglobina, na população geral, tenha distribuição normal com média de 11.9 g/dl e desvio-padrão de 1.5 g/dl. Estes valores foram extraídos de Gilbertson et al. (2008).
A população, caso fosse acessível, teria a seguinte distribuição (fi_Normal.R):
Populacao:
media populacional: 11.9
d.p. populacional: 1.5O fato de sabermos tratar-se de uma distribuição normal, completamente definida por sua média e desvio-padrão, traz a conveniência de aproveitarmos suas propriedades para saber, por exemplo, que aproximadamente 95% dos indivíduos desta população têm hemoglobina entre 8.96 e 14.84 g/dl (perto do valor da média mais ou menos 2 desvios-padrão, que correspondem a 8.9 e 14.9 g/dl).
Como é possível conhecer o valor da média populacional se não forem mensurados todos os indivíduos da população? Seria possível estimar esta média a partir de uma amostra?
Geramos uma amostra simulada do nível de hmoglobina de uma população
normal de 100 indivíduos (fi_Normal_e_Amostra.R):
Populacao:
media populacional: 11.9
d.p. populacional: 1.5
Amostra:
media amostral: 11.945
d.p. amostral: 1.6092|
|
A distribuição normal que representa a população está em coloração mais apagada porque não temos acesso a ela. Em linha mais fina e de cor mais intensa temos a distribuição da amostra de 100 indivíduos. Observe que a distribuição amostral não é perfeitamente normal, embora siga o perfil geral da população.
Na prática, porém, é tudo o que teríamos é uma amostra, sem qualquer
referência da população (fi_Amostra.R):
Amostra:
media amostral: 11.945
d.p. amostral: 1.6092Analiticamente:
\[\text{IC}^{95\%}(\mu) = \left[\bar{x} \pm t^{0.975}_{n-1}\dfrac{s}{\sqrt{n}}\right]\]
Computa-se (fi_IntervaloConfiancaAnalitico.R):
n = 100
qt(0.975,n-1) = 1.984
EP = s/sqrt(n) = 0.161
IC95 = 11.945 +- 1.984(0.975,99) * 1.609/sqrt(100)
IC95 = [11.626,12.264]|
Para recordar (nunca saberíamos), a média populacional é 11.9. A partir da amostra localizamos um intervalo de confiança 95% que contém a média populacional. |
Bootstrapping
A maior parte dos textos que explica o TCL parte de um experimento imaginário, no qual infinitas amostras (ou um número muito grande quando pretendemos fazer simulações) de igual tamanho \(n\) são retiradas de uma população hipotética.
Na prática este experimento não faz sentido porque, se pudéssemos retirar muitas amostras poderíamos retirar uma amostra maior e mais representativa da população.
Existe, porém, uma alternativa possível e que produz o mesmo resultado para que, com base em um única amostra, possamos chegar a conclusões que possam ser extrapoladas para a população de onde ela proveio: bootstrapping.
Reamostragem por bootstrapping, com base no comportamento do teorema central do limite, serve para calcular erro-padrão e intervalo de confiança.
Faremos várias reamostragens com reposição de 100 elementos sobre a única amostra de 100 elementos que temos para trabalhar.
O bootstrapping, sobre a amostra obtida da população fictícia, resulta em:
implementado com fi_boot01.R
Amostra unica (n = 100):
media amostral: 11.945
d.p. amostral: 1.6092
Amostras:10000 com n = 100
media das medias amostrais: 11.9435
media dos desvios padrao das medias amostrais: 1.5956|
|
O gráfico mostra linhas representando as 10000 reamostras de 100 indivíduos (o gráfico só exibe 500 linhas das 10000 geradas para melhor visibilidade) às quais a distribuição da amostra original foi sobreposta. Como acontecia entre a amostra e a população, também as reamostras seguem o perfil geral da amostra única que tomamos como ponto de partida. Na parte alta do gráfico vemos a média e os desvios-padrão da amostra original, e a média das médias e a média dos desvios-padrão das 10000 reamostras.
Estas últimas médias (das médias e dos desvios-padrão das 10000
reamostras) foram obtidas do bootstrapping e seus valores podem
ser encontrados: cada uma das 10000 reamostras tem sua média e
desvio-padrão, respectivamente guardadas nos vetores
boot.medias e boot.sds (contendo 10000
elementos cada um dos vetores).
Observe os próximos códigos em R para entender quais vetores e quais valores foram calculados:
boot.medias (10000 elementos): 11.82 11.95 11.96 12.17 12.26 11.76 11.92 11.89 11.66 11.71 12.12 11.88 11.98 11.78 11.79 12.15 11.87 12.06 12.09 12.07 12.13 11.94 12.13 12.03 12.08 12 11.94 11.82 11.97 12.07 12.04 11.86 12.07 11.81 11.9 12.01 12.24 12.04 11.95 12.13 12.23 12.06 12.06 12.11 11.89 11.77 11.85 12.09 12.07 12.01 ... 11.98 11.87 11.7 12 11.89 11.79 11.95 12.15 11.9 11.78 11.96 12.17 11.95 12.05 11.91 12.12 12.14 11.87 11.85 11.95 12.18 11.78 11.72 12.01 12.14 12.21 11.99 12.03 12.04 11.97 11.78 12.11 11.93 12.12 11.92 12.06 12.02 11.74 11.8 11.97 12.04 11.84 11.88 11.7 12 12.08 11.99 11.74 11.94 11.89 11.78
boot.sds (10000 elementos): 1.66 1.68 1.66 1.59 1.85 1.62 1.69 1.66 1.59 1.4 1.68 1.46 1.62 1.61 1.57 1.68 1.57 1.72 1.49 1.64 1.87 1.63 1.63 1.73 1.44 1.53 1.6 1.68 1.75 1.29 1.56 1.47 1.71 1.77 1.62 1.57 1.56 1.59 1.49 1.8 1.62 1.75 1.52 1.68 1.48 1.49 1.68 1.66 1.59 1.87 ... 1.53 1.37 1.57 1.65 1.53 1.54 1.72 1.47 1.48 1.53 1.5 1.65 1.7 1.73 1.42 1.43 1.57 1.35 1.57 1.62 1.67 1.63 1.67 1.69 1.47 1.53 1.68 1.86 1.61 1.64 1.86 1.5 1.64 1.5 1.62 1.68 1.51 1.67 1.5 1.56 1.52 1.64 1.9 1.79 1.51 1.7 1.65 1.48 1.71 1.75 1.74implementado com fi_exibevetores.R
Amostras:10000 com n = 100
media das medias amostrais: 11.9435
media dos desvios padrao das medias amostrais: 1.5956|
|
Adicionalmente, considerando apenas boot.medias, podemos
calcular mais um desvio-padrão:
desvio padrao das medias amostrais: 0.1593|
|
O erro-padrão da média, estimado aqui pelo desvio-padrão das médias amostrais obtidas por muitas reamostragens, é denotado como \(\text{EP}_{\mu}\).
A distribuição de boot.medias pode ser expressa
graficamente, mantendo-se a mesma escala do eixo \(x\) que utilizamos para a distribuição dos
valores amostrais:
A distribuição destas médias amostrais é mais alta e mais estreita (a área de uma distribuição de probabilidades é mantida em 1) porque o desvio-padrão das médias amostrais é o \(\text{EP}_{\mu}\) estimado. Observe que este erro-padrão é cerca de \(1 \big/ 10\) da média dos desvios-padrão das amostras (\(s\), que é estimador do desvio-padrão populacional, \(\sigma\)). Não é uma coincidência \(n=100\) e \(\sqrt{n}=10\).
Para melhor visualização, exibimos a mesma curva da distribuição das médias amostrais ajustando a escala do eixo \(x\), adicionando a referência amostral e populacional (esta última apenas para conferência; lembre-se que jamais a teremos na prática). Temos:
implementado com fi_boot01c.R
Amostra unica (n = 100):
media amostral: 11.945
d.p. amostral: 1.6092
Amostras:10000 com n = 100
media das medias amostrais: 11.9435
media dos d.p. amostrais: 1.5956
d.p. das medias amostrais: 0.1593 (erro padrao da media)
IC95(mu): [11.636,12.26]O círculo sólido (na mesma coordenada \(x\) que o círculo preto) e as faixas que aparecem na mesma cor da distribuição na parte alta do gráfico correspondem, respectivamente, à média das médias amostrais, \(\bar{\bar{x}}\), e aos erros padrão destas médias (de \(\pm1\) a \(\pm3\) \(\text{EP}_{\mu}\)). Na parte de baixo do gráfico adicionamos uma linha horizontal com halteres em preto que representa o intervalo de confiança 95%: seus limites indicam 95% da área sob a curva da distribuição das médias das reamostras com um círculo sólido preto na posição de \(\bar{\bar{x}}\) que corresponde à estimativa pontual da média populacional \(\mu\); observe que os limites dos halteres estão aproximadamente a \(\mu \pm2 \text{EP}_{\mu}\). Adicionamos, também na parte de baixo do gráfico, um círculo sólido na posição da média da amostra original e um círculo branco na posição da média populacional para referência.
Colocando em outras palavras, uma vez que confiamos na única amostra que obtivemos e supomos que esta amostra representa bem a população de onde foi retirada, utilizamos o bootstrapping para criar um número grande de variantes desta amostra (neste exemplo, 10000 reamostras). Cada uma destas reamostras tem sua média, \(\bar{x}\). A distribuição obtida mostra que é mais provável que a maioria dos valores de \(\bar{x}\) sejam próximas da média de nossa amostra única e que valores de \(\bar{x}\) mais afastados são progressivamente mais improváveis, de tal forma que a distribuição dos 10000 valores de \(\bar{x}\) é aproximadamente normal. Encontramos, então, a faixa de valores que contém 95% destes \(\bar{x}\). Esta faixa de valores é o intervalo em que, com confiança de 95%, esperamos encontrar a média populacional. Não conhecendo a média populacional (como acontece em uma amostra obtida na prática) poderíamos dizer que ela está localizada entre 11.636 e 12.26 com 95% de confiança. Esta média populacional está representada no gráfico como um círculo branco apenas para mostrar o sucesso deste procedimento com bootstrapping, pois o intervalo de confiança 95% estimado a englobou.
O bootstrapping serve para situações com amostras menores.
Vamos repetir este procedimento com uma amostra de \(n=9\) da mesma população com média de Hb de
11.6 g/dl e desvio-padrão de 1.5 g/dl. Com 100000 reamostragens,
obteríamos o seguinte (código disponível em fi_boot02.R):
Amostra unica (n = 9):
media amostral: 12.4556
d.p. amostral: 1.9501
Amostras:1e+05 com n = 9
media das medias amostrais: 12.4589
media dos d.p. amostrais: 1.8137
d.p. das medias amostrais: 0.6119 (erro padrao da media)
IC95: [11.2556,13.6444]Observe o que aconteceu com o intervalo de confiança 95% [11.2556, 13.6444]. Com uma amostra menor, a distribuição normal das médias das reamostras é mais larga (\(\text{EP}_{\mu}\) é aproximadamente \(1 \big/ \sqrt{9}\) do desvio-padrão amostral, que por sua vez é estimador do desvio-padrão populacional) e, consequentemente, a faixa de valores que afirmamos conter a média populacional com confiança de 95% é mais larga (com \(n=100\) o IC95 era [11.636, 12.26]).
|
O TCL não vale apenas para estimar a média populacional. Podemos
aproveitar o mesmo bootstrapping para estimar o desvio-padrão
populacional. A distribuição dos desvios-padrão das reamostras e os
cálculos para o intervalo de confiança são (código disponível em
Aqui estimamos o \(\text{EP}_{\sigma}\) e obtivemos o intervalo de confiança 95% do \(\sigma\). Encontramos que o desvio-padrão populacional deve estar localizado entre 1.1616 e 2.342. Como estamos em uma situação simulada, sabemos que a população de onde retiramos a amostra tinha desvio-padrão igual a 1.5 g/dl, valor este que está dentro do intervalo e a estimativa por bootstrapping foi bem sucedida. |
Bootstrapping (variável populacional não normal)
No exemplo acima tiramos amostras (com \(n=100\) e \(n=9\)) de uma população cuja variável Hb tinha distribuição normal. A distribuição das médias reamostrais obtida por bootstrapping tinha, igualmente, distribuição normal.
Para verificar o comportamento do bootstrapping em situação menos favorável, simularemos com uma população com variável que não tem distribuição normal.
Em primeiro lugar ativamos algumas funções de apoio que desenvolvemos e, então, criamos um vetor com uma população fictícia. Vamos imaginar que pretendemos medir o valor do colesterol LDL e, embora não saibamos, esta população tem indivíduos misturados: hipo, normo e hipercolesterolêmicos.
Observe o código em R (implementado em fi_criapopsample.R). Criaremos a
população artificialmente e, em seguida, retiramos uma única amostra
desta população com tamanho \(n=36\):
obtendo-se:
Populacao:
media populacional: 98.0423
d.p. populacional: 41.5684
media amostral: 98.5392
d.p. amostral: 38.8217Nesta situação, a amostra captura informação (com imperfeição) da população. Assim como é a população, esta amostra não tem distribuição normal. Repare que os “picos” e “vales” existentes na distribuição populacional foram incompletamente representados, mas a média e desvio-padrão amostrais são similares àqueles da população.
Lembramos novamente: na prática temos tão somente esta amostra, sem a referência populacional:
media amostral: 98.5392
d.p. amostral: 38.8217Com o bootstrapping (código em fi_bootstrapping.R) obtemos:
Amostra unica (n = 36):
media amostral: 98.5392
d.p. amostral: 38.8217
Amostras:10000 com n = 36
media das medias amostrais: 98.606
media dos d.p. amostrais: 38.0609
erro padrao da media: 6.5011
IC95(media): [86.1206,111.6244]Surpreendentemente, a distribuição das médias das reamostras continua sendo aproximadamente normal (apesar da amostra não o ser) com média igual à da amostra original e desvio-padrão dado por \(\text{EP}=s/\sqrt{n}\) (ou próximo a estes valores porque aqui estamos simulando).
Exibindo novamente o último gráfico, adicionando o intervalo de
confiança e uma distribuição normal (linha pontilhada) para comparação,
podemos observar melhor a distribuição das médias amostrais (código em
fi_bootstrapping2.R):
Neste caso temos o desvio-padrão da amostra as=38.8217 e
o desvio-padrão das 10000 médias das reamostras epm=6.5011.
Podemos conferir que as/epm=5.9715587 é valor próximo a
sqrt(36)=6. A concordância talvez fosse melhor se usássemos
um bootstrapping com mais reamostras; aqui exemplificamos com
10000 reamostras, mas é recomendado que se aplique entre 100000 (=1e+05)
e 1000000 (=1e+06) reamostras para fins de decisão estatística.
|
Interessantemente, o TCL sob bootstrapping também funciona
para o desvio-padrão neste caso de distribuição não normal da variável
populacional (código em
O intervalo de confiança 95% estimado incluiu, como pode-se verificar, o desvio-padrão populacional. |
O bootstrapping é um procedimento sempre possível com o uso de computadores, capturando o experimento imaginário de obter inúmeras amostras de uma mesma população e reconstituindo propriedades da população inalcançável de onde a amostra proveio.
|
Voltando à figura apresentada no início de texto, agora é possível entender seus elementos:
|
DescTools::MeanCI: estimação de intervalo de
confiança
alfa <- 0.05
alfa.Bonf <- alfa/length(unique(Dados$Sexo))
# IC95% por DescTools::MeanCI
print(DescTools::MeanCI(Dados.F$MCT, conf.level=1-alfa.Bonf, na.rm=TRUE),digits=4) mean lwr.ci upr.ci
57.63 56.28 58.98 mean lwr.ci upr.ci
71.59 70.05 73.13 # Os mesmos IC95% por t.test
print(t.test(Dados.F$MCT, conf.level=1-alfa.Bonf)$conf.int[1:2],digits=4)[1] 56.28 58.98[1] 70.05 73.13# IC95% por DescTools::MeanCI com desvio-padrao conhecido
print(DescTools::MeanCI(Dados.F$MCT, sd=10, conf.level=1-alfa.Bonf, na.rm=TRUE),digits=4) mean lwr.ci upr.ci
57.63 56.15 59.11 mean lwr.ci upr.ci
71.59 69.81 73.37 # IC95% unilateral por DescTools::MeanCI
# round(DescTools::MeanCI(Dados.F$MCT, sides="left", na.rm=TRUE),1)
# round(DescTools::MeanCI(Dados.M$MCT, sides="left", na.rm=TRUE),1)
# round(DescTools::MeanCI(Dados.F$MCT, sides="right", na.rm=TRUE),1)
# round(DescTools::MeanCI(Dados.M$MCT, sides="right", na.rm=TRUE),1)
# IC95% estratificado por Sexo
print(tapply(Dados$MCT,
Dados$Sexo,
FUN=DescTools::MeanCI,
conf.level=1-alfa.Bonf,
na.rm=TRUE), digits=4)$F
mean lwr.ci upr.ci
57.63 56.28 58.98
$M
mean lwr.ci upr.ci
71.59 70.05 73.13 # IC95% por bootstrapping para sexo feminino
print(DescTools::MeanCI(Dados.F$MCT,
method="boot",
type="perc",
R=1e4,
sides="two.sided",
conf.level=1-alfa.Bonf,
na.rm=TRUE), digits=4) mean lwr.ci upr.ci
57.63 56.32 59.04 # IC95% de MCT, Estatura e IMC para sexo feminino
print(do.call("rbind", lapply(Dados.F[, c("MCT","Estatura","IMC")],
FUN=DescTools::MeanCI,
conf.level=1-alfa.Bonf,
na.rm=TRUE)), digits=4) mean lwr.ci upr.ci
MCT 57.63 56.28 58.98
Estatura 164.34 163.39 165.29
IMC 21.24 20.82 21.65print(t(sapply(Dados.F[,c("MCT","Estatura","IMC")],
FUN=DescTools::MeanCI,
conf.level=1-alfa.Bonf,
na.rm=TRUE)), digits=4) mean lwr.ci upr.ci
MCT 57.63 56.28 58.98
Estatura 164.34 163.39 165.29
IMC 21.24 20.82 21.65Gráfico de intervalo de confiança de 95% de média: delineamento entre participantes
Obtido com fi_IntervalosDadosEntre.R:
ID Ano Turma Sexo Mao TipoSang ABO AtivFisica Sedentarismo MCT
1 2 B M D A+ A baixa_intensidade N 68
2 2 B M D A+ A atualmente_inativo S 82
3 1 B M D A+ A media_intensidade N 79
4 1 B F D A+ A media_intensidade N 49
5 1 B F D O- O atualmente_inativo S 52
6 3 A M D A+ A alta_intensidade N 73
Estatura Rh IMC
173 + 22.72
175 + 26.78
172 + 26.70
160 + 19.14
164 - 19.33
175 + 23.84
Gráfico de intervalo de confiança de 95% de média: delineamento intraparticipantes
Obtido com fi_IntervalosDadosIntra.R:
subject condition value
1 pretest 59.4
2 pretest 46.4
3 pretest 46.0
4 pretest 49.0
5 pretest 32.5
6 pretest 45.2
7 pretest 60.3
8 pretest 54.3
9 pretest 45.4
10 pretest 38.9
1 posttest 64.5
2 posttest 52.4
3 posttest 49.7
4 posttest 48.7
5 posttest 37.4
6 posttest 49.5
7 posttest 59.9
8 posttest 54.1
9 posttest 49.6
10 posttest 48.5 condition N value valueNormed sd se ci
posttest 10 51.43 51.43 2.262 0.7154 1.921
pretest 10 47.74 47.74 2.262 0.7154 1.921
condition N value sd se ci
pretest 10 47.74 8.599 2.719 7.301
posttest 10 51.43 7.254 2.294 6.159
Two within-subjects variables Subject RoundMono SquareMono RoundColor SquareColor
1 1 41 40 41 37
2 2 57 56 56 53
3 3 52 53 53 50
4 4 49 47 47 47
5 5 47 48 48 47
6 6 37 34 35 36
7 7 47 50 47 46
8 8 41 40 38 40
9 9 48 47 49 45
10 10 37 35 36 35
11 11 32 31 31 33
12 12 47 42 42 42 Subject Time Shape ColorScheme
1 41 Round Monochromatic
2 57 Round Monochromatic
3 52 Round Monochromatic
4 49 Round Monochromatic
5 47 Round Monochromatic
6 37 Round Monochromatic Shape ColorScheme N Time TimeNormed sd se ci
Round Colored 12 43.58 43.58 1.212 0.3500 1.043
Round Monochromatic 12 44.58 44.58 1.331 0.3844 1.146
Square Colored 12 42.58 42.58 1.462 0.4219 1.258
Square Monochromatic 12 43.58 43.58 1.261 0.3641 1.085
HH::CIplot: simulação de intervalo de confiança
Simulação de intervalo de confiança em Rpsychologist
DescTools::MedianCI: IC de mediana
alfa.Bonf <- alfa/nlevels(Dados$Sexo) # Correção de Bonferroni
print(DescTools::MedianCI(Dados.F$MCT, conf.level=1-alfa.Bonf, na.rm=TRUE), digits=4)median lwr.ci upr.ci
56 55 58
attr(,"conf.level")
[1] 0.9792median lwr.ci upr.ci
72 70 73
attr(,"conf.level")
[1] 0.9799DescTools::VarCI: IC de desvio-padrão
alfa.Bonf <- alfa/nlevels(Dados$Sexo) # Correção de Bonferroni
print(sqrt(DescTools::VarCI(Dados.F$MCT, conf.level=1-alfa.Bonf, na.rm=TRUE)), digits=4) var lwr.ci upr.ci
9.053 8.192 10.107 var lwr.ci upr.ci
12.09 11.09 13.28 Intervalo de confiança pivotal
Uma forma mais robusta é o intervalo de confiança pivotal (Chihara & Hesterberg, 2019).
O método pivotal com teste t consiste em construir um intervalo de confiança para a média utilizando uma estatística cuja distribuição não depende de parâmetros desconhecidos.
Para dados \(x_1,\dots,x_n\) com média amostral \(\bar{x}\) e desvio-padrão amostral \(x\), define-se a estatística pivotal:
\[ t = \dfrac{\bar{x} - \mu}{s/\sqrt{n}} \]
que, sob normalidade, segue distribuição t de Student com \(n-1\) graus de liberdade.
No procedimento bootstrap pivotal:
- Calcula-se \(\bar{x}\) e \(s\) da amostra observada.
- Gera-se um grande número de reamostragens com reposição.
- Para cada reamostragem \(b\), calcula-se
\[ t_b = \dfrac{\bar{x}_b - \bar{x}}{s_b / \sqrt{n}} \]
sendo que \(\bar{x}_b\) e \(s_b\) são a média e o desvio-padrão da reamostragem.
Obtêm-se os quantis empíricos \(q^{\alpha/2}_b\) e \(q^{1-\alpha/2}_b\) da distribuição de \(t_b\).
O intervalo de confiança \(1-\alpha\) para \(\mu\) é dado por
\[ \text{IC}^{(1-\alpha)\%}_b(\mu)=\left[\bar{x} - q^{\alpha/2}_b \dfrac{s}{\sqrt{n}}, \;\bar{x} + q^{1-\alpha/2}_b \dfrac{s}{\sqrt{n}}\right] \]
Diferencia-se do IC clássico porque não se usa a distribuição t teórica, mas sim a distribuição empírica de \(t_b\) estimada a partir dos dados, tornando o método mais robusto a desvios da normalidade e amostras pequenas.
Voltamos ao exemplo do colesterol LDL e simulamos com simulaIC95pivotal.R utilizando a
população simulada anteriormente. A execução pode ser feita no ambiente
do RStudio com o botão [Source] ou com o seguinte comando na
Console:
Média populacional obtida=98
Desvio padrão populacional obtido=42
Média amostral obtida=98.539
Desvio padrão amostral obtido=38.822
Assimetria=0.604
ICBootPivotal95(media pop):
[ 84.00595 , 110.6758 ] (em relação à média: -14.53329 12.13654 )
IC95 da media populacional: bootstrapping por percentil
[ 86.43313 , 111.273 ], (em relação à média: -12.1061 12.7338 )
IC95 da media populacional: teste t para um grupo:
[ 85.40386 , 111.6746 ], (em relação à média: -13.13537 13.13537 )
(esta é a estimativa tradicional paramétrica, centrada na média)Neste gráfico e nos seguintes não apresentamos mais a distribuição
das médias amostrais (já sabemos que serão aproximadamente normais), mas
apenas os intervalos de confiança computados por três métodos
diferentes: o bootstrapping percentílico, o
bootstrapping pivotal (veja o código para seu cálculo) e a
maneira tradicional por equação com base na distribuição \(t\) (aproveitando a função nativa
t.teste que fornece o intervalo de confiança sem que
precisemos implementá-lo manualmente).
Note que a média populacional igual a 98.04 (linha
vertical pontilhada) está dentro do intervalo de confiança pivotal de
95% [ 84.01, 110.68] (o mesmo ocorreu com os outros dois
tipos de intervalo). Os intervalos obtidos por bootstrapping
pivotal podem ser assimétricos (capturando a assimetria da amostra e da
variável populacional), enquanto a estimativa tradicional do intervalo
de confiança da média é sempre simétrico. No caso dos intervalos de
confiança do desvio-padrão, por bootstrapping ou por fórmula,
assimetrias costumam aparecer.
Uma condição
Em outro exemplo estes procedimentos, como não poderia deixar de ser, funcionam da mesma forma para populações cuja variável tem distribuição aproximadamente normal.
Supondo que a estatura de homens adultos tem distribuição normal com média 171 cm e desvio-padrão de 7 cm e que tenhamos feito uma amostra, podemos verificar a plausibilidade da amostra ser oriunda desta população. A hipótese nula é:
\[ H_0: \mu = 171 \]
e alternativa:
\[ H_1: \mu \ne 171 \]
Testamos o procedimento executando simulaIC95pivotal2.R, obtendo:
Média populacional = 171
Desvio padrão populacional = 7
Média amostral obtida=171.781
Desvio padrão amostral obtido=6.486
Assimetria=-0.164
ICBootPivotal95(media pop):
169.4189 174.2576 (em relação à média: -2.361658 2.47695 )
IC95 da media populacional: bootstrapping por percentil
169.4806 174.0414 (em relação à média: -2.300055 2.260756 )
IC95 da media populacional: teste t para um grupo:
169.3588 174.2024 (em relação à média: -2.421816 2.421816 )
(esta é a estimativa tradicional paramétrica, centrada na média)Novamente, média populacional igual a 171 (linha
vertical pontilhada) está dentro dos intervalos de confiança 95%
estimados pelos três métodos apresentados.
Duas condições independentes
Os intervalos de confiança também podem ser utilizados para distinguir duas populações (Chihara & Hesterberg, 2019).
Por exemplo, imagine que não soubéssemos que os homens tendem a ser mais altos que as mulheres e que estaturas têm distribuição normal. Suponha que, em dada população, os homens têm \(\mu_M=175\) e \(\sigma_M=8\) e as mulheres \(\mu_M=168\) e \(\sigma_M=6\).
Sem ter acesso à população, simulamos com simulaIC95pivotal_MF.R a retirada de
uma amostra de cada grupo com \(n=30\)
para decidir se as estaturas coincidem, utilizando um teste da hipótese
nula. No gráfico aparece a distribuição normal populacional apenas para
referência (lembre que nunca a veríamos na prática), evidenciando que as
amostras, embora próximas em termos de média e dispersão, são
imperfeitas.
Caso a estatura dos indivíduos de sexo masculino e feminino fossem iguais, esperaríamos:
\[ H_0: \mu_M - \mu_F = 0 \]
A hipótese alternativa é:
\[ H_1: \mu_M - \mu_F \ne 0 \] Executamos o bootstrapping para cada amostra em separado a cada iteração, calculamos as respectivas médias e armazenamos a diferença entre elas. A distribuição das diferenças entre as médias, como se vê no gráfico abaixo, também produz uma distribuição aproximadamente normal (o TCL também vale nesta situação), com média igual à diferença das médias (amostrais ou populacionais) e erro-padrão da diferença entre as médias.
Masculino
Média populacional=175
Desvio padrão populacional=8
Feminino
Média populacional=168
Desvio padrão populacional=6
Masculino (n=30)
Média amostral obtida=174.545
Desvio padrão amostral obtido=8.998
Assimetria=-0.186
Feminino (n=30)
Média amostral obtida=167.817
Desvio padrão amostral obtido=5.81
Assimetria=0.352
ICBootPercentilico95(dif. medias pop):
2.892144 10.45776 (em relação à diferença das médias: -3.835895 3.729721 )
ICBootPivotal95(dif. medias pop):
2.84948 10.74495 (em relação à diferença das médias: -3.878559 4.016914 )
IC95 da dif. das medias populacionais
(teste t para condições independentes):
2.799474 10.6566 (em relação à diferença das médias: -3.928565 3.928565 )
(esta é a estimativa tradicional paramétrica, centrada na média)
Observe que as estimativas dos intervalos de confiança de 95% da diferença das médias populacionais entre indivíduos masculinos e femininos não incluiram o valor zero, esperado para quando não houvesse diferença de estatura. Portanto, rejeitamos a igualdade e ficamos com a hipótese alternativa.
Como fizemos a diferença das médias amostrais \(\bar{x}_M - \bar{x}_F\) e o intervalo de confiança de 95% está acima do valor nulo, então assuminos, com base neste teste, que a média de estatura dos indivíduos do sexo masculino é maior do que a dos indivíduos do sexo femininos.
Duas condições dependentes
Uma outra situação ocorre com condições dependentes. Suponha que mulheres façam parte de um estudo para avaliar a eficácia de uma dieta. Suas massas corporais totais (MCT), então, são obtidas antes e depois da dieta. Analisa-se a diferença de MCT (adotaremos “depois”-“antes”).
Aqui, a hipótese nula é
\[ H_0: \mu_D = 0 \]
e alternativa é
\[ H_1: \mu_D \ne 0 \]
onde \(\mu_D\) é a diferença média de peso. Caso a dieta funcione, a diferença média de peso “depois”-“antes” não pode ser nula. Caso seja negativa, a dieta promoveu emagrecimento; caso seja positivo, engorda.
Neste exemplo, aproximando de uma situação que podemos encontrar na
prática e diferentemente dos exemplos anteriores, não simulamos qualquer
referência populacional. O arquivo Dieta.xlsx contém as medidas de massa
corporal total em kg de 30 mulheres em dois momentos (Antes e Depois da
dieta) e a diferença de Depois-Antes em g. Observe que há mulheres que
perderam e outras que ganharam peso com a dieta.
Antes Depois Dif
66.3 66.5 200
57.0 55.3 -1700
51.9 52.4 500
54.8 53.7 -1100
74.7 74.1 -600
77.6 77.4 -200
63.7 64.0 300
59.4 60.2 800
55.2 53.7 -1500
66.5 66.4 -100
63.5 63.3 -200
68.8 66.9 -1900
64.2 63.5 -700
73.4 74.0 600
72.0 70.4 -1600
64.2 63.9 -300
78.0 77.8 -200
66.8 67.0 200
58.0 56.5 -1500
67.7 68.1 400
67.2 66.3 -900
70.1 70.2 100
62.7 62.0 -700
60.3 58.5 -1800
74.1 73.9 -200
58.1 58.3 200
57.6 58.2 600
67.9 66.5 -1400
55.2 55.0 -200
70.0 69.3 -700Simulamos com simulaIC95pivotal_Dieta.R. O primeiro
gráfico mostra a distribuição das massas corporais totais nos dois
momentos considerados; são, visualmente similares.
O segundo gráfico mostra a diferença de massa corporal total, em
gramas, sobre a qual é feito o teste. A curva é a distribuição amostral,
anotada na terceira coluna do data frame,
dt_dieta$Dif que, embora mostre um emagrecimento de -453g
em média, 10 mulheres ganharam peso enquanto outras 20 perderam peso.
Precisamos avaliar se esta dieta funciona.
Os intervalos de confiança no topo do gráfico foram obtidos pelo teste \(t\) tradicional ou por bootstrapping (veja o código para entender sua aplicação neste caso).
n=30
Antes
Média amostral obtida=64.897kg
Desvio padrão amostral obtido=7.105kg
Assimetria=0.051
Depois
Média amostral obtida=64.443kg
Desvio padrão amostral obtido=7.203kg
Assimetria=0.062
ICBootPercentilico95(diferença pop):
-740 -173.3333 (em relação à diferença média: -286.6667 280 )
ICBootPivotal95(diferença pop):
-738.3927 -132.6996 (em relação à diferença média: -285.0593 320.6338 )
IC95 da media populacional: teste t para um grupo:
-753.4606 -153.2061 (em relação à diferença média: -300.1273 300.1273 )
(esta é a estimativa tradicional paramétrica, centrada na média)
O intervalo de confiança de 95% da diferença de MCT não contém o valor nulo e o intervalo é negativo. Rejeita-se a hipótese nula e, portanto, o experimento sugere que houve redução da MCT com a dieta.
Proporção
Em epidemiologia, a prevalência (proporção) das doenças na população é uma quantidade fundamental.
Prevalência é a proporção de indivíduos doentes em determinado momento e, portanto, pode ser vista como a probabilidade de encontrarmos um indivíduo doente através de sorteio simples.
Doença é um evento raro na população geral.
Prevalência é probabilidade (e.g., 2% de gripe em agosto/2018):
- pode ser medida em uma coorte fechada ou em uma população aberta
- geralmente usada em estudos transversais
- varia entre 0 e 1
|
Não confunda: Incidência é taxa (e.g., casos por 100 mil habitantes por ano):
Razão de Riscos (Risk Ratio (RR)) é razão entre duas probabilidades:
|
van Belle, 2008
É muito comum existirem itens com respostas dicotômicas, por exemplo quando classifica-se um indivíduo em portador ou não de uma doença, transtorno, trauma ou síndrome. Uma das formas de lidar com isto é utilizar a proporção destes eventos de interesse.
O modelo de variável dicotômica é o jogo de cara ou coroa com uma moeda. Assim, se uma moeda bem balanceda for jogada 100 vezes, esperamos obter 50 coroas. Em outras palavras, a proporção de coroas é 0.5 ou 50%, que também pode ser vista como a probabilidade de ocorrência de coroa, ou como a média esperada de coroas.
A distribuição de probabilidades de uma variável dicotômica é a de Bernoulli (um evento com probabilidade \(p\)). A soma de vários eventos deste tipo resulta em uma distribuição binomial. Assim, a média da ocorrência de \(n\) eventos é \(np\) e o desvio-padrão é \(\sqrt{np(1-p)}\).
Na área da saúde \(p\) corresponde também à prevalência. Como vimos, os valores de \(p\) costumam ser pequenos, da ordem de 10% ou (muito) menos. Em estudos, como sempre, usamos amostras e, portanto, precisamos calcular o erro-padrão e o intervalo de confiança da proporção populacional. O erro-padrão é dado por:
\[ \text{EP} = \sqrt{{\hat{p}(1-\hat{p})\over{n}}} \]
O intervalo de confiança de proporção de Wald é o mais popular. Para o IC de 95%, calcula-se com:
\[ \begin{align} \text{IC}^{95\%}_{\text{Wald}}(p) &= \left[\hat{p} \pm 1.96 \times \text{EP}\right]\\ \text{IC}^{95\%}_{\text{Wald}}(p) &= \left[\hat{p} \pm 1.96 \sqrt{{\hat{p}(1-\hat{p})\over{n}}}\right] \end{align} \]
A fórmula de Wald, porém, tem restrições. Para utilizar a distribuição normal padrão, onde \(z_{\alpha/2}=1.96, \alpha=0.05\), o tamanho da amostra deveria ser de, pelo menos, 30 indivíduos. Além disto, proporções estão limitadas ao intervalo \([0,1]\) e esta fórmula não garante que seus limites inferior e superior estejam contidos neste intervalo.
Uma forma mais robusta é o Intervalo de Confiança de Proporção de Wilson. Para 95% é dado por:
\[ \text{IC}^{95\%}_{\text{Wilson}}(p) = \left[\dfrac{{{\hat{p} + \dfrac{1.96^2}{2n} \pm 1.96 \sqrt{\dfrac{\hat{p}(1-\hat{p})}{n} +\dfrac{1.96^2}{4n^2}}}}}{1 + {\dfrac{1.96^2}{n}}}\right] \]
Wilson, 1927
Esta fórmula vale para qualquer \(n\) e também costuma manter os limites inferior e superior dentro do intervalo [0,1].
Wallis (2013) traz um exemplo verificando a proporção do uso das palavras shall ou will na língua inglesa ao longo dos anos.
Podemos comparar o desempenho das duas fórmulas com as seguintes figuras:
- com o intervalo de Wald:
- com o intervalo de Wilson:
Wallis, 2013
| Se quiser saber mais, consulte Robert, 2020. |
Uma condição: teste bilateral
Simulamos com demo_binomial.R a distribuição de
ocorrências de um evento com prevalência populacional de 10%, do qual
retiramos uma amostra de 30 indivíduos. As variantes dos intervalos de
confiança de Wald e Wilson são denominados no gráfico a seguir por teste
t e Binomial, respectivamente.
Amostra:
2 sucessos em 30 tentativas
media: 0.0667
desvio-padrão: 0.2537
erro-padrao: 0.0463
assimetria: 3.3021
Intervalos:
IC95 bootstrapping:
[0.0000,0.1667]
IC95 teste binomial:
[0.0082,0.2207]
IC95 teste t:
[-0.0281,0.1614] Observe que o teste \(t\), não concebido para proporções, pode gerar valores inadequados. Note, também, que os intervalos de confiança obtidos por bootstrapping e pelo teste binomial são assimétricos.
É importante enfatizar que o valor da prevalência de 10%, representada como uma linha horizontal pontilhada no gráfico, só está disponível porque o simulamos. Na prática, não saberemos a prevalência populacional: aqui está representada para mostrar que os intervalos de confiança estimados a partir da única amostra que obtivemos contém (e, portanto, estimaram corretamente) o valor populacional.
Uma condição: teste unilateral
Há vezes em que, para estimar a prevalência de uma doença na
população, queremos saber qual é seu valor mínimo ou máximo. Para tanto,
generalizamos este código para utilizar o parâmetro side
com um dos três valores:
two.sidedque é a situação anterior (default),greaterque estima a prevalência mínima,lessque estima a prevalência máxima
Estas são estimativas do mínimo e máximo da média populacional de onde a amostra proveio (Wonnacott & Wonnacott, 1990, p. 317-8).
Então, a estimativa da prevalência mínima é o limite mínimo do intervalo obtido por:
Amostra:
2 sucessos em 30 tentativas
media: 0.0667
desvio-padrão: 0.2537
erro-padrao: 0.0463
assimetria: 3.3021
Intervalos:
IC95 bootstrapping:
[0.0000,NA]
IC95 teste binomial:
[0.0120,NA]
IC95 teste t:
[-0.0120,NA] Observe novamente o teste \(t\) gerando probabilidades negativas. Além disto, como estamos procuramos a estimativa mínima da prevalência, o valor superior não interessa: e o código está adequado para auxiliar o usuário substituindo o valor (ainda que apareça) por NA e representando o intervalo com uma seta direcional para sugerir que prossegue para cima (até 1, no caso de proporções).
A estimativa da prevalência máxima é o limite máximo do intervalo obtida por:
Amostra:
2 sucessos em 30 tentativas
media: 0.0667
desvio-padrão: 0.2537
erro-padrao: 0.0463
assimetria: 3.3021
Intervalos:
IC95 bootstrapping:
[NA,0.1333]
IC95 teste binomial:
[NA,0.1953]
IC95 teste t:
[NA,0.1454] Novamente, agora para o limite máximo, o código foi preparado para auxiliar a análise.
Nestes três exemplos, fixamos a amostra com set.seed. A
estimativa pontual desta amostra subestimou a prevalência de nossa
população artificial (encontrou 2 sucessos em 30 tentativas, em vez de 3
que esperaríamos para a prevalência de 10% que usamos para criar a
população). No entanto, estimando o intervalo de confiança, a decisão
passa a ser sobre a faixa de prevalência dentro da qual deve estar a
prevalência populacional; como são dados simulados, sabemos que o
intervalo estimado está correto porque englobou 10%.
|
Há outras funções em outros pacotes que também dão informações interessantes sobre proporções. Há vários outros métodos possíveis, disponíveis em outros pacotes.
Por exemplo, Com o código disponível em
Compare os intervalos obtidos por Para as estimativas unilaterais implementadas neste código, basta
passar um dos três valores disponíveis em
O equivalente a
|
Duas condições independentes
Suponha que uma dieta foi instituída para pessoas dos dois sexos, 80 pessoas do sexo feminino (F) e 70 pessoas do sexo masculino (M). Após algum tempo, o pesquisador verifica quantos tiveram sucesso em perder peso: 48 (60%) mulheres e 56 (80%) homens conseguiram emagrecer. Como verificar se a proporção de sucesso é igual ou diferente para os dois sexos?
As hipóteses são:
\[ \begin{cases} H_0: \pi_F - \pi_M = 0 \\ H_1: \pi_F - \pi_M \ne 0 \end{cases}\\ \alpha = 5\% \]
sendo que \(\pi\) é a proporção populacional da condição (sucesso em emagrecer).
O código R para este tipo de teste utiliza a função
DescTools::BinomDiffCI:
sucessos.F <- 48
n.F <- 80
sucessos.M <- 56
n.M <- 70
xci <- DescTools::BinomDiffCI(sucessos.F, n.F,
sucessos.M, n.M)
print(data.frame(xci), digits=2, row.names=FALSE) est lwr.ci upr.ci
-0.2 -0.34 -0.052Observe que o valor nulo não está contido no intervalo de confiança de 95% da diferença de proporções de emagrecimento (pela ordem) das mulheres em relação aos homens. Como o intervalo está abaixo de zero, significa que a proporção de mulheres que teve sucesso em emagrecer é menor do que a dos homens.
|
A função Não encontrei na documentação nesta data (setembro/2020) qual é o
método default. Quando isto acontece, pode-se usar o R como
laboratório para descobrir. Implementamos uma generalização do código
anterior em
Comparando a saída com a anterior, vemos que o método usado como
default foi sides também está
disponível, caso precise.
|
Três ou mais condições independentes
Em um estudo verificou-se a curva de aprendizado de neurologistas,
medindo-se a proporção de erros cometidos em punções lombares. A
expectativa é de decréscimo da proporção de erros ao longo de três anos
de treinamento (R1 a R3). Como referência, o estudo incluiu os médicos
especialistas da área (assistentes). A principal função que escolhemos
para resolver esta situação é DescTools::BinomCI para o
teste omnibus e fmsb::pairwise.fisher.test para o
teste post hoc.
Se os intervalos de confiança com correção de Bonferroni dois a dois não se sobrepõem, então há diferença significante. Se há sobreposição, o valor p como correção de Bonferroni pode ser usado para analisar a significância par a par.
Os dados estão simulados em demo_binomialCondicoes.R:
Proporcao de erros de puncao:
R1: 24 em 61 = 39.34%
R2: 38 em 76 = 50%
R3: 8 em 55 = 14.55%
Assistente: 1 em 15 = 6.67%
Intervalos de confianca 95% Bonferroni:
est lwr.ci upr.ci
R1 0.393 0.2542 0.55
R2 0.500 0.3623 0.64
R3 0.145 0.0634 0.30
Assistente 0.067 0.0083 0.38
Teste omnibus:
Results of Hypothesis Test
--------------------------
Alternative Hypothesis: two.sided
Test Name: 4-sample test for equality of proportions without continuity correction
Estimated Parameter(s): prop 1 = 0.393
prop 2 = 0.500
prop 3 = 0.145
prop 4 = 0.067
Data: errospuncao out of tentativaspuncao
Test Statistic: X-squared = 24
Test Statistic Parameter: df = 3
P-value: 3e-05
Teste post hoc:
R1 R2 R3
R2 1.000 NA NA
R3 0.021 0.0002 NA
Assistente 0.094 0.0093 1
Três ou mais condições dependentes (multinomial)
Três marcas de chocolate foram disponibilizadas para degustação às
cegas por 79 voluntários. Cada voluntário experimentou os três
chocolates e então informou seu preferido. O chocolate A teve 23
(29.1%), o chocolate B 12 (15.2%) e o chocolate C 44 (55.7%) das
preferências. As marcas, então foram reveladas, respectivamente Best
Cocoa, Dream Brown e Wonka. A principal função que escolhemos para
resolver esta situação é DescTools::MultinomCI (default:
sisonglaz) para o teste omnibus e
RVAideMemoire::multinomial.multcomp para o teste post
hoc.
Se os intervalos de confiança com correção de Bonferroni dois a dois não se sobrepõem, então há diferença significante. Se há sobreposição, o valor p como correção de Bonferroni pode ser usado para analisar a significância par a par.
A análise está em demo_binomialChocolate.R:
est lwr.ci upr.ci
Dream Brown 0.15 0.025 0.28
Best Cocoa 0.29 0.165 0.42
Wonka 0.56 0.430 0.69
----------------------------
Teste omnibus:
Results of Hypothesis Test
--------------------------
Alternative Hypothesis:
Test Name: Chi-squared test for given probabilities with simulated p-value
(based on 1e+05 replicates)
Data: preferido
Test Statistic: X-squared = 20.07595
Test Statistic Parameter: df = NA
P-value: 2.99997e-05
Infere-se, para a populacao, que a preferencia por pelo menos um dos chocolates difere de algum outro.
----------------------------
Analise par a par:
Dream Brown Best Cocoa
Best Cocoa 2.7e-01 NA
Wonka 6.3e-05 0.042
Infere-se, para a populacao, entre:
- Best Cocoa e Dream Brown: nao ha diferenca de preferencia entre os chocolates
- Wonka e Dream Brown: ha diferenca de preferencia entre os chocolates
- Wonka e Best Cocoa: ha diferenca de preferencia entre os chocolatesComo no exemplo anterior, o teste estatístico é aplicado em duas fases: (1) aplicamos um qui-quadrado para avaliar globalmente (teste omnibus) se há alguma diferença entre os chocolates e (2) havendo, compara-se os chocolates dois a dois. A diferença, neste caso, está nas funções escolhidas porque há dependência sobre as escolhas dos chocolates: ao escolher um chocolate, cada indivíduo deixa de escolher os outros dois, criando uma correlação negativa entre as escolhas e as recusas.
A hipótese nula, par a par, é que a proporção de preferência pelos
chocolates é a mesma. Adotando \(\alpha=0.05\), conclui-se que não se
rejeita a hipótese nula na comparação entre Dream Brown e Best Cocoa
(\(p\)=0.2685932), mas que
o Wonka é diferente dos outros dois (respectivamente \(p\)=0.0000626297 e \(p\)=0.04180163). Observando as
proporções, Wonka é o preferido.
Intervalo de credibilidade (bayesiano)
Suposições para validade do intervalo de credibilidade:
- Delineamento entre participantes: observações independentes
- Variável intervalar
- Distribuição normal da variável
- O parâmetro média da distribuição normal tem distribuição normal
Na inferência frequentista/ clássica/ de Neyman-Pearson, o intervalo de confiança de proporção da média é um intervalo que contêm este parâmetro com 95% de confiança (https://rpsychologist.com/d3/ci/).
“Um intervalo de confiança frequentista de 95% significa que, com um grande número de amostras repetidas, 95% de tais intervalos de confiança calculados incluem o valor verdadeiro do parâmetro. Em termos frequentistas, o parâmetro é fixo (não pode ser considerado como tendo uma distribuição de valores possíveis) e o intervalo de confiança é probabilístico (pois depende da amostra aleatória).”
Na inferência bayesiana, o intervalo de credibilidade de 95% da média é um intervalo que contêm este parâmetro com 95% de probabilidade (https://rpsychologist.com/d3/bayes/).
O preço do intervalo de credibilidade bayesiano é que sua interpretação depende da distribuição prévia (prior) que é subjetiva.
“Um intervalo de credibilidade de 95% é um intervalo dentro do qual um valor de parâmetro com distribuição de probabilidade está com 95% de probabilidade.”
Intervalo de credibilidade de 95% da média populacional:
\[ \Large \text{IC}_{\text{bayes}}^{0.95}(\mu)=\left[\frac{n_0\,\mu_0+n\,\bar{x}}{n_0+n} \pm t^{0.975}_{n-1}\,\frac{s}{\sqrt{n_0+n}}\right] \]
Sendo que:
- \(\mu \sim \text{normal}\left(\mu_0,\sigma_0\right)\)
- \(n_0 = \dfrac{s}{\sigma_0}\)
Exemplo: MCT Feminino
\[ \Large \text{IC}^{0.95}_{\text{bayes}}(\mu)=\left[\frac{1.82\times55+230\times57.6}{1.82+230} \pm 1.97\,\frac{9.1}{\sqrt{1.82+230}}\right] \]
Sendo que:
- \(t^{0.975}_{229}=\) 1.9702867
- \(\mu \sim \text{normal}\left(55,5\right)\)
- \(n_0 = \dfrac{9.1}{5} = 1.82\)
\[ \Large \text{IC}^{0.95}_{\text{bayes}}(\mu)=\left[56.4,58.8\right] \]
\[ \Large \text{IC}^{0.95}_{\text{freq}}(\mu)=\left[56.5,58.8\right] \]
\[ \Large \text{IC}^{0.95}_{\text{boot}}(\mu)=\left[56.6, 58.9\right] \]
Se o tamanho da amostra é grande, i.e., \(n>30\), os três intervalos são aproxidamentamente iguais.
A distribuição não-informativa da média populacional \(\mu\) ocorre quando o desvio-padrão \(\sigma_0\) é grande, acarretatando valor baixo de \(n_0\).
Se \(n\gg n_{0}\), o intervalo de confiança é aproximadamente igual ao intervalo de credibilidade.
“Para nós, a mensagem principal do nosso artigo é a seguinte. Os intervalos de confiança frequentistas podem ser interpretados como uma aproximação razoável de um intervalo de credibilidade bayesiano (com a priori não informativa). Isso é reconfortante para aqueles que lutam por a interpretação formalmente correta dos intervalos de confiança frequentistas.”
Albers et al., 2018, p. 6
Wonnacott & Wonnacott, 1969
“Ponto 11: Fatores de Bayes frequentemente concordam com os valores p […] Os resultados entre fatores de Bayes e testes clássicos frequentemente concordam entre si.”
Tendeiro & Kiers, 2019, p. 789
Referências:
- Wonnacott, TH & Wonnacott, RJ (1969) Introductory Statistics. NJ: Wiley.
- Wonnacott, TH & Wonnacott, RJ (1981) Estatística aplicada à Economia e à Administração. RJ: LTC.
- Wonnacott, TH & Wonnacott, RJ (1990) Introductory Statistics for Business and Economics. NJ: Wiley.
- Albers, CJ et al. (2018) Credible confidence: A pragmatic view on the frequentist vs Bayesian debate. Collabra: Psychology, 4(1): 31. DOI: https://doi.org/10.1525/collabra.149; Scripts R: https://osf.io/dgfht/.
- Tendeiro, JN & Kiers, HAL (2019) A review of issues about null hypothesis Bayesian testing. Psychological Methods 24(6): 774–95. https://doi.org/10.1037/met0000221.
- R-Bloggers: Confidence vs. Credibility Intervals: https://www.r-bloggers.com/2014/11/confidence-vs-credibility-intervals/
- https://en.wikipedia.org/wiki/Credible_interval
- https://rpsychologist.com/d3/ci/
- https://rpsychologist.com/d3/bayes/
Vídeos
- Incidence and Prevalence - Everything you need to know: https://www.youtube.com/watch?v=cTp_ONVVrh8
- Test sensitivity and specificity made easy: https://www.youtube.com/watch?v=8J_i2C4elnk
- Sensitivity and Specificity – Advanced: https://www.youtube.com/watch?v=eC1K_2DR9Yg
- Positive Predictive Value - The role of specificity: https://www.youtube.com/watch?v=z4yauM3GqfI
- Positive Predictive Value - The role of prevalence: https://www.youtube.com/watch?v=QqgJHryKOSU
Apêndice: Intervalos de Predição
Seja \(Y = f(X_1, X_2, \ldots, X_k)\), \(k\in \mathbb{N}^{\ast}\), uma função de variáveis aleatórias. Portanto, \(Y\) também é uma variável aleatória.
Se \(k=1\), \(Y\) tem esperança e variância, exceto a distribuição de Cauchy (t de Student com um grau de liberdade) que não tem esperança e variância e t de Student com dois graus de liberdade que tem esperança, mas não tem variância.
Se \(k>1\) e as variáveis da função são independentes, então é possível, em geral, determinar esperança e variância exatas por meio das transformadas de Mellin e de Laplace mesmo desconhecendo a distribuição de \(Y\), mas conhecendo as distribuições de \(X_1, X_2, \ldots, X_k\).
A esperança e variância aproximadas podem ser obtidas por expansão de Taylor.
Se \(k>1\) e as variáveis da função são dependentes, então não é possível determinar esperança e variância exatas por meio das transformadas de Mellin e de Laplace. Para algumas funções específicas de variáveis aleatórias dependentes é possível a determinação exata da esperança e da variância. Se as variáveis aleatórias são dependentes e as correlações são conhecidas, a esperança e variância aproximadas podem ser obtidas por expansão de Taylor.
A obtenção da esperança e variância de função de variáveis aleatórias é importante para a determinação do coeficiente de variação de Pearson (desvio-padrão/esperança) e para o intervalo de predição por meio da desigualdade de Chebychev, desigualdade de Selberg e desigualdade de Gauss-Camp-Meidell.
Desigualdade de Chebychev
A desigualdade de Chebychev fornece um limite para a probabilidade de uma variável aleatória \(X\) se desviar da sua média \(\mu\) por mais de \(k\) vezes o desvio-padrão \(\sigma\). Isso é útil porque a desigualdade é válida para qualquer distribuição de probabilidade com variância finita, independentemente da forma dessa distribuição.
Forma da Desigualdade
Para uma variável aleatória \(X\) com média \(\mu\) e variância \(\sigma^2\), a desigualdade de Chebychev pode ser expressa da seguinte forma:
\[ P(|X - \mu| \geq k\sigma) \leq \dfrac{1}{k^2} \]
Podemos reescrever essa expressão para evitar o uso do módulo, indicando um intervalo de predição simétrico em torno da média:
\[ P(\mu - k\sigma \leq X \leq \mu + k\sigma) \geq 1 - \dfrac{1}{k^2} \]
Interpretação
Esta forma indica que a probabilidade de \(X\) estar dentro do intervalo \(\mu \pm k\sigma\) é, no mínimo, \(1 - \dfrac{1}{k^2}\). Em outras palavras, pelo menos \(1 - \dfrac{1}{k^2}\) da distribuição da variável \(X\) está contida nesse intervalo.
Exemplos
Para valores específicos de \(k\):
Para \(k = 2\): \[P(\mu - 2\sigma \leq X \leq \mu + 2\sigma) \geq 1 - \frac{1}{4} = 0.75\] Ou seja, pelo menos 75% dos valores de \(X\) estão dentro de \(2\) desvios-padrão da média.
Para \(k = 3\): \[P(\mu - 3\sigma \leq X \leq \mu + 3\sigma) \geq 1 - \frac{1}{9} \approx 0.89\] Pelo menos 89% dos valores de \(X\) estão dentro de \(3\) desvios-padrão da média.
Utilidade
A desigualdade de Chebychev é particularmente útil porque:
- Independência da Distribuição: Ela não assume que a distribuição de \(X\) segue qualquer forma específica, como a normal. Isso a torna aplicável em uma ampla gama de situações.
- Estimativa Conservadora: Fornece uma estimativa conservadora que pode ser utilizada mesmo quando pouco se sabe sobre a forma da distribuição.
- Aplicação Geral: Pode ser aplicada em diversas áreas, como controle de qualidade, análise de risco e qualquer situação onde a variância é conhecida e se deseja estimar a concentração dos valores em torno da média.
| \(k\) | Mínima porcentagem dentro de \(k\) desvios-padrão da média |
|---|---|
| 1 | 0% |
| \(\sqrt{2}\) | 50% |
| 1.5 | 55.56% |
| 2 | 75% |
| \(2\sqrt{2}\) | 87.5% |
| 3 | 88.89% |
| 4 | 93.75% |
| 5 | 96% |
| 6 | 97.22% |
| 7 | 97.96% |
| 8 | 98.44% |
| 9 | 98.77% |
| 10 | 99% |
Desigualdade de Selberg
A desigualdade de Selberg é uma generalização da desigualdade de Chebyshev para intervalos arbitrários. Suponha que \(X\) seja uma variável aleatória com média \(\mu\) e variância \(\sigma^2\). A desigualdade de Selberg estabelece que, se $ $, então:
\[ P(X \in [\mu - \alpha, \mu + \beta]) \geq \begin{cases} \dfrac{\alpha^2}{\alpha^2 + \sigma^2} & \text{se } \alpha(\beta - \alpha) \geq 2\sigma^2 \\ \dfrac{4\alpha\beta - 4\sigma^2}{(\alpha + \beta)^2} & \text{se } 2\alpha\beta \geq 2\sigma^2 \geq \alpha(\beta - \alpha) \\ 0 & \text{se } \sigma^2 \geq \alpha\beta \end{cases} \]
Quando \(\alpha = \beta\), esta desigualdade se reduz à desigualdade de Chebyshev. Esses limites são conhecidos por serem os melhores possíveis.
A desigualdade de Selberg fornece uma maneira de estimar a probabilidade de uma variável aleatória \(X\) cair dentro de um intervalo específico e assimétrico em torno de sua média, levando em consideração a variância da variável. Vamos aplicar essa desigualdade em um exemplo prático em biologia.
Exemplo Prático em Biologia
Suponha que estamos estudando o comprimento das asas de uma espécie de pássaro. Vamos dizer que o comprimento das asas, \(X\), é uma variável aleatória com média \(\mu = 15\) cm e variância \(\sigma^2 = 4\) cm².
Queremos saber a probabilidade de que o comprimento das asas esteja dentro do intervalo [13 cm, 18 cm]. Aqui, temos:
- \(\mu = 15\)
- \(\alpha = \mu - 13 = 15 - 13 = 2\)
- \(\beta = 18 - \mu = 18 - 15 = 3\)
- \(\sigma^2 = 4\)
Vamos verificar as condições para aplicar a desigualdade de Selberg:
- \(\alpha(\beta - \alpha) = 2(3 - 2) = 2\)
- \(2\sigma^2 = 2 \times 4 = 8\)
- \(2\alpha\beta = 2 \times 2 \times 3 = 12\)
- \(\alpha\beta = 2 \times 3 = 6\)
Agora, aplicamos a desigualdade de Selberg:
- \(\alpha(\beta - \alpha) = 2\) e \(2\sigma^2 = 8\). Portanto, a condição \(\alpha(\beta - \alpha) \geq 2\sigma^2\) não é satisfeita.
- \(2\alpha\beta = 12 \geq 8 \geq 2\), e a condição \(2\alpha\beta \geq 2\sigma^2 \geq \alpha(\beta - \alpha)\) é satisfeita.
Portanto, usamos a segunda fórmula da desigualdade de Selberg:
\[ P(X \in [13, 18]) \geq \dfrac{4\alpha\beta - 4\sigma^2}{(\alpha + \beta)^2} \]
Substituindo os valores:
\[ \begin{align} P(X \in [13, 18]) &\geq \dfrac{4 \times 2 \times 3 - 4 \times 4}{(2 + 3)^2}\\ &\geq \dfrac{24 - 16}{25}\\ &\geq \dfrac{8}{25}\\ P(X \in [13, 18]) &\geq 0.32 \end{align} \]
Assim, a probabilidade de que o comprimento das asas dos pássaros esteja dentro do intervalo de 13 cm a 18 cm é pelo menos 32%.
Interpretação
Essa aplicação da desigualdade de Selberg em biologia nos permite estimar a probabilidade de uma característica biológica (como o comprimento das asas de pássaros) cair dentro de um intervalo específico, considerando a variância da medida.
Desigualdade de Gauss-Camp-Meidell
A desigualdade é uma variação da desigualdade de Chebyshev que pode ser usada para distribuição unimodal assimétrica de variável intervalar.
Desigualdade de Gauss-Camp-Meidell
Para uma distribuição unimodal assimétrica, a desigualdade de Gauss-Camp-Meidell estabelece que um intervalo simétrico centrado na média, com semi-amplitude \(k\sigma\), sendo que \(k\) é dado por:
\[ k=\gamma + \dfrac{2}{3} \sqrt{\frac{1 - \gamma^2}{1 - \pi}} \]
contém pelo menos \(\pi\)% dos valores de uma variável intervalar unimodal assimétrica.
A fórmula do intervalo de predição é:
\[ IP_\pi(X) = \left[ \mu \pm k\sigma \right] \]
onde:
- \(\sigma\): desvio-padrão
- \(\mu\): média
- \(\beta_1=\dfrac{\mu_3}{\sigma^{3/2}}\): indice de assimetria
- \(\beta_2=\dfrac{\mu_4}{\sigma^{4}}\): índice de curtose
- \(\zeta=\mu-\dfrac{\sigma}{2}\dfrac{\beta_1\left(\beta_2+3\right)}{5\beta_2-6\beta_1^2-9}\): moda (Park & Sharp-Bete, 1990, p. 421)
- \(\gamma = \dfrac{|\mu-\zeta|}{\sigma}\): coeficiente de assimetria
- \(\pi\) é a porcentagem mínima de valores contidos no intervalo.
Exemplo
Suponhamos que estamos estudando a massa corporal total de uma espécie de peixe em um lago.
Dados do Exemplo:
- Média (\(\mu\)): 26 gramas
- Variância (\(\sigma^2\)): 88 gramas²
- Terceiro momento central (\(\mu_3\)): 131 gramas³
- Quarto momento central (\(\mu_4\)): 15520 gramas⁴
Cálculos Iniciais:
O desvio-padrão (\(\sigma\)) é a raiz quadrada da variância:
\[ \sigma = \sqrt{88} \approx 9.38 \]
Índice de Assimetria (\(\beta_1\)):
\[ \beta_1 = \dfrac{\mu_3}{\sigma^3} = \dfrac{131}{9.38^3}\approx 0.159 \]
Índice de Curtose (\(\beta_2\)):
\[ \beta_2 = \dfrac{\mu_4}{\sigma^4} = \dfrac{15520}{9.38^4}\approx 2.00 \]
Cálculo da Moda (\(\zeta\)):
\[ \zeta = \mu - \dfrac{\sigma}{2} \dfrac{\beta_1 (\beta_2 + 3)}{5 \beta_2 - 6 \beta_1^2 - 9} \approx 21.60 \]
Cálculo do Coeficiente de Assimetria (\(\gamma\)):
\[ \gamma = \dfrac{|\mu - \zeta|}{\sigma} = \dfrac{|26 - 21.60|}{9.38} = \dfrac{4.40}{9.38} \approx 0.469 \]
Cálculo de \(k\) para \(\pi = 0.75\):
\[ k = \gamma + \dfrac{2}{3} \sqrt{\dfrac{1 - \gamma^2}{1 - \pi}} \]
Substituindo os valores:
\[ k = 0.469 + \dfrac{2}{3} \sqrt{\dfrac{1 - 0.469^2}{1 - 0.75}}\approx 1.65 \]
Intervalo de Predição (IP) para 75%:
\[ IP_{75}(X) = \left[ \mu \pm k\sigma \right] = \left[ 26 \pm 1.65 \times 9.38 \right] = \left[ 10.53, 41.47 \right] \]
Comparação com a Desigualdade de Chebyshev:
Para a desigualdade de Chebyshev, com \(p = 75\%\) (\(p = 0.75\)):
\[ \begin{align} 1 - \dfrac{1}{k^2} &= 0.75 \\ k &= 2 \end{align} \]
O intervalo de Chebyshev é:
\[ IP_{75}(X) = \left[ \mu \pm 2\sigma \right] = \left[ 26 \pm 2 \times 9.38 \right] \]
Portanto, o intervalo é \(\left[ 7.24, 44.76 \right]\).
Conclusão:
Intervalo de Gauss-Camp-Meidell para 75% dos valores: [10.53, 41.47]
Intervalo de Chebyshev para 75% dos valores: [7.24, 44.76]
O intervalo de Gauss-Camp-Meidell é mais estreito do que o intervalo de Chebyshev, refletindo uma estimativa mais precisa para distribuições unimodais assimétricas. Este exemplo biológico demonstra como a desigualdade de Gauss-Camp-Meidell pode fornecer intervalos de predição mais precisos para dados de massa de peixes, considerando a assimetria da distribuição.
Apêndice: Intervalo de referência
Os intervalos de referência constituem a base dos exames laboratoriais e auxiliam o médico a diferenciar indivíduos saudáveis de indivíduos com doença.
Um exemplo clássico é o sódio sérico em adultos. Um laboratório pode adotar, por exemplo, o seguinte intervalo de referência:
\[ RI_{95%}=[135,;145]\ \text{mmol/L} \]
Isso significa que aproximadamente 95% dos indivíduos da população de referência apresentam concentrações de sódio sérico entre 135 e 145 mmol/L, considerando o método analítico e a população utilizados pelo laboratório.
Por exemplo,
\[ Na^+=140\ \text{mmol/L} \]
está dentro do intervalo de referência, enquanto
\[ Na^+=150\ \text{mmol/L} \]
está acima do limite superior.
O intervalo de referência não deve ser interpretado como um limite diagnóstico. Seus valores podem variar entre laboratórios, métodos analíticos e populações de referência.
Conforme Horn & Pesce (2003), O problema enfrentado pelo profissional do laboratório é estabelecer intervalos de referência para populações saudáveis. Devem ser considerados os efeitos da idade, raça, exercício físico intenso, dieta ou condições de saúde não normais sobre a estimativa do intervalo de referência da população saudável.
O problema é duplo. Primeiro, o intervalo de referência é determinado a partir das amostras disponíveis em determinado laboratório. Essas amostras podem ou não ser constituídas exclusivamente por indivíduos saudáveis. Segundo, o tamanho da amostra pode não ser suficiente para estimar os limites do intervalo de referência com um grau razoável de precisão.
Como consequência, pode ocorrer aumento no número de decisões incorretas, levando a maiores custos, investigações desnecessárias e riscos à segurança do paciente.
Conforme Doyle Y Bunch (2023),
“RESUMO
Os resultados de exames laboratoriais clínicos, isoladamente, têm pouco valor para o diagnóstico, tratamento e monitoramento das condições de saúde; portanto, é necessário um ponto de corte clinicamente acionável ou um intervalo de referência que forneça contexto para a interpretação dos resultados. Os profissionais de saúde baseiam seus diagnósticos, tratamentos subsequentes e solicitações de novos exames nesses referenciais. Entretanto, podem não estar cientes das limitações inerentes à definição e à obtenção dos intervalos de referência.
Os profissionais de laboratório são responsáveis por fornecer os intervalos de referência apresentados juntamente com os resultados. Contudo, o estabelecimento e a verificação desses intervalos por métodos diretos convencionais são dificultados por limitações de recursos ou por características demográficas específicas dos pacientes. Para promover práticas padronizadas e adequadas para intervalos de referência, diversas sociedades científicas internacionais estão elaborando ativamente diretrizes e buscando financiamento para apoiar essas iniciativas.
Várias colaborações multicêntricas nacionais e internacionais demonstram a possibilidade de combinar recursos para realizar grandes estudos de intervalos de referência por métodos diretos. Entretanto, nem todos os grupos demográficos são igualmente acessíveis. Métodos indiretos inovadores constituem alternativas atraentes, pois utilizam métodos computacionais para definir distribuições e intervalos de referência a partir de conjuntos de dados mistos, contendo resultados de exames de pacientes com e sem condições patológicas.
Na busca por intervalos de referência mais precisos e personalizados, os intervalos de referência individualizados têm demonstrado ser mais úteis que os intervalos baseados na população para detectar alterações clinicamente significativas de analitos em um paciente, alterações que poderiam passar despercebidas quando se utilizam intervalos populacionais mais amplos.
Além disso, intervalos de referência contínuos podem fornecer faixas mais precisas do que partições baseadas em idade, especialmente para indivíduos próximos aos limites das faixas etárias estabelecidas. Nesta revisão, são discutidas as vantagens e desvantagens das diferentes abordagens para determinação de intervalos de referência, bem como os avanços nos estudos não convencionais desses intervalos. “.
Em aplicações clínicas e laboratoriais é comum estabelecer uma faixa de valores considerada compatível com uma determinada população de referência. Embora expressões como normal range, reference range, intervalo probabilístico e intervalo de predição sejam por vezes utilizadas de maneira pouco rigorosa, esses conceitos não são equivalentes.
Neste texto são considerados três tipos de intervalos: o intervalo probabilístico obtido pela desigualdade de Chebyshev, o intervalo de referência central baseado em quantis e o intervalo de maior densidade, ou Highest Density Interval (HDI).
Intervalo probabilístico pela desigualdade de Chebyshev
Seja uma variável aleatória \(X\), com média finita \(\mu\) e variância finita \(\sigma^2\). A desigualdade de Chebyshev estabelece que
\[ P(|X-\mu|<k\sigma)\geq 1-\frac{1}{k^2} \]
para \(k>1\).
Consequentemente, para obter uma cobertura mínima de \(1-\alpha\), escolhe-se
\[ k=\frac{1}{\sqrt{\alpha}} \]
Assim,
\[ P\left( \mu-\frac{\sigma}{\sqrt{\alpha}} < X < \mu+\frac{\sigma}{\sqrt{\alpha}} \right) \geq 1-\alpha \]
Para uma cobertura mínima de 95%,
\[ \alpha=0.05 \]
e
\[ k=\frac{1}{\sqrt{0.05}} =\sqrt{20} \approx 4.472 \]
Portanto,
\[ I_C= \left[ \mu-4.472\sigma,\; \mu+4.472\sigma \right] \]
A principal vantagem da desigualdade de Chebyshev é não exigir uma forma específica para a distribuição de \(X\). A distribuição pode ser normal, assimétrica ou assumir outras formas, desde que média e variância sejam finitas.
Entretanto, a cobertura fornecida por Chebyshev é uma garantia mínima:
\[ P(X\in I_C)\geq0.95 \]
e não necessariamente
\[ P(X\in I_C)=0.95 \]
Por essa razão, o intervalo de Chebyshev tende a ser bastante conservador.
Intervalo de referência central
No laboratório clínico, um intervalo de referência é construído a partir da distribuição dos resultados em uma população de referência adequadamente definida.
Uma construção não paramétrica natural para um intervalo de referência de 95% é
\[ RI_{95\%} = \left[ Q_{0.025},Q_{0.975} \right] \]
em que \(Q_p\) representa o quantil de ordem \(p\).
Esse intervalo contém os 95% centrais da distribuição, deixando
\[ P(X<Q_{0.025})=0.025 \]
e
\[ P(X>Q_{0.975})=0.025 \]
Assim,
\[ P(Q_{0.025}\leq X\leq Q_{0.975})=0.95 \]
A definição não exige normalidade. É, portanto, especialmente conveniente para resultados laboratoriais cujas distribuições podem apresentar assimetria.
Quando
\[ X\sim N(\mu,\sigma^2) \]
tem-se aproximadamente
\[ Q_{0.025}=\mu-1.96\sigma \]
e
\[ Q_{0.975}=\mu+1.96\sigma \]
Logo,
\[ RI_{95\%} \approx \left[ \mu-1.96\sigma,\; \mu+1.96\sigma \right] \]
Nesse caso,
\[ P(X\in RI_{95\%})=0.95 \]
Isso contrasta com Chebyshev, que fornece apenas uma cobertura mínima de 95%.
Intervalo de maior densidade
Uma terceira possibilidade é definir uma região que contenha 95% da probabilidade e que privilegie os valores de maior densidade.
O Highest Density Interval de probabilidade \(1-\alpha\) pode ser caracterizado por
\[ HDI_{1-\alpha} = \{x:f(x)\geq c\} \]
em que \(f(x)\) é a função densidade de probabilidade e \(c\) é escolhido de modo que
\[ P(X\in HDI_{1-\alpha})=1-\alpha \]
Para uma distribuição contínua e unimodal, o HDI também pode ser interpretado como o intervalo de menor comprimento que satisfaz
\[ P(a\leq X\leq b)=1-\alpha \]
Ou seja,
\[ HDI_{1-\alpha} = \operatorname*{arg\,min}_{a,b}(b-a) \]
sujeito a
\[ P(a\leq X\leq b)=1-\alpha. \]
Portanto, o intervalo de referência central e o HDI utilizam critérios diferentes.
O intervalo central exige probabilidades iguais nas duas caudas:
\[ P(X<LI)=P(X>LS)=\frac{\alpha}{2} \]
O HDI, por sua vez, procura concentrar a probabilidade na região de maior densidade, sem exigir probabilidades iguais nas duas caudas.
Distribuições simétricas
Para uma distribuição normal,
\[ X\sim N(\mu,\sigma^2) \]
a simetria implica que o intervalo central de 95% também corresponde ao intervalo de menor comprimento contendo 95% da probabilidade.
Assim,
\[ HDI_{95\%} = RI_{95\%} = \left[ \mu-1.96\sigma,\; \mu+1.96\sigma \right] \]
Portanto, em distribuições normais, a distinção entre os dois critérios não produz diferença prática nos limites populacionais.
Isso não ocorre necessariamente em distribuições assimétricas.
Distribuições assimétricas
Considere uma variável com distribuição Gamma. Como a distribuição apresenta assimetria positiva, o intervalo central de 95% é definido por
\[ RI_{95\%} = \left[ Q_{0.025},Q_{0.975} \right] \]
Esse intervalo necessariamente mantém 2.5% da probabilidade em cada cauda.
O HDI, entretanto, pode deslocar seus limites para a esquerda, incluindo regiões de maior densidade e excluindo uma proporção maior da longa cauda direita.
Consequentemente,
\[ HDI_{95\%} \neq RI_{95\%} \]
e, tipicamente,
\[ \operatorname{largura}(HDI_{95\%}) < \operatorname{largura}(RI_{95\%}) \]
Por exemplo, para
\[ X\sim\operatorname{Gamma}(2,2) \]
obtêm-se aproximadamente
\[ RI_{95\%} = \left[ 0.484,\;11.143 \right] \]
com largura
\[ 11.143-0.484 = 10.659 \]
O HDI correspondente é aproximadamente
\[ HDI_{95\%} = \left[ 0.085,\;9.530 \right] \]
com largura
\[ 9.530-0.085 = 9.446 \]
Ambos contêm 95% da probabilidade, mas são determinados por critérios distintos.
Distribuição uniforme
A distribuição uniforme apresenta uma situação particular.
Se
\[ X\sim U(a,b) \]
a densidade é constante em todo o suporte:
\[ f(x)=\frac{1}{b-a}, \qquad a\leq x\leq b \]
O intervalo central de 95% é perfeitamente definido:
\[ RI_{95\%} = \left[ Q_{0.025},Q_{0.975} \right] \]
Para
\[ X\sim U(0,10) \]
tem-se
\[ RI_{95\%} = \left[ 0.25,\;9.75 \right] \]
Sua largura é
\[ 9.75-0.25=9.5 \]
Entretanto, não existe um único HDI de 95%, pois todos os pontos do suporte possuem a mesma densidade.
Qualquer intervalo
\[ [a^\ast,b^\ast]\subset[0,10] \]
com comprimento
\[ b^\ast-a^\ast=9.5 \]
contém 95% da probabilidade e possui a mesma densidade.
Portanto, para uma distribuição uniforme, o HDI não é único.
Comparação entre os três critérios
Os três intervalos respondem a questões diferentes.
O intervalo de Chebyshev procura fornecer uma garantia mínima:
\[ P(X\in I_C)\geq1-\alpha \]
O intervalo de referência central procura conter uma proporção determinada da população, distribuindo igualmente a probabilidade excluída entre as duas caudas:
\[ P(X<LI)=P(X>LS)=\frac{\alpha}{2} \]
O HDI procura a região de maior densidade contendo determinada probabilidade:
\[ P(X\in HDI)=1-\alpha \]
Para 95%, pode-se resumir:
| Intervalo | Critério |
|---|---|
| Chebyshev | cobertura mínima de 95% |
| Reference interval | 95% central da distribuição |
| HDI | região de maior densidade contendo 95% |
Uma diferença fundamental é que
\[ P(X\in I_C)\geq0.95 \]
para Chebyshev, enquanto, idealmente,
\[ P(X\in RI_{95\%})=0.95 \]
e
\[ P(X\in HDI_{95\%})=0.95 \]
Implicações para exames clínicos
Para resultados de exames clínicos, o intervalo de referência central apresenta uma interpretação particularmente simples:
\[ RI_{95\%} = \left[ Q_{0.025},Q_{0.975} \right] \]
Ele descreve os 95% centrais dos valores observados na população de referência.
Isso significa que, mesmo em uma população de referência adequadamente selecionada, aproximadamente 5% dos indivíduos apresentarão resultados fora do intervalo:
\[ P(X\notin RI_{95\%})=0.05 \]
Consequentemente, estar fora do intervalo de referência não é equivalente a apresentar doença. Da mesma forma, estar dentro do intervalo não implica necessariamente ausência de doença.
O intervalo de referência é uma propriedade da distribuição dos resultados na população de referência, e não um limite diagnóstico.
O HDI fornece uma descrição alternativa da região em que os valores estão mais concentrados. Essa distinção torna-se relevante principalmente em distribuições assimétricas.
Chebyshev, por sua vez, fornece uma garantia probabilística muito mais geral, mas normalmente produz intervalos excessivamente amplos para utilização como intervalo de referência clínico.
# ============================================================
# Comparação de intervalos probabilísticos
# Distribuições: Uniforme, Normal e Gamma assimétrica
#
# Intervalos comparados:
# 1. Chebyshev
# 2. Intervalo de referência central de 95%
# 3. HDI empírico de 95%
# ============================================================
set.seed(123)
# ------------------------------------------------------------
# Parâmetros gerais
# ------------------------------------------------------------
n <- 1e6
prob <- 0.95
alpha <- 1 - prob
# ------------------------------------------------------------
# 1. Intervalo de Chebyshev
# ------------------------------------------------------------
intervalo_chebyshev <- function(x, prob = 0.95) {
alpha <- 1 - prob
k <- 1 / sqrt(alpha)
LI <- mean(x) - k * sd(x)
LS <- mean(x) + k * sd(x)
c(
LI = LI,
LS = LS
)
}
# ------------------------------------------------------------
# 2. Intervalo de referência central
#
# Para 95%:
#
# [Q_0.025 ; Q_0.975]
# ------------------------------------------------------------
intervalo_referencia <- function(x, prob = 0.95) {
alpha <- 1 - prob
q <- stats::quantile(
x,
probs = c(
alpha / 2,
1 - alpha / 2
),
names = FALSE,
type = 7
)
c(
LI = q[1],
LS = q[2]
)
}
# ------------------------------------------------------------
# 3. HDI empírico
#
# Procura o intervalo de menor comprimento que contém
# aproximadamente a proporção "prob" das observações.
#
# Adequado para distribuições unimodais.
# ------------------------------------------------------------
intervalo_hdi <- function(x, prob = 0.95) {
x <- sort(x)
n <- length(x)
m <- floor(prob * n)
largura <- x[(m + 1):n] -
x[1:(n - m)]
i <- which.min(largura)
c(
LI = x[i],
LS = x[i + m]
)
}
# ------------------------------------------------------------
# Função para calcular cobertura empírica
# ------------------------------------------------------------
cobertura <- function(x, intervalo) {
mean(
x >= intervalo["LI"] &
x <= intervalo["LS"]
)
}
# ------------------------------------------------------------
# Função para comparar os três intervalos
# ------------------------------------------------------------
comparar_intervalos <- function(x,
distribuicao,
prob = 0.95) {
cheb <- intervalo_chebyshev(
x,
prob
)
ref <- intervalo_referencia(
x,
prob
)
hdi <- intervalo_hdi(
x,
prob
)
data.frame(
Distribuicao =
rep(distribuicao, 3),
Intervalo =
c(
"Chebyshev",
"Reference interval",
"HDI"
),
LI =
c(
cheb["LI"],
ref["LI"],
hdi["LI"]
),
LS =
c(
cheb["LS"],
ref["LS"],
hdi["LS"]
),
Largura =
c(
cheb["LS"] - cheb["LI"],
ref["LS"] - ref["LI"],
hdi["LS"] - hdi["LI"]
),
Cobertura =
c(
cobertura(x, cheb),
cobertura(x, ref),
cobertura(x, hdi)
),
row.names = NULL,
check.names = FALSE
)
}
# ============================================================
# GERAÇÃO DAS DISTRIBUIÇÕES
# ============================================================
# ------------------------------------------------------------
# Distribuição uniforme
#
# X ~ U(0,10)
# ------------------------------------------------------------
x_uniforme <- stats::runif(
n,
min = 0,
max = 10
)
# ------------------------------------------------------------
# Distribuição normal
#
# X ~ N(5, 2^2)
# ------------------------------------------------------------
x_normal <- stats::rnorm(
n,
mean = 5,
sd = 2
)
# ------------------------------------------------------------
# Distribuição Gamma assimétrica
#
# X ~ Gamma(shape = 2, scale = 2)
#
# média = shape * scale = 4
# variância = shape * scale^2 = 8
# ------------------------------------------------------------
x_gamma <- stats::rgamma(
n,
shape = 2,
scale = 2
)
# ============================================================
# CÁLCULO DOS INTERVALOS
# ============================================================
resultado_uniforme <-
comparar_intervalos(
x_uniforme,
"Uniforme",
prob
)
resultado_normal <-
comparar_intervalos(
x_normal,
"Normal",
prob
)
resultado_gamma <-
comparar_intervalos(
x_gamma,
"Gamma",
prob
)
resultado <- rbind(
resultado_uniforme,
resultado_normal,
resultado_gamma
)
# ------------------------------------------------------------
# Resultados
# ------------------------------------------------------------
resultado[, c(
"LI",
"LS",
"Largura",
"Cobertura"
)] <-
round(
resultado[, c(
"LI",
"LS",
"Largura",
"Cobertura"
)],
3
)
print(
resultado,
row.names = FALSE
) Distribuicao Intervalo LI LS Largura Cobertura
Uniforme Chebyshev -7.921 17.911 25.832 1.000
Uniforme Reference interval 0.248 9.749 9.502 0.950
Uniforme HDI 0.126 9.623 9.497 0.950
Normal Chebyshev -3.937 13.929 17.866 1.000
Normal Reference interval 1.081 8.908 7.827 0.950
Normal HDI 1.092 8.919 7.827 0.950
Gamma Chebyshev -8.639 16.635 25.275 0.998
Gamma Reference interval 0.485 11.131 10.646 0.950
Gamma HDI 0.084 9.530 9.445 0.950# ============================================================
# GRÁFICOS
# ============================================================
# ------------------------------------------------------------
# Função gráfica
# ------------------------------------------------------------
plot_intervalos <- function(x,
titulo,
prob = 0.95) {
cheb <-
intervalo_chebyshev(
x,
prob
)
ref <-
intervalo_referencia(
x,
prob
)
hdi <-
intervalo_hdi(
x,
prob
)
d <- stats::density(x)
plot(
d$x,
d$y,
type = "l",
lwd = 2,
main = titulo,
xlab = "Valor",
ylab = "Densidade"
)
ymax <- max(d$y)
# Chebyshev
segments(
x0 = cheb["LI"],
y0 = ymax * 0.90,
x1 = cheb["LS"],
y1 = ymax * 0.90,
lwd = 5
)
# Reference interval
segments(
x0 = ref["LI"],
y0 = ymax * 0.80,
x1 = ref["LS"],
y1 = ymax * 0.80,
lwd = 5,
lty = 2
)
# HDI
segments(
x0 = hdi["LI"],
y0 = ymax * 0.70,
x1 = hdi["LS"],
y1 = ymax * 0.70,
lwd = 5,
lty = 3
)
legend(
"bottom",
legend = c(
"Chebyshev",
"Reference interval",
"HDI"
),
lwd = c(5, 5, 5),
lty = c(1, 2, 3),
bty = "n"
)
}
# ------------------------------------------------------------
# Três gráficos
# ------------------------------------------------------------
par(
mfrow = c(1, 3),
mar = c(4, 4, 3, 1)
)
plot_intervalos(
x_uniforme,
"Uniforme U(0,10)",
prob
)
plot_intervalos(
x_normal,
"Normal N(5, 2²)",
prob
)
plot_intervalos(
x_gamma,
"Gamma(2,2)",
prob
)
par(
mfrow = c(1, 1)
)
# ============================================================
# VALORES TEÓRICOS PARA COMPARAÇÃO
# ============================================================
# ------------------------------------------------------------
# NORMAL
# ------------------------------------------------------------
mu_normal <- 5
sigma_normal <- 2
normal_chebyshev <-
c(
LI =
mu_normal -
sigma_normal / sqrt(alpha),
LS =
mu_normal +
sigma_normal / sqrt(alpha)
)
normal_reference <-
stats::qnorm(
c(
alpha / 2,
1 - alpha / 2
),
mean = mu_normal,
sd = sigma_normal
)
# Para distribuição normal:
# HDI = intervalo central
normal_hdi <-
normal_reference
cat(
"\n\nNORMAL - valores teóricos\n"
)
NORMAL - valores teóricos LI LS
Chebyshev -3.944272 13.944272
Reference 1.080072 8.919928
HDI 1.080072 8.919928# ------------------------------------------------------------
# UNIFORME
# ------------------------------------------------------------
a <- 0
b <- 10
mu_uniforme <-
(a + b) / 2
sigma_uniforme <-
(b - a) / sqrt(12)
uniforme_chebyshev <-
c(
LI =
mu_uniforme -
sigma_uniforme /
sqrt(alpha),
LS =
mu_uniforme +
sigma_uniforme /
sqrt(alpha)
)
uniforme_reference <-
stats::qunif(
c(
alpha / 2,
1 - alpha / 2
),
min = a,
max = b
)
cat(
"\n\nUNIFORME - valores teóricos\n"
)
UNIFORME - valores teóricos LI LS
Chebyshev -7.909944 17.90994
Reference 0.250000 9.75000
Na distribuição uniforme o HDI não é único.cat(
"Qualquer intervalo interno ao suporte com comprimento",
prob * (b - a),
"possui probabilidade",
prob,
".\n"
)Qualquer intervalo interno ao suporte com comprimento 9.5 possui probabilidade 0.95 .# ------------------------------------------------------------
# GAMMA
# ------------------------------------------------------------
shape <- 2
scale <- 2
mu_gamma <-
shape * scale
sigma_gamma <-
sqrt(
shape * scale^2
)
gamma_chebyshev <-
c(
LI =
mu_gamma -
sigma_gamma /
sqrt(alpha),
LS =
mu_gamma +
sigma_gamma /
sqrt(alpha)
)
gamma_reference <-
stats::qgamma(
c(
alpha / 2,
1 - alpha / 2
),
shape = shape,
scale = scale
)
cat(
"\n\nGAMMA - valores teóricos\n"
)
GAMMA - valores teóricos LI LS
Chebyshev -8.6491106 16.64911
Reference 0.4844186 11.14329# ============================================================
# RESUMO CONCEITUAL
# ============================================================
cat(
"\n\n",
"Chebyshev:\n",
"P(X no intervalo) >= 0.95\n\n",
"Reference interval:\n",
"[Q0.025 ; Q0.975]\n\n",
"HDI:\n",
"intervalo de menor comprimento contendo aproximadamente 95%\n"
)
Chebyshev:
P(X no intervalo) >= 0.95
Reference interval:
[Q0.025 ; Q0.975]
HDI:
intervalo de menor comprimento contendo aproximadamente 95%Referências
- Albers, CJ et al. (2018) Credible confidence: A pragmatic view on the frequentist vs Bayesian debate. Collabra: Psychology, 4(1): 31. DOI: https://doi.org/10.1525/collabra.149; Código R: https://osf.io/dgfht/.
- Auvin S, Irwin J, Abi-Aad, Battersby A (2018) The Problem of Rarity: Estimation of Prevalence in Rare Disease. Value in Health 21(5): 501-07. https://doi.org/10.1016/j.jval.2018.03.002
- Banjanovic, ES & Osborne, JW (2016) Confidence intervals for effect sizes: Applying bootstrap resampling. Practical Assessment, Research & Evaluation 21(5).
- Chihara LM & Hesterberg TC (2019) Mathematical statistics with resampling and R. 2nd ed. NJ: Wiley.
- Clopper C & Pearson E (1934). The Use of Confidence or Fiducial Limits Illustrated in the Case of the Binomial. Biometrika 26(4): 404-13. doi:10.2307/2331986
- Cumming, G & Finch, S (2005) Inference by eye: Confidence intervals and how to read pictures of data. American Psychologist 60(2): 170-80.
- DANCEY, CP & REIDY, J (2019) Estatística sem Matemática para Psicologia. 7ª ed. Porto Alegre: Penso.
- Doyle, K & Bunch, DR (2023) Reference intervals: Past, present, and future. Critical Reviews in Clinical Laboratory Sciences 60(6): 466-482. doi:10.1080/10408363.2023.2196746.
- Gilbertson D, Ebben J, Foley R, Weinhandl E, Bradbury B, Collins A (2008) Hemoglobin Level Variability: Associations with Mortality. Clinical journal of the American Society of Nephrology CJASN 3: 133-8. 10.2215/CJN.01610407.
- Horn PS, Pesce AJ. (2003) Reference intervals: an update. Clinica Chimica Acta. 334(1–2): 5-23. DOI: 10.1016/S0009-8981(03)00133-5.
- Krzywinski, M & Altman, N (2013) Importance of being uncertain. Nature Methods 10(9): 809-10.
- Puth, M et al. (2015) On the variety of methods for calculating confidence intervals by bootstrapping. Journal of Animal Ecology 84: 892-7.
- Rahme E & Joseph L (2002) Estimating the prevalence of a rare disease: adjusted maximum likelihood. The Statistician 47(1):149-58. https://doi.org/10.1111/1467-9884.00120
- Ranstam, J (2009) Sampling uncertainty in medical research. Osteoarthritis and Cartilage 17:1416-9.
- Roach, KE (2001) A clinician’s guide to specification and sampling. Journal of Orthopaedic & Sports Physical Therapy 31(12): 753-8.
- Robert D (2020) Five Confidence Intervals for Proportions That You Should Know About. Disponível em https://towardsdatascience.com/five-confidence-intervals-for-proportions-that-you-should-know-about-7ff5484c024f
- Siqueira JO, Vieira JE, Silveira PSP (2017) Os testes estatísticos. In: Luiz Marciano Cangiani; Maria José Carvalho Carmona; Marcelo Luiz Abramides Torres; Carlos Othon Batos; David Ferez; Enis Donizetti Silva; Leonardo Duarte; Maria Angela Tardelli. (Org.). Tratado de Anestesiologia. 8ed. Rio de Janeiro: Atheneu 2: 3863-78.
- Tendeiro, JN & Kiers, HAL (2019) A review of issues about null hypothesis Bayesian testing. Psychological Methods 24(6): 774–95. https://doi.org/10.1037/met0000221.
- van Belle G (2008) Statistical rules of thumb. 2nd ed. NJ: Wiley.
- Wallis SA (2013) Binomial confidence intervals and contingency tests: mathematical fundamentals and the evaluation of alternative methods. Journal of Quantitative Linguistics 20(3): 178-208.
- Weissgerber TL, Milic NM, Winham SJ, Garovic VD (2015) Beyond bar and line graphs: Time for a new data presentation paradigm. PLoS Biol 13(4): e1002128. doi:10.1371/journal.
- White, C (1953) Sampling in medical research. BMJ 12: 1284-8.
- Wilson EB (1927) Probable inference, the law of succession, and statistical inference. Journal of the American Statistical Association 22(158): 209–12.
- Wonnacott, TH & Wonnacott, RJ (1969) Introductory Statistics. NJ: Wiley. * Wonnacott, TH & Wonnacott, RJ (1981) Estatística aplicada à Economia e à Administração. RJ: LTC.
- Wonnacott, TH & Wonnacott, RJ (1990) Introductory Statistics for Business and Economics. NJ: Wiley.
- Zygmont, CS (2023) Managing the assumption of normality within the general linear model with small samples: Guidelines for researchers regarding if, when and how. The Quantitative Methods for Psychology 19(4), 302–32. https://doi.org/10.20982/tqmp.19.4.p302