Pacotes

  • Pacote lme4: Linear Mixed model com efeitos aleatórios

Exemplos com efeitos mistos no site

require(lme4)
require(tidyverse)
require(afex)
require(DT)
require(jtools)
require(car)
require(redres)
require(lsmeans)
  • Conjuntos de dados do pacote lme4
data(package = "lme4")

  • Exemplo 1: Dados de penicilina:

    • Seis amostras de penicilina foram testadas usando o método de placa de Bacillus subtilis em cada uma das 24 placas. A resposta é o diâmetro (mm) da zona de inibição de crescimento do organismo.
    • Os dados são descritos em Davies e Goldsmith (1972) como provenientes de uma investigação para avaliar a variabilidade entre as amostras de penicilina pelo método de B. subtilis. Neste método de teste, um meio de ágar nutriente inoculado a granel é vertido em uma placa de Petri de aproximadamente 90 mm. diâmetro, conhecido como placa. Quando o meio está curado, seis pequenos cilindros ocos ou potes (cerca de 4 mm de diâmetro) são cimentados na superfície em intervalos igualmente espaçados. Algumas gotas das soluções de penicilina a serem comparadas são colocadas nos respectivos cilindros, e toda a placa é colocada em uma incubadora por um determinado tempo. A penicilina se difunde dos potes para o ágar, e isso produz uma zona de inibição circular clara de crescimento dos organismos, que pode ser facilmente medida. O diâmetro da zona está relacionado de forma conhecida à concentração de penicilina na solução.

Acesso ao banco de dados e visualização em tabela

  • No banco de dados Penicillin,
    • A variável diameter é quantitativa;
    • As variáveis plate (indivíduos) & sample (amostras) são qualitativas.
summary(Penicillin)
##     diameter         plate     sample
##  Min.   :18.00   a      :  6   A:24  
##  1st Qu.:22.00   b      :  6   B:24  
##  Median :23.00   c      :  6   C:24  
##  Mean   :22.97   d      :  6   D:24  
##  3rd Qu.:24.00   e      :  6   E:24  
##  Max.   :27.00   f      :  6   F:24  
##                  (Other):108
datatable(Penicillin,options = list(pageLength = 5))
  • Gráficos de Boxplot dos diâmetros, por amostras coletadas.
    • Verificamos que há variação nas respostas (diâmetros) ao longo das amostras.
Penicillin  %>% 
  ggplot(aes(x=sample,y=diameter)) +
         geom_boxplot()

  • Gráfico dos perfis:
    • Para ligar os pontos por linhas é necessário escrever as.numeric()
    • Verificamos grande variabilidade nos perfis de cada unidade amostral (plate)
  ggplot(Penicillin,aes (x = as.numeric(sample), y = diameter, color = plate)) +
  geom_point() +
  geom_line() +
  theme(legend.position = "none")

Gráficos de dispersão de diameter versus sample - em facetas para cada indivíduo (plate).

  • Para traçar a linha de tendência é necessário escrever as.numeric() para os valores no eixo x, embora no nosso modelo (abaixo) será tratada como as.factor ou fator, ou variável qualitativa/classe.

  • A pergunta é:

    • cada indivíduo tem seu próprio “slope” ou inclinação, ou será que as linhas só diferem pelo intercepto?
ggplot(Penicillin, aes(x=as.numeric(sample), diameter))+
  geom_point()+
  geom_smooth(method="lm")+
  facet_wrap(~plate)
## `geom_smooth()` using formula 'y ~ x'

Modelo linear de efeitos mistos:

  • Modelo I (mais simples):
    • Efeito fixo (sample) + Efeito aleatório no intercepto (plate)
