MAN5728: Tópicos Essenciais da Bioestatística
Módulo 2: Intervalo de Confiança
07 julho 2025 18:43h
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))options(warn=-1)
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))
source("summarySEwithin2.R")Material
- HTML de R Markdown 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 (missing = NA):
- 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
Dados <- 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 vs. Intervalo de confiança
Intervalo de predição da variável
\[\Large 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
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. > Desigualdade de Chebyshev
“Teste” de normalidade
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(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()`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)
| group | n | mean | 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 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
tapply(Dados$MCT,
Dados$Sexo,
FUN=DescTools::MeanCI,
na.rm=TRUE)$F
mean lwr.ci upr.ci
57.63043 56.45426 58.80661
$M
mean lwr.ci upr.ci
71.58974 70.24329 72.93620 # 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.9 # 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))
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))`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()`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)
Anexando pacote: 'data.table'Os seguintes objetos são mascarados por 'package:reshape2':
dcast, meltO seguinte objeto é mascarado por 'package:DescTools':
%like%O seguinte objeto é mascarado por 'package:HH':
transpose
Anexando pacote: 'dplyr'Os seguintes objetos são mascarados por 'package:data.table':
between, first, lastOs seguintes objetos são mascarados por 'package:plyr':
arrange, count, desc, failwith, id, mutate, rename, summarise,
summarizeO seguinte objeto é mascarado por 'package:gridExtra':
combineO seguinte objeto é mascarado por 'package:MASS':
selectOs seguintes objetos são mascarados por 'package:stats':
filter, lagOs seguintes objetos são mascarados por 'package:base':
intersect, setdiff, setequal, union 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()`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()`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()
# 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=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")`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
Média populacional = 175desvio-padrão populacional = 8n = 30 
Proporcao dos 100000 IC95(media pop) que contem a media poppulacional 175 = 0.95set.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)\n",
"que contem a media pop ",
mp, " = ",
round(pIC95,2), sep="")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 Distribuição normal (gaussiana)
- Distribuição de probabilidade de variável contínua ilimitada bilateralmente \(X\)
- \(X \sim normal(\mu,\sigma)\): variável aleatória \(X\) representa valores contínuos reais
- \(X\in\mathbb{R}\)
- Média de \(X\): \(\mu\)
- Desvio-padrão de \(X\): \(\sigma\)
- O desvio-padrão, \(\sigma\), é a distância entre a média, \(\mu\), e o ponto de inflexão da curva da normal
- A curva da normal é unimodal
- A curva da normal é simétrica
q1 <- round(qnorm(p=alfa/2, mean=120, sd=10),2)
q2 <- round(qnorm(p=1-alfa/2, mean=120, sd=10),2)
DescTools::PlotProbDist(breaks=c(80,q1,q2,160),
function(x) dnorm(x, mean=120, sd=10),
blab=c(q1,q2),
alab=c("",1-alfa,""),
xlim=c(80,160),
main="normal(120,10)",
col=c("black", "black"),
density=c(0,20,0))
Variável com distribuição normal
Variável bruta ou com transformação potência de Tukey (e.g., logaritmo, raiz quadrada)
Somatométrica
- Massa corporal total
- Comprimento total/ estatura
- Comprimento de garra, presa, unha, pelo
- Medidas fisiológicas: pressão sanguínea de humanos adultos
Psicológica
- Escore de teste cognitivo: QI
Física
- A soma de pequenos erros independentes de mensuração
Figure 1. Line graphs show smoothed weighted frequency distribution, median, and 90th percentile of systolic pressures for populations 18 to 29 years and 60 to 74 years of age in the United States, 1960 to 1991. (Source: Centers for Disease Control and Prevention, National Center for Health Statistics.); https://www.ahajournals.org/doi/10.1161/01.HYP.26.1.60
Distribuição normal padrão
Exemplo de distribuição de normal: 1
Pressão arterial sistólica de jovens saudáveis
A pressão arterial sistólica (PAS) medida em milímetro de mercúrio (mmHg) em jovens gozando de boa saúde tem distribuição \(normal(120,10)\).
Siqueira & Tibúrcio, 2011, p. 2010-1
X ~ Normal(120;10)
P(X < 90) = 0.001349898
P(X < 100) = 0.02275013
P(X < 110) = 0.1586553
P(X < 120) = 0.5
P(X < 130) = 0.8413447
P(X < 140) = 0.9772499
P(X < 150) = 0.9986501
P(110 < X < 130) = 0.6826895
P(100 < X < 140) = 0.9544997
P(100.4 < X < 139.6) = 0.9500042
P(90 < X < 150) = 0.9973002
P(113.25 < X < 126.75) = 0.5003242
sink("exemplo05_normal.txt")
m <- 120
dp <- 10
cat("X ~ Normal(",m,";",dp,")\n\n", sep="")
cat("P(X < ",m-3*dp,") = ",pnorm(q=m-3*dp,mean=m,sd=dp),"\n", sep="")
cat("P(X < ",m-2*dp,") = ",pnorm(q=m-2*dp,mean=m,sd=dp),"\n", sep="")
cat("P(X < ",m-dp,") = ",pnorm(q=m-dp,mean=m,sd=dp),"\n", sep="")
cat("P(X < ",m,") = ",pnorm(q=m,mean=m,sd=dp),"\n", sep="")
cat("P(X < ",m+dp,") = ",pnorm(q=m+dp,mean=m,sd=dp),"\n", sep="")
cat("P(X < ",m+2*dp,") = ",pnorm(q=m+2*dp,mean=m,sd=dp),"\n", sep="")
cat("P(X < ",m+3*dp,") = ",pnorm(q=m+3*dp,mean=m,sd=dp),"\n", sep="")
cat("\nP(",m-dp," < X < ",m+dp,") = ",
pnorm(q=m+dp,mean=m,sd=dp)-pnorm(q=m-dp,mean=m,sd=dp),"\n", sep="")
cat("P(",m-2*dp," < X < ",m+2*dp,") = ",
pnorm(q=m+2*dp,mean=m,sd=dp)-pnorm(q=m-2*dp,mean=m,sd=dp),"\n", sep="")
qnorm(p=0.975,mean=0,sd=1)
cat("P(",m-1.96*dp," < X < ",m+1.96*dp,") = ",
pnorm(q=m+1.96*dp,mean=m,sd=dp)-pnorm(q=m-1.96*dp,mean=m,sd=dp),"\n", sep="")
cat("P(",m-3*dp," < X < ",m+3*dp,") = ",
pnorm(q=m+3*dp,mean=m,sd=dp)-pnorm(q=m-3*dp,mean=m,sd=dp),"\n", sep="")
qnorm(p=0.75,mean=0,sd=1)
cat("P(",m-0.675*dp," < X < ",m+0.675*dp,") = ",
pnorm(q=m+0.675*dp,mean=m,sd=dp)-pnorm(q=m-0.675*dp,mean=m,sd=dp),"\n", sep="")
sink()
png("exemplo05_normal.png")
rn <- rnorm(n=1e6,mean=m,sd=dp)
plot(density(rn),xlab="PAS(mmHg)",ylab="densidade",
main=paste0("Distribuicao normal(",m,";",dp,")"))
dev.off()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\)
[1] 0.02275013[1] 0.02275013DescTools::PlotProbDist(breaks=c(80,140,160),
function(x) dnorm(x, mean=120, sd=10),
blab=c(140),
alab=c("",round(p,3)),
xlim=c(80,160),
main="normal(120,10)",
col=c("black", "black"),
density=c(0,20))
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\)
[1] 0.9544997DescTools::PlotProbDist(breaks=c(80,100,140,160),
function(x) dnorm(x, mean=120, sd=10),
blab=c(100,140),
alab=c("",round(p,4),""),
xlim=c(80,160),
main="normal(120,10)",
col=c("black", "black"),
density=c(0,20,0))
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\)
[1] 100.4004[1] 139.5996q1 <- round(qnorm(p=alfa/2, mean=120, sd=10),2)
q2 <- round(qnorm(p=1-alfa/2, mean=120, sd=10),2)
DescTools::PlotProbDist(breaks=c(80,q1,q2,160),
function(x) dnorm(x, mean=120, sd=10),
blab=c(q1,q2),
alab=c("",1-alfa,""),
xlim=c(80,160),
main="normal(120,10)",
col=c("black", "black"),
density=c(0,20,0))
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
mean lwr.ci upr.ci
57.63043 56.45426 58.80661 mean lwr.ci upr.ci
71.58974 70.24329 72.93620 item group1 vars n mean sd median trimmed mad min max
X11 1 F 1 230 57.63043 9.052874 56 56.750 7.4130 41 120
X12 2 M 1 312 71.58974 12.087235 72 70.876 11.1195 43 125
range skew kurtosis se
X11 79 2.4309659 12.039133 0.5969288
X12 82 0.8963175 2.370362 0.6843049m.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")
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 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.
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 58.9 # 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.2 # 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 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:
\[\Large 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 normal\left(\mu_0,\sigma_0\right)\)
- \(n_0 = \frac{s}{\sigma_0}\)
Exemplo: MCT Feminino
\[\Large 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 normal\left(55,5\right)\)
- \(n_0 = \frac{9.1}{5} = 1.82\)
\[\Large IC^{0.95}_{\text{bayes}}(\mu)=\left[56.4,58.8\right]\]
\[\Large IC^{0.95}_{\text{freq}}(\mu)=\left[56.5,58.8\right]\]
\[\Large 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
“Point 11: Bayes Factors Often Agree With p Values […] Results between Bayes factors and classical tests often agree with each other.”
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/
Script 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
options(warn=-1)
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 <- readxl::read_excel(path="Biometria_FMUSP.xlsx",
sheet="dados",
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$AtivFisica <- factor(Dados$AtivFisica)
Dados$Sedentarismo <- factor(Dados$Sedentarismo)
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")
result <- MVN::mvn(data=subset(Dados,
select=c(Sexo, MCT)),
subset="Sexo",
mvnTest="hz",
univariateTest="SW")
result$univariateNormality
# 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)