MAN5728: Tópicos Essenciais da Bioestatística

Módulo 3: Testagem de Hipótese Nula
Testes qui-quadrado, t e ANOVA

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

Pacotes R

options(warn=-1)
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(MBESS, warn.conflicts=FALSE))
suppressMessages(library(RVAideMemoire, warn.conflicts=FALSE))
suppressMessages(library(corrplot, warn.conflicts=FALSE))
suppressMessages(library(gplots, warn.conflicts=FALSE))
suppressMessages(library(RcmdrMisc, warn.conflicts=FALSE))
suppressMessages(library(rcompanion, warn.conflicts=FALSE))
suppressMessages(library(Rmisc, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(lsr, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(rstatix, warn.conflicts=FALSE))
suppressMessages(library(jmv, warn.conflicts=FALSE))
source("summarySEwithin2.R")

Material

• HTML de R Markdown em RPubs
• Prof. José Siqueira: ResearchGate

R

• RStudio Cloud

Conteúdo

• Teste de hipótese nula: conceito
• Teste qui-quadrado
• Testes paramétricos mais utilizados na área da saúde:
  • teste t
  • ANOVA

Ler planilha

Os dados da planilha Excel Biometria_FMUSP.xlsx foram coletados pelos docentes do curso de Medicina da FMUSP dos estudantes do segundo ano de uma mesma disciplina em três anos consecutivos.

As variáveis do arquivo são:

• ID: idenficador do(a) estudante
  
• Ano da coleta dos dados: 1, 2, 3
  
• Turma: A, B
  
• Sexo: Feminino, Masculino
• Mao: Destro, Canhoto, Ambidestro
  
• TipoSang: A+, A-, …
• ABO: A, B, AB, O
• Rh: +, -
• AtivFisica: nível de atividade física
• Sedentarismo: Não, Sim
  
• MCT: massa corporal total (kg)
• Estatura: cm

Dados <- readxl::read_excel(path="Biometria_FMUSP.xlsx",
                            sheet="dados",
                            na=c("NA","na","Na","nA"))
Dados$MCT[Dados$MCT==max(Dados$MCT,na.rm=TRUE)] <- NA
Dados$Estatura[Dados$Estatura==min(Dados$Estatura,na.rm=TRUE)] <- NA
Dados$ID <- factor(Dados$ID)
Dados$Ano <- factor(Dados$Ano)
Dados$Turma <- factor(Dados$Turma)
Dados$Sexo <- factor(Dados$Sexo)
Dados$Mao <- factor(Dados$Mao)
Dados$TipoSang <- factor(Dados$TipoSang)
Dados$ABO <- factor(Dados$ABO)
Dados$Rh <- factor(Dados$Rh)
Dados$AtivFisica <- factor(Dados$AtivFisica,
                           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.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")

Teste de hipótese: conceito

Teste z bilateral para uma condição: não rejeição da hipótese nula

A distribuição da estatura do homem brasileiro de 19 anos de 2016 é normal com desvio-padrão \(\sigma\) = 7 cm.

Testar se a média populacional \(\mu\) é igual a 177 cm hipotetizada.

Quatro participantes deste grupo de 2020 tiveram suas estaturas medidas, cujos valores em centímetro são: 169, 174, 175, 186.

Teste

Média populacional hipotetizada: \(\mu\) = 177 cm

Suposições

1. \(n = 4\) observações independentes: 169, 174, 175, 186

2. Se \(n<30\), estatura tem distribuição normal; se \(n\ge30\), estatura tem distribuição qualquer (Teorema Central do Limite)

3. Desvio-padrão \(\sigma\) = 7 cm conhecido

4. Teste bilateral

5. Nível de confiança adotado de 95% (ou nível de significância de 5%)

Hipóteses

\(H_{0}: μ = 177\)

\(H_{1}: μ \ne 177\)

ou equivalentemente

\(H_{1}: μ < 177 \; \text{ou} \; μ > 177\)

Estatísticas

\(\bar{x} = \dfrac{175+186+169+174}{4} = 176\)

\(ep = \dfrac{\sigma}{\sqrt{n}} = \dfrac{7}{\sqrt{4}}= 3.5\)

\(\text{IC}^{95\%}(\mu) = [176 \pm 1.96\times3.5] = [169.14, 182.86]\)

Estatística de teste \(z = \dfrac{\bar{x}-177}{ep}=\dfrac{176−177}{3.5}=−0.29\)

Estatística de tamanho de efeito \(d=\dfrac{|\bar{x}-177|}{\sigma}=\dfrac{|176-177|}{7}=0.14\) (muito pequeno)

Decisão

Critério do valor crítico: Como \(|z| = 0.29 < 1.96\), não rejeitar \(H_{0}\)

Critério do \(IC\): Como \(\text{IC}^{95\%}\) contém 177, não rejeitar \(H_{0}\)

Critério do valor p: Como o valor p bilateral = 0.77 é maior que 0.05, não rejeitar \(H_{0}\)

alfa <- 0.05
estatura <- c(169, 174, 175, 186)
shapiro.test(estatura)

    Shapiro-Wilk normality test

data:  estatura
W = 0.91181, p-value = 0.492

DescTools::ZTest(x=estatura, 
                 sd_pop=7,
                 mu=177,
                 alternative="two.sided",
                 conf.level=1-alfa,
                 na.rm=TRUE)

    One Sample z-test

data:  estatura
z = -0.28571, Std. Dev. Population = 7, p-value = 0.7751
alternative hypothesis: true mean is not equal to 177
95 percent confidence interval:
 169.1401 182.8599
sample estimates:
mean of x 
      176 

effectsize::interpret_cohens_d(d=0.14, rules="cohen1988")

[1] "very small"
(Rules: cohen1988)

Teste z bilateral para uma condição: rejeição da hipótese nula

A distribuição da estatura do homem brasileiro de 19 anos de 2016 é normal com desvio-padrão \(\sigma\) = 7 cm.

Testar se a média populacional \(\mu\) é igual a 177 cm hipotetizada.

Mil participantes deste grupo de 2020 tiveram suas estaturas medidas, cujos valores em centímetro são: 169, 174, 175, 186, …

Teste

Média populacional hipotetizada: \(\mu\) = 177 cm

Suposições

1. \(n = 1000\) observações independentes: 169, 174, 175, 186, …

2. Se \(n<30\), estatura tem distribuição normal; se \(n\ge30\), estatura tem distribuição qualquer (Teorema Central do Limite)

3. Desvio-padrão \(\sigma\) = 7 cm conhecido

4. Teste bilateral

5. Nível de confiança adotado de 95% (ou nível de significância de 5%)

Hipóteses

\(H_{0}: μ = 177\)

\(H_{1}: μ \ne 177\)

Estatísticas

\(\bar{x} = \dfrac{175+186+169+174+\cdots}{1000} = 176\)

\(ep = \dfrac{\sigma}{\sqrt{n}} = \dfrac{7}{\sqrt{1000}}= 0.22\)

\(\text{IC}^{95\%}(\mu) = [176 \pm 1.96\times 0.22] = [175.57, 176.43]\)

Estatística de teste \(z = \dfrac{\bar{x}-177}{ep}=\dfrac{176−177}{0.22}=−4.52\)

Estatística de tamanho de efeito \(d=\dfrac{|\bar{x}-177|}{\sigma}=\dfrac{|176-177|}{7}=0.14\) (muito pequeno)

Decisão

Critério do valor crítico: Como \(|z| = 4.52 > 1.96\), rejeitar \(H_{0}\)

Critério do \(IC\): Como \(\text{IC}^{95\%}\) não contém 177, rejeitar \(H_{0}\)

Critério do valor p: Como o valor p bilateral = \(6.2\times10^{-6}\) é menor que 0.05, rejeitar \(H_{0}\)

alfa <- 0.05
set.seed(123)
s <- rnorm(999, 175.9, 7)
estatura <- c(s, 176*1000-999*mean(s))
shapiro.test(estatura)

    Shapiro-Wilk normality test

data:  estatura
W = 0.99837, p-value = 0.4712

DescTools::ZTest(estatura, 
                 sd_pop=7,
                 mu=177,
                 alternative="two.sided",
                 conf.level=1-alfa,
                 na.rm=TRUE)

    One Sample z-test

data:  estatura
z = -4.5175, Std. Dev. Population = 7, p-value = 6.256e-06
alternative hypothesis: true mean is not equal to 177
95 percent confidence interval:
 175.5661 176.4339
sample estimates:
mean of x 
      176 

effectsize::interpret_cohens_d(d=0.14, rules="cohen1988")

[1] "very small"
(Rules: cohen1988)

“[…] there is always a large enough sample size \(n\) for which any simple null hypothesis \(H_{0}: \mu=\mu_{0}\) will be rejected by a frequentist \(\alpha\)-significance level test.”

Spanos, 2014, p. 646

Análise estatística inferencial

• Significância estatística: valor p (existência de efeito)
• Significância prática: \(d\) de Cohen, \(\eta^{2}\) de Cohen e V de Cramer (tamanho de efeito)

Erro amostral

“Constitui um dos problemas enfrentados quando conduzimos uma pesquisa o fato de não sabermos qual é o padrão existente na população de interesse.

De fato, o motivo de realizarmos a pesquisa é, em primeiro lugar, determinar esse padrão.

Você precisa estar ciente de que, algumas vezes, devido ao erro amostral, obteremos padrões nas amostras que não refletem de forma acurada a população de onde as amostras foram retiradas.

Assim, precisamos de um algum meio para avaliar a probabilidade de que a amostra selecionada seja um retrato fiel da população.

Os testes estatísticos nos auxiliam nesta decisão, mas isso ocorre de uma forma não de todo intuitiva.”

Dancey & Reidy, 2019

Efeito populacional

• Associação/ concordância entre variáveis na população
• Diferença entre condições na população

Hipóteses nula e alternativa

• Hipótese nula/ \(H_{0}\)
  • A hipótese nula sempre declara que NÃO existe efeito na população.
• Hipótese de pesquisa/ alternativa/ \(H_{1}\)/ \(H_{a}\)
  • A hipótese de pesquisa declara que existe efeito na população.

Teste estatístico de hipótese nula

“Em Estatística, assim como em Economia, Psicologia, Medicina e Direito, o problema de testagem de hipótese nula envolve a mediação de objetivos conflitantes.

Há uma analogia legal interessante: no julgamento de um crime de homicídio, pede-se ao júri que decida entre a hipótese nula \(H_{0}\): O réu não é culpado e a hipótese alternativa \(H_{1}\): o réu é culpado.

O sistema judiciário não é perfeito.

Comete-se um erro do tipo I condenando-se um inocente. A probabilidade do erro do tipo I é chamado de nível de significância do teste (\(\alpha\)).

O nível de confiança do teste (probabilidade de absolver o inocente) é \(1 – \alpha\).

O erro do tipo II consiste em absolver um culpado. A probabilidade do erro do tipo II é o \(\beta\).

O poder do teste (probabilidade de condenar um culpado) é igual a \(1 – \beta\).

A advertência do juiz ao júri, de que o crime “a culpa pelo crime de homicídio deve ser provada além de qualquer dúvida razoável” significa que \(\alpha\) deve ser muito pequena.

Tem havido muitas reformas legais (e.g., limitar o poder da polícia para obter confissão) elaboradas a fim de reduzir \(\alpha\).

Porém, essas mesmas reformas têm contribuído para aumentar \(\beta\).

Não há meios de reduzir \(\alpha\) a 0 (impossibilidade total de condenar um inocente) sem elevar \(\beta\) a 1 (tender a libertar todos os culpados, invalidando o julgamento).

A única maneira de reduzir \(\alpha\) e \(\beta\) concomitantemente é aumentar a evidência, i.e., aumentar o tamanho da amostra.”

Wonnacott & Wonnacott, 1981

“Quais das seguintes situações representam um erro do tipo I e quais um erro do tipo II?

1. Você verificou, em um estudo, a existência de um relacionamento entre a quantidade de chá ingerida por dia e a quantidade de dinheiro ganho na loteria. Você conclui que, para ganhar na loteria, deve beber muito chá.
   • Solução: \(H_{0}\): Não há relacionamento entre quantidade de chá ingerida e quantidade de dinheiro ganho na loteria. Um estudo concluiu que essa relação existe, pois \(p<0.05\). Na realidade essa relação não existe! Portanto, \(H_{0}\) foi rejeitada, sendo que ela é verdadeira, i.e., o erro do tipo I foi cometido.
2. Você verificou, em um estudo, que não existe diferença entre as velocidades das tartarugas e dos guepardos. Você conclui que as tartarugas são tão rápidas quanto os guepardos.
   • Solução: \(H_{0}\): Não há diferença entre as velocidades das tartarugas e guepardos. Um estudo concluiu que essa relação não existe, pois \(p>0.05\). Na realidade essa diferença existe! Portanto, \(H_{0}\) não foi rejeitada, sendo que ela é falsa, i.e., o erro do tipo II foi cometido.
3. Você verificou, em um estudo, que existe um relacionamento entre estilo de vida e renda anual. No entanto, em virtude de a probabilidade associada com o relacionamento ser de \(0.5\), você conclui que não existe relacionamento entre as duas variáveis.
   • Solução: \(H_{0}\): Não há relacionamento entre estilo de vida e renda anual. Um estudo concluiu que essa relação não existe, pois \(p>0.05\); Na realidade esse relacionamento existe! Portanto, \(H_{0}\) não foi rejeitada sendo que ela é falsa, i.e., o erro do tipo II foi cometido.”

Dancey & Reidy, 2019

Teste de hipótese nula em ação

Understanding Statistical Power and Significance Testing: An interactive visualization

Decisão estatística vs. Estado da natureza

Estatística de teste

Wikipedia: Test statistic

Valor p

“O valor p é a probabilidade de que a estatística de teste seja igual ou mais extrema que o valor observado na direção prevista pela hipótese alternativa (\(H_{1}\)), presumindo que a hipótese nula (\(H_{0}\)) é verdadeira.”

Agresti & Finlay, 2012, p. 171

“Mínimo \(\alpha\) que rejeita \(H_{0}\).”

Wonnacott & Wonnacott, 1981

“Valores p e métodos de testagem de hipótese nula são frequentemente mal utilizados em pesquisas clínicas.”

Mark et al., 2016, p. 1048

Dancey & Reidy, 2019

Por que estabelecer \(\alpha = 0.05\)?

“Não recomendamos abandonar o valor p, nem reduzir arbitrariamente o limiar de significância [\(\alpha = 0.05\)].”

Di Leo & Sardanelli, 2020, p. 6

Testes unilateral e bilateral

• Quando a direção do relacionamento ou da diferença (efeito) é especificada, então o teste é unilateral; caso contrário, é bilateral.
• Em geral (mas nem sempre), se você tiver obtido um valor p para um teste bilateral e quiser saber o valor mínimo correspondente para o teste unilateral, então: \(p_{\text{unilateral}} ≥ p_{\text{bilateral}}/2\)
• Observe que o que deve ser dobrado ou dividido por 2 não é a estatística de teste (e.g.: \(z\) ou \(t\)).
• \(H_{0}: \mu \ge \mu_{0}\) vs. \(H_{1}: \mu < \mu_{0}\) é equivalente a \(H_{0}: \mu = \mu_{0}\) vs. \(H_{1}: \mu < \mu_{0}\)
  • Demonstração em Gatás, 1978, p. 220-3.

“Testes unicaudais raramente devem ser usados para pesquisa básica ou aplicada em ecologia, comportamento animal ou qualquer outra ciência.”

Lombardi & Hurlbert, 2009, p. 447

“Testes unicaudais também se mostram úteis em uma área de interesse estatístico chamada ‘não inferioridade’ e testes de ‘equivalência’.”

Lombardi & Hurlbert, 2009, p. 462

Teste de hipótese nula & Intervalo de confiança

• Os critérios de valor p e de intervalo de confiança são equivalentes.

Críticas contra o teste de hipótese nula

“Não recomendamos abandonar o valor p, nem reduzir arbitrariamente o limiar de significância [\(\alpha = 0.05\)].”

Di Leo & Sardanelli, 2020, p. 6

“A testagem da hipótese nula é a abordagem dominante na Psicologia e Medicina.

Apesar das críticas à testagem da hipótese nula, isso não significa que tal abordagem deve ser abandonada completamente.

Ao invés disso, devemos ter um entendimento completo de seu significado para podermos nos beneficiar desta tecnologia da decisão.

Além do valor p, é importante usar o intervalo de confiança e de tamanho de efeito.”

Dancey & Reidy, 2019

Significância prática

“Mesmo efeitos muito pequenos poderão apresentar significância estatística quando o tamanho da amostra for bem grande.

Para determinar a significância prática, a melhor abordagem consiste em obter uma medida do tamanho do efeito, sendo que essa medida não depende do tamanho da amostra.

E.g.: o coeficiente de correlação de Pearson amostral mede a intensidade da associação linear entre duas variáveis quantitativas e não depende do tamanho da amostra.”

Dancey & Reidy, 2019

Interpretação errônea do valor p

• Muitos pesquisadores sem experiência em estatística (e mesmo aqueles com alguma) equiparam o valor p com o tamanho do efeito, i.e., quanto menor o valor p, mais forte seria, por exemplo, o relacionamento entre duas variáveis; talvez, de fato, quanto mais forte o relacionamento, mais baixo o valor p, mas não significa que isso necessariamente ocorrerá.
• O valor p não é a probabilidade de que a hipótese nula seja verdadeira; de fato, não sabemos qual é a probabilidade de que a hipótese nula seja verdadeira.
• \(1 – p\) não é a probabilidade de que a hipótese alternativa seja verdadeira; de fato, não sabemos qual é a probabilidade de que a hipótese alternativa seja verdadeira.

Dancey & Reidy, 2019

Replicação

“A replicação é uma das pedras angulares da ciência. Se você observa um fenômeno uma vez, então pode ter sido por acaso; se o observa duas, três ou mais vezes, pode estar começando a aprender algo sobre o fenômeno estudado.

Se o seu estudo foi o primeiro neste assunto, é sensato que você trate os resultados com certo grau de cautela.”

Dancey & Reidy, 2019

“Infelizmente, nenhuma ferramenta ou técnica estatística pode garantir o acesso ao caminho mais curto para a verdade.

Experimentos repetidos, como Fisher reconheceu anos atrás, podem ser a melhor maneira de domar grande parte da bagunça.

Se isso não for viável, e muitas vezes não é, a medicina ainda pode realizar muito fazendo uso pragmático e bem fundamentado das evidências que possui.”

Mark et al., 2016, p. 1054

Planejamento do estudo: falando grego

Coelho et al. 2008, Tabela 4.1

“A análise de poder é uma das ferramentas mais fundamentais que os pesquisadores podem usar ao planejar estudos.”

Perugini et al., 2018

“Resumo

Antecedentes:

A crença continua generalizada de que os estudos de pesquisa médica devem ter poder estatístico de pelo menos 80% para serem cientificamente sólidos, e os revisores geralmente questionam se o poder é alto o suficiente.

Discussão:

Este requisito e os métodos para atendê-lo apresentam falhas graves.

Notavelmente, a verdadeira natureza de como o tamanho da amostra influencia o valor científico ou prático projetado de um estudo impede qualquer designação abrangente de <80% de poder como “inadequada”.

Além disso, os cálculos padrão são inerentemente não confiáveis, e focar apenas no poder negligencia os resultados mais importantes de um estudo concluído: estimativas e intervalos de confiança.

As convenções atuais prejudicam o processo de pesquisa de várias maneiras: promovendo má interpretação de estudos concluídos, corroendo a integridade científica, dando poder arbitrário aos revisores, inibindo a inovação, pervertendo padrões éticos, desperdiçando esforços e desperdiçando dinheiro.

A pesquisa médica se beneficiaria de abordagens alternativas, incluindo o valor estabelecido dos métodos de informação, escolhas simples baseadas no custo ou viabilidade que foram recentemente justificadas, análises de sensibilidade que examinam uma série significativa de descobertas possíveis e seguindo estudos análogos anteriores.

Para promover abordagens mais racionais, o treinamento em pesquisa deve abranger as questões apresentadas aqui, os revisores devem ser extremamente cuidadosos antes de levantar questões de tamanho de amostra “inadequado” e relatórios de estudos concluídos não devem discutir poder.

Resumo:

Convenções e expectativas comuns em relação ao tamanho da amostra são profundamente falhas, causam sérios danos ao processo de pesquisa e devem ser substituídas por alternativas mais racionais.”

Bacchetti, 2010, p. 8

\(p>0.05\) e análise de poder retrospectivo

“Evite o erro “ausência de evidência = evidência de ausência” [ou aceitar \(H_{0}\) se \(p>0.05\)].

Se um teste de hipótese tem um valor p acima do nível \(\alpha\) escolhido, não é correto inferir que não houve efeito ou diferença entre os grupos.”

Althouse et al., 2021, p. e73

“Cálculos de poder post hoc usando o tamanho do efeito observado são uma simples transformação do valor p do estudo e nunca devem ser usados porque não respondem a uma pergunta importante.

Se um cálculo de poder não foi realizado, então pode ser aceitável relatar um “cálculo de poder” ilustrando que o tamanho da amostra disponível foi suficiente para abordar a questão de interesse do estudo.

Observe que isso deve ser baseado em algum tamanho mínimo de efeito de interesse para a inferência primária, não nos dados observados.”

Althouse et al., 2021, p. e87

“Repetimos pontos relatados por outros autores que o uso do poder estimado para interpretar resultados é logicamente inconsistente e redundante com o valor p, e esse poder estimado pode ser tendencioso e altamente não confiável.”

Gerard et al., 1998, p. 802

Análise de poder retrospectivo (plug in) de ANOVA unifatorial independente balanceada (Gerard et al., 1998)

k <- 3 # numero de condicoes independentes
n <- 20 # numero de observacoes independentes de cada condicao
alfa <- 0.05
dfn <- k - 1
dfd <- k*(n - 1)
Fcrt <- qf(1-alfa, dfn, dfd)
Fobs <- Fcrt*0.99 
eta2 <- dfn*Fobs/(dfn*Fobs+dfd)
eta2lims <- MBESS::ci.pvaf(Fobs, dfn, dfd, k*n, 1-alfa)
f2 <- eta2/(1-eta2)
f2.ll <- eta2lims$Lower.Limit.Proportion.of.Variance.Accounted.for/
  (1-eta2lims$Lower.Limit.Proportion.of.Variance.Accounted.for)
f2.ul <- eta2lims$Upper.Limit.Proportion.of.Variance.Accounted.for/
  (1-eta2lims$Upper.Limit.Proportion.of.Variance.Accounted.for)
ncp.p <- dfd*f2 # ou dfn*Fobs
ncp.p.ll <- dfd*f2.ll 
ncp.p.ul <- dfd*f2.ul 
cat(paste("n =", k*n, "\tk =", k, "\tFcrt =", round(Fcrt,2),"\tFobs =", round(Fobs,2),"\n"))

n = 60  k = 3   Fcrt = 3.16     Fobs = 3.13 

poder.p <- 1-pf(Fcrt,dfn, dfd, ncp.p)
cat(paste("Estimativa pontual do poder retrospectivo =", round(poder.p,3),"\n"))

Estimativa pontual do poder retrospectivo = 0.579 

poder.p.ll <- 1-pf(Fcrt,dfn, dfd, ncp.p.ll)
cat(paste("Limite inferior do IC95% do poder retrospectivo =", round(poder.p.ll,3),"\n"))

Limite inferior do IC95% do poder retrospectivo = 0.05 

poder.p.ul <- 1-pf(Fcrt,dfn, dfd, ncp.p.ul)
cat(paste("Limite superior do IC95% do poder retrospectivo =", round(poder.p.ul,3),"\n"))

Limite superior do IC95% do poder retrospectivo = 0.967 

Teste qui-quadrado

Os testes baseados na estatística qui-quadrado podem ser divididos em dois tipos

1. Teste de aderência uma variável nominal com duas ou mais categorias: teste de frequências hipotetizadas;

2. Teste de independência entre duas variáveis nominais com duas ou mais categorias.

Teste qui-quadrado de aderência por bootstrapping

O teste de aderência nos permite descobrir se um conjunto de frequências observadas difere de um outro conjunto de frequências hipotetizadas.

Suposição

Independência das observações: cada observação da unidade observacional deve ser alocada em apenas um categoria da variável nominal.

Siegel & Castellan Jr, 1988, p. 46

Distribuição do tipo sanguíneo no Brasil

Wikipedia: Distribuição do tipo sanguíneo por país

• ABO_cluster.R

dados$Cluster: 1

[1] “Indígenas” “Bascos” “Belgas”
[4] “Brasileiros” “Povo_san” “Neerlandeses”
[7] “Ingleses” “Georgianos” “Islandeses”
[10] “Irlandeses” “Italianos_Milão” “Māori”
[13] “Escoceses” “Sul-africanos” “Americanos_brancos”

• ABORh_cluster.R

dados$Cluster: 5

[1] “Austrália” “Brasil” “Canadá” “Grécia”

[5] “Itália” “Líbano” “Lituânia” “Países Baixos”

[9] “Nova Zelândia” “Somália” “Sri Lanka” “Reino Unido”

tipsang.freq <- as.numeric(table(Dados$TipoSang))
tipsang.prob <- c(0.08, 0.34, 0.005, 0.025,
                  0.02, 0.08, 0.09, 0.36)
tslab <- c("A-", "A+", "AB-", "AB+", 
           "B-", "B+", "O-", "O+")
plot(tipsang.prob, type="p", lty=1, 
     lwd=1, pch=2, col="black", 
     xaxt = "n", xlab='Tipo Sanguineo', ylab="Probabilidade")
axis(1, at=1:length(tslab), labels=tslab)
lines(tipsang.freq/sum(tipsang.freq), pch=21, 
      type="p", lty=1, lwd=1, col="black")
legend("top", legend=c("Prob.Hipotet","Prob.Observ"), 
       pch=c(2,21), bty="n")

cat(tslab,"\n")

A- A+ AB- AB+ B- B+ O- O+ 

print(tipsang.freq)

[1]  49 171   2  40  34  65  45 132

print(round(tipsang.freq/sum(tipsang.freq), 3))

[1] 0.091 0.318 0.004 0.074 0.063 0.121 0.084 0.245

print(tipsang.prob)

[1] 0.080 0.340 0.005 0.025 0.020 0.080 0.090 0.360

round(ICM <- DescTools::MultinomCI(x=tipsang.freq,                                   conf.level=1-alfa/length(unique(Dados$TipoSang))),4)

        est lwr.ci upr.ci
[1,] 0.0911 0.0353 0.1494
[2,] 0.3178 0.2621 0.3761
[3,] 0.0037 0.0000 0.0620
[4,] 0.0743 0.0186 0.1326
[5,] 0.0632 0.0074 0.1215
[6,] 0.1208 0.0651 0.1791
[7,] 0.0836 0.0279 0.1419
[8,] 0.2454 0.1896 0.3037

ICM <- as.data.frame(ICM)
print(ICM, digits=3)

      est  lwr.ci upr.ci
1 0.09108 0.03532  0.149
2 0.31784 0.26208  0.376
3 0.00372 0.00000  0.062
4 0.07435 0.01859  0.133
5 0.06320 0.00743  0.121
6 0.12082 0.06506  0.179
7 0.08364 0.02788  0.142
8 0.24535 0.18959  0.304

ICM <- cbind(ICM,sort(unique(Dados$TipoSang)))
names(ICM) <- c("est", "lwr.ci", "upr.ci", "ts")
ggplot2::ggplot(ICM, 
                ggplot2::aes(x=ts,
                             y=est)) + 
  ggplot2::geom_errorbar(ggplot2::aes(ymin=lwr.ci, 
                                      ymax=upr.ci), 
                         width=.5) +
  ggplot2::geom_point(size=1,
                      shape=19) +
  ggplot2::xlab("Tipo Sanguineo") +
  ggplot2::ylab("Proporção") +
  ggplot2::ggtitle("I95% Bonferroni de MCT") +
  ggplot2::theme_bw()

RVAideMemoire::multinomial.multcomp(tipsang.freq, p.method = "bonferroni")

    Pairwise comparisons using exact binomial tests 

data:  tipsang.freq 

    2       34      40      45      49      65      132  
34  5.4e-07 -       -       -       -       -       -    
40  1.2e-08 1.000   -       -       -       -       -    
45  4.5e-10 1.000   1.000   -       -       -       -    
49  3.3e-11 1.000   1.000   1.000   -       -       -    
65  8.6e-16 0.067   0.525   1.000   1.000   -       -    
132 < 2e-16 2.2e-13 3.1e-11 1.1e-09 1.5e-08 5.8e-05 -    
171 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.4e-15 9.9e-11 0.808

P value adjustment method: bonferroni 

out <- chisq.test(x=tipsang.freq, 
                  p=tipsang.prob,
                  simulate.p.value=TRUE, 
                  B=1e6) 
df <- chisq.test(tipsang.freq, p=tipsang.prob)$parameter
X2 <- out$statistic
p <- out$p.value
n <- sum(chisq.test(tipsang.freq, p=tipsang.prob)$observed)
cat("X^2 = ", X2 , " df = ",  df, " p = ", p, "\n", sep="")

X^2 = 135.472 df = 7 p = 9.99999e-07

cat("Heuristica de significancia: X^2/df > 2\n", round(X2/df,2), "\n")

Heuristica de significancia: X^2/df > 2
 19.35 

cat("Residuos estandardizados\n", round(out$stdres,2), "\n")

Residuos estandardizados
 0.95 -1.08 -0.42 7.33 7.16 3.49 -0.52 -5.54 

V <- sqrt((X2/n)/df)
cat("V de Cramer = ", round(V,2), "\n")

V de Cramer =  0.19 

# cat("Grau de nao aderencia:\n")
print(effectsize::interpret_cramers_v(V))

X-squared 
  "small" 
(Rules: funder2019)

Teste qui-quadrado de independência por bootstrapping

O teste qui-quadrado de independência de Pearson permite que se descubra se existe um relacionamento ou efeito de interação entre duas variáveis nominais com duas ou mais categorias.

Suposição

Cada observação da unidade observacional deve ser alocada em apenas um categoria de cada uma das duas variáveis nominais.

Siegel & Castellan Jr, 1988, p. 111

tabela <- xtabs(~TipoSang+Sexo, data=Dados)
tabela

        Sexo
TipoSang  F  M
     A-  17 32
     A+  74 97
     AB-  2  0
     AB+ 12 28
     B-  14 20
     B+  28 37
     O-  14 31
     O+  67 65

gplots::balloonplot(t(tabela), main ="", xlab ="", ylab="",
                    label=TRUE, show.margins=TRUE,
                    dotcolor="gray")

res <- chisq.test(tabela,
                  simulate.p.value=TRUE, 
                  B=1e6)
res

    Pearson's Chi-squared test with simulated p-value (based on 1e+06
    replicates)

data:  tabela
X-squared = 12.638, df = NA, p-value = 0.07442

df <- chisq.test(tabela)$parameter
cat("Graus de liberdade (df) = ", df, "\n")

Graus de liberdade (df) =  7 

x2 <- res$statistic
n <- sum(res$observed)
r <- nrow(tabela) 
c <- ncol(tabela) 
df <- (r-1)*(c-1) 
cat("Heuristica de significancia: X^2/df > 2\n", round(x2/df,2), "\n")

Heuristica de significancia: X^2/df > 2
 1.81 

STAR <- res$stdres 
cat("MCSTAR\n") # Moment-Corrected STandardized Adjusted Residual

MCSTAR

print(MCSTAR <- STAR/(sqrt((1-1/r)*(1-1/c))), digits=3) 

        Sexo
TipoSang      F      M
     A-  -1.726  1.726
     A+   0.434 -0.434
     AB-  2.498 -2.498
     AB+ -2.490  2.490
     B-  -0.222  0.222
     B+   0.184 -0.184
     O-  -2.416  2.416
     O+   3.390 -3.390

alfa <- 0.05
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
cat("|MCSTAR critico| (alfaBonferroni=5%/",df,") =",zcrit,"\n\n")

|MCSTAR critico| (alfaBonferroni=5%/ 7 ) = 2.69011 

pSTAR <- (1-pnorm(abs(MCSTAR)))*2
pSTAR <- abs(MCSTAR)<abs(zcrit)
pSTAR <- pSTAR*1
if(prod(pSTAR)!=1)
  {
    corrplot::corrplot(MCSTAR, is.corr = FALSE, 
                       insig="label_sig",
                       p.mat=pSTAR,
                       pch="*", 
                       pch.cex=3, 
                       pch.col="yellow")  
  } else
  {
    pSTAR <- pSTAR*0
    corrplot::corrplot(MCSTAR, 
                       is.corr = FALSE, 
                       p.mat=pSTAR)  
  }

V <- DescTools::CramerV(tabela)
cat("V de Cramer = ", round(V,2), "\n")

V de Cramer =  0.15 

effectsize::interpret_cramers_v(V)

[1] "small"
(Rules: funder2019)

Teste qui-quadrado de independência

1. Delineamento entre participantes (observações independentes)

   • \(H_{0}\): independência ou ausência de efeito de interação entre as duas variáveis nominais

2. Suposições do teste qui-quadrado de Pearson

   • Valor-p assintótico: até \(20\%\) das caselas com frequência esperada menor que 5 e todas as caselas com frequência esperada maior que 1
   • Valor-p exato ou por bootstrapping: sem restrições

3. Número de graus de liberdade (df) = \((r-1)(c-1)\)

4. Número de dimensões da tabela de contingência: \(dim = mínimo\{r,c\} - 1\)

5. MCSTAR (Moment-Corrected STandardized Adjusted Residual)

   • Escore-z da casela que indica quanto ela se afasta em unidade de desvio-padrão da hipótese nula de independência
   • Adequado se a frequência absoluta da casela é pelo menos 12 e se o número de dimensões é igual a 1 ou 2; se o número de dimensões é 3 ou mais, usar Análise de Correspondência Simples
   • MCSTAR positivo (negativo): valor da casela acima (abaixo) do valor esperado indica associação positiva (negativa) entre as categorias das duas variáveis nominais
   • A associação entre as categorias é estatisticamente significante se \(|MCSTAR| > z_{crítico}\) com correção de Bonferroni (nível de significância ajustado=\(0.05/df\))

6. \(\phi^{2}\): grau de associação entre as duas variáveis nominais: varia entre \(-dim\) e \(dim\)

7. V de Cramér: medida de tamanho de efeito (correlação implícita absoluta)

   • Valor adimensional que varia entre 0 e 1
   • Não depende do tamanho da amostra e do número de dimensões da tabela de contingência
     
   • Quanto maior, maior o grau de dependência ou do efeito de interação entre as duas variáveis nominais

8. \(\eta^{2}=V^{2}\): medida de tamanho de efeito: porcentagem da variância compartilhada entre as duas variáveis nominais

Teste t para uma condição

alfa <- 0.05
estatura <- c(169, 174, 175, 186)
shapiro.test(estatura)

    Shapiro-Wilk normality test

data:  estatura
W = 0.91181, p-value = 0.492

fit <- t.test(x=estatura, 
       mu=177,
       alternative="two.sided",
       conf.level=1-alfa)
fit

    One Sample t-test

data:  estatura
t = -0.27915, df = 3, p-value = 0.7983
alternative hypothesis: true mean is not equal to 177
95 percent confidence interval:
 164.5993 187.4007
sample estimates:
mean of x 
      176 

d <- lsr::cohensD(estatura, mu=177)
cat("d de Cohen =", round(d,2), "\n")

d de Cohen = 0.14 

effectsize::interpret_cohens_d(d=d, rules="cohen1988")

[1] "very small"
(Rules: cohen1988)

df <- fit$parameter
q1 <- round(qt(p=alfa/2, df=df),2)
q2 <- round(qt(p=1-alfa/2, df=df),2)
t <- round(abs(fit$statistic),2)
DescTools::PlotProbDist(breaks=c(-4.5,q1,-t,t,q2,4.5), 
                        function(x) dt(x, df=df), 
                        blab=c(q1,-t,t,q2), 
                        alab=c("RejH0","p/2","","p/2","RejH0"),
                        xlim=c(-4.5,4.5),
                        main=paste0("t(",df,")"),
                        col=c("darkred","black",
                              "black","black","darkred"), 
                        density=c(20,20,0,20,20))

Número de graus de liberdade

Tamanho do amostra é o tamanho nominal do amostra.

O número de graus de liberdade (gl) constitui o tamanho efetivo da amostra.

Eisenhauer, 2008

O número de graus de liberdade (gl) é a quantidade necessária de unidades experimentais para que o estimador da variância do estimador do efeito do teste (diferença padronizada entre a média amostral e a média hipotetizada), usando todas as medidas independentes da amostra, seja não-viesado.

\[E\left(S^{2}\right)=\sigma^{2}, \; \text{se} \; S^{2}=\dfrac{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}}{n-1} \]

Wonnacott & Wonnacott, 1990, p. 262

Teste t de Welch para duas condições independentes

1. Delineamento entre participantes: independência entre as unidades experimentais

2. Uma VD intervalar

3. Uma VI nominal dicotômica (fator entre participantes)

4. Normalidade da VD nas duas condições independentes, se o tamanho total da amostra é menor que 30; se maior ou igual a 30, a normalidade não é requerida (Teorema Central do Limite)

5. O teste t de Welch robusto à heterocedasticidade é o default da função stats::t.test do R

alfa <- 0.05
print(psych::describeBy(Dados$MCT, group=Dados$Sexo, mat=2))

    item group1 vars   n     mean        sd median trimmed     mad min max
X11    1      F    1 230 57.63043  9.052874     56  56.750  7.4130  41 120
X12    2      M    1 312 71.58974 12.087235     72  70.876 11.1195  43 125
    range      skew  kurtosis        se
X11    79 2.4309659 12.039133 0.5969288
X12    82 0.8963175  2.370362 0.6843049

boxplot(MCT ~ Sexo, 
        data=Dados, 
        ylab="MCT (kg)", 
        xlab="Sexo")

gplots::plotmeans(MCT ~ Sexo, 
                  data=Dados, 
                  p=1-alfa/2,
                  xlab="Instructor", 
                  ylab="MCT (kg)", 
                  main="IC95% Bonferroni da média de MCT",
                  barcol="black",
                  connect=FALSE)

print(RcmdrMisc::normalityTest(MCT ~ Sexo, data=Dados))

 --------
 Sexo = F 

    Shapiro-Wilk normality test

data:  MCT
W = 0.83245, p-value = 4.94e-15

 --------
 Sexo = M 

    Shapiro-Wilk normality test

data:  MCT
W = 0.95412, p-value = 2.614e-08

 --------

 p-values adjusted by the Holm method:
  unadjusted adjusted  
F 4.9401e-15 9.8801e-15
M 2.6139e-08 2.6139e-08
NULL

cat("IC95% das medias brutas estimadas de MCT\n")

IC95% das medias brutas estimadas de MCT

print(Rmisc::group.CI(MCT ~ Sexo, data=Dados, ci=1-alfa/2))

  Sexo MCT.upper MCT.mean MCT.lower
1    F  58.97725 57.63043  56.28362
2    M  73.13101 71.58974  70.04848

fit <- t.test(MCT ~ Sexo, 
              data=Dados)
fit

    Welch Two Sample t-test

data:  MCT by Sexo
t = -15.372, df = 539.86, p-value < 2.2e-16
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
 -15.74310 -12.17552
sample estimates:
mean in group F mean in group M 
       57.63043        71.58974 

d <- lsr::cohensD(MCT ~ Sexo, 
                  data=Dados,
                  method="unequal")
cat("d de Cohen =", round(d,2), "\n")

d de Cohen = 1.31 

effectsize::interpret_cohens_d(d)

[1] "large"
(Rules: cohen1988)

df <- fit$parameter
t <- fit$statistic # estatistica de teste t
F <- t^2 # estatistica de teste F de Fisher
eta2 <- F/(F+df)
cat("eta^2 = R^2 = ",round(eta2,2),"\n")

eta^2 = R^2 =  0.3 

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

      t 
"large" 
(Rules: cohen1992)

cat("IC95% Bonferroni das medias marginais estimadas de MCT\n")

IC95% Bonferroni das medias marginais estimadas de MCT

Sum <- rcompanion::groupwiseMean(MCT ~ Sexo, 
                                 data=Dados, 
                                 conf=1-alfa/2, 
                                 digits=3, 
                                 traditional=FALSE, 
                                 percentile=TRUE,
                                 na.rm=TRUE)
print(Sum)

  Sexo   n Mean Conf.level Percentile.lower Percentile.upper
1    F 230 57.6      0.975             56.4             59.0
2    M 312 71.6      0.975             70.1             73.2

ggplot2::ggplot(Sum, 
                ggplot2::aes(x = Sexo,y = Mean)) +
  ggplot2::geom_errorbar(ggplot2::aes(ymin = Percentile.lower,
                    ymax = Percentile.upper),width=0.05,size=0.5)+
  ggplot2::geom_point(shape=15,size=4) + 
  ggplot2::theme_bw() +
  ggplot2::theme(axis.title=ggplot2::element_text(face="bold")) +
  ggplot2::ggtitle("IC95 Bonferroni de media marginal de MCT (kg)") 

\(\eta^{2}\) de Cohen

\(\eta^{2}\) de Cohen é:

1. Uma estatística de tamanho de efeito que
   • Não depende do tamanho da amostra
   • Varia entre 0 e 1
   • É adimensional
2. A proporção da variância da VD explicada pela VI
3. Uma estatística igual ao coeficiente de determinação \(R^{2}\)
4. O coeficiente de correlação de Pearson implícita ao quadrado entre a VD e a VI do modelo linear geral (GLM).

Teste t para duas condições dependentes

1. Delineamento intraparticipantes: independência entre as unidades experimentais

2. Uma VD intervalar

3. Uma VI nominal dicotômica (fator intraparticipantes)

4. Normalidade da diferença da VD nas duas condições dependentes, se o tamanho total da amostra é menor que 30; se maior ou igual a 30, a normalidade não é requerida (Teorema Central do Limite)

5. O teste t relacionado é executado por meio da função stats::t.test do R

Programa de Educação Nutricional: ingestão de sódio

No exemplo a seguir, o nutricionista faz com que seus mesmos alunos do Programa de Educacao Nutricional mantenham diários do que comem por uma semana e, em seguida, calculem a ingestão diária de sódio em miligramas.

Como a turma recebeu os mesmos programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma entre duas semanas consecutivas.

rcompanion.org de Salvatore S. Mangiafico

“A quantidade diaria de sodio disponível para consumo nos domicilios brasileiros foi de 4,7 g para ingestao diaria de 2.000 kcal, mantendo-se mais de duas vezes superior ao limite recomendado de ingestao desse nutriente.”

SARNO et al. (2013) Estimativa de consumo de sodio pela populacao brasileira 2008-9. Rev Saude Publica 47(3):571-8.

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

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()

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

Anexando pacote: 'data.table'

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

    dcast, melt

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:plyr':

    arrange, count, desc, failwith, id, mutate, rename, summarise,
    summarize

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

    recode

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

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 

effectsize::interpret_cohens_d(d)

[1] "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 

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

      t 
"small" 
(Rules: cohen1992)

ANOVA unifatorial independente de Welch

1. Delineamento entre participantes: independência entre as unidades experimentais

2. Uma VD intervalar

3. Uma VI nominal politômica (fator entre participantes)

4. Normalidade da VD nas três ou mais condições independentes, se o tamanho total da amostra é menor que 30; se maior ou igual a 30, a normalidade não é requerida (Teorema Central do Limite)

5. Teste omnibus: ANOVA de Welch robusta à heterocedasticidade é teste default da função stats::oneway.tests do R

6. Teste post hoc: testes de Games-Howell para a par robustos à heterocedasticidade; função rstatix::games_howell_test do R

Teste omnibus

• \(H_{0}\): As médias populacionais da VD nas condições são iguais

• \(H_{1}\): Pelo menos duas médias populacionais da VD nas condições são diferentes

Teste post hoc

• Testes de Games-Howell par a par.

ANOVA de Welch: MCT~AtivFisica (Masculino)

# MCT~AtivFisica Masculino
print(psych::describeBy(Dados.M$MCT, 
                  Dados.M$AtivFisica, 
                  mat=2, 
                  digits=2))

    item             group1 vars   n  mean    sd median trimmed   mad min max
X11    1     sempre_inativo    1  35 67.91 10.03     65   67.28  7.41  51  95
X12    2 atualmente_inativo    1  89 70.42 14.03     70   69.47 11.86  43 125
X13    3  baixa_intensidade    1  30 73.63 10.37     74   72.92  4.45  54 102
X14    4  media_intensidade    1  50 70.06  9.24     70   69.53  7.41  48  94
X15    5   alta_intensidade    1 108 73.89 12.22     73   73.26 10.38  50 120
    range skew kurtosis   se
X11    44 0.72    -0.05 1.69
X12    82 1.19     3.19 1.49
X13    48 0.63     1.17 1.89
X14    46 0.35     0.09 1.31
X15    70 0.66     1.23 1.18

boxplot(MCT~AtivFisica, data=Dados.M)

print(RcmdrMisc::normalityTest(MCT~AtivFisica, data=Dados.M))

 --------
 AtivFisica = sempre_inativo 

    Shapiro-Wilk normality test

data:  MCT
W = 0.94878, p-value = 0.1039

 --------
 AtivFisica = atualmente_inativo 

    Shapiro-Wilk normality test

data:  MCT
W = 0.92088, p-value = 4.365e-05

 --------
 AtivFisica = baixa_intensidade 

    Shapiro-Wilk normality test

data:  MCT
W = 0.87753, p-value = 0.002476

 --------
 AtivFisica = media_intensidade 

    Shapiro-Wilk normality test

data:  MCT
W = 0.97292, p-value = 0.3032

 --------
 AtivFisica = alta_intensidade 

    Shapiro-Wilk normality test

data:  MCT
W = 0.96681, p-value = 0.008436

 --------

 p-values adjusted by the Holm method:
                   unadjusted adjusted  
sempre_inativo     0.1039402  0.20788038
atualmente_inativo 4.3652e-05 0.00021826
baixa_intensidade  0.0024758  0.00990304
media_intensidade  0.3031693  0.30316931
alta_intensidade   0.0084365  0.02530941
NULL

s <- Rmisc::summarySE(Dados.M, 
                      measurevar="MCT", 
                      groupvars=c("AtivFisica"),
                      conf.interval = 1-alfa/(length(unique(Dados.M$AtivFisica))),
                      na.rm=TRUE)

ggplot2::ggplot(s, 
                ggplot2::aes(x=AtivFisica, 
                             y=MCT)) + 
  ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-ci, 
                                      ymax=MCT+ci), 
                         width=.3) +
  ggplot2::geom_point(size=3,
                      shape=21) +
  ggplot2::xlab("Nivel de Atividade Fisica") +
  ggplot2::ylab("MCT (kg)") +
  ggplot2::ggtitle("IC95% Bonferroni de MCT Masculino") +
  ggplot2::ylim(62,80) +  
  ggplot2::theme_bw() +
  ggplot2::theme(legend.justification=c(1,0),
                 legend.position=c(1,0))

