MAN5728: Tópicos Essenciais da Bioestatística
Módulo 1: Estatística Descritiva
07 julho 2025 18:18h
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(finalfit, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(ggpubr, warn.conflicts=FALSE))
suppressMessages(library(HH, warn.conflicts=FALSE))Material
- HTML de R Markdown em
RPubs - Prof. José Siqueira: ResearchGate
Módulo 1: Parte 1
- Tipos de variáveis e apresentação dos dados.
- Estatística descritiva e apresentação dos resultados por meio de gráficos e tabelas.
- Medidas de tendência central e de dispersão.
Ler planilha
Os dados da planilha 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
- AtivFisica: nível de atividade física
- Sedentarismo: Não, Sim
- MCT: massa corporal total (kg)
- Estatura: cm
Visualização dos dados
# A tibble: 6 × 12
ID Ano Turma Sexo Mao TipoSang ABO AtivFisica Sedentarismo MCT
<dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
1 1 2 B M D A+ A baixa_intensi… N 68
2 2 2 B M D A+ A atualmente_in… S 82
3 3 1 B M D A+ A media_intensi… N 79
4 4 1 B F D A+ A media_intensi… N 49
5 5 1 B F D O- O atualmente_in… S 52
6 6 3 A M D A+ A alta_intensid… N 73
# ℹ 2 more variables: Estatura <dbl>, Rh <chr># A tibble: 6 × 12
ID Ano Turma Sexo Mao TipoSang ABO AtivFisica Sedentarismo MCT
<dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
1 544 3 B F C <NA> <NA> alta_intensid… N 54
2 545 3 B F D B+ B alta_intensid… N 47
3 546 3 B F D B+ B sempre_inativo S 48
4 547 3 B M A A+ A alta_intensid… N 70
5 548 3 B M D O- O atualmente_in… S 67
6 549 3 B F D A+ A media_intensi… N 45
# ℹ 2 more variables: Estatura <dbl>, Rh <chr>Dados
| ID | Ano | Turma | Sexo | Mao | TipoSang | ABO | AtivFisica | Sedentarismo | MCT | Estatura | Rh |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2 | B | M | D | A+ | A | baixa_intensidade | N | 68 | 173 | + |
| 2 | 2 | B | M | D | A+ | A | atualmente_inativo | S | 82 | 175 | + |
| 3 | 1 | B | M | D | A+ | A | media_intensidade | N | 79 | 172 | + |
| 4 | 1 | B | F | D | A+ | A | media_intensidade | N | 49 | 160 | + |
| 5 | 1 | B | F | D | O- | O | atualmente_inativo | S | 52 | 164 | - |
| 6 | 3 | A | M | D | A+ | A | alta_intensidade | N | 73 | 175 | + |
| 7 | 3 | A | M | D | A+ | A | alta_intensidade | N | 73 | 175 | + |
| 8 | 2 | B | M | D | B+ | B | media_intensidade | N | 60 | 173 | + |
| 9 | 3 | B | M | D | O+ | O | alta_intensidade | N | 75 | 182 | + |
| 10 | 3 | B | M | D | O+ | O | alta_intensidade | N | 75 | 182 | + |
Estrutura dos dados
tibble [549 × 12] (S3: tbl_df/tbl/data.frame)
$ ID : num [1:549] 1 2 3 4 5 6 7 8 9 10 ...
$ Ano : num [1:549] 2 2 1 1 1 3 3 2 3 3 ...
$ Turma : chr [1:549] "B" "B" "B" "B" ...
$ Sexo : chr [1:549] "M" "M" "M" "F" ...
$ Mao : chr [1:549] "D" "D" "D" "D" ...
$ TipoSang : chr [1:549] "A+" "A+" "A+" "A+" ...
$ ABO : chr [1:549] "A" "A" "A" "A" ...
$ AtivFisica : chr [1:549] "baixa_intensidade" "atualmente_inativo" "media_intensidade" "media_intensidade" ...
$ Sedentarismo: chr [1:549] "N" "S" "N" "N" ...
$ MCT : num [1:549] 68 82 79 49 52 73 73 60 75 75 ...
$ Estatura : num [1:549] 173 175 172 160 164 175 175 173 182 182 ...
$ Rh : chr [1:549] "+" "+" "+" "+" ...Converter variável alfanumérica para fator (nominal)
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,
levels=c("sempre_inativo",
"atualmente_inativo",
"baixa_intensidade",
"media_intensidade",
"alta_intensidade"),
labels=c("sempre_inativo",
"atualmente_inativo",
"baixa_intensidade",
"media_intensidade",
"alta_intensidade"))
Dados$Sedentarismo <- factor(Dados$Sedentarismo)
Dados$Rh <- factor(Dados$Rh)
str(Dados)tibble [549 × 12] (S3: tbl_df/tbl/data.frame)
$ ID : Factor w/ 549 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Ano : Factor w/ 3 levels "1","2","3": 2 2 1 1 1 3 3 2 3 3 ...
$ Turma : Factor w/ 2 levels "A","B": 2 2 2 2 2 1 1 2 2 2 ...
$ Sexo : Factor w/ 2 levels "F","M": 2 2 2 1 1 2 2 2 2 2 ...
$ Mao : Factor w/ 3 levels "A","C","D": 3 3 3 3 3 3 3 3 3 3 ...
$ TipoSang : Factor w/ 8 levels "A-","A+","AB-",..: 2 2 2 2 7 2 2 6 8 8 ...
$ ABO : Factor w/ 4 levels "A","AB","B","O": 1 1 1 1 4 1 1 3 4 4 ...
$ AtivFisica : Factor w/ 5 levels "sempre_inativo",..: 3 2 4 4 2 5 5 4 5 5 ...
$ Sedentarismo: Factor w/ 2 levels "N","S": 1 2 1 1 2 1 1 1 1 1 ...
$ MCT : num [1:549] 68 82 79 49 52 73 73 60 75 75 ...
$ Estatura : num [1:549] 173 175 172 160 164 175 175 173 182 182 ...
$ Rh : Factor w/ 2 levels "-","+": 2 2 2 2 1 2 2 2 2 2 ...summary
Ano Turma Sexo Mao TipoSang ABO
1:171 A:285 F:232 A : 28 A+ :171 A :220
2:155 B:264 M:317 C : 37 O+ :132 AB : 42
3:223 D :442 B+ : 65 B : 99
NA's: 42 A- : 49 O :177
O- : 45 NA's: 11
(Other): 76
NA's : 11
AtivFisica Sedentarismo MCT Estatura
sempre_inativo : 71 N:316 Min. : 41.00 Min. :120.0
atualmente_inativo:152 S:233 1st Qu.: 56.00 1st Qu.:164.0
baixa_intensidade : 59 Median : 64.00 Median :172.0
media_intensidade : 99 Mean : 66.76 Mean :170.8
alta_intensidade :168 3rd Qu.: 73.00 3rd Qu.:177.0
Max. :658.00 Max. :195.0
NA's :6 NA's :6
Rh
- :130
+ :408
NA's: 11
Vídeo: Clima “atípico” em Angola deixa menina do tempo preocupada… mínimas mais altas que as máximas!
Exclusão de dados de MCT e Estatura com erros tipográficos de MCT e estatura
Dados$MCT[Dados$MCT==max(Dados$MCT,na.rm=TRUE)] <- NA
Dados$Estatura[Dados$Estatura==min(Dados$Estatura,na.rm=TRUE)] <- NA
summary(subset(Dados, select=c(MCT, Estatura))) MCT Estatura
Min. : 41.00 Min. :153.0
1st Qu.: 56.00 1st Qu.:164.0
Median : 64.00 Median :172.0
Mean : 65.67 Mean :170.9
3rd Qu.: 73.00 3rd Qu.:177.0
Max. :125.00 Max. :195.0
NA's :7 NA's :7 Análise de dados faltantes
n.total <- nrow(Dados)
n.completo <- nrow(na.omit(Dados))
n.incompleto <- n.total - n.completo
cat("Numero de casos total = ", n.total, "\n", sep="")Numero de casos total = 549cat("Numero de casos completos = ", n.completo,
" (",round(100*n.completo/n.total,2),"%)\n", sep="")Numero de casos completos = 488 (88.89%)cat("Numero de casos incompletos = ", n.incompleto,
" (",round(100*n.incompleto/n.total,2),"%)\n", sep="")Numero de casos incompletos = 61 (11.11%)obs.falt <- sum(is.na(Dados))
obs.valid <- sum(!is.na(Dados))
obs.tot <- obs.falt + obs.valid
cat("Numero de observacoes validas = ", obs.valid,
" (",round(100*obs.valid/obs.tot,2),"%)\n", sep="")Numero de observacoes validas = 6499 (98.65%)cat("Numero de observacoes faltantes = ", obs.falt,
" (",round(100*obs.falt/obs.tot,2),"%)\n", sep="")Numero de observacoes faltantes = 89 (1.35%)Análise de dados faltantes com
finalfit::missing_pattern
Ano Turma Sexo AtivFisica Sedentarismo MCT Estatura TipoSang ABO Rh Mao
488 1 1 1 1 1 1 1 1 1 1 1 0
42 1 1 1 1 1 1 1 1 1 1 0 1
11 1 1 1 1 1 1 1 0 0 0 1 3
1 1 1 1 1 1 1 0 1 1 1 1 1
1 1 1 1 1 1 0 1 1 1 1 1 1
6 1 1 1 1 1 0 0 1 1 1 1 2
0 0 0 0 0 7 7 11 11 11 42 89psych::describe e table
# item name ,item number, nvalid, mean, sd,
# median, mad (Median Absolute Deviation),
# min, max, skew, kurtosis, se
print(psych::describe(subset(Dados, select=c(MCT, Estatura)))) vars n mean sd median trimmed mad min max range skew kurtosis
MCT 1 542 65.67 12.90 64 64.58 13.34 41 125 84 1.04 2.04
Estatura 2 542 170.94 8.92 172 170.80 10.38 153 195 42 0.12 -0.66
se
MCT 0.55
Estatura 0.38
A- A+ AB- AB+ B- B+ O- O+
49 171 2 40 34 65 45 132
A- A+ AB- AB+ B- B+ O- O+
0.09 0.32 0.00 0.07 0.06 0.12 0.08 0.25 tapply
$F
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
41.00 52.00 56.00 57.63 62.00 120.00 2
$M
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
43.00 63.00 72.00 71.59 78.00 125.00 5 psych::describeBy
print(psych::describeBy(subset(Dados,select=c(MCT, Estatura)),
group=list(Dados$Sexo),
mat=TRUE,
digits=2)) item group1 vars n mean sd median trimmed mad min max range
MCT1 1 F 1 230 57.63 9.05 56 56.75 7.41 41 120 79
MCT2 2 M 1 312 71.59 12.09 72 70.88 11.12 43 125 82
Estatura1 3 F 2 229 164.34 6.38 163 163.94 5.93 153 186 33
Estatura2 4 M 2 313 175.77 7.26 175 175.88 7.41 155 195 40
skew kurtosis se
MCT1 2.43 12.04 0.60
MCT2 0.90 2.37 0.68
Estatura1 0.61 0.02 0.42
Estatura2 -0.13 0.21 0.41print(psych::describeBy(subset(Dados,select=c(MCT, Estatura)),
group=list(Dados$Sexo, Dados$ABO),
mat=TRUE,
digits=2)) item group1 group2 vars n mean sd median trimmed mad min max
MCT1 1 F A 1 89 57.79 9.32 57.0 57.12 7.41 43 120
MCT2 2 M A 1 126 73.21 13.06 72.0 72.16 10.38 50 125
MCT3 3 F AB 1 14 56.93 5.69 56.5 57.00 6.67 48 65
MCT4 4 M AB 1 27 71.93 11.33 72.0 71.65 8.90 50 102
MCT5 5 F B 1 42 55.26 7.02 53.0 54.44 5.19 45 76
MCT6 6 M B 1 56 68.38 10.48 65.0 67.61 7.41 49 95
MCT7 7 F O 1 81 59.10 10.07 57.0 57.98 7.41 41 100
MCT8 8 M O 1 96 71.38 11.48 73.0 71.12 10.38 43 120
Estatura1 9 F A 2 88 164.59 5.83 164.0 164.39 5.93 153 181
Estatura2 10 M A 2 126 175.17 6.79 175.0 175.08 7.41 160 192
Estatura3 11 F AB 2 14 165.21 5.67 164.0 165.08 5.93 156 176
Estatura4 12 M AB 2 27 175.41 8.95 177.0 175.83 8.90 155 195
Estatura5 13 F B 2 42 163.43 6.84 162.5 162.76 6.67 153 185
Estatura6 14 M B 2 57 175.26 7.73 175.0 175.23 5.93 160 195
Estatura7 15 F O 2 81 164.33 6.90 163.0 163.80 7.41 154 186
Estatura8 16 M O 2 96 176.98 7.12 177.0 177.14 7.41 156 193
range skew kurtosis se
MCT1 77 3.29 19.99 0.99
MCT2 75 1.12 2.65 1.16
MCT3 17 0.06 -1.50 1.52
MCT4 52 0.47 0.64 2.18
MCT5 31 1.03 0.52 1.08
MCT6 46 0.69 0.08 1.40
MCT7 59 1.69 4.67 1.12
MCT8 77 0.58 2.26 1.17
Estatura1 28 0.39 -0.53 0.62
Estatura2 32 0.09 -0.35 0.61
Estatura3 20 0.25 -1.05 1.52
Estatura4 40 -0.35 0.23 1.72
Estatura5 32 0.91 0.72 1.06
Estatura6 35 0.04 -0.23 1.02
Estatura7 32 0.69 -0.06 0.77
Estatura8 37 -0.36 0.81 0.73print(psych::describeBy(subset(Dados, select=c(MCT, Estatura)),
group=list(Dados$ABO, Dados$Sexo),
mat=TRUE,
digits=2)) item group1 group2 vars n mean sd median trimmed mad min max
MCT1 1 A F 1 89 57.79 9.32 57.0 57.12 7.41 43 120
MCT2 2 AB F 1 14 56.93 5.69 56.5 57.00 6.67 48 65
MCT3 3 B F 1 42 55.26 7.02 53.0 54.44 5.19 45 76
MCT4 4 O F 1 81 59.10 10.07 57.0 57.98 7.41 41 100
MCT5 5 A M 1 126 73.21 13.06 72.0 72.16 10.38 50 125
MCT6 6 AB M 1 27 71.93 11.33 72.0 71.65 8.90 50 102
MCT7 7 B M 1 56 68.38 10.48 65.0 67.61 7.41 49 95
MCT8 8 O M 1 96 71.38 11.48 73.0 71.12 10.38 43 120
Estatura1 9 A F 2 88 164.59 5.83 164.0 164.39 5.93 153 181
Estatura2 10 AB F 2 14 165.21 5.67 164.0 165.08 5.93 156 176
Estatura3 11 B F 2 42 163.43 6.84 162.5 162.76 6.67 153 185
Estatura4 12 O F 2 81 164.33 6.90 163.0 163.80 7.41 154 186
Estatura5 13 A M 2 126 175.17 6.79 175.0 175.08 7.41 160 192
Estatura6 14 AB M 2 27 175.41 8.95 177.0 175.83 8.90 155 195
Estatura7 15 B M 2 57 175.26 7.73 175.0 175.23 5.93 160 195
Estatura8 16 O M 2 96 176.98 7.12 177.0 177.14 7.41 156 193
range skew kurtosis se
MCT1 77 3.29 19.99 0.99
MCT2 17 0.06 -1.50 1.52
MCT3 31 1.03 0.52 1.08
MCT4 59 1.69 4.67 1.12
MCT5 75 1.12 2.65 1.16
MCT6 52 0.47 0.64 2.18
MCT7 46 0.69 0.08 1.40
MCT8 77 0.58 2.26 1.17
Estatura1 28 0.39 -0.53 0.62
Estatura2 20 0.25 -1.05 1.52
Estatura3 32 0.91 0.72 1.06
Estatura4 32 0.69 -0.06 0.77
Estatura5 32 0.09 -0.35 0.61
Estatura6 40 -0.35 0.23 1.72
Estatura7 35 0.04 -0.23 1.02
Estatura8 37 -0.36 0.81 0.73aggregate
Group.1 MCT Estatura
1 F 57.63043 164.3406
2 M 71.58974 175.7732print(aggregate(subset(Dados, select=c(MCT, Estatura)),
by=list(Dados$Sexo, Dados$ABO),
FUN=mean,
na.rm=TRUE)) Group.1 Group.2 MCT Estatura
1 F A 57.78652 164.5909
2 M A 73.20635 175.1746
3 F AB 56.92857 165.2143
4 M AB 71.92593 175.4074
5 F B 55.26190 163.4286
6 M B 68.37500 175.2632
7 F O 59.09877 164.3333
8 M O 71.37500 176.9792print(aggregate(subset(Dados, select=c(MCT, Estatura)),
by=list(Dados$ABO, Dados$Sexo),
FUN=mean,
na.rm=TRUE)) Group.1 Group.2 MCT Estatura
1 A F 57.78652 164.5909
2 AB F 56.92857 165.2143
3 B F 55.26190 163.4286
4 O F 59.09877 164.3333
5 A M 73.20635 175.1746
6 AB M 71.92593 175.4074
7 B M 68.37500 175.2632
8 O M 71.37500 176.9792aggregate com fórmula
Sexo MCT
1 F 56
2 M 72 Sexo ABO MCT
1 F A 57.0
2 M A 72.0
3 F AB 56.5
4 M AB 72.0
5 F B 53.0
6 M B 65.0
7 F O 57.0
8 M O 73.0Média e desvio-padrão: Parte 1
- Higgins JPT, Li T, Deeks JJ (editors). Chapter 6: Choosing effect
measures and computing estimates of effect. In: Higgins JPT, Thomas J,
Chandler J, Cumpston M, Li T, Page MJ, Welch VA (editors). Cochrane
Handbook for Systematic Reviews of Interventions version 6.3 (updated
February 2022). Cochrane, 2022. Available from
www.training.cochrane.org/handbook.
- Table 6.5.a em https://training.cochrane.org/handbook/current/chapter-06
- Suponha que não temos acesso aos dados brutos de MCT. Os dados disponíveis são apresentados na tabela abaixo:
| group | n | mean | sd |
|---|---|---|---|
| F | 230 | 57.63 | 9.05 |
| M | 312 | 71.59 | 12.09 |
- Quais são os valores da média e desvio-padrão amostrais dos 542 estudantes?
Media = 65.67dp <- sqrt(((230-1)*(9.05^2) + (312-1)*(12.09^2) +
(230*312/542)*((57.63-71.59)^2))/(542-1))
cat("Desvio-padrao = ", round(dp,2), " > ", max(9.05,12.09), "\n", sep="")Desvio-padrao = 12.9 > 12.09Média e desvio-padrão: Parte 2
MCT
n mean sd 542 65.67 12.9
| group | n | mean | sd |
|---|---|---|---|
| F | 230 | 57.63 | 9.05 |
| M | 312 | 71.59 | 12.09 |
Media = 65.67dp <- sqrt(((230-1)*9.05^2 + (312-1)*12.09^2 +
(230*312/542)*(57.63-71.59)^2)/(542-1))
cat("Desvio-padrao = ", round(dp,2), "\n", sep="")Desvio-padrao = 12.9dp <- sqrt((230*9.05^2 + 312*12.09^2 +
230*(57.63-65.67)^2 + 312*(71.59-65.67)^2)/542)
cat("Desvio-padrao = ", round(dp,2), "\n", sep="")Desvio-padrao = 12.9\[\bar{x} = \dfrac{n_F \, \bar{x}_F + n_M \, \bar{x}_M}{n_F + n_M}\]
\[\begin{align} s&=\sqrt{\dfrac{(n_F-1)\,s_F^2+(n_M-1)\,s_M^2+\dfrac{n_F\;n_M}{n_F+n_M}\left(\bar{x}_F-\bar{x}_M\right)^2}{n_F+n_M-1}}\\ s&=\sqrt{\dfrac{n_F\,s_F^2+n_M\,s_M^2+n_F\,\left(\bar{x}_F-\bar{x}\right)^2+n_M\,\left(\bar{x}_M-\bar{x}\right)^2}{n_F+n_M}}\\ s&=\sqrt{\dfrac{n_F\,\left(s_F^2+\left(\bar{x}_F-\bar{x}\right)^2\right)+n_M\,\left(s_M^2+\left(\bar{x}_M-\bar{x}\right)^2\right)}{n_F+n_M}} \end{align}\]
“Dotplot” com plot & table
Gráfico de setores
F M
232 317 
A AB B O
220 42 99 177 
Tabela de contingência 2x2
Sexo
Sedentarismo F M
N 133 183
S 99 134Sedentarismo
N S
316 233 Sexo
F M
232 317 Sexo
Sedentarismo F M
N 0.24 0.33
S 0.18 0.24 Sexo
Sedentarismo F M
N 0.42 0.58
S 0.42 0.58 Sexo
Sedentarismo F M
N 0.57 0.58
S 0.43 0.42
barplot(Sedentarismo.Sexo,
beside=TRUE,
legend.text=rownames(Sedentarismo.Sexo),
ylab="Freq",
xlab="Sexo x Sedentarismo")
barplot(proportions(Sedentarismo.Sexo),
beside=TRUE,
legend.text=rownames(Sedentarismo.Sexo),
ylab="Freq",
xlab="Sexo x Sedentarismo")
sexo.ABO.freq <- as.data.frame(table(Dados$Sexo, Dados$ABO))
names(sexo.ABO.freq) <- c("Sexo", "ABO", "Freq")
ggpubr::ggbarplot(sexo.ABO.freq,
x="ABO",
y="Freq",
color="Sexo",
palette=c("gray", "black"),
order=c("A", "B", "AB", "O"),
width=.7)
ggpubr::ggbarplot(sexo.ABO.freq,
x="ABO",
y="Freq",
color="Sexo",
palette=c("gray", "black"),
order=c("A", "B", "AB", "O"),
position = ggplot2::position_dodge(),
width=.7)
Tabela de contingência 2x4
Sexo
ABO F M
A 91 129
AB 14 28
B 42 57
O 81 96 Sexo
ABO F M
A 0.17 0.24
AB 0.03 0.05
B 0.08 0.11
O 0.15 0.18 Sexo
ABO F M
A 0.41 0.59
AB 0.33 0.67
B 0.42 0.58
O 0.46 0.54 Sexo
ABO F M
A 0.40 0.42
AB 0.06 0.09
B 0.18 0.18
O 0.36 0.31
barplot(proportions(ABO.Sexo),
beside=TRUE,
legend.text=rownames(ABO.Sexo),
ylab="Freq",
xlab="Sexo x ABO")
Tabela multivariada: mais de duas variáveis categóricas
, , ABO = A
Sedentarismo
Sexo N S
F 51 40
M 83 46
, , ABO = AB
Sedentarismo
Sexo N S
F 9 5
M 12 16
, , ABO = B
Sedentarismo
Sexo N S
F 26 16
M 27 30
, , ABO = O
Sedentarismo
Sexo N S
F 45 36
M 56 40 Sexo F M
Sedentarismo N S N S
ABO
A 51 40 83 46
AB 9 5 12 16
B 26 16 27 30
O 45 36 56 40
Boxplot
ggpubr::ggboxplot(data=Dados,
x="Sexo",
y="MCT",
add="",
orientation="horizontal",
width=.7,
order=c("F", "M"))
ggpubr::ggboxplot(data=Dados,
x="Sexo",
y="MCT",
add="",
orientation="horizontal",
width=.7,
order=c("M", "F"))
ggpubr::ggboxplot(data=Dados,
x="Sexo",
y="MCT",
add="jitter",
orientation="horizontal",
width=.7,
order=c("F", "M"))
ggpubr::ggboxplot(data=Dados,
x="ABO",
y="MCT",
add="",
orientation="horizontal",
width=.7,
select=c("A", "B", "O"),
order=c("O", "B", "A"))
Bagplot: boxplot bidimensional
bgp.F <- DescTools::PlotBag(Dados.F$Estatura,
Dados.F$MCT,
main=paste("Feminino"),
xlab="Estatura (cm)",
ylab="Massa Corporal Total (kg)",
na.rm = TRUE,
show.bagpoints=FALSE,
show.looppoints=FALSE,
show.whiskers=FALSE,
col.loophull = "white",
col.looppoints = "black",
col.baghull = "white",
col.bagpoints = "black",
cex=1)[1] "Warning: NA elements have been removed!!" x y
1 165 100
2 165 100
3 158 78for (o in 1:nrow(outliers.F))
{
r.F <- which(Dados.F$Estatura==outliers.F$x[o] &
Dados.F$MCT==outliers.F$y[o])
text(outliers.F$x[o],outliers.F$y[o], r.F, pos=1, cex=0.7)
}
bgp.M <- DescTools::PlotBag(Dados.M$Estatura,
Dados.M$MCT,
main=paste("Masculino"),
xlab="Estatura (cm)",
ylab="Massa Corporal Total (kg)",
na.rm = TRUE,
show.bagpoints=FALSE,
show.looppoints=FALSE,
show.whiskers=FALSE,
col.loophull = "white",
col.looppoints = "black",
col.baghull = "white",
col.bagpoints = "black",
cex=1)[1] "Warning: NA elements have been removed!!" x y
1 165 94
2 193 120
3 184 125
4 176 125for (o in 1:nrow(outliers.M))
{
r.M <- which(Dados.M$Estatura==outliers.M$x[o] &
Dados.M$MCT==outliers.M$y[o])
text(outliers.M$x[o], outliers.M$y[o], r.M, pos=1, cex=0.7)
}
Histograma
Um histograma é uma representação gráfica da distribuição de um conjunto de dados. Ele é utilizado para mostrar a frequência com que diferentes intervalos de valores ocorrem em um conjunto de dados. A estrutura básica de um histograma inclui:
- Eixo X (Horizontal): Este eixo representa os intervalos (ou “bins”) nos quais os dados são agrupados. Cada intervalo abrange uma faixa específica de valores.
- Eixo Y (Vertical): Este eixo mostra a frequência (ou contagem) de dados que se enquadram em cada intervalo.
Para construir um histograma: - Primeiro, os dados são divididos em intervalos não sobrepostos. - Em seguida, contam-se quantos dados caem em cada intervalo. - Finalmente, retângulos (ou barras) são desenhados para cada intervalo, onde a altura de cada barra é proporcional à frequência de dados naquele intervalo.
Silveira, PSP & Siqueira, JO (2022) Histogram lies about distribution shape and Pearson’s coefficient of variation lies about relative variability. The Quantitative Methods for Psychology 18(1). DOI 10.20982/tqmp.18.1.p091
How to Lie with Histograms
Gráfico de densidade
Gráfico de distribuição de probabilidade: ECDF plot
Gráfico de dispersão com sunflowerplot
sunflowerplot(MCT~Estatura,
data=Dados.F,
rotate=TRUE,
pch=1,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8)
sunflowerplot(MCT~Estatura,
data=Dados.M,
rotate=TRUE,
pch=2,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8,
add=FALSE)
sunflowerplot(MCT~Estatura,
data=Dados.F,
rotate=TRUE,
pch=1,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8)
sunflowerplot(MCT~Estatura,
data=Dados.M,
rotate=TRUE,
pch=2,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8,
add=TRUE)
Gráfico de dispersão com car::scatterplot
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=FALSE,
smooth=FALSE,
ellipse=FALSE,
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=FALSE,
smooth=FALSE,
boxplots=TRUE,
ellipse=list(levels=c(0.68),
robust=TRUE,
fill=FALSE,
fill.alpha=0.2),
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=FALSE,
smooth=FALSE,
ellipse=list(levels=c(0.68,0.95),
robust=TRUE,
fill=FALSE,
fill.alpha=0.2),
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=TRUE,
smooth=FALSE,
ellipse=FALSE,
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
Gráfico matricial
car::scatterplotMatrix(Dados[,c("Estatura",
"MCT")],
groups=Dados$Sexo,
regLine=TRUE,
smooth=FALSE,
boxplots=TRUE,
by.groups=TRUE,
ellipse=list(levels=c(0.5),
robust=TRUE,
fill=FALSE),
grid=FALSE,
col=c("#666666","#888888","#cccccc"),
cex=0.5,
cex.labels=1,
row1attop=TRUE)
GGally::ggpairs(subset(Dados,
select=-c(ID,Ano,Turma,Mao,TipoSang,
ABO,AtivFisica)),
ggplot2::aes(colour=Sexo))`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Gráficos avançados
- ggplot2
- Módulo V: Gráficos com ggplot2 por Bruno Lucian & Claudia Jakelline
- ggplot2 - Essentials
- The R Graph Gallery
- Advanced plotting by Jeff Powell and Remko Duursma
- Chapter 19. Advanced graphics with ggplot2
- Top 50 ggplot2 Visualizations - The Master List (With Full R Code)
Teste de distribuição normal
result <- try(MVN::mvn(data=subset(Dados,
select=c(Sexo, MCT, Estatura)),
subset="Sexo",
mvn_test = "hz",
univariate_test = "SW",
multivariate_outlier_method ="quan"))
print(result$multivariate_normality) Group Test Statistic p.value MVN
1 F Henze-Zirkler 3.507 <0.001 ✗ Not normal
2 M Henze-Zirkler 1.365 0.008 ✗ Not normal Group Test Variable Statistic p.value Normality
1 F Shapiro-Wilk MCT 0.905 <0.001 ✗ Not normal
2 F Shapiro-Wilk Estatura 0.963 <0.001 ✗ Not normal
3 M Shapiro-Wilk MCT 0.954 <0.001 ✗ Not normal
4 M Shapiro-Wilk Estatura 0.991 0.045 ✗ Not normal Group Observation Mahalanobis.Distance
1 F 5 82.681
2 F 6 82.681
3 F 128 29.788
4 F 222 20.620
5 F 223 20.620
6 F 7 18.758
7 F 4 17.027
8 F 98 15.497
9 F 71 12.204
10 F 28 11.319
11 F 135 10.173
12 F 30 9.369
13 F 31 9.369
14 F 139 9.061
15 F 106 7.966
16 F 80 7.965
17 M 283 41.064
18 M 266 31.922
19 M 44 22.582
20 M 13 18.941
21 M 24 16.292
22 M 25 16.292
23 M 153 13.602
24 M 100 11.649
25 M 224 11.589
26 M 165 11.374
27 M 166 11.374
28 M 256 10.317
29 M 257 10.317
30 M 148 10.241
31 M 259 9.553
32 M 35 8.815
33 M 290 8.609
34 M 222 8.494
35 M 54 8.387
36 M 142 8.293
37 M 104 8.081Código R completo
library(readxl)
library(finalfit)
library(car)
library(GGally)
library(ggplot2)
library(psych)
library(ggpubr)
# http://www.hiercourse.com/docs/advanced_plotting.html
# http://ecor.ib.usp.br/doku.php?id=03_apostila:start
# https://docs.ufpr.br/~aanjos/CE231/web/apostila.html#Q1-1-2
# https://pt.wikipedia.org/wiki/Distribui%C3%A7%C3%A3o_do_tipo_sangu%C3%ADneo_por_pa%C3%ADs
# Descriptive Statistics and Graphics
## http://www.sthda.com/english/wiki/descriptive-statistics-and-graphics
# Plotting means and error bars (ggplot2)
## http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
Dados <- readxl::read_excel(path="Biometria_FMUSP.xlsx",
sheet="dados",
na="NA")
# Preparacao dos dados
head(Dados)
tail(Dados)
print(str(Dados))
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)
print(str(Dados))
print(summary(subset(Dados, select=-ID)))
Dados$MCT[Dados$MCT==max(Dados$MCT,na.rm=TRUE)] <- NA
Dados$Estatura[Dados$Estatura==min(Dados$Estatura,na.rm=TRUE)] <- NA
print(summary(subset(Dados, select=c(MCT, Estatura))))
Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")
# Missing Value Analysis
print(sapply(Dados,function(x){sum(!is.na(x))}))
print(sapply(Dados,function(x){sum(is.na(x))}))
n.total <- nrow(Dados)
n.completo <- nrow(na.omit(Dados))
n.incompleto <- n.total - n.completo
cat("Numero de casos total = ", n.total, "\n", sep="")
cat("Numero de casos completos = ", n.completo,
" (",round(100*n.completo/n.total,2),"%)\n", sep="")
cat("Numero de casos incompletos = ", n.incompleto,
" (",round(100*n.incompleto/n.total,2),"%)\n", sep="")
obs.falt <- sum(is.na(Dados))
obs.valid <- sum(!is.na(Dados))
obs.tot <- obs.falt + obs.valid
cat("Numero de observacoes validas = ", obs.valid,
" (",round(100*obs.valid/obs.tot,2),"%)\n", sep="")
cat("Numero de observacoes faltantes = ", obs.falt,
" (",round(100*obs.falt/obs.tot,2),"%)\n", sep="")
print(finalfit::missing_pattern(subset(Dados,select=-ID)))
# Analise descritiva numerica
# mean,median,25th and 75th quartiles,min,max
tapply(Dados$MCT,Dados$Sexo,summary)
# Tukey min,lower-hinge, median,upper-hinge,max
tapply(Dados$MCT,Dados$Sexo,fivenum)
# item name ,item number, nvalid, mean, sd,
# median, mad (Median Absolute Deviation),
# min, max, skew, kurtosis, se
psych::describe(subset(Dados, select=c(MCT, Estatura)))
psych::describeBy(subset(Dados,select=c(MCT, Estatura)),
group=list(Dados$Sexo),
mat=TRUE,
digits=2)
psych::describeBy(subset(Dados,select=c(MCT, Estatura)),
group=list(Dados$Sexo, Dados$ABO),
mat=TRUE,
digits=2)
psych::describeBy(subset(Dados, select=c(MCT, Estatura)),
group=list(Dados$ABO, Dados$Sexo),
mat=TRUE,
digits=2)
agreg <- aggregate(subset(Dados, select=c(MCT, Estatura)),
by=list(Dados$Sexo),
FUN=mean,
na.rm=TRUE)
cat("\nMean\n")
print(agreg)
agreg <- aggregate(subset(Dados, select=c(MCT, Estatura)),
by=list(Dados$Sexo, Dados$ABO),
FUN=mean,
na.rm=TRUE)
cat("\nMean\n")
print(agreg)
agreg <- aggregate(subset(Dados, select=c(MCT, Estatura)),
by=list(Dados$ABO, Dados$Sexo),
FUN=mean,
na.rm=TRUE)
cat("\nMean\n")
print(agreg)
agreg <- aggregate(MCT ~ Sexo,
FUN=mean,
na.rm=TRUE,
data=Dados)
cat("\nMean\n")
print(agreg)
agreg <- aggregate(MCT ~ Sexo + ABO,
FUN=mean,
na.rm=TRUE,
data=Dados)
cat("\nMean\n")
print(agreg)
round(table(Dados$TipoSang)/sum(table(Dados$TipoSang),
na.rm=TRUE),2)
# Analise descritiva grafica
# Frequencia de cada observacao
plot(Dados$Sexo, xlab="Sexo", ylab="Freq")
plot(Dados$TipoSang, xlab="Tipo Sanguineo", ylab="Freq")
plot(Dados$Sedentarismo, xlab="Sedentarismo", ylab="Freq")
plot(table(Dados$MCT), xlab="Massa corporal total (kg)", ylab="Freq")
plot(table(Dados$Estatura), xlab="Estatura (cm)", ylab="Freq")
# Grafico de setores
table(Dados$Sexo)
pie(table(Dados$Sexo),
xlab="Sexo")
table(Dados$ABO)
pie(table(Dados$ABO),
xlab="ABO")
# Tabela de contingencia
print(Sedentarismo.Sexo <- xtabs(~Sedentarismo+Sexo, data=Dados))
margin.table(Sedentarismo.Sexo,1)
margin.table(Sedentarismo.Sexo,2)
round(proportions(Sedentarismo.Sexo),4)
round(proportions(Sedentarismo.Sexo,1),4)
round(proportions(Sedentarismo.Sexo,2),4)
plot(Dados$Sedentarismo~Dados$Sexo, xlab="Sexo", ylab="Sedentarismo")
mosaicplot(~Sexo+Sedentarismo, data=Dados, color=FALSE)
barplot(Sedentarismo.Sexo,
beside=TRUE,
legend.text=rownames(Sedentarismo.Sexo),
ylab="Freq",
xlab="Sexo x Sedentarismo")
barplot(proportions(Sedentarismo.Sexo),
beside=TRUE,
legend.text=rownames(Sedentarismo.Sexo),
ylab="Freq",
xlab="Sexo x Sedentarismo")
print(ABO.Sexo <- xtabs(~ABO+Sexo, data=Dados))
round(proportions(ABO.Sexo),4)
round(proportions(ABO.Sexo,1),4)
round(proportions(ABO.Sexo,2),4)
mosaicplot(~Sexo+ABO, data=Dados, color=FALSE)
barplot(ABO.Sexo,
beside=TRUE,
legend.text=rownames(ABO.Sexo),
ylab="Freq",
xlab="Sexo x ABO")
barplot(proportions(ABO.Sexo),
beside=TRUE,
legend.text=rownames(ABO.Sexo),
ylab="Freq",
xlab="Sexo x ABO")
sexo.ABO.freq <- as.data.frame(table(Dados$Sexo, Dados$ABO))
names(sexo.ABO.freq) <- c("Sexo", "ABO", "Freq")
ggpubr::ggbarplot(sexo.ABO.freq,
x="ABO",
y="Freq",
color="Sexo",
palette=c("gray", "black"),
order=c("A", "B", "AB", "O"),
width=.7)
ggpubr::ggbarplot(sexo.ABO.freq,
x="ABO",
y="Freq",
color="Sexo",
palette=c("gray", "black"),
order=c("A", "B", "AB", "O"),
position = ggplot2::position_dodge(),
width=.7)
# mean +- se
ggpubr::ggbarplot(data=Dados,
x="Sexo",
y="MCT",
add="mean_se",
error.plot="errorbar",
width=.5,
order=c("F", "M"),
ylim=c(55,75),
main="mean +- se")
ggpubr::ggbarplot(data=Dados,
x="Sexo",
y="MCT",
add="mean_se",
error.plot="errorbar",
width=.5,
order=c("M", "F"),
ylim=c(55,75),
main="mean +- se")
MCT.ci <- tapply(Dados$MCT,
Dados$Sexo,
gmodels::ci,
na.rm=TRUE)
print(MCT.ci)
# mean +- sd
ggpubr::ggbarplot(data=Dados,
x="Sexo",
y="MCT",
add="mean_sd",
error.plot="errorbar",
width=.5,
order=c("M", "F"),
ylim=c(45,85),
main="mean +- sd")
MCT.pi <- tapply(Dados$MCT,
Dados$Sexo,
sd,
na.rm=TRUE)
cat("Limite inferior do intervalo de predicao de MCT Feminino = ",
mean(Dados.F$MCT,na.rm=TRUE) - as.numeric(MCT.pi[1]), "\n", sep="")
cat("Limite superior do intervalo de predicao de MCT Feminino = ",
mean(Dados.F$MCT,na.rm=TRUE) + as.numeric(MCT.pi[1]), "\n", sep="")
cat("Limite inferior do intervalo de predicao de MCT Masculino = ",
mean(Dados.M$MCT,na.rm=TRUE) - as.numeric(MCT.pi[2]), "\n", sep="")
cat("Limite superior do intervalo de predicao de MCT Masculino = ",
mean(Dados.M$MCT,na.rm=TRUE) + as.numeric(MCT.pi[2]), "\n", sep="")
# Multiway tables: More than two categorical variables
print(xtabs(~Sexo + Sedentarismo + ABO, data=Dados))
mosaicplot(xtabs(~Sexo + Sedentarismo + ABO, data=Dados))
ftable(Sexo + Sedentarismo ~ ABO, data=Dados)
# boxplot
boxplot(Dados$MCT, horizontal=TRUE,
xlab="MCT (kg)")
rug(jitter(Dados$MCT))
boxplot(MCT~Sexo, data=Dados, horizontal=TRUE,
xlab="MCT (kg)")
boxplot(MCT~Sexo+Sedentarismo, data=Dados, horizontal=TRUE,
xlab="MCT (kg)")
boxplot(MCT~Sexo+Sedentarismo+Mao, data=Dados, horizontal=TRUE,
xlab="MCT (kg)", cex=0.6)
ggpubr::ggboxplot(data=Dados,
y="MCT")
ggpubr::ggboxplot(data=Dados,
x="Sexo",
y="MCT",
add="",
orientation="horizontal",
width=.7,
order=c("F", "M"))
ggpubr::ggboxplot(data=Dados,
x="Sexo",
y="MCT",
add="",
orientation="horizontal",
width=.7,
order=c("M", "F"))
ggpubr::ggboxplot(data=Dados,
x="Sexo",
y="MCT",
add="jitter",
orientation="horizontal",
width=.7,
order=c("F", "M"))
ggpubr::ggboxplot(data=Dados,
x="ABO",
y="MCT",
add="",
orientation="horizontal",
width=.7,
select=c("A", "B", "O"),
order=c("O", "B", "A"))
# ECDF plot
ggpubr::ggecdf(data=Dados,
x="MCT",
linetype="Sexo")
# histograma
hist(Dados$MCT)
rug(jitter(Dados$MCT))
# bagplot
bgp.F <- DescTools::PlotBag(Dados.F$Estatura,
Dados.F$MCT,
main=paste("Feminino"),
xlab="Estatura (cm)",
ylab="Massa Corporal Total (kg)",
na.rm = TRUE,
show.bagpoints=FALSE,
show.looppoints=FALSE,
show.whiskers=FALSE,
col.loophull = "white",
col.looppoints = "black",
col.baghull = "white",
col.bagpoints = "black",
cex=1)
print(outliers.F <- as.data.frame(bgp.F$pxy.outlier))
for (o in 1:nrow(outliers.F))
{
r.F <- which(Dados.F$Estatura==outliers.F$x[o] &
Dados.F$MCT==outliers.F$y[o])
text(outliers.F$x[o],outliers.F$y[o], r.F, pos=1, cex=0.7)
}
bgp.M <- DescTools::PlotBag(Dados.M$Estatura,
Dados.M$MCT,
main=paste("Masculino"),
xlab="Estatura (cm)",
ylab="Massa Corporal Total (kg)",
na.rm = TRUE,
show.bagpoints=FALSE,
show.looppoints=FALSE,
show.whiskers=FALSE,
col.loophull = "white",
col.looppoints = "black",
col.baghull = "white",
col.bagpoints = "black",
cex=1)
print(outliers.M <- as.data.frame(bgp.M$pxy.outlier))
for (o in 1:nrow(outliers.M))
{
r.M <- which(Dados.M$Estatura==outliers.M$x[o] &
Dados.M$MCT==outliers.M$y[o])
text(outliers.M$x[o], outliers.M$y[o], r.M, pos=1, cex=0.7)
}
# Comparacao de medias
source("From_RcmdrMisc_plotMeans.R")
From_RcmdrMisc_plotMeans(Dados$MCT,
Dados$Sexo,
col="black",
connect=TRUE)
xtabs(~Sexo+Ano, data=Dados)
round(proportions(xtabs(~Sexo+Ano, data=Dados),2),4)
From_RcmdrMisc_plotMeans(as.numeric(Dados$Sexo)-1,
factor(Dados$Ano),
col="black",
connect=FALSE,
ylab="Proporcao de masculino")
From_RcmdrMisc_plotMeans(Dados$MCT,
Dados$Sedentarismo,
col="black",
connect=FALSE)
From_RcmdrMisc_plotMeans(Dados$MCT,
Dados$Sexo,
Dados$Sedentarismo,
col="black",
connect=FALSE)
str(Dados.F$Sedentarismo)
cor.test(Dados.F$MCT, as.numeric(Dados.F$Sedentarismo),
method = "spearman")
cor.test(Dados.M$MCT, as.numeric(Dados.M$Sedentarismo),
method = "spearman")
# Comparacao de distribuicoes
car::densityPlot(~MCT, data=Dados)
car::densityPlot(MCT~Sexo, data=Dados,
col=c("black","black"))
car::densityPlot(MCT~Sedentarismo, data=Dados.F,
col=c("black","black"),
main="Feminino")
car::densityPlot(MCT~Sedentarismo, data=Dados.M,
col=c("black","black"),
main="Masculino")
car::densityPlot(MCT~ABO, data=Dados.F,
main="Feminino")
car::densityPlot(MCT~ABO, data=Dados.M,
main="Masculino")
# Grafico de dispersao
sunflowerplot(MCT~Estatura,
data=Dados.F,
rotate=TRUE,
pch=1,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8)
sunflowerplot(MCT~Estatura,
data=Dados.M,
rotate=TRUE,
pch=2,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8,
add=FALSE)
sunflowerplot(MCT~Estatura,
data=Dados.F,
rotate=TRUE,
pch=1,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8)
sunflowerplot(MCT~Estatura,
data=Dados.M,
rotate=TRUE,
pch=2,
size=.1,
col="black",
seg.col="black",
seg.lwd=.8,
add=TRUE)
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=FALSE,
smooth=FALSE,
ellipse=FALSE,
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=FALSE,
smooth=FALSE,
boxplots=TRUE,
ellipse=list(levels=c(0.68),
robust=TRUE,
fill=FALSE,
fill.alpha=0.2),
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=FALSE,
smooth=FALSE,
ellipse=list(levels=c(0.68,0.95),
robust=TRUE,
fill=FALSE,
fill.alpha=0.2),
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
car::scatterplot(MCT~Estatura,
group=Dados$Sexo,
regLine=TRUE,
smooth=FALSE,
ellipse=FALSE,
grid=FALSE,
col="black",
xlim=c(145,200),
ylim=c(30,130),
data=Dados)
# Grafico matricial
car::scatterplotMatrix(Dados[,c("Estatura",
"MCT")],
groups=Dados$Sexo,
regLine=TRUE,
smooth=FALSE,
boxplots=TRUE,
by.groups=TRUE,
ellipse=list(levels=c(0.5),
robust=TRUE,
fill=FALSE),
grid=FALSE,
col=c("#666666","#888888","#cccccc"),
cex=0.5,
cex.labels=1,
row1attop=TRUE)
GGally::ggpairs(subset(Dados,
select=-c(ID,Ano,Turma,Mao,TipoSang,
ABO,AtivFisica,IMC)),
ggplot2::aes(colour=Sexo))
# Teste de normalidade
result <- MVN::mvn(data=subset(Dados,
select=c(Sexo, MCT, Estatura)),
subset="Sexo",
mvnTest="hz",
univariateTest="SW")
print(result$multivariateNormality)
print(result$univariateNormality)
# Binomial
k <- 0:10
plot(k,dbinom(x=k,size=10,prob=.45),
type='h',
main='Distribuicao binomial(n=10, p=0.3)',
ylab='Probabilidade',
xlab ='Numero de sucessos',
lwd=3)
# Poisson
k <- 0:10
plot(k,dpois(x=k,lambda=2),
type='h',
main='Distribuicao de Poisson(lambda=2)',
ylab='Probabilidade',
xlab ='x',
lwd=3)
# Normal
DescTools::PlotProbDist(breaks=c(-10, -1, 12),
function(x) dnorm(x, mean=2, sd=2),
blab=c(-1),
alab=NA,
xlim=c(-7,10),
main="normal(2,2)",
col=c("black", "black"),
density=c(20, 0))
p <- pnorm(q=140, mean=120, sd=10, lower.tail=FALSE)
DescTools::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))
p <- pnorm(q=140, mean=120, sd=10)-pnorm(q=100, mean=120, sd=10)
DescTools::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))
q1 <- round(qnorm(p=0.025, mean=120, sd=10),2)
q2 <- round(qnorm(p=0.975, 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("",0.95,""),
xlim=c(80,160),
main="normal(120,10)",
col=c("black", "black"),
density=c(0,20,0))
# plot t-distribution
DescTools::PlotProbDist(breaks=c(-6, -2.3, 1.5, 6),
function(x) dt(x, df=8),
blab=c(-2.3, 1.5),
xlim=c(-4,4),
alab=NA,
col=c("black", "black"),
main="t(8)",
density=c(10, 0))
# same for Chi-square
DescTools::PlotProbDist(breaks=c(0, 15, 35),
function(x) dchisq(x, df=8),
blab=c(15),
alab=NA,
xlim=c(0, 30),
main=expression(paste(chi^2,"(8)")),
col=c("black", "black"),
density=c(0, 20))Módulo 1: Parte 2
- Descrição dos princípios básicos de probabilidade e as formas mais comuns de distribuição de probabilidade.
- Conceitos de distribuição normal, [Poisson e] binomial.
Distribuição de probabilidade
- Distribuição de probabilidade é uma maneira estatística de expressar a incerteza dos eventos de uma variável.
- A variável pode ter dois ou mais eventos mutuamente excludentes.
- O conjunto de eventos pode ser discreto ou contínuo, finito ou infinito.
- As probabilidades dos eventos variam de 0 (exclusive) a 1 (exclusive).
- A soma das probabilidades dos eventos é igual a 1.
- E.g.: Variável discreta finita
- Sexo de recém-nascido: \(P(Masc) = 0.51\) e \(P(Fem) = 0.49\).
- E.g.: Variável discreta infinita
- Número de recém-nascidos (vivos) por dia no estado de São Paulo: \(média = 2000/dia\).
- E.g.: Variável contínua
- MCT de recém-nascido: mínimo = 1300 g, máximo = 4300 g, média = 3200 g e desvio-padrão = 400 g.
Distribuição binomial
- Distribuição de probabilidade de variável discreta finita \(X\)
- \(n=1,2,3,...\): número de realizações independentes do experimento
- Cada repetição do experimento produz um de dois resultados possíveis, denominados genericamente de sucesso (S) e fracasso (F)
- A probabilidade de sucesso, \(P(S) = p\), é constante em cada repetição.
- A probabilidade de fracasso, \(P(F) = 1 – p\), é constante em cada repetição.
- \(X \sim binomial(n, p)\) : variável aleatória \(X\) com distribuição binomial representando o número de sucessos em \(n\) repetições do experimento; \(X=0,...,n\)
- Média de \(X\): \(n p\)
- Desvio-padrão de \(X\): \(\sqrt{n p (1 – p)}\)
k <- 0:10
plot(k,dbinom(x=k,size=10,prob=.45),
type='h',
main='Distribuicao binomial(n=10, p=0.3)',
ylab='Probabilidade',
xlab ='Numero de sucessos',
lwd=3)
Exemplo de distribuição binomial
Fator Rh
“Além da classificação comum de sangue nos grupos A, B, AB e O, é importante a subdivisão de acordo com uma substância chamada fator Rhesus (Rh), que pode ser positivo (Rh+) ou negativo (Rh-). Podem ocorrer reações de incompatibilidade em transfusões de sangue. Por exemplo, um indivíduo Rh- só deve receber transfusão de sangue Rh-. Caso receba sangue Rh+, haverá sua sensibilização e a formação de anticorpos anti-Rh. Além disso, pode acontecer incompatibilidade Rh entre o sangue materno e o fetal. Aproximadamente 85% da população são Rh+, os outros 15% são Rh-.’’
Siqueira & Tibúrcio, 2011, p. 198-9
Siqueira, AL & Tibúrcio, JD (2011) Estatística na Área de Saúde: conceitos, metodologia, aplicações e prática computacional. BH: Coopmed.
- \(n = 3\) é o número de pacientes que realizarão um transplante.
- \(X\) é o número de pacientes com Rh- que realizarão um transplante (sucesso), variando de 0 e 3.
- \(P(Rh^-) = 0.15\)
- \(X \sim binomial(3, 0.15)\)
X ~ binomial (3, 0.15)
Media de X = 0.45
Desvio-padrao de X = 0.6184658
P(X = 0) = 0.614125
P(X = 1) = 0.325125
P(X = 2) = 0.057375
P(X = 3) = 0.003375
P(X = 0) + P(X = 1) + P(X = 2) + P(X = 3) = 1
P(X = 0) + P(X = 1) = 0.93925
P(X = 2) + P(X = 3) = 0.06075 
sink("exemplo01_binomial.txt")
n <- 3
p <- 0.15
cat("X ~ binomial (",n,", ",p,")\n",sep="")
m <- n*p
dp <- sqrt(n*p*(1-p))
cat("Media de X =",m,"\n")
cat("Desvio-padrao de X =",dp,"\n\n")
px0 <- dbinom(x=0,size=n,prob=p)
px1 <- dbinom(x=1,size=n,prob=p)
px2 <- dbinom(x=2,size=n,prob=p)
px3 <- dbinom(x=3,size=n,prob=p)
cat("P(X = 0) =",px0,"\n")
cat("P(X = 1) =",px1,"\n")
cat("P(X = 2) =",px2,"\n")
cat("P(X = 3) =",px3,"\n")
cat("\nP(X = 0) + P(X = 1) + P(X = 2) + P(X = 3) =",
px0+px1+px2+px3, "\n")
# lower.tail: logical; if TRUE (default), probabilities are P[X <= x],
# otherwise, P[X > x]
px01 <- pbinom(q=1,size=n,prob=p,lower.tail=TRUE)
cat("P(X = 0) + P(X = 1) =", px01, "\n")
px23 <- pbinom(q=1,size=n,prob=p,lower.tail=FALSE)
cat("P(X = 2) + P(X = 3) =", px23, "\n")
sink()
png("exemplo01_binomial.png")
X <- 0:n
probs <- dbinom(x=X,size=n,prob=p)
plot(X,
probs,
type="h",
xlim=c(-1,n+1),
ylim=c(0,1.1*max(probs)),
main=paste0("Distribuicao binomial (",n,", ",p,")"),
lwd=1,
col="black",
ylab="Probabilidade")
points(X,probs,pch=16,cex=1,col="black")
dev.off()Distribuição binomial: Questão 1
Qual é a probabilidade que nenhum paciente que realizará transplante seja Rh-?
R.: \(P(X=0) = 0.61\)
[1] 0.614125Distribuição binomial: Questão 2
Qual é a probabilidade que pelo menos um paciente que realizará transplante seja Rh-?
R.: \(1 – P(X=0) = 1 – 0.61 = 0.39\)
[1] 0.385875Distribuição binomial: Questão 3
Qual é a probabilidade que todos os pacientes que realizarão transplante sejam Rh-?
R.: \(P(X=3) = 0.003\)
[1] 0.003375Distribuição de Poisson
- Distribuição de probabilidade de variável discreta infinita \(X\)
- A distribuição de Poisson surge de um experimento de contagem de \(X\) ocorrências de um evento de interesse por determinado período de tempo, área ou volume.
- \(X \sim Poisson(\lambda)\): variável aleatória \(X\) representando o número de ocorrências do evento com distribuição de Poisson de taxa média de ocorrências \(\lambda\)
- \(X=0,1,2,...\)
- Média de \(X\): \(\lambda\)
- Desvio-padrão de \(X\): \(\sqrt{\lambda}\)
- Suposições
- Modela evento raro;
- As condições permanecem estáveis no decorrer do tempo ou espaço, i.e., a taxa média de ocorrências,\(\lambda\), é constante;
- Intervalos de tempo ou espaço disjuntos são independentes, i.e., a informação sobre o número de ocorrências num intervalo nada revela sobre o número de ocorrências em outro intervalo;
- A probabilidade de ocorrência de um evento em um certo intervalo é a mesma para todos os demais intervalos de tempo;
- A probabilidade de ocorrência dos eventos é proporcional ao tamanho do intervalo;
- Em uma porção infinitesimal do intervalo, a probabilidade de mais de uma ocorrência do evento é desprezível.
k <- 0:10
plot(k,dpois(x=k,lambda=2),
type='h',
main='Distribuicao de Poisson(lambda=2)',
ylab='Probabilidade',
xlab ='x',
lwd=3)
Exemplo de distribuição de Poisson: 1
Número de atendimentos completos em pronto-atendimento
O número de pacientes que têm atendimento completo num pronto-socorro de uma pequena cidade durante uma madrugada, \(X\), tem distribuição de Poisson com taxa média \(\lambda=3\).
Portanto,
- \(X \sim Poisson(3)\)
- Média de \(X\): \(3\)
- Desvio-padrão de \(X\): \(\sqrt{3} \approx 1.732\)
Se forem consideradas duas madrugadas, a taxa média de atendimentos completos é \(6 = 2 \times 3\).
Num mês espera-se que a taxa média de atendimentos completos durante a madrugada seja \(90 = 30 \times 3\).
Siqueira & Tibúrcio, 2011, p. 201-2
X ~ Poisson(lambda = 3 )
Media de X = 3
Desvio-padrao de X = 1.732051
X probs probsacum
1 0 0.049787 0.049787
2 1 0.149361 0.199148
3 2 0.224042 0.423190
4 3 0.224042 0.647232
5 4 0.168031 0.815263
6 5 0.100819 0.916082
7 6 0.050409 0.966491
8 7 0.021604 0.988095
9 8 0.008102 0.996197
10 9 0.002701 0.998898
11 10 0.000810 0.999708
12 11 0.000221 0.999929
P(X=0) + P(X=1) + ... + P(X=10) = 0.9999286
P(X=0) + P(X=1) = 0.1991483
P(X=2) + P(X=3) + ... = 0.8008517 
sink("exemplo02_Poisson.txt")
lambda <- 3
m <- lambda
dp <- sqrt(lambda)
X <- 0:(m+5*dp)
cat("X ~ Poisson(lambda =", lambda,")\n")
cat("Media de X =",m,"\n")
cat("Desvio-padrao de X =",dp,"\n\n")
probs <- dpois(x=X,lambda=lambda)
probsacum <- ppois(q=X,lambda=lambda)
print(tbl <- round(data.frame(X,probs,probsacum),6))
cat("\nP(X=0) + P(X=1) + ... + P(X=10) =",sum(probs),"\n")
# lower.tail: logical; if TRUE (default), probabilities are P[X <= x],
# otherwise, P[X > x]
px01 <- ppois(q=1,lambda=lambda,lower.tail=TRUE)
cat("P(X=0) + P(X=1) =", px01, "\n")
px23inf <- ppois(q=1,lambda=lambda,lower.tail=FALSE)
cat("P(X=2) + P(X=3) + ... =", px23inf, "\n")
sink()
png("exemplo02_Poisson.png")
plot(X,probs,type="h",
xlim=c(-1,m+5*dp),ylim=c(0,1.1*max(probs)),
main=paste("Distribuicao de Poisson(",lambda,")"),
lwd=1,col="black",ylab="Probabilidade")
points(X,probs,pch=16,cex=1,col="black")
dev.off()Distribuição de Poisson: Questão 1
Qual é a probabilidade que nenhum paciente tenha atendimento completo durante uma madrugada no pronto-atendimento dessa pequena cidade?
R.: \(P(X=0) = 0.05\)
[1] 0.04978707Distribuição de Poisson: Questão 2
Qual é a probabilidade que pelo menos um paciente tenha atendimento completo durante uma madrugada no pronto-atendimento dessa pequena cidade?
R.: \(1 – P(X=0) = 1 – 0.05 = 0.95\)
[1] 0.9502129Distribuição de Poisson: Questão 3
Qual é a probabilidade que mais de dez pacientes tenham atendimento completo durante uma madrugada no pronto-atendimento dessa pequena cidade?
R.: \(P(X>10) = 0.0003\)
[1] 0.000292337[1] 0.000292337Exemplo de distribuição de Poisson: 2
Exemplo: Número de consultas médicas anual de plano de saúde
Siqueira & Tibúrcio, 2011, p. 202
O número de consultas médicas anual de um associado de um plano de saúde é finito.
No entanto, uma aproximação pode ser feita supondo que o número de consultas anual de um associado pode ser infinito.
Num plano de saúde com 5694 filiados, ao fim de um ano foram realizadas 13098 consultas, conforme tabela a seguir:
| Número de consultas anuais de um associado | Frequência |
|---|---|
| 0 | 589 |
| 1 | 1274 |
| 2 | 1542 |
| 3 | 1144 |
| 4 | 663 |
| 5 | 304 |
| 6 | 126 |
| 7 | 39 |
| 8 | 10 |
| 9 | 3 |
\[\begin{align} \lambda&=\dfrac{0\times589+1\times1274+\cdots+9\times3}{589+1274+\cdots+3}\\ \lambda&=2.3 \;\text{consultas por associado por ano} \end{align}\]
O valor 2.3 é o número médio de consultas anual de um associado típico deste plano de saúde.
X ~ Poisson(lambda = 2.3 )
Media de X = 2.3
Desvio-padrao de X = 1.516575
X probs probsacum
1 0 0.100259 0.100259
2 1 0.230595 0.330854
3 2 0.265185 0.596039
4 3 0.203308 0.799347
5 4 0.116902 0.916249
6 5 0.053775 0.970024
7 6 0.020614 0.990638
8 7 0.006773 0.997411
9 8 0.001947 0.999358
10 9 0.000498 0.999856
P(X=0) + P(X=1) + ... + P(X=10) = 0.9998561
P(X=0) + P(X=1) = 0.3308542
P(X=2) + P(X=3) + ... = 0.6691458 
sink("exemplo03_Poisson.txt")
lambda <- 2.3
m <- lambda
dp <- sqrt(lambda)
X <- 0:(m+5*dp)
cat("X ~ Poisson(lambda =", lambda,")\n")
cat("Media de X =",m,"\n")
cat("Desvio-padrao de X =",dp,"\n\n")
probs <- dpois(x=X,lambda=lambda)
probsacum <- ppois(q=X,lambda=lambda)
print(tbl <- round(data.frame(X,probs,probsacum),6))
cat("\nP(X=0) + P(X=1) + ... + P(X=10) =",sum(probs),"\n")
# lower.tail: logical; if TRUE (default), probabilities are P[X <= x],
# otherwise, P[X > x]
px01 <- ppois(q=1,lambda=lambda,lower.tail=TRUE)
cat("P(X=0) + P(X=1) =", px01, "\n")
px23inf <- ppois(q=1,lambda=lambda,lower.tail=FALSE)
cat("P(X=2) + P(X=3) + ... =", px23inf, "\n")
sink()
png("exemplo03_Poisson.png")
plot(X,probs,type="h",
xlim=c(-1,(m+5*dp)),ylim=c(0,1.1*max(probs)),
main=paste("Distribuicao de Poisson(",lambda,")"),
lwd=1,col="black",ylab="Probabilidade")
points(X,probs,pch=16,cex=1,col="black")
dev.off()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\Re\)
- 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
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
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=0.025, mean=120, sd=10),2)
q2 <- round(qnorm(p=0.975, 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("",0.95,""),
xlim=c(80,160),
main="normal(120,10)",
col=c("black", "black"),
density=c(0,20,0))
Relações entre distribuições de probabilidade
Relações da distribuição normal com outras distribuições discretas e contínuas.
Aproximação da binomial pela Poisson
\(X \sim binomial(n,p)\)
Se \(n > 20\) e \(p<0.05\), então a binomial é aproximadamente igual à Poisson\((np)\).
X ~ binomial(10, 0.1)
Media de X = 1
Desvio-padrao de X = 0.95
Binomial(10, 0.1) diferente de Poisson(1)
X probbinom probPois
1 0 0.348678 0.367879
2 1 0.387420 0.367879
3 2 0.193710 0.183940
4 3 0.057396 0.061313
5 4 0.011160 0.015328
6 5 0.001488 0.003066
7 6 0.000138 0.000511
8 7 0.000009 0.000073
9 8 0.000000 0.000009
10 9 0.000000 0.000001
11 10 0.000000 0.000000
X ~ binomial(30, 0.03)
Media de X = 0.9
Desvio-padrao de X = 0.93
Binomial(30, 0.03) semelhante a Poisson(0.9)
X probbinom probPois
1 0 0.401007 0.406570
2 1 0.372068 0.365913
3 2 0.166855 0.164661
4 3 0.048164 0.049398
5 4 0.010055 0.011115
6 5 0.001617 0.002001
7 6 0.000208 0.000300
8 7 0.000022 0.000039
9 8 0.000002 0.000004
10 9 0.000000 0.000000
11 10 0.000000 0.000000
12 11 0.000000 0.000000
13 12 0.000000 0.000000
14 13 0.000000 0.000000
15 14 0.000000 0.000000
16 15 0.000000 0.000000
17 16 0.000000 0.000000
18 17 0.000000 0.000000
19 18 0.000000 0.000000
20 19 0.000000 0.000000
21 20 0.000000 0.000000
22 21 0.000000 0.000000
23 22 0.000000 0.000000
24 23 0.000000 0.000000
25 24 0.000000 0.000000
26 25 0.000000 0.000000
27 26 0.000000 0.000000
28 27 0.000000 0.000000
29 28 0.000000 0.000000
30 29 0.000000 0.000000
31 30 0.000000 0.000000
sink("exemplo09_binomial_aproxPois_n30p3.txt")
n <- 30
p <- 0.03
cat("X ~ binomial(",n,", ",p,")\n", sep="")
m <- n*p
dp <- sqrt(n*p*(1-p))
cat("Media de X = ",m,"\n", sep="")
cat("Desvio-padrao de X = ",round(dp,2),"\n", sep="")
X <- 0:n
if (n >= 20 && p <= 0.05)
{cat("Binomial(",n,", ",p,") semelhante a Poisson(",m,")\n\n", sep="")} else
{cat("Binomial(",n,", ",p,") diferente de Poisson(",m,")\n\n", sep="")}
probbinom <- dbinom(x=X,size=n,prob=p)
probPois <- dpois(x=X,lambda=m)
print(round((tbl <- data.frame(X,probbinom,probPois)),6))
sink()
png("exemplo09_binomial_aproxPois_n30p3.png")
plot(X,probbinom,type="h",
main=paste0("Distribuicao binomial(",n,", ",p,"): circulo\n\n"),
xlab="X", ylab="probabilidade",
xlim=c(-1,m+5*dp), ylim=c(0,1.1*max(c(max(probbinom),max(probPois)))),
lwd=1,col="black")
points(X,probbinom,pch=16,cex=1,col="black")
par(new=TRUE)
plot(X,probPois,type="h",
main=paste0("Distribuicao de Poisson(",m,"): quadrado"),
xlab="X", ylab="probabilidade",
xlim=c(-1,m+5*dp), ylim=c(0,1.1*max(c(max(probbinom),max(probPois)))),
lwd=1,col="black")
points(X,probPois,pch=12,cex=1,col="black")
dev.off()Aproximação da binomial pela Poisson: Questão 1
Anestésico
Uma em cada mil pessoas que utilizam determinado anestésico sofre uma reação negativa (choque).
Num total de 500 cirurgias em que se empregou esse anestésico, a probabilidade de que uma pessoa sofra a reação negativa é:
R.: \(0.3\)
Como \(n = 500 > 20\) e \(p = 1/1000 \le 0.05\), então a \(binomial(500,1/1000)\) é aproximadamente igual à Poisson\((500\times(1/1000))\).
[1] 0.303493[1] 0.3032653Arango, 2012, p. 187
ARANGO, HG (2012) Bioestatística. 3a ed. RJ: Guanabara Koogan.
Aproximação da binomial pela Poisson: Questão 2
Anestésico
Uma em cada mil pessoas que utilizam determinado anestésico sofre uma reação negativa (choque).
Num total de 500 cirurgias em que se empregou esse anestésico a probabilidade de que nenhuma pessoa sofra a reação negativa é:
R.: \(0.61\)
Como \(n = 500 > 20\) e \(p = 1/1000 \le 0.05\), então a \(binomial(500,1/1000)\) é aproximadamente igual à Poisson\((500\times(1/1000))\).
[1] 0.6063789[1] 0.6065307Aproximação da binomial pela Poisson: Questão 3
Anestésico
Uma em cada mil pessoas que utilizam determinado anestésico sofre uma reação negativa (choque).
Num total de 500 cirurgias em que se empregou esse anestésico a probabilidade de que mais de uma pessoa sofra a reação negativa é:
R.: \(0.09\)
Como \(n = 500 > 20\) e \(p = 1/1000 \le 0.05\), então a \(binomial(500,1/1000)\) é aproximadamente igual à Poisson\((500\times(1/1000))\).
[1] 0.09012809[1] 0.09020401Aproximação da binomial pela normal
- \(X \sim binomial(n, p)\)
- Se
- \(np > 5\) e
- \(n(1-p) > 5\) e
- \({\frac{1}{\sqrt{n}} \left|\sqrt{\frac{1-p}{p}}-\sqrt{\frac{p}{1-p}}\right|} < 0.3\),
- \(binomial(n, p)\) é aproximadamente igual à \(normal\left(np, \sqrt{np(1-p)}\right)\)
X ~ binomial(6, 0.7)
Media de X = 4.2
Desvio-padrao de X = 1.12
Binomial(6, 0.7) diferente da Normal(4.2, 1.12)
X probbinom probnorm
1 0 0.000729 0.000324
2 1 0.010206 0.006109
3 2 0.059535 0.052072
4 3 0.185220 0.200704
5 4 0.324135 0.349809
6 5 0.302526 0.275694
7 6 0.117649 0.098253
X ~ binomial(30, 0.7)
Media de X = 21
Desvio-padrao de X = 2.51
Binomial(30, 0.7) semelhante a Normal(21, 2.51)
X probbinom probnorm
1 0 0.000000 0.000000
2 1 0.000000 0.000000
3 2 0.000000 0.000000
4 3 0.000000 0.000000
5 4 0.000000 0.000000
6 5 0.000000 0.000000
7 6 0.000000 0.000000
8 7 0.000000 0.000000
9 8 0.000001 0.000000
10 9 0.000006 0.000002
11 10 0.000030 0.000011
12 11 0.000126 0.000057
13 12 0.000464 0.000257
14 13 0.001498 0.000989
15 14 0.004246 0.003253
16 15 0.010567 0.009128
17 16 0.023115 0.021855
18 17 0.044418 0.044643
19 18 0.074852 0.077809
20 19 0.110308 0.115709
21 20 0.141562 0.146816
22 21 0.157291 0.158942
23 22 0.150141 0.146816
24 23 0.121854 0.115709
25 24 0.082928 0.077809
26 25 0.046440 0.044643
27 26 0.020838 0.021855
28 27 0.007203 0.009128
29 28 0.001801 0.003253
30 29 0.000290 0.000989
31 30 0.000023 0.000257
sink("exemplo07_binomial_aproxnormal_n30.txt")
n <- 30
p <- 0.7
cat("X ~ binomial(",n,", ",p,")\n", sep="")
m <- n*p
dp <- sqrt(n*p*(1-p))
cat("Media de X = ",m,"\n", sep="")
cat("Desvio-padrao de X = ",round(dp,2),"\n", sep="")
X <- 0:n
if (n*p > 5 && n*(1-p) > 5 && abs((1/sqrt(n))*(sqrt((1-p)/p)-sqrt(p/(1-p)))) < 0.3)
{cat("\nBinomial(",n,", ",p,") semelhante a Normal(",m,", ",round(dp,2),")\n\n", sep="")} else
{cat("\nBinomial(",n,", ",p,") diferente da Normal(",m,", ",round(dp,2),")\n\n", sep="")}
probbinom <- dbinom(x=X,size=n,prob=p)
probnorm <- dnorm(x=X,mean=m,sd=dp)
print(round((tbl <- data.frame(X,probbinom,probnorm)),6))
sink()
png("exemplo07_binomial_aproxnormal_n30.png")
plot(X,probbinom,type="h",
main=paste0("Distribuicao binomial(",n,", ",p,")\n\n"),
xlab="X", ylab="densidade & probabilidade",
xlim=c(-1,n+1), ylim=c(0,1.1*max(c(max(probbinom),max(probnorm)))),
lwd=1,col="black")
points(X,probbinom,pch=16,cex=1,col="black")
par(new=TRUE)
rn <- rnorm(n=1e6,mean=m,sd=dp)
plot(density(rn), xlab="X", ylab="densidade & probabilidade",
main=paste0("Distribuicao normal(",m,", ",round(dp,2),")"),
xlim=c(-1,n+1), ylim=c(0,1.1*max(c(max(probbinom),max(probnorm)))))
dev.off()Aproximação da Poisson pela normal
\(X \sim Poisson(\lambda)\)
Se \(\lambda > 10\), então a binomial é aproximadamente igual à \(normal\left(\lambda, \sqrt{\lambda}\right)\).
X ~ Poisson(lambda = 5)
Media de X = 5
Desvio-padrao de X = 2.24
Poisson(5) diferente da normal(5, 2.24)
X probPois probnorm
1 0 0.006738 0.014645
2 1 0.033690 0.036021
3 2 0.084224 0.072537
4 3 0.140374 0.119593
5 4 0.175467 0.161434
6 5 0.175467 0.178412
7 6 0.146223 0.161434
8 7 0.104445 0.119593
9 8 0.065278 0.072537
10 9 0.036266 0.036021
11 10 0.018133 0.014645
12 11 0.008242 0.004875
X ~ Poisson(lambda = 50)
Media de X = 50
Desvio-padrao de X = 7.07
Poisson(50) semelhante a normal(50, 7.07)
X probPois probnorm
1 0 0.000000 0.000000
2 1 0.000000 0.000000
3 2 0.000000 0.000000
4 3 0.000000 0.000000
5 4 0.000000 0.000000
6 5 0.000000 0.000000
7 6 0.000000 0.000000
8 7 0.000000 0.000000
9 8 0.000000 0.000000
10 9 0.000000 0.000000
11 10 0.000000 0.000000
12 11 0.000000 0.000000
13 12 0.000000 0.000000
14 13 0.000000 0.000000
15 14 0.000000 0.000000
16 15 0.000000 0.000000
17 16 0.000000 0.000001
18 17 0.000000 0.000001
19 18 0.000000 0.000002
20 19 0.000000 0.000004
21 20 0.000001 0.000007
22 21 0.000002 0.000013
23 22 0.000004 0.000022
24 23 0.000009 0.000038
25 24 0.000019 0.000065
26 25 0.000037 0.000109
27 26 0.000071 0.000178
28 27 0.000132 0.000284
29 28 0.000236 0.000446
30 29 0.000406 0.000686
31 30 0.000677 0.001033
32 31 0.001092 0.001526
33 32 0.001707 0.002210
34 33 0.002586 0.003136
35 34 0.003803 0.004361
36 35 0.005432 0.005947
37 36 0.007545 0.007947
38 37 0.010196 0.010410
39 38 0.013416 0.013367
40 39 0.017200 0.016824
41 40 0.021500 0.020755
42 41 0.026219 0.025098
43 42 0.031213 0.029749
44 43 0.036294 0.034564
45 44 0.041244 0.039362
46 45 0.045826 0.043939
47 46 0.049811 0.048077
48 47 0.052991 0.051563
49 48 0.055199 0.054207
50 49 0.056325 0.055858
51 50 0.056325 0.056419
52 51 0.055221 0.055858
53 52 0.053097 0.054207
54 53 0.050091 0.051563
55 54 0.046381 0.048077
56 55 0.042164 0.043939
57 56 0.037647 0.039362
58 57 0.033023 0.034564
59 58 0.028468 0.029749
60 59 0.024126 0.025098
61 60 0.020105 0.020755
62 61 0.016479 0.016824
63 62 0.013290 0.013367
64 63 0.010547 0.010410
65 64 0.008240 0.007947
66 65 0.006339 0.005947
67 66 0.004802 0.004361
68 67 0.003584 0.003136
69 68 0.002635 0.002210
70 69 0.001909 0.001526
71 70 0.001364 0.001033
72 71 0.000960 0.000686
sink("exemplo08_Poisson_aproxnormal_n50.txt")
lambda <- 50
m <- lambda
dp <- sqrt(lambda)
X <- 0:(m+3*dp)
cat("X ~ Poisson(lambda = ", lambda,")\n", sep="")
cat("Media de X = ",m,"\n", sep="")
cat("Desvio-padrao de X = ",round(dp,2),"\n\n", sep="")
if (lambda > 10)
{cat("Poisson(",lambda,") semelhante a normal(",m,", ",round(dp,2),")\n\n", sep="")} else
{cat("Poisson(",lambda,") diferente da normal(",m,", ",round(dp,2),")\n\n", sep="")}
probPois <- dpois(x=X,lambda=lambda)
probnorm <- dnorm(x=X,mean=m,sd=dp)
print(round((tbl <- data.frame(X,probPois,probnorm)),6))
sink()
png("exemplo08_Poisson_aproxnormal_n50.png")
plot(X,probPois,type="h",
xlim=c(-1,m+3*dp),ylim=c(0,1.1*max(probPois)),
xlab="X", ylab="densidade & probabilidade",
main=paste0("Distribuicao de Poisson(",lambda,")\n\n"),
lwd=1,col="black")
points(X,probPois,pch=16,cex=1,col="black")
par(new=TRUE)
rn <- rnorm(n=1e6,mean=m,sd=dp)
plot(density(rn),
xlab="X", ylab="densidade & probabilidade",
main=paste0("Distribuicao normal(",m,", ",round(dp,2),")"),
xlim=c(-1,m+3*dp), ylim=c(0,1.1*max(probPois)))
dev.off()Distribuição normal e t em R
HH::normal.and.t.dist(std.dev=2, Use.alpha.left=TRUE, deg.free=6,
gxbar.max=.20, polygon.density=10)
HH::normal.and.t.dist(std.dev=2, Use.alpha.left=FALSE, deg.free=6,
gxbar.max=.20, polygon.density=10,
mu.H1=2, Use.mu.H1=TRUE,
obs.mean=2.5, Use.obs.mean=TRUE, xmin=-7)
HH::norm.setup(xlim=c(75,125), mean=100, se=5)
HH::norm.curve(100, 5, 100+5*(1.645))
HH::norm.observed(112, (112-100)/5)
HH::norm.outline("dnorm", 112, par()$usr[2], 100, 5)
HH::norm.setup(xlim=c(-3, 6))
HH::norm.curve(critical.values=1.645, mean=1.645+1.281552, col='green',
shade="left", axis.name="z1")
HH::norm.curve(critical.values=1.645, col='red')
HH::norm.setup(xlim=c(-6, 12), se=2)
HH::norm.curve(critical.values=2*1.645, se=2, mean=2*(1.645+1.281552),
col='green', shade="left", axis.name="z1")
HH::norm.curve(critical.values=2*1.645, se=2, mean=0,
col='red', shade="right")
par(mfrow=c(2,1))
HH::norm.setup()
HH::norm.curve()
mtext("norm.setup(); norm.curve()", side=1, line=5)
HH::norm.setup(n=1)
HH::norm.curve(n=1)
mtext("norm.setup(n=1); norm.curve(n=1)", side=1, line=5)
par(mfrow=c(1,1))
par(mfrow=c(2,2))
## naively scaled,
## areas under the curve are numerically the same but visually different
HH::norm.setup(n=1)
HH::norm.curve(n=1)
HH::norm.observed(1.2, 1.2/(1/sqrt(1)))
HH::norm.setup(n=2)
HH::norm.curve(n=2)
HH::norm.observed(1.2, 1.2/(1/sqrt(2)))
HH::norm.setup(n=4)
HH::norm.curve(n=4)
HH::norm.observed(1.2, 1.2/(1/sqrt(4)))
HH::norm.setup(n=10)
HH::norm.curve(n=10)
HH::norm.observed(1.2, 1.2/(1/sqrt(10)))
mtext("areas under the curve are numerically the same but visually different",
side=3, outer=TRUE)
## scaled so all areas under the curve are numerically and visually the same
HH::norm.setup(n=1, ylim=c(0,1.3))
HH::norm.curve(n=1)
HH::norm.observed(1.2, 1.2/(1/sqrt(1)))
HH::norm.setup(n=2, ylim=c(0,1.3))
HH::norm.curve(n=2)
HH::norm.observed(1.2, 1.2/(1/sqrt(2)))
HH::norm.setup(n=4, ylim=c(0,1.3))
HH::norm.curve(n=4)
HH::norm.observed(1.2, 1.2/(1/sqrt(4)))
HH::norm.setup(n=10, ylim=c(0,1.3))
HH::norm.curve(n=10)
HH::norm.observed(1.2, 1.2/(1/sqrt(10)))
mtext("all areas under the curve are numerically and visually the same",
side=3, outer=TRUE)
par(mfrow=c(1,1))
## t distribution
mu.H0 <- 16
se.val <- .4
df.val <- 10
crit.val <- mu.H0 - qt(.95, df.val) * se.val
mu.alt <- 15
obs.mean <- 14.8
alt.t <- (mu.alt - crit.val) / se.val
HH::norm.setup(xlim=c(12, 19), se=se.val, df.t=df.val)
HH::norm.curve(critical.values=crit.val, se=se.val,
df.t=df.val, mean=mu.alt,
col='green', shade="left", axis.name="t1")
HH::norm.curve(critical.values=crit.val, se=se.val,
df.t=df.val, mean=mu.H0,
col='gray', shade="right")
HH::norm.observed(obs.mean, (obs.mean-mu.H0)/se.val)
## normal
HH::norm.setup(xlim=c(12, 19), se=se.val)
HH::norm.curve(critical.values=crit.val, se=se.val, mean=mu.alt,
col='green', shade="left", axis.name="z1")
HH::norm.curve(critical.values=crit.val, se=se.val, mean=mu.H0,
col='gray', shade="right")
HH::norm.observed(obs.mean, (obs.mean-mu.H0)/se.val)
## normal and t
HH::norm.setup(xlim=c(12, 19), se=se.val, main="t(6) and normal")
HH::norm.curve(critical.values=15.5, se=se.val, mean=16.3,
col='gray', shade="right")
HH::norm.curve(critical.values=15.5, se.val, df.t=6, mean=14.7,
col='green', shade="left", axis.name="t1", second.axis.label.line=4)
HH::norm.curve(critical.values=15.5, se=se.val, mean=16.3,
col='gray', shade="none")
HH::norm.setup(xlim=c(12, 19), se=se.val, main="t(6) and normal")
HH::norm.curve(critical.values=15.5, se=se.val, mean=15.5,
col='gray', shade="right")
HH::norm.curve(critical.values=15.5, se=se.val, df.t=6, mean=15.5,
col='green', shade="left", axis.name="t1", second.axis.label.line=4)
HH::norm.curve(critical.values=15.5, se=se.val, mean=15.5,
col='gray', shade="none")