5.1 Introdução


Nos experimentos descritos nos Capítulos 2 a 4, o objetivo era comparar a resposta entre diferentes níveis de variáveis ou fatores controláveis, a fim de prever qual poderia ser a resposta no futuro em níveis específicos dos fatores ou recomendar o melhor nível do fator para uso futuro.

Outro propósito da experimentação é estudar fontes de variabilidade na resposta. Por exemplo, sabe-se que o colesterol, a glicemia e outros testes de diagnóstico feitos por médicos variam de acordo com os procedimentos e o equipamento utilizado para fazer as medições. Experimentos podem ser conduzidos para descobrir quanto da variabilidade se deve ao equipamento e quanto se deve ao procedimento.

Simbolicamente \(\sigma^2_T = \sigma^2_p + \sigma^2_e\), onde \(\sigma^2_T\) é a variância total, \(\sigma^2_p\) e \(\sigma^2_e\) são as parcelas do total devidas ao procedimento e ao equipamento, respectivamente. \(\sigma^2_p\) e \(\sigma^2_e\) são chamados de componentes de variância.

Um experimento pode ser conduzido para coletar dados para que os componentes de variância possam ser estimados. Neste tipo de experimento, pode não haver qualquer interesse na diferença nas leituras médias de diagnóstico entre pedaços específicos de equipamentos porque há e continuará a haver, muitos em uso.

Existem pelo menos três razões para a realização de experimentos para estudar as fontes de variabilidade. Em alguns casos, o objectivo pode ser descritivo e os componentes de variância têm valor em si mesmos. Uma segunda razão para quantificar as fontes de variabilidade é obter informações sobre como reduzir a variância da resposta. Uma terceira razão para estudar as fontes de variabilidade é estimular ideias sobre as causas da variabilidade que poderiam ser testadas em experiências futuras.

Dois exemplos de onde os componentes de variância são úteis como medidas descritivas são na genética e em testes educacionais e psicológicos (ver Searle et al., 1992). Na criação de vacas leiteiras, a variabilidade na produção de leite pode ser particionada na quantidade devida à progenitora (\(\sigma^2_s\)) e à filha (\(\sigma^2_d\)), ou seja, \(\sigma^2_T = \sigma^2_s + \sigma^2_d\). A razão \(h = 4\sigma^2_s /(\sigma^2_s + \sigma^2_d)\) é chamada de herdabilidade e é altamente importante para os produtores de leite.

Em testes psicológicos e educacionais, a variabilidade nas pontuações dos testes pode ser dividida na variabilidade de pessoa para pessoa (\(\sigma^2_p\)) e nas pontuações repetidas dos testes para a mesma pessoa (\(\sigma^2_r\)), ou seja, \(\sigma^2_T = \sigma^2_p + \sigma^2_r\). Neste caso, \(\sigma^2_p /(\sigma^2_p + \sigma^2_r)\) é chamada de correlação intraclasse e valores elevados dela implicam confiabilidade do procedimento de teste.

No controle de qualidade industrial, há necessidade de reduzir a variabilidade nas medições durante o processo das principais características do produto e do processo. Se não for possível medir com precisão, não há esperança de controlar ou melhorar a qualidade.

A variabilidade total de medição dentro de uma planta pode ser atribuída ao equipamento de medição ou medidor e ao operador ou inspetor que faz a medição, ou seja, \(\sigma_T^2 = \sigma_g^2 + \sigma_o^2\). Para reduzir a variabilidade da medição, a gestão precisa saber onde concentrar os seus esforços.

Se a maior proporção da variabilidade for \(\sigma_g^2\), talvez seja necessário fazer esforços para recalibrar equipamentos de medição ou comprar mais equipamentos novos, mais precisos e consistentes. Por outro lado, se a principal fonte de variabilidade de medição for \(\sigma_o^2\), talvez um melhor treinamento dos operadores-inspetores possa resolver o problema.

Em alguns casos, um investigador gostaria de realizar uma experiência como as descritas nos Capítulos 2–4 para comparar a resposta média causada por diferentes níveis de factores controláveis; no entanto, ele ou ela tem um conhecimento tão limitado sobre o mecanismo em estudo que é difícil formular hipóteses sobre quais fatores ou níveis de fatores estudar. Nesta situação, determinar as fontes de variabilidade na resposta pode suscitar ideias sobre quais os factores que seriam mais rentáveis para estudar.

Por exemplo, saber se a maior parte da variabilidade em um processo industrial é de lote para lote ou dentro de um lote daria uma ideia sobre se os fatores que poderiam variar dentro de um lote ou os fatores que poderiam variar de um lote para outro, deve ser estudado em experimentos de otimização.

Nos casos em que um experimento fatorial ou um experimento fatorial foi conduzido e nada foi considerado significativo, Leitnaker and Cooper (2005) sugerem que um estudo de amostragem de acompanhamento, para classificar as fontes de variabilidade na resposta, pode explicar por que não foram encontrados fatores significativos.


5.2 Fatores aleatórios e experimentos de amostragem aleatória


Quando o objetivo da experimentação é estudar diferenças na resposta média causadas por diferenças nos níveis dos fatores, como os experimentos descritos nos Capítulos 2 a 4, os fatores do experimento são chamados de fatores fixos. Os níveis desses fatores são selecionados especificamente pelo experimentador. Por outro lado, quando o objetivo da experimentação é estudar a variância causada pela mudança nos níveis de um fator, o fator é chamado de fator aleatório.

Por exemplo, no modelo \(y_{ij} = \mu + \tau_i + \epsilon_{ij}\) do Capítulo 2, o fator \(\tau_i\) seria considerado um fator fixo. Embora não o tenhamos chamado de fator antes, o termo \(\epsilon_{ij}\) neste modelo seria considerado um fator aleatório. \(\epsilon_{ij}\) representa o efeito da \(j\)-ésima unidade experimental dentro do \(i\)-ésimo nível do fator de tratamento. No experimento de crescimento de pão descrito no Capítulo 2, a unidade experimental foi o pão. Pães replicados, isto é, \(j = 1,\cdots, 4\) foram usados neste experimento para que \(\sigma^2\), a variância das unidades experimentais, pudesse ser estimada e usada para julgar a significância do fator fixo, ou seja, o tempo de subida.

Como o objetivo de incluir múltiplos níveis ou pães em cada tempo de subida era estimar \(\sigma^2\), \(\epsilon_{ij}\) é considerado um fator aleatório.

Enquanto os níveis dos fatores fixos são escolhidos especificamente pelo experimentador, os níveis dos fatores aleatórios são apenas amostras de níveis possíveis que poderiam ter sido utilizados. Por exemplo, no experimento de crescimento do pão, o experimentador escolheu 35, 40 e 45 minutos para estudar os níveis de tempo de crescimento.

No entanto, os quatro pães replicados utilizados para cada tempo de fermentação representam apenas uma amostra dos pães que poderiam ter sido utilizados na experiência. Por esse motivo, os experimentos usados para estudar variâncias podem ser considerados experimentos de amostragem aleatória ou RSE, uma vez que os níveis dos fatores são apenas uma amostra de níveis possíveis.

Por exemplo, considere os dados da Tabela 5.1 modelados após a pesquisa internacional de apolipoproteínas (Apo) realizada em 1984-1985 apresentada em Henderson et al., 1987. Apo A–I é conhecida por ajudar a eliminar o colesterol das artérias. No entanto, o genótipo apo produz valores preditivos fracos no rastreio de aterosclerose clinicamente definida. Isto pode ser em parte devido à dificuldade de medição.

\[ \begin{array}{ccccc}\hline \mbox{Laboratório:} & \mbox{A} & \mbox{B} & \mbox{C} & \mbox{D} \\\hline & 1.195 & 1.155 & 1.021 & 1.163 \\ & 1.144 & 1.173 & 1.037 & 1.171 \\ & 1.167 & 1.171 & 1.022 & 1.182 \\ & 1.249 & 1.175 & 1.064 & 1.184 \\ & 1.177 & 1.153 & 1.094 & 1.175 \\ & 1.217 & 1.139 & 0.992 & 1.134 \\ & 1.187 & 1.185 & 1.072 & 1.169 \\ & & 1.144 & & 1.136 \\\hline \end{array} \] Tabela 5.1: Concentrações medidas de Apo A–I por código de laboratório.

O objetivo do estudo foi examinar e quantificar a variação entre laboratórios no que diz respeito à medição de um material de referência internacional para apo A-I e B. Várias medições das concentrações relativas de preparações do material de referência foram feitas por 28 laboratórios selecionados. A Tabela 5.1 mostra dados de quatro laboratórios representativos.

O modelo para os dados pode ser escrito como: \[\begin{equation*} \tag{5.1} y_{ij}=\mu+t_i+\epsilon_{ij}, \end{equation*}\] onde \(y_{ij}\) é a \(j\)-ésima medição da concentração de apo A-I no \(i\)-ésimo laboratório, \(\mu\) é a concentração média geral medida, \(t_i\) é o efeito do laboratório e \(\epsilon_{ij}\) é o efeito da \(j\)-ésima medição no \(i\)-ésimo laboratório. Não houve interesse em comparar as concentrações medidas de apo A-I entre os laboratórios específicos do estudo, uma vez que podem ser consideradas uma amostra aleatória ou amostra representativa, de vários laboratórios ao redor do mundo.

Como o objetivo de incluir vários laboratórios era estimar o componente de variância \(\sigma_t^2\), na concentração medida devido ao laboratório, \(t_i\) pode ser considerado um fator aleatório e o experimento pode ser pensado como um experimento de amostragem. Observe que a letra romana \(t_i\) foi utilizada para representar o efeito de região aleatória, enquanto no modelo (2.2) do Capítulo 2 a letra grega \(\tau_i\) foi utilizada para representar o fator de tratamento fixo.

As medições replicadas de apo A-I feitas em cada laboratório são apenas amostras das medições possíveis que poderiam ser feitas naquele laboratório. Como foram feitas múltiplas medições em cada laboratório para estimar \(\sigma^2\), \(\epsilon_{ij}\) também pode ser considerado um fator aleatório. Exceto pelo \(\epsilon\) usado para representar a unidade experimental aleatória, a convenção neste texto será usar letras romanas para representar fatores aleatórios e letras gregas para representar fatores fixos.

As suposições usuais em relação ao modelo (5.1) são que os efeitos aleatórios, \(t_i\) e \(\epsilon_{ij}\), são independentes e normalmente distribuídos com médias zero e variâncias iguais às componentes de variância \(\sigma_t^2\) e \(\sigma^2\), respectivamente. Como a variância de uma soma de variáveis aleatórias independentes é a soma das variâncias, \(\sigma_y^2 =\sigma_t^2 + \sigma^2\). Dados de experimentos de amostragem podem ser usados para particionar a variância na resposta, \(\sigma_y^2\), nos dois componentes de variância que a compõem.


5.3 Projetos de amostragem de um fator


Experimentos de amostragem de um fator podem ser usados para particionar a variabilidade em duas fontes. Como exemplo de partição da variância em duas fontes, considere os experimentos com helicópteros de papel descritos no Exercício 2.2. Nesse exercício poderão ter havido algumas diferenças na descrição da unidade experimental entre os alunos.

Alguns podem argumentar que a unidade experimental era uma folha de papel da qual foi cortado o desenho de um helicóptero. Outros podem argumentar que foi o teste ou as condições do ar, no momento em que o helicóptero foi lançado e cronometrado. Se a primeira definição fosse usada, então os experimentos replicados consistiriam em fabricar, lançar e cronometrar vários helicópteros do mesmo projeto uma vez, sendo cada um feito a partir de uma folha de papel diferente. Se a segunda definição fosse usada, os experimentos replicados consistiriam em lançar e cronometrar repetidamente um helicóptero. Uma maneira prática de decidir como definir a unidade experimental seria particionar a variabilidade nos tempos de queda em variabilidade de helicóptero para helicóptero e variabilidade entre lançamentos repetidos do mesmo helicóptero.

Se uma parte substancial da variabilidade ocorresse entre helicópteros do mesmo projeto cortados de diferentes pedaços de papel, haveria razão para fabricar vários helicópteros para réplicas. Se, por outro lado, toda a variabilidade estivesse entre os tempos de lançamento do mesmo helicóptero, não haveria razão para fabricar vários helicópteros do mesmo projeto para réplicas. Nesse caso, lançamentos repetidos do mesmo helicóptero poderiam ser considerados réplicas.

A variabilidade nos tempos de queda do helicóptero pode ser dividida na variabilidade de helicóptero para helicóptero e na variabilidade do helicóptero usando um experimento de amostragem de um fator. Para fazer isso, primeiro selecione aleatoriamente seis folhas de papel. De cada uma dessas folhas, corte e dobre um helicóptero padrão com largura de corpo = 4,25”, comprimento de cauda = 4,0” e comprimento de asa = 6,5”. Solte e cronometre cada um dos seis helicópteros três vezes cada, de acordo com uma ordem aleatória como a criada na Seção 2.2.

Como segundo exemplo de partição da variabilidade em duas fontes, considere o seguinte exemplo apresentado por Davies (1949). Um fabricante de corante queria saber se havia uma contribuição apreciável para a variabilidade nos rendimentos de cor do corante devido à qualidade do lote de ácido intermediário utilizado.

Lotes intermediários de ácido foram produzidos em uma etapa e, em uma etapa posterior, o ácido foi utilizado para produzir o corante. O objetivo era manter o rendimento da cor do corante consistentemente alto. Se a maior parte da variação foi causada por diferenças entre os lotes ácidos intermediários, então os esforços de melhoria deveriam ser concentrados no processo que produz os lotes ácidos.

Se a maior parte da variação ocorreu nas preparações do corante feitas a partir do mesmo lote ácido, os esforços de melhoria deveriam ser concentrados na etapa do processo de produção do corante. Um experimento de amostragem foi realizado em que seis amostras representativas do ácido intermediário H foram retiradas da etapa do processo de fabricação que o produz. De cada amostra de ácido foram feitas em laboratório cinco preparações do corante Naftaleno 12B, representativas das preparações que poderiam ser feitas com cada amostra.

Os dados do experimento de amostragem são mostrados na Tabela 5.2. Os rendimentos são dados em gramas de cor padrão.

\[ \begin{array}{ccccccc}\hline \mbox{Amostra do ácido H} & 1 & 2 & 3 & 4 & 5 & 6 \\\hline & 1440 & 1490 & 1510 & 1440 & 1515 & 1445 \\ \mbox{Rendimentos individuais} & 1440 & 1495 & 1550 & 1445 & 1595 & 1450 \\ \mbox{em gramas de cor padrão} & 1520 & 1540 & 1560 & 1465 & 1625 & 1455 \\ & 1545 & 1555 & 1595 & 1545 & 1630 & 1480 \\ & 1580 & 1560 & 1605 & 1595 & 1635 & 1520 \\\hline \end{array} \] Tabela 5.2: Rendimentos de negro de Naftaleno 12B.


5.4 Estimando componentes de variância


O modelo (5.1) pode ser expresso em termos matriciais como: \[\begin{equation*} \tag{5.2} \pmb{y}=\pmb{X}\beta+\epsilon \end{equation*}\] onde \(\beta^\top =(\mu,t^\top)\) é \(t^\top\) é o vetor de efeitos aleatórios. As suposições de independência e normalidade podem ser expressas como: \[\begin{equation*} \tag{5.3} \begin{pmatrix} t \\ \epsilon \end{pmatrix} \sim \mbox{MVN}\begin{pmatrix} \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_t^2 \mbox{I}_t & 0 \\ 0 & \sigma^2 \mbox{I}_n \end{pmatrix}\end{pmatrix}, \end{equation*}\] sendo que MVN representa a distribuição normal multivariada, \(\mbox{I}_t\) é uma matriz identidade \(t\times t\), \(t\) é o número de níveis do fator aleatório \(t_i\) e \(n = tr\) é o número total de experimentos ou execuções.

Com as suposições de independência e normalidade, existem vários métodos para estimar componentes de variância a partir dos dados obtidos em experimentos de amostragem.


5.4.1 Estimadores do método dos momentos


R.A. Fisher foi o primeiro a mostrar como estimar os componentes de variância a partir de uma análise de variância. Para isso, utiliza-se o método dos momentos. O modelo para o experimento de amostragem unifatorial é dado pela Equação (5.1), e a tabela de análise de variância é idêntica à Tabela 2.2 no Capítulo 2. O teste \(F\) na tabela de análise de variância pode ser usado para testar a hipótese nula \(H_0 ∶ \sigma_t^2 = 0\) contra a alternativa \(H_1 ∶ \sigma_t^2 > 0\).