m1 <- lmer(diameter ~ sample + (1 | plate), Penicillin)
m1
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: diameter ~ sample + (1 | plate)
##    Data: Penicillin
## REML criterion at convergence: 308.2793
## Random effects:
##  Groups   Name        Std.Dev.
##  plate    (Intercept) 0.8467  
##  Residual             0.5499  
## Number of obs: 144, groups:  plate, 24
## Fixed Effects:
## (Intercept)      sampleB      sampleC      sampleD      sampleE      sampleF  
##      25.167       -3.208       -0.250       -2.292       -2.208       -5.208

  • A fim de fazer algumas comparações entre os diferentes modelos:

    • O critério de informação de Akaike (AIC) e o critério de informação Bayesiano (BIC) são medidas de qualidade relativa de um modelo estatístico para um determinado conjunto de dados. O cálculo de AIC e BIC são dados por: \[ \begin{array}{lll} \text{AIC} &=& −2\log(L)+2\text{p};\\ \text{BIC} &=& −2\log(L)+ \text{p} \log(n), \end{array} \] onde L é o valor maximizado da função de verossimilhança, p é o número de parâmetros no modelo e n é o número total de observações no banco de dados.

    • Dado um conjunto de modelos candidatos para os dados, o modelo preferido é aquele com os valores mínimos de AIC e BIC. Esses critérios não apenas recompensam a adequação, mas também incluem uma penalidade (eles se tornam maiores quando incluímos mais parâmetros no modelo).

Referências:

  • SAKAMOTO, Y.; ISHIGURO, M.; KITAGAWA, G. Akaike information criterion statistics .
  1. ed. [S.l.]: Springer Netherlands, 1986. v. 1.
  • SCHWARZ, G. Estimating the dimension of a model. The Annals of Statistics , Institute of Mathematical Statistics, v. 6, n. 2, p. 461–464, 1978.

Exibindo o AIC & BIC do modelo de efeitos mistos

AIC(m1)
## [1] 324.2793
BIC(m1)
## [1] 348.0378

Utilizando o pacote emmeans para realizar as comparações múltiplas por método de Tukey

  • emmeans: Estimated marginal means no modelo linear de efeitos mistos
    • Verificamos que as amostras diferem entre si, exceto nas comparações A-C & D-E
require(emmeans)
emmeans(m1, list(pairwise ~ sample), adjust = "tukey")
## $`emmeans of sample`
##  sample emmean    SE   df lower.CL upper.CL
##  A        25.2 0.206 39.7     24.6     25.7
##  B        22.0 0.206 39.7     21.4     22.5
##  C        24.9 0.206 39.7     24.3     25.5
##  D        22.9 0.206 39.7     22.3     23.4
##  E        23.0 0.206 39.7     22.4     23.5
##  F        20.0 0.206 39.7     19.4     20.5
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 6 estimates 
## 
## $`pairwise differences of sample`
##  contrast estimate    SE  df t.ratio p.value
##  A - B      3.2083 0.159 115  20.210 <.0001 
##  A - C      0.2500 0.159 115   1.575 0.6167 
##  A - D      2.2917 0.159 115  14.436 <.0001 
##  A - E      2.2083 0.159 115  13.911 <.0001 
##  A - F      5.2083 0.159 115  32.809 <.0001 
##  B - C     -2.9583 0.159 115 -18.635 <.0001 
##  B - D     -0.9167 0.159 115  -5.774 <.0001 
##  B - E     -1.0000 0.159 115  -6.299 <.0001 
##  B - F      2.0000 0.159 115  12.598 <.0001 
##  C - D      2.0417 0.159 115  12.861 <.0001 
##  C - E      1.9583 0.159 115  12.336 <.0001 
##  C - F      4.9583 0.159 115  31.234 <.0001 
##  D - E     -0.0833 0.159 115  -0.525 0.9951 
##  D - F      2.9167 0.159 115  18.373 <.0001 
##  E - F      3.0000 0.159 115  18.898 <.0001 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 6 estimates

Explorando o pacote redres:

  • Github ;

  • redres é um pacote R desenvolvido para ajudar no diagnóstico de modelos lineares mistos ajustados usando a função lmer do pacote lme4. Destina-se a complementar o pacote lme4.

  • calcula os resíduos de um modelo lmer. Temos os seguintes tipos de resíduos:

    • “pearson_cond”: Pearson conditional residuals ou reíduos condicionais de Pearson;
    • “pearson_mar”: Pearson marginal residuals ou resíduos marginais de Pearson;
    • “raw_cond”: raw conditional residuals ou resíduos condicionais
    • “raw_mar”: raw marginal residuals ou resíduos marginais
    • “std_cond”: studentized conditional residuals ou resíduos condicionais estudentizados;
    • “std_mar”: studentized marginal residuals ou resíduos marginais estudentizados
    • “ranef”: efeitos aleatórios dos indivíduos