cat("Teste omnibus")

Teste omnibus

fit.M <- oneway.test(MCT~AtivFisica, 
                     data=Dados.M, 
                     na.action=na.omit)
print(fit.M)

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

data:  MCT and AtivFisica
F = 2.8463, num df = 4.00, denom df = 111.12, p-value = 0.0273

cat("Teste post hoc")

Teste post hoc

ph.M <- rstatix::games_howell_test(MCT~AtivFisica, 
                                   data=Dados.M)
print(as.data.frame(ph.M), digits=4)

   .y.             group1             group2 estimate conf.low conf.high p.adj
1  MCT     sempre_inativo atualmente_inativo   2.5014   -3.781     8.784 0.801
2  MCT     sempre_inativo  baixa_intensidade   5.7190   -1.425    12.863 0.176
3  MCT     sempre_inativo  media_intensidade   2.1457   -3.848     8.140 0.853
4  MCT     sempre_inativo   alta_intensidade   5.9746    0.197    11.752 0.039
5  MCT atualmente_inativo  baixa_intensidade   3.2176   -3.532     9.967 0.670
6  MCT atualmente_inativo  media_intensidade  -0.3557   -5.831     5.119 1.000
7  MCT atualmente_inativo   alta_intensidade   3.4732   -1.753     8.699 0.359
8  MCT  baixa_intensidade  media_intensidade  -3.5733  -10.059     2.912 0.533
9  MCT  baixa_intensidade   alta_intensidade   0.2556   -6.036     6.547 1.000
10 MCT  media_intensidade   alta_intensidade   3.8289   -1.038     8.695 0.195
   p.adj.signif
