27 Efeitos aleatórios

27.1 Componentes de variância

27.1.1 Introdução

O problema que ilustra este tópico consiste em encontrar valores "padrão"(ou de referência) para a teores de elementos químicos Para isto, amostras de referência supostamente de teores iguais foram enviadas a diferentes laboratórios nos quais determinações de teores foram feitas com replicações.

Como exemplo, considere os dados dos teores medidos de um único elemento mostrados a seguir e que podem ser obtidos em http://www.leg.ufpr.br/~paulojus/dados/MgO.xls

Lab   MgO  
  A  1.86  
  A  1.88  
  A  1.86  
  B  2.00  
  B  2.00  
  B  1.99  
  B  2.02  
  B  2.01  
  C  1.84  
  C  1.83  
  C  1.83  
  D  1.64  
  D  1.73  
  D  1.68  
  E  0.28  
  E  0.31  
  E  0.68  
  F  1.88  
  F  1.87  
  F  1.86  
  G  1.87  
  G  1.87  
  G  1.86  
  H  1.85  
  H  1.86  
  H  1.85

  > require(gdata)
  > mgo <- read.xls("MgO.xls")
  > head(mgo)
  > str(mgo)
  > summary(mgo)

Pode-se identificar duas fontes de variação nos valores medidos, uma devido à variabilidade entre laboratórios e outra devida à variabilidade das replicações feitas nos laboratórios. O objetivo é encontrar um valor "característico"para as amostras, que seria dado por alguma "média"adequada, associada a uma medida de variabilidade desta média, por exemplo dada por um intervalo de confiança. Além disto deseja-se estimar os "componentes de variância", isto é, medidas da variabilidade entre e dentro de laboratórios.

Nos resultados a seguir ajustamos um modelo ajustado com a função lme() do pacote nlme.

  > require(nlme)

O summary() mostra um resumo dos resultados mais importantes do ajuste do modelo incluindo as estimativas da média (Fixed Effects: Intercept) e dos desvios padrão entre laboratórios (Random Effects: Intercept) e das replicatas (Random Effects: Residual).

  > mgo.lme <- lme(MgO ~ 1, random = ~1 | Lab, mgo)
  > summary(mgo.lme)

  Linear mixed-effects model fit by REML
   Data: mgo
          AIC      BIC   logLik
    -13.69303 -10.0364 9.846515
  
  Random effects:
   Formula: ~1 | Lab
          (Intercept)   Residual
  StdDev:   0.5112672 0.07620438
  
  Fixed effects: MgO ~ 1
                 Value Std.Error DF  t-value p-value
  (Intercept) 1.675204 0.1813949 18 9.235126       0
  
  Standardized Within-Group Residuals:
          Min          Q1         Med          Q3         Max
  -2.00166549 -0.06901512 -0.02752389  0.10150856  3.24737682
  
  Number of Observations: 26
  Number of Groups: 8

O intervalo de confiança para média e as estimativas das variâncias e desvios padrões podem ser obtidos como mostrado a seguir.

  > intervals(mgo.lme, which = "fixed")

  Approximate 95% confidence intervals
  
   Fixed effects:
                 lower     est.    upper
  (Intercept) 1.294108 1.675205 2.056301
  attr(,"label")
  [1] "Fixed effects:"

  > VarCorr(mgo.lme)

  Lab = pdLogChol(1)
              Variance    StdDev
  (Intercept) 0.261394113 0.51126716
  Residual    0.005807107 0.07620438

27.1.2 Avaliando o ajuste e qualidade dos dados

Os resultados mostrados anteriormente devem ser vistos apenas como uma ilustração dos comandos básicos para obtenção dos resultados. Entretanto, não deve-se tomar os resultados obtidos como corretos ou definitivos pois uma análise criteriosa deve verificar anomalias dos dados e adequação a pressupostos do modelo.