Resíduos condicionais: (raw_cond)

  • Os resíduos são aproximadamente normais, independentes e com variância constante
    • O gráfico de normalidade exibe o intervalo “envelope” de 95%.
rc_resids <- compute_redres(m1)
## Loading required namespace: testthat
plot(rc_resids,main="resíduos condicionais versus índices",ylim=c(-2,2))
abline(h=2, lty=2)
abline(h=-2, lty=2)

plot_resqq(m1)

shapiro.test(rc_resids)
## 
##  Shapiro-Wilk normality test
## 
## data:  rc_resids
## W = 0.9886, p-value = 0.2878

Efeitos aleatórios dos indivíduos (ranef)

  • Os efeitos aleatórios são aproximadamente normais, independentes e com variância constante
    • O gráfico de normalidade exibe o intervalo “envelope” de 95%.
# creates normal quantile plots for each random effect
random <- ranef(m1)
aleatorio = random[["plate"]][["(Intercept)"]]
plot(aleatorio,main="efeitos aleatórios versus índices",ylim=c(-2,2))
abline(h=2, lty=3)
abline(h=-2, lty=3)
abline(h=0,lty=3,col=4)

plot_ranef(m1)

shapiro.test(aleatorio)
## 
##  Shapiro-Wilk normality test
## 
## data:  aleatorio
## W = 0.95425, p-value = 0.334
aleatorio
##  [1]  0.80454708  0.80454708  0.18167192  0.33739071  0.02595313 -0.44120324
##  [7] -1.37551598  0.80454708 -0.75264082 -0.75264082  0.96026587  0.49310950
## [13]  1.42742224  0.49310950  0.96026587  0.02595313 -0.28548445 -0.28548445
## [19] -1.37551598  0.96026587 -0.90835961 -0.28548445 -0.59692203 -1.21979719

Coloca os resíduos no mesmo banco de dados

  • os valores preditos por “sample” são os mesmos para todos os “plates”;
  • os valores observados de cada indivíduo se distanciam dos preditos, através de um resíduo e um termo aleatório;
  • pendente colocar o termo aleatório nesta tabela.
preditos=fitted(m1)
resids = data.frame(Penicillin, preditos, rc_resids)
datatable(resids,options = list(pageLength = 5))

Cria aplicativo no shiny com os resultados

#launch_redres(m1)

Modelo linear de efeitos mistos:

  • Modelo II (mais sofisticado)
    • Efeito fixo (sample) de forma quantitativa
    • Efeito aleatório no intercepto (plate) + Inclinação da Reta
    • Algoritmo de estimação deu problema quando fizemos as mudanças!
temp = Penicillin %>%
      mutate(sample= as.numeric(sample))
m12 <- lmer(diameter ~ sample + (sample | plate), temp)
## boundary (singular) fit: see ?isSingular
m12
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: diameter ~ sample + (sample | plate)
##    Data: temp
## REML criterion at convergence: 545.8233
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  plate    (Intercept) 0.58581      
##           sample      0.01151  1.00
##  Residual             1.50123      
## Number of obs: 144, groups:  plate, 24
## Fixed Effects:
## (Intercept)       sample  
##     25.4806      -0.7167  
## convergence code 0; 0 optimizer warnings; 1 lme4 warnings

Exibindo o AIC & BIC dos dois modelos

  • O AIC & BIC do segundo modelo ficaram piores!
AIC(m1,m12)
##     df      AIC
## m1   8 324.2793
## m12  6 557.8233
BIC(m1,m12)
##     df      BIC
## m1   8 348.0378
## m12  6 575.6422
  • Exemplo 2: Dados de Tempos de reação em um estudo de privação de sono:
    • O tempo médio de reação por dia para indivíduos em um estudo de privação de sono. No dia 0, os indivíduos tiveram sua quantidade normal de sono. A partir daquela noite, eles estavam restritos a 3 horas de sono por noite. As observações representam o tempo médio de reação em uma série de testes dados a cada dia para cada sujeito.