1            ns
2            ns
3            ns
4             *
5            ns
6            ns
7            ns
8            ns
9            ns
10           ns

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

Analise de significancia pratica: tamanho de efeito

F <- as.numeric(fit.M$statistic)
dfn <- as.numeric(fit.M$parameter[1])
dfd <- as.numeric(fit.M$parameter[2])
eta2 <- dfn*F/(dfn*F+dfd)
cat("eta^2 = R^2 = ", round(eta2,4), "\n\n")

eta^2 = R^2 =  0.0929 

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

[1] "small"
(Rules: cohen1992)

ANOVA unifatorial independente de Welch sem dados brutos

nA <- 35
nB <- 89
nC <- 30
nD <- 50
nE <- 108
mediaA <- 67.91
mediaB <- 70.42
mediaC <- 73.63
mediaD <- 70.06
mediaE <- 73.89
dpA <- 10.03
dpB <- 14.04
dpC <- 10.37
dpD <- 9.24
dpE <- 12.22
n <- c(nA,nB,nC,nD, nE)
media <- c(mediaA,mediaB,mediaC,mediaD, mediaE)
dp <- c(dpA,dpB,dpC,dpD, dpE)
J <- length(n)
balanc <- max(dp)/min(dp)
var <- dp^2
w <- n/var
U <- sum(w)
X_til <- as.numeric(w%*%media/U)
gl_num <- J - 1
A <- as.numeric(w%*%((media - X_til)^2)/(J-1))
B <- as.numeric((2*(J - 2)/(J^2 - 1))*(((1 - w/U)/(n-1))%*%(1 - w/U)))
gl_denom <- 1/(((3/2)/(J-2))*B)
F <- A/(1+B)
p <- pf(F,gl_num,gl_denom,lower.tail=FALSE)
razao <- max(dp)/min(dp)
eta2 <- as.numeric(gl_num*F/(gl_num*F + gl_denom))
razao <- max(dp)/min(dp) 
es <- effectsize::interpret_eta_squared(eta2, rules="cohen1992")
cat("\nmax(dp)/min(dp) = ", razao,"")

