Previsão Bayesiana de valores genéticos de touros por meio do modelo auto-regressivo para dados em painel

Bayesian forecasting of sires breeding values using autoregressive panel data model

ZOOTECNIA E TECNOLOGIA E INSPEÇÃO DE PRODUTOS DE ORIGEM ANIMAL

F.F. Silva^I; T. Sáfadi^II; J.A. Muniz^II; L.H. Aquino^II; G.B. Mourão^III

^IDepartamento de Informática - UFV, Av. P.H. Rolfs, s/n, 36570-000 - Viçosa, MG

^IIDepartamento de Ciências Exatas - UFLA - Lavras, MG

^IIIDepartamento de Ciências Exatas - ESALQ - Piracicaba, SP

RESUMO

Inferência Bayesiana do modelo auto-regressivo foi aplicado para dados em painel a dados reais de DEP´s de touros da raça Nelore publicadas entre 2000 e 2005. Na análise foram consideradas a função de verossimilhança exata e a obtenção de distribuições preditivas de dados futuros. Verificou-se ser recomendável o agrupamento dos animais em grupos homogêneos de acordo com a acurácia. Constatou-se também, que em média, a eficiência de previsão dos valores de DEP´s para um ano futuro é próxima de 80%.

Palavras-chave: gado Nelore, série temporal, método MCMC, DEP

ABSTRACT

A Bayesian inference of autoregressive panel data model was applied to real data of Nelore sires Expected Progenie Difference (EPD) during a five-year period (2000-2005). The exact likelihood function and predictive distributions of future observations were considered. The results indicated the importance of sires grouping in homogeneous groups according to accuracy and showed forecast efficiency for EPD values in a future year around 80%.

Keywords: Nelore cattle, time serie, MCMC method, EPD

INTRODUÇÃO

Os estudos de dados em estrutura de painel estão essencialmente orientados para estudar a heterogeneidade relativa a diferentes indivíduos que apresentam uma trajetória longitudinal. Esse fato permite identificar determinados aspectos, como tendências e sazonalidades, que são de difícil quantificação usando apenas modelos de regressão individuais.

Uma classe de modelos de regressão bastante retratada em estudos de dados em painel é a de séries temporais. Um modelo de série temporal é utilizado quando se dispõe de uma seqüência de dados obtidos em intervalos regulares de tempo durante um período específico. Dentre esses modelos, destaca-se o auto-regressivo, pois ele se aplica a diversas situações práticas e geralmente apresenta boa qualidade de ajuste comparado aos modelos mais complexos. De acordo com Morettin e Toloi (2004), ao se estudarem esses modelos, primeiramente há necessidade de se compreender o comportamento da série, mediante a obtenção de estimativas para os parâmetros e, depois, realizar previsões fundamentadas nessas estimativas.

