Usando o pacote NLME

Uma outra possível e elegante solução no R para este problema é utilizar a função lme do pacote nlme.

Note que a abordagem do problema por este pacote é um pouco diferente da forma apresentada no curso por se tratar de uma ferramente geral para modelos com efeitos aleatórios. Entretanto todos os elementos relevantes da análise estão incluídos nos resultados.

Vamos a seguir ver os comandos necessários e comentar os resultados.

Inicialmente temos que carregar o pacote nlme com o comando require.

A seguir criamos uma variável para indicar o efeito aleatório que neste exemplo chamamos de ex06$fa utilizando a função getGroups.

Feito isto podemos rodar a função lme que faz o ajuste do modelo.

> require(nlme)
[1] TRUE
> ex06$fa <- getGroups(ex06, ~ 1|forn/lot, level=2)
> ex06.av <- lme(resp ~ forn, random=~1|fa, data=ex06)
> ex06.av
Linear mixed-effects model fit by REML
  Data: ex06 
  Log-restricted-likelihood: -71.42198
  Fixed: resp ~ forn 
(Intercept)       forn2       forn3 
 -0.4166667   0.7500000   1.5833333 

Random effects:
 Formula: ~1 | fa
        (Intercept) Residual
StdDev:    1.307561 1.624483

Number of Observations: 36
Number of Groups: 12

Este modelo tem a variável forn como efeito fixo e a variável lot como efeito aleatório com o componente de variância $\sigma^2_{lote}$.

Além disto temos a variância residual$\sigma^2$.

A saída acima mostra as estimativas destes componentes da variância como sendo $\hat{\sigma}^2_{lote} = (1.307)^2=1.71$ e $\hat{\sigma}^2 = (1.624)^2=2.64$.

O comando anova vai mostrar a análise de variância com apenas os efeitos principais. O fato do programa não incluir o efeito aleatório de lotes na saída não causa problema algum. O comando intervals mostra os intervalos de confiança para os componentes de variância. Portanto para verificar a significância do efeito de lotes basta ver se o intervalo para este componente de variância exclui o valor $0$, o que é o caso neste exemplo conforme vamos abaixo.

> anova(ex06.av)
            numDF denDF   F-value p-value
(Intercept)     1    24 0.6043242  0.4445
forn            2     9 0.9690643  0.4158
> intervals(ex06.av)
Approximate 95% confidence intervals

 Fixed effects:
                 lower       est.    upper
(Intercept) -2.0772277 -0.4166667 1.243894
forn2       -1.8239746  0.7500000 3.323975
forn3       -0.9906412  1.5833333 4.157308

 Random Effects:
  Level: fa 
                    lower     est.    upper
sd((Intercept)) 0.6397003 1.307561 2.672682

 Within-group standard error:
   lower     est.    upper 
1.224202 1.624483 2.155644

Finalmente uma versão mais detalhada dos resultados pode ser obtida com o comando summary.

> summary(ex06.av)
Linear mixed-effects model fit by REML
 Data: ex06 
       AIC      BIC    logLik
  152.8440 160.3265 -71.42198

Random effects:
 Formula: ~1 | fa
        (Intercept) Residual
StdDev:    1.307561 1.624483

Fixed effects: resp ~ forn 
                 Value Std.Error DF    t-value p-value
(Intercept) -0.4166667 0.8045749 24 -0.5178718  0.6093
forn2        0.7500000 1.1378407  9  0.6591432  0.5263
forn3        1.5833333 1.1378407  9  1.3915246  0.1975
 Correlation: 
      (Intr) forn2 
forn2 -0.707       
forn3 -0.707  0.500

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.4751376 -0.7500844  0.0812409  0.7060895  1.8720268 

Number of Observations: 36
Number of Groups: 12

O próximo passo da seria fazer uma análise dos resíduos para verificar os pressupostos, semalhante ao que foi feito no experimento fatorial da aula anterior. Vamos deixar isto por conta do leitor.