Os gráficos de resíduos das figuras 69 e 70 mostram observações discrepantes e pode-se detectar que ocorrem nos dados do Laboratório E. No primeiro desses todos os resíduos são visualizados conjuntamente, enquanto que no segundo usa-se gráficos condicionais do sistema gráfico fornecido pelo pacote lattice para separar os resíduos de cada laboratório.


  > print(plot(mgo.lme))

PIC

Figura 69: Gráfico de resíduos do modelo ajustado



  > print(plot(mgo.lme, resid(.) ~ fitted(.) | Lab, abline = 0))

PIC

Figura 70: Gráfico de resíduos para cada laboratório do modelo ajustado.


A observação de valor 0.68 do laboratório E é bastante diferente das demais replicatas deste laboratório (0.28 e 0.31), sendo que este dado também foi considerado suspeito pela fonte dos dados. Uma possível alternativa é, em acordo com o responsável pelos dados, optar por remover este dado da análise o que pode ser feito com o comando a seguir.

  > mgo1 <- subset(mgo, !(Lab == "E" & MgO > 0.6))
  > dim(mgo1)
  [1] 25  2

O modelo ajustado assume que os dados possuem distribuição normal e os gráficos de perfil de verossimilhança do parâmetro da transformação Box-Cox na figura 71 mostram que, excluindo-se o dado atípico, a transformação não é necessária.

  > require(MASS)
  > with(mgo, boxcox(MgO ~ Lab, lam = seq(1.5, 5.5, len = 200)))
  > with(mgo1, boxcox(MgO ~ Lab, lam = seq(0, 3, len = 200)))


PIC
Figura 71: Perfis de verossimilhança do parâmetro da transformação Box-Cox na presença e ausência do ponto atípico do Laboratório E.


O modelo ajustado com o novo conjunto de dados apresenta resultados diferentes do anterior, reduzindo a estimativa de variância entre as replicatas.

  > mgo1.lme <- lme(MgO ~ 1, random = ~1 | Lab, mgo1)
  > summary(mgo1.lme)
  Linear mixed-effects model fit by REML
   Data: mgo1
          AIC       BIC   logLik
    -59.08601 -55.55185 32.54301
  
  Random effects:
   Formula: ~1 | Lab
          (Intercept)   Residual
  StdDev:   0.5577551 0.01831692
  
  Fixed effects: MgO ~ 1
                 Value Std.Error DF  t-value p-value
  (Intercept) 1.659078 0.1972321 17 8.411808       0
  
  Standardized Within-Group Residuals:
         Min         Q1        Med         Q3        Max
  -2.3652780 -0.3598894 -0.1781699  0.3673809  2.5482108
  
  Number of Observations: 25
  Number of Groups: 8
  > intervals(mgo1.lme, which = "fixed")
  Approximate 95% confidence intervals
  
   Fixed effects:
                 lower     est.    upper
  (Intercept) 1.242955 1.659078 2.075202
  attr(,"label")
  [1] "Fixed effects:"
  > VarCorr(mgo1.lme)
  Lab = pdLogChol(1)
              Variance     StdDev
  (Intercept) 0.3110907386 0.55775509
  Residual    0.0003355097 0.01831692


  > print(plot(mgo1.lme, resid(., type = "p") ~ fitted(.) | Lab,
  +     abline = 0))

PIC

Figura 72: Gráfico de resíduos para cada laboratório do modelo ajustado.


Além disto, nota-se que na verdade todas as observações do Laboratório E parecem atípicas com valores inferiores aos obtidos nos demais laboratórios. Poderia-se então considerar ainda remover todas as observações deste laboratório.

  > mgo2 <- subset(mgo, Lab != "E")
  > dim(mgo2)
  [1] 23  2
  > mgo2.lme <- lme(MgO ~ 1, random = ~1 | Lab, mgo2)
  > summary(mgo2.lme)
  Linear mixed-effects model fit by REML
   Data: mgo2
          AIC       BIC   logLik
    -78.17204 -74.89891 42.08602
  
  Random effects:
   Formula: ~1 | Lab
          (Intercept)   Residual
  StdDev:  0.09324064 0.01811513
  
  Fixed effects: MgO ~ 1
                 Value  Std.Error DF  t-value p-value
  (Intercept) 1.854012 0.03545001 16 52.29932       0
  
  Standardized Within-Group Residuals:
         Min         Q1        Med         Q3        Max
  -2.5091805 -0.3302089 -0.1587727  0.3606918  2.4590419
  
  Number of Observations: 23
  Number of Groups: 7
  > intervals(mgo2.lme, which = "fixed")
  Approximate 95% confidence intervals
  
   Fixed effects:
                 lower     est.    upper
  (Intercept) 1.778861 1.854012 1.929162
  attr(,"label")
  [1] "Fixed effects:"
  > VarCorr(mgo2.lme)
  Lab = pdLogChol(1)
              Variance    StdDev
  (Intercept) 0.008693816 0.09324064
  Residual    0.000328158 0.01811513

