> print(plot(mgo.lme))
O problema que ilustra este tópico consiste em encontrar valores "padrão"(ou de referência) para os 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/\(\sim \)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(readxl)
Carregando pacotes exigidos: readxl
> mgo <- as.data.frame(read_excel("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)
Carregando pacotes exigidos: 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.63106 -9.974431 9.815529 Random effects: Formula: ~1 | Lab (Intercept) Residual StdDev: 0.5112631 0.07633555 Fixed effects: MgO ~ 1 Value Std.Error DF t-value p-value (Intercept) 1.675227 0.1813956 18 9.235211 0 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.99864232 -0.06883532 -0.02227551 0.10285379 3.24138039 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.294129 1.675227 2.056325
> VarCorr(mgo.lme)
Lab = pdLogChol(1) Variance StdDev (Intercept) 0.261389947 0.51126309 Residual 0.005827116 0.07633555
Os resultados mostrados anteriormente devem ser vistos apenas como uma ilustração dos comandos básicos para obtenção dos resultados. Entretanto, não se deve 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 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))
> print(plot(mgo.lme, resid(.) ~ fitted(.) | Lab, abline = 0))
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))) > #with(mgo2, boxcox(MgO~Lab, lam=seq(1.5,10, len=200)))
Carregando pacotes exigidos: MASS
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 -58.04521 -54.51105 32.02261 Random effects: Formula: ~1 | Lab (Intercept) Residual StdDev: 0.5577469 0.01888633 Fixed effects: MgO ~ 1 Value Std.Error DF t-value p-value (Intercept) 1.659103 0.1972315 17 8.411958 0 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.4333720 -0.3487900 -0.1439774 0.3565135 2.5808360 Number of Observations: 25 Number of Groups: 8
> intervals(mgo1.lme, which="fixed")
Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 1.242981 1.659103 2.075225
> VarCorr(mgo1.lme)
Lab = pdLogChol(1) Variance StdDev (Intercept) 0.3110815938 0.55774689 Residual 0.0003566936 0.01888633
> print(plot(mgo1.lme, resid(., type = "p") ~ fitted(.) | Lab, abline = 0))
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 -77.11573 -73.8426 41.55786 Random effects: Formula: ~1 | Lab (Intercept) Residual StdDev: 0.09317201 0.01872537 Fixed effects: MgO ~ 1 Value Std.Error DF t-value p-value (Intercept) 1.854044 0.03543846 16 52.31728 0 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.5760112 -0.3054222 -0.1399831 0.3484244 2.4812983 Number of Observations: 23 Number of Groups: 7
> intervals(mgo2.lme, which="fixed")
Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 1.778918 1.854044 1.92917
> VarCorr(mgo2.lme)
Lab = pdLogChol(1) Variance StdDev (Intercept) 0.0086810227 0.09317201 Residual 0.0003506395 0.01872537
Os resultados são substancialmente diferentes e a decisão de exclusão ou não dos dados deste Laboratório deve ser cuidadosamente investigada dentro do contexto destes dados e em conjunto com especialista da área.
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: \[Y_{ij} = \mu + \varepsilon _i + \epsilon _{ij} , \] em que \(y_{ij}\) são valores observados na \(j\)-ésima medida feita no \(i\)-ésimo laboratório, \(\mu \) é o valor real do elemento na amostra padrão, \(\varepsilon _i \sim N(0, \sigma ^2_{\varepsilon })\) é o efeito aleatório do \(i\)-ésimo laboratório e \(\sigma ^2_{\varepsilon }\) que representa a variabilidade de medidas fornecidas por diferentes laboratórios (entre laboratórios) e \(\epsilon _{ij} \sim N(0, \sigma ^2_{\epsilon })\) é o termo associado à \(j\)-ésima medida feita no \(i\)-ésimo laboratório e \(\sigma ^2_{\epsilon }\) é a variabilidade das medidas de replicatas dentro dos laboratórios.
O problema então consiste em estimar \(\mu \) e a variância associada à esta estimativa, que por sua vez está associada aos valores dos parâmetros de variância do modelo \(\sigma ^2_{\varepsilon }\) e \(\sigma ^2_{\epsilon }\). 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 as observações têm distribuição normal \[Y \sim N(\onevec \mu , V) , \] em que \(\onevec \) é um vetor unitário de dimensão igual ao número de observações \(n\) e \(V\) é a matriz de variâncias e covariâncias das observações com elementos dados por: \({\rm Var}(Y_{i,j}) = \sigma ^2_{\varepsilon } + \sigma ^2_{\epsilon }\), a variância de cada observação individual; \({\rm Cov}(Y_{i,j}, Y_{i,j^\prime }) = \sigma ^2_{\varepsilon }\) 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 verossimilhança de uma projeção do vetor dos dados no espaço complementar ao definido pela parte fixa do modelo. Tipicamente, os estimadores de \(\sigma ^2_{\varepsilon }\) e \(\sigma ^2_{\epsilon }\) 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: \begin {align} \hat {\mu } &= (\onevec ^\prime \hat {V}^{-1} \onevec )^{-1} \onevec ^\prime \hat {V}^{-1} y \label {eq:est-media}\\ \hat {{\rm Var}}(\hat {\mu }) &= (\onevec ^\prime \hat {V}^{-1} \onevec )^{-1} \label {eq:est-var} \end {align}
em que \(\hat {V}\) é a matriz de variâncias e covariâncias estimada das observações obtida a partir das estimativas \(\hat {\sigma }^2_{\varepsilon }\) e \(\hat {\sigma }^2_{\epsilon }\).
No exemplo em questão são estes os estimadores utilizados para obter as estimativas mostradas na seção anterior (ver o resultado de summary(mgo1.lme)), 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.659103
> VarCorr(mgo1.lme)[,1]
(Intercept) Residual "0.3110815938" "0.0003566936"
O intervalo de confiança para média pode então ser obtido por: \[ \hat {\mu } \pm t_{1-\alpha /2, n-1} \sqrt {\hat {{\rm Var}}(\hat {\mu })} , \] 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ção intervals.lme().
> mgo1.lme$varFix
(Intercept) (Intercept) 0.03890025
> with(mgo1.lme, coefficients$fixed + qt(c(0.025, 0.975), df=fixDF$X) * sqrt(varFix))
Warning in qt(c(0.025, 0.975), df = fixDF$X) * sqrt(varFix): Reciclar uma array de comprimento 1 na aritmética de vetor-array foi descontinuado. Em vez disso, use c() ou as.vector().
[1] 1.242981 2.075225
> intervals(mgo1.lme, which="fixed")
Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 1.242981 1.659103 2.075225
Para uma observação individual o intervalo é dado por \[ y \pm t_{1-\alpha /2, n-1} \sqrt {\sigma ^2_{\varepsilon } + \sigma ^2_{\epsilon }} ; \] e as estimativas \(\hat {\sigma }^2_{\varepsilon }\) e \(\hat {\sigma }^2_{\epsilon }\) podem ser obtidas da seguinte forma.
> vcomp <- as.numeric(VarCorr(mgo1.lme)[,1]) > vcomp
[1] 0.3110815938 0.0003566936
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 pela expressão a seguir e calculado como mostrado nas linhas de comando. \[ \rho = \frac {\sigma ^2_{\varepsilon }}{\sigma ^2_{\varepsilon }+\sigma ^2_{\epsilon }}. \]
> vcomp[1]/sum(vcomp)
[1] 0.9988547
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)
Carregando pacotes exigidos: lme4
Carregando pacotes exigidos: Matrix
Anexando pacote: 'lme4'
O seguinte objeto é mascarado por 'package:nlme': lmList
O comando para se obter uma análise equivalente à anterior é mostrado a seguir. Os resultados são apresentados de forma diferente, porém os elementos são equivalentes.
> mgo1.lmer <- lmer(MgO ~ 1 + (1|Lab), mgo1) > summary(mgo1.lmer)
Linear mixed model fit by REML ['lmerMod'] Formula: MgO ~ 1 + (1 | Lab) Data: mgo1 REML criterion at convergence: -64 Scaled residuals: Min 1Q Median 3Q Max -2.4334 -0.3488 -0.1440 0.3565 2.5808 Random effects: Groups Name Variance Std.Dev. Lab (Intercept) 0.3110816 0.55775 Residual 0.0003567 0.01889 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 ['lmerMod'] Formula: MgO ~ 1 + (1 | Lab) Data: mgo1 AIC BIC logLik -2*log(L) df.resid -59.5 -55.9 32.8 -65.5 22 Scaled residuals: Min 1Q Median 3Q Max -2.4333 -0.3482 -0.1434 0.3570 2.5809 Random effects: Groups Name Variance Std.Dev. Lab (Intercept) 0.2721697 0.52170 Residual 0.0003567 0.01889 Number of obs: 25, groups: Lab, 8 Fixed effects: Estimate Std. Error t value (Intercept) 1.6591 0.1845 8.993