Prática R: Fundamentos de Inferência Estatística

MSP4112: Estatística Aplicada à Medicina

Bastão de Asclépio & Distribuição Normal

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))

suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(gplots, warn.conflicts=FALSE))
suppressMessages(library(kableExtra, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))

Material

• HTML de R Markdown em RPubs

Introdução

Esta prática é relativa aos capítulos 3 (Estatística Descritiva) e 4 (Probabilidade, Amostragem e Distribuições) do livro adotado:

• Dancey, CP & Reidy, J (2019) Estatística sem Matemática para Psicologia. 7a ed. Porto Alegre: Penso.

Prepare-se

SEMPRE crie um projeto em uma pasta separada em seu computador.

Baixe os arquivos Biometria_FMUSP.rds e NewDrug.rds para seu computador:

• https://drive.google.com/drive/folders/1dV4df6I1A4UAKgNXvztS1A8IznwPt3ut?usp=sharing

Atividade 1

ler Biometria_FMUSP.rds

Leia o arquivo Biometria_FMUSP.rds.

Dados <- readRDS("Biometria_FMUSP.rds")

estatística descritiva por sexo

Considerando a variável Sexo, há dois grupos neste dataframe, masculino e feminino.

Exiba a estatística descritiva numérica (média, mediana, desvio-padrão, etc.) das variáveis intervalares (neste caso, MCT, Estatura e IMC) por sexo.

print(psych::describeBy(Estatura ~ Sexo, data=data.frame(Dados),
                        mat=1, digits=2))

          item group1 vars   n   mean   sd median trimmed  mad min max range
Estatura1    1      F    1 229 164.34 6.38    163  163.94 5.93 153 186    33
Estatura2    2      M    1 313 175.77 7.26    175  175.88 7.41 155 195    40
           skew kurtosis   se
Estatura1  0.61     0.02 0.42
Estatura2 -0.13     0.21 0.41

print(psych::describeBy(MCT ~ Sexo, data=data.frame(Dados), 
                        mat=1, digits=2))

     item group1 vars   n  mean    sd median trimmed   mad min max range skew
MCT1    1      F    1 230 57.63  9.05     56   56.75  7.41  41 120    79 2.43
MCT2    2      M    1 312 71.59 12.09     72   70.88 11.12  43 125    82 0.90
     kurtosis   se
MCT1    12.04 0.60
MCT2     2.37 0.68

print(psych::describeBy(IMC ~ Sexo, data=data.frame(Dados), 
                        mat=1, digits=2))

     item group1 vars   n  mean   sd median trimmed  mad   min   max range skew
IMC1    1      F    1 229 21.24 2.76  20.70   20.95 1.88 16.28 36.73 20.45 2.09
IMC2    2      M    1 312 23.12 3.29  22.86   22.89 3.06 14.53 40.35 25.82 1.10
     kurtosis   se
IMC1     8.13 0.18
IMC2     3.16 0.19

gráficos

Exiba a distribuição dos valores das variáveis intervalares MCT, Estatura e IMC por sexo.

Para comparação, o melhor é que ambos os grupos apareçam no mesmo gráfico.

boxplot(Estatura ~ Sexo, data=Dados)

boxplot(MCT ~ Sexo, data=Dados)

boxplot(IMC ~ Sexo, data=Dados)

car::densityPlot(Estatura ~ Sexo, data=Dados)

car::densityPlot(MCT ~ Sexo, data=Dados)

car::densityPlot(IMC ~ Sexo, data=Dados)

comparação de médias populacionais

Verifique se as distribuições de MCT, Estatura e IMC em cada grupo têm distribuição normal.

Calcule os intervalos de confiança de 95% da média populacional de MCT, Estatura e IMC em cada sexo e verifique se são diferentes.

Consulte as funções disponíveis no pacote DescTools e as funções MVN::mvn e gplots::plotmeans.

# https://en.wikipedia.org/wiki/Shapiro%E2%80%93Wilk_test
tapply(Dados$IMC, Dados$Sexo, shapiro.test)

$F

    Shapiro-Wilk normality test

data:  X[[i]]
W = 0.85432, p-value = 6.495e-14


$M

    Shapiro-Wilk normality test

data:  X[[i]]
W = 0.94337, p-value = 1.44e-09

tapply(Dados$MCT, Dados$Sexo, shapiro.test)

$F

    Shapiro-Wilk normality test

data:  X[[i]]
W = 0.83245, p-value = 4.94e-15


$M

    Shapiro-Wilk normality test

data:  X[[i]]
W = 0.95412, p-value = 2.614e-08

tapply(Dados$Estatura, Dados$Sexo, shapiro.test)

$F

    Shapiro-Wilk normality test

data:  X[[i]]
W = 0.96265, p-value = 1.053e-05


$M

    Shapiro-Wilk normality test

data:  X[[i]]
W = 0.99062, p-value = 0.04289

alfa <- 0.05
result <- MVN::mvn(data=subset(Dados, 
                               select=c(Sexo, 
                                        IMC, MCT, Estatura)), 
                   subset="Sexo",
                   alpha=alfa,
                   univariate_test="SW")

