16
Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 55 MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM MODELO LINEAR GENERALIZADO MISTO José Airton Rodrigues NUNES 1 Augusto Ramalho de MORAIS 2 Júlio Sílvio de Sousa BUENO FILHO 2 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]. 2 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].

MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Embed Size (px)

Citation preview

Page 1: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 55

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]. 2 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].

Page 2: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 56

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

Page 3: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 57

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 é

( ) π=RE

βπ

πη X=��

���

−=

1log

(1)

em que: R é o vetor de proporções observadas, de dimensões n x 1, tal que se admite:

i

iii m

mBinR

),(~

π, com i = 1,...,n;

η é 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:

Page 4: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 58

( ) ( )

( ) ( ) ( ) ( )ijkijkijkijkijk

ijkijk

ijkijkijk

ijkijk

xxxxM

xxM

211222110

22110

1log:2

1log:1

ααααπ

πη

αααπ

πη

+++=��

��

−=

++=��

��

−=

(2)

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

jiij

ijij βαµ

ππ

η ++=��

��

−=

1log (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

( ) ( ) ( )[ ]rlrrlrD ,ˆ,2ˆ, ππ −=

( ) ( )���

����

����

−−

−+���

����

�=

− i

ii

i

ii

n

ii

rr

rrmrD

πππ

ˆ11

log1ˆ

log2ˆ,1

(4)

Por conseguinte, realizou-se o teste da adequação do modelo corrente com p

parâmetros linearmente independentes por meio da estatística 2( ; )n 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

Page 5: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 59

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

ˆ ˆ( )' ( )ˆ Wn p

φ − −=−

z zη ηη ηη ηη η (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 extra-binomial 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

variáveis aleatórias i

iiii m

mBinuR

),(~|

π, com função de densidade discreta dada por

3 GILMOUR, S. G. 2002, Lavras - MG, comunicação pessoal.

Page 6: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 60

( )iiRii urfiiduRui

|.~||

( ) ( )��

��

��

��

���

����

�+

��

�−+

��

��

−=

ii

ii

i

iiiiiR rm

mrmurf

iuilog1log

1logexp|

ππ

(6)

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

( )

�uX�

uR

+==��

���

=

ηπ

ππ

1log

|E (7)

em que: R|u é o vetor de proporções condicionalmente independentes, de dimensões n x 1, tal que

se admite i

iiii m

mBinuR

),(~|

π, correspondente a (6);

η é 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

( ) ( ) uGurmruQVP ii

ii

ii

1'

21

1log1

log,, −−���

�−+��

����

−= π

πππ (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

( )πη −∆+= ry* (9)

tem-se que esta pode ser aproximada por um modelo linear misto (Henderson et al., 1959) da forma

** eXy +Ζ+= β

Page 7: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 61

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

��

���

�=

��

� −

+ −

*'

*'

'

'

'

'

ˆ

ˆ 1

WyX

WyZ

WZX

GWZZ

WXX

WZZu

β (10)

em que, para um MLGM para dados binomiais, tem-se que

( )���

���

−=

��

���

��

���

���

����

����

−∂∂=

���

���

∂∂

=∆iii

i diagdiagdiagπππ

πηη

η11

1log

11

1

( )[ ] ( )[ ]( ) ( ){ }iii

ii

iiiii

i

i mdiagm

diaguyVdiagW ππππππ

ππ

−=��

���

��

���

−−

=��

��

��

��

���

����

∂∂

= − 11

1|

21

2

.

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 quase-logverossimilhanç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:

( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) kijkijkijkijkkijk

kijkijkkijk

bxxxxbx

bxxbx

++++==Μ

+++==Μ

211222110,

22110,

|1Prlogit:2

|1Prlogit:1

αααα

ααα (11)

em que:

Page 8: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 62

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:

( ) kjikjiij bb +++=== βαµβαη ,,|1Prlog (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 ( 2σ̂ ) conforme descrito por Shall (1991). O teste para a

superdispersão presente foi realizado verificando se os valores de 2σ̂ 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

( ) ( )( )Xrn

uZXzWuXz−

−−Ζ−−=ˆˆ'ˆˆ

ˆ 2 ββσ (13)

em que, ( )Xr é 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 ( )0023,0Pr => D , 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 ( ) 5871,2705.0;172 =X ,

Page 9: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 63

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

GENMOD – SAS GLM – R ALGORITMO – MLG Variável ˆ ( )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.

Page 10: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 64

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:

( ) ( )( )( ) ( )( )ijkijk

ijkijk

xx

xxM

22110

22110

ˆˆˆexp1

ˆˆˆexpˆ:1

αααααα

π+

+

+++

=

( )( ) ( ) ( ) ( )( )ijkijkijkijk

ijkijkijkijk

xxxx

xxxxM

211222110

)(2)(112)(22)(110

ˆˆˆˆexp1

expˆ:2

αααααααα

π++++

+++=

����

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 (σ2

p). 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

Page 11: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 65

Tabela 2 - Estimativas dos parâmetros para os MLGM’s descritos em (11) por três algoritmos para MLGM

ALGORITMOS Glimmix – SAS glmmPQL – R ALGORITMO - MLGM Variável

ˆ ( )EPβ ˆ ( )EPβ ˆ ( )EPβ

Modelo com efeitos principais Intercepto -0,38 (0,1821) -0,38 (0,1643) -0,38 (0,1821) Variedade -0,36 (0,2284) -0,34 (0,2085) -0,36 (0,2284)

Extrato raiz 1,02 (0,2236) 1,02 (0,2030) 1,02 (0,2236) 2ˆ pσ 0,1242 0,0857 0,1242 2σ̂ 0,9805 0,9871 0,9805

Modelo com interação Intercepto -0,54 (0,1904) -0,54 (0,1656) -0,54 (0,1904) Variedade 0,08 (0,3081) 0,10 (0,2746) 0,08 (0,3081)

Extrato raiz 1,34 (0,2699) 1,33 (0,2347) 1,34 (0,2699) Interação -0,83 (0,4296) -0,81 (0,3817) -0,83 (0,4296)

2ˆ pσ 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:

( ) ( )( )( ) ( )( )kijkijk

kijkijkijk

bxx

bxxM ˆˆˆˆexp1

ˆˆˆˆexpˆ:1

22110

22110

++++

+++=

αααααα

π

( ) ( ) ( ) ( )( )( ) ( ) ( ) ( )( )kijkijkijkijk

kijkijkijkijkijk

bxxxx

bxxxxM ˆˆˆˆˆexp1

ˆˆˆˆˆexpˆ:2

211222110

211222110

+++++

++++=

αααα

ααααπ

Page 12: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 66

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.

Page 13: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 67

Tabela 3 - Estimativas dos parâmetros para o modelo binomial padrão (3) por três algoritmos para ajuste do MLG

ALGORITMOS GENMOD – SAS GLM – R ALGORITMO – MLG Variável

ˆ ( )EPβ ˆ ( )EPβ ˆ ( )EPβ Intercepto -2,43 (0,1919) 4,64 (0,2810) 1,00 (0,0502) Local 01 4,61 (0,2502) 0 1,65 (0,1630) Local 02 4,20 (0,2317) -0,42 (0,2461) 1,23 (1,1445) Local 03 3,37 (0,2014) -1,24 (0,2194) 0,40 (0,1119) Local 04 3,66 (0,2131) -0,95 (0,2287) 0,69 (0,1244) Local 05 0 -4,61 (0,2500) -2,99 (0,1439) Período (4 sem.) 2,45 (0,2341) 0 1,99 (0,1684) Período (7 sem.) 0,28 (0,1640) -2,17 (0,2381) -0,18 (0,1118) Período (8 sem.) 0,12 (0,1648) -2,33 (0,2456) -0,34 (0,1144) Período (11 sem.) 0 -2,45 (0,2338) -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

( )( )ji

jiij βαµ

βαµπ ˆˆˆexp1

ˆˆˆexpˆ

+++

++= ,

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.

Page 14: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 68

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

( )( )ijji

ijjiij

b

bˆˆˆˆexp1

ˆˆˆˆexpˆ

++++

+++=

βαµ

βαµπ

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

ALGORITMOS Glimmix – SAS GLM – R Algoritmo – MLGM Variável

ˆ ( )EPβ ˆ ( )EPβ ˆ ( )EPβ

Intercepto -2,87 (0,5877) 4,36 (0,4367) 0,96 (0,1392) Local 01 4,61 (0,6401) 0 1,46 (0,4065) Local 02 4,42 (0,6404) -0,25 (0,4375) 1,26 (0,4058) Local 03 3,61 (0,6230) -1,07 (0,4190) 0,45 (0,3888) Local 04 4,09 (0,6433) -0,69 (0,4303) 0,94 (0,4081) Local 05 0 -4,54 (0,4321) -3,15 (0,4043) Período (4 sem.) 2,59 (0,5996) 0 1,91 (0,3827) Período (7 sem.) 0,59 (0,5534) -2,01 (0,4109) -0,09 (0,3410) Período (8 sem.) 0,48 (0,5502) -2,11 (0,4080) -0,19 (0,3378) Período (11 sem.) 0 -2,47 (0,4131) -0,67 (0,3460)

2ˆ pσ 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

Page 15: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 69

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

Page 16: MODELAGEM DA SUPERDISPERSÃO EM DADOS POR UM …jaguar.fcav.unesp.br/RME/fasciculos/v22/v22_n1/A5_JAirton.pdf · binomial efetuou-se a modelagem dos dados correspondentes pelo uso

Rev. Mat. Estat., São Paulo, v. 22, n.1, p.55-70, 2004 70

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. SASSystem 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.