Acesso ao banco de dados e visualização em tabela

  • No banco de dados sleepstudy,
    • As variáveis Reaction & Days são quantitativas;
    • A variável Subject (indivíduos) é qualitativa.
summary(sleepstudy)
##     Reaction          Days        Subject   
##  Min.   :194.3   Min.   :0.0   308    : 10  
##  1st Qu.:255.4   1st Qu.:2.0   309    : 10  
##  Median :288.7   Median :4.5   310    : 10  
##  Mean   :298.5   Mean   :4.5   330    : 10  
##  3rd Qu.:336.8   3rd Qu.:7.0   331    : 10  
##  Max.   :466.4   Max.   :9.0   332    : 10  
##                                (Other):120
datatable(sleepstudy,options = list(pageLength = 5))
  • Gráficos de Boxplot dos diâmetros, por amostras coletadas.
    • Verificamos que os valores (de tempo de reação no teste) em geral são maiores ao longo dos dias, devido à privação de sono.
#Criando um banco de dados "temp" para não modificar o original "sleepstudy"
temp = sleepstudy %>% 
  mutate(Days = as.factor(Days))  

temp %>% 
  ggplot(aes(x=Days,y=Reaction)) +
         geom_boxplot()

  • Gráfico dos perfis:
    • Verificamos grande variabilidade nos perfis de cada unidade amostral (Subject), sendo que em geral as respostas crescem ao longo dos dias.
  ggplot(temp,aes (x = as.numeric(Days), y = Reaction, color = Subject)) +
  geom_point() +
  geom_line() +
  theme(legend.position = "none")

Gráficos de dispersão de Reaction versus Days - em facetas para cada indivíduo (Subject).

ggplot(temp, aes(x=as.numeric(Days), Reaction))+
  geom_point()+
  geom_smooth(method="lm")+
  facet_wrap(~Subject)
## `geom_smooth()` using formula 'y ~ x'

Modelo linear de efeitos mistos:

  • Modelo I (mais simples)
    • Efeito fixo (Days) + Efeito aleatório no intercepto (Subject)
m2 <- lmer(Reaction ~ Days + (1 | Subject), temp)
m2
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: temp
## REML criterion at convergence: 1729.493
## Random effects:
##  Groups   Name        Std.Dev.
##  Subject  (Intercept) 37.09   
##  Residual             31.43   
## Number of obs: 180, groups:  Subject, 18
## Fixed Effects:
## (Intercept)        Days1        Days2        Days3        Days4        Days5  
##     256.652        7.844        8.710       26.340       31.998       51.867  
##       Days6        Days7        Days8        Days9  
##      55.526       62.099       79.978       94.199

Exibindo o AIC & BIC do modelo de efeitos mistos

AIC(m2)
## [1] 1753.493
BIC(m2)
## [1] 1791.808

Comparações múltiplas pelo método de Tukey

  • Verificamos que as amostras diferem entre si, exceto nas comparações A-C & D-E
  • As seguintes comparações foram significativas ao nível de 5%:
    • Dias 0, 1 e 2: diferem dos dias 5 a 9;
    • Dia 3: difere dos dias 7 a 9;
    • Dia 4: difere dos dias 8 a 9;
    • Dias 5 e 6: diferem do dia 9;