Os resultados são substancialmente diferentes e a decisão de exclusão on não dos dados deste Laboratório deve ser cuidadosamente investigada dentro do contexto destes dados e conjunto com especialista da área.

27.1.3 Fundamentos

Assumindo que efeitos aleatórios podem ser usados para descrever o efeito de laboratórios, podemos descrever os teores por um modelo de efeitos aleatórios:

Yij = μ + εi + ϵij,

em que yij são valores observados na j-ésima medida feita no i-ésimo laboratório, μ é o valor real do elemento na amostra padrão, εi ~ N(0ε2) é o efeito aleatório do i-ésimo laboratório e σ ε2 que representa a variabilidade de medidas fornecidas por diferentes laboratórios (entre laboratórios) e ϵij ~ N(0ϵ2) é o termo associado à j-ésima medida feita no i-ésimo laboratório e σ ϵ2 é a variabilidade das medidas de replicatas dentro dos laboratórios.

O problema então consiste em estimar μ e a variância associada à esta estimativa, que por sua vez está associada aos valores dos parâmetros de variância do modelo σε2 e σ ϵ2. Esses últimos parâmetros são chamados de componentes de variância. Diferentes métodos de estimação são propostos na literatura tais como estimadores de momentos baseados na análise de variância, estimadores minque (estimadores de norma quadrática mínima), estimadores de máxima verossimilhança e máxima verossimilhança restrita.

Sob o modelo assumido os observações tem distribuição normal

Y ~ N (1lμ,V ),
em que 1l é um vetor unitário de dimensão igual ao numero de observações n e V é a matriz de variâncias e covariâncias das observações com elementos dados por: V ar(Y i,j) = σε2 + σ ϵ2, a variância de cada observação individual; Cov(Y i,j,Y i,j) = σε2 a covariância entre observações diferentes do mesmo laboratório, e os demais elementos são nulos. No caso balanceado, isto é, igual número de replicatas nos diferentes laboratórios, a matriz V pode ser obtida por um produto de Kronecker simples entre matrizes diagonais e unitárias multiplicadas pelos componentes de variância.

Considerando os recursos computacionais atualmente disponíveis e as propriedades dos diferentes estimadores, nossa preferência é pelo uso de estimadores de máxima verossimilhança restrita. Estes estimadores são obtidos maximizando-se a função de verossimilhaça de uma projeção do vetor dos dados no espaço complementar os definido pela parte fixa do modelo. Tipicamente, os estimadores de σε2 e σ ϵ2 são obtidos por maximização numérica de tal função e o estimador do parâmetro de interesse e sua variância são então obtidos por:

ˆμ = (1lˆ
V-11l)-11l ˆ
V-1y (6)
V ˆar(ˆμ) = (1lˆV-11l)-1 (7)
em que ˆV é a matrix de variâncias e covariâncias estimada das observações obtida a partir das estimativas ˆσε2 e σˆ ϵ2.

