Intervalo de Confiança em R
MAN5728: Tópicos Essenciais da
Bioestatística
MÓDULO 7 - ESTATÍSTICA APLICADA I:
FUNDAMENTOS E
TESTES ESTATÍSTICOS
23 março 2026 21:36h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(HH, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(Rmisc, warn.conflicts=FALSE))Material
- HTML de Rmarkdown em
RPubs - Prof. José Siqueira: ResearchGate
Conteúdo
- Métodos básicos de estatística inferencial
- Intervalo de confiança.
Ler planilha
Os dados da planilha Excel Biometria_FMUSP.xlsx foram
coletados pelos docentes do curso de Medicina da FMUSP dos estudantes do
segundo ano de uma mesma disciplina em três anos consecutivos.
As variáveis do arquivo são:
ID: idenficador do(a) estudanteAno: Ano da coleta dos dados: 1, 2, 3
Turma: A, BSexo: Feminino, Masculino
Mao: Destro, Canhoto, Ambidestro
TipoSang: A+, A-, …
ABO: A, B, AB, ORh: +, -AtivFisica: nível de atividade física
- `Sedentarismo`` : Não, Sim
MCT: massa corporal total (kg)Estatura: Estatura (cm)
Dados <- data.frame(readxl::read_excel(path="Biometria_FMUSP.xlsx",
sheet="dados",
na=c("NA","Na","nA","na")))
Dados$MCT[Dados$MCT==max(Dados$MCT,na.rm=TRUE)] <- NA
Dados$Estatura[Dados$Estatura==min(Dados$Estatura,na.rm=TRUE)] <- NA
Dados$ID <- factor(Dados$ID)
Dados$Ano <- factor(Dados$Ano)
Dados$Turma <- factor(Dados$Turma)
Dados$Sexo <- factor(Dados$Sexo)
Dados$Mao <- factor(Dados$Mao)
Dados$TipoSang <- factor(Dados$TipoSang)
Dados$ABO <- factor(Dados$ABO)
Dados$Rh <- factor(Dados$Rh)
Dados$AtivFisica <- factor(Dados$AtivFisica)
Dados$Sedentarismo <- factor(Dados$Sedentarismo)
Dados$IMC <- Dados$MCT/(Dados$Estatura/100)^2
Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")
alfa <- 0.05Amostragem e inferência estatística
- 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.
Como seria o mundo se tivesse 100 pessoas?
- LORD, A (2016) The World as 100 People: A Visual Guide to 7 Billion Humans. Smith Street Books.
- Vídeo: The world as 100 people
Amostragem estatística
Amostragem aleatória simples (SRS)
Amostragem sistemática
Amostragem estratificada
Amostragem por conglomerado
Vídeos sobre Amostragem
Intervalo de predição e Intervalo de confiança
Intervalo de predição da variável
\[ \Large \text{IP}^{\le75\%}(X) = \left[\bar{x} \pm 2\,s\right] \]
- \(X\): variável intervalar
- \(\bar{x}\): média amostral
- \(s\): desvio-padrão amostral
Descrição da amostra
Desigualdade de Chebyshev: Wikipedia
Independentemente da distribuição de MCT por condição, 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.
“Teste” de normalidade
Distribuição normal: Wikipedia
A distribuição normal contém aproximadamente 95% dos valores da amostra no intervalo simétrico centrado na média com amplitude 4 desvios-padrão.
# summarySE provides the standard deviation, standard error of the mean, and a (default 95%) confidence interval
s <- Rmisc::summarySE(Dados,
measurevar="MCT",
groupvars=c("Sexo","Sedentarismo"),
na.rm=TRUE)
pd <- ggplot2::position_dodge(0.5) # move them .05 to the left and right
# Intervalo de predicao >= 75%: Standard deviation (sd)
k75 <- 2
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
colour=Sedentarismo)) +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-k75*sd,
ymax=MCT+k75*sd),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de predição de >= 75% de MCT") +
ggplot2::ylim(30,100) +
ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(),
legend.justification=c(1,0),
legend.position=c(1,0))`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
# 75% Error bar: sd
k75 <- 2
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
# group=Sedentarismo,
colour=Sedentarismo)) +
ggplot2::geom_bar(position=ggplot2::position_dodge(1),
stat="identity",
width=.8,
fill="white") +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-k75*sd,
ymax=MCT+k75*sd),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de predição >= 75% de MCT") +
ggplot2::scale_y_continuous(breaks=0:100*10) +
ggplot2::coord_cartesian(ylim=c(30,100)) +
ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
Intervalo de confiança da média
Wikipédia: Intervalo de confiança
Referências
- 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.
- Banjanovic, ES & Osborne, JW (2016) Confidence intervals for effect sizes: Applying bootstrap resampling. Practical Assessment, Research & Evaluation 21(5).
- 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.
Estimação de intervalo de confiança
- MCT
- sd: desvio-padrão
- se: erro-padrão (desvio-padrão da média)
| Sexo | \(n\) | média | sd | se |
|---|---|---|---|---|
| F | 230 | 57.6 | 9.1 | 0.6 |
| M | 312 | 71.6 | 12.1 | 0.7 |
# IC95% por formula
## MCT Feminino
n.F <- sum(!is.na(Dados.F$MCT))
m.F <- mean(Dados.F$MCT, na.rm=TRUE)
sd.F <- sd(Dados.F$MCT, na.rm=TRUE)
se.F <- sd.F/sqrt(n.F)
cat("Erro-padrao da media de MCT Feminino = ",
round(se.F,1), "\n", sep="")Erro-padrao da media de MCT Feminino = 0.6me.F <- qt(p=1-alfa/2, df=n.F-1)*se.F
cat("Margem de erro da media de MCT Feminino = ",
round(me.F,1), "\n", sep="")Margem de erro da media de MCT Feminino = 1.2IC95(media de MCT Feminino) = [ 56.5 58.8 ]gl.F <- n.F-1
qt1.F <- round(qt(p=alfa/2, df=gl.F),2)
qt2.F <- round(qt(p=1-alfa/2, df=gl.F),2)
DescTools::PlotProbDist(breaks=c(-5, qt1.F, 0, qt2.F, 5),
function(x) dt(x, df=gl.F),
blab=c(qt1.F, "", qt2.F),
xlim=c(-5,5),
alab=c(alfa/2,"","", alfa/2),
col=c("black", "black"),
main=paste0("MCT Feminino: t(",gl.F,")"),
density=c(20, 0, 0, 20))
## MCT Masculino
n.M <- sum(!is.na(Dados.M$MCT))
m.M <- mean(Dados.M$MCT, na.rm=TRUE)
sd.M <- sd(Dados.M$MCT, na.rm=TRUE)
se.M <- sd.M/sqrt(n.M)
cat("Erro-padrao da media de MCT Masculino = ",
round(se.M,1), "\n", sep="")Erro-padrao da media de MCT Masculino = 0.7me.M <- qt(p=1-alfa/2, df=n.M-1)*se.M
cat("Margem de erro da media de MCT Masculino = ",
round(me.M,1), "\n", sep="")Margem de erro da media de MCT Masculino = 1.3IC95(media de MCT Masculino) = [ 70.2 72.9 ]gl.M <- n.M-1
qt1.M <- round(qt(p=alfa/2, df=gl.M),2)
qt2.M <- round(qt(p=1-alfa/2, df=gl.M),2)
DescTools::PlotProbDist(breaks=c(-5, qt1.M, 0, qt2.M, 5),
function(x) dt(x, df=gl.M),
blab=c(qt1.M, "", qt2.M),
xlim=c(-5,5),
alab=c(alfa/2,"","", alfa/2),
col=c("black", "black"),
main=paste0("MCT Masculino: t(",gl.M,")"),
density=c(20, 0, 0, 20))
Fórmula do intervalo de confiança da média
\[ \Large \text{IC}^{95\%}(\mu) = \left[\bar{x} \pm t^{0.975}_{n-1}\frac{s}{\sqrt{n}}\right] \]
- \(\mu\): média populacional
- \(n\): tamanho da amostra
- \(\bar{x}\): média amostral
- \(s\): desvio-padrão amostral
Quantil da distribuição t: \(t^{0.975}_{n-1} =\)
qt(p=.975, df=n-1):
- \(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
DescTools::MeanCI: estimação de intervalo de
confiança
mean lwr.ci upr.ci
57.6 56.5 58.8 mean lwr.ci upr.ci
71.6 70.2 72.9 [1] 56.5 58.8[1] 70.2 72.9# IC95% por DescTools::MeanCI com desvio-padrao conhecido
round(DescTools::MeanCI(Dados.F$MCT, sd=10, na.rm=TRUE), 1) mean lwr.ci upr.ci
57.6 56.3 58.9 mean lwr.ci upr.ci
71.6 70.0 73.1 # 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,
na.rm=TRUE), digits=4)$F
mean lwr.ci upr.ci
57.63 56.45 58.81
$M
mean lwr.ci upr.ci
71.59 70.24 72.94 # IC95% por bootstrapping para sexo feminino
round(DescTools::MeanCI(Dados.F$MCT,
method="boot",
type="perc",
R=1e4,
na.rm=TRUE), 1) mean lwr.ci upr.ci
57.6 56.5 58.8 # IC95% de MCT, Estatura e IMC para sexo feminino
round(do.call("rbind", lapply(Dados.F[, c("MCT","Estatura","IMC")],
FUN=DescTools::MeanCI,
na.rm=TRUE)), 1) mean lwr.ci upr.ci
MCT 57.6 56.5 58.8
Estatura 164.3 163.5 165.2
IMC 21.2 20.9 21.6 mean lwr.ci upr.ci
MCT 57.6 56.5 58.8
Estatura 164.3 163.5 165.2
IMC 21.2 20.9 21.6Gráfico de intervalo de confiança de 95% de média: delineamento entre participantes
# summarySE provides the standard deviation, standard error of the mean, and a (default 95%) confidence interval
s <- Rmisc::summarySE(Dados,
measurevar="MCT",
groupvars=c("Sexo"),
na.rm=TRUE)
# Intervalo de confianca de 95%: Standard error of the mean (se)
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT)) +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-ci,
ymax=MCT+ci),
width=.3) +
ggplot2::geom_point(size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de confiança de 95% de MCT") +
ggplot2::ylim(50,80) +
ggplot2::theme_bw() +
ggplot2::theme(legend.justification=c(1,0),
legend.position=c(1,0),
panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())
s <- Rmisc::summarySE(Dados,
measurevar="MCT",
groupvars=c("Sexo","Sedentarismo"),
na.rm=TRUE)
# Intervalo de confianca de 95%: Standard error of the mean (se)
pd <- ggplot2::position_dodge(0.5) # move them .05 to the left and right
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
colour=Sedentarismo)) +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-ci,
ymax=MCT+ci),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de confiança de 95% de MCT") +
ggplot2::ylim(50,80) +
ggplot2::theme_bw() +
ggplot2::theme(legend.justification=c(1,0),
legend.position=c(1,0),
panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
# 95% Error bar: se
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
# group=Sedentarismo,
colour=Sedentarismo)) +
ggplot2::geom_bar(position=ggplot2::position_dodge(1),
stat="identity",
width=.8,
fill="white") +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-ci,
ymax=MCT+ci),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de confiança de 95% de MCT") +
ggplot2::scale_y_continuous(breaks=0:100*10) +
ggplot2::coord_cartesian(ylim=c(50,80)) +
ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
Gráfico de intervalo de confiança de 95% de média: delineamento intraparticipantes
alfa <- 0.05
# Within-subject CI
dfw <- read.table(header=TRUE,
text='
subject pretest posttest
1 59.4 64.5
2 46.4 52.4
3 46.0 49.7
4 49.0 48.7
5 32.5 37.4
6 45.2 49.5
7 60.3 59.9
8 54.3 54.1
9 45.4 49.6
10 38.9 48.5
')
# Treat subject ID as a factor
dfw$subject <- factor(dfw$subject)
# Convert to long format
dfw_long <- reshape2::melt(dfw,
id.vars="subject",
measure.vars=c("pretest","posttest"),
variable.name="condition")
print(dfw_long) 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.5dfwc <- summarySEwithin2(dfw_long,
measurevar="value",
withinvars="condition",
idvar="subject",
na.rm=TRUE,
conf.interval=1-alfa/2)
print(dfwc) 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# Make the graph with the 95% confidence interval (Bonferroni)
ggplot2::ggplot(dfwc,
ggplot2::aes(x=condition,
y=value)) +
ggplot2::geom_line() +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=value-ci,
ymax=value+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylim(40,60) +
ggplot2::ylab("Value") +
ggplot2::ggtitle("Within-subject CI95 Bonferroni") +
ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
# Instead of summarySEwithin, use summarySE, which treats condition as though it were a between-subjects variable
dfwc_between <- Rmisc::summarySE(data=dfw_long,
measurevar="value",
groupvars="condition",
na.rm=FALSE,
conf.interval=1-alfa/2)
print(dfwc_between) 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# Show the between-S CI's in red, and the within-S CI's in black
ggplot2::ggplot(dfwc_between,
ggplot2::aes(x=condition,
y=value)) +
ggplot2::geom_line() +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=value-ci,
ymax=value+ci),
colour="darkgray",
data=dfwc_between) +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=value-ci,
ymax=value+ci),
colour="black",
data=dfwc) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylim(40, 60) +
ggplot2::ylab("Value") +
ggplot2::ggtitle("Within vs. Between-subject CI95 Bonferroni") +
ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
# Use a consistent y range
ymax <- max(dfw_long$value)
ymin <- min(dfw_long$value)
# Plot the individuals
ggplot2::ggplot(dfw_long,
ggplot2::aes(x=condition,
y=value,
colour=subject,
group=subject)) +
ggplot2::geom_line() +
ggplot2::geom_point(shape=21,
fill="white") +
ggplot2::ylim(ymin,ymax) +
ggplot2::ylab("Value") +
ggplot2::ggtitle("Subject profile") +
ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())
# Two within-subjects variables
data <- read.table(header=TRUE,
text='
Subject RoundMono SquareMono RoundColor SquareColor
1 41 40 41 37
2 57 56 56 53
3 52 53 53 50
4 49 47 47 47
5 47 48 48 47
6 37 34 35 36
7 47 50 47 46
8 41 40 38 40
9 48 47 49 45
10 37 35 36 35
11 32 31 31 33
12 47 42 42 42
')
# Convert it to long format
data_long <- reshape2::melt(data=data,
id.var="Subject",
measure.vars=c("RoundMono",
"SquareMono",
"RoundColor",
"SquareColor"),
variable.name="Condition")
names(data_long)[names(data_long)=="value"] <- "Time"
# Split Condition column into Shape and ColorScheme
data_long$Shape <- NA
data_long$Shape[grepl("^Round", data_long$Condition)] <- "Round"
data_long$Shape[grepl("^Square", data_long$Condition)] <- "Square"
data_long$Shape <- factor(data_long$Shape)
data_long$ColorScheme <- NA
data_long$ColorScheme[grepl("Mono$", data_long$Condition)] <- "Monochromatic"
data_long$ColorScheme[grepl("Color$", data_long$Condition)] <- "Colored"
data_long$ColorScheme <- factor(data_long$ColorScheme, levels=c("Monochromatic","Colored"))
# Remove the Condition column now
data_long$Condition <- NULL
# Look at first few rows
print(head(data_long)) 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 Monochromaticdatac <- summarySEwithin2(data_long,
measurevar="Time",
withinvars=c("Shape",
"ColorScheme"),
idvar="Subject",
na.rm=TRUE)
print(datac) Shape ColorScheme N Time TimeNormed sd se ci
1 Round Colored 12 43.58333 43.58333 1.212311 0.3499639 0.7702654
2 Round Monochromatic 12 44.58333 44.58333 1.331438 0.3843531 0.8459554
3 Square Colored 12 42.58333 42.58333 1.461630 0.4219364 0.9286757
4 Square Monochromatic 12 43.58333 43.58333 1.261312 0.3641095 0.8013997# Within-subject Error Bar 95%
ggplot2::ggplot(datac,
ggplot2::aes(x=Shape,
y=Time,
fill=ColorScheme)) +
ggplot2::geom_bar(position=ggplot2::position_dodge(.9),
colour="black",
stat="identity") +
ggplot2::geom_errorbar(position=ggplot2::position_dodge(.9),
width=.25,
ggplot2::aes(ymin=Time-ci,
ymax=Time+ci)) +
ggplot2::coord_cartesian(ylim=c(40,46)) +
ggplot2::scale_fill_manual(values=c("#CCCCCC",
"#FFFFFF")) +
ggplot2::scale_y_continuous(breaks=seq(1:100)) +
ggplot2::theme_bw() +
ggplot2::ggtitle("Within-subject CI95") +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())
# Within-subject CI95
ggplot2::ggplot(datac,
ggplot2::aes(x=Shape,
y=Time,
colour=ColorScheme)) +
ggplot2::geom_errorbar(position=ggplot2::position_dodge(.9),
width=.25,
ggplot2::aes(ymin=Time-ci,
ymax=Time+ci)) +
ggplot2::geom_line(position=ggplot2::position_dodge(.9)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white",
position=ggplot2::position_dodge(.9)) +
ggplot2::coord_cartesian(ylim=c(41,46)) +
ggplot2::scale_fill_manual(values=c("#CCCCCC",
"#FFFFFF")) +
ggplot2::scale_y_continuous(breaks=seq(1:100)) +
ggplot2::theme_bw() +
ggplot2::ggtitle("Within-subject CI95") +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank())`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
HH::CIplot: simulação de intervalo de confiança
Simulação de intervalo de confiança em Rpsychologist
DescTools::MedianCI: IC de mediana
median lwr.ci upr.ci
56 55 58
attr(,"conf.level")
[1] 0.9592875median lwr.ci upr.ci
72 70 73
attr(,"conf.level")
[1] 0.9526369DescTools::VarCI: IC de desvio-padrão
var lwr.ci upr.ci
9.052874 8.294251 9.965432 var lwr.ci upr.ci
12.08723 11.20739 13.11815 Validade do intervalo de confiança
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.
mean lwr.ci upr.ci
57.6 56.5 58.8 mean lwr.ci upr.ci
71.6 70.2 72.9 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.
Os intervalos de confiança da média populacional por fórmula é válida, 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 reamostragem (bootstrapping). TCL e bootstrapping não necessitam da suposição de normalidade para a validade do intervalo de confiança.
Teorema Central do Limite
O Teorema Central do Limite (TCL) para a média amostral assume:
- Delineamento entre participantes: observações independentes
- Variável intervalar
- Distribuição qualquer da variável
- Amostra com pelo menos 30 observações
TCL: 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.
Fonte: Wonnacott & Wonnacott (1969)
Wikipedia: Central Limit Theorem (CLT)
- Referências
- Krzywinski, M & Altman, N (2013) Importance of being uncertain. Nature Methods 10(9): 809-10.
- 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.
Reamostragem estatística
Bootstrapping
Suposições para validade do intervalo de confiança por bootstrapping:
- Variável intervalar
- Observações independentes
Referências
- Carpenter, J & Bithell, John (2000) Bootstrap confidence intervals: when, which, what? A practical guide for medical statisticians. Statistics in Medicine 19(9): 1141-64. https://doi.org/10.1002/(SICI)1097-0258(20000515)19:9%3C1141::AID-SIM479%3E3.0.CO;2-F
- Johnson, RW (2001) An introduction to the bootstrap. Teaching Statistics 23(2): 49-54. https://onlinelibrary.wiley.com/doi/10.1111/1467-9639.00050
# IC95% por bootstrapping da média para sexo feminino
round(DescTools::MeanCI(Dados.F$MCT,
method="boot",
type="bca",
R=1e5,
na.rm=TRUE), 1) mean lwr.ci upr.ci
57.6 56.6 59.0 # IC95% por bootstrapping da média para sexo masculino
round(DescTools::MeanCI(Dados.M$MCT,
method="boot",
type="bca",
R=1e5,
na.rm=TRUE), 1) mean lwr.ci upr.ci
71.6 70.3 73.0 # IC95% por bootstrapping da mediana para sexo feminino
round(DescTools::MedianCI(Dados.F$MCT,
method="boot",
R=1e5,
na.rm=TRUE), 1) median lwr.ci upr.ci
56 55 58 # IC95% por bootstrapping da mediana para sexo masculino
round(DescTools::MedianCI(Dados.M$MCT,
method="boot",
R=1e5,
na.rm=TRUE), 1) median lwr.ci upr.ci
72 70 73 # IC95% por bootstrapping do desvio-padrao para sexo feminino
round(sqrt(DescTools::VarCI(Dados.F$MCT,
method="perc",
R=1e5,
na.rm=TRUE)),1) var lwr.ci upr.ci
9.1 7.0 11.3 # IC95% por bootstrapping do desvio-padrao para sexo masculino
round(sqrt(DescTools::VarCI(Dados.M$MCT,
method="perc",
R=1e5,
na.rm=TRUE)),1) var lwr.ci upr.ci
12.1 10.7 13.5 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/
Código R
# "Módulo 2"
# José O Siqueira siqueira@usp.br
# Paulo S P Silveira silveira@usp.br
# "Tópicos Essenciais da Bioestatística"
# 2022
# Descriptive Statistics and Graphics
## http://www.sthda.com/english/wiki/descriptive-statistics-and-graphics
# Be Awesome in ggplot2: A Practical Guide to be Highly Effective - R software and data visualization
## http://www.sthda.com/english/wiki/be-awesome-in-ggplot2-a-practical-guide-to-be-highly-effective-r-software-and-data-visualization
# Plotting means and error bars (ggplot2)
## http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
# Carregar pacotes
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(HH, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(Rmisc, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(ggpubr, warn.conflicts=FALSE))
source("summarySEwithin2.R")
# Ler planilha
Dados <- data.frame(readxl::read_excel(path="Biometria_FMUSP.xlsx",
sheet="dados",
na=c("NA","Na","nA","na")))
Dados$MCT[Dados$MCT==max(Dados$MCT,na.rm=TRUE)] <- NA
Dados$Estatura[Dados$Estatura==min(Dados$Estatura,na.rm=TRUE)] <- NA
Dados$ID <- factor(Dados$ID)
Dados$Ano <- factor(Dados$Ano)
Dados$Turma <- factor(Dados$Turma)
Dados$Sexo <- factor(Dados$Sexo)
Dados$Mao <- factor(Dados$Mao)
Dados$TipoSang <- factor(Dados$TipoSang)
Dados$ABO <- factor(Dados$ABO)
Dados$Rh <- factor(Dados$Rh)
Dados$AtivFisica <- factor(Dados$AtivFisica)
Dados$Sedentarismo <- factor(Dados$Sedentarismo)
Dados$IMC <- Dados$MCT/(Dados$Estatura/100)^2
Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")
# Estimação de intervalo de confiança de MCT
# group | n | mean | sd| se|
# :----:|----:|------:|-----:|---:|
# F | 230 | 57.6 | 9.1| 0.6|
# M | 312 | 71.6 | 12.1| 0.7|
# IC95% por formula
alfa <- 0.05
## MCT Feminino
n.F <- sum(!is.na(Dados.F$MCT))
m.F <- mean(Dados.F$MCT, na.rm=TRUE)
sd.F <- sd(Dados.F$MCT, na.rm=TRUE)
se.F <- sd.F/sqrt(n.F)
cat("Erro-padrao da media de MCT Feminino = ",
round(se.F,1), "\n", sep="")
me.F <- qt(p=1-alfa/2, df=n.F-1)*se.F
cat("Margem de erro da media de MCT Feminino = ",
round(me.F,1), "\n", sep="")
ic.F <- m.F+c(-1,1)*me.F
cat("IC95(media de MCT Feminino) = [",
round(ic.F,1), "]\n")
## MCT Masculino
n.M <- sum(!is.na(Dados.M$MCT))
m.M <- mean(Dados.M$MCT, na.rm=TRUE)
sd.M <- sd(Dados.M$MCT, na.rm=TRUE)
se.M <- sd.M/sqrt(n.M)
cat("Erro-padrao da media de MCT Masculino = ",
round(se.M,1), "\n", sep="")
me.M <- qt(p=1-alfa/2, df=n.M-1)*se.M
cat("Margem de erro da media de MCT Masculino = ",
round(me.M,1), "\n", sep="")
ic.M <- m.M+c(-1,1)*me.M
cat("IC95(media de MCT Masculino) = [",
round(ic.M,1), "]\n")
# plot t-distribution
gl.F <- n.F-1
qt1.F <- round(qt(p=alfa/2, df=gl.F),2)
qt2.F <- round(qt(p=1-alfa/2, df=gl.F),2)
DescTools::PlotProbDist(breaks=c(-5, qt1.F, 0, qt2.F, 5),
function(x) dt(x, df=gl.F),
blab=c(qt1.F, "", qt2.F),
xlim=c(-5,5),
alab=c(alfa/2,"","", alfa/2),
col=c("black", "black"),
main=paste0("MCT Feminino: t(",gl.F,")"),
density=c(20, 0, 0, 20))
gl.M <- n.M-1
qt1.M <- round(qt(p=alfa/2, df=gl.M),2)
qt2.M <- round(qt(p=1-alfa/2, df=gl.M),2)
DescTools::PlotProbDist(breaks=c(-5, qt1.M, 0, qt2.M, 5),
function(x) dt(x, df=gl.M),
blab=c(qt1.M, "", qt2.M),
xlim=c(-5,5),
alab=c(alfa/2,"","", alfa/2),
col=c("black", "black"),
main=paste0("MCT Masculino: t(",gl.M,")"),
density=c(20, 0, 0, 20))
# `DescTools::MeanCI`: estimação de intervalo de confiança
# IC95% por DescTools::MeanCI
round(DescTools::MeanCI(Dados.F$MCT, na.rm=TRUE),1)
round(DescTools::MeanCI(Dados.M$MCT, na.rm=TRUE),1)
# Os mesmos IC95% por t.test
round(t.test(Dados.F$MCT)$conf.int[1:2],1)
round(t.test(Dados.M$MCT)$conf.int[1:2],1)
# IC95% por DescTools::MeanCI com desvio-padrao conhecido
round(DescTools::MeanCI(Dados.F$MCT, sd=10, na.rm=TRUE), 1)
round(DescTools::MeanCI(Dados.M$MCT, sd=14, na.rm=TRUE), 1)
# 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,
na.rm=TRUE)
# IC95% por bootstrapping para sexo feminino
round(DescTools::MeanCI(Dados.F$MCT,
method="boot",
type="perc",
R=1e4,
na.rm=TRUE), 1)
# IC95% de MCT, Estatura e IMC para sexo feminino
round(do.call("rbind", lapply(Dados.F[, 10:12],
FUN=DescTools::MeanCI,
na.rm=TRUE)), 1)
t(round(sapply(Dados.F[,c("MCT","Estatura","IMC")],
FUN=MeanCI,
na.rm=TRUE), 1))
# Cookbook for R: Plotting means and error bars (ggplot2)
# http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
# summarySE provides the standard deviation, standard error of the mean, and a (default 95%) confidence interval
s <- Rmisc::summarySE(Dados,
measurevar="MCT",
groupvars=c("Sexo","Sedentarismo"),
na.rm=TRUE)
# Intervalo de confianca de 95%: Standard error of the mean (se)
pd <- ggplot2::position_dodge(0.5) # move them .05 to the left and right
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
colour=Sedentarismo)) +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-ci,
ymax=MCT+ci),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de confiança de 95% de MCT") +
ggplot2::ylim(50,80) +
ggplot2::theme_bw() +
ggplot2::theme(legend.justification=c(1,0),
legend.position=c(1,0))
# 95% Error bar: se
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
# group=Sedentarismo,
colour=Sedentarismo)) +
ggplot2::geom_bar(position=ggplot2::position_dodge(1),
stat="identity",
width=.8,
fill="white") +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-ci,
ymax=MCT+ci),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de confiança de 95% de MCT") +
ggplot2::scale_y_continuous(breaks=0:100*10) +
ggplot2::coord_cartesian(ylim=c(50,80)) +
ggplot2::theme_bw()
# Intervalo de predicao >= 75%: Standard deviation (sd)
k75 <- 2
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
colour=Sedentarismo)) +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-k75*sd,
ymax=MCT+k75*sd),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de predição de 95% de MCT") +
ggplot2::ylim(30,100) +
ggplot2::theme_bw() +
ggplot2::theme(legend.justification=c(1,0),
legend.position=c(1,0))
# 95% Error bar: sd
k75 <- 2
ggplot2::ggplot(s,
ggplot2::aes(x=Sexo,
y=MCT,
# group=Sedentarismo,
colour=Sedentarismo)) +
ggplot2::geom_bar(position=ggplot2::position_dodge(1),
stat="identity",
width=.8,
fill="white") +
ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-k75*sd,
ymax=MCT+k75*sd),
width=.3,
position=pd) +
ggplot2::geom_line(position=pd) +
ggplot2::geom_point(position=pd,
size=3,
shape=21) +
ggplot2::xlab("Sexo") +
ggplot2::ylab("MCT (kg)") +
ggplot2::ggtitle("Intervalo de predição >= 75% de MCT") +
ggplot2::scale_y_continuous(breaks=0:100*10) +
ggplot2::coord_cartesian(ylim=c(30,100)) +
ggplot2::theme_bw()
# Within-subject CI
dfw <- read.table(header=TRUE,
text='
subject pretest posttest
1 59.4 64.5
2 46.4 52.4
3 46.0 49.7
4 49.0 48.7
5 32.5 37.4
6 45.2 49.5
7 60.3 59.9
8 54.3 54.1
9 45.4 49.6
10 38.9 48.5
')
# Treat subject ID as a factor
dfw$subject <- factor(dfw$subject)
# Convert to long format
dfw_long <- reshape2::melt(dfw,
id.vars="subject",
measure.vars=c("pretest","posttest"),
variable.name="condition")
dfw_long
#> 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
#>
dfwc <- summarySEwithin2(dfw_long,
measurevar="value",
withinvars="condition",
idvar="subject",
na.rm=TRUE,
conf.interval=.95)
dfwc
#> condition N value value_norm sd se ci
#> 1 posttest 10 51.43 51.43 2.262361 0.7154214 1.618396
#> 2 pretest 10 47.74 47.74 2.262361 0.7154214 1.618396
#>
# Make the graph with the 95% confidence interval
ggplot2::ggplot(dfwc,
ggplot2::aes(x=condition,
y=value)) +
ggplot2::geom_line() +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=value-ci,
ymax=value+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylim(40,60) +
ggplot2::ylab("Value") +
ggplot2::ggtitle("Within-subject CI95") +
ggplot2::theme_bw()
# Instead of summarySEwithin, use summarySE, which treats condition as though it were a between-subjects variable
dfwc_between <- Rmisc::summarySE(data=dfw_long,
measurevar="value",
groupvars="condition",
na.rm=FALSE,
conf.interval=.95)
dfwc_between
#> condition N value sd se ci
#> 1 pretest 10 47.74 8.598992 2.719240 6.151348
#> 2 posttest 10 51.43 7.253972 2.293907 5.189179
# Show the between-S CI's in red, and the within-S CI's in black
ggplot2::ggplot(dfwc_between,
ggplot2::aes(x=condition,
y=value)) +
ggplot2::geom_line() +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=value-ci,
ymax=value+ci),
colour="gray",
data=dfwc_between) +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=value-ci,
ymax=value+ci),
colour="black",
data=dfwc) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylim(40, 60) +
ggplot2::ylab("Value") +
ggplot2::ggtitle("Within vs. Between-subject CI95") +
ggplot2::theme_bw()
# Use a consistent y range
ymax <- max(dfw_long$value)
ymin <- min(dfw_long$value)
# Plot the individuals
ggplot2::ggplot(dfw_long,
ggplot2::aes(x=condition,
y=value,
colour=subject,
group=subject)) +
ggplot2::geom_line() +
ggplot2::geom_point(shape=21,
fill="white") +
ggplot2::ylim(ymin,ymax) +
ggplot2::ylab("Value") +
ggplot2::ggtitle("Subject profile") +
ggplot2::theme_classic()
# Two within-subjects variables
data <- read.table(header=TRUE,
text='
Subject RoundMono SquareMono RoundColor SquareColor
1 41 40 41 37
2 57 56 56 53
3 52 53 53 50
4 49 47 47 47
5 47 48 48 47
6 37 34 35 36
7 47 50 47 46
8 41 40 38 40
9 48 47 49 45
10 37 35 36 35
11 32 31 31 33
12 47 42 42 42
')
# Convert it to long format
data_long <- reshape2::melt(data=data,
id.var="Subject",
measure.vars=c("RoundMono",
"SquareMono",
"RoundColor",
"SquareColor"),
variable.name="Condition")
names(data_long)[names(data_long)=="value"] <- "Time"
# Split Condition column into Shape and ColorScheme
data_long$Shape <- NA
data_long$Shape[grepl("^Round", data_long$Condition)] <- "Round"
data_long$Shape[grepl("^Square", data_long$Condition)] <- "Square"
data_long$Shape <- factor(data_long$Shape)
data_long$ColorScheme <- NA
data_long$ColorScheme[grepl("Mono$", data_long$Condition)] <- "Monochromatic"
data_long$ColorScheme[grepl("Color$", data_long$Condition)] <- "Colored"
data_long$ColorScheme <- factor(data_long$ColorScheme, levels=c("Monochromatic","Colored"))
# Remove the Condition column now
data_long$Condition <- NULL
# Look at first few rows
head(data_long)
#> 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
datac <- summarySEwithin2(data_long,
measurevar="Time",
withinvars=c("Shape",
"ColorScheme"),
idvar="Subject",
na.rm=TRUE)
datac
#> Shape ColorScheme N Time Time_norm sd se ci
#> 1 Round Colored 12 43.58333 43.58333 1.212311 0.3499639 0.7702654
#> 2 Round Monochromatic 12 44.58333 44.58333 1.331438 0.3843531 0.8459554
#> 3 Square Colored 12 42.58333 42.58333 1.461630 0.4219364 0.9286757
#> 4 Square Monochromatic 12 43.58333 43.58333 1.261312 0.3641095 0.8013997
# Within-subject Error Bar 95%
ggplot2::ggplot(datac,
ggplot2::aes(x=Shape,
y=Time,
fill=ColorScheme)) +
ggplot2::geom_bar(position=position_dodge(.9),
colour="black",
stat="identity") +
ggplot2::geom_errorbar(position=position_dodge(.9),
width=.25,
ggplot2::aes(ymin=Time-ci,
ymax=Time+ci)) +
ggplot2::coord_cartesian(ylim=c(40,46)) +
ggplot2::scale_fill_manual(values=c("#CCCCCC",
"#FFFFFF")) +
ggplot2::scale_y_continuous(breaks=seq(1:100)) +
ggplot2::theme_bw() +
ggplot2::ggtitle("Within-subject CI95")
# Within-subject CI95
ggplot2::ggplot(datac,
ggplot2::aes(x=Shape,
y=Time,
colour=ColorScheme)) +
ggplot2::geom_errorbar(position=position_dodge(.9),
width=.25,
ggplot2::aes(ymin=Time-ci,
ymax=Time+ci)) +
ggplot2::geom_line(position=position_dodge(.9)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white",
position=position_dodge(.9)) +
ggplot2::coord_cartesian(ylim=c(41,46)) +
ggplot2::scale_fill_manual(values=c("#CCCCCC",
"#FFFFFF")) +
ggplot2::scale_y_continuous(breaks=seq(1:100)) +
ggplot2::theme_bw() +
ggplot2::ggtitle("Within-subject CI95")
# `HH::CIplot`: simulação de intervalo de confiança
HH::CIplot(n.intervals=1e2,
n.per.row=40,
pop.mean=175,
pop.sd=10,
conf.level=0.95)
HH::CIplot(n.intervals=1e4,
n.per.row=40,
pop.mean=175,
pop.sd=10,
conf.level=0.95)
set.seed(9342)
mp <- 175 # media populacional
dpp <- 8 # desvio-padrao populacional
cat ("Média populacional = ",mp,"\n",sep="")
cat ("Desvio padrão populacional = ",dpp,"\n",sep="")
n <- 30 # tamanho da amostra
cat("n =", n, "\n")
alfa <- 0.05 # nivel de significancia
plot(x = c(mp-dpp, mp+dpp), y = c(1, 100), type = "n",
xlab = paste("Estatura masculina (média = ",mp," cm)",sep=""),
ylab = "",
axes=FALSE)
abline(v = mp, col = "red") # linha vertical da media populacional
# sampling
counter <- 0 # set counter to 0
B <- 1e5
for (i in 1:B)
{
x <- rnorm(n, mp, dpp) # amostra aleatoria normal de tamanho n
s <- sd(x)
L <- mean(x) - qt(1-alfa/2,n-1)*s/sqrt(n) # lower limit
U <- mean(x) + qt(1-alfa/2,n-1)*s/sqrt(n) # upper limit
if (L < mp && mp < U) # verifica se a media pop esta contida em IC95
{
color <- "#222222" # cinza
counter <- counter + 1 # Se sim, soma 1 no contador
espessura <- 0.5
} else
{
color <- "#1965B0" # azul cobalto
espessura <- 2
}
if (i <= 100) # plota os primeiros IC95%
segments(L, i, U, i, col=color, lwd=espessura)
}
pIC95 <- counter/B # % de IC95 que contem a media populacional.
cat("\nProporcao dos ",
format(B, scientific = FALSE),
" IC95(media pop) que contem a media poppulacional ",
mp, " = ",
round(pIC95,2), sep="")
# `DescTools::MedianCI`: IC de mediana
DescTools::MedianCI(Dados.F$MCT, na.rm=TRUE)
DescTools::MedianCI(Dados.M$MCT, na.rm=TRUE)
# `DescTools::VarCI`: IC de desvio-padrão
sqrt(DescTools::VarCI(Dados.F$MCT, na.rm=TRUE))
sqrt(DescTools::VarCI(Dados.M$MCT, na.rm=TRUE))
# Variável com distribuição normal
# Exemplo de distribuição de normal: 1
## Pressão arterial sistólica de jovens saudáveis
source("exemplo05_normal.R")
# Distribuição normal: Questão 1
# Qual é a probabilidade que nenhum paciente tenha atendimento completo durante uma madrugada no pronto-atendimento dessa pequena cidade?
# R.: $P(X > 140) = 0.023$
pnorm(q=140, mean=120, sd=10, lower.tail=FALSE)
# OU
1-pnorm(q=140, mean=120, sd=10)
# Distribuição normal: Questão 2
# Qual é a probabilidade dos valores de PAS de jovens sadios estarem entre 100 e 140 mmHg?
# R.: $P(100 < X < 140) = 0.9545$
pnorm(q=140, mean=120, sd=10)-pnorm(q=100, mean=120, sd=10)
# Distribuição normal: Questão 3
# Quais são os limites de um intervalo simétrico em relação à média que engloba 95% dos valores de PAS de jovens sadios?
# R.: $P(100.4 < X < 139.6) = 0.95$
qnorm(p=0.025, mean=120, sd=10)
qnorm(p=0.975, mean=120, sd=10)
DescTools::MeanCI(Dados.F$MCT, na.rm=TRUE)
DescTools::MeanCI(Dados.M$MCT, na.rm=TRUE)
describeBy(Dados$MCT,group=Dados$Sexo,mat=2)
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")
m.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")
print(by(Dados$MCT, Dados$Sexo, shapiro.test))
# Teorema Central do Limite
# IC95% por bootstrapping da média para sexo feminino
round(DescTools::MeanCI(Dados.F$MCT,
method="boot",
type="perc",
R=1e5,
na.rm=TRUE), 1)
# IC95% por bootstrapping da média para sexo masculino
round(DescTools::MeanCI(Dados.M$MCT,
method="boot",
type="perc",
R=1e5,
na.rm=TRUE), 1)
# IC95% por bootstrapping da mediana para sexo feminino
round(DescTools::MedianCI(Dados.F$MCT,
method="boot",
R=1e5,
na.rm=TRUE), 1)
# IC95% por bootstrapping da mediana para sexo masculino
round(DescTools::MedianCI(Dados.M$MCT,
method="boot",
R=1e5,
na.rm=TRUE), 1)
# IC95% por bootstrapping do desvio-padrao para sexo feminino
round(sqrt(DescTools::VarCI(Dados.F$MCT,
method="perc",
R=1e5,
na.rm=TRUE)),1)
# IC95% por bootstrapping do desvio-padrao para sexo masculino
round(sqrt(DescTools::VarCI(Dados.M$MCT,
method="perc",
R=1e5,
na.rm=TRUE)),1)