Exemplos Resolvidos - Ex. III
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 3 Velocidade de transporte mucociliar relativa (m/s) - Singer & Andrade (2000)
- Banco de dados no formato long table:
- A variável Concentração é qualitativa ordinal (codificada em 0,1,8,16,32,64);
- A variável tempo (5,10,15,20,25,30) é quantitativa discreta;
- A variável Palato é 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 & concentração, exceto para concentração igual a Zero. Sendo assim, pode-se fazer a primeira análise com filtro Concentração > 0, para proceder com Anova de medidas repetidas.
temp %>%
group_by (Concentracao,tempo) %>%
summarise(n=n()) %>%
ungroup() #volta a ser dataframe## `summarise()` regrouping output by 'Concentracao' (override with `.groups` argument)## # A tibble: 41 x 3
## Concentracao tempo n
## <int> <dbl> <int>
## 1 0 5 6
## 2 0 10 6
## 3 0 15 6
## 4 0 20 6
## 5 0 25 6
## 6 0 30 6
## 7 1 5 10
## 8 1 10 10
## 9 1 15 10
## 10 1 20 10
## # ... with 31 more rows- Gráfico de perfis
ggplot(temp, aes(x=tempo, valores,color = as.factor(Palato)))+
geom_point()+
geom_line() +
facet_wrap(~Concentracao) +
theme(legend.position = "none")
- Construção do Boxplot
temp = temp %>%
filter(Concentracao > 0)
temp %>%
ggplot(aes(x=as.factor(tempo),y=valores)) +
geom_boxplot()+
facet_wrap(vars(Concentracao),nrow=3)
- Gráfico de Interação
interaction.plot(x.factor = temp$tempo,
trace.factor = temp$Concentracao,
response = temp$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 = temp %>%
mutate(tempo = as.factor(tempo),
Concentracao = as.factor(Concentracao),
Palato = as.factor(Palato))
require(afex)
fit = aov_ez("Palato","valores",temp,between=c("Concentracao"),within=c("tempo"))## Contrasts set to contr.sum for the following variables: Concentracaosummary(fit)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 267.217 1 23.7155 45 507.0416 < 2.2e-16 ***
## Concentracao 25.428 4 23.7155 45 12.0625 9.594e-07 ***
## tempo 1.173 6 6.3452 270 8.3180 2.791e-08 ***
## Concentracao:tempo 2.459 24 6.3452 270 4.3594 6.964e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## tempo 0.013114 1.2258e-28
## Concentracao:tempo 0.013114 1.2258e-28
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## tempo 0.35354 0.0003581 ***
## Concentracao:tempo 0.35354 0.0001209 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## tempo 0.3714965 2.739429e-04
## Concentracao:tempo 0.3714965 8.579929e-05- Análise de resíduos:
res=residuals(fit)
hist(res)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.97922, p-value = 6.102e-05qqnorm(res)
qqline(res)
plot(fitted(fit),res)
Problemas, alternativas & desafios (em aberto):
- Transformação de Box-Cox para medidas repetidas:
- artigo sobre Box Cox (ainda pendente implementar no R).
- Uma alternativa simplificada com transformação de Box Cox:
- Ajustar anova sem medidas repetidas
fit_anova = aov(valores ~ Concentracao + tempo,data=temp)
a=boxCox(fit_anova)
maximo=cbind(a$x,a$y)
datatable(maximo)maximo[which.max(maximo[,2])]## [1] 0.5454545#0.5454545 é o valor de lambda- Agora voltamos a Anova com medidas repetidas, só que com a variável transformada:
lambda = 0.5454545
temp = temp %>%
mutate(valores_t = (valores^lambda-1)/lambda)
fit2 = aov_ez("Palato","valores_t",temp,between=c("Concentracao"),within=c("tempo"))## Contrasts set to contr.sum for the following variables: Concentracaosummary(fit2)##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 11.2154 1 26.7332 45 18.8789 7.844e-05 ***
## Concentracao 29.7166 4 26.7332 45 12.5055 6.382e-07 ***
## tempo 1.6329 6 7.2539 270 10.1297 4.194e-10 ***
## Concentracao:tempo 3.2391 24 7.2539 270 5.0235 7.318e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## tempo 0.015492 2.9839e-27
## Concentracao:tempo 0.015492 2.9839e-27
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## tempo 0.34506 8.616e-05 ***
## Concentracao:tempo 0.34506 2.730e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## tempo 0.3619767 6.247077e-05
## Concentracao:tempo 0.3619767 1.832749e-05- Análise de resíduos já na abordagem de dados transformados
res=residuals(fit2)
hist(res)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.99751, p-value = 0.8809qqnorm(res)
qqline(res)
plot(fitted(fit2),res)
- 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(Palato,Concentracao,tempo)
tabela = temp %>%
mutate(preditos_t = fitted(fit2),
preditos = (fitted(fit2)*lambda +1)^(1/lambda),
residuos = res)
datatable(tabela)- Por fim, as comparações de Tukey:
comp = lsmeans(fit2,specs = c("Concentracao","tempo"))
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%")- Conclusões:
- A metodologia de Anova com medidas repetidas aplicada aos dados de velocidade de transporte mucociliar relativa apresentou resultados satisfatórios. A transformação de Box - Cox na variável resposta que são os valores de velocidade propiciou a análise de resíduos da Anova serem normais e com variância constante. Com respeito a esfericidade de Mauchly, esta não foi satisfeita sendo assim aplicou-se as correções de Greenhouse-Geisser e Huynh-Feldt. Concluímos que a velocidade varia de acordo com o tempo e concentração, sendo significativa a interação entre tempo e concentração. O teste de comparações múltiplas de Tukey mostra as comparações que são significativas estatisticamente.