Em recentes estudos envolvendo análise de modelos de regressão e de séries temporais (Barreto e Andrade, 2004; Silva et al., 2005; Silva et al., 2006) a metodologia Bayesiana foi utilizada com sucesso, pois sua característica de considerar todos os parâmetros como variáveis aleatórias, reduziu substancialmente o número de estimativas viesadas e possibilitou a utilização de um número menor de observações, pois os conceitos probabilísticos envolvidos independem do número de dados utilizados. Outro fato que merece destaque é que na metodologia Bayesiana a estimação por intervalo é realizada de forma direta, considerando a incerteza existente sobre todos os parâmetros simultaneamente, sendo, portanto, a estimação por intervalo mais precisa em relação àquela apresentada pela metodologia freqüentista, que, em geral, utiliza variâncias assintóticas e aproximadas (O'Hagan, 1994).

A inferência Bayesiana consiste de uma informação a priori, dos dados amostrais e do cálculo da densidade a posteriori dos parâmetros. A informação a priori é dada pela densidade de probabilidade, que expressa o conhecimento do pesquisador sobre os parâmetros a serem estimados. Os dados representam amostra aleatória de uma população com densidade, e são utilizados na análise Bayesiana por meio da função de verossimilhança f, que representa a distribuição conjunta desses dados.

Portanto, a partir do momento em que se opta por uma distribuição a priori, obtém-se a função de verossimilhança, torna-se possível, por meio do Teorema de Bayes, representado pela fórmula (1), obter a distribuição densidade a posteriori de , de forma que qualquer conclusão a seu respeito é realizada a partir desta distribuição.

sendo Y={Y₁, Y₂,..., Y_n}. O denominador, chamado de constante de integração, não depende de 0. Portanto, tem-se a expressão,P(θ | Y) ∝ L (Y | θ )P (θ), que pode ser entendida como: posteriori ∝verossimilhança x priori.

A distribuição a posteriori de um parâmetro contém toda a informação probabilística a seu respeito. Assim, toda a inferência sobre o parâmetro é realizada por meio desta distribuição. Segundo Rosa (1998), para se inferir com relação a qualquer elemento de q, a distribuição a posteriori conjunta dos parâmetros, P(θ |Y) deve ser integrada em relação a todos os outros elementos que a constituem. Assim, se o interesse do pesquisador concentra-se sobre determinado conjunto de θ, por exemplo 0₁, tem-se a necessidade da obtenção da distribuição p(q₁ |Y), denominada de distribuição marginal, que é dada por:

A integração da distribuição conjunta a posteriori para a obtenção das distribuições marginais geralmente não é analítica, necessitando de algoritmos iterativos especializados como o Gibbs Sampler e o Metropolis-Hastings, os quais são denominados de algoritmos Markov Chain-Monte Carlo (MCMC). Portanto, para a utilização desses algoritmos, é necessário que se obtenha, a partir da distribuição a posteriori, um conjunto de distribuições chamadas de distribuições condicionais completas.

Ao se utilizar o modelo auto-regressivo de ordem p, que apresenta termos de defasagem temporal, verifica-se uma redução na quantidade de observações que compõem a função de verossimilhança quando esta se apresenta condicionada às p primeiras observações (Box et al., 1994), uma vez que essas não são consideradas na análise. Ao se tratar de dados em painel, esse problema é ampliado, pois a condicionalidade é inerente a todos os indivíduos, o que acarreta uma perda considerável de informações. Uma maneira de solucionar o problema é a designação de uma função de verossimilhança exata, ou incondicional, pois apesar de proporcionar uma maior complexidade ao processo de análise, permite considerar todas as observações (Dufour e Kiviet, 1998).

O objetivo central da análise Bayesiana de Séries Temporais é gerar mecanismos robustos de previsão. Uma observação futura é descrita por uma distribuição condicionada aos dados passados, denominada distribuição preditiva (Barreto e Andrade, 2004). Quando se dispõe da distribuição conjunta das previsões futuras pode-se facilmente obter a distribuição desses valores futuros. Por exemplo, pode-se estar interessado em realizar previsões a respeito de valores assumidos por uma determinada variável ao longo do ano de 2007, baseando-se em dados observados durante o ano de 2006.

Segundo Box et al. (1994) a grande utilização de modelos auto-regressivos é justificada pelas suas diversas aplicações em situações reais. Uma possível aplicação, considerando dados em painel, está relacionada com modelagem e previsão de valores genéticos, também denominados de diferenças esperadas nas progênies (DEP). Essa aplicação é possível porque as DEP's são calculadas e publicadas anualmente em sumários, e ao verificar que um mesmo touro tem seus valores genéticos calculados ao longo de vários anos, tem-se então uma série temporal. Ao se considerar vários touros, pode-se assumir a caracterização de um arquivo de dados estruturado em forma de painel.

Ao se promover seleção, além de ser importante ter uma idéia do que se espera com o processo, ou seja, a predição do resultado da seleção, é importante também obter uma estimativa da segurança que se terá ao se adotar determinado procedimento (Carneiro et al., 2001). Essa segurança é definida pela acurácia da seleção, e em um sumário de touros, têm-se indivíduos cujas DEP´s foram obtidas sob diferentes graus de acurácia, o que exige um agrupamento desses indivíduos em classes homogêneas para que essas possam ser consideradas as populações requeridas na análise de dados em painel.

O presente trabalho teve como objetivo apresentar a inferência Bayesiana do modelo auto-regressivo de ordem p, AR(p), para dados em painel, aplicada a dados de DEP´s calculadas durante os anos de 2000 a 2005 para a característica ganho de peso pós-desmama de touros da raça Nelore. Tal metodologia foi utilizada com o intuito de realizar previsões desses valores para anos futuros de acordo com o grau de acurácia das DEP´s.

MATERIAL E MÉTODOS

O modelo auto-regressivo, denotado por AR(p), em que p é o número de parâmetros considerados, é apresentado pela fórmula (3):

em que:

Y_t é o valor atual de um processo estocástico, cujos valores já assumidos no passado são dados por Y_t-1, Y_t-2, ... , Y_t-p;

Φ₁, Φ₂,..., Φ_p são os parâmetros do modelo, denominados de parâmetros de auto-regressão;

e_t é o resíduo associado ao modelo, também denominado de ruído branco,

Na presença de dados em estrutura de painel, o modelo AR(p) indicado em (3) continua apresentando as mesmas características, porém a ele deve ser acrescentado um índice i referente ao indivíduo, como em (4):

em que: i=1,2,...m; j=1,2,...,p e t=1,2,...,n_i.

De acordo com esta notação, têm-se m indivíduos, cada um com n_i observações longitudinais, indicando que cada indivíduo i pode apresentar um número específico de observações. Neste estudo optou-se pela utilização do modelo AR(2), ou seja p=2, portanto o modelo contempla dois parâmetros por indivíduo. Em relação aos resíduos, assume-se que eles são amostras aleatórias provenientes de uma mesma distribuição, portanto

A função de verossimilhança exata, de acordo com o modelo apresentado em (4), é a seguinte:

em que:

Os vetores Φ e Y₁, considerados nas expressões, apresentam, respectivamente, as dimensões mp x m(n-p) e m(n-p) x 1.

Optou-se pela utilização de uma priori independente t-Student multivariada - Gama Inversa, respectivamente para os parâmetros Φ e , a qual é dada por:

Depois de obtidas as distribuições conjuntas a posteriori, por meio do teorema de Bayes, torna-se necessário apresentar as distribuições condicionais completas a posteriori para Φ e , sendo estas respectivamente:

Nota-se que a distribuição condicional completa para o parâmetro é dada por uma distribuição gama-inversa, ou seja, ela apresenta uma forma conhecida, portanto passível ao uso do algoritmo Gibbs Sampler. O mesmo não acontece para a distribuição condicional do parâmetro Φ, que não apresenta uma forma definida, requerendo nessa situação o uso do algoritmo Metropolis-Hastings.

Os algoritmos Gibbs Sampler e Metropolis-Hastings foram implementados matricialmente no software estatístico R (R Development Core Team, 2006 <tttp://www.R-project.org>). A função mnormt (multivariate Normal and t-Student distributions) foi utilizada para a geração de números aleatórios implícita nos algoritmos MCMC. Estes algoritmos foram implementados considerando uma cadeia de 20.000 iterações, de forma que as 10.000 primeiras foram eliminadas para reduzir o efeito dos valores iniciais (burn-in). Para avaliar suas convergências foram utilizados os critérios de Gelman e Rubin (1992) e Raftery e Lewis (1992) mediante o pacote Bayesian Output Analysis (BOA) do R.

Para realizar a previsão de um dado futuro, utilizou-se a teoria de distribuição preditiva descrita por Heckman e Leamer (2001), que é dada por:

em que: I é uma matriz quadrada de ordem mp x mp. Assim, o conjunto de valores gerados para esta distribuição normal multivariada, provenientes de cada q iteração dos algoritmos Metropolis-Hastings e Gibbs Sampler, constituêm a distribuição preditiva para um dado futuro, cuja estimativa, , é representada pela média dessa distribuição.

Neste estudo foram utilizados dados cedidos pelo Grupo de Melhoramento Animal da Faculdade de Zootecnia e Engenharia de Alimentos da USP, GMA/FZEA, com sede em Pirassununga, SP. Estes correspondem aos sumários de touros Nelore compreendidos entre os anos de 2000 e 2005. Constam no arquivo dados de DEP´s para a característica ganho de peso entre a desmama (205 dias de idade) e o sobreano (550 dias de idade) de 117 touros, divididos em três grupos, conforme o valor da acurácia referente a última observação, isto é, aos dados de 2005. Os grupos foram divididos em: acurácia baixa (0 a 40%), com 31 touros , acurácia média (41 a 60%), com 63 touros e acurácia alta (acima de 60%), com a freqüência de 23 touros. Considerou-se também a análise de todos os animais, sem a separação em grupos de acordo com a acurácia. Dessa forma, utilizaram-se, portanto, quatro arquivos de dados. Em todos eles, a última observação, referente ao ano de 2005, foi suprimida visando à comparação entre os valores estimados por meio da distribuição preditiva e os verdadeiros valores observados.

RESULTADOS E DISCUSSÃO

Um resumo dos resultados, expressos como porcentagem de significância dos parâmetros Φ₁ e Φ₂, obtidos no ajuste do modelo AR(2) aos dados individuais de DEP's considerando cada estruturada de dados determinada pelos valores da acurácia é apresentado na Tab. 1.

É possível inferir que as séries referentes à população de indivíduos com baixa acurácia apresentam um comportamento diferenciado das demais, pois a porcentagem de significância dos parâmetros foi menor. Nota-se também, que o parâmetro Φ₂ apresentou porcentagem de significância um pouco menor que aquela obtida para o parâmetro Φ₁ nas estruturas de dados referentes à acurácia média e geral, indicando que as séries de DEP's de alguns indivíduos não apresentam um comportamento auto-regressivo de segunda ordem, podendo essas serem descritas por uma ordem inferior, AR(1), ou ainda se apresentarem como um processo estacionário, que indica a ausência de autocorrelação.

Para complementar estes resultados, foram realizadas análises das funções de autocorrelação e autocorrelação parcial para as séries dos indivíduos cujos parâmetros Φ₁ e Φ₂ não se mostraram significativos. Por meio destas verificou-se que na estrutura relacionada com baixa acurácia, a função de autocorrelação parcial apresentou lag's significativos para alguns indivíduos, indicando que modelos de ordem superior, AR(3), devem ser considerados em 38% dos casos. Além disso, também foram identificados indivíduos (62%) cujas séries já são caracterizadas como um processo estacionário. A variedade de processos identificados nesta estrutura de dados não foi observada nos conjuntos relacionados com a acurácia média e alta, pois o número de indivíduos cujas séries demonstraram seguir um modelo AR(3) foram, respectivamente, 23,8% e 11,1%, observando-se, portanto, uma quase totalidade de processos estacionários.

Essa discussão evidencia o fato de a acurácia baixa proporcionar maior variabilidade entre as séries, o que justifica a separação dos indivíduos em grupos. Porém, não foram observadas evidências a respeito de diferenças na variabilidade do comportamento de séries representadas pelas estruturas de acurácia média e alta, o que pode indicar uma possível união desses indivíduos em uma só estrutura.

O ponto mais marcante em relação aos resultados Tab. 1 está relacionado com a identificação de indivíduos cujas densidades preditivas serão apresentadas, pois, não há razão em discutir o processo de previsão de valores futuros de séries individuais que não apresentam um processo auto-regressivo de primeira e segunda ordem, uma vez que a densidade preditiva está fundamentada no modelo AR(2) para dados em painel. Assim, serão discutidas sob o aspecto preditivo apenas séries cujos parâmetros mostraram-se significativos.

Nas Tab. 2, 3 e 4 são apresentados os verdadeiros valores da última observação (y₆), correspondentes as DEP's do ano de 2005, os quais não foram considerados na análise para possibilitar a comparação entre valores observados e estimados por via das densidades preditivas.

Os resultados permitem afirmar que a metodologia utilizada para realizar previsões de dados futuros individuais com base na obtenção das distribuições preditivas foi eficiente, uma vez que a porcentagem de intervalos de credibilidade que continham os verdadeiros valores das DEP´s referentes ao ano de 2005 variou entre 77,8% e 85,7%. Estes resultados também conduzem a discussões relacionadas com as classes de acurácia consideradas, pois nota-se que a capacidade preditiva do modelo é menor quando a acurácia é baixa, 77,8%, e para as acurácias médias e altas, respectivamente 83,3% e 85,7%. Este fato também é explicado mediante estimativas obtidas para a variância do erro () em cada uma das estruturas de dados consideradas, as quais são mostradas na Tab. 5.

Realmente nota-se que a variância do erro foi maior na estrutura de dados com baixa acurácia, e essa revelação apenas vem a confirmar discussões anteriores relacionadas com a maior variabilidade entre o comportamento das séries dos indivíduos considerados na amostra referente a esta classe de acurácia.

CONCLUSÕES

A análise de dados de DEP´s de touros Nelore, em sucessivos anos, resultou em boas estimativas para o valor de um dado futuro, obtidas a partir da distribuição preditiva. Esse processo possibilitou a detecção de diferenças nos comportamentos de séries referentes a indivíduos pertencentes a grupos distintos de acurácia, indicando que este fator deve ser levado em consideração ao predizer DEP´s em anos futuros. As diferenças relatadas não impõem barreiras para a não adoção da metodologia em dados provenientes de animais com baixa acurácia, visto que a porcentagem de acerto na previsão de um dado futuro para essa classe foi próxima a 80%.

Recebido em 27 de dezembro de 2006

Aceito em 5 de setembro de 2008

E-mail: [e-mail]

• BARRETO, G.; ANDRADE, M.G. Robust Bayesian approach for AR(p) models applied to streamflow forecasting. J. Appl. Stat. Sci, v.12, p.269-292, 2004.
• BOX, G.E.P.; JENKINS, G.M.; G.C. REINSEL (Eds). Time series analysis: forecasting and control. 3.ed. San Francisco: Holden-Day, 1994. 500p.
• CARNEIRO, A.P.S. ; TORRES, R.A. ; EUCLYDES, R.F. et al. Efeito da conexidade dos dados sobre a acurácia dos testes de progênie e performance. Rev. Bras. de Zootec, v.30, p.342-347, 2001.
• DUFOUR, J.M.; KIVIET, J.F. Exact inference methods for first-order autoregressive distributed lag models.Econometrica, v.66, p.79-104, 1998.
• GELMAN, A.; RUBIN, D.B. Inference from iterative simulation using multiple sequence. Stat. Sci., v.7, p.457-511, 1992.
• HECKMAN, J.; LEAMER, E. (Eds). Handbook of Econometrics Amsterdam: Elsevier, 2001. p.744.
• MORETTIN, P.A.; TOLOI, C.M.C. (Eds). Análise de séries temporais São Paulo: Edgard Blucher, 2004. 535p.
• O'HAGAN, A. (Ed). Kendall's advanced theory of statistics: Bayesian Inference. New York: Edward Arnold Press, 1994. 402p.
• RAFTERY, A.E.; LEWIS, S. How many interations in the Gibbs sampler? In: BERNARDO, J.M. (Ed). Bayesian Statistics Oxford University, 1992. p.763-773.
• ROSA, G.J.M. Análise Bayesiana de modelos mistos robustos via amostrador de Gibbs 1998. 57f. Tese (Doutorado) - Universidade de São Paulo, Piracicaba.
• SILVA, F.F.; MUNIZ, J.A.; AQUINO, L.H. et al. Abordagem Bayesiana da curva de lactação de cabras Saanen de primeira e segunda ordem de parto. Pesq. Agropec. Brás, v.40, p.27-33, 2005.
• SILVA, N.M.A.; MUNIZ, J.A.; SILVA, F.F. et al. Estudo de parâmetros de crescimento de bezerros Nelore por meio de um modelo de regressão linear: uma abordagem Bayesiana. Cien. Anim. Bras, v.7, p.57-65, 2006.

Previsão Bayesiana de valores genéticos de touros por meio do modelo auto-regressivo para dados em painel

Bayesian forecasting of sires breeding values using autoregressive panel data model

Datas de Publicação

• Publicação nesta coleção
  03 Dez 2008
• Data do Fascículo
  Out 2008

Histórico

• Aceito
  05 Set 2008
• Recebido
  27 Dez 2006