max(dp)/min(dp) =  1.519481 

cat("\nAnalise de significancia estatistica (ANOVA unifatorial de Welch): valor-p")

Analise de significancia estatistica (ANOVA unifatorial de Welch): valor-p

cat("\n F(",gl_num,",",gl_denom,") = ",F,", p = ",p,"",sep="")

    F(4,111.1203) = 2.847271, p = 0.02726022

cat("\nAnalise de significancia pratica: tamanho de efeito")

Analise de significancia pratica: tamanho de efeito

cat("\n eta^2 = R^2 = ",eta2," (",es[1],")",sep="")

    eta^2 = R^2 = 0.09296497 (small)

Referências

• AGRESTI, A & FINLAY, B (2012) Métodos estatísticos para as Ciências Sociais. Porto Alegre: PENSO.
• ALTHOUSE, AD et al. (2021) Recommendations for statistical reporting in cardiovascular medicine: A special report from the American Heart Association. Circulation 144(4): e70-e91.
• BACCHETTI, P (2010) Current sample size conventions: Flaws, harms, and alternatives. BMC Med 8: 17. https://doi.org/10.1186/1741-7015-8-17.
• COELHO, JP et al. (2008) Inferência estatística: com utilização do SPSS e G*Power. Lisboa: Sílabo.
• DANCEY, C & REIDY, J (2019) Estatística sem Matemática para Psicologia. 7ª ed. Porto Alegre: Penso.
• DI LEO, G & SARDANELLI, F (2020) Statistical significance: p value, 0.05 threshold, and applications to radiomics-reasons for a conservative approach. Eur Radiol Exp 4: 18. https://doi.org/10.1186/s41747-020-0145-y.
• EISENHAUER, JG (2008) Degrees of freedom. Teaching Statistics 30(3): 75-8.
• García-pérez, MA & Núñez-antón, V (2003) Cellwise Residual Analysis in Two-Way Contingency Tables. Educational and Psychological Measurement, 63(5), 825–839. https://doi.org/10.1177/0013164403251280.
• GATÁS, RR (1978) Elementos de probabilidade e inferência. SP: Atlas.
• GERARD, PD & SMITH, DR & WEERRAKKODY, G (1998) Limits of retrospective power analysis. The Journal of Wildlife Management 62(2): 801–7. https://doi.org/10.2307/3802357.
• LOMBARDI, CM & HURLBERT, SH (2009) Misprescription and misuse of one-tailed tests. Austral Ecology 34: 447-68.
• MARK, DB et al. (2016) Understanding the Role of P Values and Hypothesis Tests in Clinical Research. JAMA Cardiol. 1(9): 1048–54. doi:10.1001/jamacardio.2016.3312.
• PERUGINI, M et al. (2018) A practical primer to power analysis for simple experimental designs. International Review of Social Psychology 31(1): 1–23, DOI: https://doi.org/10.5334/irsp.181. Scripts R e Rmd: https://github.com/mcfanda/primerPowerIRSP.
• SIEGEL, S & CASTELLAN Jr, NJ (1988) Nonparametric statistics for the behavioral sciences. 2ª ed. NY: McGraw-Hill.
• SPANOS, A (2014) Recurring controversies about P values and confidence intervals revisited. Ecology 95(3):645-51. doi: 10.1890/13-1291.1. PMID: 24804448.
• STEPHENS, LJ (2009) Statistics in Psychology. NY: McGraw-Hill, Schaum’s Outline Series.
• WONNACOTT, T & WONNACOTT, R (1981) Estatística aplicada à Economia e à Administração. RJ: LTC.
• WONNACOTT, T & WONNACOTT, R (1990) Introductory Statistics for Business and Economics. 4ª ed. NJ: Wiley.