No exemplo em questão são estes os estimadores utilizados para obter as estimativas mostradas na Sessão anterior (ver o resultado de summary(mgo1.lme) e com valores mostrados novamente a seguir.

  > names(mgo1.lme)

   [1] "modelStruct"  "dims"         "contrasts"    "coefficients" "varFix"
   [6] "sigma"        "apVar"        "logLik"       "numIter"      "groups"
  [11] "call"         "terms"        "method"       "fitted"       "residuals"
  [16] "fixDF"        "na.action"    "data"

  > mgo1.lme$coeff$fixed

  (Intercept)
     1.659078

  > VarCorr(mgo1.lme)[, 1]

     (Intercept)       Residual
  "0.3110907386" "0.0003355097"

O intervalo de confiança para média pode então ser obtido por:

              ∘ --------
ˆμ ± t1-α∕2,n- 1  Vˆar(ˆμ),

. Nos comandos a seguir mostramos a obtenção do intervalo segundo cálculos dessa expressão e a equivalência com o informado pela funçãom intervals.lme().

  > mgo1.lme$varFix

              (Intercept)
  (Intercept)   0.0389005

  > with(mgo1.lme, coefficients$fixed + qt(c(0.025, 0.975), df = fixDF$X) *
  +     sqrt(varFix))

  [1] 1.242955 2.075202

  > intervals(mgo1.lme, which = "fixed")

  Approximate 95% confidence intervals
  
   Fixed effects:
                 lower     est.    upper
  (Intercept) 1.242955 1.659078 2.075202
  attr(,"label")
  [1] "Fixed effects:"

Para uma observação individual o intervalo é dado por

             ∘  --------
y ± t1-α∕2,n-1   σ2ε + σ2ϵ;

e as estimativas ˆσε2 e ˆσ ϵ2 podem obtidas da seguinte forma.

  > vcomp <- as.numeric(VarCorr(mgo1.lme)[, 1])
  > vcomp

  [1] 0.3110907386 0.0003355097

O coeficiente de correlação intraclasse reflete a relação entre a variabilidade das observações dentro dos laboratórios em relação a variabilidade total. É definido ela expressão a seguir e calculado como mostrado nas linhas de comando.

       σ2ε
ρ = -2----2-.
    σε + σϵ

  > vcomp[1]/sum(vcomp)

  [1] 0.9989227

27.1.4 Alternativas de código

O pacote lme4 reimplementa algumas funcionalidades do nlme onde o modelo é definido indicando os termos aleatórios entre parênteses na fórmula e eliminando o uso do argumento random.

  > require(lme4)

O comando para se obter uma análise equivalente à anterior é mostrado a seguir. Os resultados são apresentados de forma diferente, prém os elementos são equivalentes.

  > mgo1.lmer <- lmer(MgO ~ 1 + (1 | Lab), mgo1)
  > summary(mgo1.lmer)

  Linear mixed model fit by REML
  Formula: MgO ~ 1 + (1 | Lab)
     Data: mgo1
      AIC    BIC logLik deviance REMLdev
   -59.09 -55.43  32.54   -66.52  -65.09
  Random effects:
   Groups   Name        Variance   Std.Dev.
   Lab      (Intercept) 0.31109074 0.557755
   Residual             0.00033551 0.018317
  Number of obs: 25, groups: Lab, 8
  
  Fixed effects:
              Estimate Std. Error t value
  (Intercept)   1.6591     0.1972   8.412

A opção padrão é o ajuste por máxima verossimilhança restrita. Estimativas de máxima verossimilhança podem ser obtidas usando o argumento REML=FALSE.

  > mgo1.lmer.ml <- lmer(MgO ~ 1 + (1 | Lab), mgo1, REML = FALSE)
  > summary(mgo1.lmer.ml)

  Linear mixed model fit by maximum likelihood
  Formula: MgO ~ 1 + (1 | Lab)
     Data: mgo1
      AIC    BIC logLik deviance REMLdev
   -60.56 -56.91  33.28   -66.56  -65.04
  Random effects:
   Groups   Name        Variance   Std.Dev.
   Lab      (Intercept) 0.27217959 0.521708
   Residual             0.00033552 0.018317
  Number of obs: 25, groups: Lab, 8
  
  Fixed effects:
              Estimate Std. Error t value
  (Intercept)   1.6591     0.1845   8.993