MODELAGEM DA SUPERDISPERSÃO EM DADOS
POR UM MODELO LINEAR GENERALIZADO MISTO
José Airton Rodrigues NUNES1
Augusto Ramalho de MORAIS2
Júlio Sílvio de Sousa BUENO FILHO2
RESUMO: A modelagem de dados binomiais, em muitos casos, constitui-se numa tarefa difícil
uma vez que estes, não raramente, apresentam variação extra-binomial. O objetivo deste trabalho
constituiu-se na aplicação de modelos lineares generalizados mistos (MLGM’s), em dois
conjuntos de dados da literatura, com o intuito de acomodar de forma satisfatória a
superdispersão presente. As estimativas dos parâmetros foram derivadas a partir da função de
quase-logverossimilhança penalizada conjunta utilizando função de ligação logística por meio
dos seguintes algoritmos: macro glimmix do programa SAS 8.10, glmmPQL (Biblioteca MASS)
do programa R 1.6.1 e algoritmo MLGM implementado no R 1.6.1. Primeiramente os conjuntos
de dados foram submetidos ao ajuste pelo modelo binomial padrão, obtendo-se fortes evidências
de superdispersão em ambos. Os resultados do ajuste pelos MLG’s com incorporação de efeito
aleatório de parcela mostraram que os modelos sugeridos conseguiram explicar satisfatoriamente
a variabilidade presente nos dados. Observou-se ainda que as proporções ajustadas pelo
algoritmo implementado no R 1.6.1 foram iguais às obtidas pelo SAS 8.10 para os conjuntos de
dados. Em virtude dos resultados apresentados pode-se aferir que o algoritmo implementado
constitui-se um procedimento bastante confiável para ajuste de MLGM’s para dados binomiais.
PALAVRAS-CHAVE: Algoritmo; variação extra-binomial; quase-verossimilhança penalizada;
desvio; ensaio de germinação.
1 Introdução
Na pesquisa agronômica, o pesquisador se depara, não raramente, com situações nas
quais os dados obtidos apresentam distribuição binomial. Nestes casos, as pressuposições
básicas requeridas para aplicação da metodologia da análise de variância associada ao
teste F, técnica normalmente utilizada, são violadas. Negligenciando estas restrições, o
pesquisador incorrerá em elevadas taxas de erro para os testes de hipóteses realizados e
em inferências pouco confiáveis.
Para tornar válida a utilização deste método estatístico, os pesquisadores têm
adotado a mudança adequada da escala da variável aleatória por meio de transformações
nestes dados.
1
Departamento de Biologia, Universidade Federal de Lavras - UFLA, Caixa postal 37, CEP: 37200-000, Lavras,
MG, Brasil. E-mail: [email protected].
Departamento de Ciências Exatas, Universidade Federal de Lavras - UFLA, Caixa postal 37, CEP: 37200-000,
Lavras, MG, Brasil. E-mail: [email protected] / [email protected].
2
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
55
Com a introdução dos modelos lineares generalizados, os problemas com escalas
foram grandemente reduzidos (McCullagh e Nelder, 1989). Trata-se de uma extensão dos
modelos lineares, desenvolvida por Nelder e Wedderburn (1972), para dados não
normalmente distribuídos. Esta metodologia motiva-se no fato que os efeitos sistemáticos
são linearizados por uma transformação adequada dos valores esperados, permitindo aos
valores ajustados variarem dentro da amplitude real das respostas.
Não obstante o uso deste método estatístico, dados binomiais apresentam, em muitas
ocasiões práticas, uma variância nas respostas superior à variância nominal da distribuição
binomial comportada pelo modelo, denominada de variação extra-binomial ou
superdispersão. Vários autores têm mencionado a importância de considerar a presença
deste fenômeno na modelagem e, com isso, têm sugerido várias formas de lidar com este
problema prático.
A incorporação de efeitos aleatórios no preditor linear, os modelos lineares
generalizados mistos (MLGM´s), têm se mostrado numa técnica de grande utilidade e
aplicabilidade na área das ciências biológicas para acomodação da superdispersão. Esta se
fundamenta numa extensão da teoria dos modelos mistos para dados com distribuições
pertencentes à família exponencial, assumindo-se uma distribuição particular para os
efeitos aleatórios.
O objetivo deste trabalho foi a aplicação de modelos lineares generalizados mistos
com função de ligação logística baseado numa quase-logverossimilhança penalizada
conjunta, em dois ensaios envolvendo dados binomiais superdispersos, com o intuito de
acomodar de forma satisfatória a variabilidade extra presente. Um importante objetivo
adicional foi o de ajustar estes modelos com um algoritmo implementado em ambiente R
1.6.1 e comparar os resultados aos de algoritmos já implementados nos programas
estatísticos SAS 8.10 (macro GLIMMIX ), e R 1.6.1 (procedimento glmmPQL).
2 Material e métodos
2.1 Material
No presente trabalho foram analisados dois conjuntos de dados presentes na
literatura com respostas binomiais com fins metodológicos de aplicação. A descrição
desses conjuntos de dados é feita a seguir.
Exemplo 1: Ensaio de germinação
Crowder (1978), citado por Breslow e Clayton (1993), apresenta dados sobre a
proporção de sementes germinadas observadas (ri) em 21 bandejas provenientes de um
experimento conduzido num delineamento inteiramente casualizado com os tratamentos
repetidos, dispostos num esquema fatorial cruzado 2 x 2, sendo duas espécies de sementes
(Orobanche aegyptiaco 75 (= 0), Orobanche aegyptiaco 73 (= 1)) e dois diferentes meios
para germinação provenientes de extratos de raiz (feijão (= 0), pepino (= 1)). As sementes
de Orobanche foram colocadas em diluições 1/125 dos referidos extratos.
Exemplo 2: Ensaio de sobrevivência
Manly (1978), citado por Hinde e Demétrio (1998), retrata um experimento sobre a
sobrevivência de ovos de truta. O experimento constitui-se de caixas contendo
56
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
quantidades variáveis de ovos, que foram colocadas aleatoriamente em cinco diferentes
locais num fluxo e a quatro diferentes períodos (4, 7, 8 e 11 semanas). Uma caixa de cada
local foi amostrada e os números de ovos sobreviventes foram contados.
3 Métodos
3.1 Modelo linear generalizado binomial
Tendo os dados descritos proporções de sucessos (ri) em amostras de tamanho (mi),
admitiu-se, primeiramente, o modelo binomial para explicar os dados apresentados.
Assim, conforme a teoria de modelos lineares generalizados considerando efeitos fixos,
descrita extensivamente em McCullagh e Nelder (1989), Dobson (1990) e Demétrio
(2001), optou-se pelo uso da função de ligação logística (canônica) como uma
transformação do valor esperado que se deseja modelar como uma combinação linear nos
parâmetros, conforme estrutura dos experimentos, devido à interpretação mais simples
como o logaritmo da razão de chances. Assim, o MLG ajustado é
E (R ) = π
π
η = log
= Xβ
(1)
1−π
em que:
R é o vetor de proporções observadas, de dimensões n x 1, tal que se admite:
Ri ~
Bin( mi , π i )
, com i = 1,...,n;
mi
η é o vetor dos preditores lineares, de dimensões n x 1;
X é a matriz do delineamento, de dimensões n x p;
β é o vetor de p parâmetros desconhecidos do preditor linear do modelo, de dimensões
p x 1.
Para efetuar o ajuste dos modelos generalizados propostos para os conjuntos de
dados descritos foram utilizados o procedimento GENMOD do programa SAS 8.10,
comando GLM do programa R 1.6.1 (Ihaka e Gentleman, 1996) e um algoritmo MLG
implementado no programa R 1.6.1. O algoritmo implementado no R utiliza o
procedimento de quadrados mínimos reponderados iterativamente para obtenção das
estimativas dos parâmetros com uso do algoritmo de Escore de Fisher (Demétrio, 2001).
Exemplo 1: Ensaio de germinação
A proporção de sementes germinadas (rijk) em bandejas com mijk sementes foram
modeladas admitindo-se dois possíveis modelos, “M1” com efeitos principais e “M2”
incluindo o efeito da interação:
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
57
M 1 : η ijk = log
M 2 : η ijk = log
π ijk
1 − π ijk
π ijk
1 − π ijk
= α 0 + α 1 x1(ijk ) + α 2 x 2 (ijk )
(2)
= α 0 + α 1 x1(ijk ) + α 2 x 2 (ijk ) + α 12 x1(ijk ) x 2(ijk )
em que:
ηijk = o preditor linear correspondente à observação rijk, k = 1,...,nij;
x1(ijk) = 1 se a variedade de semente i foi O. aegyptiaco 73, 0 se O. aegyptiaco 75;
x2(ijk) = 1 se o extrato de raiz j foi pepino, 0 se feijão;
α = representa os efeitos fixos associados com variedade de semente, extrato de raiz e
interação.
Exemplo 2: Ensaio de sobrevivência
As proporções de ovos sobreviventes (rij), obtidas a partir da contagem em caixas
com mij ovos, foram analisados conforme MLG descrito por
ηij = log
π ij
1 − π ij
= µ + αi + β j
(3)
em que:
ηij = o preditor linear correspondente à observação rij, i = 1,...,5 e j = 1,...,4;
µ = constante associada a todas as observações;
αi = o efeito do local i;
βj = o efeito do período j.
4 Critérios de verificação de ajuste
4.1
“Deviance” binomial residual
Para verificação do ajuste dos modelos descritos utilizou-se a “deviance” como
medida de ajustamento. A “deviance” binomial residual é dada por
D(r , πˆ ) = 2[l (r , r ) − l (πˆ , r )]
D(r , πˆ ) =
n
i −1
2mi ri log
ri
1 − ri
+ (1 − ri ) log
πˆ i
1 − πˆ i
(4)
Por conseguinte, realizou-se o teste da adequação do modelo corrente com p
parâmetros linearmente independentes por meio da estatística
χ (2n− p;α ) ; se a “deviance”
binomial residual (4) fosse inferior a este quantil superior, não existia evidências para
supor que o acréscimo de parâmetros no modelo forneceria ganhos no ajuste; caso
58
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
contrário, o modelo proposto não comportaria a variação presente nos dados, a qual em
muitos casos se deve à variação extra-binomial presente.
4.2 Estimativa do parâmetro de dispersão
A adequação do modelo linear generalizado binomial foi complementarmente
verificada pela estimação do parâmetro de escala (φ); dentre vários estimadores
encontrados na literatura, utilizou-se o estimador de momentos da estatística X2
generalizada de Pearson, dado por
φˆ =
W ( z − ηˆ )
( z − ηˆ ) '
n− p
(5)
sendo z o vetor de variáveis dependentes ajustadas, η
η̂ o vetor de estimativas dos
preditores lineares, n o número de observações e p o número de parâmetros linearmente
independentes.
A estimativa (5) foi utilizada para checagem de superdispersão, indicada por valores
calculados significativamente maiores 1.
4.3 Modelo linear generalizado binomial misto
A aplicação do modelo binomial padrão em dados de proporções superdispersos é
limitada pelo fato de a variância ser automaticamente determinada estando a média
especificada (Hinde e Demétrio, 1998). Assim, evidenciada a presença de variação extrabinomial efetuou-se a modelagem dos dados correspondentes pelo uso da abordagem que
utiliza a incorporação de efeitos aleatórios no preditor linear, ou seja, por um modelo
linear generalizado misto descrito em McCulloch e Searle (2001).
Os modelos lineares generalizados mistos (MLGM) surgem em experimentos
aleatorizados em que se observam dados discretos da mesma forma como surgem os
modelos lineares Gauss-Markov normal (MLGMN) para uma variável contínua.
Assumindo aditividade entre efeitos de unidade experimental (UE) e de tratamentos, o
modelo hierárquico produzido é aproximadamente normal-normal (nos níveis de parcelas
e dados) e a convolução de duas normais gera uma nova distribuição normal. No caso dos
dados discretos, a convolução da aproximação normal para o efeito subjacente de parcelas
com a distribuição dos dados observados (contagens ou proporções) gera um MLGM na
forma de modelo hierárquico (Normal-Poisson ou Normal-Binomial, respectivamente).3
Seja r1, r2,..., rn um conjunto de realizações condicionalmente independentes de
Bin(mi , π i )
variáveis aleatórias Ri | ui ~
, com função de densidade discreta dada por
mi
3
GILMOUR, S. G. 2002, Lavras - MG, comunicação pessoal.
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
59
Ri | ui ~ iid . f Ri|u (ri | ui )
f Ri|u (ri | u i ) = exp mi ri log
i
πi
1−πi
+ log(1 − π i ) + log
mi
(6)
mi ri
O MLGM para análise destas proporções condicionais foi construído utilizando-se
uma função de ligação logística, resultando no seguinte modelo proposto
π = E (R | u )
log
π
1−π
=η = X + u
(7)
em que:
R|u é o vetor de proporções condicionalmente independentes, de dimensões n x 1, tal que
Bin (mi , π i )
se admite Ri | u i ~
, correspondente a (6);
mi
η é o vetor dos preditores lineares, de dimensões n x 1;
X é a matriz do delineamento referente aos efeitos fixos, de dimensões n x p;
Z é a matriz do delineamento referente aos efeitos aleatórios, de dimensões n x q;
β é o vetor de efeitos fixos desconhecidos do preditor linear do modelo, de dimensões
p x 1;
u é o vetor de efeitos aleatórios desconhecidos do preditor linear do modelo, de dimensões
q x 1, assumindo-se u ~ Nq(0, G);
Assim, tendo-se admitido o modelo (7), submeteu-se os dados binomiais a uma
análise pela teoria dos MLGM’s conforme descrito nos trabalhos de Gilmour et al. (1985),
Shall (1991) e Breslow e Clayton (1993), assumindo-se os efeitos aleatórios normalmente
distribuídos.
As estimativas dos vetores de parâmetros β e u foram obtidos a partir da
maximização da função de quase-logverossimilhança penalizada conjunta (QVP) dada por
QVP (π , u , r ) =
i
mi ri log
πi
1
+ log(1 − π i ) − u 'G −1u
1− πi
2
(8)
As diferenciações realizadas resultam num sistema de equações pouco tratável;
porém, com base na variável dependente ajustada dada por
y* = η + ∆(r − π )
(9)
tem-se que esta pode ser aproximada por um modelo linear misto (Henderson et al., 1959)
da forma
y* = Xβ + Ζ + e*
60
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
em que:
e* = ∆(r - π) ; Cov(e*) = W -1 ; Cov(u) = G;
E(Y*) = Xβ ; Cov(Y*) = W -1 + ZGZ’.
Logo, aqueles autores sugerem o uso de algoritmo iterativo similar, em que um
modelo generalizado misto é ajustado para encontrar soluções para estimativas dos
vetores de parâmetros β e u por
βˆ
uˆ
=
X'
WX
X'
WZ
− X'
Wy*
WZ +G −1
Z'
WZ Z '
(10)
Z'
Wy*
em que, para um MLGM para dados binomiais, tem-se que
πi
∂η1
∂
log
= diag
1− π i
∂η1
∂η1
∆ = diag
W = diag
∂π i
∂π i
2
[V ( y i | u i )]−1
= diag
= diag
m i [π i (1 − π i )]2
π i (1 − π i )
1
π i (1 − π i )
= diag{mi π i (1 − π i )} .
Os MLGM’s foram ajustados por algoritmos já implementados nos pacotes
estatísticos SAS 8.10, por meio da macro GLIMMIX associada ao PROC MIXED, e R
1.6.1, pelo comando glmmPQL, e ainda por um algoritmo implementado em ambiente R
1.6.1, sendo as soluções obtidas a partir da maximização da função de quaselogverossimilhança penalizada conjunta (8).
O algoritmo implementado no R 1.6.1 utiliza o algoritmo de Newton-Raphson no
passo do modelo generalizado e o algoritmo EM no passo de maximização para estimar os
componentes de dispersão referente aos efeitos aleatórios.
O algoritmo utilizado pelo SAS 8.10, assim como o algoritmo implementado no R
1.6.1 constituem algoritmos de maximização conjunta descritos em Shall (1991) e
Breslow e Clayton (1993) e obtêm as estimativas dos parâmetros por meio da resolução
do sistema de equações dos MLGM’s (9), enquanto que as estimativas obtidas pelo
comando glmmPQL do programa R 1.6.1, derivadas de (8), utiliza o procedimento de
quadrados mínimos reponderados iterativamente.
Exemplo 1: Ensaio de germinação
As proporções condicionalmente independentes de sementes germinadas (rijk |bk)
em bandejas com mijk sementes foram modeladas admitindo-se dois possíveis modelos,
“M1” com efeitos principais e “M2” incluindo o efeito da interação:
(
)
Μ 2 : logit (Pr = 1 | xijk , bk ) = α 0 + α1 x1(ijk ) + α 2 x 2(ijk ) + α12 x1(ijk ) x 2 (ijk ) + bk
Μ1 : logit Pr = 1 | xijk , bk = α 0 + α1 x1(ijk ) + α 2 x 2(ijk ) + bk
(11)
em que:
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
61
bk = o efeito aleatório associado à parcela k, considerando bk ~ N(0,σ2p).
Exemplo 2: Ensaio de sobrevivência
As proporções condicionalmente independentes de ovos sobreviventes no local i e
período j (rij|bk) em caixas com mij ovos, foram modeladas admitindo-se um MLGM por:
(
)
ηij = log Pr = 1 | α i , β j , bk = µ + α i + β j + bk
(12)
no qual bk é o efeito aleatório da parcela k, considerando bk ~ N(0,σ2p).
4.4 Verificação de ajuste
Verificou-se o ajuste do MLGM’s propostos por meio do cálculo da estimativa do
componente de extra-dispersão ( σ̂ ) conforme descrito por Shall (1991). O teste para a
2
superdispersão presente foi realizado verificando se os valores de σ̂ eram
significativamente maiores que 1, respectivamente (Shall, 1991). O estimador sugerido foi
calculado de forma análoga a (5) e aproximado pelo algoritmo implementado por
2
σˆ 2 =
(z − Xβˆ − Ζuˆ )'W (z − Xβˆ − Zuˆ )
n − r(X )
(13)
em que, r ( X ) é o posto da matriz de efeitos fixos.
5 Resultados e discussão
Exemplo 1: Ensaio de germinação
Primeiramente foram ajustados os modelos de regressão logística ordinários para
estes dados. A Tabela 1 apresenta as estimativas de máxima verossimilhança dos
coeficientes de regressão logística para os modelos ajustados, de (2), para as 21
proporções de sementes germinadas por meio de programas executados nos pacotes
estatísticos SAS 8.10 e R 1.6.1. As estimativas dos parâmetros obtidas foram idênticas
para as três análises do MLG proposto.
Verifica-se que para o modelo de efeitos principais “M1”, de (2), a estimativa de φ,
de (5), foi de 2,1284, excedendo o valor assumido para o modelo (φ = 1) , evidenciando
presença de variação extra não comportada pelo modelo sugerido. O valor da “deviance”
residual, de (4), para o modelo “M1” com 18 graus de liberdade, foi de 39,6859,
tendo (Pr > D = 0,0023) , rejeitando-se a hipótese de que o modelo está bem ajustado e
somando forte evidência de variação extra não modelada ou má especificação do modelo.
A partir da evidência do efeito da interação dos fatores presente, conforme mostrado
na Figura 1, ajustou-se o modelo “M2”, de (2), incluindo o efeito da interação da
variedade de semente e tipo de raiz. Para este modelo, a “deviance” residual, de (4), com
17 graus de liberdade, foi de 33,2778, a qual superou o valor do X 2 (17;0.05) = 27,5871 ,
62
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
rejeitando-se a hipótese de ajustamento do modelo. A estimativa de φ, de (5), para o
modelo sugerido foi igual a 1,8618, excedendo o valor assumido para o modelo (φ=1),
evidenciando-se presença de superdispersão.
Tabela 1 - Estimativas dos parâmetros para os modelos de regressão logística descritos em
(2) por três algoritmos para MLG
ALGORITMOS
Variável
GENMOD – SAS
GLM – R
ALGORITMO – MLG
βˆ ( EP)
βˆ ( EP)
βˆ ( EP)
Modelo com efeitos principais
Intercepto
-0,43 (0,1137)
-0,43 (0,1137)
-0,43 (0,1137)
Variedade
-0,27 (0,1547)
-0,27 (0,1547)
-0,27 (0,1547)
Extrato raiz
1,07 (0,1442)
1,07 (0,1442)
1,07 (0,1442)
φˆ
2,1284
2,1284
2,1284
Modelo com interação
Intercepto
-0,56 (0.1260)
-0,56 (0,1260)
-0,56 (0,1260)
Variedade
0,15 (0,2232)
0,15 (0,2232)
0,15 (0,2232)
Extrato raiz
1,32 (0,1775)
1,32 (0,1775)
1,32 (0,1775)
Interação
-0,78(0,3064)
-0,78 (0,3064)
-0,78 (0,3064)
φˆ
1,8618
1,8618
1,8618
FIGURA 1 - Representação gráfica da proporção média de sementes germinadas em função das
variedades O. aegyptiaco 73 e 75 e dos extratos de raízes de feijão e pepino.
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
63
A Figura 2 apresenta as proporções observadas e correspondentes valores ajustados
pelo modelo “M2”, de (2), obtidos pela inversa da função de ligação logística
correspondente, respectivamente:
M 1 : πˆ =
M 2 : πˆ =
(
(
)
exp αˆ 0 + αˆ1 x1(ijk )+αˆ 2 x2(ijk )
1 + exp αˆ 0 + αˆ1 x1(ijk )+αˆ 2 x2(ijk )
(
)
)
exp α 0 + α1 x1(ijk ) + α 2 x 2(ijk ) + α12 x1(ijk ) x 2(ijk )
1 + exp αˆ 0 + αˆ1 x1(ijk ) + αˆ 2 x 2(ijk ) + αˆ12 x1(ijk ) x 2(ijk )
(
)
Crowder (1978), citado por Breslow e Calyton (1993), observou que existe uma
variação dentro de parcelas que excede à predita pelos modelos acima ajustados,
caracterizando o fenômeno da superdispersão.
Com base na evidência anunciada, partiu-se para o ajuste de modelos com inclusão
do efeito aleatório de parcela no preditor linear (10) como forma de modelar a variação
extra-binomial por um MLGM. Devido à dificuldade na convergência do algoritmo,
optou-se por re-expressar os ensaios binomiais na forma de ensaios de Bernoulli (1’s e
0’s) (Littell et al., 1996).
A Tabela 2 apresenta as estimativas de máxima quase-verossimilhança penalizada
dos coeficientes de regressão logística obtidas pela macro Glimmix-SAS 8.10, pelo
algoritmo MLGM implementado na linguagem R 1.6.1 e pelo comando glmmPQL do R
1.6.1 para os efeitos fixos dos modelos generalizados mistos com ausência e presença da
interação dos fatores (11) e as estimativas dos componentes de dispersão referentes ao
parâmetro de extra-dispersão e ao efeito aleatório da parcela.
Observa-se que as estimativas obtidas pela macro Glimmix do SAS 8.10 e algoritmo
MLGM apresentam os mesmos resultados para ambos os modelos, com e sem interação,
denotando uma boa caracterização do algoritmo para fins de ajuste de um MLGM via
quase-verossimilhança penalizada conjunta, tendo assumido pressuposição distribucional
de normalidade para o efeito aleatório presente.
Os resultados obtidos pelo comando glmmPQL do software R 1.6.1 apresentaram
resultados similares com relação aos dois outros procedimentos de análise, por utilizar
uma implementação diferente a partir da QVP (8). Porém, as estimativas obtidas pelo
glmmPQL foram bem próximas das estimativas obtidas por Breslow e Clayton (1993) e
Hinde e Demétrio (1998) via máxima verossimilhança, demonstrando uma robustez dos
estimadores de QVP. Observa-se, ainda, que os erros padrões das estimativas obtidas pelo
glmmPQL do R 1.6.1 foram menores do que os obtidos pelo glimmix – SAS 8.10 e pelo
algoritmo MLGM implementado.
Para todos os algoritmos de análise foram calculadas as estimativas do parâmetro
referente ao componente de extra-dispersão (σ2), de (13), para fins de verificação da
acomodação da superdispersão (Shall, 1991), e a estimativa do parâmetro relacionado
com a variância do efeito aleatório de parcela (σ2p).
Para o modelo “M1”, de (11), as estimativas alcançadas para o algoritmo MLGM e
SAS 8.10 foram iguais a 0,9805 e 0,1242, enquanto as obtidas no R 1.6.1 foram de 0,9871
e 0,0857, respectivamente, para os componentes de extra-dispersão e componentes de
variância da parcela. Como os valores dos componentes de extra-dispersão estão bem
64
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
Tabela 2 - Estimativas dos parâmetros para os MLGM’s descritos em (11) por três
algoritmos para MLGM
ALGORITMOS
glmmPQL – R
ALGORITMO - MLGM
Glimmix – SAS
Variável
βˆ ( EP)
Intercepto
Variedade
Extrato raiz
βˆ ( EP)
βˆ ( EP)
Modelo com efeitos principais
-0,38 (0,1821)
-0,38 (0,1643)
-0,36 (0,2284)
-0,34 (0,2085)
1,02 (0,2236)
1,02 (0,2030)
-0,38 (0,1821)
-0,36 (0,2284)
1,02 (0,2236)
σˆ 2p
0,1242
0,0857
0,1242
σ̂
0,9805
0,9871
0,9805
2
Intercepto
Variedade
Extrato raiz
Interação
Modelo com interação
-0,54 (0,1904)
-0,54 (0,1656)
0,08 (0,3081)
0,10 (0,2746)
1,34 (0,2699)
1,33 (0,2347)
-0,83 (0,4296)
-0,81 (0,3817)
-0,54 (0,1904)
0,08 (0,3081)
1,34 (0,2699)
-0,83 (0,4296)
σˆ 2p
0,09780
0,0551
0,09780
σ̂ 2
0,9855
0,9899
0,9855
próximos a 1, valor admitido e fixo para este no modelo binomial, indicam que a variância
está consistente com a distribuição assumida, isto é, os dados não fornecem evidência de
superdispersão para o modelo assumido.
Para o modelo “M2”, de (11), considerando o efeito de interação, as estimativas
alcançadas para os componentes de extra-dispersão e de variância da parcela pelo
algoritmo MLGM e Glimmix-SAS 8.10 foram similares e iguais a 0,9855 e 0,0978,
enquanto as obtidas no R foram de 0,9899 e 0,0551, respectivamente. A inclusão do efeito
da interação, a qual foi significativa no modelo, aproximou a estimativa de σ2 a 1,
evidenciando melhor acomodação da variação extra-binomial.
A Figura 2 apresenta as proporções observadas e correspondentes valores ajustados
pelo modelo “M2”, de (11), a partir da inversa da função de ligação logística por meio das
estimativas obtidas pelo glimmix SAS 8.10 e algoritmo MLGM (A) e pelo comando
glmmPQL do programa R 1.6.1(B), conforme apresentado abaixo:
M 1 : πˆ ijk =
M 2 : πˆ ijk =
(
(
exp αˆ 0 + αˆ1 x1(ijk ) + αˆ 2 x2(ijk ) + bˆk
(
)
1 + exp αˆ 0 + αˆ1 x1(ijk ) + αˆ 2 x2(ijk ) + bˆk
)
)
exp αˆ 0 + αˆ1 x1(ijk ) + αˆ 2 x2(ijk ) + αˆ12 x1(ijk ) x2 (ijk ) + bˆk
1 + exp αˆ + αˆ x
+ αˆ x
+ αˆ x
x
+ bˆ
(0
1 1(ijk )
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
2 2(ijk )
12 1(ijk ) 2 (ijk )
k
)
65
A Figura 2 mostra que os valores das proporções ajustadas pelo MLGM “M2”, de
(11), em ambos algoritmos utilizados, ajustou-se de forma mais adequada aos dados
observados.
FIGURA 2 - Proporções de sementes germinadas e correspondentes proporções ajustadas pelo
modelo MLG “M2”, de (2), e pelo MLGM “M2”, de (11): SAS 8.10 e algoritmo
implementado (A) e comando glmmPQL do programa R 1.6.1 (B).
Exemplo 2: Ensaio de sobrevivência
A Tabela 3 apresenta algumas de muitas soluções possíveis para as estimativas dos
parâmetros para o modelo binomial via máxima verossimilhança, de (3), para as 20
proporções de ovos sobreviventes em função das variáveis classificatórias local e período,
por meio de programas executados nos pacotes estatísticos SAS 8.10 e R 1.6.1.
Verifica-se que a estimativa de φ, obtida de (5), para as análises realizadas, foi de
5,3303, excedendo o valor assumido para o modelo (φ = 1), evidenciando presença
de superdispersão. O valor da “deviance” residual, de (4), para o modelo com 12 graus de
liberdade, foi 64,495, sendo superior ao quantil do qui-quadrado com 12 graus de
liberdade para o nível de significância adotado (α = 0,05), rejeitado-se a hipótese de
ajustamento do modelo.
66
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
Tabela 3 - Estimativas dos parâmetros para o modelo binomial padrão (3) por três
algoritmos para ajuste do MLG
Variável
Intercepto
Local 01
Local 02
Local 03
Local 04
Local 05
Período (4 sem.)
Período (7 sem.)
Período (8 sem.)
Período (11 sem.)
φˆ
ALGORITMOS
GLM – R
ALGORITMO – MLG
GENMOD – SAS
βˆ ( EP)
βˆ ( EP)
βˆ ( EP)
-2,43 (0,1919)
4,61 (0,2502)
4,20 (0,2317)
3,37 (0,2014)
3,66 (0,2131)
0
2,45 (0,2341)
0,28 (0,1640)
0,12 (0,1648)
0
4,64 (0,2810)
0
-0,42 (0,2461)
-1,24 (0,2194)
-0,95 (0,2287)
-4,61 (0,2500)
0
-2,17 (0,2381)
-2,33 (0,2456)
-2,45 (0,2338)
1,00 (0,0502)
1,65 (0,1630)
1,23 (1,1445)
0,40 (0,1119)
0,69 (0,1244)
-2,99 (0,1439)
1,99 (0,1684)
-0,18 (0,1118)
-0,34 (0,1144)
-0,46 (0,1029)
5,3303
5,3303
5,3303
A Figura 3 mostra as proporções de ovos de truta observados e ajustados pelo
modelo (3), obtidos pela inversa da função ligadora denotada por
πˆij =
(
)
exp µˆ + αˆ i + βˆ j
1 + exp µˆ + αˆ + βˆ
(
i
j
),
para as 20 unidades experimentais. Denota-se que o modelo ajustado não descreve de
forma satisfatória os dados deste experimento.
FIGURA 3 - Proporções de ovos de truta sobreviventes observados e ajustados pelo modelo (3)
para as 20 parcelas.
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
67
Em virtude da inadequação do modelo proposto (3) em explicar os dados, propôs-se,
analogamente, modelar esta variabilidade por meio de um MLG com inclusão do efeito
aleatório da parcela experimental, conforme (12).
Com base neste modelo (12), as estimativas de máxima quase-verossimilhança
penalizada dos parâmetros para os efeitos fixos de local e período resultaram em iguais
proporções ajustadas para as análises realizadas pela macro glimmix do SAS 8.10 e
algoritmo MLGM e estas foram próximas às obtidas por meio das estimativas
provenientes pelo comando glmmPQL do programa R 1.6.1, conforme mostrado na
Tabela 4.
A Figura 4 apresenta a plotagem das proporções de ovos de truta sobreviventes
observados e ajustados pelo modelo (12), obtidos analogamente pela inversa da função de
ligação
πˆ ij =
(
)
exp µˆ + αˆ i + βˆ j + bˆij
1 + exp µˆ + αˆ + βˆ + bˆ
(
i
j
ij
)
com substituição pelas estimativas dos parâmetros correspondentes ao modelo. Nota-se
que o modelo misto proposto explica de forma bastante satisfatória a variabilidade
inerente aos dados experimentais.
Tabela 4 - Estimativas dos parâmetros para o MLGM descrito em (12) por três algoritmos
para MLGM
Variável
Glimmix – SAS
ALGORITMOS
GLM – R
Algoritmo – MLGM
βˆ ( EP)
βˆ ( EP)
βˆ ( EP)
Intercepto
Local 01
Local 02
Local 03
Local 04
Local 05
Período (4 sem.)
Período (7 sem.)
Período (8 sem.)
Período (11 sem.)
-2,87 (0,5877)
4,61 (0,6401)
4,42 (0,6404)
3,61 (0,6230)
4,09 (0,6433)
0
2,59 (0,5996)
0,59 (0,5534)
0,48 (0,5502)
0
4,36 (0,4367)
0
-0,25 (0,4375)
-1,07 (0,4190)
-0,69 (0,4303)
-4,54 (0,4321)
0
-2,01 (0,4109)
-2,11 (0,4080)
-2,47 (0,4131)
0,96 (0,1392)
1,46 (0,4065)
1,26 (0,4058)
0,45 (0,3888)
0,94 (0,4081)
-3,15 (0,4043)
1,91 (0,3827)
-0,09 (0,3410)
-0,19 (0,3378)
-0,67 (0,3460)
σˆ 2p
0,6550
0,2330
0,6550
σ̂ 2
0,9097
0,9392
0,9097
Estimaram-se, ainda, os componentes de dispersão ou variância relativos à parcela e
à variação extra (13) iguais a 0,655 e 0,9097, respectivamente, para as análises realizadas
no SAS 8.10 e algoritmo MLGM. A análise oriunda do glmmPQL do programa R 1.6.1
resultou em estimativas equivalentes a 0,2330 e 0,9392 para estas componentes. Nota-se
que a estimativa da componente de variância para o efeito aleatório de parcela para as
68
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
análises diferem substancialmente devido à diferente implementação utilizada a partir de
(8), utilizada pelo comando glmmPQL.
Com base no critério sugerido por Shall (1991), aproximado por (13), tem-se que o
valor estimado para a componente de extra-dispersão foi, relativamente, próximo a 1, o
que nos permite aferir que o modelo misto comporta de forma satisfatória a variação
extra-binomial presente.
FIGURA 4 - Proporções de ovos de truta sobreviventes observados e ajustados pelo modelo (12)
para as 20 parcelas.
Conclusões
A aplicação dos modelos lineares generalizados mistos para fins de acomodação da
variação extra-binomial presente nos dados referentes aos ensaios de germinação e
sobrevivência considerados constitui uma técnica bastante adequada;
O algoritmo implementado em ambiente R 1.6.1 mostrou-se consistente com
softwares de renomada precisão, nos dois exemplos utilizados para o ajuste de modelos
lineares generalizados mistos com ligação logística e pressuposição de normalidade para o
efeito aleatório presente.
NUNES, J.A.R.; MORAIS, A.R. de; BUENO FILHO, J.S. de S. Modelling overdispersion
in binomial data by a generalized linear mixed model. Rev. Mat. Estat., São Paulo, v.22,
n.1, p.55-70, 2004.
ABSTRACT: In agronomic research, experimental data with binomial distribution are found
frequently. In model selection, the presence of extra-binomial variation needs to be considered.
The use of generalized linear mixed models (GLMM) is a flexible methodology to analyze
proportions as well as an elegant way to model such overdispersion. The objective of this study
was to apply this methodology in two experiments involving binomial data with the purpose to
model the extra-binomial variability to an appropriate and valid inference. Furthermore, it
aimed at implementing an algorithm in R environment to fit GLMM with the logit link function in
order to obtain maxima penalized quasi-likelihood estimators assuming normally distributed
random effects. The estimates of the dispersion parameters using standard binomial fit to both
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004
69
data sets were greater than one, giving strong evidence of overdispersion. Analysis using a
generalized linear model with a random effect to plots was performed in SAS 8.1 by the
“glimmix” macro, R 1.6.1 through the “glmmPQL” command (Library MASS) and the
implemented algorithm in R 1.6.1 environment. Results showed that the suggested models
explained the present variability in a very acceptable way. Estimates of the components of extra
dispersion close to one provided no indication of overdispersion. The estimates achieved by the
implemented algorithm for fitting GLMM were similar to the estimates of the SAS program,
denoting that it can be used to fit mixed models for binomial data using the logit link function.
KEYWORDS: Algorithm; extra-binomial variation; penalized quasi-likelihood; deviance;
germination trial.
Referências
BRESLOW, N.E.; CLAYTON, D.G. Approximate inference in generalized linear mixed
models. J. Am. Stat. Assoc., Washington, v.88, n.421, p.9-25, 1993.
DEMÉTRIO, C.G.B. Modelos lineares generalizados em experimentação agronômica. In:
REUNIÃO ANUAL DA RBRAS, 46.; SEAGRO, 9., 2001. Piracicaba. Resumos...
Piracicaba: ESALQ/USP, 2001. 113 p.
DOBSON, A.J. An introduction to generalized linear models. London: Chapman and
Hall, 1990. 173 p.
GILMOUR, A.R.; ANDERSON, R.D.; RAE, A.L. The analysis of binomial data by a
generalized linear mixed model. Biometrika, London, v.72, n.3, p.593-9, 1985.
HENDERSON, C.R.; KEMPTHORNE, O.; SEARLE, S.R.; KROSIGK, C.M.. The
estimation of environmental and genetic trends from records subject to culling.
Biometrics, London, v.15, n.1, p.192-218, 1959.
HINDE, J.P.; DEMÉTRIO, C.G.B. Overdispersion: models and estimation. Comp. Stat.
Data Anal., v.27, n.2, p.151-70, 1998.
IHAKA, R.; GENTLEMAN, R.R: A language for data analysis and graphics. J. Comp.
Graphical Stat., Alexandria, v.5, n.3, p.299-314, 1996.
LITTLEL, R.C.; MILLEKEN, G.A.; STROUP, W.W.; WOLFINGER, R.D. SASSystem
for mixed models. Cary: SAS Institute, 1996. 633 p.
McCULLAGH, P.; NELDER, J. A. Generalized linear model. 2.ed. London: Chapman
and Hall, 1989. 511 p.
McCULLOCH, C. E.; SEARLE, S. R. Generalized, linear, and mixed models. New York:
E. Willey-Interscience, 2001. 324 p.
NELDER, J.A.; WEDDERBURN, R. W.M.Generalized linear model. J. R. Stat. Soc. A,
London, v.135, n.3, p.370-84, 1972.
SHALL, R. Estimation in generalized linear models with random effects. Biometrika,
London, v.78, n.4, p.719-27, 1991.
Recebido em 23.04.2003.
Aprovado após revisão em 01.11.2003.
70
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004