ANOVA Unifatorial em R
MSP4112: Estatística Aplicada à Medicina
05 outubro 2026 11:30h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(apaTables, warn.conflicts=FALSE))
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(Valores p incorretos)
- testes post hoc:
- Fisher sem dados individualizados (homocedástico)
- teste omnibus:
demo_ANOVA_Independente_Unifatorial_Fisher_SemDadosindividualizados.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(Valores p incorretos)
- testes post hoc:
- Welch sem dados individualizados
- 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(Valores p incorretos)
- 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
História
ANOVA
A análise de variância (ANOVA) foi desenvolvida por Ronald A. Fisher, principalmente durante seu trabalho na Rothamsted Experimental Station, na Inglaterra. Seu desenvolvimento ocorreu gradualmente entre o final da década de 1910 e a década de 1920.
1918 — Fisher introduziu o conceito de variância e ideias de decomposição da variabilidade em seu trabalho sobre genética mendeliana. Ainda não se tratava da ANOVA experimental como a conhecemos atualmente.
1923 — surgiu um dos primeiros exemplos de análise de variância aplicada a um experimento agrícola.
1924 — Fisher apresentou sistematicamente o método no International Mathematical Congress, em Toronto. O trabalho correspondente foi publicado em 1928.
1925 — a ANOVA foi apresentada de forma ampla no livro Statistical Methods for Research Workers, que teve grande influência na disseminação do método.
Fisher posteriormente integrou a ANOVA aos princípios de randomização, replicação e delineamento experimental, contribuindo decisivamente para a formação da estatística experimental moderna.
Ronald A. Fisher
Ronald A. Fisher não iniciou sua carreira como professor universitário. Após concluir sua formação em Matemática na University of Cambridge, trabalhou inicialmente como professor de matemática e física no ensino secundário.
Em 1919, aos 29 anos, ingressou na Rothamsted Experimental Station como estatístico. Foi nesse período, entre 1919 e 1933, que desenvolveu grande parte de suas contribuições fundamentais para a Estatística, incluindo a análise de variância, o delineamento de experimentos e importantes métodos de inferência estatística.
Sua carreira como professor universitário começou apenas em 1933, aos 43 anos, quando assumiu a posição de Galton Professor of Eugenics na University College London.
Em 1943, transferiu-se para a University of Cambridge, onde se tornou Arthur Balfour Professor of Genetics, permanecendo nessa posição até sua aposentadoria, em 1957.
Assim, Fisher já era um estatístico de grande reconhecimento científico antes de ocupar sua primeira cátedra universitária.
Cronologia:
- 1912 — formação em Matemática na University of Cambridge
- 1913–1919 — professor de matemática e física no ensino secundário
- 1919–1933 — estatístico na Rothamsted Experimental Station
- 1933–1943 — Galton Professor of Eugenics, University College London
- 1943–1957 — Arthur Balfour Professor of Genetics, University of Cambridge
- 1957 — aposentadoria
Ronald Fisher, fertilizantes e experimentação agrícola
Ronald A. Fisher ingressou na Rothamsted Experimental Station em 1919, aos 29 anos. Rothamsted já realizava experimentos agrícolas de longa duração desde o século XIX, incluindo estudos sobre os efeitos de diferentes fertilizantes na produtividade das culturas. O experimento de Broadbalk com trigo, por exemplo, havia começado em 1843.
O problema enfrentado pelos pesquisadores era distinguir o verdadeiro efeito de um tratamento da grande variabilidade natural existente entre parcelas agrícolas. Diferenças de solo, clima, umidade e outras condições podiam produzir diferenças de produtividade mesmo quando os tratamentos eram iguais.
Foi nesse contexto que Fisher desenvolveu e sistematizou métodos como análise de variância (ANOVA), randomização, replicação, blocos e delineamentos fatoriais. A partir de 1926, os novos princípios estatísticos desenvolvidos por Fisher passaram a ser utilizados de maneira sistemática nos experimentos anuais de Rothamsted.
Entre seus trabalhos estão experimentos destinados a quantificar a resposta das culturas a diferentes doses e tipos de fertilizantes. Em um estudo publicado com Tim Eden em 1929, por exemplo, foram analisados experimentos realizados entre 1925 e 1927 sobre a resposta da batata a diferentes quantidades de nitrogênio e potássio. O objetivo não era apenas verificar se o fertilizante funcionava, mas estimar de maneira confiável o efeito das diferentes doses e formulações.
Relação com a revolução dos fertilizantes nitrogenados
Esse desenvolvimento estatístico ocorreu justamente em um período de profunda transformação da agricultura. Em 1909, Fritz Haber demonstrou a síntese de amônia a partir do nitrogênio atmosférico e do hidrogênio. Carl Bosch tornou posteriormente o processo viável em escala industrial. O processo Haber–Bosch permitiu produzir grandes quantidades de amônia sintética, matéria-prima para fertilizantes nitrogenados.
Antes disso, a agricultura dependia principalmente do nitrogênio proveniente de esterco, leguminosas, depósitos naturais de nitrato e outras fontes relativamente limitadas. A produção industrial de amônia tornou possível fornecer nitrogênio às culturas em uma escala sem precedentes. Seu uso agrícola cresceu principalmente depois da Segunda Guerra Mundial e tornou-se um dos fatores centrais do grande aumento da produtividade agrícola no século XX.
É importante distinguir as contribuições:
- Haber e Bosch tornaram possível a produção industrial de grandes quantidades de fertilizantes nitrogenados
- Rothamsted forneceu o ambiente experimental para estudar quantitativamente seus efeitos sobre as culturas
- Fisher desenvolveu os métodos estatísticos que permitiram separar os efeitos reais dos tratamentos da variabilidade natural dos experimentos agrícolas
Assim, a expansão dos fertilizantes nitrogenados e o desenvolvimento da estatística experimental ocorreram aproximadamente no mesmo período e foram complementares: tornou-se possível produzir novos insumos agrícolas em grande escala e, simultaneamente, surgiu uma metodologia científica rigorosa para determinar se funcionavam, quanto funcionavam e em quais condições.
Eden & Fisher (1929) analisam experimentalmente a resposta da batata a diferentes aplicações de nitrogênio e potássio e documenta diretamente a atuação de Fisher em experimentos com fertilizantes em Rothamsted.
Modelo
As variáveis envolvidas na ANOVA são:
- VD (variável dependente) quantitativa, geralmente em escala intervalar ou de razão
- VI (variável independente) nominal com três ou mais níveis, denominada fator
Na ANOVA independente, assume-se que:
- a VD apresente distribuição aproximadamente normal em cada condição
- as variâncias da VD sejam homogêneas entre as condições
- as observações sejam independentes
Na ANOVA unifatorial (one-way ANOVA) existe apenas um fator.
A ANOVA unifatorial pode ser:
- independente: fator entre participantes
- relacionada: fator intraparticipantes
A ANOVA independente unifatorial (one-way analysis of variance) é uma extensão do teste t para grupos independentes quando existem três ou mais condições.
O teste \(F\) apresentado na tabela da ANOVA testa o efeito do fator fixo sobre a VD. Portanto, rigorosamente, não se testa o fator em si, mas o efeito associado ao fator.
Fatores fixos e aleatórios na ANOVA
Na terminologia tradicional da ANOVA, fala-se em fatores fixos e fatores aleatórios. Mais rigorosamente, porém, a distinção refere-se à maneira como os efeitos associados aos níveis do fator são tratados no modelo.
Portanto, o teste \(F\) apresentado na tabela da ANOVA não testa o fator como entidade abstrata, mas uma hipótese sobre os efeitos associados aos seus níveis.
Fator fixo
Um fator fixo possui níveis especificamente definidos pelo pesquisador, e a inferência é dirigida a esses níveis ou às diferenças entre seus efeitos.
O fator fixo pode ser:
- de interesse: quando seu efeito constitui o objetivo principal da análise
- de controle: quando é incluído para controlar uma fonte sistemática de variação e melhorar a estimação ou o teste dos efeitos de interesse
Para um fator fixo com \(a\) níveis, uma hipótese típica é
\[ H_0:\mu_1=\mu_2=\cdots=\mu_a \]
ou, equivalentemente, sob uma parametrização por efeitos,
\[ H_0:\alpha_1=\alpha_2=\cdots=\alpha_a=0 \]
O teste \(F\) avalia, portanto, uma hipótese sobre os efeitos do fator fixo sobre a variável dependente.
Fator aleatório
Um fator aleatório possui níveis considerados como uma amostra de uma população maior de níveis possíveis.
Nesse caso, o interesse não está primariamente nas diferenças entre os níveis específicos observados, mas na variabilidade dos efeitos associada à população de níveis do fator.
Se
\[ u_j\sim N(0,\sigma_u^2) \]
o parâmetro de interesse é, em geral, o componente de variância \(\sigma_u^2\), e não cada efeito \(u_j\) individualmente.
Exemplos usuais de fatores aleatórios incluem participantes, blocos, centros, avaliadores e unidades experimentais quando seus níveis representam uma população mais ampla de níveis possíveis.
Na ANOVA clássica, fatores aleatórios também são fundamentais para determinar os quadrados médios esperados e, consequentemente, o denominador apropriado de determinados testes \(F\).
Assim, de forma resumida:
- fator fixo → o interesse recai principalmente sobre os efeitos dos níveis especificados
- fator aleatório → o interesse recai principalmente sobre a variabilidade dos efeitos entre níveis
- teste \(F\) → compara quadrados médios construídos de acordo com a estrutura fixa, aleatória ou mista do delineamento
Portanto, é conveniente falar em fator fixo e fator aleatório, desde que se entenda que essa classificação descreve, em última análise, como os efeitos associados ao fator são modelados e sobre qual população se pretende fazer inferência.
Suposições
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:
- População muito grande, idealmente infinita, de modo que não ocorra deriva genética
- Acasalamento aleatório (panmixia)
- Ausência de mutação
- Ausência de migração, isto é, de fluxo gênico
- Ausência de 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: \(p=p_{\text{A}}\),
- selecionando para o alelo A.
cycle fA fB fAA fAB fBB
1 0.000 1.000 0.000 0.000 1.000
100 0.074 0.926 0.006 0.136 0.858
200 0.304 0.696 0.100 0.408 0.492
300 0.680 0.320 0.466 0.428 0.106
400 0.860 0.140 0.740 0.240 0.020
500 0.927 0.073 0.858 0.138 0.004
600 0.954 0.046 0.908 0.092 0.000
700 0.972 0.028 0.944 0.056 0.000
800 0.983 0.017 0.966 0.034 0.000
900 0.982 0.018 0.964 0.036 0.000
1000 0.985 0.015 0.970 0.030 0.000
1100 0.977 0.023 0.954 0.046 0.000
1200 0.967 0.033 0.934 0.066 0.000
1300 0.976 0.024 0.952 0.048 0.000
1400 0.980 0.020 0.960 0.040 0.000
1500 0.971 0.029 0.942 0.058 0.000
1600 0.984 0.016 0.968 0.032 0.000
1700 0.985 0.015 0.970 0.030 0.000
1800 0.986 0.014 0.972 0.028 0.000
1900 0.992 0.008 0.984 0.016 0.000
2000 0.992 0.008 0.984 0.016 0.000
2100 0.977 0.023 0.954 0.046 0.000
2200 0.987 0.013 0.974 0.026 0.000
2300 0.987 0.013 0.974 0.026 0.000
2400 0.987 0.013 0.974 0.026 0.000
2500 0.957 0.043 0.916 0.082 0.002
2600 0.975 0.025 0.950 0.050 0.000
2700 0.978 0.022 0.956 0.044 0.000
2800 0.975 0.025 0.950 0.050 0.000
2900 0.974 0.026 0.948 0.052 0.000
3000 0.977 0.023 0.956 0.042 0.002 
3000 0.977 0.023 0.956 0.042 0.002 
As condições matemáticas de população muito grande ou infinita, ausência de mutação e ausência de seleção contribuem para a manutenção do equilíbrio de Hardy–Weinberg. Entretanto, populações biológicas reais não satisfazem rigorosamente essas condições. Ainda assim, quando os desvios são suficientemente pequenos, as frequências genotípicas previstas pelo modelo podem constituir uma aproximação aceitável, permitindo sua aplicação a 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 conceito de variância
Para uma variável quantitativa intervalar, a variância pode ser interpretada como:
- o quadrado do desvio-padrão;
- a covariância da variável com ela mesma;
- uma medida da dispersão dos valores em torno da média.
Quando uma amostra é particionada em duas ou mais condições ou grupos, a variância total também pode ser interpretada como uma composição da variabilidade dentro dos grupos e da variabilidade entre os grupos.
- 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.
Suponha que não temos acesso aos dados individualizados de massa corporal total (MCT) dos estudantes do curso de Medicina da FMUSP. Estão disponíveis apenas o tamanho amostral, a média e o desvio-padrão amostral da MCT em cada grupo de sexo:
| Sexo | \(n\) | média | desvio-padrão |
|---|---|---|---|
| F | 230 | 57.63 | 9.05 |
| M | 312 | 71.59 | 12.09 |
Quais são a média e o desvio-padrão amostrais da MCT para o conjunto dos 542 estudantes?
A média total é a média ponderada das médias dos grupos:
\[ \bar{x} = \dfrac{ n_F\,\bar{x}_F+n_M\,\bar{x}_M }{ n_F+n_M } \]
A variância amostral total é:
\[ s^2 = \dfrac{ (n_F-1)s_F^2+(n_M-1)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-1 } \]
| \(n\) | média | desvio-padrão |
|---|---|---|
| 542 | 65.67 | 12.90 |
nF <- 230
nM <- 312
mF <- 57.63
mM <- 71.59
sF <- 9.05
sM <- 12.09
n <- nF + nM
m <- (nF*mF + nM*mM)/n
cat("Média = ", round(m, 2), "\n", sep="")Média = 65.67dp <- sqrt(
((nF - 1)*sF^2 +
(nM - 1)*sM^2 +
nF*(mF - m)^2 +
nM*(mM - m)^2) /
(n - 1)
)
cat("Desvio-padrão = ", round(dp, 2), "\n", sep="")Desvio-padrão = 12.9A variabilidade total pode ser decomposta em:
\[ \text{SQ}_{\text{total}} = \text{SQ}_{\text{dentro}} + \text{SQ}_{\text{entre}} \]
A variabilidade dentro dos grupos é:
\[ \text{SQ}_{\text{dentro}} = (n_F-1)s_F^2+(n_M-1)s_M^2 \]
A variabilidade entre os grupos é:
\[ \text{SQ}_{\text{entre}} = n_F\left(\bar{x}_F-\bar{x}\right)^2+ n_M\left(\bar{x}_M-\bar{x}\right)^2 \]
Dividindo todos os termos por \(n_F+n_M-1\):
\[ s^2 = s_{\text{dentro}}^2+ s_{\text{entre}}^2 \]
sendo que
\[ s_{\text{dentro}}^2 = \dfrac{ (n_F-1)s_F^2+(n_M-1)s_M^2 }{ n_F+n_M-1 } \]
e
\[ 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-1 } \]
A média harmônica dos tamanhos dos dois grupos é
\[ n_H=\dfrac{2n_Fn_M}{n_F+n_M} \]
Assim,
\[ n_F\left(\bar{x}_F-\bar{x}\right)^2+ n_M\left(\bar{x}_M-\bar{x}\right)^2 = \dfrac{n_H}{2} \left(\bar{x}_F-\bar{x}_M\right)^2 \] de modo que
\[ s_{\text{entre}}^2 = \dfrac{ \dfrac{n_H}{2}\left(\bar{x}_F-\bar{x}_M\right)^2 }{ n_F+n_M-1 } \]
Portanto, a variância amostral total pode ser escrita como
\[ s^2= \dfrac{ (n_F-1)s_F^2+(n_M-1)s_M^2+ \dfrac{n_H}{2} \left(\bar{x}_F-\bar{x}_M\right)^2 }{ n_F+n_M-1 } \]
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 e 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\)
sendo que \(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, \(\text{df}_A=g-1\), e do efeito residual, \(\text{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. \(\text{df}_A\) quantifica a complexidade do modelo (GLM univariado) e \(\text{df}_E\) quantifica o tamanho de amostra efetivo, sendo que \(n\) é o tamanho de amostra nominal.
\[
\begin{align}
F &\underset{H_0}{\sim} F_{\text{df}_A, \, \text{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 para fator com três grupos
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.
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?
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 1100 Brendon McGuirk Melissa
20 20 20 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
alpha <- 0.05
g <- nlevels(Dados$Instructor)
alphaBonf <- alpha/g
gplots::plotmeans(Sodium ~ Instructor,
connect=FALSE,
barcol="black",
p=1-alphaBonf,
main="IC95% Bonferroni",
data=Dados)
t_out <- stats::pairwise.t.test( # testes t de Student
Dados$Sodium,
Dados$Instructor,
p.adjust.method = "none",
pool.sd = FALSE,
var.equal = TRUE
)
print(t_out)
Pairwise comparisons using t tests with non-pooled SD
data: Dados$Sodium and Dados$Instructor
Brendon McGuirk
McGuirk 0.4477 -
Melissa 0.0043 0.0100
P value adjustment method: none 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.
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 entre participantes 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.
Hipó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
NormalityRegistered S3 methods overwritten by 'car':
method from
hist.boot FSA
confint.boot FSA
--------
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, \(\text{df}_A=g-1\), e do efeito residual, \(\text{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. \(\text{df}_A\) quantifica a complexidade do modelo (GLM univariado) e \(\text{df}_E\) quantifica o tamanho de amostra efetivo, sendo que \(n\) é o tamanho de amostra nominal.
\[ \begin{align} F &\underset{H_0}{\sim} F_{\text{df}_A, \, \text{df}_E}\\ &\underset{H_0}{\sim} F_{g-1, \, \sum_{i=1}^{g}{(n_i-1)}}\\ &\underset{H_0}{\sim} F_{g-1, \, (n-1)-(g-1)}\\ F &\underset{H_0}{\sim} F_{g-1, \, n-g}\\ \end{align} \]
Se \(\text{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 \(\text{df}_E < 30\), pode-se recorrer ao argumento de que a distribuição normal 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).
Algumas das funções principais são:
fit1 <- lm(Sodium ~ Instructor, data=Dados)
car::Anova(fit1)
oneway.test(Sodium ~ Instructor,
var.equal=TRUE,
data=Dados)
pairwise.t.test(
x = Dados$Sodium,
g = Dados$Instructor,
pool.sd = TRUE,
p.adjust.method = "holm")
fit2 <- onewaytests::onewaytests(Sodium ~ Instructor,
method="aov",
alpha=alfa,
na.rm=TRUE,
verbose=TRUE,
data=Dados)
onewaytests::paircomp(fit2, # Valores p incorretos
adjust.method="holm")
jmv::ANOVA(dep = "Sodium",
factors = "Instructor",
modelTes = TRUE,
postHoc= ~ Instructor,
emMeans = ~ Instructor,
ciWidthEmm = 100*(1-alfa),
emmPlots = TRUE,
emmTables = TRUE,
postHocCorr = "holm",
effectSize = "eta",
data = Dados) 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
ANOVAAnova 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.006235R^2 = eta^2 = 0.16eta2 <- effectsize::eta_squared(anv,
partial=FALSE,
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 | 95% CI
--------------------------------------
Instructor | 0.1632 | [0.0166, 0.3265]# gráfico do test F
reg <- summary(fit)
fobs <- reg$fstatistic[1]
dfn <- reg$fstatistic[2]
dfd <- reg$fstatistic[3]
fc <- qf(1-alfa, dfn, dfd)
p <- 1-pf(fobs, dfn, dfd)
if (p < 1e-4)
{
p <- sprintf("%.2e",p)
} else
{
p <- sprintf("%.4f",p)
}
f <- seq(0,1.4*max(fc,fobs),length.out=300)
densf <- df(f, dfn, dfd)
plot(f, densf,
main="One-way ANOVA: omnibus test",
xlab="F", ylab="Density",
lwd=2, type="l")
abline(v=fc, lty=3)
abline(v=fobs, lty=4)
legend("topright",
c("H0: igualdade das médias populacionais",
paste("F(",dfn,",",dfd,",",1-alfa,") = ",round(fc,3),sep=""),
paste("F(",dfn,",",dfd,") = ",round(fobs,3),"\n",
"p = ",p,sep="")
),
lwd=c(2,1,1), lty=c(1,3,4),
cex=0.8)
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: DunnettEMM.contrast <- emmeans::contrast(object = emm,
method = "trt.vs.ctrl",
ref = "Brendon",
adjust="holm",
level=1-alfa)
print(EMM.contrast) 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: ConsecutivoEMM.contrast <- emmeans::contrast(object=emm,
method="consec",
adjust="holm",
level=1-alfa)
print(EMM.contrast) 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 
# oneway.test: Fisher
fit1 <- oneway.test(Sodium ~ Instructor,
var.equal=TRUE,
data=Dados)
print(fit1, digits=5)
One-way analysis of means
data: Sodium and Instructor
F = 5.56, num df = 2, denom df = 57, p-value = 0.0062print(pairwise.t.test(
x = Dados$Sodium,
g = Dados$Instructor,
pool.sd = TRUE,
p.adjust.method = "holm"))
Pairwise comparisons using t tests with pooled SD
data: Dados$Sodium and Dados$Instructor
Brendon McGuirk
McGuirk 0.4228 -
Melissa 0.0066 0.0396
P value adjustment method: holm
Análise de significância prática: tamanho de efeitoeta2 <- effectsize::eta_squared(fit1,
partial=FALSE,
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
Eta2 | 95% CI
-------------------------
0.1632 | [0.0166, 0.3265]# onewaytests::onewaytests: Fisher
fit2 <- onewaytests::onewaytests(Sodium ~ Instructor,
method="aov",
alpha=alfa,
na.rm=TRUE,
verbose=TRUE,
data=Dados)
One-Way Analysis of Variance
data: Sodium and Instructor
F = 5.5579, num df = 2, denom df = 57, p-value = 0.006235
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
----------------------------------------------------- # jmv::ANOVA: Fisher
print(jmv::ANOVA(dep = "Sodium",
factors = "Instructor",
modelTes = TRUE,
postHoc= ~ Instructor,
emMeans = ~ Instructor,
ciWidthEmm = 100*(1-alfa),
emmPlots = TRUE,
emmTables = TRUE,
postHocCorr = "holm",
effectSize = "eta",
data = Dados))
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 0.1631905
Residuals 1487812.5 57 26101.97
──────────────────────────────────────────────────────────────────────────────────────────────
POST HOC TESTS
Post Hoc Comparisons - Instructor
──────────────────────────────────────────────────────────────────────────────────────────────────────
Instructor Instructor Mean Difference SE df t p-holm
──────────────────────────────────────────────────────────────────────────────────────────────────────
Brendon - McGuirk 41.25000 51.09009 57.00000 0.8073973 0.4227942
- Melissa 163.75000 51.09009 57.00000 3.2051225 0.0066393
McGuirk - Melissa 122.50000 51.09009 57.00000 2.3977252 0.0395910
──────────────────────────────────────────────────────────────────────────────────────────────────────
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
A ANOVA independente unifatorial de Welch é uma alternativa à ANOVA clássica de Fisher quando a hipótese de homogeneidade das variâncias entre as condições não é plausível.
A hipótese nula continua sendo a igualdade das médias populacionais:
\[ H_0:\mu_1=\mu_2=\cdots=\mu_g \]
Entretanto, diferentemente da ANOVA de Fisher, o teste de Welch não exige igualdade das variâncias populacionais. As médias dos grupos são comparadas utilizando pesos que dependem das variâncias amostrais e dos tamanhos amostrais de cada condição.
Para a condição \(i\), o peso é proporcional a
\[ w_i=\dfrac{n_i}{s_i^2} \]
Assim, grupos com maior precisão, isto é, maior tamanho amostral e menor variância, recebem maior peso na comparação das médias.
A ANOVA de Welch é especialmente útil quando há heterocedasticidade, principalmente se os tamanhos amostrais também forem diferentes entre as condições.
ANOVA independente unifatorial de Welch demo_ANOVA1f_indep_Welch_sodio.R pode
ser executada pelas funções jmv::anovaOneW(),
stats::oneway.test() ou
rstatix::games_howell_test().
oneway.test(
Sodium ~ Instructor,
var.equal = FALSE,
data = Dados
)
rstatix::welch_anova_test(Sodium ~ Instructor,
data=Dados)
rstatix::games_howell_test(
Sodium ~ Instructor,
conf.level = 1 - alfa,
detailed = FALSE,
data = Dados
)
fit <- onewaytests::onewaytests(Sodium ~ Instructor,
method="welch",
alpha=alfa,
na.rm=TRUE,
verbose=TRUE,
data=Dados)
onewaytests::paircomp(fit, # Valores p incorretos
adjust.method = "holm")
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# oneway.test: Welch
fit <- oneway.test(Sodium ~ Instructor,
var.equal=FALSE,
data=Dados)
print(fit, digits=5)
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.0066out <- rstatix::games_howell_test(Sodium ~ Instructor,
conf.level=1-alfa,
detailed=FALSE,
data=Dados)
print(data.frame(out), digits=2) .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 efeitoeta2 <- effectsize::eta_squared(fit,
partial=FALSE,
generalized=FALSE,
ci=1-alfa,
alternative="two.sided",
verbose=TRUE)`var.equal = FALSE` - effect size is an approximation.# Effect Size for ANOVA
Eta2 | 95% CI
-------------------------
0.2357 | [0.0251, 0.4351]# jmv::anovaOneW: Welch
print(jmv::anovaOneW(
data = Dados,
dep = "Sodium",
group = "Instructor",
welch = TRUE,
phMethod = "gamesHowell",
phTest = TRUE,
phFlag = TRUE,
desc = FALSE
))
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# rstatix::welch_anova_test: Welch
fit2 <- rstatix::welch_anova_test(Sodium ~ Instructor,
data=Dados)
print(data.frame(fit2), digits=4) .y. n statistic DFn DFd p method
1 Sodium 60 5.77 2 37.4 0.007 Welch ANOVAout <- rstatix::games_howell_test(Sodium ~ Instructor,
data=Dados,
conf.level=1-alfa,
detailed=FALSE)
print(data.frame(out), digits=2) .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 *# onewaytests::onewaytests: Welch
fit3 <- onewaytests::onewaytests(Sodium ~ Instructor,
method="welch",
alpha=alfa,
na.rm=TRUE,
verbose=TRUE,
data=Dados)
Welch's Heteroscedastic F Test
data: Sodium and Instructor
F = 5.7655, num df = 2.000, denom df = 37.396, p-value = 0.006569
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
----------------------------------------------------- Quando o teste global de Welch é significante, comparações múltiplas entre pares de médias podem ser realizadas pelo procedimento de Games-Howell, que também não pressupõe homogeneidade das variâncias.
A conclusão, neste exemplo, é semelhante à obtida com a ANOVA de Fisher: o programa adotado por Melissa apresentou ingestão diária média de sódio significantemente menor que a dos programas dos outros dois instrutores.
A matemática da ANOVA 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.006569No caso balanceado e homocedástico, considere o caso em que todos os grupos têm o mesmo tamanho e a mesma variância amostral:
\[ n_i=m \qquad\text{e}\qquad s_i^2=s^2 \]
Nesse caso, todos os pesos são iguais e
\[ \frac{w_i}{w}=\frac{1}{g} \]
Substituindo na fórmula de Welch–Satterthwaite:
\[ \nu = \frac{g^2-1} {3g\dfrac{\left(1-1/g\right)^2}{m-1}} = \frac{g(g+1)(m-1)} {3(g-1)} \]
Como o tamanho amostral total é
\[ N=gm \]
temos
\[ N-g=g(m-1) \]
e, portanto,
\[ \boxed{ \nu_{\max} = \frac{(g+1)(N-g)}{3(g-1)} } \]
Assim, mesmo quando as amostras são perfeitamente balanceadas e apresentam variâncias amostrais iguais, o \(\text{df}\) do denominador do teste de Welch não é, em geral, igual ao \(\text{df}\) residual da ANOVA clássica, que é
\[ N-g \]
Por exemplo, para \(g=3\) grupos e \(N=120\) observações,
\[ \nu_{\max} = \frac{4}{6}(120-3) = 78 \]
Portanto,
\[ \text{df}_{\text{Welch}}\leq 78 \]
enquanto a ANOVA clássica teria
\[ \text{df}_{\text{resíduo}}=120-3=117 \]
No extremo oposto, suponha que um único grupo concentre praticamente todo o peso:
\[ \frac{w_1}{w}\to1 \]
enquanto, para os demais grupos,
\[ \frac{w_i}{w}\to0 \]
Se esses grupos tiverem o menor tamanho amostral possível, \(n_i=2\), então
\[ \nu_{\min} = \frac{g^2-1}{3(g-1)} = \boxed{\frac{g+1}{3}} \]
Para \(g=3\),
\[ \nu_{\min} = \frac{4}{3} \approx1.33 \]
Portanto, para três grupos e \(N=120\), os graus de liberdade do denominador podem, teoricamente, variar aproximadamente no intervalo
\[ \boxed{ 1.33 < \nu \leq 78 }. \]
O valor \(117=N-g\) é o número de graus de liberdade residual da ANOVA clássica, e não o número de graus de liberdade do teste de Welch.
Interpretação do resultado:
O teste de Welch resultou em
\[ F(2,37.396)=5.7655, \qquad p=0.006569 \]
Como \(p<0.05\), rejeita-se a hipótese nula \(H_0:\mu_1=\mu_2=\mu_3\).
Há evidência estatística de que as médias populacionais não são iguais entre os três grupos.
O número de graus de liberdade do numerador, \(\text{df}_1=2\), corresponde a \(g-1=3-1=2\).
O número de graus de liberdade do denominador, \(\text{df}_2=37.396\), é obtido pela aproximação de Welch–Satterthwaite e, por isso, não precisa ser inteiro.
Esse teste não pressupõe igualdade das variâncias entre os grupos. Assim, a conclusão sobre a diferença entre as médias permanece apropriada mesmo na presença de heterocedasticidade.
O resultado global indica apenas que pelo menos uma média difere das demais. Para identificar quais grupos diferem entre si, é necessário realizar comparações múltiplas compatíveis com heterocedasticidade, como o teste de Games–Howell.
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 Fisher-White
A ANOVA independente unifatorial clássica de Fisher testa a hipótese de igualdade das médias populacionais entre \(g\) condições independentes:
\[ H_0:\mu_1=\mu_2=\cdots=\mu_g \]
O teste clássico pressupõe homocedasticidade, isto é, igualdade das variâncias populacionais entre as condições. Quando essa suposição é questionável, pode-se utilizar uma versão robusta à heterocedasticidade baseada no estimador de White da matriz de covariâncias dos coeficientes do modelo linear.
O modelo continua sendo
\[ Y_{ij}=\beta_0+\beta_i+\varepsilon_{ij} \]
mas a matriz de covariâncias dos estimadores dos coeficientes é
corrigida para heterocedasticidade. A correção HC2 ajusta
os resíduos também em função das alavancagens das observações.
No R, o procedimento pode ser realizado ajustando inicialmente o modelo linear e aplicando posteriormente o teste \(F\) robusto:
A hipótese testada permanece sendo a igualdade das médias:
\[ H_0:\mu_1=\mu_2=\cdots=\mu_g \]
A diferença em relação à ANOVA de Fisher clássica está no cálculo da variância dos estimadores e, consequentemente, da estatística de teste. A correção de White-HC2 torna a inferência robusta à heterocedasticidade, mas não elimina a necessidade de independência entre as unidades experimentais.
demo_ANOVA1f_indep_FisherWhite.R
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
A ANOVA independente unifatorial por bootstrapping utiliza reamostragem para aproximar a distribuição amostral da estatística de teste, reduzindo a dependência de aproximações assintóticas da ANOVA clássica.
A hipótese nula continua sendo
\[ H_0:\mu_1=\mu_2=\cdots=\mu_g \]
No bootstrap, novas amostras são geradas repetidamente a partir dos dados observados ou dos resíduos do modelo. Para cada reamostra, calcula-se novamente a estatística de interesse. A distribuição das estatísticas obtidas nas \(B\) reamostras é então utilizada para realizar a inferência.
A função lmboot::ANOVA.boot() permite utilizar
bootstrap dos resíduos ou wild bootstrap. No
bootstrap dos resíduos, utilizado abaixo, os resíduos do modelo
são reamostrados e adicionados aos valores ajustados para gerar novas
respostas.
demo_ANOVA1f_indep_Bootstrap_sodio.R
implementa o procedimento por meio da função
lmboot::ANOVA.boot(). Essa função realiza o teste global,
mas não fornece diretamente comparações post hoc.
lmboot::ANOVA.boot(
Sodium ~ Instructor,
data = Dados,
B = B,
type = "residual",
seed = 123,
keep.boot.resp = FALSE
)Instructor
Brendon McGuirk Melissa
20 20 20
Bootstrap One-Way ANOVA
F(2, 57) = 5.56
p bootstrap = 0.0052
Bootstrap samples = 10000
Effect size analysis
eta^2 = 0.16O argumento wild.dist não é necessário quando
type="residual". Ele somente é utilizado quando se
especifica
O bootstrap dos resíduos é mais diretamente compatível com a hipótese de homocedasticidade do modelo. Quando há preocupação com heterocedasticidade, uma alternativa é o wild bootstrap:
lmboot::ANOVA.boot(
Sodium ~ Instructor,
data = Dados,
B = B,
type = "wild",
wild.dist = "normal",
seed = 123,
keep.boot.resp = FALSE
)Outra possibilidade é utilizar o pacote onewaytests. O
método "gtb" corresponde a um teste generalizado
equivalente ao bootstrap paramétrico para comparação de médias
sob variâncias possivelmente diferentes.
fit <- onewaytests::onewaytests(
Sodium ~ Instructor,
method = "gtb",
alpha = alfa,
na.rm = TRUE,
N = B,
verbose = TRUE,
data = Dados
)
onewaytests::paircomp( # valores p incorretos
fit,
adjust.method = "holm"
)Nesse caso, onewaytests::paircomp() realiza as
comparações dois a dois e ajusta os valores de \(p\) pelo método especificado, aqui
Holm.
Uma terceira alternativa é o pacote WRS2, que implementa
procedimentos robustos por bootstrap para médias aparadas.
demo_WRS2_ANOVA_independente.R
A função WRS2::t1waybt() executa o teste global e
fornece, entre outras medidas, a proporção de variância explicada. O
argumento tr determina a proporção aparada em cada cauda.
Com tr=0, são utilizadas as médias aritméticas usuais.
A função WRS2::mcppb20() realiza comparações múltiplas
por bootstrap:
Embora mcppb20() aceite tr=0, o
procedimento foi originalmente formulado para médias aparadas e a
própria função produz uma advertência quando tr difere de
0.20. Portanto, seu uso com tr=0 deve ser interpretado com
cautela.
Instructor
Brendon McGuirk Melissa
20 20 20
Generalized test equivalent to parametric bootstrap
Generalized Test Equivalent to Parametric Bootstrap Test (size close to
intended)
data: Sodium and Instructor
p-value = 0.006664
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 ANOVAWarning in WRS2::t1waybt(Sodium ~ Instructor, data = Dados, tr = 0.2, nboot =
B): Some bootstrap estimates of the test statistic could not be computed.Call:
WRS2::t1waybt(formula = Sodium ~ Instructor, data = Dados, tr = 0.2,
nboot = B)
Effective number of bootstrap samples was 99990.
Test statistic: 3.7723
p-value: 0.04828
Variance explained: 0.283
Effect size: 0.532
Call:
WRS2::mcppb20(formula = Sodium ~ Instructor, data = Dados, tr = 0.2,
nboot = B)
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.02193As diferentes implementações de bootstrap não são
estatisticamente idênticas. lmboot::ANOVA.boot() com
type="residual" reamostra resíduos do modelo;
lmboot::ANOVA.boot() com type="wild" utiliza
wild bootstrap; onewaytests com
method="gtb" utiliza um teste generalizado equivalente ao
bootstrap paramétrico; e WRS2::t1waybt() utiliza
um procedimento bootstrap robusto para médias aparadas.
Consequentemente, pequenas diferenças entre os valores de \(p\) obtidos por esses procedimentos são esperadas.
Neste exemplo, os métodos conduzem à mesma conclusão geral de que existem diferenças nas médias de ingestão diária de sódio entre os programas dos instrutores. As comparações dois a dois devem ser interpretadas utilizando os valores de \(p\) ou intervalos de confiança já corrigidos para multiplicidade pelo procedimento adotado.
Comparação dos resultados
Os três procedimentos rejeitam a hipótese de igualdade entre os grupos ao nível de significância de 5%, mas testam parâmetros diferentes.
O lmboot::ANOVA.boot() e o
onewaytests::onewaytests(method="gtb") comparam as médias
aritméticas e produzem resultados muito semelhantes:
\[ p=0.0052 \]
e
\[ p=0.0067 \]
As comparações múltiplas indicam que Melissa apresenta ingestão média de sódio menor que Brendon e McGuirk, enquanto Brendon e McGuirk não diferem entre si.
O WRS2::t1waybt() com tr=0.2 compara médias
aparadas em 20% e apresenta resultado global mais próximo do limite de
5%:
\[ p=0.0483 \]
Nas comparações múltiplas robustas, apenas Brendon e Melissa diferem após o controle da multiplicidade; a diferença entre McGuirk e Melissa não permanece significante.
Assim, a conclusão sobre Melissa é mais forte quando se comparam médias aritméticas do que quando se utilizam médias aparadas em 20%.
Comparativo (ANOVA independente)
Para facilitar a comparação das variantes de ANOVA independente unifatorial (omnibus) apresentadas acima:
| ANOVA | Estatística | \(\text{df}\) | \(p\) | Função |
|---|---|---|---|---|
| Fisher | \(F=5.56\) | 2, 57 | 0.0062 | jmv::ANOVA |
| Fisher | \(F=5.56\) | 2, 57 | 0.0062 | lm, car::Anova |
| Welch | \(F=5.77\) | 2, 37.4 | 0.0066 | oneway.test |
| Welch | \(F=5.77\) | 2, 37.4 | 0.0066 | rstatix::welch_anova_test |
| Welch | \(F=5.77\) | 2, 37.4 | 0.0066 | onewaytests::onewaytests |
| Fisher-White | \(F=5.87\) | 2, 57 | 0.0048 | lm, car::Anova |
| Bootstrap residual | \(F=5.56\) | 2, 57 | 0.0052 | lmboot::ANOVA.boot |
| Bootstrap paramétrico (GTB) | — | — | 0.0067 | onewaytests::onewaytests |
| Bootstrap robusto, média aparada 20% | \(T=3.7723\) | — | 0.0483 | WRS2::t1waybt |
Os métodos de Fisher, Welch, Fisher-White, bootstrap residual e bootstrap paramétrico apresentam resultados muito semelhantes e rejeitam claramente a hipótese de igualdade das médias ao nível de 5%.
O resultado de WRS2::t1waybt() não é diretamente
comparável aos demais, pois com tr=0.2 testa a igualdade
das médias aparadas em 20%. Nesse caso, a evidência contra \(H_0\) é consideravelmente menor, com
\[ p=0.0483 \]
Portanto, a conclusão de diferença entre os grupos é robusta entre os métodos baseados nas médias aritméticas, enquanto se torna limítrofe quando a localização dos grupos é representada pelas médias aparadas em 20%.
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 paramétrico
(
onewaytests::onewaytests+onewaytests::paircomp, ajuste de Holm)
| Comparação | \(p\) ajustado |
|---|---|
| Brendon - McGuirk | 0.4482 |
| Brendon - Melissa | 0.0134 |
| McGuirk - Melissa | 0.0200 |
Ao nível de significância de 5%, Melissa difere de Brendon e de McGuirk, enquanto Brendon e McGuirk não diferem entre si.
- Bootstrap robusto para médias aparadas em 20%
(
WRS2::t1waybt+WRS2::mcppb20)
| Comparação | \(p\) |
|---|---|
| Brendon - McGuirk | 0.3691 |
| Brendon - Melissa | 0.0047 |
| McGuirk - Melissa | 0.0219 |
Os valores de \(p\) apresentados por
WRS2::mcppb20() não são diretamente ajustados para
comparações múltiplas. A decisão deve considerar os intervalos de
confiança simultâneos fornecidos pela função. Nesse procedimento, apenas
Brendon e Melissa apresentam diferença significante após o controle da
multiplicidade.
ANOVA independente SEM dados individualizados
É 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 individualizados 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 individualizados 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 individualizados\n")
ANOVA independente unifatorial de Welch sem dados individualizados# 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 suposições para validade do teste F da ANOVA relacionada unifatorial são:
- 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
Dados
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 1100GLM e GLMM univariados omnibus
Considere os dois modelos:
# GLMM
fit1 <- lmerTest::lmer(Sodium ~ Instructor + (1 | Student),
data = Dados)
car::Anova(fit1,
test.statistic = "F")
# GLM
fit2 <- lm(Sodium ~ Instructor + Student,
data = Dados)
car::Anova(fit2,
type = 2)No GLMM, o efeito da unidade experimental é representado por um intercepto aleatório:
\[ Y_{ij}=\mu+\alpha_j+b_i+\varepsilon_{ij} \]
com
\[ b_i\sim N(0,\sigma_b^2) \]
No GLM, o efeito da unidade experimental é introduzido como efeito fixo:
\[ Y_{ij}=\mu+\alpha_j+\beta_i+\varepsilon_{ij} \]
Os dois modelos podem produzir o mesmo teste omnibus para
Instructor quando o delineamento é completo e perfeitamente
balanceado, isto é, quando:
- cada unidade experimental é observada em todas as \(q\) condições dependentes
- todas as unidades experimentais possuem o mesmo número de observações
- não há valores ausentes
- o GLMM contém apenas intercepto aleatório por unidade experimental
- a variância residual é homogênea entre as condições
- não há inclinações aleatórias
- a estrutura de covariância intraunidade é compatível com a estrutura de intercepto aleatório
- o teste de interesse é o efeito fixo de
Instructor
Nessas condições, o efeito da unidade experimental é removido da comparação entre as condições tanto no GLM quanto no GLMM, e os dois modelos utilizam essencialmente a mesma informação para testar
\[ H_0:\mu_1=\mu_2=\cdots=\mu_q \]
No GLM perfeitamente balanceado,
\[ \text{df}_E=(n-1)(q-1) \]
No GLMM perfeitamente balanceado, os graus de liberdade obtidos por Satterthwaite ou Kenward–Roger podem coincidir com esse valor.
A equivalência deixa de ser garantida quando há dados ausentes, números diferentes de observações por unidade experimental, heterogeneidade de variâncias, estruturas de covariância mais complexas ou inclinações aleatórias.
Assim, em delineamentos completos e perfeitamente balanceados, o GLM
com Student como efeito fixo e o GLMM com
(1 | Student) podem fornecer resultados omnibus
equivalentes para o efeito de Instructor. Em dados
desbalanceados ou com estrutura de dependência mais complexa, o GLMM é
mais geral.
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# ============================================================
# GLM: Student como efeito fixo
# ============================================================
cat("\nGLM\n")
GLMmodelo.glm <- stats::lm(
Sodium ~ Instructor + Student,
data = Dados
)
# ANOVA tipo II
print(
car::Anova(
modelo.glm,
type = 2
),
digits = 4
)Anova 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# ============================================================
# GLMM: Student como intercepto aleatório
# ============================================================
cat("\nGLMM\n")
GLMMRegistered S3 method overwritten by 'lme4':
method from
na.action.merMod car # car::Anova: teste F
print(
car::Anova(
modelo.glmm,
type = 2,
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 ' ' 1# Kenward-Roger
# Requer o pacote pbkrtest
print(
stats::anova(
modelo.glmm,
type = 2,
ddf = "Kenward-Roger"
),
digits = 4
)Type II 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 II 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 <0.0001
McGuirk - Melissa 122.5 21.8 38 67.9 177.1 5.624 <0.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. Graus de liberdade
O 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_q\), a estatística de teste F segue uma
distribuição F de Fisher-Snedecor com graus de liberdade \(\text{df}_A=q-1\) e \(\text{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. Essa expressão pressupõe que o modelo aditivo
y ~ Condition + Subject tenha o posto esperado, isto é, que
o delineamento seja conectado.
\[ \begin{align} F &\underset{H_0}{\sim} F_{\text{df}_A, \, \text{df}_E}\\ &\underset{H_0}{\sim} F_{q-1, \, (\sum_{i=1}^{n}{q_i}-1)-((n-1)+(q-1))}\\ &\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} \]
Como \(r=n\bar q\), sendo que \(\bar q=\frac{1}{n}\sum_{i=1}^{n}{q_i}\), temos exatamente \(\text{df}_E=n(\bar q-1)-(q-1)\). Se \(n(\bar q-1)\gg(q-1)\), então \(\text{df}_E\approx n(\bar q-1)\). Se os dados são perfeitamente balanceados, \(\bar q=q\) e, para \(n\) suficientemente grande, \(\text{df}_E\approx n(q-1)\).
Se os dados são completos (perfeitamente balanceados):
Modelo Linear Geral (GLM): lm:
- \(\text{df}_E = (n-1)(q-1) = (nq-1)-((n-1)+(q-1))\) exatamente. No estudo do exemplo, \(\text{df}_E = (20-1)(3-1) = (20\times 3-1)-((20-1)+(3-1))=38\).
Modelo Linear Misto (LMM): lmerTest::lmer:
- \(\text{df}_E\) é aproximado por Satterthwaite ou Kenward–Roger. Não há, em geral, uma expressão algébrica exata para \(\text{df}_E\) equivalente a \((n-1)(q-1)\). No estudo do exemplo, o valor de \(\text{df}_E\) deve ser obtido diretamente do modelo ajustado. Em dados perfeitamente balanceados, as aproximações podem coincidir com \((n-1)(q-1)\).
Para o Modelo Linear Geral com dados completos e perfeitamente balanceados:
\[ \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 \(\text{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 \(\text{df}_2\) no modelo
y ~ Condition + Subject:
- \(\text{df}_{\text{total}} = r - 1\)
- \(\text{df}_{\text{subject}} = n - 1\)
- \(\text{df}_{\text{condition}} = q - 1\)
- \(\text{df}_{\text{residual}} = r - n - q + 1\)
Observe que:
\[ \text{df}_{\text{residual}} = \text{df}_{\text{total}} - (\text{df}_{\text{subject}}+\text{df}_{\text{condition}}) \]
Caso balanceado (cada unidade experimental em todas as \(q\) condições dependentes), \(r=nq\) e então:
\[ \begin{align} \text{df}_{\text{residual}} &= (r-1)-(n-1)-(q-1)\\ &= (nq-1)-(n-1)-(q-1)\\ \text{df}_{\text{residual}} &= (n-1)(q-1) \end{align} \]
DADOS BALANCEADOSn = 20 q = 3 r = 60 df total = 59 df Subject = 19 df Condition = 2 df residual: decomposição = 38 df residual: r-n-q+1 = 38 df residual: (n-1)(q-1) = 38 MODELO LINEAR GERAL: DADOS BALANCEADOSAnova Table (Type II tests)
Response: measure
Sum Sq Df F value Pr(>F)
Condition 148.37 2 10.19 0.0002863 ***
Subject 1464.86 19 10.59 6.505e-10 ***
Residuals 276.65 38
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df residual fornecido por stats::lm = 38 Posto da matriz X = 22 n + q - 1 = 22 r - rank(X) = 38 Todas as fórmulas coincidem exatamente.Sob H0: F ~ F(2, 38)DADOS DESBALANCEADOSn = 20 q = 3 r = 55 q médio = 2.75 df residual pela decomposição = 33 df residual por r-n-q+1 = 33 df residual por n(qbar-1)-(q-1) = 33 df aproximado por n(qbar-1) = 35 MODELO LINEAR GERAL: DADOS DESBALANCEADOSAnova Table (Type II tests)
Response: measure
Sum Sq Df F value Pr(>F)
Condition 180.94 2 14.337 3.302e-05 ***
Subject 1286.80 19 10.732 3.447e-09 ***
Residuals 208.24 33
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df residual fornecido por stats::lm = 33 Posto da matriz X = 22 n + q - 1 = 22 r - rank(X) = 33 MODELO LINEAR MISTO: DADOS BALANCEADOSSatterthwaiteType III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Condition 148.37 74.186 2 38 10.19 0.0002863 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Kenward-RogerType III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Condition 148.37 74.186 2 38 10.19 0.0002863 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MODELO LINEAR MISTO: DADOS DESBALANCEADOSSatterthwaiteType III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Condition 184.15 92.076 2 33.426 14.602 2.77e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Kenward-RogerType III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Condition 183.74 91.871 2 33.382 14.569 2.829e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANOVA 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 na condição
Melissa, conforme a planilha Nutricao3rmdesb.xlsx.
alfa <- 0.05
Dados <- data.frame(
readxl::read_excel("Nutricao3rmdesb.xlsx")
)
Dados$Student <- factor(Dados$Student)
Dados$Instructor <- factor(
Dados$Instructor,
labels = c("Brendon", "McGuirk", "Melissa")
)
saveRDS(Dados, "Nutricao3rmdesb.rds")Com a exclusão de uma observação, o delineamento deixa de ser perfeitamente balanceado. Nesse caso, o GLMM é particularmente útil, pois permite utilizar todas as observações disponíveis sem exigir que cada unidade experimental tenha medidas em todas as condições.
Utilizamos o código R implementado em demo_ANOVA1f_dep_desbalanc_sodio.R.
fit1 <- lmerTest::lmer(
Sodium ~ Instructor + (1 | Student),
data = Dados
)
car::Anova(
fit1,
test.statistic = "F"
)
fit2 <- stats::lm(
Sodium ~ Instructor + Student,
data = Dados
)
car::Anova(
fit2,
type = 2
)- GLM e GLMM omnibus
source("summarySEwithin2.R")
Dados <- readRDS("Nutricao3rmdesb.rds")
Dados$Student <- factor(Dados$Student)
Dados$Instructor <- factor(Dados$Instructor)
print(
stats::xtabs(~ Instructor + Student, data = Dados)
) 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 1alfa <- 0.05
g <- nlevels(Dados$Instructor)
# ============================================================
# INTERVALOS DE CONFIANÇA INTRAUNIDADES
# Bonferroni para os g intervalos das médias
# ============================================================
alfaBonf <- alfa/g
ic <- summarySEwithin2(
Dados,
measurevar = "Sodium",
withinvars = "Instructor",
idvar = "Student",
na.rm = TRUE,
conf.interval = 1 - alfaBonf
)
print(ic, digits = 3) 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.0grf <- ggplot2::ggplot(
ic,
ggplot2::aes(
x = Instructor,
y = Sodium
)
) +
ggplot2::geom_errorbar(
width = 0.1,
ggplot2::aes(
ymin = Sodium - ci,
ymax = Sodium + ci
)
) +
ggplot2::geom_point(
shape = 21,
size = 3,
fill = "white"
) +
ggplot2::ylab("Sodium") +
ggplot2::ggtitle(
paste0(
g,
" Instructors: Sodium\n",
"Within-subject CI95% Bonferroni"
)
) +
ggplot2::theme_bw()
print(grf)
# ============================================================
# GLM
# Student como efeito fixo
# ============================================================
cat("\nGLM\n")
GLMmodelo.glm <- stats::lm(
Sodium ~ Instructor + Student,
data = Dados
)
anv.glm <- car::Anova(
modelo.glm,
type = 2
)
print(anv.glm, digits = 4)Anova Table (Type II tests)
Response: Sodium
Sum Sq Df F value Pr(>F)
Instructor 256495 2 27.55 4.71e-08 ***
Student 1262858 19 14.28 1.09e-11 ***
Residuals 172255 37
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ============================================================
# GLMM
# Student como intercepto aleatório
# ============================================================
cat("\nGLMM\n")
GLMMmodelo.glmm <- lmerTest::lmer(
Sodium ~ Instructor + (1 | Student),
data = Dados
)
# ============================================================
# TESTE OMNIBUS GLMM
# Kenward-Roger
# ============================================================
cat("\nANOVA\n")
ANOVAAnalysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Sodium
F Df Df.res Pr(>F)
Instructor 27.36 2 37.06 5.04e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ============================================================
# REGRESSÃO GLMM
# ============================================================
cat("\nRegression\n")
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# ============================================================
# TAMANHO DE EFEITO GLMM
# ============================================================
cat("\nEffect size analysis\n")
Effect size analysiseta2 <- effectsize::eta_squared(
anv.glmm,
partial = TRUE,
generalized = FALSE,
ci = 1 - alfa,
alternative = "two.sided",
verbose = TRUE
)
print(eta2, digits = 4)# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI
----------------------------------------------
Instructor | 0.5962 | [0.3729, 0.7253]# ============================================================
# MÉDIAS MARGINAIS ESTIMADAS
# ============================================================
cat("\nEstimated marginal means\n")
Estimated marginal meansemm <- emmeans::emmeans(
modelo.glmm,
specs = ~ Instructor,
lmer.df = "satterthwaite",
lmerTest.limit = nrow(Dados)
)
print(
summary(
emm,
infer = c(TRUE, FALSE),
level = 1 - alfa
)
) 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 # ============================================================
# TESTES POST HOC
# ============================================================
cat("\nPost hoc tests\n")
Post hoc testscontrastes <- emmeans::contrast(
emm,
method = "pairwise"
)
# Valores p ajustados por Holm
print(
summary(
contrastes,
infer = c(FALSE, TRUE),
adjust = "holm"
)
) contrast estimate SE df t.ratio p.value
Brendon - McGuirk 41.2 21.6 37.0 1.912 0.0636
Brendon - Melissa 157.5 22.0 37.1 7.166 <0.0001
McGuirk - Melissa 116.3 22.0 37.1 5.289 <0.0001
Degrees-of-freedom method: satterthwaite
P value adjustment: holm method for 3 tests # Intervalos de confiança simultâneos de Bonferroni
print(
stats::confint(
contrastes,
level = 1 - alfa,
adjust = "bonferroni"
)
) contrast estimate SE df lower.CL upper.CL
Brendon - McGuirk 41.2 21.6 37.0 -12.8 95.3
Brendon - Melissa 157.5 22.0 37.1 102.4 212.6
McGuirk - Melissa 116.3 22.0 37.1 61.1 171.4
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates # ============================================================
# GRÁFICOS DAS MÉDIAS MARGINAIS E DOS CONTRASTES
# ============================================================
print(
plot(
emm,
colors = "black"
) +
ggplot2::theme_bw()
)
# ============================================================
# COMPACT LETTER DISPLAY
# ============================================================
print(
multcomp::cld(
object = emm,
level = 1 - alfa,
adjust = "holm",
Letters = letters,
alpha = alfa
)
) 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. # ============================================================
# COMPARAÇÃO FINAL: GLM versus GLMM
# ============================================================
tab.glm <- as.data.frame(anv.glm)
tab.glmm <- as.data.frame(anv.glmm)
F.glm <- tab.glm["Instructor", "F value"]
df1.glm <- tab.glm["Instructor", "Df"]
df2.glm <- stats::df.residual(modelo.glm)
p.glm <- tab.glm["Instructor", "Pr(>F)"]
F.glmm <- tab.glmm["Instructor", "F"]
df1.glmm <- tab.glmm["Instructor", "Df"]
df2.glmm <- tab.glmm["Instructor", "Df.res"]
p.glmm <- tab.glmm["Instructor", "Pr(>F)"]
Comparacao <- data.frame(
"Modelo" = c("GLM", "GLMM"),
"Student" = c("Efeito fixo", "Intercepto aleatório"),
"F" = c(F.glm, F.glmm),
"df1" = c(df1.glm, df1.glmm),
"df2" = c(df2.glm, df2.glmm),
"p" = c(p.glm, p.glmm),
check.names = FALSE
)
cat("\nCOMPARAÇÃO GLM versus GLMM\n")
COMPARAÇÃO GLM versus GLMM Modelo Student F df1 df2 p
GLM Efeito fixo 27.55 2 37.00 4.714e-08
GLMM Intercepto aleatório 27.36 2 37.06 5.043e-08No GLM, Student é tratado como efeito fixo e, com dados
desbalanceados, utilizamos car::Anova(type=2) para que o
teste de Instructor não dependa da ordem dos termos no
modelo.
No GLMM, Student é tratado como intercepto aleatório. Os
graus de liberdade do denominador são aproximados por Kenward–Roger ou
Satterthwaite, permitindo acomodar o desbalanceamento decorrente da
observação ausente.
|
A ANOVA relacionada unifatorial desbalanceada feita com o pacote
|
Comparação entre GLM e GLMM com dados desbalanceados
Com uma observação ausente, o delineamento deixa de ser perfeitamente balanceado e a equivalência exata entre GLM e GLMM deixa de ser garantida.
No GLM,
Student é tratado como efeito fixo e o número de graus
de liberdade residual é exato:
\[ \text{df}_E=r-n-q+1=59-20-3+1=37 \]
O teste omnibus fornece
\[ F(2,37)=27.55 \qquad p<0.0001 \]
No GLMM,
Student é tratado como efeito aleatório e os graus de
liberdade do denominador são aproximados por Satterthwaite ou
Kenward–Roger. Neste exemplo, Kenward–Roger fornece
\[ F(2,37.06)=27.36 \qquad p<0.0001 \]
Os resultados são muito próximos:
| Modelo | Student |
\(F\) | \(\text{df}_1\) | \(\text{df}_2\) | \(p\) |
|---|---|---|---|---|---|
| GLM | efeito fixo | 27.55 | 2 | 37.00 | \(4.714\times10^{-8}\) |
| GLMM | intercepto aleatório | 27.36 | 2 | 37.06 | \(5.043\times10^{-8}\) |
A perda de uma única observação produz, portanto, apenas uma pequena diferença entre os resultados dos dois modelos. No GLM, \(\text{df}_2=37\) é determinado exatamente pela decomposição dos graus de liberdade; no GLMM, \(\text{df}_2=37.06\) é aproximado por Kenward–Roger.
O GLMM também estima separadamente a variabilidade entre estudantes e a variabilidade residual:
\[ \hat\sigma^2_{\text{Student}}=20846 \qquad \hat\sigma^2_E=4654 \]
Ambos os modelos utilizam as 59 observações disponíveis. Entretanto, o GLMM é mais geral porque modela explicitamente a dependência intraestudante por meio do efeito aleatório e permite extensões para estruturas de dependência mais complexas.
O teste omnibus indica diferença entre os instrutores nos dois modelos. No GLMM, as comparações post hoc com correção de Holm mostram que Melissa apresenta ingestão média de sódio menor que Brendon e McGuirk, enquanto Brendon e McGuirk não diferem significantemente.
Assim, com desbalanceamento pequeno, GLM e GLMM produzem resultados praticamente equivalentes neste exemplo, embora essa equivalência não seja garantida para desbalanceamentos maiores ou estruturas de covariância mais complexas.
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
------------------------------------------------
factor(Year) | 0.8702 | [0.7167, 0.9254]
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 <0.0001
Year2 - Year5 -47.57 5.01 18 -62.4 -32.73 -9.499 <0.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 <0.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 <0.0001
Year5 - Year2 47.6 5.01 18 9.499 <0.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 <0.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
------------------------------------------------
factor(Year) | 0.8614 | [0.6980, 0.9203]
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 <0.0001
Year2 - Year5 -33.9 3.55 24 -44.1 -23.65 -9.537 <0.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 <0.0001
Year5 - Year2 33.9 3.55 24 9.537 <0.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 <0.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
Lógica estatística da estatística de teste F
Considere um delineamento entre participantes com \(g\geq3\) condições independentes e perfeitamente balanceado, com \(n\) unidades experimentais em cada condição. Portanto, o tamanho amostral total é
\[ N=gn \]
A ANOVA unifatorial compara duas formas de estimar a variância populacional \(\sigma^2\) da variável dependente.
A ideia básica é estimar \(\sigma^2\) por dois métodos diferentes: um que não depende da veracidade de \(H_0\) e outro que depende.
Se todas as condições tiverem a mesma média populacional,
\[ H_0:\mu_1=\mu_2=\cdots=\mu_g \]
as duas estimativas de \(\sigma^2\) deverão ser próximas. Se as médias populacionais diferirem, a estimativa baseada na variabilidade entre as médias amostrais tenderá a ser maior.
Uma primeira forma de estimar \(\sigma^2\), independentemente da veracidade de \(H_0\), consiste em calcular a variância amostral da variável dependente em cada condição,
\[ S_1^2,S_2^2,\ldots,S_g^2 \]
e, como o delineamento é perfeitamente balanceado, tomar a média dessas variâncias:
\[ \hat{\sigma}_{\text{dentro}}^2 = \overline{S^2} = \dfrac{1}{g}\sum_{i=1}^{g}S_i^2 \]
Essa estimativa não depende da veracidade de \(H_0\) e reflete a variabilidade amostral dentro de cada condição independente.
Considere agora as \(g\) médias amostrais,
\[ \bar X_1,\bar X_2,\ldots,\bar X_g \]
Como cada condição contém \(n\) unidades experimentais independentes,
\[ \operatorname{Var}(\bar X_i) = \dfrac{\sigma^2}{n} \]
Se \(H_0\) for verdadeira, todas as médias amostrais estimam a mesma média populacional \(\mu\). Assim,
\[ E(\bar X_i)=\mu \]
e
\[ \operatorname{Var}(\bar X_i)=\dfrac{\sigma^2}{n} \]
A variância amostral das \(g\) médias é
\[ S_{\bar X}^2 = \dfrac{ \sum_{i=1}^{g} \left(\bar X_i-\bar{\bar X}\right)^2 }{ g-1 } \]
onde
\[ \bar{\bar X} = \dfrac{1}{g}\sum_{i=1}^{g}\bar X_i \]
Sob \(H_0\), \(S_{\bar X}^2\) estima \(\sigma^2/n\). Multiplicando por \(n\), obtemos uma segunda estimativa de \(\sigma^2\):
\[ \hat{\sigma}_{\text{entre}}^2 = nS_{\bar X}^2 \]
Essa estimativa depende da veracidade de \(H_0\). Se as médias populacionais forem diferentes, a variabilidade observada entre as médias amostrais contém tanto a variabilidade amostral quanto a variabilidade decorrente das diferenças entre as médias populacionais. Consequentemente, \(nS_{\bar X}^2\) tende a ficar inflacionada.
A estatística de teste \(F\) compara essas duas estimativas:
\[ F = \dfrac{ \hat{\sigma}_{\text{entre}}^2 }{ \hat{\sigma}_{\text{dentro}}^2 } \]
Para um delineamento perfeitamente balanceado,
\[ F = \dfrac{ nS_{\bar X}^2 }{ \overline{S^2} } \]
ou, de forma particularmente simples,
\[ F = n\; \dfrac{ \text{variância das médias} }{ \text{média das variâncias} } \]
Portanto:
- \(\hat{\sigma}_{\text{dentro}}^2\) estima \(\sigma^2\) independentemente da veracidade de \(H_0\)
- \(\hat{\sigma}_{\text{entre}}^2\) estima \(\sigma^2\) quando \(H_0\) é verdadeira
- sob \(H_0\), as duas quantidades estimam o mesmo parâmetro \(\sigma^2\)
- quando \(H_0\) é falsa, \(\hat{\sigma}_{\text{entre}}^2\) tende a incorporar também a variabilidade entre as médias populacionais
Escrevendo diretamente em função das observações,
\[ F = n\; \dfrac{ \dfrac{ \sum_{i=1}^{g} \left(\bar X_i-\bar{\bar X}\right)^2 }{ g-1 } }{ \dfrac{ \sum_{i=1}^{g} \sum_{j=1}^{n} \left(X_{ij}-\bar X_i\right)^2 }{ g(n-1) } } \]
Sob \(H_0\) e assumindo normalidade, independência e homocedasticidade,
\[ F \sim F_{g-1,\,g(n-1)} \]
Os graus de liberdade são, portanto,
\[ \text{df}_{\text{entre}}=g-1 \]
e
\[ \text{df}_{\text{dentro}}=g(n-1)=N-g \]
Definições auxiliares:
\[ \begin{align} \bar X_i &= \dfrac{1}{n} \sum_{j=1}^{n}X_{ij} \qquad i=1,2,\ldots,g \\ \bar{\bar X} &= \dfrac{1}{g} \sum_{i=1}^{g}\bar X_i \end{align} \]
com
\[ n=n_1=n_2=\cdots=n_g \]
No delineamento perfeitamente balanceado, \(\bar{\bar X}\) coincide com a média geral das \(N=gn\) observações.
O Teorema do Limite Central não é necessário para mostrar que
\[ \operatorname{Var}(\bar X_i)=\dfrac{\sigma^2}{n} \]
Essa relação é exata para observações independentes com variância finita. O Teorema do Limite Central justifica aproximadamente a normalidade das médias amostrais quando \(n\) é suficientemente grande. Já a distribuição \(F\) exata da ANOVA clássica decorre da hipótese de normalidade da variável dependente em cada condição.
Aplicação em Nutricao3.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.
Dados <- readRDS("Nutricao3.rds")
print(FSA::Summarize(Sodium ~ Instructor,
data=Dados), digits=2, row.names=FALSE) Instructor n mean sd min Q1 median Q3 max
Brendon 20 1288 194 950 1150 1300 1400 1700
McGuirk 20 1246 142 1000 1144 1212 1350 1525
Melissa 20 1124 143 900 1006 1112 1231 1400Anova 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.00624
ANOVA results using Sodium as the dependent variable
Predictor SS df MS F p partial_eta2 CI_95_partial_eta2
Instructor 290145.83 2 145072.92 5.56 .006 .16 [.02, .32]
Error 1487812.50 57 26101.97 [NA, NA]
Note: Values in square brackets indicate the bounds of the 95% confidence interval for partial eta-squared
General Linear Model (OLS Regression)
Model Fit:
F(2, 57) = 5.56, p = 0.006 **
R² = 0.16319 (Adjusted R² = 0.13383)
Unstandardized Coefficients:
Outcome Variable: Sodium
N = 60
──────────────────────────────────────────────────────────────────────────────
b S.E. t p [95% CI of b] VIF
──────────────────────────────────────────────────────────────────────────────
(Intercept) 1287.500 (36.126) 35.639 <.001 *** [1215.159, 1359.841]
InstructorMcGuirk -41.250 (51.090) -0.807 .423 [-143.556, 61.056]
InstructorMelissa -163.750 (51.090) -3.205 .002 ** [-266.056, -61.444]
──────────────────────────────────────────────────────────────────────────────Registered S3 methods overwritten by 'MuMIn':
method from
nobs.multinom broom
nobs.fitdistr broomStandardized Coefficients (β):
Outcome Variable: Sodium
N = 60
──────────────────────────────────────────────────────────────────────────────────────
β S.E. t p [95% CI of β] r(partial) r(part)
──────────────────────────────────────────────────────────────────────────────────────
InstructorMcGuirk -0.113 (0.140) -0.807 .423 [-0.393, 0.167] -0.106 -0.098
InstructorMelissa -0.448 (0.140) -3.205 .002 ** [-0.729, -0.168] -0.391 -0.388
──────────────────────────────────────────────────────────────────────────────────────
NULLTeste t de Student é 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.} \]
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.
- Eden, T & Fisher, RA (1929) Studies in Crop Variation: VI. Experiments on the Response of the Potato to Potash and Nitrogen. Journal of Agricultural Science 19: 201–13.
- 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.
- Schmider, E et al. (2010) Is it really robust? Reinvestigating the robustness of ANOVA against violations of the normal distribution assumption. Methodology: European Journal of Research Methods for the Behavioral and Social Sciences 6(4): 147–51. https://doi.org/10.1027/1614-2241/a000016
- 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.