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

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

options(warn=-1)
suppressMessages(library(bootES, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(effsize, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lsr, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(sjPlot, warn.conflicts=FALSE))
source("summarySEwithin2.R")

Material

HTML de Rmarkdown em RPubs
- Prática R: Teste t

Introdução

Esta prática é relativa ao capítulo 7 do livro-texto:

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

Teste t

A função stats::t.test realiza os seguintes testes de médias populacionais:

Teste t para uma condição,
Duas condições independentes
- t de Student (homocedástico)
- t de Welch-Satterthwaite (heterocedástico) (default)
Duas condições dependentes
- t relacionado ou para medidas repetidas

O planejamento do teste t pode ser realizado por meio da função pwr::pwr.t.test.

Uma condição

Teste t bilateral: Estatura

Ler Biometria_FMUSP.rds.

Dados <- data.frame(readRDS("Biometria_FMUSP.rds"))

A distribuição da estatura do homem brasileiro de 19 anos de 2016 é normal com média igual a \(\mu = 177~cm\). Assumir que o desvio-padrão populacional é desconhecido.

Verificar se os estudantes de medicina de 2015-2018 têm estatura média igual a da população masculina brasileira de 19 anos de 2016.

Testet_1condicao.R

Teste t unilateral: Estatura

Verificar se os estudantes de medicina de 2015-2018 têm estatura média menor do que a da população masculina brasileira de 19 anos de 2016.

Testet_1condicao.R

Duas condições independentes

Teste t independente: Estatura e Sexo

Testar a igualdade das médias populacionais das estaturas masculina e feminina. Assumir que os desvios-padrão populacionais da VD nas duas condições são desconhecidos.

\[ \begin{cases} H_0: \mu_{\text{Estatura}}^{\text{M}} = \mu_{\text{Estatura}}^{\text{F}}\\ H_1: \mu_{\text{Estatura}}^{\text{M}} \ne \mu_{\text{Estatura}}^{\text{F}} \end{cases}\\ \alpha=0.05 \]

Estimar e classificar a significância prática da estimativa pontual da estatística de tamanho de efeito d de Cohen e obter os limites inferior e superior de seu intervalo de confiança de 95%.

Testet_independente.R

Duas condições dependentes

Teste t relacionado

Ler NewDrug.rds.

Dados <- readRDS("NewDrug.rds")
Dados <- data.frame(Dados$wide)

As medidas pulse e resp foram observadas em 12 participantes em três semanas consecutivas, sendo 6 participantes do grupo Placebo e 6 participantes do grupo New Drug.

Testar a igualdade das médias populacionais de pulse1 e pulse3 no grupo New Drug.

\[ \begin{cases} H_0: \mu_{\text{pulse1}}^{\text{New Drug}} = \mu_{\text{pulse3}}^{\text{New Drug}}\\ H_1: \mu_{\text{pulse1}}^{\text{New Drug}} \ne \mu_{\text{pulse3}}^{\text{New Drug}} \end{cases}\\ \alpha=0.05 \]

Estimar e classificar a significância prática da estimativa pontual da estatística de tamanho de efeito d de Cohen e obter os limites inferior e superior de seu intervalo de confiança de 95%.

Testet_relacionado.R

Teste t relacionado com dados no formato wide

Testet_relacionado_wide.R

Tabela <- ("
Participante Semana1 Semana2
a            1200    1100
b            1400    1200
c            1350    1250
d             950    1050
e            1400    1200
f            1150    1250
g            1300    1350
h            1325    1350
i            1425    1325
j            1500    1525
k            1250    1225
l            1150    1125
")
Dados <- read.table(textConnection(Tabela),header=TRUE)
Dados$Participante <- factor(Dados$Participante)
dif <- Dados$Semana2 - Dados$Semana1
alfa <- 0.05
boxplot(dif, ylab="Sodio Depois-Antes")

car::densityPlot(dif)

shapiro.test(dif)


    Shapiro-Wilk normality test

data:  dif
W = 0.92523, p-value = 0.3322

Dados_long <- reshape2::melt(Dados,
                             id.vars="Participante",
                             measure.vars=c("Semana1","Semana2"),
                             variable.name="Semana")
names(Dados_long) <- c("Participante", "Semana", "Sodio")
print(Dados_long)

   Participante  Semana Sodio
1             a Semana1  1200
2             b Semana1  1400
3             c Semana1  1350
4             d Semana1   950
5             e Semana1  1400
6             f Semana1  1150
7             g Semana1  1300
8             h Semana1  1325
9             i Semana1  1425
10            j Semana1  1500
11            k Semana1  1250
12            l Semana1  1150
13            a Semana2  1100
14            b Semana2  1200
15            c Semana2  1250
16            d Semana2  1050
17            e Semana2  1200
18            f Semana2  1250
19            g Semana2  1350
20            h Semana2  1350
21            i Semana2  1325
22            j Semana2  1525
23            k Semana2  1225
24            l Semana2  1125

print(psych::describeBy(Sodio~Semana, data=Dados_long, mat=2, digits=2))

       item  group1 vars  n    mean     sd median trimmed    mad  min  max
Sodio1    1 Semana1    1 12 1283.33 152.38 1312.5  1295.0 148.26  950 1500
Sodio2    2 Semana2    1 12 1245.83 129.61 1237.5  1237.5 148.26 1050 1525
       range  skew kurtosis    se
Sodio1   550 -0.61    -0.52 43.99
Sodio2   475  0.46    -0.49 37.42

print(psych::describe(dif))

   vars  n  mean     sd median trimmed    mad  min max range  skew kurtosis
X1    1 12 -37.5 103.63    -25     -35 111.19 -200 100   300 -0.22    -1.38
      se
X1 29.91

g1 <- ggplot2::ggplot(Dados_long, 
                      ggplot2::aes(x=Semana, 
                                   y=Sodio, 
                                   colour=Participante, 
                                   group=Participante)) +
  ggplot2::geom_line() + 
  ggplot2::geom_point(shape=21, 
                      fill="white") + 
  ggplot2::ylab("Sodio (mg/dia)") +
  ggplot2::ggtitle("Perfil do Participante") +
  ggplot2::theme_bw()
g1

Dados_long_ICIP <- summarySEwithin2(Dados_long, 
                                    measurevar="Sodio", 
                                    withinvars="Semana",
                                    idvar="Participante", 
                                    na.rm=TRUE, 
                                    conf.interval=1-alfa/2)


Anexando pacote: 'data.table'

O seguinte objeto é mascarado por 'package:DescTools':

    %like%


Anexando pacote: 'dplyr'

Os seguintes objetos são mascarados por 'package:data.table':

    between, first, last

Os seguintes objetos são mascarados por 'package:stats':

    filter, lag

Os seguintes objetos são mascarados por 'package:base':

    intersect, setdiff, setequal, union

print(Dados_long_ICIP)

   Semana  N    Sodio SodioNormed       sd       se       ci
1 Semana1 12 1283.333    1283.333 73.27563 21.15285 54.85131
2 Semana2 12 1245.833    1245.833 73.27563 21.15285 54.85131

g2 <- ggplot2::ggplot(Dados_long_ICIP, 
                      ggplot2::aes(x=Semana, 
                                   y=Sodio)) +
  ggplot2::geom_errorbar(width=.1, 
                         ggplot2::aes(ymin=Sodio-ci, 
                                      ymax=Sodio+ci)) +
  ggplot2::geom_point(shape=21, 
                      size=3, 
                      fill="white") +
  ggplot2::ylab("Sodio (mg/dia)") +
  ggplot2::ggtitle("IC95% Bonferroni intraparticipantes") +
  ggplot2::theme_bw()
g2

cat("Analise de significancia estatistica: valor p\n")

Analise de significancia estatistica: valor p

print(res <- t.test(dif,mu=0))


    One Sample t-test

data:  dif
t = -1.2536, df = 11, p-value = 0.236
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -103.3417   28.3417
sample estimates:
mean of x 
    -37.5

cat("Analise de significancia pratica: tamanho de efeito\n")

Analise de significancia pratica: tamanho de efeito

d <- lsr::cohensD(x=Dados$Semana2, 
                  y=Dados$Semana1,
                  method="paired")
cat("d de Cohen =", round(d,2), "\n")

d de Cohen = 0.36

es <- effectsize::interpret_cohens_d(d)
names(es) <- c("Tamanho de efeito: estimativa pontual")
print(es)

Tamanho de efeito: estimativa pontual 
                              "small" 
(Rules: cohen1988)

t <- res$statistic
F <- t^2
df <- res$parameter
eta2 <- F/(F+df)
cat("eta^2 = R^2 = ", eta2, "\n\n")

eta^2 = R^2 =  0.125

print(effectsize::interpret_eta_squared(eta2, rules="cohen1992"))

      t 
"small" 
(Rules: cohen1992)

Planejamento de teste t frequentista

Todo teste de hipótese nula frequentista consiste num sistema fechado entre o nível de significância, poder prospectivo, tamanho de efeito populacional e tamanho de amostra, i.e., se três deles são estipulados, o quarto elemento pode ser calculado.

G*Power consegue calcular cada um deles (https://www.psychologie.hhu.de/arbeitsgruppen/allgemeine-psychologie-und-arbeitspsychologie/gpower).

Teste t de Student frequentista

A função pwr::pwr.t.test foi adotada para o planejamento (seus resultados coincidem com os obtidos pelo G*Power). Recebe três dos seguintes quatro parâmetros, calculando o que não for fornecido (NULL):

d … tamanho de efeito populacional
power … poder prospectivo
sig.level … nível de significância
n … tamanho da amostra (cada condição)

Note que, nesta função, o tamanho de efeito populacional é dado pelo d de Cohen; d é distância padronizada entre \(H_0\) e \(H_1\). Conforme Cohen (1988), as magnitudes populacionais de significância prática pequena, intermediária e grande são, respectivamente 0.2, 0.5 e 0.8.

Os valores de poder prospectivo comumente usados são 0.8 e 0.9. Esta é a probabilidade de rejeitar corretamente \(H_0\).

O nível de confiança (\(1-\alpha\)) é a probabilidade de não rejeitar corretamente \(H_0\) (por isso, o nível de significância \(\alpha\) é a probabilidade de rejeitar incorretamente \(H_0\)).

Qual é o tamanho de amostra (\(n\)) para cada grupo se o nível de significância é igual a 0.05, poder prospectivo igual a 0.9 e o tamanho de efeito é intermediário?
Quais são os valores de poder prospectivo do teste t de Student frequentista para cada um dos três níveis de tamanho de efeito, nível de significância igual a 0.05 e 5 participantes em cada grupo?

Teste t de Welch frequentista

A função pwr::pwr.t.test realiza o planejamento do teste t de Student (homocedástico).

O código da função que calcula os tamanhos de amostra por grupo do teste t de Welch é de autoria de Jan & Shieh (2011).

Referências

Jan SL & Shieh G (2011) Optimal sample sizes for Welch’s test under various allocation and cost considerations. Behav Res Methods 43(4):1014-22. doi: 10.3758/s13428-011-0095-7.
Power calculation for two-sample Welch’s t test: https://stats.stackexchange.com/questions/413815/power-calculation-for-two-sample-welchs-t-test.

Duas condições independentes: tamanho de amostra muito grande

Teste t independente: idade de óbito e sexo

Testet_independente_DOBR2020.R

Ler o arquivo DOBR2020.rds.

Dados <- readRDS("DOBR2020.rds")
# sjPlot::view_df(Dados)

Questões

Realizar o teste t para comparar as médias populacionais de idade de óbito de adultos (18+) dos sexos feminino e masculino.
Interpretar o valor p e d Cohen. Qual faz mais sentido: valor p ou d de Cohen?

boxplot(data=Dados, 
        IdadeObito~SEXO)

print(psych::describeBy(IdadeObito~SEXO,
                        mat=2,
                        digits=6,
                        data=Dados[Dados$IdadeObito >= 18,]))

            item    group1 vars      n     mean       sd median  trimmed
IdadeObito1    1  Ignorado    1     89 69.22472 19.69123     73 70.53425
IdadeObito2    2 Masculino    1 842099 65.04460 18.55398     68 66.39839
IdadeObito3    3  Feminino    1 661045 72.18622 16.79651     75 73.60812
                mad min max range      skew  kurtosis       se
IdadeObito1 14.8260  18 120   102 -0.622505  0.218859 2.087266
IdadeObito2 17.7912  18 125   107 -0.580035 -0.266404 0.020219
IdadeObito3 16.3086  18 120   102 -0.739004  0.223625 0.020659

alfa <- 0.05
gplots::plotmeans(data=Dados,
                  subset=(SEXO=="Masculino" | SEXO=="Feminino") 
                         & IdadeObito >= 18,
                  IdadeObito~SEXO,
                  p=1-alfa/2,
                  connect=FALSE,
                  barcol="black",
                  main="IC95% Bonferroni")

Registered S3 method overwritten by 'gplots':
  method         from     
  reorder.factor DescTools

tt <- t.test(data=Dados,
             subset=(SEXO=="Masculino" | SEXO=="Feminino")  
                    & IdadeObito >= 18,
             IdadeObito~SEXO,
             conf.level=1-alfa)
print(tt)


    Welch Two Sample t-test

data:  IdadeObito by SEXO
t = -247.06, df = 1473015, p-value < 2.2e-16
alternative hypothesis: true difference in means between group Masculino and group Feminino is not equal to 0
95 percent confidence interval:
 -7.198278 -7.084966
sample estimates:
mean in group Masculino  mean in group Feminino 
               65.04460                72.18622

print(try(shapiro.test(Dados[Dados$SEXO=="Masculino"
                             & Dados$IdadeObito >= 18,
                             "IdadeObito"])))

Error in shapiro.test(Dados[Dados$SEXO == "Masculino" & Dados$IdadeObito >=  : 
  tamanho da amostra deve ser entre 3 e 5000
[1] "Error in shapiro.test(Dados[Dados$SEXO == \"Masculino\" & Dados$IdadeObito >=  : \n  tamanho da amostra deve ser entre 3 e 5000\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in shapiro.test(Dados[Dados$SEXO == "Masculino" & Dados$IdadeObito >=     18, "IdadeObito"]): tamanho da amostra deve ser entre 3 e 5000>

print(try(shapiro.test(Dados[Dados$SEXO=="Masculino"
                             & Dados$IdadeObito >= 18,
                             "IdadeObito"])))

Error in shapiro.test(Dados[Dados$SEXO == "Masculino" & Dados$IdadeObito >=  : 
  tamanho da amostra deve ser entre 3 e 5000
[1] "Error in shapiro.test(Dados[Dados$SEXO == \"Masculino\" & Dados$IdadeObito >=  : \n  tamanho da amostra deve ser entre 3 e 5000\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in shapiro.test(Dados[Dados$SEXO == "Masculino" & Dados$IdadeObito >=     18, "IdadeObito"]): tamanho da amostra deve ser entre 3 e 5000>

hmc <- fligner.test(IdadeObito~SEXO, data = Dados)
print(hmc)


    Fligner-Killeen test of homogeneity of variances

data:  IdadeObito by SEXO
Fligner-Killeen:med chi-squared = 4260.2, df = 2, p-value < 2.2e-16

dC <- psych::cohen.d(data=subset(Dados, 
                                 subset=(SEXO=="Masculino" | 
                                         SEXO=="Feminino") 
                                        & IdadeObito >= 18),
                     IdadeObito~SEXO,
                     alpha=alfa)
cat("IC95%(d de Cohen) =", dC$cohen.d, "\n")

IC95%(d de Cohen) = 0.3979067 0.4011593 0.4044117

dCp <- dC$cohen.d[2]
names(dC) <- c("Effect size")
es <- effectsize::interpret_cohens_d(dCp, # estimativa pontual
                                     rules = "sawilowsky2009")
print(es)

[1] "small"
(Rules: sawilowsky2009)

Dadoss <- subset(Dados, 
                 subset=(SEXO=="Masculino" | 
                         SEXO=="Feminino") 
                 & IdadeObito >= 18)
Dadoss$SEXO <- droplevels(Dadoss$SEXO)
boxplot(data=Dadoss, 
        IdadeObito~SEXO)

dC <- effsize::cohen.d(IdadeObito~SEXO, 
                       data=Dadoss, 
                       paired=FALSE,  
                       na.rm=TRUE) # Classificação de Cohen (1988)
print(dC)


Cohen's d

d estimate: -0.401159 (small)
95 percent confidence interval:
     lower      upper 
-0.4044115 -0.3979065

es <- effectsize::interpret_cohens_d(d=dC$estimate, rules="sawilowsky2009")
names(es) <- c("Tamanho de efeito: estimativa pontual")
print(es)

Tamanho de efeito: estimativa pontual 
                              "small" 
(Rules: sawilowsky2009)

fit <- oneway.test(data=Dadoss,
                   IdadeObito~SEXO)
print(fit)


    One-way analysis of means (not assuming equal variances)

data:  IdadeObito and SEXO
F = 61039, num df = 1, denom df = 1473015, p-value < 2.2e-16

eta2 <- effectsize::eta_squared(fit,
                                ci = 1-alfa,
                                alternative = "two.sided")

For one-way between subjects designs, partial eta squared is equivalent
  to eta squared. Returning eta squared.

`var.equal = FALSE` - effect size is an approximation.

print(eta2, digits = 6)

# Effect Size for ANOVA

Eta2     |               95% CI
-------------------------------
0.039789 | [0.039179, 0.040404]

es <- effectsize::interpret_eta_squared(eta2$Eta2,
                                        rules="cohen1992")
names(es) <- c("Tamanho de efeito: estimativa pontual")
print(es)

Tamanho de efeito: estimativa pontual 
                              "small" 
(Rules: cohen1992)

Prática R: Teste t

MSP4112: Estatística Aplicada à Medicina

José Siqueira (siqueira@usp.br)
Paulo Silveira (silveira@usp.br)

15 agosto 2025 13:55h

Material

Introdução

Teste t

Uma condição

Teste t bilateral: Estatura

Teste t unilateral: Estatura

Duas condições independentes

Teste t independente: Estatura e Sexo

Duas condições dependentes

Teste t relacionado

Teste t relacionado com dados no formato wide

Planejamento de teste t frequentista

Teste t de Student frequentista

Teste t de Welch frequentista

Referências

Duas condições independentes: tamanho de amostra muito grande

Teste t independente: idade de óbito e sexo

Questões

Prática R: Teste t

MSP4112: Estatística Aplicada à Medicina

José Siqueira (siqueira@usp.br) Paulo Silveira (silveira@usp.br)

15 agosto 2025 13:55h

Material

Introdução

Teste t

Uma condição

Teste t bilateral: Estatura

Teste t unilateral: Estatura

Duas condições independentes

Teste t independente: Estatura e Sexo

Duas condições dependentes

Teste t relacionado

Teste t relacionado com dados no formato wide

Planejamento de teste t frequentista

Teste t de Student frequentista

Teste t de Welch frequentista

Referências

Duas condições independentes: tamanho de amostra muito grande

Teste t independente: idade de óbito e sexo

Questões

José Siqueira (siqueira@usp.br)
Paulo Silveira (silveira@usp.br)