Exemplos Resolvidos - Ex. IV
Lia Hanna Martins Morita
2020
require(DT) #para usar o datatable
require(tidyverse) #entre outros, a função pivot_longer
require(dplyr)
require(lme4) #efeitos mistos
require(jtools)
require(car)
require(redres)
require(lsmeans)
require(lmerTest)
require(afex)
require(car)
require(tidyverse)
require(rstatix)
require(gridExtra)
require(DT)
require(lsmeans)Exemplo 4 Dados de pH da placa bacteriana dentária - Grande, Oliveira, Singer, Santos & Nicolau (1998)
- Banco de dados no formato long table:
- A variável Periodo é qualitativa nominal (Ante & Depois);
- A variável tempo (0,5,10,15,30,45,60) é quantitativa discreta;
- A variável Voluntário é qualitativa nominal (codificada em 1,2,3,…)
- Disposição do banco de dados em longtable para análise;
O banco de dados é balanceado com respeito ao tempo & Periodo:
dados_ex4 %>%
group_by (Periodo,tempo) %>%
summarise(n=n()) %>%
ungroup() #volta a ser dataframe## `summarise()` regrouping output by 'Periodo' (override with `.groups` argument)## # A tibble: 14 x 3
## Periodo tempo n
## <chr> <dbl> <int>
## 1 Antes 0 21
## 2 Antes 5 21
## 3 Antes 10 21
## 4 Antes 15 21
## 5 Antes 30 21
## 6 Antes 45 21
## 7 Antes 60 21
## 8 Depois 0 21
## 9 Depois 5 21
## 10 Depois 10 21
## 11 Depois 15 21
## 12 Depois 30 21
## 13 Depois 45 21
## 14 Depois 60 21- Gráfico de perfis
ggplot(dados_ex4, aes(x=tempo, valores,color = as.factor(Voluntario)))+
geom_point()+
geom_line() +
facet_wrap(~Periodo) +
theme(legend.position = "none")
- Construção do Boxplot
dados_ex4 %>%
ggplot(aes(x=Periodo,y=valores,color=as.factor(tempo))) +
geom_boxplot()
- Gráfico de Interação
interaction.plot(x.factor = dados_ex4$tempo,
trace.factor = dados_ex4$Periodo,
response = dados_ex4$valores,
fun = mean,
type="b",
col=c("black","red","blue","orange","green"), ### 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 I - Anova com medidas repetidas:
- Nesta abordagem, o tempo é tratado como fator.
temp = dados_ex4 %>%
mutate(tempo = as.factor(tempo))
require(afex)
fit = aov_ez("Voluntario","valores",temp,within=c("Periodo","tempo"))
summary(fit)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 10053.2 1 6.9676 20 28856.993 < 2.2e-16 ***
## Periodo 4.3 1 7.2264 20 11.931 0.002508 **
## tempo 34.0 6 2.1238 120 320.197 < 2.2e-16 ***
## Periodo:tempo 0.9 6 1.1755 120 16.118 1.555e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## tempo 0.0018916 0.0000e+00
## Periodo:tempo 0.0111472 5.9872e-09
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## tempo 0.29699 < 2.2e-16 ***
## Periodo:tempo 0.33768 6.535e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## tempo 0.3240411 1.069897e-24
## Periodo:tempo 0.3759892 2.334663e-06- Análise de resíduos:
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.98989, p-value = 0.03935qqnorm(res)
qqline(res)
plot(fitted(fit),res)## Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data.
## fitted(..., append = TRUE) will return data and fitted values.
- Uma alternativa para o problema com transformação de Box Cox:
- Ajustar anova sem medidas repetidas;
- Encontrar a transformação ótima para a anova ajustada;
- A função boxCox do pacote car não é implementada para o anova com medidas repetidas - mas a transformação indicada já auxiliará a resover o problema de normalidade que tivemos anteriormente.
fit_anova = aov(valores ~ Periodo + tempo,data=temp)
a=boxCox(fit_anova)
maximo=cbind(a$x,a$y)
datatable(maximo)maximo[which.max(maximo[,2])]## [1] 2#2 é o valor de lambda- Agora voltamos a Anova com medidas repetidas, só que com a variável transformada:
lambda = 2
temp = temp %>%
mutate(valores_t = (valores^lambda-1)/lambda)
fit2 = aov_ez("Voluntario","valores_t",temp,within=c("Periodo","tempo"))
summary(fit2)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 81933 1 229.893 20 7127.949 < 2.2e-16 ***
## Periodo 144 1 244.974 20 11.753 0.002662 **
## tempo 1145 6 65.890 120 347.443 < 2.2e-16 ***
## Periodo:tempo 26 6 35.388 120 14.632 1.742e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## tempo 0.003244 1.0000e-12
## Periodo:tempo 0.019601 2.5869e-07
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## tempo 0.32135 < 2.2e-16 ***
## Periodo:tempo 0.37474 6.172e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## tempo 0.3549695 1.499429e-27
## Periodo:tempo 0.4245376 1.829330e-06- Análise de resíduos:
res=residuals(fit2)## 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.99079, p-value = 0.062qqnorm(res)
qqline(res)
plot(fitted(fit2),res)## Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data.
## fitted(..., append = TRUE) will return data and fitted values.
- Um conjunto de dados (dataframe) com os valores preditos
- A transformação na variável resposta permite calcular valores preditos: basta utilizar a transformação inversa no final da análise
- Os intervalos de confiança para os preditos são possíveis porém mais complicados de formalizar - envolve transformação de variáveis & método do Jacobiano (derivadas parciais na análise);
- Ao colocar os preditos & resíduos na tabela, deve-se verificar a relação de correspondência entre o banco de dados original, e os preditos e resíduos do objeto fit2 (modelo ajustado pela Anova).
temp = temp %>%
arrange(Voluntario,Periodo,tempo)
tabela = temp %>%
mutate(preditos_t = fitted.values(fit2),
preditos = (fitted.values(fit2)*lambda +1)^(1/lambda),
residuos = res)## Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data.
## fitted(..., append = TRUE) will return data and fitted values.
## Data was changed during ANOVA calculation. Thus, fitted values cannot be added to original data.
## fitted(..., append = TRUE) will return data and fitted values.datatable(tabela)- Por fim, as comparações de Tukey:
comp = lsmeans(fit2,specs = c("tempo","Periodo"))
contrastes=contrast(comp,method="pairwise")
datatable(data.frame(contrastes),list(pageLength = 5),cap="comparações múltiplas de Tukey")final = contrastes %>%
data.frame() %>%
filter(p.value <= 0.05) %>%
arrange(desc(contrast))
datatable(final, list(pageLength = 5),cap="comparações múltiplas - significativas ao nível 5%")