Anova com medidas repetidas
Lia Hanna Martins Morita
2020
Materiais didáticos & exemplos
Conjuntos de dados
Pacotes
- Pacote AFEX:Analysis of Factorial Experiments;
- Pacote car: Companion to Applied Regression;
- Pacote psycho: ferramentas para área de psicologia blog & github do projeto ;
- Pacote lsmeans: Para testes post-hoc - assim como o TukeyHSD() do pacote stats, só que aplica-se para ANOVA com medidas repetidas.
require(afex)
require(car)
require(tidyverse)
require(rstatix)
require(gridExtra)
require(DT)
require(psycho)
require(lsmeans)1 Exemplo 1: dados de atividade aerobica (fictícios)
https://statistics.laerd.com/statistical-guides/sphericity-statistical-guide.php - é um exemplo de Anova com medidas repetidas com somente um fator intra individuos, ou em inglês, one way within ANOVA; - Medidas de capacidade aeróbica (ml/min/kg) em três instantes (T1, T2 & T3) em 6 indivíduos.
aerobic=read.table("exemplo1_anova_repetidas.txt",sep=";",header=TRUE)
aerobic## id t1 t2 t3
## 1 1 45 50 55
## 2 2 42 42 45
## 3 3 36 41 43
## 4 4 39 35 40
## 5 5 51 55 59
## 6 6 44 49 56Descritivas pacotes: tidyverse, rstatix & DT
require(tidyverse)
require(rstatix)
require(DT)descritivas=aerobic %>%
get_summary_stats(t1,t2,t3, type = "mean_sd")
datatable(descritivas)- colocar a numeração dos indivíduos como um fator, para não haver confusão nos cálculos da ANOVA.
temp=aerobic %>%
mutate(id=as.factor(id)) %>%
pivot_longer(!id,names_to="tempo",values_to="valores")
datatable(temp)temp %>%
ggplot(aes(x=tempo,y=valores)) +
geom_boxplot()
Gráfico de perfis
- Perfis dos individuos
- foi utilizada as.numeric(as.factor()) para ser possível traçar os segmentos de reta nos intervalos de tempo.
ggplot(temp,aes (x = as.numeric(as.factor(tempo)), y = valores, color = id)) +
geom_point() +
geom_line() +
theme(legend.position = "none")
Anova com medidas repetidas, juntamente com teste de esfericidade de Mauchly.
- As correções de Greenhouse-Geisser and Huynh-Feldt servem para contornar problemas e exibe um teste com correção, no caso em que a suposição de esfericidade é violada;
- Conclusões: os valores de atividade aeróbica variam significativamente entre as avaliações (valor-p= 0.001886);
- Pelo teste de esfericidade de Mauchy, não rejeitamos a hipótese nula: \(H_0:\) mesma variancia para todos os tempos e mesma correlação entre os diferentes tempos.
fit = aov_ez("id","valores",temp,within=c("tempo"))
summary(fit)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 37996 1 658.28 5 288.602 1.293e-05 ***
## tempo 143 2 57.22 10 12.534 0.001886 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## tempo 0.43353 0.18795
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## tempo 0.63838 0.008985 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## tempo 0.760165 0.005284458Outro modo: porém não exibe o resultado do teste de esfericidade de Mauchly
- procede da mesma forma que a Anova One way, acrescentando o termo Error(id/tempo): para as medidas repetidas dos indivíduos na variável “tempo”:
fit2 = aov(valores ~ tempo + Error(id/tempo), data=temp)
summary(fit2)##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 5 658.3 131.7
##
## Error: id:tempo
## Df Sum Sq Mean Sq F value Pr(>F)
## tempo 2 143.44 71.72 12.53 0.00189 **
## Residuals 10 57.22 5.72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Análise de residuos:
- Note que a forma de chamar os resíduos é diferente (não é fit$residuals):
- Os resíduos são aproximadamente normais com média igual a zero e variância constante.
- Construção de um banco de dados agregando os valores preditos e resíduos da ANOVA.
res=residuals(fit)
#hist(res)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.94586, p-value = 0.3638pred=fitted(fit)
qqnorm(res)
qqline(res)
plot(pred,res,main="Gráfico de resíduos versus predito")
pred## 1 2 3 4 5 6 7 8
## 42.83333 45.33333 49.66667 42.83333 45.33333 49.66667 42.83333 45.33333
## 9 10 11 12 13 14 15 16
## 49.66667 42.83333 45.33333 49.66667 42.83333 45.33333 49.66667 42.83333
## 17 18
## 45.33333 49.66667res## 1 2 3 4 5 6
## 2.1666667 4.6666667 5.3333333 -0.8333333 -3.3333333 -4.6666667
## 7 8 9 10 11 12
## -6.8333333 -4.3333333 -6.6666667 -3.8333333 -10.3333333 -9.6666667
## 13 14 15 16 17 18
## 8.1666667 9.6666667 9.3333333 1.1666667 3.6666667 6.3333333mean(res)## [1] 0length(res)## [1] 18length(pred)## [1] 18a=cbind(temp,pred,res)
names(a)=c("id","tempo","valores","predito","residuos")
datatable(a)- a função lsmeans, juntamente com contrastes, permite efetuar as comparações múltiplas
- Verificamos que ao nível de significância \(\alpha\)=5%:
- os valores das avaliações 1 e 2 não diferem entre si de forma significativa, e os valores da avaliação 3 difere dos demais.
comp = lsmeans(fit,specs = c("tempo"))
comp## tempo lsmean SE df lower.CL upper.CL
## t1 42.8 2.82 5.88 35.9 49.8
## t2 45.3 2.82 5.88 38.4 52.3
## t3 49.7 2.82 5.88 42.7 56.6
##
## Warning: EMMs are biased unless design is perfectly balanced
## Confidence level used: 0.95contrast(comp,method="pairwise")## contrast estimate SE df t.ratio p.value
## t1 - t2 -2.50 1.38 10 -1.810 0.2156
## t1 - t3 -6.83 1.38 10 -4.948 0.0015
## t2 - t3 -4.33 1.38 10 -3.138 0.0261
##
## P value adjustment: tukey method for comparing a family of 3 estimatesExemplo 2: dados de biomassa (dados reais)
- é um exemplo de Anova com medidas repetidas com um fator entre grupos e um fator intra individuo: one way between & one way within Anova; O artigo exibe os dados detalhados
- Artigo ;
Descrição do experimento: O experimento consistiu em testar dois diferentes meios de cultura para o crescimento aeróbico da levedura Meyerozyma guilliermondii. As condições de agitação e temperatura foram idênticas. Ambas as fermentações nos dois diferentes meios foram conduzidas em triplicata e a quantidade de biomassa foi determinada por espectrofotometria de absorção molecular (g/l).
dados_artigo=read.table("exemplo2_anova_repetidas.txt",sep=";",header=T)
datatable(dados_artigo)Disposição em dados longitudinais
organizando o conjunto de dados
- colocar o tempo como fator;
- dispor em dados longitudinais - com a função pivot_longer;
- esta disposição difere dos dados multivariados em que pode-se usar a função pivot_wider.
temp = dados_artigo %>%
mutate(Tempo=as.factor(Tempo))
names(temp)=c("tempo","1","2")
temp = temp %>%
pivot_longer(!tempo, names_to = "meios", values_to = "valores",values_drop_na = TRUE)
temp=temp[order(temp$meios),]criando uma coluna de invididuo
- Falta inserir a coluna do indivíduo (as unidades são leveduras) para prosseguir com a análise;
- 1 triplicata pra cada tratamento e são dois tratamentos => 6 leveduras sendo medidos ao longo do tempo
triplicata_grupo1 = seq(1,3,1) #são tres leveduras para cada tratamento
invididuos_grupo1 = rep(triplicata_grupo1,7) #são 7 tempos
triplicata_grupo2 = seq(4,6,1) #a mesma levedura não fica nos dois tratamentos
individuos_grupo2 = rep(triplicata_grupo2,7)
individuos = as.factor(c(invididuos_grupo1,individuos_grupo2))
temp = cbind(individuos,temp)
temp## individuos tempo meios valores
## 1 1 0 1 0.04531428
## 2 2 0 1 0.05150598
## 3 3 0 1 0.05782735
## 4 1 6 1 0.25529718
## 5 2 6 1 0.24846289
## 6 3 6 1 0.24996831
## 7 1 24 1 4.63745973
## 8 2 24 1 4.66825945
## 9 3 24 1 4.69937670
## 10 1 30 1 5.12132367
## 11 2 30 1 5.18323901
## 12 3 30 1 5.15220196
## 13 1 48 1 6.11829714
## 14 2 48 1 6.04449834
## 15 3 48 1 6.08130752
## 16 1 72 1 7.30006064
## 17 2 72 1 7.38729706
## 18 3 72 1 7.34355278
## 19 1 96 1 8.88692174
## 20 2 96 1 8.91372769
## 21 3 96 1 8.86021068
## 22 4 0 2 0.03073367
## 23 5 0 2 0.03663051
## 24 6 0 2 0.04264485
## 25 4 6 2 0.19056737
## 26 5 6 2 0.19619761
## 27 6 6 2 0.19746317
## 28 4 24 2 5.36119470
## 29 5 24 2 5.40039456
## 30 6 24 2 5.44011022
## 31 4 30 2 6.08792411
## 32 5 30 2 6.16065624
## 33 6 30 2 6.12418064
## 34 4 48 2 8.09991714
## 35 5 48 2 8.19697142
## 36 6 48 2 8.14828824
## 37 4 72 2 9.18715857
## 38 5 72 2 9.29936683
## 39 6 72 2 9.24305413
## 40 4 96 2 9.38463177
## 41 5 96 2 9.41326887
## 42 6 96 2 9.35610294ordenando por individuo
- Para melhor visualização dos dados.
temp=temp[order(temp$individuos),]
datatable(temp)Descritivas
- Criando um banco de dados com as descritivas:
temp_res=temp %>%
group_by(tempo,meios) %>%
get_summary_stats(valores, type = "mean_sd")
temp_res = temp_res[order(temp_res$tempo),]
datatable(temp_res)- Vizualização dos dados com a construção do box plot;
- são somente 3 unidades por tempo & tratamento, então há pouca variabilidade nos boxplots;
temp %>%
ggplot(aes(x=tempo,y=valores,color=meios)) +
geom_boxplot()
Gráfico de perfis
- Há pouca variabilidade nos perfis;
- as.numeric(tempo) para ligar os “pontinhos” dos perfis.
plot1 <- temp %>% filter(meios == "1") %>%
ggplot(aes (x = as.numeric(tempo), y = valores, color = individuos))+
geom_point() +
geom_line() +
ggtitle("Meio de cultura 1") +
theme(legend.position = "none")
plot2 <- temp %>% filter(meios == "2") %>%
ggplot(aes (x = as.numeric(tempo), y = valores, color = individuos))+
geom_point() +
geom_line() +
ggtitle("Meio de cultura 2") +
theme(legend.position = "none")
grid.arrange(plot1,plot2, nrow =2)
Gráfico de interação: mostra a relação entre as duas variáveis na interação.
O gráfico sugere que existe interação entre tempo e meio de cultura os segmentos de reta não são paralelos, sugerindo que a cultura age de forma diferente para cada intervalo de tempo; o próximo passo é aplicar a Análise de Variância (ANOVA) com medidas repetidas.
interaction.plot(x.factor = temp$tempo,
trace.factor = temp$meios,
response = temp$valores,
fun = mean,
type="b",
col=c("black","red"), ### Colors for levels of trace var.
pch=c(19, 17), ### Symbols for levels of trace var.
fixed=TRUE, ### Order by factor order in data
leg.bty = "o",
ylab="Valores",
xlab="Tempo")
Análise de Medidas repetidas com a função aov_ez do pacote afex
- Fator entre grupos é o meio de cultura;
- Fator intraindividuo é o tempo.
- Embora a tabela ANOVA nos mostre os resultados do tratamento (meios) e tempo, além da interação entre eles, houve problemas na estimação;
- Não exibe o teste e esfericidade de Mauchly e nem o fator de correção.
- Deve-se pensar em alternativas: Transformação de dados ou métodos não paramétricos.
fit <- aov_ez("individuos","valores",temp,between=c("meios"),within=c("tempo"))## Converting to factor: meios## Contrasts set to contr.sum for the following variables: meiossummary(fit)## Warning in summary.Anova.mlm(object$Anova, multivariate = FALSE): Singular error SSP matrix:
## non-sphericity test and corrections not available##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1079.24 1 0.010676 4 404373.5 3.669e-11 ***
## meios 7.97 1 0.010676 4 2984.8 6.720e-07 ***
## tempo 481.08 6 0.019718 24 97590.8 < 2.2e-16 ***
## meios:tempo 6.45 6 0.019718 24 1309.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Análise de resíduos
- Desvio de normalidade e resíduos não têm variância constante;
res=residuals(fit)
hist(res)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.93583, p-value = 0.02048pred=fitted(fit)
plot(pred,res,main="resíduos versus predito")
qqnorm(res)
qqline(res)
Exemplo 3: dados de emotions (dados reais)
- diretamente do pacote psycho Descrição do experimento: Vários participantes são expostos a imagens com fotos neutras ou negativas, foram 48 ensaios por participante. Durante cada ensaio, o participante tinha que dar uma nota para sua valência emocional (Subjective_Valence: pontuação de -100 a 100). Além disso, 20 minutos após a tarefa, o participante tinha que dizer se lembrava da foto.
data(emotion)
temp = emotion %>%
select(Participant_ID,
Participant_Sex,
Emotion_Condition,
Subjective_Valence,
Recall)
summary(temp)## Participant_ID Participant_Sex Emotion_Condition Subjective_Valence
## 10S : 48 Female:720 Negative:456 Min. :-100.000
## 11S : 48 Male :192 Neutral :456 1st Qu.: -65.104
## 12S : 48 Median : -2.604
## 13S : 48 Mean : -18.900
## 14S : 48 3rd Qu.: 7.000
## 15S : 48 Max. : 100.000
## (Other):624
## Recall
## Mode :logical
## FALSE:600
## TRUE :312
##
##
##
## datatable(temp)Descritivas
- Vamos considerar as variáveis “Emotion_Condition” e “Participant_Sex”;
- Sendo assim, temos uma ANOVA com medidas repetidas, sendo um fator entre grupos (Sexo), e um fator intra indivíduos (Emotion_Condition).
temp_res=temp %>%
group_by(Participant_Sex,Emotion_Condition) %>%
get_summary_stats(Subjective_Valence, type = "mean_sd")
datatable(temp_res)- Exibe os boxplots por Sexo e por tipo de foto: negativa ou neutra;
temp %>%
ggplot(aes(x=Participant_Sex,y=Subjective_Valence,color=Emotion_Condition)) +
geom_boxplot()
Gráfico de perfis
- Perfis dos individuos
plot1 <- temp %>% filter(Participant_Sex == "Female") %>%
ggplot(aes (x = as.numeric(as.factor(Emotion_Condition)), y = Subjective_Valence, color = Participant_ID))+
geom_point() +
geom_line() +
theme(legend.position = "none") +
ggtitle("Feminino")
plot2 <- temp %>% filter(Participant_Sex == "Male") %>%
ggplot(aes (x = as.numeric(as.factor(Emotion_Condition)), y = Subjective_Valence, color = Participant_ID))+
geom_point() +
geom_line() +
theme(legend.position = "none") +
ggtitle("Masculino")
grid.arrange(plot1,plot2, nrow =2)
Gráfico de interação: Mostra a relação entre as duas variáveis na interação.
O gráfico sugere que existe interação entre tipo de foto e Sexo do participante - os segmentos de reta não são paralelos, sugerindo que a característica da foto age de forma diferente para cada sexo. o próximo passo é aplicar a Análise de Variância (ANOVA) com medidas repetidas.
interaction.plot(x.factor = temp$Emotion_Condition,
trace.factor = temp$Participant_Sex,
response = temp$Subjective_Valence,
fun = mean,
type="b",
col=c("black","red"), ### Colors for levels of trace var.
pch=c(19, 17), ### Symbols for levels of trace var.
fixed=TRUE, ### Order by factor order in data
leg.bty = "o",
ylab="Subjective_Valence",
xlab="Emotion")
fit = aov_ez("Participant_ID","Subjective_Valence",temp,between=c("Participant_Sex"),within=c("Emotion_Condition"))## Warning: More than one observation per cell, aggregating the data using mean
## (i.e, fun_aggregate = mean)!## Contrasts set to contr.sum for the following variables: Participant_Sexsummary(fit)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value
## (Intercept) 7829.3 1 4684.9 17 28.4101
## Participant_Sex 126.6 1 4684.9 17 0.4593
## Emotion_Condition 29331.8 1 3044.9 17 163.7645
## Participant_Sex:Emotion_Condition 854.0 1 3044.9 17 4.7680
## Pr(>F)
## (Intercept) 5.525e-05 ***
## Participant_Sex 0.5071
## Emotion_Condition 3.739e-10 ***
## Participant_Sex:Emotion_Condition 0.0433 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Outro modo: não exibe o teste de esfericidade de Mauchly
fit2 = aov(Subjective_Valence ~ Emotion_Condition*Participant_Sex+ Error(Participant_ID/Emotion_Condition), data=temp)
summary(fit2)##
## Error: Participant_ID
## Df Sum Sq Mean Sq F value Pr(>F)
## Participant_Sex 1 3038 3038 0.459 0.507
## Residuals 17 112437 6614
##
## Error: Participant_ID:Emotion_Condition
## Df Sum Sq Mean Sq F value Pr(>F)
## Emotion_Condition 1 1278417 1278417 297.401 3.32e-12 ***
## Emotion_Condition:Participant_Sex 1 20496 20496 4.768 0.0433 *
## Residuals 17 73077 4299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 874 935646 1070Verificação dos resíduos
- Verificamos desvios de normalidade
res=residuals(fit)## Data was changed during ANOVA calculation. Thus, residuals cannot be added to original data.
## residuals(..., append = TRUE) will return data and residuals.hist(res)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.93612, p-value = 0.03139pred=fitted(fit)## Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data.
## fitted(..., append = TRUE) will return data and fitted values.qqnorm(res)
qqline(res)
Banco de dados agregando os resíduos e preditos
a=cbind(temp,pred,res)## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discardednames(a)=c("ID","Sex","Emotion_Condition","Subjective_Valence","Recall","Predito","Residuos")
datatable(a)Comparações Múltiplas com teste de Tukey
- Há problema pois o experimento não é perfeitamente balanceado;
- Alternativa: Modelo de regressão com efeito aleatório do indivíduo;
- pacote lme4: ajuste do modelo linear com medidas repetidas;
comp = lsmeans(fit,specs = c("Participant_Sex","Emotion_Condition"))
comp## Participant_Sex Emotion_Condition lsmean SE df lower.CL upper.CL
## Female Negative -61.0 4.91 27.1 -71.10 -50.9
## Male Negative -44.9 6.37 33.6 -57.87 -32.0
## Female Neutral 18.8 4.91 27.1 8.67 28.8
## Male Neutral 11.6 6.37 33.6 -1.35 24.5
##
## Warning: EMMs are biased unless design is perfectly balanced
## Confidence level used: 0.95contrast(comp,method="pairwise")## contrast estimate SE df t.ratio p.value
## Female Negative - Male Negative -16.10 8.48 32.5 -1.898 0.2487
## Female Negative - Female Neutral -79.78 4.89 17.0 -16.325 <.0001
## Female Negative - Male Neutral -72.62 8.48 32.5 -8.559 <.0001
## Male Negative - Female Neutral -63.67 8.48 32.5 -7.504 <.0001
## Male Negative - Male Neutral -56.52 9.46 17.0 -5.973 0.0001
## Female Neutral - Male Neutral 7.15 8.48 32.5 0.843 0.8336
##
## P value adjustment: tukey method for comparing a family of 4 estimatesVerifica-se mais pessoas no sexo feminino
temp %>%
distinct(Participant_ID,Participant_Sex)## Participant_ID Participant_Sex
## 1 1S Female
## 2 2S Female
## 3 3S Female
## 4 4S Female
## 5 5S Female
## 6 6S Female
## 7 7S Female
## 8 8S Female
## 9 9S Female
## 10 10S Female
## 11 11S Female
## 12 12S Female
## 13 13S Female
## 14 14S Female
## 15 15S Female
## 16 16S Male
## 17 17S Male
## 18 18S Male
## 19 19S Male2 Exemplo 1: dados de atividade aerobica (fictícios)
https://statistics.laerd.com/statistical-guides/sphericity-statistical-guide.php - é um exemplo de Anova com medidas repetidas com somente um fator intra individuos, ou em inglês, one way within ANOVA; - Medidas de capacidade aeróbica (ml/min/kg) em três instantes (T1, T2 & T3) em 6 indivíduos.
aerobic=read.table("exemplo1_anova_repetidas.txt",sep=";",header=TRUE)
aerobic## id t1 t2 t3
## 1 1 45 50 55
## 2 2 42 42 45
## 3 3 36 41 43
## 4 4 39 35 40
## 5 5 51 55 59
## 6 6 44 49 56Descritivas
descritivas=aerobic %>%
get_summary_stats(t1,t2,t3, type = "mean_sd")
datatable(descritivas)- colocar a numeração dos indivíduos como um fator, para não haver confusão nos cálculos da ANOVA.
temp=aerobic %>%
mutate(id=as.factor(id)) %>%
pivot_longer(!id,names_to="tempo",values_to="valores")
datatable(temp)temp %>%
ggplot(aes(x=tempo,y=valores)) +
geom_boxplot()
Gráfico de perfis
- Perfis dos individuos
- foi utilizada as.numeric(as.factor()) para ser possível traçar os segmentos de reta nos intervalos de tempo.
ggplot(temp,aes (x = as.numeric(as.factor(tempo)), y = valores, color = id)) +
geom_point() +
geom_line() +
theme(legend.position = "none")
Anova com medidas repetidas, juntamente com teste de esfericidade de Mauchly.
- As correções de Greenhouse-Geisser and Huynh-Feldt servem para contornar problemas e exibe um teste com correção, no caso em que a suposição de esfericidade é violada;
- Conclusões: os valores de atividade aeróbica variam significativamente entre as avaliações (valor-p= 0.001886);
- Pelo teste de esfericidade de Mauchy, não rejeitamos a hipótese nula: \(H_0:\) mesma variancia para todos os tempos e mesma correlação entre os diferentes tempos.
fit = aov_ez("id","valores",temp,within=c("tempo"))
summary(fit)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 37996 1 658.28 5 288.602 1.293e-05 ***
## tempo 143 2 57.22 10 12.534 0.001886 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## tempo 0.43353 0.18795
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## tempo 0.63838 0.008985 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## tempo 0.760165 0.005284458Outro modo: porém não exibe o resultado do teste de esfericidade de Mauchly
- procede da mesma forma que a Anova One way, acrescentando o termo Error(id/tempo): para as medidas repetidas dos indivíduos na variável “tempo”:
fit2 = aov(valores ~ tempo + Error(id/tempo), data=temp)
summary(fit2)##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 5 658.3 131.7
##
## Error: id:tempo
## Df Sum Sq Mean Sq F value Pr(>F)
## tempo 2 143.44 71.72 12.53 0.00189 **
## Residuals 10 57.22 5.72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Análise de residuos:
- Note que a forma de chamar os resíduos é diferente (não é fit$residuals):
- Os resíduos são aproximadamente normais com média igual a zero e variância constante.
- Construção de um banco de dados agregando os valores preditos e resíduos da ANOVA.
res=residuals(fit)
#hist(res)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.94586, p-value = 0.3638pred=fitted(fit)
qqnorm(res)
qqline(res)
plot(pred,res,main="Gráfico de resíduos versus predito")
pred## 1 2 3 4 5 6 7 8
## 42.83333 45.33333 49.66667 42.83333 45.33333 49.66667 42.83333 45.33333
## 9 10 11 12 13 14 15 16
## 49.66667 42.83333 45.33333 49.66667 42.83333 45.33333 49.66667 42.83333
## 17 18
## 45.33333 49.66667res## 1 2 3 4 5 6
## 2.1666667 4.6666667 5.3333333 -0.8333333 -3.3333333 -4.6666667
## 7 8 9 10 11 12
## -6.8333333 -4.3333333 -6.6666667 -3.8333333 -10.3333333 -9.6666667
## 13 14 15 16 17 18
## 8.1666667 9.6666667 9.3333333 1.1666667 3.6666667 6.3333333mean(res)## [1] 0length(res)## [1] 18length(pred)## [1] 18a=cbind(temp,pred,res)
names(a)=c("id","tempo","valores","predito","residuos")
datatable(a)- a função lsmeans, juntamente com contrastes, permite efetuar as comparações múltiplas
- Verificamos que ao nível de significância \(\alpha\)=5%:
- os valores das avaliações 1 e 2 não diferem entre si de forma significativa, e os valores da avaliação 3 difere dos demais.
comp = lsmeans(fit,specs = c("tempo"))
comp## tempo lsmean SE df lower.CL upper.CL
## t1 42.8 2.82 5.88 35.9 49.8
## t2 45.3 2.82 5.88 38.4 52.3
## t3 49.7 2.82 5.88 42.7 56.6
##
## Warning: EMMs are biased unless design is perfectly balanced
## Confidence level used: 0.95contrast(comp,method="pairwise")## contrast estimate SE df t.ratio p.value
## t1 - t2 -2.50 1.38 10 -1.810 0.2156
## t1 - t3 -6.83 1.38 10 -4.948 0.0015
## t2 - t3 -4.33 1.38 10 -3.138 0.0261
##
## P value adjustment: tukey method for comparing a family of 3 estimates