Para estimar os componentes de variância da ANOVA, os quadrados médios são igualados aos seus valores esperados e as equações simultâneas são resolvidas. Quando o fator de tratamento é um fator aleatório, como no modelo (5.1), e há um número igual de repetições em cada nível, o termo \(msT\) na tabela ANOVA segue uma distribuição que é um múltiplo da distribuição qui-quadrado central.

Isso é diferente do modelo com fator de tratamento de efeito fixo. Nesse caso, descrito na Seção 2.3.5, a distribuição do \(msT\) era um qui-quadrado não central. Em ambos os casos, o \(msT\) da ANOVA pode ser representado como a forma quadrática \(\pmb{y}^\top \pmb{A}\pmb{y}\) e no modelo de efeitos aleatórios (5.1), Hartley (1967) mostrou que seu valor esperado pode ser escrito como \(\sigma^2 + c\sigma_t^2\), onde \(c = \sum_i \pmb{x}_i^\top \pmb{A}\pmb{x}_i\) e \(\pmb{x}_i\) é um indicador do \(i\)-ésimo nível de \(t_i\), essa é a (\(i\) + 1)-ésima coluna na matriz \(\pmb{X}\) conforme mostrado na Equação (2.6).

Quando há um número igual de repetições \(r\), em cada nível do fator aleatório, o coeficiente no quadrado médio esperado simplifica para \(c = r\). O método dos estimadores de momentos para os componentes de variância pode ser utilizado quando houver um número igual de repetições em cada nível do fator aleatório conforme mostrado na Tabela 5.2, ou um número desigual conforme mostrado na Tabela 5.1.

Para o caso de réplicas iguais, as estimativas revelam-se estimadores não viciados uniformemente melhores, mas para o caso desigual, não existem estimadores que sejam uniformemente melhores, ver Searle et al., 1992.

Para ilustrar a estimativa dos componentes de variância usando R, o código para abrir o conjunto de dados contendo dados do experimento de amostragem de medição de Apo, mostrado na Tabela 5.1, juntamente com a chamada da função aov para produzir a ANOVA são mostrados a continuação.

library(daewr)
mod1 = aov( conc ~ lab, data = Apo )
sm1 = summary(mod1)
sm1
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## lab          3 0.09223 0.03074   42.11 4.01e-10 ***
## Residuals   26 0.01898 0.00073                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Neste caso, a tabela ANOVA não é apenas impressa, mas também armazenada no objeto sm1. Ao fazer isso, os quadrados médios do tratamento e do erro podem ser recuperados para cálculos posteriores.

O primeiro passo no cálculo do coeficiente quadrático médio esperado \(c\) usando o método de Hartley (1967) é extrair os indicadores de cada nível do fator de tratamento \(\pmb{x}_i\), da matriz do modelo. O código R para fazer isso é mostrado abaixo.

X = as.matrix( model.matrix(mod1) )
labB = X[,2]
labC = X[,3]
labD = X[,4]
labA = X[,1]-(labB+labC+labD)


A próxima etapa é calcular uma ANOVA usando cada vetor indicador como resposta, extrair o quadrado médio do tratamento de cada ANOVA e, finalmente, somar os resultados para obter o coeficiente \(c\) conforme mostrado abaixo.

s1 = summary(aov (labA ~ Apo$lab ))
x1 = as.matrix( s1[[1]][1,3] )
s2 = summary(aov( labB ~ Apo$lab ))
x2 = as.matrix( s2[[1]][1,3] )
s3 = summary(aov(labC ~ Apo$lab ))
x3 = as.matrix( s3[[1]][1,3] )
s4 = summary(aov(labD ~ Apo$lab ))
x4 = as.matrix( s4[[1]][1,3] )
c = x1+x2+x3+x4


Conforme mostrado na Seção 2.4.3, o valor esperado do \(msE\) na tabela ANOVA é \(\sigma^2\), portanto, duas equações simultâneas, nos dois componentes de variância desconhecidos \(\sigma^2\) e \(\sigma_t^2\), podem ser criadas igualando os quadrados médios do modelo e do erro a seus valores esperados.

O código R para recuperar os quadrados médios e imprimir os resultados é mostrado abaixo.

sigma2 = as.matrix( sm1[[1]][2,3] )
mslab = as.matrix( sm1[[1]][1,3] )
cat(" Mean Square for Lab = ",mslab,"\n"," Mean Square for Error = ", sigma2,"\n",
    " Expected Mean Square for Lab","\n", "Var(error)+",c,"Var(Lab)","\n")
##  Mean Square for Lab =  0.03074443 
##   Mean Square for Error =  0.0007301573 
##   Expected Mean Square for Lab 
##  Var(error)+ 7.488889 Var(Lab)


A solução para as equações que igualam os quadrados médios aos seus valores esperados é mostrada abaixo. \[ 0.03074443 = \sigma^2 + 7.488889 \sigma_t^2\cdot \] Sabemos que \(\widehat{\sigma}^2\)=0.0007301573, então \[ \widehat{\sigma}_t^2 = (0.03074443-0.0007301573)/ 7.488889 = 0.004007841\cdot \]

Essas equações podem ser resolvidas usando R conforme mostrado no código abaixo.

sigma2t = (mslab - sigma2) / c
cat("Method of Moments Variance Component Estimates","\n", "Var(error)=",sigma2, 
    "\n","Var(Lab)=",sigma2t,"\n")
## Method of Moments Variance Component Estimates 
##  Var(error)= 0.0007301573 
##  Var(Lab)= 0.00400784


Como \(\widehat{\sigma}^2 < \widehat{\sigma}_t^2\), essas estimativas mostram que há muita ou mais variabilidade entre laboratórios do que dentro de um laboratório.


5.4.2 Estimativas de intervalo de confiança


Quando as suposições de normalidade se aplicam e há um número igual de réplicas para cada nível do fator aleatório no modelo (5.1), existem estimativas de intervalo exato para \(\sigma^2\), \(\sigma_t^2 /(\sigma^2 + \sigma_t^2)\), \(\sigma^2 /(\sigma^2 +\sigma_t^2)\) e \(\sigma_t^2 /\sigma^2\) com base nas distribuições dos quadrados médios.

A função aov em R não produz esses intervalos de confiança, mas eles podem ser calculados manualmente a partir das estatísticas produzidas na tabela ANOVA. Se houverem \(i = 1,\cdots,T\) níveis do fator aleatório e \(j = 1,\cdots,r\) replica em cada nível do fator aleatório, então a Tabela 5.3, retirada de Searle et al. (1992), contém fórmulas exatas para intervalos de confiança nas linhas 1, 3-5 e uma fórmula aproximada na linha 2.

Tabela 5.3: Intervalos de confiança para funções de componentes de variância em modelo aleatório de um fator com replicação igual.

Como exemplo de cálculo dos intervalos de confiança, considere os dados do experimento de amostragem mostrado na Tabela 5.2 para estudar a variabilidade nos rendimentos de cor dos corantes. Da tabela ANOVA \(ssT = 56.358\), \(ssE = 58.830\), \(T = 6\), \(r = 5\) e \(F = msT/msE = 4.59847\). O percentil superior 0.975 da distribuição qui-quadrado com 24 graus de liberdade é 39.36.

Isso pode ser obtido usando a função R qchisq como qchisq(0.975,24). O percentil superior 0.975 da distribuição \(F\) com 5 e 24 graus de liberdade é 3.15482. Isso pode ser obtido usando a função R qf como qf(0.975,5,24). Os demais percentis da distribuição \(F\)-Fischer e do qui-quadrado podem ser obtidos de forma semelhante.

O intervalo de confiança exato de 95% para \(\sigma^2\) é dado por: \[ \left(\dfrac{ssE}{\chi_{24,0.975}^2},\dfrac{ssE}{\chi_{24,0.025}^2}\right) = \left(\dfrac{58830}{39.36},\dfrac{58830}{12.4}\right) = (1494.51,4744.35)\cdot \]

O intervalo de confiança aproximado de 90% para \(\sigma_t^2\) é dado por: \[ \left(\dfrac{ssT \big(1-F_{5,24,0.975}/F\big)}{r \, \chi_{5,0.975}^2},\dfrac{ssT\big(1-F_{5,24,0.025}/F \big) }{r \, \chi_{5,0.025}^2}\right) = \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\ = \left(\dfrac{56358(1-3.5482/4.60)}{5(12.8325)},\dfrac{56385(1-0.15929/4.60)}{5(0.83121)}\right) = (275.9551,13097.168) \] e o intervalo de confiança exato de 95% em \(\sigma_t^2/(\sigma_t^2 +\sigma^2)\) é dado por: \[ \left(\dfrac{F/F_{5,24,0.975}-1}{r + F/F_{5,24,0.975} -1},\dfrac{F/F_{5,24,0.025}-1 \big) }{r + F/F_{5,24,0.025}-1}\right) = \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\ = \left(\dfrac{4.59847/3.15482-1}{5+4.59847/3.15482-1},\dfrac{4.59847/0.15929-1}{5+4.59847/0.15929-1}\right) = (0.0838,0.8479)\cdot \]

Observe que a estimativa intervalar de \(\sigma_t^2\) é muito mais ampla que a estimativa intervalar de \(\sigma^2\). Isso ocorre porque \(ssT\) possui apenas 5 graus de liberdade na ANOVA, enquanto \(ssE\) possui 24 graus de liberdade.


5.4.3 Estimadores de máxima verossimilhança e REML


Embora os estimadores do método dos momentos sejam melhores estimadores uniformemente não viciados, eles têm uma propriedade infeliz. Quando \(msT\) for menor que \(msE\) na análise de variância, o estimador de \(\sigma_t^2\) será negativo. Isto pode acontecer com bastante frequência se \(\sigma_t^2/\sigma^2\leq 0.10\), e houver menos de \(T = 10\) níveis do fator aleatório \(t_i\).

A máxima verossimilhança (ML) e a máxima verossimilhança reduzida ou restrita (REML) são métodos preferidos de estimação que evitam esse problema. REML é uma adaptação da técnica de máxima verossimilhança que maximiza parte da verossimilhança. O fato dos estimadores de máxima verossimilhança não poderem ficar fora do seu espaço de parâmetros impede que ambos os métodos ML e REML obtenham estimativas negativas de \(\sigma_t^2\).

Para entender como funcionam a máxima verossimilhança e o REML, consideraremos o caso de replicação igual. Dado o modelo e as suposições nas Equações (5.2) e (5.3), a distribuição de \(\pmb{y}\) pode ser escrita como: \[\begin{equation*} \tag{5.4} \pmb{y}\sim MVN(\mu\pmb{1},\pmb{V}), \end{equation*}\] onde \(\pmb{V}\) é uma matriz bloco diagonal com \(T\) blocos de \((\sigma_t^2\pmb{J}_r+\sigma^2\pmb{1}_r)\) ao longo da diagonal.

A função de verossimilhança é \[\begin{equation*} \tag{5.5} L(\mu,\pmb{V} \, | \, \pmb{y}) =\dfrac{1}{(2\pi)^{n/2}|\pmb{V}|^{1/2}}\exp\left(-\dfrac{1}{2}\big(\pmb{y}-\mu\pmb{1}_n\big)^\top \pmb{V}^{-1}\big(\pmb{y}-\mu\pmb{1}_n\big) \right)\cdot \end{equation*}\] Para o caso de replicação igual, isso pode ser simplificado para: \[\begin{equation*} \tag{5.6} L(\mu,\sigma^2,\lambda \, | \, \pmb{y}) =\dfrac{1}{(2\pi)^{n/2}\sigma^{2[\frac{1}{2}n]}\lambda^{T/2}}\exp\left(-\dfrac{1}{2}\Big(\frac{ssE}{\sigma^2}+\frac{ssT}{\lambda}+\frac{(\overline{y}_{\cdot\cdot}-\mu)^2}{\lambda/n} \Big)\right), \end{equation*}\] onde \(\lambda=\sigma^2+r\sigma_t^2\). As estimativas de máxima verossimilhança são obtidas maximizando a verossimilhança em relação a \(\mu\), \(\sigma^2\) e \(\lambda\).

Os estimadores REML de \(\sigma^2\) e \(\sigma_t^2\) são obtidos pela maximização da verossimilhança \[\begin{equation*} \tag{5.7} L(\sigma^2,\sigma_t^2 \, | \, ssT,ssE) = \dfrac{L(\mu,\sigma^2,\lambda \, | \, \pmb{y})}{L(\mu \, | \, \overline{y}_{\cdot\cdot})} \end{equation*}\] a qual pode ser obtida fatorando \(L(\mu, \sigma^2, \lambda \, | \, \pmb{y})\) usando o fato de que \(ssE\) e \(ssT\) são independentes de \(\overline{y}_{\cdot\cdot}\). A maximização pode ser feita analiticamente para o caso de replicação igual conforme mostrado por Searle et al. (1992) e pode ser feita numericamente para o caso desequilibrado.

Uma propriedade desejável das estimativas REML é que elas sejam iguais ao método dos momentos, estimativas de análise de variância, quando há replicação igual em cada nível do fator aleatório e \(msT > msE\). Os estimadores de máxima verossimilhança e os estimadores REML são calculados pelo pacote R lme4 usando uma solução numérica para os casos desequilibrados e balanceados.

Para ilustrar os estimadores REML, considere o seguinte caso. Um fabricante de misturas para sopas secas embaladas estava enfrentando uma variabilidade excessiva nos pesos das embalagens de um componente da mistura para sopas secas chamado “intermix”. A mistura é uma mistura de ingredientes saborosos, como óleo vegetal, sal e assim por diante.

Muita mistura em um pacote de sopa confere um sabor muito forte e pouca mistura dá um sabor muito fraco. Foi um processo de duas etapas para fazer a mistura de sopa embalada. O primeiro passo foi fazer uma grande quantidade de sopa e secar na secadora rotativa. Em seguida, o lote de sopa seca foi colocado em um misturador, onde a mistura foi adicionada através das portas à medida que era misturada.

Em seguida, foi embalado em sacos lacrados de peso uniforme. Havia vários fatores que poderiam ser alterados na primeira etapa, produção do lote de sopa seca, e vários fatores que poderiam ser alterados na segunda etapa, adição da mistura e mistura, que poderiam afetar a variabilidade do peso do mistura em cada saco selado.

Um experimento fatorial deveria ser planejado para descobrir quais fatores afetavam a variabilidade nos pesos da mistura. Para determinar quais fatores incluir em um experimento fatorial, um primeiro passo razoável seria particionar a variabilidade no peso da mistura entre a variabilidade de lote de sopa para lote de sopa e a variabilidade dentro de um lote causada pelo processo de mistura e adição de mistura para a sopa seca. Se houvesse pouca variabilidade de lote para lote, o experimento só precisaria considerar fatores envolvidos na etapa de mistura.

Para particionar a variabilidade nos pesos das embalagens entre lotes e dentro dos lotes, foi realizada uma experiência de amostragem. Uma amostra aleatória de quatro lotes foi selecionada ao longo de um mês de produção. De cada lote, três amostras de 1 lb da mistura de sopa acabada foram retiradas do tanque de mistura enquanto ela era embalada. Uma libra avoirdupois ou libra internacional, abreviatura: “lb”, ou por vezes # nos Estados Unidos, é a unidade de massa equivalente a exatamente 0.45359237 quilogramas ou 453.59237 gramas. Também a libra é chamada de pound avoirdupois ou simplesmente pound.

O peso da mistura foi determinado para cada amostra. Os resultados do experimento de amostragem são mostrados na Tabela 5.4. A função lmer no pacote R lme4 pode calcular as estimativas REML de componentes de variância. O código para abrir o conjunto de dados para os dados da mistura de sopa seca e chamar a função lmer para produzir as estimativas REML de \(\sigma_t^2\) e \(\sigma^2\) é mostrado logo a seguir e os resultados são mostrados abaixo dos comandos.

\[ \begin{array}{cccc}\hline \mbox{Lote} \, \setminus \mbox{Pesos} & & & \\\hline 1 & 0.52 & 2.94 & 2.03 \\ 2 & 4.59 & 1.26 & 2.78 \\ 3 & 2.87 & 1.77 & 2.68 \\ 4 & 1.38 & 1.57 & 4.10 \\\hline \end{array} \] Tabela 5.4: Variabilidade nos pesos da mistura de sopa seca.