Script R completo

  cat(readLines("Modulo3.R"), sep="\n")

options(warn=-1)   
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(MBESS, warn.conflicts=FALSE))
suppressMessages(library(RVAideMemoire, warn.conflicts=FALSE))
suppressMessages(library(corrplot, warn.conflicts=FALSE))
suppressMessages(library(gplots, warn.conflicts=FALSE))
suppressMessages(library(RcmdrMisc, warn.conflicts=FALSE))
suppressMessages(library(rcompanion, warn.conflicts=FALSE))
suppressMessages(library(Rmisc, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(lsr, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(rstatix, warn.conflicts=FALSE))
suppressMessages(library(jmv, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8")
Sys.setlocale("LC_ALL", "pt_BR.UTF-8")

Dados <- readxl::read_excel(path="Biometria_FMUSP.xlsx",
                            sheet="dados",
                            na=c("NA","na","Na","nA"))
Dados <- dplyr::mutate_if(Dados, is.character, as.factor)
categAF <- c("sempre_inativo",
            "atualmente_inativo",
            "baixa_intensidade",
            "media_intensidade",
            "alta_intensidade")
Dados$AtivFisica <- factor(Dados$AtivFisica,
                           levels=categAF,
                           labels=categAF)
categSexo <- c("M","F")
Dados$Sexo <- factor(Dados$Sexo,
                     levels=categSexo,
                     labels=categSexo)
Dados$MCT[Dados$MCT==658] <- NA
Dados$Estatura[Dados$Estatura==120] <- NA
# https://en.wikipedia.org/wiki/Body_mass_index
Dados$IMC.categ <- cut(Dados$IMC, 
                       breaks=c(0,18.5,25,Inf),
                       labels=c("sub","normal","sobre"),
                       ordered_result=TRUE)
Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")

alfa <- 0.05
estatura <- c(169, 174, 175, 186)
shapiro.test(estatura)
DescTools::ZTest(x=estatura, 
                 sd_pop=7,
                 mu=177,
                 alternative="two.sided",
                 conf.level=1-alfa,
                 na.rm=TRUE)
effectsize::interpret_cohens_d(d=0.14, rules="cohen1988")

alfa <- 0.05
q1 <- round(qnorm(p=alfa/2, mean=0, sd=1),2)
q2 <- round(qnorm(p=1-alfa/2, mean=0, sd=1),2)
DescTools::PlotProbDist(breaks=c(-3,q1,-0.29,0.29,q2,3), 
                        function(x) dnorm(x, mean=0, sd=1), 
                        blab=c(q1,-0.29,0.29,q2), 
                        alab=c("RejH0","p/2","","p/2","RejH0"),
                        xlim=c(-3,3),
                        main="normal(0,1)",                                              col=c("darkred","black",                                                "black","black","darkred"), 
                        density=c(20,20,0,20,20))

alfa <- 0.05
set.seed(123)
s <- rnorm(999, 175.9, 7)
estatura <- c(s, 176*1000-999*mean(s))
shapiro.test(estatura)
DescTools::ZTest(estatura, 
                 sd_pop=7,
                 mu=177,
                 alternative="two.sided",
                 conf.level=1-alfa,
                 na.rm=TRUE)
effectsize::interpret_cohens_d(d=0.14, rules="cohen1988")

alfa <- 0.05
q1 <- round(qnorm(p=alfa/2, mean=0, sd=1),2)
q2 <- round(qnorm(p=1-alfa/2, mean=0, sd=1),2)
DescTools::PlotProbDist(breaks=c(-5,-4.52,q1,q2,4.52,5), 
                        function(x) dnorm(x, mean=0, sd=1), 
                        blab=c(-4.52,q1,q2,4.52), 
                        alab=c("","RejH0","","RejH0",""),
                        xlim=c(-5,5),
                        main="normal(0,1)",                                              col=c("darkred","darkred",                                             "black","darkred","darkred"),                              density=c(20,20,0,20,20))

# Análise de poder retrospectivo (_plug in_) de ANOVA unifatorial independente balanceada (Gerard et al., 1998)

k <- 3 # numero de condicoes independentes
n <- 20 # numero de observacoes independentes de cada condicao
alfa <- 0.05
dfn <- k - 1
dfd <- k*(n - 1)
Fcrt <- qf(1-alfa, dfn, dfd)
Fobs <- Fcrt*0.99 
eta2 <- dfn*Fobs/(dfn*Fobs+dfd)
eta2lims <- MBESS::ci.pvaf(Fobs, dfn, dfd, k*n, 1-alfa)
f2 <- eta2/(1-eta2)
f2.ll <- eta2lims$Lower.Limit.Proportion.of.Variance.Accounted.for/
  (1-eta2lims$Lower.Limit.Proportion.of.Variance.Accounted.for)
f2.ul <- eta2lims$Upper.Limit.Proportion.of.Variance.Accounted.for/
  (1-eta2lims$Upper.Limit.Proportion.of.Variance.Accounted.for)
ncp.p <- dfd*f2 # ou dfn*Fobs
ncp.p.ll <- dfd*f2.ll 
ncp.p.ul <- dfd*f2.ul 
cat(paste("n =", k*n, "\tk =", k, "\tFcrt =", round(Fcrt,2),"\tFobs =", round(Fobs,2),"\n"))
poder.p <- 1-pf(Fcrt,dfn, dfd, ncp.p)
cat(paste("Estimativa pontual do poder retrospectivo =", round(poder.p,3),"\n"))
poder.p.ll <- 1-pf(Fcrt,dfn, dfd, ncp.p.ll)
cat(paste("Limite inferior do IC95% do poder retrospectivo =", round(poder.p.ll,3),"\n"))
poder.p.ul <- 1-pf(Fcrt,dfn, dfd, ncp.p.ul)
cat(paste("Limite superior do IC95% do poder retrospectivo =", round(poder.p.ul,3),"\n"))

tipsang.freq <- as.numeric(table(Dados$TipoSang))
tipsang.prob <- c(0.08, 0.34, 0.005, 0.025,
                  0.02, 0.08, 0.09, 0.36)
tslab <- c("A-", "A+", "AB-", "AB+", 
           "B-", "B+", "O-", "O+")
plot(tipsang.prob, type="p", lty=1, 
     lwd=1, pch=2, col="black", 
     xaxt = "n", xlab='Tipo Sanguineo', ylab="Probabilidade")
axis(1, at=1:length(tslab), labels=tslab)
lines(tipsang.freq/sum(tipsang.freq), pch=21, 
      type="p", lty=1, lwd=1, col="black")
legend("top", legend=c("Prob.Hipotet","Prob.Observ"), 
       pch=c(2,21), bty="n")
cat(tslab,"\n")
print(tipsang.freq)
print(round(tipsang.freq/sum(tipsang.freq), 3))
print(tipsang.prob)
out <- chisq.test(x=tipsang.freq, 
                  p=tipsang.prob,
                  simulate.p.value=TRUE, 
                  B=1e4) 
df <- chisq.test(tipsang.freq, p=tipsang.prob)$parameter
X2 <- out$statistic
p <- out$p.value
n <- sum(chisq.test(tipsang.freq, p=tipsang.prob)$observed)
cat("X^2 = ", X2 , " df = ",  df, " p = ", p, "\n", sep="")
cat("Heuristica de significancia: X^2/df > 2\n", round(X2/df,2), "\n")
cat("Residuos estandardizados\n", round(out$stdres,2), "\n")
V <- sqrt((X2/n)/df)
cat("V de Cramer = ", round(V,2), "\n")
cat("Grau de não aderência:\n")
print(effectsize::interpret_cramers_v(V))
RVAideMemoire::multinomial.multcomp(tipsang.freq, p.method="holm")

tabela <- xtabs(~TipoSang+Sexo, data=Dados)
tabela
gplots::balloonplot(t(tabela), main ="", xlab ="", ylab="",
                    label=TRUE, show.margins=TRUE,
                    dotcolor="gray")
res <- chisq.test(tabela,
                  simulate.p.value=TRUE, 
                  B=1e4)
res
df <- chisq.test(tabela)$parameter
cat("Graus de liberdade (df) = ", df, "\n")
x2 <- res$statistic
n <- sum(res$observed)
r <- nrow(tabela) 
c <- ncol(tabela) 
df <- (r-1)*(c-1) 
cat("Heuristica de significancia: X^2/df > 2\n", round(x2/df,2), "\n")
STAR <- res$stdres 
cat("MCSTAR\n") # Moment-Corrected STandardized Adjusted Residual
print(MCSTAR <- STAR/(sqrt((1-1/r)*(1-1/c)))) 
alfa <- 0.05
alfaBonf <- alfa/df
zcrit <- abs(qnorm(alfaBonf/2))
cat("|MCSTAR critico| (alfaBonferroni=5%/",df,") =",zcrit,"\n\n")
pSTAR <- (1-pnorm(abs(MCSTAR)))*2
pSTAR <- abs(MCSTAR)<abs(zcrit)
pSTAR <- pSTAR*1
if(prod(pSTAR)!=1)
{
  corrplot::corrplot(MCSTAR, is.corr = FALSE, 
                     insig="label_sig",
                     p.mat=pSTAR,
                     pch="*", 
                     pch.cex=3, 
                     pch.col="brown")  
} else
{
  pSTAR <- pSTAR*0
  corrplot::corrplot(MCSTAR, 
                     is.corr = FALSE, 
                     p.mat=pSTAR)  
}
V <- DescTools::CramerV(tabela)
cat("V de Cramer = ", round(V,2), "\n")
effectsize::interpret_cramers_v(V)

alfa <- 0.05
estatura <- c(169, 174, 175, 186)
shapiro.test(estatura)
fit <- t.test(x=estatura, 
       mu=177,
       alternative="two.sided",
       conf.level=1-alfa)
fit
d <- lsr::cohensD(estatura, mu=177)
cat("d de Cohen =", round(d,2), "\n")
effectsize::interpret_cohens_d(d=d, rules="cohen1988")
df <- fit$parameter
q1 <- round(qt(p=alfa/2, df=df),2)
q2 <- round(qt(p=1-alfa/2, df=df),2)
t <- round(abs(fit$statistic),2)
DescTools::PlotProbDist(breaks=c(-4.5,q1,-t,t,q2,4.5), 
                        function(x) dt(x, df=df), 
                        blab=c(q1,-t,t,q2), 
                        alab=c("RejH0","p/2","","p/2","RejH0"),
                        xlim=c(-4.5,4.5),
                        main=paste0("t(",df,")"),
                        col=c("darkred","black",
                              "black","black","darkred"), 
                        density=c(20,20,0,20,20))

alfa <- 0.05
psych::describeBy(Dados$MCT, group=Dados$Sexo, mat=2)
boxplot(MCT ~ Sexo, 
        data=Dados, 
        ylab="MCT (kg)", 
        xlab="Sexo")
car::densityPlot(MCT ~ Sexo, 
                 data=Dados, 
                 col=c("black","black"),
                 rug=FALSE)
gplots::plotmeans(MCT ~ Sexo, 
                  data=Dados, 
                  p=1-alfa/2,
                  xlab="Instructor", 
                  ylab="MCT (mg/semana)", 
                  main="IC95% Bonferroni da média de MCT",
                  barcol="black",
                  connect=FALSE)
RcmdrMisc::normalityTest(MCT ~ Sexo, data=Dados)
cat("IC95% das medias brutas estimadas de MCT\n")
Rmisc::group.CI(MCT ~ Sexo, data=Dados, ci=1-alfa/2)
fit <- t.test(MCT ~ Sexo, 
              data=Dados)
fit
d <- lsr::cohensD(MCT ~ Sexo, 
                  data=Dados,
                  method="unequal")
cat("d de Cohen =", round(d,2), "\n")
effectsize::interpret_cohens_d(d)
df <- fit$parameter
t <- fit$statistic # estatistica de teste t
F <- t^2 # estatistica de teste F de Fisher
eta2 <- F/(F+df)
cat("eta^2 = R^2 = ",round(eta2,2),"\n")
effectsize::interpret_eta_squared(eta2, rules="cohen1992")
cat("IC95% Bonferroni das medias marginais estimadas de MCT\n")
Sum <- rcompanion::groupwiseMean(MCT ~ Sexo, 
                                 data=Dados, 
                                 conf=1-alfa/2, 
                                 digits=3, 
                                 traditional=FALSE, 
                                 percentile=TRUE,
                                 na.rm=TRUE)
print(Sum)
ggplot2::ggplot(Sum, 
                ggplot2::aes(x = Sexo,y = Mean)) +
  ggplot2::geom_errorbar(ggplot2::aes(ymin = Percentile.lower,
                    ymax = Percentile.upper),width=0.05,size=0.5)+
  ggplot2::geom_point(shape=15,size=4) + 
  ggplot2::theme_bw() +
  ggplot2::theme(axis.title=ggplot2::element_text(face="bold")) +
  ggplot2::ggtitle("IC95 Bonferroni de media marginal de MCT (kg)") 

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)
Dados_long <- reshape2::melt(Dados,
                             id.vars="Participante",
                             measure.vars=c("Semana1","Semana2"),
                             variable.name="Semana")
names(Dados_long) <- c("Participante", "Semana", "Sodio")
Dados_long
psych::describeBy(Sodio~Semana, data=Dados_long, mat=2, digits=2)
psych::describe(dif)

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()

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

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()

cat("Analise de significancia estatistica: valor p\n")
print(res <- t.test(dif,mu=0))
cat("Analise de significancia pratica: tamanho de efeito\n")
d <- lsr::cohensD(x=Dados$Semana2, 
                  y=Dados$Semana1,
                  method="paired")
cat("d de Cohen =", round(d,2), "\n")
effectsize::interpret_cohens_d(d)
t <- res$statistic
F <- t^2
df <- res$parameter
eta2 <- F/(F+df)
cat("eta^2 = R^2 = ", eta2, "\n\n")
effectsize::interpret_eta_squared(eta2, rules="cohen1992")

# ANOVA de Welch: MCT~AtivFisica (Masculino)

psych::describeBy(Dados.M$MCT, 
                  Dados.M$AtivFisica, 
                  mat=2, 
                  digits=2)
boxplot(MCT~AtivFisica, data=Dados.M)
car::densityPlot(MCT~AtivFisica, data=Dados.M)
RcmdrMisc::normalityTest(MCT~AtivFisica, data=Dados.M)
plot.design(MCT~AtivFisica, data=Dados.M)
s <- Rmisc::summarySE(Dados.M, 
                      measurevar="MCT", 
                      groupvars=c("AtivFisica"),
                      conf.interval = 1-alfa/(length(unique(Dados.M$AtivFisica))),
                      na.rm=
                        TRUE)
pd <- ggplot2::position_dodge(0.5) 
ggplot2::ggplot(s, 
                ggplot2::aes(x=AtivFisica, 
                             y=MCT)) + 
  ggplot2::geom_errorbar(ggplot2::aes(ymin=MCT-ci, 
                                      ymax=MCT+ci), 
                         width=.3, 
                         position=pd) +
  ggplot2::geom_line(position=pd) +
  ggplot2::geom_point(position=pd, 
                      size=3,
                      shape=21) +
  ggplot2::xlab("Nivel de Atividade Fisica") +
  ggplot2::ylab("MCT (kg)") +
  ggplot2::ggtitle("IC95% Bonferroni de MCT Masculino") +
  ggplot2::ylim(62,80) +  
  ggplot2::theme_bw() +
  ggplot2::theme(legend.justification=c(1,0),
                 legend.position=c(1,0))

cat("Teste omnibus")
fit.M <- oneway.test(MCT~AtivFisica, 
                     data=Dados.M, 
                     na.action=na.omit)
print(fit.M)
cat("Teste post hoc")
ph.M <- rstatix::games_howell_test(MCT~AtivFisica, 
                                   data=Dados.M)
print(as.data.frame(ph.M))
cat("Analise de significancia pratica: tamanho de efeito\n")
F <- as.numeric(fit.M$statistic)
dfn <- as.numeric(fit.M$parameter[1])
dfd <- as.numeric(fit.M$parameter[2])
eta2 <- dfn*F/(dfn*F+dfd)
effectsize::interpret_eta_squared(eta2, rules="cohen1992")