ANOVA Unifatorial em R
MSP4112: Estatística Aplicada à Medicina
01 outubro 2025 21:54h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(bootES, warn.conflicts=FALSE))
suppressMessages(library(bruceR, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts = FALSE))
suppressMessages(library(data.table, warn.conflicts = FALSE))
suppressMessages(library(dplyr, warn.conflicts = FALSE))
suppressMessages(library(effectsize, warn.conflicts = FALSE))
suppressMessages(library(emmeans, warn.conflicts=FALSE))
suppressMessages(library(FSA, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts = FALSE))
suppressMessages(library(gplots, warn.conflicts = FALSE))
suppressMessages(library(HH, warn.conflicts = FALSE))
suppressMessages(library(jmv, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lmboot, warn.conflicts=FALSE))
suppressMessages(library(lmerTest, warn.conflicts=FALSE))
suppressMessages(library(MuMIn, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(onewaytests, warn.conflicts=FALSE))
suppressMessages(library(pbkrtest, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(pwr, warn.conflicts=FALSE))
suppressMessages(library(RcmdrMisc, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(Rmisc, warn.conflicts=FALSE))
suppressMessages(library(rstatix, warn.conflicts=FALSE))
suppressMessages(library(sandwich, warn.conflicts=FALSE))
suppressMessages(library(WRS2, warn.conflicts=FALSE))source("eiras.OneWayANOVA.Fisher.WORD.R")
source("eiras.OneWayANOVA.Welch.WORD.R")
source("summarySEwithin2.R")Material
- HTML de Rmarkdown em
RPubs
Objetivos
- reconhecer e mencionar propriedades da distribuição F.
- reconhecer as indicações e aplicar ANOVA para três ou mais condições independentes.
- reconhecer as indicações e aplicar ANOVA para três ou mais condições dependentes (dados longitudinais).
- definir hipóteses estatísticas nula e alternativa.
- executar e interpretar os testes estatísticos omnibus e post hoc.
- conceituar estatísticas de tamanho de efeito.
ANOVA unifatorial em R
- Planejamento de ANOVA unifatorial
- Independente de Fisher
pwr::pwr.anova.testWebPower::wp.anova
- Relacionada
WebPower::wp.rmanova
- ANOVA independente unifatorial
- Fisher (homocedástico)
- teste omnibus:
lm,car::Anova,effectsize::eta_squared- testes post hoc:
emmeans::emmeans,emmeans::contrast,multcomp::cld
- testes post hoc:
- teste omnibus:
oneway.test(var.equal=TRUE, ...),effectsize::eta_squared- testes post hoc:
pairwise.t.test(pool.sd=FALSE, ...)
- testes post hoc:
- teste omnibus:
jmv::ANOVA- testes post hoc:
jmv::ANOVA
- testes post hoc:
- teste omnibus:
onewaytests::onewaytests(method="aov", ...)- testes post hoc:
onewaytests::paircomp
- testes post hoc:
- Fisher sem dados brutos (homocedástico)
- teste omnibus:
demo_ANOVA_Independente_Unifatorial_Fisher_SemDadosBrutos.R,HH::anovaMean- testes post hoc: Post-Hoc Analysis with Tukey’s Test
- Fisher-White (heterocedástico)
- teste omnibus:
lm,car::Anova(white.adjust="hc2", ...)),effectsize::eta_squared- testes post hoc:
sandwich::vcovHC,emmeans::emmeans,emmeans::contrast,multcomp::cld
- testes post hoc:
- Welch (heterocedástico)
- teste omnibus:
oneway.test(var.equal=FALSE, ...)- testes post hoc:
rstatix::games_howell_test
- testes post hoc:
- teste omnibus:
jmv::anovaOneW- testes post hoc:
jmv::anovaOneW(phMethod="gamesHowell", ...)
- testes post hoc:
- teste omnibus:
rstatix::welch_anova_test- testes post hoc:
rstatix::games_howell_test
- testes post hoc:
- teste omnibus:
onewaytests::onewaytests(method="welch", ...)- testes post hoc:
onewaytests::paircomp
- testes post hoc:
- Welch sem dados brutos
- teste omnibus:
eiras.OneWayANOVA.Welch.WORD.R- testes post hoc:
eiras.OneWayANOVA.Welch.WORD.R
- testes post hoc:
- Reamostragem (bootstrapping)
- teste omnibus:
WRS2::t1waybt- testes post hoc:
WRS2::mcppb20
- testes post hoc:
- teste omnibus:
onewaytests::onewaytests(method="gtb",...)- testes post hoc:
onewaytests::paircomp
- testes post hoc:
- teste omnibus:
lmboot::ANOVA.boot- testes post hoc: ?
- Permutação (
Permutation F test)
- teste omnibus:
onewaytests::onewaytests(method="pf", ...)- testes post hoc:
onewaytests::paircomp
- testes post hoc:
- ANOVA relacionada unifatorial ou para medidas repetidas (rm)
- rmANOVA
- teste omnibus:
lmerTest::lmer,car::Anova(test.statistic="F", ...),effectsize::eta_squared- testes post hoc:
emmeans::emmeans,emmeans::contrast,multcomp::cld
- testes post hoc:
- Reamostragem (bootstrapping)
- teste omnibus:
WRS2::rmanovab- testes post hoc:
WRS2::pairdepb
- testes post hoc:
Introdução
As variáveis envolvidas no teste estatístico são:
- VD (variável dependente) intervalar com distribuição normal por condição: univariada
- VI (variável independente) nominal com três ou mais níveis (nominal politômica, fator): unifatorial, one way
ANOVA unifatorial pode ser:
- independente (fator entre participantes)
- relacionada (fator intraparticipantes)
ANOVA independente unifatorial (One-way Analysis of Variance) é uma extensão do teste t independente.
As suposições de normalidade homocedasticidade são condições suficientes e a suposição de independência é necessária para ANOVA independente unifatorial.
Exemplo de condições suficientes: lei do equilíbrio genético de Hardy-Weinberg
Ilustramos com uma analogia, recordando o a lei do equilíbrio genético de Hardy-Weinberg.
Em 1908, G. Hardy e W. Weinberg propuseram independentemente que a frequência de alelos e genótipos em uma população permanecerá constante de geração para geração se a população for estável e em equilíbrio genético.
Cinco condições são necessárias para que uma população permaneça em equilíbrio de Hardy-Weinberg:
- Uma grande população reprodutiva
- Acasalamento aleatório
- Nenhuma mudança na frequência alélica devido a mutação
- Nenhuma imigração ou emigração
- Nenhuma seleção natural
http://www.phschool.com/science/biology_place/labbench/lab8/concepts.html
Violando as condições de H-W, em simulação (implementada em demo_HardyWeinberg.R) com:
- 500 indivíduos,
- mutação alterando a frequência dos alelos,
- selecionando para o alelo A.
Observamos, por exemplo:
As condições matemáticas de população muito grande (ou infinita), sem mutação e sem seleção garantem a população no equilíbrio de Hardy-Weinberg; porém, dentro de certos limites, ainda encontramos predição aceitável para o Equilíbrio de Hardy-Weinberg e poderíamos utilizá-los em populações biológicas.
O que é independência?
Um delineamento entre participantes é aquele em que cada unidade experimental é submetida a somente uma condição experimental.
Por exemplo, três grupos com cinco participantes cada um (total de 15 participantes, nomeados de A a O) são alocados assim:
| Condição 1 | Condição 2 | Condição 3 | |
| A | F | K | |
| B | G | L | |
| C | H | M | |
| D | I | N | |
| E | J | O |
Num delineamento intraparticipantes, o mesmo participante é submetido a todas as condições experimentais. Os mesmos 15 participantes são alocados assim:
| Condição 1 | Condição 2 | Condição 3 | |
| A | A | A | |
| B | B | B | |
| C | C | C | |
| D | D | D | |
| E | E | E | |
| F | F | F | |
| G | G | G | |
| H | H | H | |
| I | I | I | |
| J | J | J | |
| K | K | K | |
| L | L | L | |
| M | M | M | |
| N | N | N | |
| O | O | O |
O que é balanceamento?
Em delineamento em que as condições independentes são balanceadas, o número de participantes em cada condição independente é aproximadamente igual. Por exemplo, com 15 participantes, poderíamos ter algo desbalanceado:
| Condição 1 | Condição 2 | Condição 3 | |
| A | H | K | |
| B | I | L | |
| C | J | M | |
| D | N | ||
| E | O | ||
| F | |||
| G |
Em um delineamento intraparticipantes com medidas repetidas no formato de arquivo largo (wide), o desbalanceamento ocorre quando pelo menos uma observação de um participante está ausente numa condição experimental dependente, por exemplo:
| Condição 1 | Condição 2 | Condição 3 | |
| A | A | A | |
| B | B | ||
| C | C | C | |
| D | D | D | |
| E | E | E | |
| F | F | F | |
| G | G | G | |
| H | H | H | |
| I | I | I | |
| J | J | ||
| K | K | K | |
| L | L | L | |
| M | M | M | |
| N | N | N | |
| O | O | O |
Se os dados de delineamento intraparticipantes estão no formato de arquivo longo (long), i.e., dados longitudinais, então o desbalanceamento não acarreta valor faltante.
Em arquivo no formato wide, basta uma observação faltante (missing value) para caracterizar o desbalanceamento. As estatísticas e modelos de efeitos fixos do R funcionam com a eliminação da unidade experimental com valor faltante. Para delineamento com dados longitudinais (arquivo long), a sugestão é usar modelo com efeito aleatório para usar todos os dados observados.
O que é variância?
Variância é o quadrado do desvio-padrão.
Variância é também a covariância da variável com ela mesma.
Variância, conceitualmente, é a quantidade de informação de variável intervalar.
Se há duas condições independentes, a terceira interpretação da variância total da variável é que ela é uma composição da variâncias entre e intra condições.
- Higgins JPT, Li T, Deeks JJ (editors). Chapter 6: Choosing effect
measures and computing estimates of effect. In: Higgins JPT, Thomas J,
Chandler J, Cumpston M, Li T, Page MJ, Welch VA (editors). Cochrane
Handbook for Systematic Reviews of Interventions version 6.3 (updated
February 2022). Cochrane, 2022. Available from
www.training.cochrane.org/handbook.
- Table 6.5.a em https://training.cochrane.org/handbook/current/chapter-06
Suponha que não temos acesso aos dados brutos de massa corporal total (MCT) dos estudantes do curso de Medicina da FMUSP. Os dados disponíveis são apresentados na tabela abaixo:
| group | n | mean | sd |
|---|---|---|---|
| F | 230 | 57.63 | 9.05 |
| M | 312 | 71.59 | 12.09 |
Quais são os valores da média e desvio-padrão amostrais dos 542 estudantes?
\[ \bar{x} = \dfrac{n_F \, \bar{x}_F + n_M \, \bar{x}_M}{n_F + n_M} \]
\[\begin{align} s^2&=\dfrac{n_F\,s_F^2+n_M\,s_M^2+n_F\,\left(\bar{x}_F-\bar{x}\right)^2+n_M\,\left(\bar{x}_M-\bar{x}\right)^2}{n_F+n_M} \end{align} \]
| n | mean | sd |
|---|---|---|
| 542 | 65.67 | 12.9 |
Media = 65.67dp <- sqrt((230*9.05^2 + 312*12.09^2 +
230*(57.63-65.67)^2 + 312*(71.59-65.67)^2)/542)
cat("Desvio-padrao = ", round(dp,2), "\n", sep="")Desvio-padrao = 12.9dp <- sqrt(((230-1)*9.05^2 + (312-1)*12.09^2 +
(230*312/542)*(57.63-71.59)^2)/(542-1))
cat("Desvio-padrao = ", round(dp,2), "\n", sep="")Desvio-padrao = 12.9A variância total pode ser decomposta em:
\[ s^2 \;=\; \dfrac{n_F\,s_F^2+n_M\,s_M^2}{n_F+n_M} \;+\; \dfrac{n_F\,\left(\bar{x}_F-\bar{x}\right)^2+n_M\,\left(\bar{x}_M-\bar{x}\right)^2}{n_F+n_M} \]
Primeiro termo: variância-dentro (within):
\[ s_{\text{dentro}}^2 = \dfrac{n_F\,s_F^2+n_M\,s_M^2}{n_F+n_M} \]
Segundo termo: variância-entre (between):
\[ s_{\text{entre}}^2 = \dfrac{n_F\,\left(\bar{x}_F-\bar{x}\right)^2+n_M\,\left(\bar{x}_M-\bar{x}\right)^2}{n_F+n_M}. \]
Assim,
\[ s^2 = s_{\text{dentro}}^2 + s_{\text{entre}}^2 \]
Teste F
O teste F na ANOVA independente unifatorial é usado para testar a hipótese nula (\(H_0\)) de que as médias populacionais de três ou mais condições independentes são iguais. A hipótese alternativa (\(H_1\)) afirma que pelo menos duas médias populacionais diferem entre si.
A distribuição da estatística de teste F tem distribuição F de Fisher-Snedecor.
Para a comparação de três ou mais médias populacionais das condições, com base no teste F, os graus de liberdade do numerador dependem do número de condições independentes e os do denominador dependem do tamanho da amostra do número de condições independentes.
A estatística de teste do efeito do fator e também omnibus da ANOVA independenten unifatorial é definida por
\[ F = \frac{\text{QM}_{\text{entre}}}{\text{QM}_{\text{dentro}}} = \frac{\dfrac{\text{SQ}_{\text{entre}}}{g-1}}{\dfrac{\text{SQ}_{\text{dentro}}}{n-g}} \]
Sendo que:
- \(\text{SQ}_{\text{entre}} = \sum_{i=1}^g
n_i(\bar X_{i} - \bar X)^2\)
- \(\text{SQ}_{\text{dentro}} = \sum_{i=1}^g \sum_{j=1}^{n_i}(X_{ij} - \bar X_{i})^2\)
com \(g\) é o número de condições independentes, \(n = \sum_{i=1}^g n_i\) é número total de unidades experimentais, \(\text{SQ}\) é soma de quadrados e \(\text{QM}\) é média da soma de quadrados (“quadrado médio”).
Sob a hipótese nula \(H_0: \mu_1 = \mu_2 = \cdots = \mu_g\), a estatística de teste F, segue a distribuição F de Fisher-Snedecor com os números de graus de liberdade do efeito do fator A, \(df_A=g-1\), e do efeito residual, \(df_E=\sum_{i=1}^{g}{(n_i-1)}=(n-1)-(g-1)=n-g\), no numerador e denominador, respectivamente, sendo que \(g\) é o número de condições independentes e \(n = \sum_{i=1}^g n_i\) é número total de unidades experimentais. \(df_A\) quantifica a complexidade do modelo (GLM univariado) e \(df_E\) quantifica o tamanho de amostra efetivo, sendo que \(n\) é o tamanho de amostra nominal.
\[
\begin{align}
F &\underset{H_0}{\sim} F_{df_A, \, df_E}\\
F &\underset{H_0}{\sim} F_{g-1, \, n-g}\\
\end{align}
\] Familiarize-se com a distribuição F, observando demo_AnimacaoF.R:
- F é um valor positivo.
- sob \(H_0\), a distribuição
F da estatística de teste F tem um parâmetro de não
centralidade (
ncp) nulo; sob \(H_1\) o parâmetro de não centralidade é positivo e é função do tamanho de efeito e do tamanho da amostra. - a distribuição F é assimétrica positiva.
- consideramos somente a cauda superior para localizar \(\alpha\).
- há dois valores para graus de liberdade: para o numerador (número de grupos - 1) e denominador (número de unidades experimentais - número de grupos).
- observe o que acontece com a distribuição sob \(H_0\) e com o valor de \(F\) crítico de 95% (\(F_c\)) à medida que os graus de liberdade aumentam.
Distribuição F em R:
O delineamento do estudo, o tipo de variável e, consequentemente, a estatística adequada mudam, mas o problema é sempre o mesmo: há incerteza para a tomada de decisão sobre a hopótese nula porque lidamos com uma evidência amostral.
A diferença, aqui, é que teremos que lidar com três ou mais condições independentes.
A hipótese nula é a igualdade das médias populacionais dos \(m\) grupos. Caso não rejeitemos \(H_0\):
Quando rejeitamos \(H_0\), concluímos que as médias populacionais não são todas iguais:
que pelo duas médias populacionais são diferentes. Por exemplo:
Teste t independente
Brendon Small e Coach McGuirk fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.
Desde que as classes receberam diferentes programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma para as duas turmas.
Aplicamos o teste t de Welch, robusto à heterocedasticidade
e default da função t.test:
Brendon Small Coach McGuirk
20 20
Instructor Student Sodium
1 Brendon Small a 1200
2 Brendon Small b 1400
3 Brendon Small c 1350
4 Brendon Small d 950
5 Brendon Small e 1400
6 Brendon Small f 1150
7 Brendon Small g 1300
8 Brendon Small h 1325
9 Brendon Small i 1425
10 Brendon Small j 1500
11 Brendon Small k 1250
12 Brendon Small l 1150
13 Brendon Small m 950
14 Brendon Small n 1150
15 Brendon Small o 1600
16 Brendon Small p 1300
17 Brendon Small q 1050
18 Brendon Small r 1300
19 Brendon Small s 1700
20 Brendon Small t 1300
21 Coach McGuirk u 1100
22 Coach McGuirk v 1200
23 Coach McGuirk w 1250
24 Coach McGuirk x 1050
25 Coach McGuirk y 1200
26 Coach McGuirk z 1250
27 Coach McGuirk aa 1350
28 Coach McGuirk ab 1350
29 Coach McGuirk ac 1325
30 Coach McGuirk ad 1525
31 Coach McGuirk ae 1225
32 Coach McGuirk af 1125
33 Coach McGuirk ag 1000
34 Coach McGuirk ah 1125
35 Coach McGuirk ai 1400
36 Coach McGuirk aj 1200
37 Coach McGuirk ak 1150
38 Coach McGuirk al 1400
39 Coach McGuirk am 1500
40 Coach McGuirk an 1200
Instructor n mean sd min Q1 median Q3 max
1 Brendon Small 20 1287.5 193.7 950 1150.0 1300.0 1400 1700
2 Coach McGuirk 20 1246.2 142.4 1000 1143.8 1212.5 1350 1525
Welch Two Sample t-test
data: Sodium by Instructor
t = 0.76722, df = 34.893, p-value = 0.4481
alternative hypothesis: true difference in means between group Brendon Small and group Coach McGuirk is not equal to 0
95 percent confidence interval:
-67.91132 150.41132
sample estimates:
mean in group Brendon Small mean in group Coach McGuirk
1287.50 1246.25 Teste t independente para fator com três grupos
Suponha que uma terceira classe junte-se ao experimento (Melissa Robins). Podemos, então, verificar as diferenças destas três condições experimentais independentes com testes t, comparando-se Brendon Small com Coach McGuirk, Brendon Small com Melissa Robins e Coach McGuirk com Melissa Robins?
Quantos testes t precisaremos de acordo com o número de grupos?
É a combinatória do número de grupos (\(m\)) dois a dois:
\[ {m \choose 2} = { {m!} \over {2! (m-2)!} } \]
- para 2 grupos: \({2 \choose 2} = { {2!} \over {2! (2-2)!} } = 1~\text{teste}\)
- para 3 grupos: \({3 \choose 2} = { {3!} \over {2! (3-2)!} } = 3~\text{testes}\)
- para 4 grupos: \({4 \choose 2} = { {4!} \over {2! (4-2)!} } = 6~\text{testes}\)
- para 6 grupos: \({6 \choose 2} = { {6!} \over {2! (6-2)!} } = 15~\text{testes}\)
Graficamente (demo_Testes_t_2a2.R):
Embora o número cresça rapidamente, podemos não nos impressionar, pois temos computadores para o trabalho repetitivo.
Há, porém, um problema mais grave: as probabilidades de erros do tipo I (\(\alpha\)) e do tipo II (\(\beta\)) acumulam. Quando escolhemos \(\alpha\), a probabilidade de rejeitar incorretamente a hipótese nula, o máximo valor que o valor p pode assumir, a probabilidade de erro do tipo I para \(m\) grupos é:
\[ P(\alpha | m) = 1-(1-\alpha)^{m \choose 2} \]
Também gostaríamos de manter a probabilidade de rejeitar corretamente a hipótese nula, o poder do teste, \(1-\beta\); para \(m\) grupos é:
\[ P(\beta | m) = (1-\beta) ^ {m \choose 2} \]
Graficamente (demo_Testes_t_AlfaPoder.R), podemos
observar o que acontece com o número crescente de pares de teste
t necessários, considerando os tradicionais \(\alpha=0.05\) e \(\beta=0.1\):
Portanto, a probabilidade de erro do tipo I cresce rapidamente para quase \(100\%\) e o poder de seus testes combinados converge para zero. Na prática, dificilmente teremos mais do que 6 grupos. Porém, ainda assim, teríamos \(\alpha \approx 53.7\%\) e poder \(1-\beta \approx 3.5\%\), valores totalmente inaceitáveis para uma boa inferência estatística.
Podemos resolver com vários testes t?
A resposta é:
NÃO!
É melhor, portanto, testar tudo simultaneamente paramanter \(\alpha\) no valor pretendido e
preservar \(1-\beta\).
Análise da variância
Esta análise, que tem o acrônimo ANOVA (do inglês, Analysis of Variance) utiliza variâncias, mas…
… as conclusões são sobre as médias populacionais dos \(m\) grupos.
Como?
Princípio do teste F da ANOVA independente unifatorial de Fisher
Consideremos que existam três grupos nos quais a VD tem distribuição normal homocedástica.
De cada um deles, retiramos uma amostra:
Não havendo problemas, as três amostras reproduzem a distribuição, média e variância das populações das quais se originaram.
A estatística de teste F da ANOVA independente unifatorial é dada por
\[ F = \dfrac{S_E^2}{S_D^2} \]
em que \(S_E^2\) é a variância entre os grupos e \(S_D^2\) é a variância dentro dos grupos.
Variância entre (between) os grupos
A variância entre os grupos presume que as médias amostrais (\(\bar{X}_m\)) refletem as respectivas médias populacionais (\(\mu_m\)):
Como sempre, não temos esta certeza e lidamos somente com a informação das amostras, \(\bar{X}_m\):
Cada média é um número e, portanto, podemos calcular a variância destes números:
Sendo assim, para \(F = \dfrac{S_E^2}{S_D^2}\), a estatística F aumenta quando a variância entre os grupos aumentar.
Variância dentro (within) dos grupos
A variância dentro dos grupos é uma medida da variância total, desconsiderando a média de cada uma das condições. Cada amostra tem sua própria distribuição (presumivelmente, reflexo da distribuição da população de onde foi extraída):
Como sempre, não temos esta certeza e lidaremos somente com a informação das amostras, \(\bar{X}_m\):
Mesclando as três distribuições, estimamos a variância dentro dos grupos, \(S_D^2\), uma medida de quanto, como um todo, a variável é dispersa:
Caso a variância em cada condição seja maior,
esta variância será refletida em \(S_D^2\):
Sendo assim, para \(F = \dfrac{S_E^2}{S_D^2}\), a estatística F diminui quando a variância dentro os grupos (ou condições) aumentar.
Comportamento de \(F = \dfrac{S_E^2}{S_D^2}\)
É fácil imaginar o comportamento da estatística F combinando-se o que pode acontecer com \(S_E^2\) e \(S_D^2\):
Desta forma, Ronald Fisher inventou uma forma de comparar médias entre várias condições, simultaneamente, utilizando somente a comparação entre duas variâncias, com o numerador refletindo a dispersão das médias e o denominador refletindo a dispersão do fenômeno em estudo.
Lógica estatística da estatística de teste F
A ANOVA unifatorial compara duas formas de estimar a variância populacional \(\sigma^2\) da variável dependente.
A ideia de base é a seguinte: Vamos estimar a variância \(\sigma^2\) da VD por dois métodos diferentes, um que não depende da veracidade de \(H_0\) e outro que sim.
Depois comparamos as duas estimativas. Se os grupos tiverem todos a mesma média (\(H_0\) verdadeiro) as duas estimativas deverão ser próximas, senão deverão diferir significativamente.
Uma forma de estimar \(\sigma^2\), sem depender da veracidade de \(H_0\), consiste em calcular para cada grupo a variância amostral (estimativa de \(\sigma^2\)) e tomar a média das várias estimativas que se obtêm.
Se pensarmos agora que as médias são todas iguais (\(H_0\) verdadeiro) estamos perante um
conjunto de \(g\) amostras todas da
mesma população.
Sabemos que:
\[ \text{Var}\left[\bar X\right] = \frac{\sigma^2_{\text{entre}}}{n} \]
e podemos obter uma “amostra” de \(g\) médias amostrais (uma para cada grupo).
Calculando a variância amostral desta “amostra” de médias amostrais temos uma estimativa de \(\sigma^2/n\).
Multiplicando por \(n\) temos uma estimativa de \(\sigma^2\).
Mas esta última estimativa só é boa se \(H_0\) for verdadeira. Senão fica muito inflacionada.
Assim, ao dividir a última estimativa pela primeira devemos obter um valor próximo de 1 se \(H_0\) for verdadeiro e muito maior que 1 caso contrário.
Suponha um delineamento entre participantes com \(g\ge3\) perfeitamente balanceado.
\[ F = \frac{\hat{\sigma}^2_{\text{entre}}}{\hat{\sigma}^2_{\text{dentro}}} = \frac{n\,\hat\sigma^2_{\bar X}}{\tfrac{1}{g}\sum_{i=1}^{g}\hat\sigma^2_{X_i}}=\frac{n S^2_{\bar X}}{\overline{S^2}} = n \;\frac{\dfrac{\sum_{i=1}^g (\bar X_i - \bar{\bar{X}})^2}{g-1}}{\dfrac{\sum_{i=1}^g S_i^2}{g}} \]
Note que a estatística de teste F é:
\[ F=n\;\dfrac{\text{variância das médias} }{\text{médias das variância}} \]
O seguinte estimador de \(\sigma^2\) não depende de \(H_0\): \(\hat{\sigma}^2_{\text{dentro}}\).
O seguinte estimador de \(\sigma^2\) depende de \(H_0\) e do Teorema do Limite Central: \(\hat{\sigma}^2_{\text{dentro}}\).
Outra forma de escrever:
\[ F = n \;\frac{\dfrac{\sum_{i=1}^g (\bar X_i - \bar X)^2}{g-1}} {\dfrac{\sum_{i=1}^g \sum_{j=1}^n (X_{ij} - \bar X_i)^2}{g(n-1)}} \;\;\sim\; F_{g-1,\,g(n-1)} \]
Definições auxiliares:
\[ \begin{align} \bar X_i &= \frac{\sum_{j=1}^n X_{ij}}{n}, \quad i=1,2,\ldots,g\\ \bar{{\bar X}} &= \frac{\sum_{i=1}^g \bar X_i}{g} \end{align} \]
com \(n = n_1 = n_2 = \cdots = n_g\).
Essa estimativa não depende da veracidade de \(H_0\) e reflete a variabilidade amostral dentro de cada condição independente.
Aplicação aos dados Nutrição3.rds
Médias e variâncias por grupo:
| Instrutor | Média | Variância |
|---|---|---|
| Brendon | 1287.50 | 37532.89 |
| McGuirk | 1246.25 | 20281.25 |
| Melissa | 1123.75 | 20491.78 |
- Variância das médias = 7253.65
- Média das variâncias = 26101.97
Estatística F (cálculo direto):
\[ F = n \; \frac{s^2_{\bar X}}{\overline{s^2}} = 20 \times \frac{7253.65}{26101.97} = 5.56 \]
Tabela ANOVA:
| Fonte | SQ | gl | QM | F | p |
|---|---|---|---|---|---|
| Instrutor | 290145.8 | 2 | 145072.9 | 5.56 | 0.0062 |
| Residual | 1487812.5 | 57 | 26102.0 |
Conclusão:
Como \(F=5.56 > F_{0.95}(2,57) \approx 3.16\), rejeitamos \(H_0\). Equivalentemente, rejeitamos \(H_0\), pois \(p=0.0062<0.05\).
Há evidência de que as médias populacionais de sódio/semana diferem significativamente entre os três instrutores.
Instructor n mean sd min Q1 median Q3 max
1 Brendon 20 1288 194 950 1150 1300 1400 1700
2 McGuirk 20 1246 142 1000 1144 1212 1350 1525
3 Melissa 20 1124 143 900 1006 1112 1231 1400Registered S3 methods overwritten by 'car':
method from
hist.boot FSA
confint.boot FSA Anova Table (Type II tests)
Response: Sodium
Sum Sq Df F value Pr(>F)
Instructor 290146 2 5.56 0.0062 **
Residuals 1487812 57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
lm(formula = Sodium ~ Instructor, data = Dados)
Residuals:
Min 1Q Median 3Q Max
-337.5 -121.2 2.5 105.9 412.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1287.5 36.1 35.64 <2e-16 ***
InstructorMcGuirk -41.3 51.1 -0.81 0.4228
InstructorMelissa -163.8 51.1 -3.21 0.0022 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 162 on 57 degrees of freedom
Multiple R-squared: 0.163, Adjusted R-squared: 0.134
F-statistic: 5.56 on 2 and 57 DF, p-value: 0.00624Teste t independente é um teste F?
- Condições baleanceada: \(n_1=n_2=n\). Portanto, \(n_{\text{total}}=2n\).
Estatística de teste t de Student:
\[ t_{2(n-1)}=\frac{\bar X_1-\bar X_2}{\sqrt{\dfrac{S_1^2+S_2^2}{n}}} \]
\[ \bar X_1=\frac{1}{n}\sum_{j=1}^{n} X_{1j},\qquad \bar X_2=\frac{1}{n}\sum_{j=1}^{n} X_{2j} \]
\[ S_i^2=\frac{\sum_{j=1}^{n}\left(X_{ij}-\bar X_i\right)^2}{\,n-1\,},\quad i=1,2 \]
Estatística de teste F:
\[ F_{1,\,2(n-1)}=\frac{\hat\sigma^2_{\text{entre}}}{\hat\sigma^2_{\text{dentro}}} =\frac{n\,\hat\sigma^2_{\bar X}}{\tfrac12\sum_{i=1}^{2}\hat\sigma^2_{X_i}} =\frac{n\displaystyle\left[\dfrac{(\bar X_1-\bar X)^2+(\bar X_2-\bar X)^2}{2-1}\right]} {\displaystyle \dfrac{\sum_{j=1}^{n}(X_{1j}-\bar X_1)^2+\sum_{j=1}^{n}(X_{2j}-\bar X_2)^2}{2(n-1)}} \]
Demonstração
\[ \begin{align} F_{1,\,2(n-1)}&=\frac{n\left[(\bar X_1-\bar X)^2+(\bar X_2-\bar X)^2\right]} {\displaystyle \frac{\sum_{j=1}^{n}(X_{1j}-\bar X_1)^2+\sum_{j=1}^{n}(X_{2j}-\bar X_2)^2}{2(n-1)}}\\ &=\frac{n\left[\bar X_1^{\,2}+\bar X_2^{\,2}-2\bar X^{\,2}\right]} {\displaystyle \frac{\sum_{j=1}^{n}(X_{1j}-\bar X_1)^2+\sum_{j=1}^{n}(X_{2j}-\bar X_2)^2}{2(n-1)}}\\ &=\frac{n\left[\tfrac12(\bar X_1-\bar X_2)^2\right]} {\displaystyle \frac{\sum_{j=1}^{n}(X_{1j}-\bar X_1)^2+\sum_{j=1}^{n}(X_{2j}-\bar X_2)^2}{2(n-1)}}\\ &=\frac{(\bar X_1-\bar X_2)^2} {\displaystyle \frac{\sum_{j=1}^{n}(X_{1j}-\bar X_1)^2+\sum_{j=1}^{n}(X_{2j}-\bar X_2)^2}{n(n-1)}}\\ &=\frac{(\bar X_1-\bar X_2)^2}{\dfrac{S_1^2+S_2^2}{n}}\\ F_{1,\,2(n-1)}&=t_{2(n-1)}^{\,2} \end{align}\\ \quad\text{Q.E.D.} \]
ANOVA independente unifatorial
A ANOVA independente unifatorial é utilizada quando os participantes são avaliados em somente uma das condições experimentais, i.e., um delineamento entre participantes.
Suposições
- independência entre as unidades experimentais,
- normalidade da VD em todas as condições,
- homocedasticidade da VD entre todas as condições.
O teste analítico pode prescindir da suposição de normalidade para amostras maiores invocando o teorema central do limite (\(n\ge30\)). Podemos testar a suposição de normalidade da VD em cada condição com tamanhos de amostra pequenos. Enfatizamos que a suposição de normalidade é uma condição suficiente para a validade do teste estatístico.
O teste por bootstrapping pode ser usado para qualquer tamanho de amostra (\(n\ge2\)), e não exige a suposição de normalidade. Portanto, exige apenas a suposição de independência entre as observações.
O teste analítico pode prescindir da suposição de homocedasticidade se as condições são balanceadas.
Normalidade
“Vale a pena observar que a MANOVA é ainda um teste válido, mesmo com modestas violações na suposição de normalidade, particularmente quando os tamanhos dos grupos são iguais e existe um número razoável de participantes em cada grupo; por “razoável” entendemos que, em um delineamento completamente entre participantes, deve haver pelo menos 12 participantes por grupo[…]”
Dancey & Reidy, 2019, p. 472
Como ANOVA é um caso particular de MANOVA, podemos transpor desta referência e assumir o número de unidades experimentais igual a pelo menos 12 em cada grupo para considerar a suposição de normalidade não necessária.
Balanceamento e Homocedasticidade
O balanceamento não precisa ser estrito:
“Consideram-se grupos de dimensão semelhante quando o quociente entre a maior dimensão e a menor for inferior a 1,5.”
Pestana & Gageiro, 2008, p. 278
“A heterocedasticidade ocorre heuristicamente se (maior desvio-padrão)/(menor desvio-padrão) é maior que 2.”
Johnson & Wichern, 2007, p. 291; Moore, 1995, apud Norusis, 2009, p. 148
Quando a maior quantidade de unidades observacionais numa condição NÃO superar 1.5 vezes a condição de menor quantidade de unidades observacionais, então a suposição de homocedasticidade não precisa ser considerada para o teste ANOVA independente unifatorial .
“Em geral, quando você tiver tamanhos de amostras iguais, esse pressuposto [homocedasticidade] não será um grande problema.”
Dancey & Reidy, 2019, p. 472-3
ANOVA independente unifatorial de Fisher para condições balanceadas é adequada na situação de heterocedasticidade populacional.
ANOVA independente unifatorial de Fisher com correção de White e a ANOVA independente unifatorial de Welch funcionam para condições desbalanceadas, pois lidam com a situação de heterocedasticidade populacional.
Exemplo: SNAP-Ed entreparticipantes com 3 grupos
Este exemplo foi obtido de Salvatore S. Mangiafico: Summary and Analysis of Extension Program Evaluation in R, disponível em https://rcompanion.org/handbook/I_05.html.
O SNAP-Ed (Supplemental Nutrition Assistance Program Education) é um programa baseado em evidências que ajuda as pessoas a terem uma vida mais saudável.
O SNAP-Ed ensina às pessoas que usam ou qualificam para o SNAP uma boa nutrição e como fazer com que o seu dinheiro de alimentação se estenda ainda mais.
Os participantes do SNAP-Ed também aprendem a ser fisicamente ativos.
Brendon Small, Coach McGuirk e Melissa Robins fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.
Desde que as classes receberam diferentes programas de educação nutricional, eles querem ver se a ingestão média de sódio é a mesma para as três turmas.
Arquivo de dados
alfa <- 0.05
Dados <- data.frame(readxl::read_excel("Nutricao3.xlsx"))
Dados$Student <- factor(Dados$Student)
Dados$Instructor <- factor(Dados$Instructor,
labels=c("Brendon", "McGuirk", "Melissa"))
print(data.frame(Dados)) Instructor Student Sodium
1 Brendon a 1200
2 Brendon b 1400
3 Brendon c 1350
4 Brendon d 950
5 Brendon e 1400
6 Brendon f 1150
7 Brendon g 1300
8 Brendon h 1325
9 Brendon i 1425
10 Brendon j 1500
11 Brendon k 1250
12 Brendon l 1150
13 Brendon m 950
14 Brendon n 1150
15 Brendon o 1600
16 Brendon p 1300
17 Brendon q 1050
18 Brendon r 1300
19 Brendon s 1700
20 Brendon t 1300
21 McGuirk u 1100
22 McGuirk v 1200
23 McGuirk w 1250
24 McGuirk x 1050
25 McGuirk y 1200
26 McGuirk z 1250
27 McGuirk aa 1350
28 McGuirk ab 1350
29 McGuirk ac 1325
30 McGuirk ad 1525
31 McGuirk ae 1225
32 McGuirk af 1125
33 McGuirk ag 1000
34 McGuirk ah 1125
35 McGuirk ai 1400
36 McGuirk aj 1200
37 McGuirk ak 1150
38 McGuirk al 1400
39 McGuirk am 1500
40 McGuirk an 1200
41 Melissa ao 900
42 Melissa ap 1100
43 Melissa aq 1150
44 Melissa ar 950
45 Melissa as 1100
46 Melissa at 1150
47 Melissa au 1250
48 Melissa av 1250
49 Melissa aw 1225
50 Melissa ax 1325
51 Melissa ay 1125
52 Melissa az 1025
53 Melissa ba 950
54 Melissa bc 925
55 Melissa bd 1200
56 Melissa be 1100
57 Melissa bf 950
58 Melissa bg 1300
59 Melissa bh 1400
60 Melissa bi 1100Hipóteses nula e alternativa
As três turmas receberam diferentes programas de educação nutricional. A ingestão diária média de sódio é a mesma (populacionalmente) para os três programas?
\[ \begin{cases} H_0: \mu_{\text{Brendon}} = \mu_{\text{McGuirk}} = \mu_{\text{Robins}}\\ H_1: \exists \,\mu_i \ne \mu_j,\quad i \ne j,\quad i,j=1,2,3 \end{cases}\\ \alpha=5\% \]
|
Esta é uma forma matemática para escrever a hipótese alternativa, lida como “Existe pelo menos alguma média \(\mu_i\), diferente de outra média \(\mu_j\), com \(i\) e \(j\) assumindo os valores \(1\), \(2\) e \(3\)”. É dizer que \(\mu_1 \ne \mu_2\) OU \(\mu_1 \ne \mu_3\) OU \(\mu_2 \ne \mu_3\). No entanto, aparece, frequentemente, como: \[ H_1: \mu_1 \ne \mu_2 \ne \mu_3 \] tentando indicar as diferenças, mas esta forma não funciona porque sugere que a rejeição de \(H_0\) implica que todos os grupos são diferentes entre si, como \(\mu_1 \ne \mu_2\) E \(\mu_1 \ne \mu_3\) E \(\mu_2 \ne \mu_3\). Basta que haja pelo menos duas médias populacionais diferentes para que se rejeite \(H_0\) com teste F da ANOVA. Em português: \[ H_1: \text{Existe pelo menos duas médias populacionais diferentes} \] |
ANOVA independente unifatorial de Fisher
Análise descritiva e teste de suposições
O arquivo de dados está no formato wide: cada linha do
arquivo contém os dados de cada unidade experimental. No caso particular
do nosso estudo, o arquivo wide Nutricao3.rds tem
\(n=20\) linhas ou unidades
experimentais distintas.
Instructor
Brendon McGuirk Melissa
20 20 20
Instructor
Student Brendon McGuirk Melissa
a 1 0 0
aa 0 1 0
ab 0 1 0
ac 0 1 0
ad 0 1 0
ae 0 1 0
af 0 1 0
ag 0 1 0
ah 0 1 0
ai 0 1 0
aj 0 1 0
ak 0 1 0
al 0 1 0
am 0 1 0
an 0 1 0
ao 0 0 1
ap 0 0 1
aq 0 0 1
ar 0 0 1
as 0 0 1
at 0 0 1
au 0 0 1
av 0 0 1
aw 0 0 1
ax 0 0 1
ay 0 0 1
az 0 0 1
b 1 0 0
ba 0 0 1
bc 0 0 1
bd 0 0 1
be 0 0 1
bf 0 0 1
bg 0 0 1
bh 0 0 1
bi 0 0 1
c 1 0 0
d 1 0 0
e 1 0 0
f 1 0 0
g 1 0 0
h 1 0 0
i 1 0 0
j 1 0 0
k 1 0 0
l 1 0 0
m 1 0 0
n 1 0 0
o 1 0 0
p 1 0 0
q 1 0 0
r 1 0 0
s 1 0 0
t 1 0 0
u 0 1 0
v 0 1 0
w 0 1 0
x 0 1 0
y 0 1 0
z 0 1 0
Instructor n mean sd min Q1 median Q3 max
1 Brendon 20 1287.5 193.7 950 1150.0 1300.0 1400.0 1700
2 McGuirk 20 1246.2 142.4 1000 1143.8 1212.5 1350.0 1525
3 Melissa 20 1123.8 143.1 900 1006.2 1112.5 1231.2 1400
boxplot
gplots::plotmeans
NULL
Normality
--------
Instructor = Brendon
Shapiro-Wilk normality test
data: Sodium
W = 0.97212, p-value = 0.7989
--------
Instructor = McGuirk
Shapiro-Wilk normality test
data: Sodium
W = 0.96772, p-value = 0.7062
--------
Instructor = Melissa
Shapiro-Wilk normality test
data: Sodium
W = 0.95835, p-value = 0.5114
--------
p-values adjusted by the Holm method:
unadjusted adjusted
Brendon 0.79894 1
McGuirk 0.70625 1
Melissa 0.51138 1
NULL
Homoscedasticity
Fligner-Killeen test of homogeneity of variances
data: Sodium by Instructor
Fligner-Killeen:med chi-squared = 0.95543, df = 2, p-value = 0.6202Testes omnibus e post-hoc
O arquivo de dados está no formato wide: cada linha do arquivo contém os dados de cada unidade experimental.
Sob a hipótese nula \(H_0: \mu_1 = \mu_2 = \cdots = \mu_g\), a estatística de teste F, segue a distribuição F de Fisher-Snedecor com os números de graus de liberdade do efeito do fator A, \(df_A=g-1\), e do efeito residual, \(df_E=\sum_{i=1}^{g}{(n_i-1)}=(n-1)-(g-1)=n-g\), no numerador e denominador, respectivamente, sendo que \(g\) é o número de condições independentes e \(n = \sum_{i=1}^g n_i\) é número total de unidades experimentais. \(df_A\) quantifica a complexidade do modelo (GLM univariado) e \(df_E\) quantifica o tamanho de amostra efetivo, sendo que \(n\) é o tamanho de amostra nominal.
\[ \begin{align} F &\underset{H_0}{\sim} F_{df_A, \, df_E}\\ F &\underset{H_0}{\sim} F_{g-1, \, \sum_{i=1}^{g}{(n_i-1)}}\\ F &\underset{H_0}{\sim} F_{g-1, \, (n-1)-(g-1)}\\ F &\underset{H_0}{\sim} F_{g-1, \, n-g}\\ \end{align} \]
Se \(df_E \ge 30\), pode-se recorrer ao Teorema Central do Limite, dispensando-se a suposição de multinormalidade para a validade do teste de hipótese nula omnibus.
Se \(df_E < 30\), pode-se recorrer ao argumento de que a distribuição multinormal em cada condição independente é uma condição suficiente para a validade do teste de hipótese nula omnibus. Nenhum teste estatístico de normalidade multivariada ou univariada é prova de normalidade, sendo possível apenas desprovar a normalidade por meio de teste estatístico. Dessa forma, mesmo em situação de amostra pequena, o teste omnibus pode ser válido.
O código R que testa a hipótese nula omnibus e realiza os
testes post-hoc da ANOVA de Fisher é (demo_ANOVA1f_indep_Fisher_sodio_Parte2.R).
modelo <-
lm(Sodium ~ Instructor, data=Dados)
car::Anova(modelo)
jmv::ANOVA(dep = “Sodium”, factors = “Instructor”, modelTes = TRUE, postHoc= ~ Instructor, emMeans = ~ Instructor, emmPlots = TRUE, emmTables = TRUE, data = Dados)
Instructor
Brendon McGuirk Melissa
20 20 20
Instructor n mean sd min Q1 median Q3 max
1 Brendon 20 1288 194 950 1150 1300 1400 1700
2 McGuirk 20 1246 142 1000 1144 1212 1350 1525
3 Melissa 20 1124 143 900 1006 1112 1231 1400
Fisher's One-way ANOVA
Statistical analysis: omnibus test
ANOVA
Anova Table (Type II tests)
Response: Sodium
Sum Sq Df F value Pr(>F)
Instructor 290146 2 5.5579 0.006235 **
Residuals 1487812 57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
lm(formula = Sodium ~ Instructor, data = Dados)
Residuals:
Min 1Q Median 3Q Max
-337.5 -121.2 2.5 105.9 412.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1287.50 36.13 35.639 < 2e-16 ***
InstructorMcGuirk -41.25 51.09 -0.807 0.42279
InstructorMelissa -163.75 51.09 -3.205 0.00221 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 161.6 on 57 degrees of freedom
Multiple R-squared: 0.1632, Adjusted R-squared: 0.1338
F-statistic: 5.558 on 2 and 57 DF, p-value: 0.006235
R^2 = eta^2 =0.16
# Effect Size for ANOVA (Type II)
Parameter | Eta2 | 95% CI | interpret
--------------------------------------------------
Instructor | 0.1632 | [0.0166, 0.3265] | large
Fisher's One-way ANOVA
Statistical analysis: Post hoc tests$emmeans
Instructor emmean SE df lower.CL upper.CL
Brendon 1288 36.1 57 1215 1360
McGuirk 1246 36.1 57 1174 1319
Melissa 1124 36.1 57 1051 1196
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Brendon - McGuirk 41.2 51.1 57 0.807 0.4228
Brendon - Melissa 163.8 51.1 57 3.205 0.0066
McGuirk - Melissa 122.5 51.1 57 2.398 0.0396
P value adjustment: holm method for 3 tests 
Instructor emmean SE df lower.CL upper.CL .group
Melissa 1124 36.1 57 1035 1213 a
McGuirk 1246 36.1 57 1157 1335 b
Brendon 1288 36.1 57 1198 1377 b
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates
P value adjustment: holm method for 3 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Teste post hoc: Dunnett contrast estimate SE df t.ratio p.value
McGuirk - Brendon -41.2 51.1 57 -0.807 0.4228
Melissa - Brendon -163.8 51.1 57 -3.205 0.0044
P value adjustment: holm method for 2 tests 
Teste post hoc: Consecutivo contrast estimate SE df t.ratio p.value
McGuirk - Brendon -41.2 51.1 57 -0.807 0.4228
Melissa - McGuirk -122.5 51.1 57 -2.398 0.0396
P value adjustment: holm method for 2 tests 
One-way analysis of means
data: Sodium and Instructor
F = 5.56, num df = 2, denom df = 57, p-value = 0.0062
Análise de significância prática: tamanho de efeito
# Effect Size for ANOVA
Eta2 | 95% CI | interpret
-------------------------------------
0.1632 | [0.0166, 0.3265] | large
One-Way Analysis of Variance (alpha = 0.05)
-------------------------------------------------------------
data : Sodium and Instructor
statistic : 5.557929
num df : 2
denom df : 57
p.value : 0.006235321
Result : Difference is statistically significant.
-------------------------------------------------------------
Holm Correction (alpha = 0.05)
-----------------------------------------------------
Level (a) Level (b) p.value No difference
1 Brendon McGuirk 0.44768959 Not reject
2 Brendon Melissa 0.01280080 Reject
3 McGuirk Melissa 0.01992282 Reject
-----------------------------------------------------
ANOVA
ANOVA - Sodium
─────────────────────────────────────────────────────────────────────────────────
Sum of Squares df Mean Square F p
─────────────────────────────────────────────────────────────────────────────────
Overall model 290145.8 2 145072.92 5.557929 0.0062353
Instructor 290145.8 2 145072.92 5.557929 0.0062353
Residuals 1487812.5 57 26101.97
─────────────────────────────────────────────────────────────────────────────────
POST HOC TESTS
Post Hoc Comparisons - Instructor
──────────────────────────────────────────────────────────────────────────────────────────────────────
Instructor Instructor Mean Difference SE df t p-tukey
──────────────────────────────────────────────────────────────────────────────────────────────────────
Brendon - McGuirk 41.25000 51.09009 57.00000 0.8073973 0.7000083
- Melissa 163.75000 51.09009 57.00000 3.2051225 0.0061739
McGuirk - Melissa 122.50000 51.09009 57.00000 2.3977252 0.0510271
──────────────────────────────────────────────────────────────────────────────────────────────────────
Note. Comparisons are based on estimated marginal means
ESTIMATED MARGINAL MEANS
INSTRUCTOR
Estimated Marginal Means - Instructor
──────────────────────────────────────────────────────────────
Instructor Mean SE Lower Upper
──────────────────────────────────────────────────────────────
Brendon 1287.500 36.12615 1215.159 1359.841
McGuirk 1246.250 36.12615 1173.909 1318.591
Melissa 1123.750 36.12615 1051.409 1196.091
────────────────────────────────────────────────────────────── 
Pr(>F) é o valor p nas saídas.
Como p value = 0.006235 omnibus da regressão
(igual ao valor p da tabela ANOVA do teste do efeito do fator
Instructor) é menor que 0.05, rejeita-se a hipótese nula
omnibus (que é a mesma hipótese nula de ausência do efeito do
fator).
Multiple R-squared \(R^2=0.1632\) é a estimativa pontual do
tamanho de efeito \(\eta^2\) de Cohen
omnibus (que é também a estimativa pontual do tamanho de efeito
do fator Instructor). Note que o intervalo de confiança
bilateral de 95% do tamanho de efeito de Instructor sobre
Sodium, \(\eta^2\), está
entre 1.66% e 32.65%.
Se \(H_0\) é rejeitada, podemos
localizar as diferenças com testes post-hoc (par-a-par
simultâneas) com os contrastes de Holm. A conclusão é que
Melissa tem média populacional diferente e menor.
Exemplo de relatório para a ANOVA independente unifatorial de Fisher
Há 60 participantes no estudo com delineamento entreparticipantes com
três condições, sendo que nenhuma delas é de controle. Os participantes
foram distribuídos aleatoriamente e balanceadamente nos três grupos. As
suposições de normalidade e homocedasticidade da VD Sodium
nos grupos foram consideradas válidas. A análise por meio da ANOVA
independente unifatorial de Fisher mostra que a menor ingestão de sódio
diária média foi observada entre estudantes submetidos ao programa
aplicado por Melissa Robins. As médias populacionais de ingestão de
sódio dos estudantes de Brendon Small e do Coach McGuirk são
semelhantes.
A análise de variância de um fator fixo entre participantes mostrou
que o efeito fixo Instructor é estatisticamente
significante, pois o teste omnibus produziu
F(2,57) = 5.56 e p = 0.0062. O tamanho do
efeito de Instructor é expresso por eta ao quadrado de
Cohen, sendo que seu valor é igual a 16,3% e seu
IC95% = [1.7%, 32.7%]. Portanto, 16% da variância da
ingesta de sódio é explicada pelo programa adotado pelos instrutores. Os
testes post hoc de Holm encontraram diferenças estatisticamente
significantes entre os programas adotados por Melissa e McGuirk e
Melissa e Brendon. Não se observou diferença estatisticamente
significante entre McGuirk e Brendon.
ANOVA independente unifatorial de Welch
- Welch (1951); Games & Howell (1976)
O arquivo de dados está no formato wide: cada linha do arquivo contém os dados de cada unidade experimental.
Sob \(H_0:\ \mu_1=\mu_2=\cdots=\mu_g\), a estatística de Welch é
\[ F_W=\frac{\displaystyle \frac{1}{g-1}\sum_{i=1}^g w_i\left(\bar X_i-\bar X_w\right)^2}{\displaystyle 1+\frac{2(g-2)}{g^2-1}\sum_{i=1}^g\frac{(1-w_i/w)^2}{n_i-1}}, \]
com \(w_i=\dfrac{n_i}{s_i^2}\), \(w=\sum_{i=1}^g w_i\) e \(\bar X_w=\dfrac{\sum_{i=1}^g w_i\bar X_i}{W}\).
A distribuição sob \(H_0\) é
\[ F_W \sim F_{g-1,\ \nu} \]
em que o numerador permanece \(g-1\) e o denominador é
\[ \nu=\frac{g^2-1}{3\displaystyle\sum_{i=1}^g\frac{(1-w_i/w)^2}{n_i-1}} \]
ANOVA de Welch F df1 df2 p
5.765 2 37.4 0.006569
One-way analysis of means (not assuming equal variances)
data: Sodium and Instructor
F = 5.7655, num df = 2.000, denom df = 37.396, p-value = 0.006569Para \(g\ge 3\), mesmo sob homocedasticidade amostral (\(s_i^2=s^2\)) e amostras balanceadas (\(n_i= m\)), a aproximação de Satterthwaite dá
\[ \nu=\frac{g^2-1}{3\sum_{i=1}^g \dfrac{(1-w_i/W)^2}{n_i-1}}= \frac{g+1}{3(g-1)}\,(n-g) \]
Máximo:
\[ \nu_{\max}=\frac{g+1}{3(g-1)}(n-g)\quad(\text{homocedástico amostral e balanceado}) \]
Para \(g=3\), \(n=120\): \(\nu_{\max}=\dfrac{4}{6}\cdot117=78\).
Mínimo (infimum teórico) sob \(n_i\ge2\) e pesos extremos (\(w_1/w\to1\), demais \(\to0\)):
\[ \nu_{\min}=\frac{g^2-1}{3(g-1)}=\frac{g+1}{3} \]
Para \(g=3\): \(\nu_{\min}=4/3\approx1.33\).
Com \(n_i\) maiores (e.g.: 4, 4, 120) e variâncias não “infinitas”, o menor \(\nu\) observado sobe para aproximadamente 4.
Para \(g=3\), \(2\) não é o limite inferior; o
infimum é \(4/3\). E \(117\) é o df clássico (\(N-g\)), não o de Welch; para Welch o máximo
é \(78\) no seu cenário.
Teórico: nu_max = 78 nu_min (infimum) = 1.333333
Simulado: nu_max ≈ 77.789 nu_min ≈ 1.333 (seed= 2985 )
N - g = 117 ; nu_max_teor = 78 ANOVA independente unifatorial de Welch demo_ANOVA1f_indep_Welch_sodio.R pode
ser executado pelas funções jmv::anovaOneW ou
oneway.test e rstatix::games_howell_test.
oneway.test(Sodium ~ Instructor, var.equal=FALSE, data=Dados)
rstatix::games_howell_test(Sodium ~ Instructor, conf.level=1-alfa, detailed=FALSE, data=Dados)
jmv::anovaOneW(data = Dados, dep = “Sodium”, group = “Instructor”, welch = TRUE, phMethod = “gamesHowell”, phTest = TRUE, phFlag = TRUE, desc = FALSE)
Instructor
Brendon McGuirk Melissa
20 20 20
ANOVA unifatorial independente de Welch
Análise de significância estatística: testes omnibus e posthoc
One-way analysis of means (not assuming equal variances)
data: Sodium and Instructor
F = 5.77, num df = 2.0, denom df = 37.4, p-value = 0.0066
.y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
1 Sodium Brendon McGuirk -41 -173 90 0.725 ns
2 Sodium Brendon Melissa -164 -296 -32 0.012 *
3 Sodium McGuirk Melissa -122 -233 -12 0.026 *
Análise de significância prática: tamanho de efeito`var.equal = FALSE` - effect size is an approximation.# Effect Size for ANOVA
Eta2 | 95% CI | interpret
-------------------------------------
0.2357 | [0.0251, 0.4351] | large
ONE-WAY ANOVA
One-Way ANOVA (Welch's)
──────────────────────────────────────────────────────
F df1 df2 p
──────────────────────────────────────────────────────
Sodium 5.765465 2 37.39637 0.0065692
──────────────────────────────────────────────────────
POST HOC TESTS
Games-Howell Post-Hoc Test – Sodium
─────────────────────────────────────────────────────────────────────
Brendon McGuirk Melissa
─────────────────────────────────────────────────────────────────────
Brendon Mean difference — 41.25000 163.7500
t-value — 0.7672235 3.040115
df — 34.89307 34.98266
p-value — 0.7253215 0.0120952
McGuirk Mean difference — 122.5000
t-value — 2.713093
df — 37.99899
p-value — 0.0263388
Melissa Mean difference —
t-value —
df —
p-value —
─────────────────────────────────────────────────────────────────────
Note. * p < .05, ** p < .01, *** p < .001
.y. n statistic DFn DFd p method
1 Sodium 60 5.77 2 37.4 0.007 Welch ANOVA
.y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
1 Sodium Brendon McGuirk -41 -173 90 0.725 ns
2 Sodium Brendon Melissa -164 -296 -32 0.012 *
3 Sodium McGuirk Melissa -122 -233 -12 0.026 *
Welch's Heteroscedastic F Test (alpha = 0.05)
-------------------------------------------------------------
data : Sodium and Instructor
statistic : 5.765465
num df : 2
denom df : 37.39637
p.value : 0.006569226
Result : Difference is statistically significant.
-------------------------------------------------------------
Holm Correction (alpha = 0.05)
-----------------------------------------------------
Level (a) Level (b) p.value No difference
1 Brendon McGuirk 0.44810920 Not reject
2 Brendon Melissa 0.01337414 Reject
3 McGuirk Melissa 0.01992301 Reject
----------------------------------------------------- A conclusão, neste exemplo, é similar à obtida com o teste de ANOVA de Fisher: o programa adotado por Melissa obteve ingestões diárias médias de sódio significantemente menores que a dos outros dois instrutores.
ANOVA independente unifatorial de Fisher-White
ANOVA independente unifatorial de Fisher-White demo_ANOVA1f_indep_FisherWhite.R.
modelo <-
lm(Sodium~Instructor, data=TH)
car::Anova(modelo, white.adjust=“hc2”)
Instructor
Brendon McGuirk Melissa
20 20 20
Fisher-White's One-way ANOVA
----------------------------------
Statistical analysis: Omnibus test
----------------------------------Coefficient covariances computed by hccm()Analysis of Deviance Table (Type II tests)
Response: Sodium
Df F Pr(>F)
Instructor 2 5.8682 0.004815 **
Residuals 57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
estimatr::lm_robust(formula = Sodium ~ Instructor, data = Dados,
se_type = "HC2")
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 1287.50 43.32 29.7205 2.185e-36 1200.8 1374.25 57
InstructorMcGuirk -41.25 53.77 -0.7672 4.461e-01 -148.9 66.41 57
InstructorMelissa -163.75 53.86 -3.0401 3.568e-03 -271.6 -55.89 57
Multiple R-squared: 0.1632 , Adjusted R-squared: 0.1338
F-statistic: 5.868 on 2 and 57 DF, p-value: 0.004815
R^2 = eta^2 omnibus = 0.16[1] "large"
(Rules: field2013)
$emmeans
Instructor emmean SE df lower.CL upper.CL
Brendon 1288 43.3 57 1201 1374
McGuirk 1246 31.8 57 1182 1310
Melissa 1124 32.0 57 1060 1188
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Brendon - McGuirk 41.2 53.8 57 0.767 0.4461
Brendon - Melissa 163.8 53.9 57 3.040 0.0107
McGuirk - Melissa 122.5 45.2 57 2.713 0.0176
P value adjustment: holm method for 3 tests 
Instructor emmean SE df lower.CL upper.CL .group
Melissa 1124 32.0 57 1045 1203 a
McGuirk 1246 31.8 57 1168 1325 b
Brendon 1288 43.3 57 1181 1394 b
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates
P value adjustment: holm method for 3 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Teste post hoc: Dunnett contrast estimate SE df t.ratio p.value
McGuirk - Brendon -41.2 53.8 57 -0.767 0.4461
Melissa - Brendon -163.8 53.9 57 -3.040 0.0071
P value adjustment: holm method for 2 tests 
Teste post hoc: Consecutivo contrast estimate SE df t.ratio p.value
McGuirk - Brendon -41.2 53.8 57 -0.767 0.4461
Melissa - McGuirk -122.5 45.2 57 -2.713 0.0176
P value adjustment: holm method for 2 tests 
A conclusão, neste exemplo, é similar à obtida com o testes de ANOVA de Fisher e de Welch: o programa adotado por Melissa obteve ingestões diárias médias de sódio significantemente menores que a dos outros dois instrutores.
ANOVA independente unifatorial por bootstrapping
demo_ANOVA1f_indep_Bootstrap_sodio.R é
a versão com bootstrapping por meio das funções
lmboot::ANOVA.boot (sem teste post hoc).
lmboot::ANOVA.boot(Sodium ~ Instructor, B=B, type=“residual”, wild.dist=“normal”, seed=123, data=TH, keep.boot.resp=FALSE)
Instructor
Brendon McGuirk Melissa
20 20 20
Bootstrap One-Way ANOVA
F(2, 57) = 5.78, p = 0.0052
(10000 bootstrap samples)
Effect size analysis
eta^2 = 0.17\nTamanho de efeito: estimativa pontual
"large"
(Rules: field2013)Uma outra possível solução (com teste post hoc) são as
funções WRS2::t1waybt e WRS2::mcppb20, e
onewaytests::onewaytests. Implementamos em demo_WRS2_ANOVA_independente.R.
onewaytests::onewaytests(Sodium ~ Instructor, method=“gtb”, alpha=alfa, na.rm=TRUE, N=B, verbose=TRUE, data=Dados)onewaytests::paircomp(fit, adjust.method = “holm”)
WRS2::t1waybt(Sodium ~ Instructor, data=Dados, tr=0, nboot=B))WRS2::mcppb20(Sodium~Instructor, data=Dados, tr=0.2, level=1-alfaBonf, nboot=B)
Instructor
Brendon McGuirk Melissa
20 20 20
Generalized Test Equivalent to Parametric Bootstrap Test (size close to intended) (alpha = 0.05)
----------------------------------------------------------------------------------------------------
data : Sodium and Instructor
p.value : 0.006664
Result : Difference is statistically significant.
----------------------------------------------------------------------------------------------------
Holm Correction (alpha = 0.05)
-----------------------------------------------------
Level (a) Level (b) p.value No difference
1 Brendon McGuirk 0.447798 Not reject
2 Brendon Melissa 0.013284 Reject
3 McGuirk Melissa 0.019962 Reject
-----------------------------------------------------
Bootstrapped independent one-Way ANOVA
Call:
WRS2::t1waybt(formula = Sodium ~ Instructor, data = Dados, tr = 0,
nboot = B)
Effective number of bootstrap samples was 100000.
Test statistic: 5.7655
p-value: 0.00686
Variance explained: 0.241
Effect size: 0.491
Call:
WRS2::mcppb20(formula = Sodium ~ Instructor, data = Dados, tr = 0.2,
nboot = B, level = 1 - alfaBonf)
psihat ci.lower ci.upper p-value
Brendon vs. McGuirk 50.0 -87.50000 177.0833 0.36933
Brendon vs. Melissa 162.5 27.08333 302.0833 0.00539
McGuirk vs. Melissa 112.5 -4.16667 241.6667 0.02193Usamos duas funções neste último código R. A primeira é um ANOVA independente unifatorial que utiliza bootstrapping e computa \(\eta^2\) (chamado de variância explicada). A segunda faz um teste post hoc, também por bootstrapping, mas que tem um problema: obriga uma média aparada de 20% que não conseguimos desligar; é provavelmente um erro de implementação, e esperamos que o autor o resolva em versões futuras.
As conclusões são quase as mesmas: os valores p não são ajustados na saída da função, então as multiplicamos por 3 (correção de Bonferroni) para o gráfico. Com isto, a diferença entre McGuirk e Melissa deixou de ser significante (decisão de acordo com o intervalo de confiança 95% que inclui a diferença igual a zero).
Comparativo (ANOVA independente)
Para facilitar a comparação das três variantes de ANOVA independente unifatorial (ominibus) feitas acima:
| ANOVA | F | df | p | Função |
|---|---|---|---|---|
| Fisher | 5.56 | 2, 57 | 0.0062 | jmv::ANOVA |
| 5.56 | 2, 57 | 0.0062 | lm, car::Anova |
|
| Welch | 5.77 | 2, 37.4 | 0.0066 | oneway.test |
| 5.77 | 2, 37.4 | 0.0066 | rstatix::welch_anova_test |
|
| 5.77 | 2, 37.4 | 0.0066 | onewaytests::onewaytests |
|
| Fisher-White | 5.87 | 2, 57 | 0.0048 | lm, car::Anova |
| Bootstrap | 5.78 | 2, 57 | 0.0052 | lmboot::ANOVA.boot |
| 0.0065 | onewaytests::onewaytests |
|||
| 5.77 | 0.0080 | WRS2::t1waybt |
Para os testes post hoc:
- Fisher (
jmv::ANOVA)
| p | |
|---|---|
| Brendon - McGuirk | 0.7000 |
| Brendon - Melissa | 0.0062 |
| McGuirk - Melissa | 0.0510 |
- Fisher (
lm:emmeans::emmeans)
| p | |
|---|---|
| Brendon - McGuirk | 0.4228 |
| Brendon - Melissa | 0.0066 |
| McGuirk - Melissa | 0.0396 |
- Welch (
oneway.test:rstatix::games_howell_test)
| p | |
|---|---|
| Brendon - McGuirk | 0.725 |
| Brendon - Melissa | 0.012 |
| McGuirk - Melissa | 0.026 |
- Welch (
rstatix::welch_anova_test:rstatix::games_howell_test)
| p | |
|---|---|
| Brendon - McGuirk | 0.725 |
| Brendon - Melissa | 0.012 |
| McGuirk - Melissa | 0.026 |
- Welch (
onewaytests::onewaytests:onewaytests::paircomp)
| p | |
|---|---|
| Brendon - McGuirk | 0.4481 |
| Brendon - Melissa | 0.0134 |
| McGuirk - Melissa | 0.0199 |
- Fisher-White (
lm:emmeans::emmeans)
| p | |
|---|---|
| Brendon - McGuirk | 0.4461 |
| Brendon - Melissa | 0.0107 |
| McGuirk - Melissa | 0.0176 |
(os mesmos valores obtidos por ANOVA de Fisher)
- Bootstrap (
lmboot::ANOVA.boot: ?)
post hoc não disponível
- Bootstrap (
onewaytests::onewaytests:onewaytests::paircomp)
| p | |
|---|---|
| Brendon - McGuirk | 0.4482 |
| Brendon - Melissa | 0.0134 |
| McGuirk - Melissa | 0.0200 |
- Bootstrap (
WRS2::t1waybt:WRS2::mcppb20)
| p | |
|---|---|
| Brendon - McGuirk | 0.3691 |
| Brendon - Melissa | 0.0047 |
| McGuirk - Melissa | 0.0219 |
ANOVA independente SEM dados brutos
É muito comum, em publicações, que somente tenhamos acesso às medidas-resumo (número de participantes, média, desvio-padrão e correlação).
Nestes casos, os códigos R acima não são utilizáveis.
Os testes F de ANOVA independente unifatorial de Fisher e de Welch sem os dados brutos e também as comparações aos pares pode ser feito com a implementação da seguinte função:
source("eiras.OneWayANOVA.Fisher.WORD.R")
n <- c(10, 10, 10, 10)
media <- c(0.2760, 0.1956, 0.1495, 0.3997)
dp <- c(0.4095, 0.1808, 0.2451, 0.2382)
names(media) <- c("A","B","C","D") # rótulos dos grupos (opcional)
res <- OneWayANOVA.Fisher.WORD(n, media, dp,
alpha = 0.05,
posthoc = TRUE)
print(res$anova, digits=4) Source SS df MS F p eta2
1 Between 0.3604 3 0.1201 1.515 0.2272 0.1121
2 Within 2.8547 36 0.0793 NA NA NA
3 Total 3.2151 39 NA NA NA NA Group.1 Group.2 Mean.1 st.dev.1 n.1 Mean.2 st.dev.2 n.2 diff SE t df
1 A B 0.28 0.41 10 0.20 0.18 10 0.080 0.13 0.64 36
2 A C 0.28 0.41 10 0.15 0.25 10 0.127 0.13 1.00 36
3 A D 0.28 0.41 10 0.40 0.24 10 -0.124 0.13 0.98 36
4 B C 0.20 0.18 10 0.15 0.25 10 0.046 0.13 0.37 36
5 B D 0.20 0.18 10 0.40 0.24 10 -0.204 0.13 1.62 36
6 C D 0.15 0.25 10 0.40 0.24 10 -0.250 0.13 1.99 36
p.raw p.adj sig. eta2
1 0.527 1.00 ns 0.0112
2 0.322 1.00 ns 0.0273
3 0.333 1.00 ns 0.0261
4 0.716 1.00 ns 0.0037
5 0.114 0.68 ns 0.0680
6 0.055 0.33 ns 0.0988
ANOVA independente unifatorial de Fisher sem dados brutos Grupo media dp n
1 A 0.28 0.41 10
2 B 0.20 0.18 10
3 C 0.15 0.25 10
4 D 0.40 0.24 10factor <- "Grupo"
mean <- "media"
sd <- "dp"
n <- "n"
alpha <- 0.05
resF <- HH::anovaMean(object=Dados$Grupo,
n=Dados$n,
ybar=Dados$media,
s=Dados$dp,
ylabel="SDB")Registered S3 method overwritten by 'HH':
method from
print.ancova WRS2Warning in ybar - (ybar %*% n)/sum(n): Reciclar uma array de comprimento 1 na aritmética de vetor-array foi descontinuado.
Em vez disso, use c() ou as.vector().Analysis of Variance Table
Response: SDB
Terms added sequentially (first to last)
Df Sum of Sq Mean Sq F value Pr(F)
SDB 3 0.3604 0.1201 1.515 0.227
Residuals 36 2.8547 0.0793 source("eiras.OneWayANOVA.Welch.WORD.R")
cat("\nANOVA independente unifatorial de Welch sem dados brutos\n")
ANOVA independente unifatorial de Welch sem dados brutos# vetores resumo
n <- c(10,10,10,10)
media <- c(0.2760,0.1956,0.1495,0.3997)
sd <- c(0.4095,0.1808,0.2451,0.2382)
names(media) <- c("A","B","C","D")
# omnibus: Welch + post-hoc: Games–Howell
res <- OneWayWelchSummary(n, media, sd, labels = names(media), posthoc = "games-howell")
print(res$omnibus_title)[1] "ANOVA independente unifatorial de Welch" Source SS df1 df2 MS F p eta2
1 Groups 6.535 3 19.47 2.178 2.039 0.1417 0.239[1] "Post hoc: Games–Howell" Group.1 Group.2 Mean.1 st.dev.1 n.1 Mean.2 st.dev.2 n.2 diff SE q
A B 0.276 0.409 10 0.196 0.181 10 0.0804 0.1416 0.568
A C 0.276 0.409 10 0.149 0.245 10 0.1265 0.1509 0.838
A D 0.276 0.409 10 0.400 0.238 10 -0.1237 0.1498 0.826
B C 0.196 0.181 10 0.149 0.245 10 0.0461 0.0963 0.479
B D 0.196 0.181 10 0.400 0.238 10 -0.2041 0.0946 2.158
C D 0.149 0.245 10 0.400 0.238 10 -0.2502 0.1081 2.315
df p conf.low conf.high sig.
12.4 0.940 -0.338 0.4987 ns
14.7 0.836 -0.309 0.5625 ns
14.5 0.841 -0.557 0.3100 ns
16.6 0.963 -0.228 0.3206 ns
16.8 0.175 -0.473 0.0651 ns
18.0 0.132 -0.556 0.0553 nsPlanejamento de ANOVA independente unifatorial de Fisher
O tamanho do efeito pode ser expresso por eta ao quadrado de Cohen (\(\eta^2\)).
Boa parte dos valores relatados por Ellis (2010) são originais de Cohen (1992).
Cohen, 1992
Do mesmo trabalho, a tabela 2 indica o planejamento dos estudos:
Esta tabela pode ser reproduzida com as funções do pacote
pwr. Por exemplo, para 3 grupos, nivel de significância de
5%, poder de 80%, com tamanho de efeito intermediário, obtemos:
Balanced one-way analysis of variance power calculation
k = 3
n = 52.3966
f = 0.25
sig.level = 0.05
power = 0.8
NOTE: n is number in each groupPara poder de 90%:
Balanced one-way analysis of variance power calculation
k = 3
n = 68.49707
f = 0.25
sig.level = 0.05
power = 0.9
NOTE: n is number in each groupO valor f de Cohen, necessário para esta função, precisa
ser calculado em função do tamanho de efeito desejado. Supondo que você
esteja mais familiarizado com \(\eta^2\) (igual ao \(R^2\)), a conversão é dada por:
\[ f = \sqrt{\dfrac{\eta^2}{1-\eta^2}} \]
Basta escolher a partir de uma das tabelas o tamanho de efeito que se deseja detectar.
Esta função pode receber os seguintes parâmetros:
pwr.anova.test(k = NULL, n = NULL, f = NULL, sig.level = 0.05,
power = NULL)
Nos exemplos acima, o valor dado em n, que não foi
fornecido na chamada da função, foi calculado.
Além deste planejamento mais habitual (achar o tamanho da amostra), podemos fornecer quaisquer quatro dos cinco valores, para obter o quinto, considerado como incógnita e calculado. E.g.: suponha que dispomos de 50 pacientes e pretendemos medir tamanho de efeito intermediário, com nível de significância de 5% e poder de 80%:
Balanced one-way analysis of variance power calculation
k = 3.268196
n = 50
f = 0.25
sig.level = 0.05
power = 0.8
NOTE: n is number in each groupNecessitaremos de 3 a 4 grupos. Na verdade, se arredondarmos para 3, o poder será um pouco menor que 80%. Quanto?
Balanced one-way analysis of variance power calculation
k = 3
n = 50
f = 0.25
sig.level = 0.05
power = 0.7795803
NOTE: n is number in each groupObteremos poder próximo a 78%.
Tradicionalmente exigia-se, para a ANOVA independente unifatorial sem reamostragem, que a distribuição populacional da variável fosse normal para cada condição e que existisse homocedasticidade entre os grupos. As suposições são:
- Independência das observações dentro e entre condições: delineamento entreparticipantes.
- Normalidade da VD se a amostra é pequena (menos de 12 observações em cada grupo).
- Homocedasticidade da VD se os grupos são desbalanceados.
Normalidade e homocedasticidade podem ser testadas. Quando não eram atendidas, alguns autores recomendam o teste não-paramétrico H de Kruskal-Wallis.
No entanto, este teste não-paramétrico não deve ser usado. O teste de Kruskal-Wallis com VD intervalar não precisa de normalidade, mas precisa de homocedasticidade para testar diferenças de pseudo-medianas; para testar diferença de médias, também precisa de simetria da VD em cada condição independente. Testes não paramétricos não são livres de suposições.
Quando o tamanho da amostra é pequeno e há heterocedasticidade ou desbalanceamento da VD, há alternativas de ANOVA independente unifatorial que podem funcionar: Fisher-White, Welch ou por reamostragem (bootstrapping).
Quando o tamanho da amostra é grande (tamanho de amostra total maior que 30 e menor grupo igual a 12 observações), o teorema central do limite permitirá que ANOVA independente unifatorial sem reamostragem funcionem independentemente da distribuição da VD em cada condição independente.
ANOVA relacionada unifatorial
ANOVA relacionada unifatorial é utilizada quando pelo menos um participante é submetido a pelo menos duas condições experimentais, i.e., um delineamento intraparticipantes. Os dados desse delineamento são chamados de dados longitudinais.
O delineamento com medidas repetidas conduz a uma ANOVA com maior poder.
Neste delineamento intraparticipantes, cada participante pode ser controle de si mesmo. Nos participantes com pelo menos duas condições experimentais, sugere-se que esta ordem de submissão, se possível, seja aleatorizada ou contrabalanceada. Para o cálculo da estatística F, a variância total da VD é particionada entre as variâncias das condições dependentes e dos participantes.
Suposições
– as diferenças dos valores das VD são independentes entre as unidades observacionais. – as diferenças dos valores das VD têm distribuição normal multivariada. – esfericidade: homocedasticidade das variâncias das diferenças das VDs (explicações estão em Mauchly’s sphericity test: Wikipedia e Ferrari, Motta-Jr & Siqueira (2017).
“Vale a pena observar que a MANOVA é ainda um teste válido, mesmo com modestas violações na suposição de normalidade, particularmente quando os tamanhos dos grupos são iguais e existe um número razoável de participantes em cada grupo; por “razoável” entendemos que [para um delineamento] completamente intraparticipantes, [deve haver] pelo menos 22 participantes ao todo.”
Dancey & Reidy, 2019, p. 472
Exemplo: SNAP-Ed intraparticipantes com 3 condições
Usaremos os mesmos dados do exemplo anterior, mas imaginando que as 20 medidas de ingestão de sódio sejam do mesmo participante, submetido aos três diferentes programas educacionais. Com 20 unidades observacionais no estudo, a normalidade multivariada das três diferenças não pode ser automaticamente assumida.
Brendon, McGuirk e Melissa fazem com que seus alunos do SNAP-Ed mantenham diários do que comem por uma semana e depois calculem a ingestão diária de sódio em miligramas.
Estudantes atenderam os diferentes programas de educação nutricional sucessivamente, cada um deles em uma ordem aleatorizada previamente. Os instrutores querem ver se a ingestão média de sódio é a mesma quando cada um dos três programas foi seguido.
Hipóteses nula e alternativa
O delineamento do estudo é diferente, mas as hipóteses são as mesmas:
\[ \begin{cases} H_0: &\mu_{\text{Brendon}} = \mu_{\text{McGuirk}} = \mu_{\text{Robins}}\\ H_1: &\exists\,\mu_i \ne \mu_j,\quad i \ne j,\quad i,j=1,2,3 \end{cases}\\ \alpha=5\% \]
ANOVA relacionada unifatorial
O arquivo de dados no formato long é Nutricao3rm.xlsx. Verifique a coluna
Student, que foi alterada para indicar o participante nos
três programas (compare com Nutricao3.xlsx, usada para ANOVA
independente unifatorial).
alfa <- 0.05
Dados <- data.frame(readxl::read_excel("Nutricao3rm.xlsx"))
Dados$Student <- factor(Dados$Student)
Dados$Instructor <- factor(Dados$Instructor,
labels=c("Brendon", "McGuirk", "Melissa"))
print(data.frame(Dados)) Instructor Student Sodium
1 Brendon a 1200
2 McGuirk a 1100
3 Melissa a 900
4 Brendon b 1400
5 McGuirk b 1200
6 Melissa b 1100
7 Brendon c 1350
8 McGuirk c 1250
9 Melissa c 1150
10 Brendon d 950
11 McGuirk d 1050
12 Melissa d 950
13 Brendon e 1400
14 McGuirk e 1200
15 Melissa e 1100
16 Brendon f 1150
17 McGuirk f 1250
18 Melissa f 1150
19 Brendon g 1300
20 McGuirk g 1350
21 Melissa g 1250
22 Brendon h 1325
23 McGuirk h 1350
24 Melissa h 1250
25 Brendon i 1425
26 McGuirk i 1325
27 Melissa i 1225
28 Brendon j 1500
29 McGuirk j 1525
30 Melissa j 1325
31 Brendon k 1250
32 McGuirk k 1225
33 Melissa k 1125
34 Brendon l 1150
35 McGuirk l 1125
36 Melissa l 1025
37 Brendon m 950
38 McGuirk m 1000
39 Melissa m 950
40 Brendon n 1150
41 McGuirk n 1125
42 Melissa n 925
43 Brendon o 1600
44 McGuirk o 1400
45 Melissa o 1200
46 Brendon p 1300
47 McGuirk p 1200
48 Melissa p 1100
49 Brendon q 1050
50 McGuirk q 1150
51 Melissa q 950
52 Brendon r 1300
53 McGuirk r 1400
54 Melissa r 1300
55 Brendon s 1700
56 McGuirk s 1500
57 Melissa s 1400
58 Brendon t 1300
59 McGuirk t 1200
60 Melissa t 1100O arquivo de dados está no formato long: cada linha do arquivo contém uma observação da medida (measure) de cada unidade experimental. Portanto, o número de linhas do arquivo long é \(r = \sum_{i=1}^{n}{q_i}\), sendo \(n\) o número de unidades experimentais e \(q_i\) o número de condições dependentes nas quais a unidade experimental \(i\) foi observada na medida (measure).
No caso particular desse estudo, com arquivo de dados perfeitamente balanceado, \(r=20\times 3=60\).
Sob a hipótese nula \(H_0: \mu_1 = \mu_2 = \cdots = \mu_g\), a estatística de teste F, segue uma distribuição F de Fisher-Snedecor com graus de liberdade \(df_A=q-1\) e \(df_E=(r-1)-((n-1)+(q-1))\) no numerador e denominador, respectivamente, sendo que \(q\) = número de condições dependentes, \(n\) é o número de unidades experimentais e \(r = \sum_{i=1}^{n}{q_i}\) é o número total de observações da variável dependente (measure) ou número de linhas do arquivo long.
\[ \begin{align} F &\underset{H_0}{\sim} F_{df_A, \, df_E}\\ F &\underset{H_0}{\sim} F_{q-1, \, (\sum_{i=1}^{n}{q_i}-1)-((n-1)+(q-1))}\\ F &\underset{H_0}{\sim} F_{q-1, \, (r-1)-((n-1)+(q-1))}\\ F &\underset{H_0}{\sim} F_{q-1, \, r-n-q+1} \end{align} \]
Se o número de unidades experimentais é muito maior do que o número de condições dependentes, i.e. (\(n\gg\bar{q}\)), sendo que \(\bar{q}=\frac{1}{n}\sum_{i=1}^{n}{q_i}\) (dados desbalanceados), \(df_E\approx n(\bar{q}-1)\). Se os dados são perfeitamente balanceados, \(df_E \approx n(q-1)\).
Se os dados são completos (perfeitamente balanceados):
Modelo Linear Geral (GLM): lm: * \(df_E = (n-1)(q-1) = (nq-1)-((n-1)+(q-1))\)
exatamente. No estudo do exemplo, \(df_E =
(20-1)(3-1) = (20\times 3-1)-((20-1)+(3-1))=38\).
Modelo Linear Geral Misto (GLMM):lmerTest::lmer: *\(df_E\) é aproximado por Satterthwaite ou
Kenward–Roger. \(df_E =
(n-1)(q-1)(nq-1)-((n-1)+(q-1))\) exatamente. No estudo do
exemplo, \(df_E = (20-1)(3-1) = (20\times
3-1)-((20-1)+(3-1))=38\).
\[ \begin{align} F &\sim F_{q-1, \; (nq-1)-((n-1)+(q-1))}\\ F &\sim F_{q-1, \; (n-1)(q-1)} \end{align} \]
Interpretação de \(df_2\):
No formato long, use \(r\) para o total de observações (número de linhas).
Decomposição geral do número de graus de liberdade de \(df_2\) no modelo
y ~ Condition + Subject:
- \(df_{\text{total}} = r − 1\)
- \(df_{\text{subject}} = n − 1\)
- \(df_{\text{condition}} = q − 1\)
- \(df_{\text{residual}} = r − n − q + 1\)
Observe que:
\[ df_{\text{residual}}= df_{\text{total}} - (df_{\text{subject}} + df_{\text{condition}}) \]
Caso balanceado (cada unidade experimental em todas as \(q\) condições dependentes), \(r = n q\) e então:
\[ \begin{align} df_{\text{residual}} &= (r - 1) − (n-1) - (q - 1) \\ &= (nq - 1) − (n-1) - (q - 1) \\ df_{\text{residual}} &=(n − 1)(q − 1) \end{align} \]
modelo <-
lmerTest::lmer(Sodium ~ Instructor + (1|Student), data=Dados)
car::Anova(modelo, test.statistic=“F”)
modelo <-
lm(Sodium ~ Instructor + Student, data=Dados)
car::Anova(modelo)
- GLM e GLMM univariados omnibus
Student
Instructor a b c d e f g h i j k l m n o p q r s t
Brendon 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
McGuirk 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Melissa 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
GLMAnova Table (Type II tests)
Response: Sodium
Sum Sq Df F value Pr(>F)
Instructor 290146 2 30.58 1.22e-08 ***
Student 1307542 19 14.51 5.37e-12 ***
Residuals 180271 38
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Response: Sodium
Df Sum Sq Mean Sq F value Pr(>F)
Instructor 2 290146 145073 30.58 1.22e-08 ***
Student 19 1307542 68818 14.51 5.37e-12 ***
Residuals 38 180271 4744
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GLMMmodelo <- lmerTest::lmer(Sodium ~ Instructor + (1|Student),
data=Dados)
print(anv <- car::Anova(modelo,
test.statistic="F"), digits=4)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Sodium
F Df Df.res Pr(>F)
Instructor 30.58 2 38 1.22e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Instructor 290146 145073 2 38 30.58 1.22e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Instructor 290146 145073 2 38 30.58 1.22e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- GLMM omnibus e post hoc
Student
Instructor a b c d e f g h i j k l m n o p q r s t
Brendon 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
McGuirk 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Melissa 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Instructor Sodium n_obs.x SodiumNormed sd n_obs.y se ci
1 Brendon 1288 20 1288 91.5 20 20.5 53.7
2 McGuirk 1246 20 1246 47.6 20 10.6 28.0
3 Melissa 1124 20 1124 59.9 20 13.4 35.2
GLMM
ANOVA
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Sodium
F Df Df.res Pr(>F)
Instructor 30.581 2 38 1.217e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RegressionLinear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Sodium ~ Instructor + (1 | Student)
Data: Dados
REML criterion at convergence: 704
Scaled residuals:
Min 1Q Median 3Q Max
-1.71211 -0.44338 -0.05566 0.54555 2.09262
Random effects:
Groups Name Variance Std.Dev.
Student (Intercept) 21358 146.14
Residual 4744 68.88
Number of obs: 60, groups: Student, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1287.50 36.13 24.37 35.639 < 2e-16 ***
InstructorMcGuirk -41.25 21.78 38.00 -1.894 0.0659 .
InstructorMelissa -163.75 21.78 38.00 -7.518 4.95e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Effect size analysis# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
----------------------------------------------------------
Instructor | 0.6168 | [0.4041, 0.7388] | large
Post hoc tests Instructor emmean SE df lower.CL upper.CL
Brendon 1288 36.1 24.4 1213 1362
McGuirk 1246 36.1 24.4 1172 1321
Melissa 1124 36.1 24.4 1049 1198
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
contrast estimate SE df lower.CL upper.CL t.ratio p.value
Brendon - McGuirk 41.2 21.8 38 -13.3 95.8 1.894 0.0659
Brendon - Melissa 163.8 21.8 38 109.2 218.3 7.518 <.0001
McGuirk - Melissa 122.5 21.8 38 67.9 177.1 5.624 <.0001
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates
P value adjustment: holm method for 3 tests 
Instructor emmean SE df lower.CL upper.CL .group
Melissa 1124 36.1 24.4 1031 1217 a
McGuirk 1246 36.1 24.4 1153 1339 b
Brendon 1288 36.1 24.4 1195 1380 b
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates
P value adjustment: holm method for 3 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. ANOVA relacionada unifatorial perfeitamente balanceada por bootstrapping
As funções WRS2::rmanovab e WRS2::pairdepb
realizam a ANOVA relacionada unifatorial perfeitamente balanceada por
reamostragem (demo_WRS2_ANOVA_dependente.R).
WRS2::rmanovab(Dados$Sodium, Dados$Instructor, tr=0, Dados$Student, nboot=B)
WRS2::pairdepb(Dados$Sodium, Dados$Instructor, tr=0, Dados$Student, nboot=B)
Student
Instructor a b c d e f g h i j k l m n o p q r s t
Brendon 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
McGuirk 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Melissa 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
"WRS" bootstrapped repeated measure one-Way ANOVACall:
WRS2::rmanovab(y = Dados$Sodium, groups = Dados$Instructor, blocks = Dados$Student,
tr = 0, nboot = B)
Test statistic: 30.5805
Critical value: 4.0915
Significant: TRUE
"WRS" post hocCall:
WRS2::pairdepb(y = Dados$Sodium, groups = Dados$Instructor, blocks = Dados$Student,
tr = 0, nboot = B)
psihat ci.lower ci.upper test crit sig
Brendon vs. McGuirk 41.25 -35.92592 118.4259 1.75616 2.71067 FALSE
Brendon vs. Melissa 163.75 82.40419 245.0958 5.41496 2.71067 TRUE
McGuirk vs. Melissa 122.50 85.28899 159.7110 8.19517 2.71067 TRUEAs conclusões permanecem as mesmas.
|
Na linguagem do modelo linear geral (GLM), na ANOVA relacionada unifatorial precisamos distinguir medida de variável dependente ou de desfecho: há apenas uma medida (measure) que, em nosso exemplo, é a quantidade de sódio ingerida. As variáveis dependentes (VD) são estas observações da medida em cada condição experimental (os programas de cada instrutor). Adicionalmente, na execução do teste, buscamos mitigar a dependência das três observações em cada indivíduo, o que é feito pelas diferenças, par-a-par, das observações entre todas as condições às quais cada indivíduo foi submetido. A rigor, as suposições de multinormalidade e de esfericidade poderiam ser feitas sobre estas diferenças. Note que, no caso da ANOVA independente unifatorial, medida e VD são idênticas: há apenas uma medida, que é a própria VD. |
ANOVA relacionada unifatorial desbalanceada
Vamos usar os mesmos dados, mas eliminar uma única observação
(removemos uma observação do estudante a de Melissa ,
planilha Nutricao3rmdesb.xlsx).
alfa <- 0.05
Dados <- data.frame(readxl::read_excel("Nutricao3rmdesb.xlsx"))
Dados$Student <- factor(Dados$Student)
Dados$Instructor <- factor(Dados$Instructor)
Dados$Instructor <- factor(Dados$Instructor,
labels=c("Brendon", "McGuirk", "Melissa"))
saveRDS(Dados, "Nutricao3rmdesb.rds")Para lidar com esta situação, utilizamos o mesmo código R utilizado
para a versão com dados perfeitamente balanceados, implementado em demo_ANOVA1f_dep_desbalanc_sodio.R. A
função lmerTest::lmer é capaz de lidar com dados
longitudinais.
modelo <-
lmerTest::lmer(Sodium ~ Instructor + (1|Student), data=Dados)
car::Anova(modelo, test.statistic=“F”)
modelo <-
lm(Sodium ~ Instructor + Student, data=Dados)
car::Anova(modelo)
- GLM e GLMM omnibus
Student
Instructor a b c d e f g h i j k l m n o p q r s t
Brendon 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
McGuirk 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Melissa 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
GLMAnova Table (Type II tests)
Response: Sodium
Sum Sq Df F value Pr(>F)
Instructor 290146 2 30.58 1.22e-08 ***
Student 1307542 19 14.51 5.37e-12 ***
Residuals 180271 38
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Response: Sodium
Df Sum Sq Mean Sq F value Pr(>F)
Instructor 2 290146 145073 30.58 1.22e-08 ***
Student 19 1307542 68818 14.51 5.37e-12 ***
Residuals 38 180271 4744
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GLMMmodelo <- lmerTest::lmer(Sodium ~ Instructor + (1|Student),
data=Dados)
print(anv <- car::Anova(modelo,
test.statistic="F"), digits=4)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Sodium
F Df Df.res Pr(>F)
Instructor 30.58 2 38 1.22e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Instructor 290146 145073 2 38 30.58 1.22e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Instructor 290146 145073 2 38 30.58 1.22e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- GLMM omnibus e post hoc
Student
Instructor a b c d e f g h i j k l m n o p q r s t
Brendon 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
McGuirk 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Melissa 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Instructor Sodium n_obs.x SodiumNormed sd n_obs.y se ci
1 Brendon 1288 20 1289 89.7 20 20.0 52.6
2 McGuirk 1246 20 1247 52.0 20 11.6 30.5
3 Melissa 1136 19 1133 57.8 19 13.3 35.0
GLMM
ANOVA
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Sodium
F Df Df.res Pr(>F)
Instructor 27.356 2 37.061 5.043e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RegressionLinear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Sodium ~ Instructor + (1 | Student)
Data: Dados
REML criterion at convergence: 691.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.7013 -0.4632 -0.0990 0.5474 2.1420
Random effects:
Groups Name Variance Std.Dev.
Student (Intercept) 20846 144.38
Residual 4654 68.22
Number of obs: 59, groups: Student, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1287.50 35.71 24.34 36.057 < 2e-16 ***
InstructorMcGuirk -41.25 21.57 37.03 -1.912 0.0636 .
InstructorMelissa -157.51 21.98 37.12 -7.166 1.68e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Effect size analysis# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
----------------------------------------------------------
Instructor | 0.5962 | [0.3729, 0.7253] | large
Post hoc tests Instructor emmean SE df lower.CL upper.CL
Brendon 1288 35.7 24.3 1214 1361
McGuirk 1246 35.7 24.3 1173 1320
Melissa 1130 36.0 24.9 1056 1204
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
contrast estimate SE df lower.CL upper.CL t.ratio p.value
Brendon - McGuirk 41.2 21.6 37.0 -12.8 95.3 1.912 0.0636
Brendon - Melissa 157.5 22.0 37.1 102.4 212.6 7.166 <.0001
McGuirk - Melissa 116.3 22.0 37.1 61.1 171.4 5.289 <.0001
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates
P value adjustment: holm method for 3 tests 
Instructor emmean SE df lower.CL upper.CL .group
Melissa 1130 36.0 24.9 1038 1222 a
McGuirk 1246 35.7 24.3 1154 1338 b
Brendon 1288 35.7 24.3 1196 1379 b
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates
P value adjustment: holm method for 3 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
|
A ANOVA relacionada unifatorial desbalanceada feita com o pacote
|
Exemplo de relatório para a ANOVA relacionada unifatorial
Vinte participantes foram selecionados para um estudo com
delineamento intrapartipantes com três condições experimentais. A ordem
de aplicação das três condições experimentais em cada participante foi
aleatorizada. As médias amostrais brutas mostram que menor ingestão
diária média de sódio foi observada entre estudantes submetidos ao
programa aplicado por Melissa Robins. A ingestão diária média de sódio
dos estudantes de Brendon e McGuirk são semelhantes. A ANOVA relacionada
unifatorial por meio do modelo linear misto geral
(lmerTest::lmer) rejeitou a hipótese nula omnibus,
pois a estatística de teste observada é F(2,37.1) = 27.369
e o valor p associado é igual a \(1.22 \times 10^{-8}\). O tamanho do efeito
do fator intraparticipantes Instructor é estimado pelo eta
ao quadrado de Cohen cujo valor indica que 62% da variância da ingesta
de sódio é explicada pelo efeito do fator fixo Instructor.
A análise post-hoc confirmou que as diferenças das médias populacionais
entre entre os programas adotados por Melissa e Brendon, e entre os
adotados por Melissa e McGuirk são estatisticamente significantes. Não
se observou diferença estatisticamente significante entre os programas
adotados por Brendon e McGuirk.
|
Para aprofundar em métodos para resolver ANOVA relacionada, há oito possibilidades discutidas em Keselman et al. (2001). |
Análise de curva de crescimento
- Johnson & Wichern (2007)
Quando a estatura de uma criança pequena é medida a cada aniversário, os pontos podem ser plotados e conectados por linhas para produzir um gráfico. Isso é um exemplo de uma curva de crescimento. Em geral, medições repetidas da mesma característica na mesma unidade observacional ou participante podem dar origem a uma curva de crescimento se for esperado um padrão crescente, decrescente ou até mesmo crescente seguido de decrescente.
Exemplo: Curva de crescimento de ursa do Alasca
O Departamento de Pesca e Caça do Alasca monitora urso-cinzento ou
pardo (Ursus arctos horribilis) com o objetivo de
manter uma população saudável. Os ursos são anestesiados com um dardo
para induzir o sono e pesados em uma balança suspensa em um tripé.
Medidas de comprimento são feitas com uma fita métrica de aço. A Tabela
1.4 fornece as massas corporais totais (Wt) em quilograma e
os comprimentos (Length) em centímetro de sete ursos fêmeas
nas idades de 2, 3, 4 e 5 anos.
Tabela 1.4 Dados de Ursas
| Ursa | Wt2 | Wt3 | Wt4 | Wt5 | Lngth2 | Lngth3 | Lngth4 | Lngth5 |
|---|---|---|---|---|---|---|---|---|
| 1 | 48 | 59 | 95 | 82 | 141 | 157 | 168 | 183 |
| 2 | 59 | 68 | 102 | 102 | 140 | 168 | 174 | 170 |
| 3 | 61 | 77 | 93 | 107 | 145 | 162 | 172 | 177 |
| 4 | 54 | 43 | 104 | 104 | 146 | 159 | 176 | 171 |
| 5 | 45 | 66 | 84 | 112 | 150 | 158 | 168 | 175 |
| 6 | 68 | 82 | 95 | 118 | 142 | 140 | 178 | 189 |
| 7 | 68 | 95 | 109 | 111 | 139 | 171 | 176 | 175 |
Primeiro, para cada ursa, plotamos os pesos em relação às idades e, em seguida, conectamos os pesos nos anos sucessivos por linhas retas. Isso dá uma aproximação da curva de crescimento para a massa corporal total.
No campo, os ursos são pesados em uma balança que indica libras.
Por ser difícil inspecionar visualmente as curvas de crescimento individuais em um gráfico combinado, as curvas individuais devem ser plotadas novamente em um arranjo onde semelhanças e diferenças sejam facilmente observadas. Algumas curvas de crescimento parecem lineares, enquanto outras são não lineares.
A ursa 6 parece ter ficado menor dos 2 aos 3 anos de idade, mas o pesquisador sabe que a medição do comprimento com uma fita de aço pode ser afetada pela postura da ursa quando sedada.
Pessoa (2008) usa a estrutura de dados de série temporal (dados equiespaçados) para plotar a curva de crescimento de cada ursa. Este não é método geral, pois os dados da curva de crescimento podem ser não equiespaçados.
Os intervalos de confiança de Mirnam (2014) são incorretos, pois são intervalos para delineamento entre participantes e sem correção de Bonferroni. Os intervalos de confiança têm que ser adequados ao delineamento intraparticipantes (i.e., com medidas repetidas).
source("summarySEwithin2.R")
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
# Growth Curve Analysis and Visualization Using R - Mirman - 2014
Year <- rep(2:5, each = 7)
Bear <- rep(1:7, 4)
Wt2 <- c(48, 59, 61, 54, 45, 68, 68)
Wt3 <- c(59, 68, 77, 43, 66, 82, 95)
Wt4 <- c(95, 102, 93, 104, 84, 95, 109)
Wt5 <- c(82, 102, 107, 104, 112, 118, 111)
Wt <- c(Wt2, Wt3, Wt4, Wt5)
Lngth2 <- c(141, 140, 145, 146, 150, 142, 139)
Lngth3 <- c(157, 168, 162, 159, 158, 140, 171)
Lngth4 <- c(168, 174, 172, 176, 168, 178, 176)
Lngth5 <- c(183, 170, 177, 171, 175, 189, 175)
Lngth <- c(Lngth2, Lngth3, Lngth4, Lngth5)
dados_ursas <- data.frame(Year, Bear, Wt, Lngth)
dados_ursas$Bear <- factor(dados_ursas$Bear)
print(car::some(dados_ursas)) Year Bear Wt Lngth
2 2 2 59 140
3 2 3 61 145
5 2 5 45 150
6 2 6 68 142
15 4 1 95 168
16 4 2 102 174
18 4 4 104 176
20 4 6 95 178
27 5 6 118 189
28 5 7 111 175 Bear
Year 1 2 3 4 5 6 7
2 1 1 1 1 1 1 1
3 1 1 1 1 1 1 1
4 1 1 1 1 1 1 1
5 1 1 1 1 1 1 1#knitr::kable(dados_ursas, caption="Dados de Ursas")
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Wt, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10) +
ggplot2::facet_wrap(~ Bear)
ggplot2::ggsave("BearWtSep.png", width=3, height=3, dpi=300)
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Lngth, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10) +
ggplot2::facet_wrap(~ Bear)
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Wt, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10)
ggplot2::ggplot(dados_ursas,
ggplot2::aes(Year, Lngth, linetype=Bear)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::theme_bw(base_size = 10)
alfa <- 0.05
q <- nlevels(factor(dados_ursas$Year))
alfaBonf <- alfa/q
ic <- summarySEwithin2(dados_ursas,
measurevar="Wt",
withinvars="Year",
idvar="Bear",
na.rm=TRUE,
conf.interval=1-alfaBonf)Automatically converting to factor: Year Year Wt n_obs.x WtNormed sd n_obs.y se ci
1 2 57.6 7 57.6 4.04 7 1.53 5.38
2 3 70.0 7 70.0 11.70 7 4.42 15.57
3 4 97.4 7 97.4 10.35 7 3.91 13.78
4 5 105.1 7 105.1 9.52 7 3.60 12.67grf <- ggplot2::ggplot(ic,
ggplot2::aes(x=Year,
y=Wt)) +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=Wt-ci,
ymax=Wt+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylab("Weight") +
ggplot2::ggtitle("7 Female Bear: Weight\nWithin-subject CI95 Bonferroni") +
ggplot2::theme_bw()
print(grf, digits=3)
ic <- summarySEwithin2(dados_ursas,
measurevar="Lngth",
withinvars="Year",
idvar="Bear",
na.rm=TRUE,
conf.interval=
1-alfaBonf)Automatically converting to factor: Year Year Lngth n_obs.x LngthNormed sd n_obs.y se ci
1 2 143 7 143 4.99 7 1.89 6.65
2 3 159 7 159 10.70 7 4.04 14.24
3 4 173 7 173 4.43 7 1.67 5.89
4 5 177 7 177 8.37 7 3.16 11.14grf <- ggplot2::ggplot(ic,
ggplot2::aes(x=Year,
y=Lngth)) +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=Lngth-ci,
ymax=Lngth+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylab("Length") +
ggplot2::ggtitle("7 Female Bear: Length\nWithin-subject CI95 Bonferroni") +
ggplot2::theme_bw()
print(grf)
GLMM: Weight
RegressionLinear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Wt ~ factor(Year) + (1 | Bear)
Data: dados_ursas
REML criterion at convergence: 190.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.4063 -0.3574 0.1039 0.5184 1.6687
Random effects:
Groups Name Variance Std.Dev.
Bear (Intercept) 53.40 7.307
Residual 87.78 9.369
Number of obs: 28, groups: Bear, 7
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 57.571 4.491 16.793 12.820 4.26e-10 ***
factor(Year)3 12.429 5.008 18.000 2.482 0.0232 *
factor(Year)4 39.857 5.008 18.000 7.959 2.64e-07 ***
factor(Year)5 47.571 5.008 18.000 9.499 1.96e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVAAnalysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Wt
F Df Df.res Pr(>F)
factor(Year) 40.224 3 18 3.464e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
factor(Year) 10593 3531 3 18 40.22 3.46e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
factor(Year) 10593 3531 3 18 40.22 3.46e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GLMAnova Table (Type II tests)
Response: Wt
Sum Sq Df F value Pr(>F)
factor(Year) 10593 3 40.224 3.46e-08 ***
Bear 1808 6 3.433 0.0194 *
Residuals 1580 18
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Response: Wt
Df Sum Sq Mean Sq F value Pr(>F)
factor(Year) 3 10593 3531 40.224 3.46e-08 ***
Bear 6 1808 301 3.433 0.0194 *
Residuals 18 1580 88
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Effect size analysiseta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alfa,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=4)# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
------------------------------------------------------------
factor(Year) | 0.8702 | [0.7167, 0.9254] | large
Post hoc testsTeste post hoc: Pairwiseemm <- emmeans::emmeans(modelo,
specs=pairwise~Year,
adjust="holm",
level=1-alfa,
lmer.df="satterthwaite",
lmerTest.limit=nrow(dados_ursas))
print(summary(emm$emmeans)) Year emmean SE df lower.CL upper.CL
2 57.6 4.49 16.8 48.1 67.1
3 70.0 4.49 16.8 60.5 79.5
4 97.4 4.49 16.8 87.9 106.9
5 105.1 4.49 16.8 95.7 114.6
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95 contrast estimate SE df lower.CL upper.CL t.ratio p.value
Year2 - Year3 -12.43 5.01 18 -27.3 2.41 -2.482 0.0463
Year2 - Year4 -39.86 5.01 18 -54.7 -25.02 -7.959 <.0001
Year2 - Year5 -47.57 5.01 18 -62.4 -32.73 -9.499 <.0001
Year3 - Year4 -27.43 5.01 18 -42.3 -12.59 -5.477 0.0001
Year3 - Year5 -35.14 5.01 18 -50.0 -20.31 -7.017 <.0001
Year4 - Year5 -7.71 5.01 18 -22.6 7.12 -1.540 0.1409
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 6 estimates
P value adjustment: holm method for 6 tests 
Year emmean SE df lower.CL upper.CL .group
2 57.6 4.49 16.8 45.0 70.1 a
3 70.0 4.49 16.8 57.4 82.6 b
4 97.4 4.49 16.8 84.9 110.0 c
5 105.1 4.49 16.8 92.6 117.7 c
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 4 estimates
P value adjustment: holm method for 6 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Teste post hoc: DunnettEMM.contrast <- emmeans::contrast(object=emm,
# by="Grupo",
method="trt.vs.ctrl",
ref="Year2",
adjust="holm",
level=1-alfa)
print(EMM.contrast) contrast estimate SE df t.ratio p.value
Year3 - Year2 12.4 5.01 18 2.482 0.0232
Year4 - Year2 39.9 5.01 18 7.959 <.0001
Year5 - Year2 47.6 5.01 18 9.499 <.0001
Degrees-of-freedom method: satterthwaite
P value adjustment: holm method for 3 tests 
Teste post hoc: ConsecutivoEMM.contrast <- emmeans::contrast(object=emm,
# by="Grupo",
method = "consec",
adjust="holm",
level=1-alfa)
print(EMM.contrast) contrast estimate SE df t.ratio p.value
Year3 - Year2 12.43 5.01 18 2.482 0.0463
Year4 - Year3 27.43 5.01 18 5.477 0.0001
Year5 - Year4 7.71 5.01 18 1.540 0.1409
Degrees-of-freedom method: satterthwaite
P value adjustment: holm method for 3 tests 
Teste post hoc: curvas polinomiais de graus 1, 2 e 3EMM.contrast <- emmeans::contrast(object=emm,
#by="Grupo",
method="poly",
adjust="holm",
level=1-alfa)
print(EMM.contrast) contrast estimate SE df t.ratio p.value
linear 170.14 15.80 18 10.743 <.0001
quadratic -4.71 7.08 18 -0.666 0.5141
cubic -34.71 15.80 18 -2.192 0.0835
Degrees-of-freedom method: satterthwaite
P value adjustment: holm method for 3 tests 
GLMM: Lengthboundary (singular) fit: see help('isSingular')
RegressionLinear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Lngth ~ factor(Year) + (1 | Bear)
Data: dados_ursas
REML criterion at convergence: 166.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.9039 -0.3818 -0.1076 0.4302 1.7854
Random effects:
Groups Name Variance Std.Dev.
Bear (Intercept) 0.00 0.000
Residual 44.11 6.641
Number of obs: 28, groups: Bear, 7
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 143.29 2.51 24.00 57.082 < 2e-16 ***
factor(Year)3 16.00 3.55 24.00 4.507 0.000145 ***
factor(Year)4 29.86 3.55 24.00 8.411 1.29e-08 ***
factor(Year)5 33.86 3.55 24.00 9.537 1.23e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
ANOVAAnalysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Lngth
F Df Df.res Pr(>F)
factor(Year) 37.304 3 18 6.204e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
factor(Year) 4936 1645 3 18 37.3 6.2e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
factor(Year) 4936 1645 3 24 37.3 3.39e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GLMmodeloGLM <- lm(Lngth ~ factor(Year) + Bear,
data=dados_ursas)
print(car::Anova(modeloGLM), digits=4)Anova Table (Type II tests)
Response: Lngth
Sum Sq Df F value Pr(>F)
factor(Year) 4936 3 28.730 4.44e-07 ***
Bear 28 6 0.081 0.997
Residuals 1031 18
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Response: Lngth
Df Sum Sq Mean Sq F value Pr(>F)
factor(Year) 3 4936 1645.4 28.730 4.44e-07 ***
Bear 6 28 4.6 0.081 0.997
Residuals 18 1031 57.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Effect size analysiseta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alfa,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=4)# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
------------------------------------------------------------
factor(Year) | 0.8614 | [0.6980, 0.9203] | large
Post hoc testsTeste post hoc: Pairwiseemm <- emmeans::emmeans(modelo,
specs=pairwise~Year,
adjust="holm",
level=1-alfa,
lmer.df="satterthwaite",
lmerTest.limit=nrow(dados_ursas))
print(summary(emm$emmeans)) Year emmean SE df lower.CL upper.CL
2 143 2.51 24 138 148
3 159 2.51 24 154 164
4 173 2.51 24 168 178
5 177 2.51 24 172 182
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95 contrast estimate SE df lower.CL upper.CL t.ratio p.value
Year2 - Year3 -16.0 3.55 24 -26.2 -5.79 -4.507 0.0004
Year2 - Year4 -29.9 3.55 24 -40.1 -19.65 -8.411 <.0001
Year2 - Year5 -33.9 3.55 24 -44.1 -23.65 -9.537 <.0001
Year3 - Year4 -13.9 3.55 24 -24.1 -3.65 -3.903 0.0013
Year3 - Year5 -17.9 3.55 24 -28.1 -7.65 -5.030 0.0002
Year4 - Year5 -4.0 3.55 24 -14.2 6.21 -1.127 0.2710
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 6 estimates
P value adjustment: holm method for 6 tests 
Year emmean SE df lower.CL upper.CL .group
2 143 2.51 24 137 150 a
3 159 2.51 24 153 166 b
4 173 2.51 24 166 180 c
5 177 2.51 24 170 184 c
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 4 estimates
P value adjustment: holm method for 6 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Teste post hoc: DunnettEMM.contrast <- emmeans::contrast(object=emm,
# by="Grupo",
method="trt.vs.ctrl",
ref="Year2",
adjust="holm",
level=1-alfa)
print(EMM.contrast) contrast estimate SE df t.ratio p.value
Year3 - Year2 16.0 3.55 24 4.507 0.0001
Year4 - Year2 29.9 3.55 24 8.411 <.0001
Year5 - Year2 33.9 3.55 24 9.537 <.0001
Degrees-of-freedom method: satterthwaite
P value adjustment: holm method for 3 tests 
Teste post hoc: ConsecutivoEMM.contrast <- emmeans::contrast(object=emm,
# by="Grupo",
method = "consec",
adjust="holm",
level=1-alfa)
print(EMM.contrast) contrast estimate SE df t.ratio p.value
Year3 - Year2 16.0 3.55 24 4.507 0.0004
Year4 - Year3 13.9 3.55 24 3.903 0.0013
Year5 - Year4 4.0 3.55 24 1.127 0.2710
Degrees-of-freedom method: satterthwaite
P value adjustment: holm method for 3 tests 
Teste post hoc: curvas polinomiais de graus 1, 2 e 3EMM.contrast <- emmeans::contrast(object=emm,
#by="Grupo",
method="poly",
adjust="holm",
level=1-alfa)
print(EMM.contrast) contrast estimate SE df t.ratio p.value
linear 115.43 11.20 24 10.282 <.0001
quadratic -12.00 5.02 24 -2.390 0.0501
cubic -7.71 11.20 24 -0.687 0.4986
Degrees-of-freedom method: satterthwaite
P value adjustment: holm method for 3 tests 
GLM multivariadomodelo <- lm(cbind(Wt,Lngth) ~ factor(Year) + Bear,
data=dados_ursas)
print(car::Anova(modelo), digits=4)
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
factor(Year) 3 1.0314 6.390 6 36 0.00012 ***
Bear 6 0.5525 1.145 12 36 0.35715
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Df Pillai approx F num Df den Df Pr(>F)
(Intercept) 1 0.9988 7262 2 17 < 2e-16 ***
factor(Year) 3 1.0314 6 6 36 0.00012 ***
Bear 6 0.5525 1 12 36 0.35715
Residuals 18
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Referências
- Cohen, J (1992) A power primer. Quantitative methods in Psychology 112(1): 155-159.
- Dancey C & Reidy J (2019) Estatística sem Matemática para Psicologia. 7a ed. Porto Alegre: PENSO.
- Ferrari, A, Motta-Junior, JC & Siqueira, JO (2017) Métodos de amostragem e análise em estudos sobre comportamento de forrageio de aves. Oecologia Australis 21(2): 191–206. https://doi.org/10.4257/oeco.2017.2102.03
- Games, PA & Howell, JF (1976) Pairwise Multiple Comparison Procedures with Unequal N’s and/or Variances: A Monte Carlo Study. Journal of Educational and Behavioral Statistics 1: 113–25. https://doi.org/10.3102/10769986001002113
- Johnson, R & Wichern, D (2007) Applied Multivariate Statistical Analysis. 6th ed. NJ: Prentice-Hall.
- Kenward, MG & Roger, JH (1997) Small sample inference for fixed effects from restricted maximum likelihood. Biometrics 53(3): 983-97.
- Keselman HJ, Algina J, Kowalchuk, R (2001) The analysis of repeated measures designs: A review. British Journal of Mathematical and Statistical Psychology 54(1).
- Moore, D (1995) The basic practice of statistics. NY: W. H. Freeman and Company.
- Norusis, M (2009) PASW Statistics 18: Statistics Procedures Companion. NJ: Prentice-Hall.
- Pestana, M & Gageiro, J (2008) Análise de dados para Ciências Sociais: a complementaridade do SPSS. Lisboa: Sílabo.
- Satterthwaite, FE (1946) An Approximate Distribution of Estimates of Variance Components. Biometrics Bulletin 2(6), 110-4.
- Welch, BL (1947) The generalization of “Student’s” problem when several different population variances are involved. Biometrika 34: 28-35.
- Welch, BL (1951) On the comparison of several mean values: An alternative approach. Biometrika 38(3-4): 330–6.