A fórmula do modelo weight ~ 1 +(1|batch) para a função lmer é diferente da declaração do modelo para as funções lm ou aov mostradas no Capítulo 2. A resposta está à esquerda do operador “~” como antes, mas o 1 à direita indica que uma média global será adequada. A barra vertical “∣” no segundo termo da fórmula do modelo separa os efeitos fixos dos efeitos aleatórios. O termo à esquerda da barra vertical é um efeito fixo, ou seja, a média geral e o termo à direita, o lote neste exemplo, é um efeito aleatório.

library(daewr)
library(lme4)
rmod2 = lmer( weight ~ 1 + (1|batch), data = soupmx)
summary(rmod2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: weight ~ 1 + (1 | batch)
##    Data: soupmx
## 
## REML criterion at convergence: 37.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.56147 -0.71722 -0.01614  0.43230  1.86604 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  batch    (Intercept) 0.00     0.000   
##  Residual             1.41     1.187   
## Number of obs: 12, groups:  batch, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   2.3742     0.3428   6.926
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')


Na saída, a tabela ANOVA familiar está faltando porque não é necessária para obter as estimativas REML. Em vez disso, a saída mostra os estimadores numéricos resultantes. Pode-se observar que \(\widehat{\sigma}_t^2\), a variabilidade entre lotes, é estimada em zero, enquanto \(\widehat{\sigma}^2\), a variabilidade causada pelo processo de mistura dentro de um lote, é estimada em 1.41. Se o método dos estimadores de momentos for usado com esses dados, deixados como exercício, o estimador de \(\sigma_t^2\) acaba sendo negativo.

As conclusões desta experiência de amostragem indicam que novas experiências para identificar as causas da variabilidade na mistura devem concentrar-se em fatores que podem ser variados na segunda etapa do processo, onde a mistura é adicionada e misturada com o lote. Esses experimentos adicionais foram realizados e serão mostrados no próximo capítulo.

Em processos que contêm mais de duas etapas, Leitnaker and Cooper (2005) mostram que experimentos de amostragem em vários estágios, a serem descritos na Seção 5.6, são muito úteis para identificar etapas do processo e fatores dentro dessas etapas que seriam bons candidatos para estudos adicionais com fatorial tipo experimentos.

Um intervalo de confiança aproximado assintótico para os componentes de variância pode ser obtido no pacote R lme4 usando o método da verossimilhança perfilada. As funções profile e confint fornecem um intervalo de confiança para a raiz quadrada, ou seja, para o desvio padrão dos componentes de variância.

O código abaixo mostra como obter os intervalos de confiança usando os dados da Tabela 5.2.

library(daewr)
library(lme4)
pr1 = profile( fm1M <- lmer( yield ~ 1 + (1| sample), data = Naph, REML = FALSE))
confint (pr1) # 95% confidence interval on sigma
##                  2.5 %     97.5 %
## .sig01        12.19854   84.06305
## .sigma        38.22998   67.65770
## (Intercept) 1486.45150 1568.54849


A quadratura do ponto final esquerdo e direito do intervalo para .sigma fornece o intervalo de confiança de 95% (1461.53–4577.56) para \(\sigma^2\), o componente de variância dos rendimentos dentro da amostra de negro de Naftaleno 12B. Isto está muito próximo do intervalo de confiança exato encontrado na Seção 5.4.2.

Uma estimativa de intervalo semelhante pode ser obtida para \(\sigma_t^2\). Os intervalos de confiança assintóticos estarão razoavelmente próximos dos intervalos de confiança exatos produzidos utilizando as fórmulas da Tabela 5.3, quando os graus de liberdade para o termo correspondente ao componente de variância que está sendo estimado for superior a 45.


5.4.4 Determinação do tamanho da amostra para estudos de amostragem de um fator


Como há dois componentes de variância sendo estimados em um plano amostral de um fator, há duas coisas a serem consideradas ao determinar o tamanho da amostra. Para estimar com precisão a variância replicada \(\sigma^2\), o importante a considerar é o número de graus de liberdade para o erro \(\nu^2 = t(r-1)\).

A precisão da estimativa do fator aleatório \(\sigma_t^2\), será sempre relativa à precisão na estimativa de \(\sigma^2\). A precisão para estimar \(\sigma^2\) pode ser expressa em termos da largura do intervalo de confiança dado na linha 1 da Tabela 5.3.

Como \(\mbox{E}(ssE) = t(r−1)\sigma^2\), a largura esperada do intervalo de confiança de 95% pode ser escrita como: \[\begin{equation*} \tag{5.8} \sigma^2\left( \dfrac{t(r-1)\times \Big( \chi^2_{t(r-1),0.975}-\chi^2_{t(r-1),0.025}\Big)}{\chi^2_{t(r-1),0.975}\times \chi^2_{t(r-1),0.025}}\right)\cdot \end{equation*}\]

Portanto, se você deseja que a meia largura do intervalo de confiança seja 50% de \(\sigma^2\), procure o número de níveis do fator aleatório \(t\), e o número de repetições \(r\), de modo que o multiplicador de \(\sigma^2\) na Equação (5.8) seja 1,0.

Isso pode ser feito usando a função qchisq no R enumerando vários casos. O exemplo abaixo mostra o cálculo do multiplicador de \(\sigma^2\) que determina a largura esperada do intervalo de confiança. A saída resultante abaixo do código mostra que qualquer combinação de \(t\) e \(r\) que resulte em \(\nu_2 = t(r−1)\) no intervalo de 36 a 38 dará a precisão desejada na estimativa de \(\sigma^2\).

nu2 = 36:44
chiu = qchisq(.975, nu2)
chil = qchisq(.025, nu2)
width = nu2 * (chiu - chil) / (chil * chiu)
halfw = width/2
data.frame(nu2, width, halfw)
##   nu2     width     halfw
## 1  36 1.0259871 0.5129936
## 2  37 1.0091269 0.5045635
## 3  38 0.9930584 0.4965292
## 4  39 0.9777224 0.4888612
## 5  40 0.9630653 0.4815327
## 6  41 0.9490392 0.4745196
## 7  42 0.9356004 0.4678002
## 8  43 0.9227095 0.4613548
## 9  44 0.9103307 0.4551654


Uma regra prática simples pode ser usada para se ter uma ideia de quantos níveis do fator aleatório \(t\), incluir no experimento de amostragem. Quando se espera que \(\sigma_t^2\) seja maior que \(\sigma^2\), \(t\) deve ser o maior possível, então \(t = \nu_2\) e \(r = 2\) seriam razoáveis.

Outra forma de determinar \(t\) e \(r\) seria considerar o poder do teste \(F\) para testar a hipótese \(H_0 ∶ \sigma_t^2 = 0\). Sob a hipótese alternativa \(H_1 ∶ \sigma^2_t > 0\), a estatística \(F = msT/msE\) segue um múltiplo de a distribuição \(F\) central e o poder ou probabilidade de exceder o limite crítico podem ser expressos por \[\begin{equation*} \tag{5.9} 1-\beta = P\left( F_{t-1,t(r-1)} > \dfrac{1}{1+r\times \rho}F_{t-1,t(r-1),\alpha} \right), \end{equation*}\] onde \(\rho=\sigma_t^2/\sigma^2\).

Novamente, os tamanhos de amostra \(t\) e \(r\) que fornecem potência adequada para alternativas especificadas podem ser determinados enumerando vários casos em R. Nesse caso, o uso das funções R qf e pf facilita isso. Por exemplo, se você quiser ter um poder maior que \(1-\beta = 0.90\) para rejeitar \(H_0 ∶ \sigma_t^2 = 0\), quando \(\rho = \sigma_t^2/\sigma^2\geq 3.0\) o código R abaixo encontrará algumas alternativas.

A saída resultante abaixo do código mostra várias combinações de \(t\) e \(r\) que resultam em potência superior a 0.90.

alpha = 0.05
rho = 3.0
t = rep(5:7, each = 3)
r = rep(2:4, 3)
nu_1 = t-1
nu_2 = t * (r - 1)
fcrit = qf( 1 - alpha, nu_1, nu_2 )
factor = 1 / ( 1 + r * rho )
plimit = factor * fcrit
power = 1 - pf( plimit, nu_1, nu_2 )
data.frame( t, r, power)
##   t r     power
## 1 5 2 0.6025330
## 2 5 3 0.8397523
## 3 5 4 0.9142402
## 4 6 2 0.6876308
## 5 6 3 0.8972133
## 6 6 4 0.9523702
## 7 7 2 0.7565926
## 8 7 3 0.9346005
## 9 7 4 0.9737459


É claro que este método não considera a precisão da estimativa de \(\sigma^2\). Pode-se considerar usar o primeiro método mostrado acima para determinar que \(\nu_2 = t(r-1)\) tenha a precisão desejada na estimativa de \(\sigma^2\); e então, com \(\nu_2 = t(r-1)\) fixado no número determinado, use o segundo método mostrado acima para determinar quão grande \(t\) deve ser para ter o poder adequado na rejeição de \(H_0 ∶ \sigma_t^2 = 0\), em um nível de significância especificado \(\alpha\), em favor de \(H_1 ∶ \sigma_t^2 > 0\) quando \(\sigma_t^2/\sigma^2\geq \rho\).


5.5 Projetos de amostragem de dois fatores


Quando o objetivo da experimentação é estudar a variância na resposta causada pela variação dos níveis de dois fatores independentes, o planejamento é semelhante ao fatorial de dois fatores apresentado no Capítulo 3.

Entretanto, nos experimentos fatoriais de dois fatores apresentados no Capítulo 3, os níveis dos fatores seriam selecionados especificamente pelo experimentador porque ele estaria interessado em comparar a resposta média entre esses níveis. No experimento de amostragem de dois fatores, por outro lado, os níveis dos fatores são apenas uma amostra aleatória ou representativa de níveis possíveis e o objetivo é determinar quanto da variância na resposta pode ser atribuída a níveis variados dos fatores. Em geral, esses projetos são chamados de experimentos de amostragem aleatória fatorial ou FRSE.

Considere os dados de exemplo apresentados na Tabela 5.5 abaixo retirados de Sower et al. (1999). Estes são dados de um estudo Gage R&R comumente realizado em departamentos de garantia de qualidade industrial.

\[ \begin{array}{cccc}\hline \mbox{Parte} \, \setminus \mbox{Operador} & 1 & 2 & 3 \\\hline 1 & 0.71 & 0.56 & 0.52 \\ & 0.69 & 0.57 & 0.54 \\\hline 2 & 0.98 & 1.03 & 1.04 \\ & 1.00 & 0.96 & 1.01 \\\hline 3 & 0.77 & 0.76 & 0.81 \\ & 0.77 & 0.76 & 0.81 \\\hline 4 & 0.86 & 0.82 & 0.82 \\ & 0.94 & 0.78 & 0.82 \\\hline 5 & 0.51 & 0.42 & 0.46 \\ & 0.51 & 0.42 & 0.49 \\\hline 6 & 0.71 & 1.00 & 1.04 \\ & 0.59 & 1.04 & 1.00 \\\hline 7 & 0.96 & 0.94 & 0.97 \\ & 0.96 & 0.91 & 0.95 \\\hline 8 & 0.86 & 0.72 & 0.78 \\ & 0.86 & 0.74 & 0.78 \\\hline 9 & 0.96 & 0.97 & 0.84 \\ & 0.96 & 0.94 & 0.81 \\\hline 10 & 0.64 & 0.56 & 1.01 \\ & 0.72 & 0.52 & 1.01 \\\hline \end{array} \] Tabela 5.5: Dados do estudo Gage R&R.

Nestes estudos o objetivo é classificar a variabilidade nas características medidas de produtos fabricados ou componentes de produtos. Supondo que o medidor ou instrumento de medição esteja devidamente calibrado, um valor medido determinado durante uma inspeção de controle de qualidade pode ser considerado uma função da verdadeira dimensão do recurso, da repetibilidade e da reprodutibilidade do medidor.

A repetibilidade do medidor é a capacidade de um único operador obter o mesmo valor de medição várias vezes usando o mesmo instrumento de medição ou medidor na mesma característica de um único componente ou peça fabricada. A reprodutibilidade do medidor é a capacidade de diferentes operadores obterem o mesmo valor medido várias vezes usando o mesmo medidor na mesma peça. Se a variabilidade nas medições causada pela repetibilidade do medidor mais a reprodutibilidade do medidor for superior a 10% da faixa de tolerância, as medições podem não ser precisas o suficiente para serem usadas no monitoramento da qualidade do produto.

O experimento de amostragem Gage R&R consiste em selecionar um conjunto de peças ou componentes fabricados que sejam representativos da variabilidade peça a peça na fabricação normal. Na Tabela 5.5 foi selecionada uma amostra de 10 partes e essas partes representam os níveis do primeiro fator no experimento de amostragem. Em seguida, é selecionada uma amostra aleatória ou representativa de inspetores.

Os inspetores ou operadores representam os níveis do segundo fator no experimento de amostragem. Finalmente, cada inspetor mede cada peça duas vezes numa ordem aleatória e os resultados são reunidos numa tabela como a Tabela 5.5. As medições replicadas representam as réplicas em cada célula.

Como cada operador ou inspetor mediu cada peça, o modelo para os dados da Tabela 5.5 pode ser escrito na forma de um modelo fatorial \[\begin{equation*} \tag{5.10} y_{ijk} = \mu + a_i + b_j + {ab}_{ij} + \epsilon_{ijk}, \end{equation*}\] onde \(y_{ijk}\) é a \(k\)-ésima medição, \(k = 1,\cdots,r\) feita pelo \(j\)-ésimo operador, \(j = 1,\cdots,b\) na \(i\)-ésima parte \(i = 1,\cdots,a\), \(a_i\) é o efeito parcial, \(b_j\) é o efeito operador ou inspetor e \({ab}_{ij}\) é o efeito de interação.

A diferença entre este modelo e o modelo (3.2) no Capítulo 3 é o fato de que os efeitos \(a_i\), \(b_j\) e \({ab}_{ij}\) são agora assumidos como variáveis aleatórias independentes e normalmente distribuídas com médias zero e variâncias \(\sigma_a^2\), \(\sigma_b^2\) e \(\sigma_{ab}^2\).

Como a variância de uma soma de variáveis aleatórias independentes é a soma das variâncias, a variância total na resposta \[ \mbox{Var}(y)= \sigma_y^2 = \sigma_a^2 + \sigma_b^2 + \sigma_{ab}^2 + \sigma^2\cdot \]

\(\sigma_a^2\) representa a parcela da variância total devido às diferenças reais nas características da peça, \(\sigma_b^2\) é a parcela da variância causada por às diferenças entre operadores, \(\sigma_{ab}^2\) é a parcela da variância causada pela interação do operador e da peça e \(\sigma^2\) é a parte da variação causada por medições replicadas ou repetibilidade do medidor.

A soma de \(\sigma_b^2 + \sigma_{ab}^2\) é a reprodutibilidade. A repetibilidade mais reprodutibilidade \(\sigma_b^2 + \sigma_{ab}^2 + \sigma^2\) é uma medida da variância atribuível ao erro de medição.


5.5.1 Estimando componentes de variância


Para o caso com igual número de réplicas por subclasse, como no exemplo da Tabela 5.5, é conveniente utilizar o método dos momentos ou estimadores REML dos componentes de variância. Ao utilizar o método dos momentos, é necessário conhecer os coeficientes dos componentes da variância nos quadrados médios esperados.

Para um planejamento amostral equilibrado de dois fatores, o valor esperado dos quadrados médios pode ser encontrado usando o método tabular de Bennett and Franklin (1954), mostrado na Tabela 5.6. As estimativas dos componentes de variância podem então ser obtidas igualando os quadrados médios na ANOVA aos quadrados médios esperados e resolvendo os componentes de variância.

\[ \begin{array}{ccc} \hline \mbox{Fonte} & \mbox{df} & \mbox{Esperança dos quadrados médios} \\\hline \mbox{A} & a-1 & \sigma^2+r\, \sigma^2_{AB}+rb \, \sigma^2_A \\ \mbox{B} & b-1 & \sigma^2+r\, \sigma^2_{AB}+ra \, \sigma^2_B \\ \mbox{AB} & (a-1)(b-1) & \sigma^2+r\, \sigma^2_{AB} \\ \mbox{Erro} & (r-1)ab & \sigma^2 \\\hline \end{array} \] Tabela 5.6: Quadrados médios esperados em projeto de amostragem de dois fatores.

O código R para produzir a tabela ANOVA é mostrado abaixo.

library(daewr)
modr1 = aov( y ~ part + oper + part:oper, data = gagerr)
summary(modr1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## part         9 1.4489 0.16099  214.18  < 2e-16 ***
## oper         2 0.0297 0.01485   19.76 3.35e-06 ***
## part:oper   18 0.4839 0.02688   35.77 1.87e-15 ***
## Residuals   30 0.0225 0.00075                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Usando os quadrados médios na tabela acima, o código R abaixo calcula as estimativas dos componentes de variância usando o método dos momentos.

sigma2 = 0.000752
sigma2po = (0.026885 - sigma2) / 2
sigma2o = (0.014852 - sigma2 - 2 * sigma2po ) / 20
sigma2p = (.160991 - sigma2 - 2 * sigma2po ) / 6
cat("Method of Moments Variance Component Estimates","\n", 
    "Var(error)=",sigma2,"\n","Var(part x oper)=",sigma2po,"\n", 
    "Var(oper)=",sigma2o,"\n","Var(part)=",sigma2p,"\n")
## Method of Moments Variance Component Estimates 
##  Var(error)= 0.000752 
##  Var(part x oper)= 0.0130665 
##  Var(oper)= -0.00060165 
##  Var(part)= 0.022351


Para este exemplo o estimador do método dos momentos de \(\sigma_b^2\) é negativo. Uma ou duas observações atípicas podem resultar em um componente de variância estimado muito grande, muito pequeno ou até negativo. Técnicas gráficas, como as descritas na Seção 5.9, muitas vezes podem revelar valores atípicos e aumentar a chance de uma interpretação útil dos dados.

Nos casos em que valores atípicos não são causa de estimativas negativas de componentes de variância, os estimadores REML evitam os valores negativos. O código R para produzir as estimativas REML é mostrado abaixo. A partir desses resultados pode-se observar que a maior parte de \[ 94.3\% = 100\times \big(0.01247/(0.01247 + 0.0007517)\big) \] do erro de medição se deve à reprodutibilidade.

Portanto, para reduzir o erro de medição, os esforços devem concentrar-se numa melhor formação dos inspetores operadores, em vez de investir em medidores mais precisos.

library(lme4)
modr2 = lmer(y ~ 1 + (1|part) + (1|oper) + (1|part:oper), data = gagerr)
summary(modr2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ 1 + (1 | part) + (1 | oper) + (1 | part:oper)
##    Data: gagerr
## 
## REML criterion at convergence: -133.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.43502 -0.36558 -0.01169  0.38978  1.94190 
## 
## Random effects:
##  Groups    Name        Variance  Std.Dev.
##  part:oper (Intercept) 0.0124650 0.11165 
##  part      (Intercept) 0.0225515 0.15017 
##  oper      (Intercept) 0.0000000 0.00000 
##  Residual              0.0007517 0.02742 
## Number of obs: 60, groups:  part:oper, 30; part, 10; oper, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   0.7982     0.0518   15.41
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')



5.5.2 Intervalos de confiança nos componentes de variância em experimentos de dois fatores


Quando há um número igual de réplicas por subclasse e as suposições de normalidade são válidas, os quadrados médios da ANOVA são distribuídos independentemente como múltiplos de variáveis aleatórias qui-quadrado. Portanto, os intervalos de confiança em qualquer quadrado médio esperado podem ser obtidos de forma semelhante à linha 1 da Tabela 5.3. Contudo, exceto para \(\mbox{E}(msE) = \sigma^2\), o valor esperado de todos os outros quadrados médios são combinações lineares de dois ou mais componentes de variância.

Embora intervalos de confiança exatos possam ser obtidos em componentes de variância individuais no experimento balanceado de dois fatores seguindo a fórmula na linha 2 da Tabela 5.3, eles não são aplicáveis a todos os experimentos ou dados não balanceados. Burdick and Graybill (1992) mostram um método de cálculo de intervalos de confiança aproximados que é de aplicação mais geral.

Sempre que um componente de variância puder ser expresso na forma \[ \delta = c_1 \mbox{E}(ms_1) − c_2 \mbox{E}(ms_2), \] onde \(ms_1\) e \(ms_2\) são quadrados médios na tabela ANOVA e \(c_1\) e \(c_2\) são positivos, o intervalo de confiança aproximado mostrado por Burdick and Graybill é aplicável.

Por exemplo, na tabela ANOVA mostrada acima para o estudo de medição R&R, \[ \mbox{E}(ms \, \mbox{part:oper}) = \sigma^2+ 2\sigma^2_{ab} \] e \[ \mbox{E}(ms \, \mbox{oper}) = \sigma^2+ 2\sigma^2_{ab} + 20\sigma^2_ b, \] portanto \[ \delta = 0.05 \times \mbox{E}(ms \, \mbox{oper}) - 0.05 \times \mbox{E}(ms \, \mbox{part:oper}) = \sigma^2_b\cdot \]

Um intervalo de confiança aproximado \(1-\alpha\) em \(\delta\) é dado por: \[\begin{equation*} \tag{5.11} \Big( \widehat{\delta}-\sqrt{V_L}, \; \widehat{\delta}+\sqrt{V_U} \Big), \end{equation*}\] onde \(\widehat{\delta} = c_1 ms_1 − c_2 ms_2\), \(\nu_1\) são os graus de liberdade para \(ms_1\), \(\nu_2\) são os graus de liberdade para \(ms_2\), \[ \begin{array}{rcl} V_L & = & G_1^2 c_1^2 ms_1^2+H_2^2 c_2^2 ms_2^2 +G_{12} c_1c_2 (ms_1)(ms_2), \\ V_U & = & H_1^2 c_1^2 ms_1^2+G_2^2 c_2^2 ms_2^2 +H_{12} c_1c_2 (ms_1)(ms_2), \\ \end{array} \] \[ G_1=1-\dfrac{1}{F_{\alpha,\nu_1,\infty}}, \qquad H_2 = \dfrac{1}{F_{1-\alpha,\nu_2,\infty}}-1, \qquad G_{12} = \dfrac{\big(F_{\alpha,\nu_1,\nu_2}-1\big)^2 -G_1^2 F^2_{\alpha,\nu_1,\nu_2} -H_2^2}{F_{\alpha,\nu_1,\nu_2}}, \] \[ H_1=\dfrac{1}{F_{1-\alpha,\nu_1,\alpha}}-1, \qquad G_2=1-\dfrac{1}{F_{\alpha,\nu_2,\alpha}} \quad \mbox{e} \quad H_{12} = \dfrac{\big(1-F_{1-\alpha,\nu_1,\nu_2}\big)^2 -H_1^2 F^2_{\alpha,\nu_1,\nu_2} -H_2^2}{F_{\alpha,\nu_1,\nu_2}}\cdot \]

Embora essas fórmulas pareçam formidáveis, elas podem ser avaliadas usando R. Uma função vci, no pacote daewr, calculará intervalos de confiança usando essas fórmulas. As entradas para esta função são confl = \(1-\alpha\), c1 = \(c_1\) , ms1 = \(ms_1\), nu1 = \(\nu_1\), c2 = \(c_2\), ms2 = \(ms_2\) e nu2 = \(\nu_2\).

Como ilustração do uso desta função, considere criar intervalos de confiança no componente de variância para o operador e no componente de variância para a interação parte por operador nos dados do estudo de medição R&R.

Para obter um intervalo de confiança no componente de variância para o operador \(\sigma_b^2\), \(ms_1\) = 0.01485 com 2 graus de liberdade, \(ms_2\) = 0.02689 com 18 graus de liberdade. A chamada da função e a saída são mostradas abaixo.

library(daewr)
options(digits = 3)
vci(confl = .90, c1 = .05, ms1 = .01485, nu1 = 2, c2 = .05, ms2 = .02689, nu2 = 18)
## delta= -0.000602  Lower Limit= -0.00158  Upper Limit= 0.00572


Isso mostra que embora o estimador do método dos momentos para \(\sigma_b^2 = −0.000602\) seja negativo, o limite superior de confiança de 90% é positivo. O intervalo de confiança de 90% para \(\sigma_{ab}^2\) é (0.008936 - 0.021895) usando a mesma fórmula de aproximação.


5.5.3 Determinando tamanhos de amostra para experimentos de amostragem de dois fatores


Em um experimento de amostragem de dois fatores, há três tamanhos de amostra a serem considerados. Primeiro, o número de níveis do fator A \(a\), segundo, o número de níveis do fator B \(b\) e, finalmente, o número de réplicas dentro de cada célula, \(r\). Os graus de liberdade para o quadrado médio replicado são \(ab(r-1)\).

Ao substituir \(t(r-1)\) por \(ab(r-1)\) na fórmula (5.8), pode ser usado para determinar o valor de \(ab(r-1)\) que resultará na largura desejada do intervalo de confiança para \(\sigma^2\). Em seguida, utilizando uma regra prática como a expressa na Secção 5.4.4, os níveis do factor A e do factor B podem ser determinados.


5.5.4 Estudos de dois fatores com replicação desigual


Quando há um número desigual de réplicas em cada subclasse, o método dos momentos ou estimadores de análise de variância dos componentes de variância são mais difíceis de calcular. Os coeficientes dos componentes de variância nos quadrados médios esperados não são mais funções inteiras do número de níveis dos fatores e do número de repetições, conforme mostrado na Tabela 5.6. Em vez disso, devem ser calculados utilizando o método de Hartley conforme descrito na Seção 5.4.1.

Além disso, as estimativas dos componentes de variância obtidas resolvendo as equações que igualam os quadrados médios aos seus valores esperados não são únicas. As estimativas serão diferentes dependendo se as somas dos quadrados sequenciais ou do tipo III são usadas nas equações. Por outro lado, os estimadores de máxima verossimilhança e os estimadores REML são únicos e podem ser preferidos nesta situação.

Como exemplo de cálculo dos componentes de variância usando o método REML, considere os dados da Tabela 5.7 retirados de um estudo de amostragem para estimar as fontes de variabilidade em um ensaio interlaboratorial de cálcio no soro sanguíneo que foi mostrado por Rao and Rao (1997).

\[ \begin{array}{ccccc}\hline \mbox{Laboratório} \, \setminus \mbox{Solulçao padrão} & 1 & 2 & 3 & 4 \\\hline \mbox{A} & 87 & 92 & 179 & 177 \\ & & 84 & 83 & 173 \\ & & 80 & 76 & 166 \\ & & & & \\ \mbox{B} & 80 & 69 & 138 & 151 \\ & & 70 & 46 & 138 \\ & & & & 132 \\ & & & & \\ \mbox{C} & 70 & 67 & 173 & 176 \\ & & 60 & 63 & 166 \\ & & 44 & 48 & \\ \hline \end{array} \] Tabela 5.7: Cálcio em soluções de soro sanguíneo com concentrações desconhecidas.

O código para abrir o quadro de dados e calcular os componentes de variância usando a função lmer no pacote lme4 é semelhante àqueles mostrados na última seção para dados balanceados. O código e os resultados da análise são mostrados na próxima página. O componente de variância estimado para a interação laboratório por solução é próximo de zero, mas o estimador do método dos momentos do mesmo componente seria negativo.

Pode-se observar que a maior parte da variabilidade no ensaio de cálcio no soro sanguíneo se deve a diferenças entre as soluções padrão e entre as análises repetidas utilizando a mesma solução no mesmo laboratório. Uma proporção muito pequena da variabilidade é causada por diferenças nos laboratórios.

library(daewr)
library(lme4)
rmod3 = lmer( calcium ~ 1 + (1|lab) + (1|sol) + (1|lab:sol), data = blood)
summary(rmod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: calcium ~ 1 + (1 | lab) + (1 | sol) + (1 | lab:sol)
##    Data: blood
## 
## REML criterion at convergence: 265
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6405 -0.5553 -0.0612  0.4173  2.3327 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  lab:sol  (Intercept) 1.27e-07 3.56e-04
##  sol      (Intercept) 1.49e+03 3.86e+01
##  lab      (Intercept) 2.80e+01 5.29e+00
##  Residual             1.05e+03 3.24e+01
## Number of obs: 27, groups:  lab:sol, 12; sol, 4; lab, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    103.2       20.7    4.99
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')


O código para abrir o conjunto de dados e calcular os componentes de variância usando a função lmer no pacote lme4 é semelhante àqueles mostrados na última seção para dados balanceados. O componente de variância estimado para a interação laboratório por solução é próximo de zero, mas o estimador do método dos momentos do mesmo componente seria negativo.

Pode-se observar que a maior parte da variabilidade no ensaio de cálcio no soro sanguíneo se deve a diferenças entre as soluções padrão e entre as análises repetidas utilizando a mesma solução no mesmo laboratório. Uma proporção muito pequena da variabilidade é causada por diferenças nos laboratórios.

Como existem apenas três laboratórios, qualquer intervalo de confiança assintótico para \(\sigma^2_a\) não seria preciso. Com replicação desigual, os quadrados médios da ANOVA não são mais distribuídos como múltiplos de variáveis aleatórias qui-quadradas, portanto os intervalos de confiança que usam a linha 2 na Tabela 5.3 também não são apropriados.

Intervalos de confiança aproximados ainda podem ser obtidos se os quadrados médios não ponderados forem substituídos pelos quadrados médios na Fórmula (5.11). Os quadrados médios não ponderados para os fatores A, B e a interação podem ser obtidos realizando uma ANOVA nas médias das células. A Tabela 5.8 é uma representação simbólica da ANOVA com média quadrada não ponderada e seus valores esperados de um planejamento amostral de dois fatores.

\[ \begin{array}{lccl}\hline \mbox{Fonte} & \mbox{df} & \mbox{Quadrados médios} & \mbox{Quadrados médios esperados} \\\hline \mbox{Fator A} & a-1 & msA_U & \sigma^2 + \widetilde{c}\sigma^2_{ab} + b \, \widetilde{c}\sigma^2_a \\ \mbox{Fator A} & b-1 & msB_U & \sigma^2 + \widetilde{c}\sigma^2_{ab} + a\, \widetilde{c}\sigma^2_b \\ \mbox{Interação AB} & (a-1)(b-1) & msAB_U & \sigma^2 + \widetilde{c}\sigma^2_{ab} \\ \mbox{Erro} & \sum_i\sum_j r_{ij}-ab & msE & \sigma^2\\\hline \end{array} \] Tabela 5.8: ANOVA simbólica com quadrados médios não ponderados.

O símbolo \(\widetilde{c}=ab/\sum_i\sum_j (1/r_{ij})\) e \(r_{ij}\) é o número de repetições na \(ij\)-ésima célula. A partir desta tabela pode-se ver que \[ \widehat{\sigma}_a^2 =\dfrac{1}{b\widetilde{c}}msA_U - \dfrac{1}{b\widetilde{c}}msAB_U, \] o qual é da forma \[ \delta = c_1\mbox{E}(ms_1) - c_2\mbox{E}(ms_2), \] mostrada na Seção 5.5.2.

Para os dados da Tabela 5.7, os quadrados médios não ponderados são obtidos pelo código R abaixo. Nesse código, a primeira instrução usa a função tapplycalcium no sangue do conjunto de dados.

Em seguida, são criados fatores com o mesmo comprimento do vetor de médias das células e a função aov é usada para fazer a tabela ANOVA de médias não ponderadas.

cellmeans = tapply( blood$calcium, list(blood$lab, blood$sol), mean)
dim(cellmeans) = NULL
Lab = factor( rep(c("A", "B", "C"), 4))
Solution = factor(rep( c( 1, 2, 3, 4), each = 3))
mod2 = aov( cellmeans ~ Lab + Solution + Lab:Solution )
summary(mod2)
##              Df Sum Sq Mean Sq
## Lab           2    826     413
## Solution      3  15035    5012
## Lab:Solution  6    625     104


Nos resultados, a média quadrada não ponderada para o fator A, \(maA_U = 413\) e a média quadrada não ponderada para o fator B, \(msAB_U = 104\). O fator \[ \begin{array}{rcl} \widetilde{c} & = & \dfrac{ab}{\sum_i\sum_j (1/r_{ij})} \\ & = & \dfrac{(3)(4)}{\Big( \frac{1}{1} + \frac{1}{3} + \frac{1}{3} + \frac{1}{3} + \frac{1}{1} + \frac{1}{2} + \frac{1}{2} + \frac{1}{3} + \frac{1}{1} + \frac{1}{3} + \frac{1}{3} + \frac{1}{2} \Big)} \, = \, 1.846, \end{array} \] e \[ c_1 \, = \, c_2 \, = \, \dfrac{1}{b\widetilde{c}} = \dfrac{1}{4(1.846)} \, = \, 0.13541\cdot \]

Assim, as entradas necessárias para a função vci calcular um intervalo de confiança de 90% para \(\sigma_a^2\) são mostradas no código abaixo.

library(daewr)
vci(confl = .90, c1 = .1354166, ms1 = 413, nu1 = 2, c2 = .1354166, ms2 = 104, nu2 = 6)
## delta= 41.8  Lower Limit= 4.14  Upper Limit= 517

e o intervalo de confiança resultante é (4.138, 516.772) muito amplo, novamente devido ao fato de haver apenas três laboratórios nos dados.


5.6 Experimentos de amostragem aninhada (NSE)


Muitos experimentos de amostragem com mais de um fator usam fatores aninhados ou um projeto hierárquico. Os níveis de um fator aninhado são fisicamente diferentes dependendo do nível do fator no qual ele está aninhado.

Esse não foi o caso dos fatores descritos na última seção. Por exemplo, no estudo de medição R&R, cada operador mediu cada peça; portanto, o número do operador foi definido exclusivamente e referia-se ao mesmo operador, independentemente de qual parte ele mediu. Chamaríamos o operador e a parte do estudo de medição R&R de fatores cruzados.

Para alterar o plano de modo que o operador seja um fator aninhado, considere um experimento onde \(n\) partes foram selecionadas e cada parte foi medida por dois operadores, embora não precisem ser os mesmos dois operadores medindo cada parte. Esta pode ser uma maneira mais conveniente de conduzir o experimento de amostragem se as peças forem selecionadas durante um longo período de tempo e os mesmos operadores nem sempre estiverem presentes para fazer as medições.

O fato de os operadores diferirem dependendo do número da peça que está sendo medida torna o operador um fator aninhado, aninhado dentro da peça. O primeiro operador que mede a primeira peça não é fisicamente a mesma pessoa que o primeiro operador que mede a peça subsequente.

Outro exemplo comum onde ocorre um projeto aninhado é quando as medições são destrutivas. Nesse caso, cada operador deve medir um conjunto diferente de peças e a parte torna-se o fator aninhado, aninhado dentro do operador, porque a primeira parte medida pelo primeiro operador não é fisicamente igual à primeira parte medida pelos operadores subsequentes.

Este tipo de plano amostral é chamado de experimento de amostragem aninhada ou NSE. Um exemplo onde já vimos fatores aninhados está no termo \(\epsilon_{ij}\) nos modelos que usamos até agora. Ele representa o efeito da \(j\)-ésima unidade experimental replicada e, como diferentes unidades experimentais são usadas para cada nível de fator ou combinação de níveis de fator, a unidade experimental é sempre aninhada em outro nível de fator ou célula no experimento.

Quando dois fatores são fatores cruzados, podemos incluir sua interação no modelo para representar a variabilidade extra causada por mudanças no nível de um fator entre os níveis do outro fator. No entanto, se um fator estiver aninhado dentro de outro fator, não podemos incluir uma interação entre eles porque o fator aninhado inclui os graus de liberdade que poderiam ser assumidos pela interação.

O modelo para um experimento aninhado de dois estágios com fator B aninhado dentro do fator A é escrito como: \[\begin{equation*} \tag{5.12} y_{ijk} = \mu + a_i + b_{(i)j} + \epsilon_{ijk}, \end{equation*}\] e se houver um número igual de réplicas \(r\), por célula, a tabela ANOVA para o modelo aninhado pode ser representada como: \[ \begin{array}{cccc} \hline \mbox{Fonte} & \mbox{df} & \mbox{Quadrados médios} & \mbox{Quadrados médios esperados} \\\hline \mbox{Fator A} & a-1 & msA & \sigma^2+ r\, \sigma_b^2 + br\, \sigma_a^2 \\ \mbox{Fator B} & a(b-1) & msB & \sigma^2+ r\, \sigma_b^2 \\ \mbox{Erro} & ab(r-1) & msE & \sigma^2 \\\hline \end{array} \] Tabela 5.9: ANOVA simbólica para experimento aninhado de dois fatores.

Aqui pode ser visto que os graus de liberdade para o fator B, \[ a(b-1) \, = \, (b-1)+(a-1)(b-1) \] é igual aos graus de liberdade de um fator cruzado mais os graus da interação AB. Os quadrados médios esperados na tabela podem novamente ser obtidos usando o método tabular de Bennett and Franklin (1954).

Projetos aninhados ou hierárquicos podem ser estendidos para incluir vários estágios ou fatores. Por exemplo, a Tabela 5.10 mostra os resultados de um estudo de amostragem agrupada em quatro estágios sobre a variabilidade das propriedades da borracha bruta, extraído de Bennett and Franklin (1954).

\[ \begin{array}{c|cc|cc|cc|cc}\hline \mbox{Fornecedor} & \mbox{A} & & \mbox{B} & & \mbox{C} & & \mbox{D} & \\\hline \mbox{Misturas de amostras} & 1 & 2 & 1 & 2 & 1 & 2 & 1 & 2 \\\hline \mbox{Lote I} & 211 & 171 & 196 & 196 & 200 & 240 & 323 & 262 \\ & 215 & 198 & 186 & 210 & 221 & 229 & 279 & 234 \\ & 197 & 268 & 190 & 156 & 198 & 217 & 251 & 249 \\\hline \mbox{Lote II} & 229 & 234 & 209 & 200 & 191 & 196 & 255 & 249 \\ & 196 & 210 & 193 & 186 & 189 & 198 & 235 & 247 \\ & 200 & 226 & 204 & 196 & 186 & 175 & 223 & 239 \\\hline \mbox{Lote III} & 204 & 225 & 204 & 174 & 211 & 196 & 228 & 262 \\ & 221 & 215 & 165 & 172 & 197 & 184 & 250 & 227 \\ & 238 & 196 & 194 & 171 & 210 & 190 & 260 & 272 \\\hline \mbox{Lote IV} & 229 & 248 & 198 & 202 & 196 & 180 & 273 & 273 \\ & 250 & 249 & 209 & 211 & 197 & 166 & 241 & 256 \\ & 238 & 249 & 221 & 204 & 186 & 172 & 221 & 230 \\\hline \end{array} \] Tabela 5.10: Módulo de elasticidade a 700% de alongamento de 96 amostras preparadas de folha de borracha defumada.

Neste estudo foi retirada uma amostra de quatro lotes de borracha de cada um dos quatro fornecedores. Como o primeiro lote obtido do primeiro fornecedor não é fisicamente igual ao primeiro lote obtido do segundo fornecedor, o lote é aninhado no fornecedor.

Em seguida, foram feitas duas misturas de amostras de cada lote e, como as duas misturas de amostras de um lote são fisicamente diferentes das misturas de amostras de qualquer outro lote, a mistura de amostras é aninhada dentro do lote. Finalmente, três réplicas de testes foram realizadas em cada mistura de amostra para determinar a elasticidade.

O modelo para os dados pode ser escrito como: \[ y_{ijkl}=\mu+a_i+b_{(i)j} + c_{(ij)k} + \epsilon_{ijkl} \] onde \(y_{ijkl}\) é a \(l\)-ésima determinação da elasticidade feita a partir da \(k\)-ésima mistura de amostras, retirada do \(j\)-ésimo lote do \(i\)-ésimo fornecedor, \(a_i\) é o efeito aleatório do fornecedor, \(b_{(i)j}\) é o efeito aleatório do lote, \(c_{(ij)k}\) é o efeito aleatório de mistura de amostra e \(\epsilon_{ijkl}\) é o efeito de determinação de replicação aleatória, \(i = 1,\cdots, 4\), \(j = 1,\cdots, 4\), \(k = 1, 2\) e \(i = 1,2, 3\).

Este modelo pode ser escrito na notação da função R aov como

mod2 = aov( elasticity ~ supplier + supplier:batch + supplier:batch:sample, data = rubber)
summary(mod2)
##                       Df Sum Sq Mean Sq F value  Pr(>F)    
## supplier               3  51990   17330   57.67 < 2e-16 ***
## supplier:batch        12  12735    1061    3.53 0.00048 ***
## supplier:batch:sample 16   5080     317    1.06 0.41365    
## Residuals             64  19233     301                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


ou na notação da função R lmer no pacote lme4 como

modr3 = lmer(elasticity ~ 1+ (1|supplier) + (1|supplier:batch) + (1|supplier:batch:sample), data = rubber)
summary(modr3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: elasticity ~ 1 + (1 | supplier) + (1 | supplier:batch) + (1 |  
##     supplier:batch:sample)
##    Data: rubber
## 
## REML criterion at convergence: 844
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.440 -0.659 -0.015  0.564  3.489 
## 
## Random effects:
##  Groups                Name        Variance Std.Dev.
##  supplier:batch:sample (Intercept)   5.66    2.38   
##  supplier:batch        (Intercept) 123.95   11.13   
##  supplier              (Intercept) 677.86   26.04   
##  Residual                          300.52   17.34   
## Number of obs: 96, groups:  
## supplier:batch:sample, 32; supplier:batch, 16; supplier, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    215.9       13.4    16.1


A notação ‘:’ indica que o lote (batch) está aninhado no fornecedor (supplier) e assim adiante. Os componentes de variância \(\sigma_a^2\), \(\sigma_b^2\), \(\sigma_c^2\) e \(\sigma^2\) podem ser estimados usando o método dos momentos ou REML.

Para aumentar a confiança nas estimativas dos componentes de variância, o número de níveis de um fator aleatório deve ser aumentado. Entretanto, em projetos hierárquicos com vários estágios, aumentar o número de níveis do fator mais alto aumenta muito o tamanho geral da amostra, mesmo que todos os outros fatores aninhados tenham apenas dois níveis. Por exemplo, no projeto mostrado na Tabela 5.10, se o número de fornecedores fosse aumentado de 4 para 20 para obter uma estimativa mais precisa de \(\sigma_a^2\), o número de determinações que teriam que ser feitas aumentaria dos 96 mostrados na Tabela 5.10 a 480.

Mesmo que o número de lotes por fornecedor e o número de misturas de amostras por lote e de determinações por mistura fossem reduzidos para 2 cada, ainda haveria \(20\times 2\times 2\times 2 = 160\) determinações. Se o estudo de amostragem tivesse sido feito desta forma haveria \(a-1 = 19\) graus de liberdade para o efeito fornecedor, \(a(b-1) = 20(2-1) = 20\) graus de liberdade para o efeito lote, \(ab(c-1) = 20\times 2 (2-1) = 40\) graus de liberdade para o efeito amostral e \(abc(r-1) = 20(2) (2) (2-1) = 80\) graus de liberdade para o efeito de replicação aleatória.

Portanto, a maioria das 160 observações são usadas para aumentar a precisão dos dois componentes de variância inferiores \(\sigma^2\) e \(\sigma_c^2\). Por esta razão, projetos hierárquicos balanceados geralmente não são recomendados se houver mais de três estágios ou fontes de variabilidade sendo estudadas. Os projetos aninhados escalonados apresentados na próxima seção permitem a conveniência dos fatores aninhados em estudos de amostragem, mas permitem que os vários componentes de variância sejam estimados com precisão mais uniforme.


5.7 Projetos aninhados escalonados


Experimentos de amostragem aninhada escalonada ou SNSE, foram desenvolvidos independentemente por Bainbridge (1965) e Prairie and Anderson (1962). Em um projeto completamente aninhado, conforme discutido na última seção, cada nível do fator superior leva a dois ou mais níveis de cada fator ou estágio seguinte.

Por outro lado, em um experimento aninhado escalonado, apenas um dos dois níveis do fator sucessor leva ao próximo estágio de dois níveis. A Figura 5.1 ilustra a diferença entre um design aninhado e um design aninhado escalonado para três estágios. Se houver \(a\) níveis do superior do fator, o experimento aninhado (nested design) requer \(4a\) observações no total, enquanto o experimento aninhado escalonado (staggered nested design) requer apenas \(3a\) observações.

Figura 5.1: Comparação entre experimentos aninhados (nested design) de 3 estágios e experimentos aninhados escalonados (staggered nested design) de 3 estágios.

A economia nas observações é multiplicada à medida que aumenta o número de estágios em um projeto aninhado. A Figura 5.2 mostra o esquema para projetos aninhados escalonados de três a seis estágios.

Figura 5.2: Experimentos aninhados escalonados para 3 a 6 estágios.

Embora as informações ou graus de liberdade disponíveis para estimar os componentes de variância em um experimento completamente aninhado estejam concentrados no nível inferior do fator, as informações são balanceadas em um experimento aninhado escalonado. A Tabela 5.11 compara a distribuição dos graus de liberdade entre um experimento aninhado escalonado e um experimento completamente aninhado, onde cada fator, exceto o nível mais alto, tem apenas dois níveis.

\[ \begin{array}{ccc}\hline & \mbox{Experimento} & \\ & \mbox{aninhado} & \mbox{Experimento} \\ & \mbox{escalonado} & \mbox{aninhado} \\ \mbox{Fonte} & \mbox{df} & \mbox{df} \\\hline \mbox{A} & a-1 & a-1 \\ \mbox{B en A} & a & a \\ \mbox{C en B} & a & 2a \\ \mbox{D en C} & a & 4a \\ \mbox{E en D} & a & 8a \\ \mbox{F en E} & a & 16a \\\hline \end{array} \] Tabela 5.11: Comparação de graus de liberdade entre experimentos aninhados e aninhados escalonados.

Embora os coeficientes quadráticos médios esperados para os experimentos aninhados escalonados possam ser determinados usando o método de Hartley (Hartley, 1967), seu método não está disponível em nenhuma função R. Felizmente, estes coeficientes foram tabulados por Nelson (1983) e são apresentados na Tabela 5.12. Eles são necessários no cálculo dos estimadores dos componentes de variância pelo método dos momentos.

\[ \begin{array}{ccl}\hline \mbox{Estágios} & \mbox{Termo} & \mbox{Quadrados médios esperados} \\\hline 3 & \mbox{A} & \sigma^2_C + (5/3)\sigma_B^2+ 2\sigma_A^2 \\ & \mbox{B} & \sigma^2_C + (4/3)\sigma_B^2 \\ & \mbox{C} & \sigma^2_C \\\hline 4 & \mbox{A} & \sigma^2_D + (3/2)\sigma_C^2+ (5/2)\sigma_B^2+ 4\sigma_A^2 \\ & \mbox{B} & \sigma^2_D + (7/6)\sigma_C^2+ (3/2)\sigma_B^2 \\ & \mbox{C} & \sigma^2_D + (4/3)\sigma_C^2 \\ & \mbox{D} & \sigma^2_D \\\hline 5 & \mbox{A} & \sigma^2_E + (7/5)\sigma_D^2+ (17/5)\sigma_B^2+ 5\sigma_A^2 \\ & \mbox{B} & \sigma^2_E + (11/10)\sigma_D^2+ (13/10)\sigma_C^2 +(8/5)\sigma_B^2 \\ & \mbox{C} & \sigma^2_E + (7/6)\sigma_D^2 +(3/2)\sigma_C^2 \\ & \mbox{D} & \sigma^2_E + (4/3)\sigma_D^2 \\ & \mbox{E} & \sigma_E^2 \\ \hline 6 & \mbox{A} & \sigma^2_F + (4/3)\sigma_E^2+ 2\sigma_D^2+ 3\sigma_C^2 + (13/12)\sigma_B^2 + 6\sigma_A^2\\ & \mbox{B} & \sigma^2_F + (14/15)\sigma_E^2+ (6/5)\sigma_D^2 +(7/5)\sigma_B^2 +(5/3)\sigma_B^2 \\ & \mbox{C} & \sigma^2_F + (11/10)\sigma_E^2 +(13/10)\sigma_D^2 +(8/5)\sigma_C^2 \\ & \mbox{D} & \sigma^2_F + (7/6)\sigma_E^2 + (3/2)\sigma_D^2 \\ & \mbox{E} & \sigma_F^2 + (4/3)\sigma_E^2 \\ & \mbox{F} & \sigma_F^2 \\\hline \end{array} \] Tabela 5.12: Coeficientes quadrados médios esperados para experimentos aninhados escalonados.

Mason et al. (1989) descreveram um estudo onde um projeto aninhado escalonado foi usado para estimar as fontes de variabilidade em um processo de polimerização contínua. Nesse processo são produzidos pellets de polietileno em lotes de cem mil libras.

Um projeto de quatro estágios foi utilizado para dividir a fonte de variabilidade na resistência à tração entre lotes, dentro dos lotes e devido ao processo de medição. Trinta lotes foram amostrados aleatoriamente. O lote representou o fator mais alto ou fonte de variabilidade A. De cada lote foram selecionadas aleatoriamente duas caixas de pellets.

Isto representou o segundo estágio ou fonte de variabilidade B. Da primeira caixa selecionada de cada lote foram feitas duas preparações para teste de resistência, mas da segunda caixa selecionada de cada lote foi feita apenas uma preparação. Isto representou o terceiro estágio ou fonte de variabilidade C. Finalmente, dois testes repetidos de resistência foram feitos a partir da primeira preparação da caixa um, enquanto apenas um teste de resistência foi feito a partir das outras três preparações.

Este esquema de amostragem está diagramado na Figura 5.3, e os dados são mostrados na Tabela 5.13 logo depois.

Figura 5.3: Diagrama do esquema de amostragem para estudo de polimerização.

\[ \begin{array}{c|ccc|c}\hline & & \mbox{Caixa 1} & & \mbox{Caixa 2} \\ & \mbox{Preparação} & \mbox{Preparação} & \mbox{Preparação} & \mbox{Preparação} \\ & 1 & 1 & 2 & 1 \\ \mbox{Lote} & \mbox{Teste 1} & \mbox{Teste 2} & \mbox{Teste 1} & \mbox{Teste 1} \\\hline 1 & 9.76 & 9.24 & 11.91 & 9.02 \\ 2 & 10.65 & 7.77 & 10.00 & 13.69 \\ 3 & 6.50 & 6.26 & 8.02 & 7.95 \\ 4 & 8.08 & 5.28 & 9.15 & 7.46 \\ 5 & 7.84 & 5.91 & 7.43 & 6.11 \\ 6 & 9.00 & 8.38 & 7.01 & 8.58 \\ 7 & 12.81 & 13.58 & 11.13 & 10.00 \\ 8 & 10.62 & 11.71 & 14.07 & 14.56 \\ 9 & 4.88 & 4.96 & 4.08 & 4.76 \\ 10 & 9.38 & 8.02 & 6.73 & 6.99 \\ 11 & 5.91 & 5.79 & 6.59 & 6.55 \\ 12 & 7.19 & 7.22 & 5.77 & 8.33 \\ 13 & 7.93 & 6.48 & 8.12 & 7.43 \\ 14 & 3.70 & 2.86 & 3.95 & 5.92 \\ 15 & 4.64 & 5.70 & 5.96 & 5.88 \\ 16 & 5.94 & 6.28 & 4.18 & 5.24 \\ 17 & 9.50 & 8.00 & 11.25 & 11.14 \\ 18 & 10.93 & 12.16 & 9.51 & 12.71 \\ 19 & 11.95 & 10.58 & 16.79 & 13.08 \\ 20 & 4.34 & 5.45 & 7.51 & 5.21 \\ 21 & 7.60 & 6.72 & 6.51 & 6.35 \\ 22 & 5.12 & 5.85 & 6.31 & 8.74 \\ 23 & 5.28 & 5.73 & 4.53 & 5.07 \\ 24 & 5.44 & 5.38 & 4.35 & 7.04 \\ 25 & 3.50 & 3.88 & 2.57 & 3.76 \\ 26 & 4.80 & 4.46 & 3.48 & 3.18 \\ 27 & 5.35 & 6.39 & 4.38 & 5.50 \\ 28 & 3.09 & 3.19 & 3.79 & 2.59 \\ 29 & 5.30 & 4.72 & 4.39 & 6.13 \\ 30 & 7.09 & 7.82 & 5.96 & 7.14 \\\hline \end{array} \] Tabela 5.13: Dados do estudo de variabilidade de resistência de polimerização.

O código R para abrir o conjunto de dados e calcular a tabela ANOVA para uso na estimativa dos componentes de variância pelo método dos momentos é mostrado abaixo.

library(daewr)
mod2 = aov(strength ~ lot + lot:box + lot:box:prep, data = polymer)
summary(mod2)
##              Df Sum Sq Mean Sq F value  Pr(>F)    
## lot          29    856   29.52   45.55 < 2e-16 ***
## lot:box      30     50    1.67    2.58 0.00577 ** 
## lot:box:prep 30     68    2.28    3.52 0.00046 ***
## Residuals    30     19    0.65                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Os estimadores do método dos momentos dos componentes de variância são encontrados igualando os quadrados médios na tabela ANOVA acima aos quadrados médios esperados na Tabela 5.12 e resolvendo. Os resultados são mostrados a continuação. Mais uma vez o procedimento do método dos momentos produz uma estimativa negativa para uma das fontes. \[ \begin{array}{rcl} \sigma_R^2 & = & 0.648 \\ \sigma_P^2 & = & (2.281-0.648)/(4/3) \, = \, 1.22475 \\ \sigma_B^2 & = & \big(1.670-(0.648+(7/6)1.22475) \big)/(3/2) \, = \, -0.27125 \\ \sigma_L^2 & = & \Big( 29.516 - \big(0.648 + (3/2)(1.22475)+(5/2)(-0.27125) \big) \Big)/4 \, = \, 6.92725 \end{array} \]

Para estimar os componentes de variância usando o método REML, a função lmer do pacote lme4 é usada conforme mostrado abaixo. Os resultados seguem o código.

O estimador REML do componente de variância para caixa dentro do lote é essencialmente zero e não o valor negativo obtido pelo método dos momentos.

library(lme4)
modr3 = lmer( strength ~ 1 + (1|lot) + (1|lot:box) + (1|lot:box:prep), data = polymer)
summary(modr3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: strength ~ 1 + (1 | lot) + (1 | lot:box) + (1 | lot:box:prep)
##    Data: polymer
## 
## REML criterion at convergence: 469
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1896 -0.4119 -0.0206  0.3826  1.7703 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  lot:box:prep (Intercept) 1.03e+00 1.014671
##  lot:box      (Intercept) 2.37e-08 0.000154
##  lot          (Intercept) 7.24e+00 2.691221
##  Residual                 6.57e-01 0.810433
## Number of obs: 120, groups:  lot:box:prep, 90; lot:box, 60; lot, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    7.221      0.509    14.2


A partir desses resultados, vemos que \[ 81\% = 100 \times (7.2427)/(7.2427+1.0296+0.6568)) \] da variação total é devida à variabilidade entre lotes, enquanto a variabilidade dentro do lote ou caixa a caixa é insignificante. Portanto, se os fabricantes quiserem diminuir a variabilidade na resistência à tração, deverão concentrar-se na redução da variação nos fatores influentes que mudam entre lotes.

Embora nenhum intervalo de confiança exato de forma fechada tenha sido desenvolvido para componentes de variância estimados a partir de dados em planejamentos aninhados escalonados, quando o número de níveis do fator mais alto é maior que 30, estimativas assintóticas aproximadas podem ser criadas usando o método da verossimilhança perfilada com o pacote R lme4 conforme discutido na Seção 5.4.3. As estimativas de intervalo bayesiano dos componentes de variância também são úteis para esses projetos (ver Lawson, 2008).

Para determinar o tamanho da amostra para um experimento aninhado escalonado, siga o procedimento descrito na Seção 5.4.4 para determinar o tamanho da amostra para estimar a variância replicada em um desenho amostral de um fator. Como os graus de liberdade para todos os fatores ou estágios no experimento aninhado escalonado são quase iguais, seguir este procedimento fornecerá aproximadamente a mesma precisão em todos os componentes de variância.


5.8 Experimentos com fatores fixos e aleatórios


Em alguns casos em que fatores de tratamento fixos estão sendo estudados em um projeto experimental como aqueles discutidos nos Capítulos 2 a 4, fatores aleatórios também são introduzidos no modelo pela forma como o experimento é conduzido.

Por exemplo, considere os dados da Tabela 5.14 que resultaram de uma experiência comparando diferentes formulações e métodos de aplicação de um pesticida nas folhas das plantas de algodão. O objetivo era aumentar a quantidade de pesticida ativo remanescente nas folhas do algodoeiro uma semana após a aplicação. O pesticida em estudo degrada-se à luz solar e um determinado aditivo na formulação retarda esse processo.

Diferentes técnicas de aplicação podem diferir na quantidade de pesticida aplicada nas folhas da planta. Os fatores de tratamento neste experimento foram duas formulações diferentes do pesticida e dois métodos de aplicação diferentes, resultando em um experimento fatorial \(2^2\).

A unidade experimental foi uma fileira de 20 pés de algodoeiro chamada parcela, porque era uma área conveniente dentro da qual a aplicação de pesticida poderia ser controlada. Oito parcelas foram selecionadas e duas foram designadas aleatoriamente para cada uma das quatro combinações de tratamento, resultando em duas repetições por combinação de tratamento. Uma semana após a aplicação, os experimentadores estavam prontos para determinar os resíduos de pesticidas remanescentes nas folhas das plantas.

No entanto, havia muito material vegetal em uma parcela inteira para ser enviado ao laboratório para análise. Para tanto, foram selecionadas de cada parcela duas amostras de folhas em quantidade conveniente para análise laboratorial de resíduos de agrotóxicos. Cada amostra foi enviada ao laboratório resultando nos dados mostrados na Tabela 5.14.

\[ \begin{array}{ccccc}\hline & \mbox{Técnica de} & & \mbox{Amostra} & \mbox{Amostra} \\ \mbox{Formulação} & \mbox{aplicação} & \mbox{Parcela} & 1 & 2 \\\hline \mbox{A} & 1 & 1 & 0.237 & 0.252 \\ \mbox{A} & 1 & 2 & 0.281 & 0.274 \\ \mbox{B} & 1 & 1 & 0.247 & 0.294 \\ \mbox{B} & 1 & 2 & 0.321 & 0.267 \\ \mbox{A} & 2 & 1 & 0.392 & 0.378 \\ \mbox{A} & 2 & 2 & 0.381 & 0.346 \\ \mbox{B} & 2 & 1 & 0.351 & 0.363 \\ \mbox{B} & 2 & 2 & 0.334 & 0.348 \\\hline \end{array} \] Tabela 5.14: Resíduos de pesticidas em plantas de algodão.

A formulação, a técnica de aplicação e sua interação são fatores fixos porque os experimentadores estavam interessados em comparar a resposta média entre os níveis desses fatores. A parcela, por outro lado, é um fator aleatório que representa diferenças nas unidades experimentais. Está aninhado nas combinações de formulação por técnica de aplicação.

Não há interesse em comparar unidades experimentais dentro de cada combinação de formulação e aplicação. Em vez disso, múltiplas parcelas por combinação de tratamento foram incluídas no projeto para que a variação causada por parcelas diferentes pudesse ser estimada e usada para avaliar a significância dos efeitos da formulação e da aplicação.

As amostras replicadas retiradas de cada parcela foram por conveniência na condução do experimento. Elas seriam classificadas como subamostras ou unidades observacionais definidas no Capítulo 1.

A forma mais simples de analisar os dados seria calcular a média das duas subamostras e proceder conforme ilustrado na Secção 3.5. Contudo, se as subamostras puderem ser consideradas independentes e for desejável incluir todos os dados (mostrados na Tabela 5.14 na análise, então um termo adicional para amostra deve ser incluído no modelo.

A amostra é outro efeito aleatório, pois não há interesse específico em comparar a resposta entre as duas amostras de cada parcela.

O modelo para os dados pode ser escrito na forma \[\begin{equation*} \tag{5.13} y_{ijkl} = \mu + \alpha_i + \beta_j + \alpha\beta_{ij} + p_{(ij)k} + \epsilon_{ijkl}, \end{equation*}\] onde \(y_{ijkl}\) é o resíduo de agrotóxico encontrado na \(l\)-ésima amostra retirada da \(k\)-ésima parcela, tratada com nível de formulação \(i\) e técnica de aplicação \(j\).

Em geral \(i = 1,\cdots,a\), \(j = 1,\cdots,b\), \(k = 1,\cdots,r\) e \(l = 1,\cdots,s\). Neste exemplo específico, \(a = 2\), \(b = 2\), \(r = 2\) e \(s = 2\), \(\alpha_i\) é o efeito da formulação, \(\beta_j\) é o efeito da aplicação, \(\alpha\beta_{ij}\) é o efeito da interação, \(p_{(ij)k}\) é o efeito aleatório da parcela e \(\epsilon_{ijkl}\) é o efeito aleatório de amostra.

O modelo pode ser escrito em notação matricial como: \[\begin{equation*} \tag{5.14} \pmb{y} = \pmb{X}\beta + \pmb{Z}\gamma + \epsilon, \end{equation*}\] onde \[ \pmb{y} = \begin{pmatrix} y_{1111} \\ y_{1112} \\ y_{1121} \\ y_{1122} \\ y_{2111} \\ y_{2112} \\ y_{2121} \\ y_{2122} \\ y_{1211} \\ y_{1212} \\ y_{1221} \\ y_{1222} \\ y_{2211} \\ y_{2212} \\ y_{2221} \\ y_{2222} \end{pmatrix}, \qquad \pmb{X} = \begin{pmatrix} 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \end{pmatrix}, \qquad \beta = \begin{pmatrix} \mu \\ \alpha_1 \\ \alpha_2 \\ \beta_1 \\ \beta_2 \\ \alpha\beta_{11} \\ \alpha\beta_{21} \\ \alpha\beta_{12} \\ \alpha\beta_{22} \end{pmatrix}, \]

\[ \pmb{Z} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}, \qquad \gamma = \begin{pmatrix} p_{(11)1} \\ p_{(11)2} \\ p_{(21)1} \\ p_{(21)2} \\ p_{(12)1} \\ p_{(12)2} \\ p_{(22)1} \\ p_{(22)2} \end{pmatrix}, \qquad \epsilon = \begin{pmatrix} \epsilon_{1111} \\ \epsilon_{1112} \\ \epsilon_{1121} \\ \epsilon_{1122} \\ \epsilon_{2111} \\ \epsilon_{2112} \\ \epsilon_{2121} \\ \epsilon_{2122} \\ \epsilon_{1211} \\ \epsilon_{1212} \\ \epsilon_{1221} \\ \epsilon_{1222} \\ \epsilon_{2211} \\ \epsilon_{2212} \\ \epsilon_{2221} \\ \epsilon_{2222} \end{pmatrix}\cdot \]

\(\widehat{\beta}\) representa o vetor de efeitos fixos, enquanto \(\gamma\) e \(\epsilon\) representa os vetores de efeitos aleatórios. A suposição de efeitos aleatórios independentes e normalmente distribuídos pode ser expressa por \[ \pmb{γ}\sim MVN\big(0, \sigma_p^2\pmb{I}_{abr}\big) \] e \[ \epsilon\sim MVN\big(0, \sigma^2\pmb{I}_{abrs}\big); \] portanto, \(\pmb{y}\sim MVN\big(\pmb{X}\beta,\pmb{V}\big)\), onde \(\pmb{V} = \pmb{Z}\big(\sigma_p^2\pmb{I}_{abr}\big)\pmb{Z}^\top + \sigma^2\pmb{I}_{abrs}\).

O estimador de mínimos quadrados de \(\beta\) seria \(\widehat{\beta} = \big(\pmb{X}^\top \pmb{V}^{-}\pmb{X}\big)^{-} \pmb{X}^\top \pmb{V}^{-}\pmb{y}\), onde \(−\) refere-se a um inverso generalizado; entretanto, \(\sigma_p^2\), \(\sigma^2\) são desconhecidos e, portanto, \(\pmb{V}\) é desconhecido.

A função R aov resolve o problema tratando \(\beta\) e \(\gamma\) como fixos. As somas dos quadrados da ANOVA para cada termo do modelo são calculadas como aquelas mostradas na Tabela 3.2, com um termo adicional chamado Resíduals, que representa as subamostras.

Os comandos para abrir o conjunto de dados com os dados da Tabela 5.14 e criar a tabela ANOVA são mostrados abaixo.

library(daewr)
mod4 = aov( residue ~ form + tech + form:tech + plot:form:tech, data = pesticide)
summary(mod4)
##                Df Sum Sq Mean Sq F value  Pr(>F)    
## form            1 0.0000  0.0000    0.04   0.846    
## tech            1 0.0323  0.0323   72.43 2.8e-05 ***
## form:tech       1 0.0022  0.0022    4.90   0.058 .  
## form:tech:plot  4 0.0023  0.0006    1.31   0.343    
## Residuals       8 0.0036  0.0004                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Os valores da estatística \(F\) e os \(p\)-valores estão incorretos nesta tabela porque usam o termo subamostra ou resíduos como denominador dos testes \(F\). Os quadrados médios esperados para este projeto balanceado podem ser encontrados usando o método tabular de Bennett and Franklin (1954) e são mostrados na Tabela 5.15.

\[ \begin{array}{ccc}\hline \mbox{Fonte} & \mbox{df} & \mbox{Quadrados médios esperados} \\\hline \mbox{A} & a-1 & \sigma^2+r\, \sigma_P^2 + srb\, \tau_A^2 \\ \mbox{B} & b-1 & \sigma^2+r\, \sigma_P^2 + sr\, \tau_{AB}^2 \\ \mbox{Caixa} & (r-1)ab & \sigma^2+r\, \sigma_P^2 \\ \mbox{Sub-amostra} & (s-1)rab & \sigma^2 \\\hline \end{array} \] Tabela 5.15: Quadrados médios esperados para experimento de dois fatores com subamostras.

Nesta tabela \(\sigma^2\) e \(\sigma_P^2\) são as variâncias dos efeitos aleatórios \(\epsilon_{ijkl}\) e \(p_{(ij)k}\) no modelo 5.13. Os termos \(\tau^2\) representam funções quadráticas dos efeitos fixos \(\alpha_i\), \(\beta_j\) e \(\alpha\beta_{ij}\).

Como todos os quadrados médios esperados para formulação, aplicação e sua interação contêm \(\sigma^2 + 2\sigma_P^2\) além de uma forma quadrática envolvendo os efeitos fixos, o quadrado médio correto a ser usado como denominador da razão \(F\) para testar esses efeitos fixos é o termo da caixa ou form:tech:plot cuja esperança é \(\sigma^2 + 2\sigma_P^2\).

O valor esperado do erro ou quadrado médio dos resíduos é \(\sigma^2\) e é muito pequeno para ser usado como denominador na razão \(F\) para testar os efeitos fixos. No entanto, é isso que a tabela ANOVA da função aov usa.

Os valores \(F\) e \(p\)-valores corretos para a ANOVA são: \[ \begin{array}{l} \mbox{form:} \;F_{1,4} = 0.00002/0.00059 = 0.03, \; p\mbox{-valor}=0.8709, \\ \mbox{tech:} \; F_{1,4} = 0.03231/0.00059 = 54.76, \; p\mbox{-valor}=0.00018, \\ \mbox{form:tech} \; F_{1,4} = 0.00219/0.00059 = 3.71, \; p\mbox{-valor}=0.1264, \end{array} \] onde os \(p\)-valores podem ser obtidos usando as funções de probabilidade no R, ou seja, 1-pf(0.03,1,4). Agora pode-se ver que o único termo significativo é a técnica de aplicação.

Em muitos casos em que as unidades observacionais são diferentes das unidades experimentais, as unidades observacionais não serão independentes. Neste exemplo, a aplicação comum do pesticida em cada parcela pode induzir uma correlação, \[ \mbox{E}(\epsilon_{ijkl}\times \epsilon_{ijkl'}) = \rho, \] entre subamostras da mesma parcela.

Mesmo que a suposição de independência seja violada, Casella (2008) mostra o teste \(F\) sobre efeitos fixos usando a média quadrática da parcela, pois o denominador ainda é válido. Neste exemplo, onde a técnica de aplicação foi considerada significativa, devemos observar as médias marginais de aplicação, uma vez que há replicação igual, para determinar qual técnica de aplicação foi a melhor.

As médias marginais podem ser expressas como \(\overline{y}_{\cdot 1\cdot\cdot}\) e \(\overline{y}_{\cdot 2\cdot\cdot}\). O valor esperado é \[ \mbox{E}\big(\overline{y}_{\cdot j\cdot\cdot}\big) = \mu+\beta_j+\overline{\alpha\beta}_{\cdot j\cdot}\cdot \] A variância \[ \mbox{Var}\big(\overline{y}_{\cdot j\cdot\cdot}\big) = \dfrac{\sigma_P^2}{ar} +\dfrac{\sigma^2}{ars}, \] e a variância da diferença em duas médias marginais seria \(2\mbox{Var}\big(\overline{y}_{\cdot j\cdot\cdot}\big)\).

No entanto, o erro padrão da diferença nas médias relatadas pela função estimable no pacote R gmodels é \(\sqrt{2\sigma^2/ars} = 0.01056\), conforme mostrado abaixo.

c1 = c(-0.5, 0.5)
mod4 = aov( residue ~ form + tech + form:tech + plot:form:tech, 
            contrasts = list( form = c1, tech = c1, plot = c1 ), data = pesticide)
efeitos = ("application effect" = c(0,0,1,0,0,0,0,0))
library(gmodels)
estimable(mod4,efeitos)
##                   Estimate Std. Error t value DF Pr(>|t|)
## (0 0 1 0 0 0 0 0)   0.0899     0.0106    8.51  8 2.79e-05


A função R aov não estima o componente de variância para \(\sigma_P^2\), pois a parcela (plot) é tratada como efeito fixo nos cálculos. Um problema semelhante ocorre com todos os erros padrão de funções estimáveis de efeitos fixos no modelo quando a função estimable em gmodels está operando em um objeto criado por aov.

Quando existem efeitos fixos e aleatórios no modelo, devido à forma como o experimento foi conduzido, chamamos de modelo misto e a função lmer no pacote R lme4 é a melhor opção para análise. Por padrão, strong>lmer usa o método REML para estimar os componentes de variância para os efeitos aleatórios no modelo e, em seguida, estima os efeitos fixos usando a fórmula \[ \widehat{\beta} = \big(\pmb{X}^\top \widehat{\pmb{V}}^{-1} \pmb{X} \big)^{-1} \pmb{X}^\top \widehat{\pmb{V}}^{-1}\pmb{y}\cdot \]

Os erros padrão corretos para funções estimáveis de efeitos fixos são produzidos por lmer. Os comandos para executar o lmer com os dados da Tabela 5.14 são mostrados abaixo.

library(lme4)
c1 = c( -0.5, 0.5 )
mod5 = lmer( residue ~ 1 + form + tech + form:tech + (1|plot:form:tech), 
             contrasts = list( form = c1, tech = c1), data = pesticide)
summary(mod5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: residue ~ 1 + form + tech + form:tech + (1 | plot:form:tech)
##    Data: pesticide
## 
## REML criterion at convergence: -51.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5362 -0.6718  0.0541  0.5771  1.7019 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  plot:form:tech (Intercept) 6.99e-05 0.00836 
##  Residual                   4.46e-04 0.02112 
## Number of obs: 16, groups:  plot:form:tech, 8
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.31656    0.00605   52.31
## form1       -0.00212    0.01210   -0.18
## tech1        0.08987    0.01210    7.43
## form1:tech1 -0.04675    0.02421   -1.93
## 
## Correlation of Fixed Effects:
##             (Intr) form1 tech1
## form1       0.000             
## tech1       0.000  0.000      
## form1:tech1 0.000  0.000 0.000


Nestes resultados são encontradas as estimativas dos componentes de variância \(\widehat{\sigma}_P^2 = 0.00006994\) e \(\widehat{\sigma}^2 = 0.000446\) e pode-se observar que o erro padrão correto para a diferença das médias das duas técnicas de aplicação é mostrado como, \(2(\widehat{\sigma}_P^2/ar + \widehat{\sigma}^2/ars) = 0.012103\).

Para projetos com mais de dois níveis para os fatores fixos, a função estimable no pacote gmodels, bem como a função lsmeans no pacote lsmeans, descrita nos Capítulos 3 e 4, produzem os erros padrão corretos das médias ao operar em um objeto criado por lmer em vez de aov. O pacote lsmeans também pode calcular comparações de pares ajustadas de Tukey das médias usando os erros padrão corretos, conforme mostrado no código abaixo.

library(lsmeans)
lsmeans(mod5, pairwise ~ tech, adjust = c("tukey"))
## $lsmeans
##  tech lsmean      SE df lower.CL upper.CL
##  1     0.272 0.00856  4    0.248    0.295
##  2     0.361 0.00856  4    0.338    0.385
## 
## Results are averaged over the levels of: form 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast      estimate     SE df t.ratio p.value
##  tech1 - tech2  -0.0899 0.0121  4  -7.430  0.0018
## 
## Results are averaged over the levels of: form 
## Degrees-of-freedom method: kenward-roger


lmer também produz os testes \(F\) corretos para os fatores fixos, conforme mostrado abaixo.

anova(mod5)
## Analysis of Variance Table
##           npar  Sum Sq Mean Sq F value
## form         1 0.00001 0.00001    0.03
## tech         1 0.02460 0.02460   55.14
## form:tech    1 0.00166 0.00166    3.73


Outro delineamento cujo modelo contém efeitos fixos e aleatórios são os de blocos aleatórios discutidos no último capítulo.

Por exemplo, considere o exemplo da Golf Magazine discutido na Seção 4.7. Nesse experimento; o fator de tratamento, altura do tee, pode ser considerado um efeito fixo, uma vez que o objetivo do experimento era determinar se existe alguma diferença na distância percorrida causada por diferentes alturas do tee.

No entanto, não houve interesse em comparar os jogadores de golfe específicos utilizados na experiência, uma vez que apenas denotavam uma amostra de jogadores de golfe representativos. Portanto, o fator de bloqueio, ou jogador de golfe, pode ser considerado um efeito aleatório. A interação entre um efeito aleatório e fixo também é definida como um efeito aleatório. Portanto, a interação entre o jogador de golfe e a altura do tee também deve ser considerada um efeito aleatório, e o modelo para o experimento de golfe pode ser escrito como: \[\begin{equation*} \tag{5.15} y_{ijk} = \mu + b_i + \tau_j + b\, \tau_{ij} + \epsilon_{ijk}, \end{equation*}\] onde \(b_i\) representa o efeito aleatório do jogador de golfe, \(\tau_j\) representa o efeito da altura fixa do tee, \(b\, \tau_{ij}\) representa o efeito de interação aleatória e \(\epsilon_{ijk}\) representa o efeito de repetição aleatória de golpes, \(i = 1,2,3\), \(j = 1,\cdots,9\) e \(k = 1,\cdots,5\).

Assumindo que os efeitos aleatórios do jogador de golfe, da interação e do golpe repetido são independentes e normalmente distribuídos com variâncias \(\sigma_b^2\), \(\sigma_{b\tau}^2\) e \(\sigma^2\), respectivamente, o valor esperado e a variância de \(y_{ijk}\) são dados por: \[ \mbox{E}(y_{ijk}) = \mu+\tau_i,\\ \mbox{Var}(y_{ijk}) = \sigma_b^2+\sigma_{b\tau}^2+ \sigma^2 \cdot \]

Os valores esperados para os quadrados médios da ANOVA são apresentados na Tabela 5.16, onde \[ Q(\tau) = \dfrac{8}{3-1}\sum_{i=1}^3 \big( \tau_i-\overline{\tau}_{\cdot}\big)^2\cdot \]

\[ \begin{array}{lcl}\\\hline \mbox{Fonte} & \mbox{df} & \mbox{Quadros médios esperados} \\\hline \mbox{Jogados de golfe} \; b_i & (9-1) =8 & \sigma^2+ 5\sigma_{b\tau}^2 + 3\sigma_b^2 \\ \mbox{Altura do tee} \; \tau_j & (3-1)=2 & \sigma^2+5\sigma_{b\tau}^2+Q(\tau) \\ \mbox{Jogador de golfe}\times \mbox{Altura do tee} \; b\, \tau_{ij} & 8\times 2 =16 & \sigma^2+5\sigma_{b\tau}^2 \\ \mbox{Repetições} \; \epsilon_{ijk} = \mbox{Erro} & (5-1)\times 16 = 64 & \sigma^2\\\hline \end{array} \] Tabela 5.16: Quadrados médios esperados para projeto de bloco randomizado com réplicas dentro de um bloco.

Esses quadrados médios esperados mostram que para testar a hipótese de não haver diferença no efeito da altura do tee, ou seja, \(H_0 ∶ \tau_1 = \tau_2 = \tau_3\), o denominador correto para a razão \(F\) seria a interação \[ \mbox{Jogador de golfe}\times \mbox{Altura do tee}, \] conforme usado na Seção 4.7, e não o quadrado médio do erro usado por padrão nas tabelas ANOVA produzidas pela função aov.

Além disso, a variância correta de uma diferença nas médias marginais ou de mínimos quadrados para diferentes alturas do tee \[ \mbox{Var}\big( \overline{y}_{i\cdot\cdot} -\overline{y}_{i'\cdot\cdot}\big)= 2\left(\dfrac{\sigma_{b\tau}^2}{45}+\dfrac{\sigma^2}{45} \right)\cdot \]

Por padrão, a função R aov não estima \(\sigma_{b\tau}^2\). Portanto, no Capítulo 4, a opção Error(id/teehgt) foi usada na chamada da função aov para obter o teste \(F\) correto para a altura do tee. Outra maneira de obter o teste \(F\) correto para a altura do tee é usar a função lmer no pacote lme4.

Além de produzir a estatística \(F\) correta, esta função calculará o erro padrão correto para as diferenças nas médias da altura do tee.


5.9 Métodos gráficos para verificar as suposições do modelo


A análise gráfica de dados de experimentos de amostragem para estudar variâncias é útil para verificar suposições dos modelos estatísticos e identificar observações atípicas que podem ter uma forte influência nas estimativas dos componentes de variância resultantes. Snee (1983) explica que outra vantagem importante da análise gráfica é que ela obriga o analista a se familiarizar mais com os dados e a pensar criticamente sobre os mecanismos para explicar a variação observada.

Ele sugere ferramentas gráficas simples, como gráficos semi-normais e gráficos gama na análise de projetos de amostragem aninhados. Gráficos normais também podem ser usados para verificar a suposição de normalidade para efeitos aleatórios.


5.9.1 Gráficos simples


Considere novamente os dados da Tabela 5.2. O modelo (5.1) é apropriado para os dados deste experimento, onde os efeitos amostrais aleatórios \(t_i\) e os erros experimentais aleatórios \(\epsilon_{ij}\), são assumidos como normalmente distribuídos.

A suposição de variância constante dos \(\epsilon_{ij}\)’s entre as diversas amostras também está implícita em \(\epsilon_{ij}\sim N(0,\sigma^2)\). Uma maneira simples de verificar essas suposições é fazer um boxplot simples dos dados, como o mostrado na Figura 2.3. Lá, variâncias não constantes ou observações discrepantes (outlier) podem ser rapidamente detectadas. A normalidade do \(\epsilon_{ij}\) também pode ser verificada com um gráfico de probabilidade normal como o mostrado na Figura 2.3.

Para o experimento de amostragem de dois fatores Gage R&R apresentado na Seção 5.5, o uso do método dos momentos resultou em uma estimativa negativa do componente de variância para \(\sigma_b^2\), a variância devida aos operadores. Essa estimativa foi negativa porque o quadrado médio da interação da peça por operador foi maior que o quadrado médio do operador.

Gráficos de interação simples como as Figuras 3.4 a 3.6 podem facilitar a explicação desse resultado. A Figura 5.4 mostra o gráfico de interação peça por operador. Uma interação é caracterizada pelo fato de que o traço da resposta média entre os níveis de um fator graficado separadamente para cada nível do outro fator não será paralelo.

Na Figura 5.4, pode-se observar que os segmentos de reta estão próximos de serem paralelos, exceto os segmentos que unem as partes 5 a 7 e 9 a 10. Isso se deve ao fato do operador 1 ter tido uma medição média muito inferior na peça número 6, do que os outros dois operadores, e o operador 3 teve uma medição muito mais elevada na peça número 10 do que os outros dois operadores.

library(daewr)
data(gagerr)
par(mar=c(4,4,1,1))
interaction.plot(gagerr$part, gagerr$oper, gagerr$y, type="l", 
  main = "Interaction Plot for part by Operator",  xlab="Part Number", ylab="Average measurement")
grid()

Figura 5.4: Gráfico de interação para peça por operador.

Quando a média quadrada da interação se deve em grande parte a uma ou duas observações, como é neste caso, pode ser sensato questionar os resultados. A Tabela 5.17 mostra as estimativas dos componentes de variância a partir dos dados do estudo Gage R&R após a eliminação das partes 6 e 10. Com essas partes eliminadas o método dos momentos não produz uma estimativa negativa do componente de variância e os resultados obtidos pelo método dos momentos e REML são bastante semelhantes.

\[ \begin{array}{lcc} \hline & \mbox{Método dos momentos} & \mbox{REML} \\ \mbox{Componente} & \mbox{Estimativas} & \mbox{Estimativas} \\\hline \mbox{part} \, \sigma_a^2 & 0.0319100 & 0.0280800 \\ \mbox{oper} \, \sigma_b^2 & 0.0008601 & 0.0008089 \\ \mbox{part}\times \mbox{oper} \, \sigma_{ab}^2 & 0.0020045 & 0.0020082 \\ \mbox{Erro} \, \sigma^2 & 0.0004062 & 0.0004063 \\\hline \end{array} \] Tabela 5.17: Comparação do método dos momentos e estimativas REML no estudo de medição R&R após eliminação das partes 6 e 10.

Numa situação como esta, se as peças ainda estiverem disponíveis, pode ser sensato pedir aos operadores que meçam novamente as peças 6 e 10 para ver se houve erros.


5.9.2 Gráficos gama e semi-normais para verificar a suposição de variância constante


Wilk et al. (1962) propuseram gráficos de probabilidade gama como ferramentas eficazes para investigar a homogeneidade de grupos de variâncias. A amostra de variâncias, com \(\nu\) graus de liberdade, segue a distribuição gama com parâmetro de forma \(\nu/2\) e, portanto, os pontos em um gráfico de probabilidade gama de variâncias amostrais que estimam a mesma constante \(\sigma^2\) devem aparecer como uma linha reta.

Muitos experimentos de amostragem aninhada e experimentos de amostragem aninhada escalonada envolvem apenas amostras de tamanho dois. Neste caso, gráficos semi-normais de desvios padrão podem ser usados da mesma maneira. Para ilustrar isto, considere novamente os dados da Tabela 5.2.

O código R abaixo calcula as variâncias dentro de cada amostra usando a função tapply, armazena-as no conjunto de dados s2 e, em seguida, usa a função qgamma para calcular os quantis gama. O gráfico das pontuações gama versus as variâncias amostrais é mostrado na Figura 5.5.

library(daewr)
data(Naph)
s2 = tapply( Naph$yield, Naph$sample, var )
os2 = sort(s2)
r = c( 1:length(s2) )
gscore =- qgamma( (r - 0.5 ) / length (s2), 2)
par(mar=c(4,4,1,1), pch=19)
plot(gscore, os2, main = "Gamma plot of within sample variances", 
     xlab = "Gamma score", ylab = "Sample Variance")
grid()

Figura 5.5: Gráfico gama de variâncias dentro da amostra a partir dos dados da Tabela 5.2.

Se as variâncias forem homogêneas, os pontos no gráfico de probabilidade gama ficarão aproximadamente ao longo de uma linha reta. Pontos individuais acima da linha à extrema direita indicariam variações maiores que o normal.

Na Figura 5.5, com apenas 6 pontos, parece que os pontos ficam bastante próximos de uma linha reta, indicando variações constantes. Para obter um exemplo onde há mais variações para representar graficamente, considere o estudo de polimerização mostrado na Tabela 5.13.

O estimador do métodos de momentos para o componente de variância caixa dentro do lote para este estudo de amostragem escalonada aninhada foi negativo. Tal como acontece com o estudo de medição R&R, um ou dois valores atípicos podem ser a causa deste resultado, e uma análise gráfica dos dados pode revelar isso.

Os quadrados médios na ANOVA apresentada na Seção 5.7 são estimativas agrupadas das variâncias para cada fonte no projeto. Se as quatro observações para o número de lote \(i\) forem denotadas por \(Y_{1i}\), \(Y_{2i}\), \(Y{3i}\) e \(Y_{4i}\) conforme mostrado abaixo, \[ \begin{array}{ccccc}\hline \mbox{Caixa} & 1 & 1 & 1 & 2 \\ \mbox{Preparação} & 1 & 1 & 2 & 1 \\\hline \mbox{Lote} & \mbox{Teste 1} & \mbox{Teste 2} & \mbox{Teste 1} & \mbox{Teste 1} \\\hline i & Y_{1i} & Y_{2i} & Y_{3i} & Y_{4i} \\ \end{array} \]

Snee (1983) mostra que as variâncias a serem agrupadas de cada fonte para criar os quadrados médios da ANOVA são dadas por: \[ \begin{array}{ll} \mbox{Fonte} & \mbox{Variância} \quad \sigma_i^2 \\\hline \mbox{Erro ou teste (prep)} & \frac{1}{2}\big(Y_{2i}-Y_{1i}\big)^2 \\ \mbox{prep (caixa)} & \frac{2}{3}\Big(Y_{3i}-\frac{1}{2}\big(Y_{1i}+Y_{2i}\big)\Big)^2 \\ \mbox{caixa} & \frac{3}{4}\Big(Y_{4i}-\frac{1}{3}\big(Y_{1i}+Y_{2i}+Y_{3i}\big)\Big)^2 \\ \end{array} \]

O código R para calcular os desvios padrão dentro de cada fonte é mostrado abaixo.

library(daewr)
data(polymer)
y = array( polymer$strength, c(4,30) )
sd1 = sqrt( (y[2,] - y[1,])**2 / 2)
sd2 = sqrt( (2/3) * ( y[3,] - (y[1,] + y[2,]) / 2)**2 )
sd3 = sqrt( (3/4) * (y[4,] - (y[1,] + y[2,] + y[3,] ) / 3 )**2)


O método do estimador de momentos do componente de variância para caixa dentro do lote foi negativo porque o quadrado médio para preparação dentro da caixa e lote foi maior que o quadrado médio para caixa. Isto sugeriria uma investigação das variações que são agrupadas para formar o quadrado médio para preparação dentro da caixa.

O código para fazer um gráfico meio normal dos desvios padrão para preparação dentro da caixa sd2 é mostrado abaixo.

osd2 = sort(sd2)
r = c( 1: length(sd2))
zscore = qnorm( ( ( r - .5 ) / length(sd2) +1 )/ 2)
par(mar=c(4,4,1,1), pch=19)
plot( zscore, osd2, main = "Half-normal plot of prep(box) standard deviations", 
      xlab = "Half Normal Score", ylab = "std. due to prep within box")
grid()

Figura 5.6: Gráfico meio normal de desvios padrão de preparação (caixa).

No gráfico mostrado na Figura 5.6 acima, vemos uma linha relativamente reta de pontos que se estende da parte inferior esquerda do gráfico, exceto que há um desvio padrão que se destaca acima da linha no lado direito.

A Tabela 5.18 mostra os desvios padrão, que foram calculados com o código R acima, dentro de cada fonte para cada lote. Nesta tabela o desvio padrão para preparação dentro da caixa é rotulado como \(s_2\), e pode-se observar que há um valor grande para \(s_2\) no lote 19, o que se deve principalmente ao alto resultado para \(Y_3 =16.79\), o primeiro teste da segundo preparação para a caixa 1.

\[ \begin{array}{cccccccc}\hline \mbox{Lote} & Y_1 & Y_2 & Y_3 & Y_4 & s_1 & s_2 & s_3 \\\hline 1 & 9.76 & 9.24 & 11.91 & 9.02 & 0.368 & 1.968 & 1.111 \\ 2 & 10.65 & 7.77 & 10.00 & 13.69 & 2.036 & 0.645 & 3.652 \\ 3 & 6.50 & 6.26 & 8.02 & 7.95 & 0.170 & 1.339 & 0.886 \\ 4 & 8.08 & 5.28 & 9.15 & 7.46 & 1.980 & 2.017 & 0.038 \\ 5 & 7.84 & 5.91 & 7.43 & 6.11 & 1.365 & 0.453 & 0.823 \\ 6 & 9.00 & 8.38 & 7.01 & 8.58 & 0.438 & 1.372 & 0.390 \\ 7 & 12.81 & 13.58 & 11.13 & 10.00 & 0.544 & 1.686 & 2.171 \\ 8 & 10.62 & 11.71 & 14.07 & 14.56 & 0.771 & 2.372 & 2.102 \\ 9 & 4.88 & 4.96 & 4.08 & 4.76 & 0.057 & 0.686 & 0.104 \\ 10 & 9.38 & 8.02 & 6.73 & 6.99 & 0.962 & 1.608 & 0.912 \\ 11 & 5.91 & 5.79 & 6.59 & 6.55 & 0.085 & 0.604 & 0.393 \\ 12 & 7.19 & 7.22 & 5.77 & 8.33 & 0.021 & 1.172 & 1.389 \\ 13 & 7.93 & 6.48 & 8.12 & 7.43 & 1.025 & 0.747 & 0.069 \\ 14 & 3.70 & 2.86 & 3.95 & 5.92 & 0.594 & 0.547 & 2.093 \\ 15 & 4.64 & 5.70 & 5.96 & 5.88 & 0.750 & 0.645 & 0.387 \\ 16 & 5.94 & 6.28 & 4.18 & 5.24 & 0.240 & 1.576 & 0.196 \\ 17 & 9.50 & 8.00 & 11.25 & 11.14 & 1.061 & 2.041 & 1.348 \\ 18 & 10.93 & 12.16 & 9.51 & 12.71 & 0.870 & 1.662 & 1.596 \\ 19 & 11.95 & 10.58 & 16.79 & 13.08 & 0.969 & 4.511 & 0.023 \\ 10 & 4.34 & 5.45 & 7.51 & 5.21 & 0.785 & 2.135 & 0.482 \\ 20 & 4.34 & 5.45 & 7.51 & 5.21 & 0.785 & 2.135 & 0.482 \\ 21 & 7.60 & 6.72 & 6.51 & 6.35 & 0.622 & 0.531 & 0.514 \\ 22 & 5.12 & 5.85 & 6.31 & 8.74 & 0.516 & 0.674 & 2.581 \\ 23 & 5.28 & 5.73 & 4.53 & 5.07 & 0.318 & 0.796 & 0.095 \\ 24 & 5.44 & 5.38 & 4.35 & 7.04 & 0.042 & 0.865 & 1.718 \\ 25 & 3.50 & 3.88 & 2.57 & 3.76 & 0.269 & 0.914 & 0.384 \\ 26 & 4.80 & 4.46 & 3.48 & 3.18 & 0.240 & 0.939 & 0.924 \\ 27 & 5.35 & 6.39 & 4.38 & 5.50 & 0.735 & 1.217 & 0.110 \\ 28 & 3.09 & 3.19 & 3.79 & 2.59 & 0.071 & 0.531 & 0.664 \\ 29 & 5.30 & 4.72 & 4.39 & 6.13 & 0.410 & 0.506 & 1.149 \\ 30 & 7.09 & 7.82 & 5.96 & 7.14 & 0.516 & 1.221 & 0.159 \\ \end{array} \] Tabela 5.18: Dados brutos para cada lote e desvios padrão calculados.

Se o lote número 19 for removido dos dados, a estimativa do método dos momentos do componente de variância para a caixa dentro do lote não é mais negativa e o método dos momentos e os estimadores REML são bastante consistentes, como pode ser visto nos resultados mostrados na Tabela 5.19.

Isto pode levar a questionar o resultado do teste 1 da preparação 2 da caixa do lote 19. Se o material ainda estivesse disponível, seria necessário repetir o teste.

\[ \begin{array}{lll} \hline & \mbox{Método dos momentos} & \mbox{REML} \\ \mbox{Componente} & \mbox{Esttimativa} & \mbox{Estimativa} \\\hline \mbox{Lote} \; \sigma_a^2 & 5.81864 & 6.09918 \\ \mbox{Caixa (Lote)} \; \sigma_b^2 & 0.13116 & 0.04279 \\ \mbox{Preparação (Caixa)} \; \sigma_c^2 & 0.76517 & 0.79604 \\ \mbox{Erro} \; \sigma^2 & 0.63794 & 0.64364 \\\hline \end{array} \] Tabela 5.19: Comparação do método dos momentos e estimativas REML para estudo de polimerização após remoção do lote 19.


5.9.3 Gráficos de probabilidade dos EBLUP’s de efeitos aleatórios


Na Seção 2.4, para verificar a suposição de que os erros experimentais seguem uma distribuição normal, foi feito um gráfico de probabilidade normal dos resíduos.

Na notação matricial, o vetor de resíduos pode ser escrito como: \[\begin{equation*} \tag{5.16} \widehat{\pmb{\epsilon}}=\pmb{y}-\pmb{X}\widehat{\pmb{\beta}}\cdot \end{equation*}\]

Quando há termos aleatórios no modelo além do termo de erro experimental, podemos verificar a normalidade desses efeitos aleatórios fazendo gráficos de probabilidade normal dos efeitos aleatórios estimados.

Uma linha reta de pontos nesse gráfico justificaria a suposição de normalidade para o efeito aleatório em questão. Esta técnica gráfica será útil na detecção de observações atípicas ou desvios da normalidade quando houver pelo menos 12 a 15 pontos no gráfico.

No modelo \[\begin{equation*} \tag{5.17} \pmb{y} = \pmb{X}{\pmb{\beta}}+\pmb{Z}\pmb{\gamma}+{\pmb{\epsilon}}, \end{equation*}\] os melhores preditores lineares não viciados empíricos (EBLUPs) dos efeitos aleatórios \(\gamma\) são dados pela equação \[\begin{equation*} \tag{5.18} \widehat{\pmb{\gamma}} = \widehat{\pmb{G}}\pmb{Z}^\top \widehat{\pmb{V}}^{-1}(\pmb{y}-\pmb{X}\widehat{\pmb{\beta}}), \end{equation*}\] onde \(\widehat{\pmb{G}}\) é a matriz de covariância de variância estimada de \(\gamma\).

A função lmer no pacote R lme4 calcula os EBLUPs e pode ser usada para modelos com todos os efeitos aleatórios, bem como para modelos com efeitos fixos e aleatórios.

Abaixo está o código R para ajustar o modelo de efeitos aleatórios aninhados para os dados de polimerização usando lmer. Os EBLUPs podem ser recuperados do objeto criado por lmer com o comando ranef().

Os pontos no gráfico normal dos EBLUPs para preparação dentro da caixa mostrada na Figura 5.7 seguem quase uma linha reta. No entanto, existem dois valores discrepantes aparentes na parte superior direita do gráfico. Ao comparar os valores destes dois EBLUPs com uma lista de dados, verifica-se que o maior EBLUP (2.2488) corresponde ao primeiro teste da segunda preparação para a caixa 1 do lote 19.

library(lme4)
modr3 = lmer( strength ~ 1 + (1|lot) + (1|lot:box) + (1|lot:box:prep), data = polymer)
ranef(modr3)$`lot:box:prep`[1:6,]
## [1] -0.305  1.225 -0.539 -1.049 -0.362  1.891
ranef(modr3)$`lot:box`[1:6,]
## [1]  2.12e-08 -1.24e-08 -3.25e-08  4.35e-08 -7.83e-09  8.29e-09
ranef(modr3)$lot[1:6,]
## [1]  2.682  3.373  0.139  0.431 -0.383  0.856
par(mar=c(4,4,1,1), pch=19)
qqnorm( ranef(modr3)$"lot:box:prep"[[1]], 
        main= "prep within box and lot", ylab="EBLUP", xlab = "Normal Score" )
grid()

Figura 5.7: Gráfico normal de efeitos aleatórios estimados para prep(box).

Este é o mesmo ponto identificado pelo gráfico meio normal na última seção e merece uma investigação mais aprofundada. O segundo maior EBLUP (1.8905) está associado ao teste da preparação 1 na caixa 2 do lote 2.

Este valor (13.69) é superior a qualquer outro valor para o lote 2 e novamente pode justificar uma investigação mais aprofundada.


5.10 Exercícios