Fundamentos de Inferência Estatística em R
MSP4112: Estatística Aplicada à Medicina
15 agosto 2025 12:43h
Bastão de Asclépio & Distribuição Normal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))options(warn=-1)
suppressMessages(library(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 do limite central
- 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
tibble [549 × 13] (S3: tbl_df/tbl/data.frame)
$ ID : Factor w/ 549 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Ano : num [1:549] 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 [1:549] 68 82 79 49 52 73 73 60 75 75 ...
$ Estatura : num [1:549] 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 [1:549] 22.7 26.8 26.7 19.1 19.3 ...
NULL[1] 232[1] 317 pos variable label col_type missing
1 1 Ano <NA> dbl 0
2 2 Turma <NA> fct 0
3 3 Sexo <NA> fct 0
4 4 Mao <NA> fct 42
5 5 TipoSang <NA> fct 11
6 6 ABO <NA> fct 11
7 7 AtivFisica <NA> fct 0
8 8 Sedentarismo <NA> fct 0
9 9 MCT <NA> dbl 7
10 10 Estatura <NA> dbl 7
11 11 Rh <NA> fct 11
12 12 IMC <NA> dbl 8
levels
1 NULL
2 A, B
3 F, M
4 A, C, D
5 A+, O+, B+, AB+, A-, O-, B-, AB-
6 A, O, B, AB
7 sempre_inativo, atualmente_inativo, baixa_intensidade, media_intensidade, alta_intensidade
8 N, S
9 NULL
10 NULL
11 -, +
12 NULL
value_labels
1 NULL
2 NULL
3 NULL
4 NULL
5 NULL
6 NULL
7 NULL
8 NULL
9 NULL
10 NULL
11 NULL
12 NULLtabela <- 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 | 3 | B | F | C | NA | NA | alta_intensidade | N | 54 | 157 | NA | 21.90758 |
| 545 | 3 | B | F | D | B+ | B | alta_intensidade | N | 47 | 157 |
|
19.06771 |
| 546 | 3 | B | F | D | B+ | B | sempre_inativo | S | 48 | 162 |
|
18.28989 |
| 547 | 3 | B | M | A | A+ | A | alta_intensidade | N | 70 | 175 |
|
22.85714 |
| 548 | 3 | B | M | D | O- | O | atualmente_inativo | S | 67 | 171 |
|
22.91303 |
| 549 | 3 | B | F | D | A+ | A | media_intensidade | N | 45 | 160 |
|
17.57812 |
ID Ano Turma Sexo Mao TipoSang
1 : 1 Min. :1.000 A:285 F:232 A : 28 A+ :171
2 : 1 1st Qu.:1.000 B:264 M:317 C : 37 O+ :132
3 : 1 Median :2.000 D :442 B+ : 65
4 : 1 Mean :2.095 NA's: 42 A- : 49
5 : 1 3rd Qu.:3.000 O- : 45
6 : 1 Max. :3.000 (Other): 76
(Other):543 NA's : 11
ABO AtivFisica Sedentarismo MCT
A :220 sempre_inativo : 71 N:316 Min. : 41.00
O :177 atualmente_inativo:152 S:233 1st Qu.: 56.00
B : 99 baixa_intensidade : 59 Median : 64.00
AB : 42 media_intensidade : 99 Mean : 65.67
NA's: 11 alta_intensidade :168 3rd Qu.: 73.00
Max. :125.00
NA's :7
Estatura Rh IMC
Min. :153.0 - :130 Min. :14.53
1st Qu.:164.0 + :408 1st Qu.:20.03
Median :172.0 NA's: 11 Median :21.79
Mean :170.9 Mean :22.32
3rd Qu.:177.0 3rd Qu.:24.07
Max. :195.0 Max. :40.35
NA's :7 NA's :8
Usando os dados de `Biometria_FMUSP.rds`, executamos [`fi_IntervaloPredicao.R`](fi_IntervaloPredicao.R){target="_blank"}:
<img src="FundamentosInferencia_files/figure-html/unnamed-chunk-8-1.png" width="672" style="display: block; margin: auto;" /><img src="FundamentosInferencia_files/figure-html/unnamed-chunk-8-2.png" width="672" style="display: block; margin: auto;" />
__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.
<img src="image/me.png" width="80%" style="display: block; margin: auto;" />
* [Margem de erro: Wikipédia](https://pt.wikipedia.org/wiki/Margem_de_erro){target="_blank"}
* [Intervalo de confiança: Wikipédia](https://pt.wikipedia.org/wiki/Intervalo_de_confian%C3%A7a){target="_blank"}
* [Vídeo: Confidence Interval Creation: YouTube](https://www.youtube.com/watch?v=vNdxtdEjBQk){target="_blank"}
* [Confidence Interval Creation: Demo](https://wise.cgu.edu/portfolio2/demo-confidence-interval-creation/){target="_blank"}
__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`](fi_ToDoMargemErro_Bonf.R){target="_blank"}, resulta:
``` output
---- 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 brutos
alpha <- alfa <- 0.05
alfa.Bonf <- alpha/length(unique(Dados$Sexo)) # Correção de Bonferroni
# Com dados brutos
CI_Fmean_t <- confintr::ci_mean(x = Dados.F$MCT,
probs = c(alfa.Bonf/2, 1-alfa.Bonf/2),
type = "t")
CI_Fmean_t
Two-sided 97.5% t confidence interval for the population mean
Sample estimate: 57.63043
Confidence interval:
1.25% 98.75%
56.28362 58.97725 CI_Mmean_t <- confintr::ci_mean(x = Dados.M$MCT,
probs = c(alfa.Bonf/2, 1-alfa.Bonf/2),
type = "t")
CI_Mmean_t
Two-sided 97.5% t confidence interval for the population mean
Sample estimate: 71.58974
Confidence interval:
1.25% 98.75%
70.04848 73.13101 mean lwr.ci upr.ci
57.63043 56.28362 58.97725 mean lwr.ci upr.ci
71.58974 70.04848 73.13101 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 brutos
# Sem dados brutos
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/2, 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
X11 1 F 1 230 57.63 9.05 56 56.75 7.41 41 120 79 2.43
X12 2 M 1 312 71.59 12.09 72 70.88 11.12 43 125 82 0.90
kurtosis se
X11 12.04 0.60
X12 2.37 0.68
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.832448
P-value: 4.940064e-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.9541195
P-value: 2.613898e-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
variable <- Dados.M$MCT
print(m <- DescTools::MeanCI(variable, na.rm=TRUE)) mean lwr.ci upr.ci
71.58974 70.24329 72.93620 var lwr.ci upr.ci
12.08723 11.20739 13.11815 invisible(out <- conf::crplot(dataset=sort(variable),
distn="norm",
alpha=alfa,
pts=FALSE,
animate=FALSE,
info=TRUE,
silent=TRUE,
# animate=TRUE,
main="Normal\nMCT Masculino"))
conf::crplot(dataset=sort(variable),
distn="norm",
alpha=alfa,
pts=FALSE,
origin=TRUE,
main="Normal\nMCT Masculino")
[1] "95% confidence region plot complete; made using 120 boundary points."# cat("\nMean\n")
# print(identical(as.numeric(out$muhat), as.numeric(m[1])))
# print(all.equal(as.numeric(out$muhat), as.numeric(m[1])))
# cat("\nSD\n")
# print(identical(as.numeric(out$sigmahat), as.numeric(s[1])))
# print(all.equal(as.numeric(out$sigmahat), as.numeric(s[1])))
# cat("\nMean CI95%: Lower Limit\n")
# print(all.equal(min(as.numeric(out$mu)), as.numeric(m[2])))
# cat("\nMean CI95%: Upper Limit\n")
# print(all.equal(max(as.numeric(out$mu)), as.numeric(m[3])))
# cat("\nSD CI95%: Lower Limit\n")
# print(all.equal(min(as.numeric(out$sigma)), as.numeric(s[2])))
# cat("\nSD CI95%: Upper Limit\n")
# print(all.equal(max(as.numeric(out$sigma)), as.numeric(s[3])))
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\)).
alpha <- alfa <- 0.05
alfa.Bonf <- alpha/length(unique(Dados$Sexo)) # Correção de Bonferroni
# Com dados brutos
R <- 1e4
seed <- 123
CI_Fmean_t <- confintr::ci_mean(x = Dados.F$MCT,
probs = c(alfa.Bonf/2, 1-alfa.Bonf/2),
type = "t")
CI_Fmean_t
Two-sided 97.5% t confidence interval for the population mean
Sample estimate: 57.63043
Confidence interval:
1.25% 98.75%
56.28362 58.97725 CI_Fmean_B <- confintr::ci_mean(x = Dados.F$MCT,
probs = c(alfa.Bonf/2, 1-alfa.Bonf/2),
type = "bootstrap",
R = R,
seed = seed)
CI_Fmean_B
Two-sided 97.5% bootstrap confidence interval for the population mean
based on 10000 bootstrap replications and the student method
Sample estimate: 57.63043
Confidence interval:
1.25% 98.75%
56.42961 59.19607 CI_Mmean_t <- confintr::ci_mean(x = Dados.M$MCT,
probs = c(alfa.Bonf/2, 1-alfa.Bonf/2),
type = "t")
CI_Mmean_t
Two-sided 97.5% t confidence interval for the population mean
Sample estimate: 71.58974
Confidence interval:
1.25% 98.75%
70.04848 73.13101 CI_Mmean_B <- confintr::ci_mean(x = Dados.M$MCT,
probs = c(alfa.Bonf/2, 1-alfa.Bonf/2),
type = "bootstrap",
R = R,
seed = seed)
CI_Mmean_B
Two-sided 97.5% bootstrap confidence interval for the population mean
based on 10000 bootstrap replications and the student method
Sample estimate: 71.58974
Confidence interval:
1.25% 98.75%
70.09580 73.21973 mean lwr.ci upr.ci
57.63043 56.28362 58.97725 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)) mean lwr.ci upr.ci
57.63043 56.25223 58.90435 mean lwr.ci upr.ci
71.58974 70.04848 73.13101 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)) mean lwr.ci upr.ci
71.58974 70.04167 73.12496 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 Cramer populacionalci_f_ncp— parâmetro de não centralidade da distribuição Fci_IQR— intervalo interquartílicoci_kurtosis— curtoseci_mad— desvio absoluto mediano (MAD)ci_mean— média populacionalci_mean_diff— diferença de médias populacionaisci_median— mediana populacionalci_median_diff— diferença de medianas populacionais entre duas amostrasci_oddsratio— odds ratioci_proportion— proporção populacionalci_quantile— quantil populacionalci_quantile_diff— diferença de quantis populacionais entre duas amostrasci_rsquared— R² populacionalci_sd— desvio-padrão populacionalci_skewness— assimetriaci_var— variância populacional
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)
DescTools::BootCI(Dados.F$MCT,
FUN=mean,
conf.level=1-alpha,
sides="two.sided",
na.rm=TRUE,
R=1e4) mean lwr.ci upr.ci
57.63043 56.45361 58.80224 set.seed(123)
DescTools::BootCI(Dados.F$MCT,
FUN=DescTools::Var,
conf.level=1-alpha,
sides="two.sided",
na.rm=TRUE,
R=1e4)DescTools::Var lwr.ci upr.ci
81.95453 42.13299 122.55889 spearman <- function(x,y) cor(x, y, method="spearman", use="p")
set.seed(123)
DescTools::BootCI(Dados.F$MCT,
Dados.F$Estatura,
FUN=spearman,
conf.level=1-alpha,
sides="two.sided",
R=1e4) spearman lwr.ci upr.ci
0.4917515 0.3865515 0.6003334 Teorema Central do Limite
O Teorema Central do Limite (TCL) assume:
- observações independentes,
- variável intervalar,
- amostra com pelo menos 30 observações.
e enuncia que:
“[…] the Central Limit Theorem […] states that as the size of the samples we select increases, the nearer to the population mean will be the mean of these sample means and the closer to normal will be the distribution of the sample means.”
traduzido para:
“[…] 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.”
Dancey & Reidy, 2019, p. 108
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}\).”
Ou, 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.
A frase define seu cálculo:
“… erro-padrão [é igual ao] desvio-padrão populacional dividido pela raiz quadrada do tamanho da amostra.”
Dancey & Reidy, 2019
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\)?
De onde veio esta 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?
Os autores tentam ilustrar o teorema com figuras:
Wonnacott & Wonnacott, 1990
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.
Fonte: 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.984217
EP = s/sqrt(n) = 0.1609152
IC95 = 11.945 +- 1.984217(0.975,99) * 1.609152/sqrt(100)
IC95 = [11.62571,12.26429]|
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:
v <- ""
v <- paste0(v,"boot.medias (",length(boot.medias)," elementos): ")
v <- paste0(v,paste(round(boot.medias[1:50],2),collapse=" "))
v <- paste0(v," ... ")
v <- paste0(v,paste(round(boot.medias[(B-50):B],2),collapse=" "))
v <- paste0(v,"\n\n")
v <- paste0(v,"boot.sds (",length(boot.sds)," elementos): ")
v <- paste0(v,paste(round(boot.sds[1:50],2),collapse=" "))
v <- paste0(v," ... ")
v <- paste0(v,paste(round(boot.sds[(B-50):B],2),collapse=" "))
cat(v)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
v <- ""
bm <- mean(boot.medias)
bs <- mean(boot.sds)
v <- paste0(v,"\nAmostras:",B," com n = ",n,"\n")
v <- paste0(v,"\tmedia das medias amostrais: ",round(bm,4),"\n")
v <- paste0(v,"\tmedia dos desvios padrao das medias amostrais: ",round(bs,4),"\n")
cat(v)
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:
v <- ""
ep <- sd(boot.medias)
v <- paste0(v,"\tdesvio padrao das medias amostrais: ",round(ep,4),"\n")
cat(v) desvio padrao das medias amostrais: 0.1593 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 de Wonnacott & Wonnacott (1990) 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
round(DescTools::MeanCI(Dados.F$MCT, conf.level=1-alfa.Bonf, na.rm=TRUE),1) mean lwr.ci upr.ci
57.6 56.3 59.0 mean lwr.ci upr.ci
71.6 70.0 73.1 [1] 56.3 59.0[1] 70.0 73.1# IC95% por DescTools::MeanCI com desvio-padrao conhecido
round(DescTools::MeanCI(Dados.F$MCT, sd=10, conf.level=1-alfa.Bonf, na.rm=TRUE), 1) mean lwr.ci upr.ci
57.6 56.2 59.1 mean lwr.ci upr.ci
71.6 69.8 73.4 # 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
tapply(Dados$MCT,
Dados$Sexo,
FUN=DescTools::MeanCI,
conf.level=1-alfa.Bonf,
na.rm=TRUE)$F
mean lwr.ci upr.ci
57.63043 56.28362 58.97725
$M
mean lwr.ci upr.ci
71.58974 70.04848 73.13101 # IC95% por bootstrapping para sexo feminino
round(DescTools::MeanCI(Dados.F$MCT,
method="boot",
type="perc",
R=1e4,
sides="two.sided",
conf.level=1-alfa.Bonf,
na.rm=TRUE), 1) mean lwr.ci upr.ci
57.6 56.3 59.0 # IC95% de MCT, Estatura e IMC para sexo feminino
round(do.call("rbind", lapply(Dados.F[, c("MCT","Estatura","IMC")],
FUN=DescTools::MeanCI,
conf.level=1-alfa.Bonf,
na.rm=TRUE)), 1) mean lwr.ci upr.ci
MCT 57.6 56.3 59.0
Estatura 164.3 163.4 165.3
IMC 21.2 20.8 21.6t(round(sapply(Dados.F[,c("MCT","Estatura","IMC")],
FUN=DescTools::MeanCI,
conf.level=1-alfa.Bonf,
na.rm=TRUE), 1)) mean lwr.ci upr.ci
MCT 57.6 56.3 59.0
Estatura 164.3 163.4 165.3
IMC 21.2 20.8 21.6Grá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 1 2 B M D A+ A baixa_intensidade N 68
2 2 2 B M D A+ A atualmente_inativo S 82
3 3 1 B M D A+ A media_intensidade N 79
4 4 1 B F D A+ A media_intensidade N 49
5 5 1 B F D O- O atualmente_inativo S 52
6 6 3 A M D A+ A alta_intensidade N 73
Estatura Rh IMC
1 173 + 22.72044
2 175 + 26.77551
3 172 + 26.70362
4 160 + 19.14062
5 164 - 19.33373
6 175 + 23.83673
Gráfico de intervalo de confiança de 95% de média: delineamento intraparticipantes
Obtido com fi_IntervalosDadosIntra.R:
subject condition value
1 1 pretest 59.4
2 2 pretest 46.4
3 3 pretest 46.0
4 4 pretest 49.0
5 5 pretest 32.5
6 6 pretest 45.2
7 7 pretest 60.3
8 8 pretest 54.3
9 9 pretest 45.4
10 10 pretest 38.9
11 1 posttest 64.5
12 2 posttest 52.4
13 3 posttest 49.7
14 4 posttest 48.7
15 5 posttest 37.4
16 6 posttest 49.5
17 7 posttest 59.9
18 8 posttest 54.1
19 9 posttest 49.6
20 10 posttest 48.5 condition N value valueNormed sd se ci
1 posttest 10 51.43 51.43 2.262361 0.7154214 1.920914
2 pretest 10 47.74 47.74 2.262361 0.7154214 1.920914
condition N value sd se ci
1 pretest 10 47.74 8.598992 2.719240 7.301189
2 posttest 10 51.43 7.253972 2.293907 6.159166
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 1 41 Round Monochromatic
2 2 57 Round Monochromatic
3 3 52 Round Monochromatic
4 4 49 Round Monochromatic
5 5 47 Round Monochromatic
6 6 37 Round Monochromatic Shape ColorScheme N Time TimeNormed sd se ci
1 Round Colored 12 43.58333 43.58333 1.212311 0.3499639 1.043198
2 Round Monochromatic 12 44.58333 44.58333 1.331438 0.3843531 1.145707
3 Square Colored 12 42.58333 42.58333 1.461630 0.4219364 1.257739
4 Square Monochromatic 12 43.58333 43.58333 1.261312 0.3641095 1.085364
HH::CIplot: simulação de intervalo de confiança
Uma outra implementação (fi_Manual.R), com 10 mil
reamostragens:
Média populacional=175
Desvio padrão populacional=8
n = 30 
Proporcao dos 10000 IC95(media pop)
que contem a media pop 175 = 0.9535Simulação de intervalo de confiança em Rpsychologist
DescTools::MedianCI: IC de mediana
alfa.Bonf <- alfa/length(unique(Dados$Sexo)) # Correção de Bonferroni
DescTools::MedianCI(Dados.F$MCT, conf.level=1-alfa.Bonf, na.rm=TRUE)median lwr.ci upr.ci
56 55 58
attr(,"conf.level")
[1] 0.9791935median lwr.ci upr.ci
72 70 73
attr(,"conf.level")
[1] 0.9798693DescTools::VarCI: IC de desvio-padrão
alfa.Bonf <- alfa/length(unique(Dados$Sexo)) # Correção de Bonferroni
sqrt(DescTools::VarCI(Dados.F$MCT, conf.level=1-alfa.Bonf, na.rm=TRUE)) var lwr.ci upr.ci
9.052874 8.192485 10.106666 var lwr.ci upr.ci
12.08723 11.08826 13.27594 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}} \]
onde \(\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}^{95\%}_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.
Dataframe dt_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
Prevalência
Em epidemiologia, a prevalência das doenças na população é uma pergunta 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):
Risco é 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:
\[ \text{IC}^{95\%}_{\text{Wald}} = \left[\hat{p} \pm 1.96 \text{EP}\right] \]
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}} = \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}}} \]
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, obtendo:
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} \]
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)
prmatrix(xci, rowlab="", quote=FALSE) est lwr.ci upr.ci
-0.2 -0.3357585 -0.05245293Observe 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:
$method
[1] "Pairwise comparison of proportions (Fisher)"
$data.name
[1] "errospuncao out of tentativaspuncao"
$p.value
R1 R2 R3
R2 1.00000000 NA NA
R3 0.02091641 0.0002035023 NA
Assistente 0.09381001 0.0093063228 1
$p.adjust.method
[1] "bonferroni"
attr(,"class")
[1] "pairwise.htest"
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:
$method
[1] "exact binomial tests"
$data.name
preferido
$p.adjust.method
[1] "bonf"
$p.value
Dream Brown Best Cocoa
Best Cocoa 0.2685932368 NA
Wonka 0.0000626297 0.04180163
attr(,"class")
[1] "pairwise.htest"
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 intervalar
- O parâmetro média 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).” https://en.wikipedia.org/wiki/Credible_interval
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.” https://en.wikipedia.org/wiki/Credible_interval
Intervalo de credibilidade de 95% da média populacional:
\[ \text{IC}_{\text{bayes}}^{95\%}(\mu)=\left[\dfrac{n_0\,\mu_0+n\,\bar{x}}{n_0+n} \pm t^{0.975}_{n-1}\,\dfrac{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
\[ \text{IC}^{95\%}_{\text{bayes}}(\mu)=\left[\dfrac{1.82\times55+230\times57.6}{1.82+230} \pm 1.97\,\dfrac{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\)
\[ \begin{align} \text{IC}^{95\%}_{\text{bayes}}(\mu)&=\left[56.4,58.8\right] \\ \text{IC}^{95\%}_{\text{freq}}(\mu)&=\left[56.5,58.8\right] \\ \text{IC}^{95\%}_{\text{boot}}(\mu)&=\left[56.6, 58.9\right] \end{align} \]
Se o tamanho da amostra é grande, i.e., \(n>30\), as amplitudes dos três intervalos de confiança de 95% 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 uma 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 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; Código 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
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.
- 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.
- 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.