Warning in MVN::mvn(data = subset(Dados, select = c(Sexo, IMC, MCT, Estatura)),
: Missing values detected in 8 rows. These rows will be removed.

print(result$univariate_normality)

  Group         Test Variable Statistic p.value    Normality
1     F Shapiro-Wilk      IMC     0.854  <0.001 ✗ Not normal
2     F Shapiro-Wilk      MCT     0.905  <0.001 ✗ Not normal
3     F Shapiro-Wilk Estatura     0.963  <0.001 ✗ Not normal
4     M Shapiro-Wilk      IMC     0.943  <0.001 ✗ Not normal
5     M Shapiro-Wilk      MCT     0.954  <0.001 ✗ Not normal
6     M Shapiro-Wilk Estatura     0.991   0.045 ✗ Not normal

# Qual é o intervalo com 95% de confiança com correção de Bonferroni 
# da média populacional da estatura masculina?

Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")
gplots::plotmeans(IMC ~ Sexo, 
                  data=Dados,
                  p=1-alfa/length(unique(Dados$Sexo)), 
                  barcol="black", 
                  connect=FALSE)

DescTools::MeanCI(Dados.M$IMC, 
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE)

    mean   lwr.ci   upr.ci 
23.11829 22.69838 23.53819 

DescTools::MeanCI(Dados.F$IMC, 
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE)

    mean   lwr.ci   upr.ci 
21.23509 20.82358 21.64660 

DescTools::MeanDiffCI(Dados.M$IMC, Dados.F$IMC, 
                      paired=FALSE,                                                   conf.level=1-alfa,
                      na.rm=TRUE)

meandiff   lwr.ci   upr.ci 
1.883195 1.370850 2.395541 

set.seed(123)
DescTools::BootCI(Dados.M$IMC,
                  FUN=mean,
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE,
                  R=1e4)

    mean   lwr.ci   upr.ci 
23.11829 22.69981 23.53620 

set.seed(123)
DescTools::BootCI(Dados.F$IMC,
                  FUN=mean,
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE,
                  R=1e4)

    mean   lwr.ci   upr.ci 
21.23509 20.82472 21.64250 

set.seed(123)
DescTools::MeanDiffCI(Dados.M$IMC,
                      Dados.F$IMC,
                      paired=FALSE,
                      na.rm=TRUE,                                                     conf.level=1-alfa,
                      method="perc",
                      R=1e4)

meandiff   lwr.ci   upr.ci 
1.883195 1.362122 2.390445 

gplots::plotmeans(MCT ~ Sexo, 
                  data=Dados,
                  p=1-alfa/length(unique(Dados$Sexo)), 
                  barcol="black", 
                  connect=FALSE)

DescTools::MeanCI(Dados.M$MCT, 
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE)

    mean   lwr.ci   upr.ci 
71.58974 70.04848 73.13101 

DescTools::MeanCI(Dados.F$MCT, 
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE)

    mean   lwr.ci   upr.ci 
57.63043 56.28362 58.97725 

DescTools::MeanDiffCI(Dados.M$MCT, 
                      Dados.F$MCT, 
                      paired=FALSE,                                                   conf.level=1-alfa,
                      na.rm=TRUE)

meandiff   lwr.ci   upr.ci 
13.95931 12.17552 15.74310 

gplots::plotmeans(Estatura ~ Sexo, 
                  data=Dados,
                  p=1-alfa/length(unique(Dados$Sexo)), 
                  barcol="black", 
                  connect=FALSE)

DescTools::MeanCI(Dados.M$Estatura, 
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE)

    mean   lwr.ci   upr.ci 
175.7732 174.8487 176.6976 

DescTools::MeanCI(Dados.F$Estatura, 
                  conf.level=1-alfa/length(unique(Dados$Sexo)),
                  na.rm=TRUE)

    mean   lwr.ci   upr.ci 
164.3406 163.3896 165.2916 

DescTools::MeanDiffCI(Dados.M$Estatura, 
                      Dados.F$Estatura, 
                      paired=FALSE,                                                   conf.level=1-alfa,
                      na.rm=TRUE)

meandiff   lwr.ci   upr.ci 
11.43255 10.27678 12.58833 

Atividade 2

Investigue se o nível médio de sedentarismo é similar nos dois sexos.

Consulte as funções disponíveis no pacote DescTools.

Uma alternativa é lembrar que proporção é média; converta a variável nominal dicotômia de sedentarismo para intervalar binária com valores 0 e 1).

alfa <- 0.05
Dados$Sedentarismo.bin[Dados$Sedentarismo=="N"] <- 0

Warning: Unknown or uninitialised column: `Sedentarismo.bin`.

Dados$Sedentarismo.bin[Dados$Sedentarismo=="S"] <- 1
xtabs(~Dados$Sedentarismo.bin)

Dados$Sedentarismo.bin
  0   1 
316 233 