emmeans(m2, list(pairwise ~ Days), adjust = "tukey")
## $`emmeans of Days`
##  Days emmean   SE df lower.CL upper.CL
##  0       257 11.5 42      223      291
##  1       264 11.5 42      231      298
##  2       265 11.5 42      232      299
##  3       283 11.5 42      249      317
##  4       289 11.5 42      255      323
##  5       309 11.5 42      275      342
##  6       312 11.5 42      278      346
##  7       319 11.5 42      285      353
##  8       337 11.5 42      303      370
##  9       351 11.5 42      317      385
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 10 estimates 
## 
## $`pairwise differences of Days`
##  contrast estimate   SE  df t.ratio p.value
##  0 - 1      -7.844 10.5 153 -0.749  0.9991 
##  0 - 2      -8.710 10.5 153 -0.831  0.9980 
##  0 - 3     -26.340 10.5 153 -2.515  0.2692 
##  0 - 4     -31.998 10.5 153 -3.055  0.0770 
##  0 - 5     -51.867 10.5 153 -4.951  0.0001 
##  0 - 6     -55.526 10.5 153 -5.301  <.0001 
##  0 - 7     -62.099 10.5 153 -5.928  <.0001 
##  0 - 8     -79.978 10.5 153 -7.635  <.0001 
##  0 - 9     -94.199 10.5 153 -8.993  <.0001 
##  1 - 2      -0.866 10.5 153 -0.083  1.0000 
##  1 - 3     -18.496 10.5 153 -1.766  0.7554 
##  1 - 4     -24.154 10.5 153 -2.306  0.3914 
##  1 - 5     -44.023 10.5 153 -4.203  0.0018 
##  1 - 6     -47.682 10.5 153 -4.552  0.0004 
##  1 - 7     -54.255 10.5 153 -5.179  <.0001 
##  1 - 8     -72.134 10.5 153 -6.886  <.0001 
##  1 - 9     -86.355 10.5 153 -8.244  <.0001 
##  2 - 3     -17.630 10.5 153 -1.683  0.8033 
##  2 - 4     -23.288 10.5 153 -2.223  0.4457 
##  2 - 5     -43.157 10.5 153 -4.120  0.0024 
##  2 - 6     -46.816 10.5 153 -4.469  0.0006 
##  2 - 7     -53.389 10.5 153 -5.097  <.0001 
##  2 - 8     -71.268 10.5 153 -6.803  <.0001 
##  2 - 9     -85.489 10.5 153 -8.161  <.0001 
##  3 - 4      -5.657 10.5 153 -0.540  0.9999 
##  3 - 5     -25.526 10.5 153 -2.437  0.3118 
##  3 - 6     -29.186 10.5 153 -2.786  0.1506 
##  3 - 7     -35.759 10.5 153 -3.414  0.0274 
##  3 - 8     -53.637 10.5 153 -5.120  <.0001 
##  3 - 9     -67.859 10.5 153 -6.478  <.0001 
##  4 - 5     -19.869 10.5 153 -1.897  0.6712 
##  4 - 6     -23.529 10.5 153 -2.246  0.4303 
##  4 - 7     -30.101 10.5 153 -2.874  0.1223 
##  4 - 8     -47.980 10.5 153 -4.580  0.0004 
##  4 - 9     -62.202 10.5 153 -5.938  <.0001 
##  5 - 6      -3.660 10.5 153 -0.349  1.0000 
##  5 - 7     -10.232 10.5 153 -0.977  0.9932 
##  5 - 8     -28.111 10.5 153 -2.684  0.1898 
##  5 - 9     -42.333 10.5 153 -4.041  0.0033 
##  6 - 7      -6.572 10.5 153 -0.627  0.9998 
##  6 - 8     -24.451 10.5 153 -2.334  0.3734 
##  6 - 9     -38.673 10.5 153 -3.692  0.0112 
##  7 - 8     -17.879 10.5 153 -1.707  0.7901 
##  7 - 9     -32.101 10.5 153 -3.064  0.0750 
##  8 - 9     -14.222 10.5 153 -1.358  0.9379 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 10 estimates

Explorando o pacote redres:

Resíduos condicionais: (raw_cond)

  • Os resíduos NÃO são aproximadamente normais, e NÂO são independentes e com variância constante
    • O gráfico de normalidade exibe o intervalo “envelope” de 95%.
rc_resids <- compute_redres(m2)
plot(rc_resids,main="resíduos condicionais versus índices",ylim=c(-2,2))
abline(h=2, lty=2)
abline(h=-2, lty=2)

plot_resqq(m2)

shapiro.test(rc_resids)
## 
##  Shapiro-Wilk normality test
## 
## data:  rc_resids
## W = 0.97274, p-value = 0.001345

Efeitos aleatórios dos indivíduos (ranef)

  • Os efeitos aleatórios são aproximadamente normais, independentes porém não parecem ter variância constante
    • O gráfico de normalidade exibe o intervalo “envelope” de 95%, que apesar de ter pontos fora deste intervalo, não rejeita-se a hipótese de normalidade pelo teste de Shapiro-Wilk.
