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 2 Tempo gasto para calçar meia(s) - Elias & Osake (1998)

  • Banco de dados:
    • A variável Sexo é qualitativa nominal;
    • A variável desempenho (antes ou depois) é quantitativa contínua - sendo apropriado disposição em longtable;
    • A variável Avaliação é qualitativa nominal;
    • É necessário criar uma coluna de indivíduo para auxiliar na análise dos dados;
    • Devido aos dados serem desbalanceados com relação ao sexo, e também dados omissos, é aconselhável utilizar o modelo de efeitos mistos ao invés de Anova com medidas repetidas.

  • Perfis do Desempenho (em segundos) por avaliação (antes ou depois)
ggplot(dados_ex2_long, aes (x = as.numeric(as.factor(Avaliacao)), y = Desempenho, color = as.factor(Ind))) +
  geom_point() +
  geom_line() +
  facet_wrap(~Sexo) +
  theme(legend.position = "none") 
## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 6 row(s) containing missing values (geom_path).

  • BoxPlot
dados_ex2_long %>% 
  ggplot(aes(x=Avaliacao,y=Desempenho,fill=Sexo)) +
         geom_boxplot()
## Warning: Removed 6 rows containing non-finite values (stat_boxplot).

Gráfico de interação:

  • Para construir este gráfico, é necessário considerar apenas as linhas em que todas as colunas estão preenchidas (função complete.cases);
  • Vale ressaltar a presença de dados omissos prejudica análise (pois alguns indivíduos são medidos apenas em um momento de avaliação);
  • O gráfico sugere que não existe interação entre Sexo e Avaliação; os segmentos de reta são quase paralelos, sugerindo que o Sexo age de forma semelhante para cada avaliação;
temp = dados_ex2_long[complete.cases(dados_ex2_long),]
datatable(temp)
interaction.plot(x.factor     = temp$Avaliacao,
                 trace.factor = temp$Sexo,
                 response     = temp$Desempenho,
                 fun = mean,
                 type="b",
                 col=c("red","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="Desempenho",
                 xlab="Avaliacao")

  • Modelo Linear de efeitos mistos

    • Avaliação sendo variável quantitativa & Sexo sendo variável qualitativa;
    • Indivíduo sendo variável qualitativa.
temp2 = temp %>%
                 mutate(Ind = as.factor(Ind))
m2 <- lmer(Desempenho ~ Avaliacao*Sexo + (1 | Ind), temp2)
m2
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: Desempenho ~ Avaliacao * Sexo + (1 | Ind)
##    Data: temp2
## REML criterion at convergence: 291.1709
## Random effects:
##  Groups   Name        Std.Dev.
##  Ind      (Intercept) 1.604   
##  Residual             1.338   
## Number of obs: 72, groups:  Ind, 39
## Fixed Effects:
##           (Intercept)        AvaliacaoDepois                  SexoM  
##                6.9386                -0.9076                 0.8204  
## AvaliacaoDepois:SexoM  
##               -0.1572

Valor-p do modelo misto:

  • Verificamos que não há interação significativa entre Sexo e Avaliação
  • Não há significância estatística de Avaliação e Sexo separadamente;
  • Para Obter os valores p:
    • Carrega-se o pacote lmerTest;
    • e faz-se o procedimento de ajustar novamente o modelo, mas agora incluindo a aproximação de graus de liberdade e valores-p com a opção REML = FALSE;
  • Saiba mais
# re-fit model
m3 <- lmer(Desempenho ~ Avaliacao*Sexo + (1 | Ind), temp2, REML = FALSE)

# get Satterthwaite-approximated degrees of freedom
a <- coef(summary(m3))[, 3]
# get approximate p-values
b <- coef(summary(m3))[, 5]
resultados=cbind(a,b)
resultados
##                              a            b
## (Intercept)           59.26664 1.533869e-24
## AvaliacaoDepois       35.36221 2.125727e-02
## SexoM                 60.78378 2.747117e-01
## AvaliacaoDepois:SexoM 35.43080 8.194470e-01
resultados=data.frame(resultados)
names(resultados)=c("coeficiente","valor-p")
resultados
##                       coeficiente      valor-p
## (Intercept)              59.26664 1.533869e-24
## AvaliacaoDepois          35.36221 2.125727e-02
## SexoM                    60.78378 2.747117e-01
## AvaliacaoDepois:SexoM    35.43080 8.194470e-01

Análise de resíduos condicionais:

  • apresentou desvio de normalidade
  • sugere-se retirar a observação que ficou fora do intervalo de 95%.
rc_resids <- compute_redres(m2)
## Loading required namespace: testthat
plot(rc_resids,main="resíduos condicionais versus índices")

plot_resqq(m2)

shapiro.test(rc_resids)
## 
##  Shapiro-Wilk normality test
## 
## data:  rc_resids
## W = 0.94311, p-value = 0.002722
  • Identificando outliers:
require(rstatix)
temp %>%
  group_by(Avaliacao,Sexo) %>%
  identify_outliers(Desempenho)
## # A tibble: 2 x 6
##   Sexo  Avaliacao   Ind Desempenho is.outlier is.extreme
##   <chr> <chr>     <dbl>      <dbl> <lgl>      <lgl>     
## 1 F     Antes        39       15.8 TRUE       TRUE      
## 2 M     Depois        6       10.4 TRUE       FALSE

Alternativas & desafios (em aberto):

  • Retirar o valor de “outlier” extremo;
  • Imputação de dados para completar os indivíduos que não tiveram valores duas vezes;
  • Retirar os indivíduos que tiveram somente uma medida.