gplots::plotmeans(Sedentarismo.bin ~ Sexo, 
                  data=Dados,
                  p=1-alfa/length(unique(Dados$Sexo)), 
                  barcol="black", 
                  connect=FALSE)

Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")
DescTools::BinomCI(x=c(sum(Dados.M$Sedentarismo.bin, na.rm=TRUE),
                       sum(Dados.F$Sedentarismo.bin, na.rm=TRUE)),
                   n=c(sum(!is.na(Dados.M$Sedentarismo.bin)),
                       sum(!is.na(Dados.F$Sedentarismo.bin))),
                   conf.level=1-alfa/length(unique(Dados$Sexo)))

              est    lwr.ci    upr.ci
x.1:n.1 0.4227129 0.3622056 0.4856318
x.2:n.2 0.4267241 0.3562528 0.5003018

DescTools::BinomDiffCI(x1=sum(Dados.M$Sedentarismo.bin, na.rm=TRUE),
                       x2=sum(Dados.F$Sedentarismo.bin, na.rm=TRUE),
                       n1=sum(!is.na(Dados.M$Sedentarismo.bin)),
                       n2=sum(!is.na(Dados.F$Sedentarismo.bin)),
                       conf.level=1-alfa)

              est      lwr.ci     upr.ci
[1,] -0.004011204 -0.08756152 0.07925565

Atividade 3

Nesta atividade, a análise é intraparticipantes, pois cada paciente teve três medidas de respiração e pulso.

Usaremos o arquivo wide armazenado em NewDrug.rds.

Suponha que as mensurações tenham sido feitas em três semanas diferentes.

Investigue se as mensurações de pulso feitas na terceira semana (por exemplo, ao final do tratamento) diferem das mensurações iniciais (na primeira semana) em cada um dos tratamentos (os que tomaram a nova droga e os que receberam placebo).

Consulte as funções disponíveis no pacote DescTools.

alfa <- 0.05
Drug <- readRDS("NewDrug.rds")
Drug <- Drug$wide
Drug$Difpulse31 <- Drug$pulse3 - Drug$pulse1
Drug.New <- subset(Drug, drug=="New Drug")
Drug.Pcb <- subset(Drug, drug=="Placebo")

gplots::plotmeans(Difpulse31 ~ drug, 
                  data=Drug,
                  p=1-alfa/length(unique(Drug$drug)), 
                  barcol="black", 
                  connect=FALSE)
abline(h=0, lty=2)

# versao analitica
DescTools::MeanCI(Drug.New$Difpulse31, 
                  conf.level=1-alfa,
                  na.rm=TRUE)

       mean      lwr.ci      upr.ci 
 0.08333333 -0.10922214  0.27588881 

DescTools::MeanDiffCI(Drug.New$pulse3, 
                      Drug.New$pulse1, 
                      conf.level=1-alfa,
                      paired=TRUE, 
                      na.rm=TRUE)

   meandiff      lwr.ci      upr.ci 
 0.08333333 -0.10922214  0.27588881 

DescTools::MeanCI(Drug.Pcb$Difpulse31, 
                  conf.level=1-alfa,
                  na.rm=TRUE)

       mean      lwr.ci      upr.ci 
 0.11666667 -0.03780608  0.27113941 

DescTools::MeanDiffCI(Drug.Pcb$pulse3, 
                      Drug.Pcb$pulse1, 
                      conf.level=1-alfa,
                      paired=TRUE, 
                      na.rm=TRUE)

   meandiff      lwr.ci      upr.ci 
 0.11666667 -0.03780608  0.27113941 

# versao bootstrapping
set.seed(123)
DescTools::BootCI(Drug.New$Difpulse31, 
                  FUN=mean, 
                  na.rm=TRUE, 
                  conf.level=1-alfa,
                  bci.method="perc",
                  R=1e5)

       mean      lwr.ci      upr.ci 
 0.08333333 -0.03333333  0.23333333 

set.seed(123)
DescTools::MeanDiffCI(Drug.New$pulse3, 
                      Drug.New$pulse1, 
                      paired=TRUE, 
                      na.rm=TRUE,
                      conf.level=1-alfa,
                      method="perc", 
                      R=1e5)

   meandiff      lwr.ci      upr.ci 
 0.08333333 -0.03333333  0.23333333 

set.seed(123)
DescTools::BootCI(Drug.Pcb$Difpulse31, 
                  FUN=mean, 
                  na.rm=TRUE, 
                  conf.level=1-alfa,
                  bci.method="perc", 
                  R=1e5)

        mean       lwr.ci       upr.ci 
1.166667e-01 2.220446e-16 2.166667e-01 

set.seed(123)
DescTools::MeanDiffCI(Drug.Pcb$pulse3, 
                      Drug.Pcb$pulse1, 
                      paired=TRUE, 
                      na.rm=TRUE,
                      conf.level=1-alfa,
                      method="perc", 
                      R=1e5)

    meandiff       lwr.ci       upr.ci 
1.166667e-01 2.220446e-16 2.166667e-01