# creates normal quantile plots for each random effect
random <- ranef(m2)
aleatorio = random[["Subject"]][["(Intercept)"]]
plot(aleatorio,main="efeitos aleatórios versus índices")
abline(h=55, lty=3)
abline(h=-55, lty=3)
abline(h=0,lty=3,col=4)

plot_ranef(m2)

shapiro.test(aleatorio)
## 
##  Shapiro-Wilk normality test
## 
## data:  aleatorio
## W = 0.94606, p-value = 0.3666
aleatorio
##  [1]  40.703420 -77.696293 -62.984327   4.397767  10.196077   8.205053
##  [7]  16.468010  -2.991081 -45.192981  72.040582 -21.154520  14.083583
## [13]  -7.846743  36.306808   7.022529  -6.350177  -3.287787  18.080083

Coloca os resíduos no mesmo banco de dados

  • os valores preditos por “Days” são os mesmos para todos os “Subjects”;
  • os valores observados de cada indivíduo se distanciam dos preditos, através de um resíduo e um termo aleatório;
    • pendente colocar o termo aleatório nesta tabela.
preditos=fitted(m2)
resids = data.frame(temp, preditos, rc_resids)
datatable(resids,options = list(pageLength = 5))

Cria aplicativo no shiny com os resultados

#launch_redres(m2)

Modelo linear de efeitos mistos:

  • Modelo II (mais sofisticado)
    • Efeito fixo (Days) de forma quantitativa
    • Efeito aleatório no intercepto (Subject) + Inclinação da Reta
m3 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
m3
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: Reaction ~ Days + (Days | Subject)
##    Data: sleepstudy
## REML criterion at convergence: 1743.628
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  Subject  (Intercept) 24.741       
##           Days         5.922   0.07
##  Residual             25.592       
## Number of obs: 180, groups:  Subject, 18
## Fixed Effects:
## (Intercept)         Days  
##      251.41        10.47

Exibindo o AIC & BIC dos dois modelos

  • O AIC do modelo II ficou bem próximo e o BIC ficou melhor!
AIC(m2,m3)
##    df      AIC
## m2 12 1753.493
## m3  6 1755.628
BIC(m2,m3)
##    df      BIC
## m2 12 1791.808
## m3  6 1774.786
  • O modelo II NÂO apresentou melhoras quanto a normalidade dos resíduos condicionais.
rc_resids <- compute_redres(m3)
plot(rc_resids,main="resíduos condicionais versus índices",ylim=c(-3,3))
abline(h=2, lty=2)
abline(h=-2, lty=2)

plot_resqq(m3)

shapiro.test(rc_resids)
## 
##  Shapiro-Wilk normality test
## 
## data:  rc_resids
## W = 0.90146, p-value = 1.408e-09

Efeitos aleatórios dos indivíduos (ranef)

  • O modelo II manteve a normalidade dos efeitos aleatórios dos indivíduos!
# creates normal quantile plots for each random effect
random <- ranef(m3)
valores = random[["Subject"]][["(Intercept)"]]
plot(valores,main="efeitos aleatórios versus índices")
abline(h=55, lty=3)
abline(h=-55, lty=3)
abline(h=0,lty=3,col=4)

plot_ranef(m3)

shapiro.test(valores)
## 
##  Shapiro-Wilk normality test
## 
## data:  valores
## W = 0.95017, p-value = 0.4277
  • Pendências no assunto Modelos Lineares com Efeitos Mistos:
    • Para o Exemplo 2 não foi verificado normalidade nos resíduos condicionais - verificar transformação nos dados ou métodos não paramétricos;
    • O pacote redres exibe vários tipos de resíduos - verificar as propriedades destes resíduos e se todos precisam exibir normalidade e homogeneidade das variâncias;
    • Explorar o efeito aleatório de diferentes formas - no intercepto + slope (inclinação da reta) - os modelos podem ser comparados através do AIC & BIC;
    • Os AICs & BICs só podem ser comparados com modelos distintos para banco de dados fixado;
    • Incluir mais covariáveis no modelo;
    • Incluir covariável Tempo de forma quantitativa (ao invés de fatores);
    • Comparar/estabelecer relações com Anova de medidas repetidas;