MATERIALES Y MÉTODOS RESULTADOS Y DISCUSIÓN CONCLUSIONES BIBLIOGRAFÍABAYESIAN PHENOTYPIC STABILITY ANALYSIS USING JEFFREYS’S PRIOR
José Miguel Cotes Torres1 y Adhemar Sanches2
1 Profesor Asistente. Universidad Nacional de Colombia, Sede Medellín. Facultad de Ciencias Agropecuarias. A.A. 1779, Medellín, Colombia. <jmcotes@unalmed.edu.co>
2 Profesor Departamento de Ciencias Exatas. Universidade Estadual Paulista. Facultade de Ciencias Agrárias e Veterinárias de Jaboticabal (FCAV/UNESP). <adhesan@fcav.unesp.br>
Recibido: Marzo 11 de 2005; aceptado: Octubre 13 de 2005.
RESUMEN
]]> Uno de los métodos utilizados para evaluar estabilidad fenotípica es el propuesto por Shukla, el cual calcula la varianza de los genotipos dentro de la interacción genotipo por ambiente, para lo cual se hace uso de la estimación de los componentes de varianza dentro del análisis de varianza combinado. El acceso al cálculo de metodologías como REML e ML permitieron trabajar con datos que presentan algún grado de desbalance, sin embargo no solucionan de una manera adecuada el problema de la estimación de componentes de varianza negativos, los cuales son asumidos como cero y redistribuidos en los demás componentes positivos. El uso de la metodología bayesiana en la estimación de componentes de varianza resuelve satisfactoriamente este problema sin afectar los demás componentes. En este trabajo, se utilizaron datos de producción comercial de papa de 10 pruebas regionales realizadas en la región andina colombiana y se utilizó la metodología bayesiana en la solución del modelo mixto para la estimación de la varianza de Shukla con base en una distribución a priori no informativa de Jeffreys. Fueron obtenidas muestras de la distribución a posteriori conjunta mediante el algoritmo Independence Chain, con un tamaño de muestreo de 1,16x105 y un burnin de 500. Los resultados muestran que en la estimación REML de componentes de varianza tres genotipos presentan componentes de varianza estimados como cero. Las estimativas bayesianas son 89,35; 377,18 y 101,12; y los respectivos intervalos de credibilidad al 95% son: (2,13 – 371,70), (35,26 – 1363,67) y (2,33 – 434,53). Finalmente con estas estimativas no se afectó la estimación de los demás componentes de varianza.Palabras claves: Varianza de Shukla, priori de Jeffreys, importance sampling.
Shukla's variance is a very useful method for the analysis of phenotypic stability, computing the genotypic variance among the genotype by environment interaction, using variance component estimation of combined analysis of variance. New methodologies like REML or ML allow work with unbalanced data but do not have a good solution for the negative variance estimate problem. In this case, the component is taken as zero and the rest of the variance is redistributed into other components with positive estimates. The use of Bayesian methodology in the variance component estimation resolves satisfactory this problem without affecting the other components. This research uses the Bayesian estimation methodology for the solution of the mixed model in order to obtain the Shukla's variance for potato production data of 10 regional trials established in the Colombian Andean Region. The noninformative Jeffreys’s prior and Independence Chain algorithm were used. The burnin period was 500 iterations and 1,16x105 generations of joint posterior distribution were obtained. The REML methodology found three Shukla's variances with zero estimates. The corresponding Bayesian estimates were 89,35; 377,18 and 101,12 and the 95% confidence intervals were (2,13 – 371,70); (35,26 – 1363,67) and (2,33 – 434,53), respectively. Finally, these estimates do not affect other variance estimate components.
Key words: Shukla’s variance, Jeffrey's prior, importance sampling.
MATERIALES Y MÉTODOS
RESULTADOS Y DISCUSIÓN
CONCLUSIONES
BIBLIOGRAFÍA
La interacción genotipo por ambiente (IGA) es una de las principales preocupaciones de los mejoradores por dos razones (Kang y Magari 1996): primero porque reduce el progreso de selección (Kang y Gorman 1989, Cooper y DeLacy 1994 y Delacy, Cooper y Basford 1996) y segundo porque hace imposible interpretar los efectos principales (debidos exclusivamente a los genotipos o al ambiente).
]]> Para el análisis de la IGA se han desarrollado varias metodologías, tales como: Análisis de Regresión (Finlay e Wilkinson 1963, Eberhart y Russell 1966 y Cruz y Regazzi 1994), Estimación de Componentes de Varianza (Shukla 1972) y Análisis Multivariado (Gauch e Zobel 1989). Debido a su facilidad de cálculo el Análisis de Regresión ha sido la metodología más utilizada destacándose los métodos de Finlay e Wilkinson 1963 y Eberhart y Russell 1966, aunque hayan recibido numerosas críticas debido a que la parte predictiva de la regresión no explica satisfactoriamente la IGA.En los estudios mas recientes las metodologías mas utilizadas son las que usan estimación de componentes de varianza, ya que ha sido posible implementar metodologías computacionales de estimación diferentes de aquellas que usan Mínimos Cuadrados, como es el caso de la metodología de máxima verosimilitud o ML (“Maximum Likelihood”), y máxima verosimilitud restricta o REML (“Restricted Maximum Likelihood”), las cuales mejoran la precisión de las estimativas.
Aunque estas metodologías presentan varia ventajas respecto a la estimación tradicional de componentes de varianza, existe el problema de que cuando por lo menos un componente de varianza es estimado como nulo, el análisis resultante es sesgado (Searle, Casella y Mcculloch 1992), ya que la variación total es redistribuida en los componentes de varianza no nulos. Así en estos casos una alternativa viable es utilizar la metodología de estimación bayesiana de componentes de varianza (Box y Tiao 1972, Searle, Casella y Mcculloch 1992 y Wolfinger y Kass 2000), usando preferiblemente distribuciones a priori no informativas (Box y Tiao 1972, Sun et al. 1996, Pauler, Wakefield y Kass 1999, y Wolfinger y Kass 2000).
La metodología bayesiana se basa principalmente en la obtención de la distribución a posteriori en función de las distribuciones a priori para cada uno de los componentes de varianza a estimar y de la función de verosimilitud. Comúnmente se define como estimador de Bayes la media o la mediana de la distribución a posteriori. Si se escoge la media se minimiza el riesgo de Bayes según la función de pérdida cuadrática, y si se elige la mediana entonces se minimiza el riesgo de Bayes de acuerdo con la función de perdida absoluta (Mood, Graybill y Boes 1974).
Cualquiera que sea el estimador de Bayes seleccionado, se necesita obtener las distribuciones marginales para cada uno de los componentes de variancia y calcular su respectiva media o mediana. Para este proceso es necesario integrar la distribución a posteriori, que frecuentemente es multidimensional y de difícil obtención por métodos analíticos (Searle, Casella y Mcculloch 1992). Por esta razón se hace necesario el uso de métodos numéricos computacionalmente intensivos como lo son los métodos de integración Monte Carlo. Estos métodos se basan en la simulación de muestras de la distribución a posteriori, las cuales son utilizadas para calcular cualquier característica de interés de las distribuciones a posteriori marginales, como lo pueden ser la media, la mediana, la moda, los intervalos de credibilidad (IC) y los intervalos de alta probabilidad (HPD), entre otras.
Para casos multidimensionales, existen algoritmos especiales que obtiene muestras (métodos de muestreo por simulación) de la distribución a posteriori, tales como (a) algoritmo de Gibbs (“Gibbs sampling”) y (b) algoritmo de Metropolis-Hastings (Gelman et al. 1995, Sorensen 1996 y Gamerman 1997). Tales algoritmos son iterativos y se basan en propiedades de las Cadenas de Markov por lo que se han denominado métodos de Monte Carlo en Cadenas de Markov, abreviadamente MCMC (“Monte Carlo Markov Chain”).
Otros métodos de integración Monte Carlo que no usan Cadenas de Markov, son los algoritmos de Cadenas Independientes y el algoritmo de Importance Sampling. Estos se basan en simular datos de una distribución que sea muy próxima de la distribución marginal a posteriori deseada, y por una determinada regla de aceptación y/o rechazo, se obtienen muestras reales de la distribución deseada. La distribución base para obtener las muestras candidatas, es escogida por ser de fácil simulación, sin embrago estos algoritmos funcionan mejor en la medida que esta distribución base se aproxima más a la distribución marginal a posteriori deseada (Ripley 1987 y Wolfinger y Kass 2000).
Así, en esta investigación se utilizó la estimación bayesiana de los componentes de varianza de estabilidad fenotípica de Shukla a través de los recursos computacionales disponibles en el procedimiento SAS/MIXED.
]]> Los datos del presente trabajo se obtuvieron de las pruebas regionales del programa de mejoramiento en papa del convenio UNIPAPA, liderado por la Facultad de Agronomía de la Universidad Nacional de Colombia sede Bogotá, las cuales se realizaron en los años de 1998 y 1999 en 10 localidades, seleccionadas estas por su facilidad de acceso y por encontrarse en sitios estratégicos de la producción de papa en el país.
Se utilizó un diseño de bloques completos al azar con tres repeticiones en cada localidad. Se utilizaron como testigos las 4 variedades comerciales Diacol Capiro, Parda Pastusa, ICA Morita e ICA Única y se evaluaron 11 genotipos promisorios. Las pruebas regionales realizadas, no tuvieron el mismo número de genotipos debido a la falta de disponibilidad de semilla. La variable evaluada fue rendimiento de papa total (ton/ha). El modelo estadístico utilizado fue (Magari y Kang 1997).
donde, 
es el vector del número total de observaciones, 
son respectivamente la matriz y el vector de efectos (parametrizados como m+bk) para los genotipos (1 £ k £g) donde g es el número de genotipos, 
son respectivamente la matriz y el vector de efectos (ai) para ambientes (1 £ i £ a) donde a es el número de ambientes, 
son respectivamente la matriz y el vector de efectos (gij) para replicaciones j dentro de i ambientes, 
son respectivamente la matriz y el vector de efectos para la interacción genotipo por ambiente del k-ésimo genotipo y 
es el vector de efectos del error experimental. Se suponen que ambientes y replicaciones son efectos aleatorios y los genotipos son efectos fijos.
Según la definición de modelos jerárquicos se define la distribución a posteriori conjunta como:
Para este modelo la distribución condicional a posteriori para los efectos fijos y aleatorios es una normal multivariada con media,
y varianza,
es proporcional al producto de, 
donde L(.) representa la función de verosimilitud integrada o restricta exclusivamente a los componentes de varianza y p(.) es la distribución a priori para cada componente de varianza.
En este estudio se utilizó la regla de Jeffreys para la obtención de la distribución a priori no informativa para cada uno de los componentes de varianza. Así, para el modelo en consideración,
es proporcional a la raíz cuadrada de la matriz de información de Fisher, esto es:
Para obtener muestras simuladas de la distribución a posteriori conjunta 
se utilizó el algoritmo de Independence Chain (Wolfinger y Kass 1987).
Finalmente se utilizó el procedimiento SAS/MIXED para la obtención de la muestras según el método de Independece Chain eliminando las primera 500 muestras simuladas (burnin=500), y obteniendo un total de 1,16 E5 muestras de la distribución a posteriori conjunta. Como estimativa de Bayes fue considerada la media de la distribución a posteriori marginal para cada componente de varianza. También fueron programadas rutinas auxiliares en SAS/IML, SAS/GPLOT y SAS/KDE, para obtener la densidad a posteriori marginal, la estimativa de Bayes y el respectivo intervalo de credibilidad para cada componente de varianza.
]]> RESULTADOS Y DISCUSIÓN
La estimación de los componentes de varianza muestra marcadas diferencias entre los métodos REML y Bayesiano. En el primer caso se obtienen estimativas nulas para la varianza de estabilidad fenotípica de Shukla para tres de los 15 genotipos evaluados (Tabla 1). También se puede observar como estas estimativas afectan las de los demás componentes de varianza. Así, por ejemplo, para el genotipo 1 la estimativa REML es cero y la estimativa bayesiana es 89,35; por otro lado para el genotipo 3 la estimativa REML es positiva y es de 41,77 y la estimativa bayesiana es de 429,38; siendo esta última mayor. Esto está de acuerdo con las limitaciones del método REML, el cual al enfrentarse a una estimación de varianzas negativas, el componente de varianza respectivo es colocado como cero y la variación total es de nuevo redistribuida en los otros componentes de varianza (Searle, Casella y Mcculloch 1992), dando como resultado una subestimación de los parámetros. Además de esto se aprecia que el genotipo 3 presenta un límite superior del intervalo de confianza al 95 % de 1,18 x 1013, el cual es un valor exagerado y poco realista para este genotipo. En la metodología bayesiana se ve como este genotipo presenta un intervalo de credibilidad al 95 % entre 67,30 y 1381,57.
Desde el punto de vista del fitomejorador, utilizando la metodología Bayesiana, es posible considerar los genotipos uno, dos y 10 como estables ya que presentan el menor valor estimado para la varianza de Shukla, lo cual puede ser fácilmente observado en la Figura 1. También se puede notar que, los genotipos más inestables serian los genotipos 6 y 8, ya que presentan una estimativa de Bayes mayor e intervalos de credibilidad más amplios. Llama la atención el hecho de que las variedades testigos, pueden considerarse como más inestables que la mayoría de los nuevos genotipos evaluados. Sin embargo, las variedades Parda Pastusa y Diacol Capiro presentan una mayor estabilidad al ser comparadas con las variedades ICA Unica e ICA Morita (Tabla 1 y Figura 2).
Figura 1. Función de densidad a posteriori marginal para cada uno de los 11 genotipos promisorios de papa
Tabla 1. Estimativas de la variancia de estabilidad fenotípica de Shukla a través de la metodología bayesiana y REML, en genotipos de papa.
Figura 2. Función de densidad a posteriori marginal para las variedades testigo de papa.
Otra de las dificultades que presenta el método REML es la realización de una prueba de hipótesis adecuada para cada uno de los componentes de varianza. Frecuentemente es calculado el error estándar de cada componente y realizada una prueba de hipótesis asumiendo una distribución Normal asintótica. Sin embargo, para que esta prueba de hipótesis sea válida se necesita un elevado número de ambientes. Así, al analizar en las Figuras 1 y 2, se deduce que la distribución marginal a posteriori para la varianza de Shukla con 10 ambientes dista mucho de una Normal asintótica. Este problema no se presenta para el componente de varianza de bloques dentro de ambientes, debido al elevado número de grados de libertad que presenta, sin embargo este componente de varianza no es de interés en la selección de genotipos. Finalmente este problema no se presenta en la metodología bayesiana, en la cual, para comparar formalmente dos componentes de varianza se puede utilizar el Factor de Bayes (Kass y Raftery 1995), o la distribución a posteriori de las diferencias o de la razón de los dos componentes de varianza involucrados (Box y Tiao 1972).
]]>La metodología bayesiana permite realizar una mejor aproximación al problema de estimación de componentes de varianza, ofreciendo varias herramientas para el análisis de estabilidad fenotípica, y evitando varios de los problemas que se presentan al utilizar otras metodologías estadísticas como REML.
Según los valores de estabilidad fenotípica obtenidos se recomiendan la selección de los genotipos de papa uno, dos y 10.
Box, G.E.P. and Tiao, G.C. 1973. Bayesian inference in statistical analysis. New York: Wiley Classics Library. 608 p. [ Links ]
Cooper, M. and Delacy, I. H. 1994. Relationships among analytical methods used to study genotypic variation and genotypeby environment interaction in plant breeding multienvironment experiments. En: Theoretical and Applied Genetics. Vol.88, no. 2; p. 561-572. [ Links ]
Cruz, C. D. e Regazzi, A. J. 1994. Modelos biométricos aplicados ñao melhoramento genético. Viçosa – Minas Gerais, Brasil: UFV-Imprensa Universitaria. 390 p. [ Links ]
Delacy, I. H., Cooper, M. and Basford, K. E. 1996. Relationships among analytical methods used to study genotype-by-environment interactions and evaluation of their impact on response to selection. p. 51-84. En: Genotype-by-environmet interaction. New York: CRC. [ Links ]
Eberhart, S. A. and Russell, W. A. 1966. Stabillity parameters for comparing varietes. En: Crop Science. Vol. 6; p. 36-40. [ Links ]
Finlay, K. W. and Wilkinson, G. N. 1963. The analysis of adaptation in a plant breeding programme. En: Australian Journal of Agriculture Research. Vol. 14, no. 6; p. 742-754. [ Links ]
Gauch, H. G. and Zobel R. W. 1989. Predictive and posdictive success of statistical analises of yield trials. En: Theoretical and Applied Genetics. Vol. 76, no. 1; p. 1-10. [ Links ]
Gamerman D. 1997. Markov chain Monte Carlo: stochastic simulation for bayesian inference. London. U.K: Chapman and Hall. 245 p. [ Links ]
Gelman, A., Carlin, J. B., Stern, H. S. and Rubin, D. B. 1995. Bayesian data analysis. Great Britain: Chapman and Hall. 526 p. [ Links ]
Kang, M. S. and Gorman D. P. 1989. Genotype x environment interactions in maize. En: Agronomy Journal. Vol. 81, no. 2; p. 662-664. [ Links ]
________ and Magari, R. 1996. New development in selecting for phenotypic stability in crop breeding. p. 1-14. En: Genotype-by-environment interaction. New York: CRC. [ Links ]
Kass, R. E. and Raftery A. E. 1995. Bayes factor. En: Journal of the American Statistical Association Vol. 90, no. 430; p. 773-796. [ Links ]
Magari R. and Kang M. S. 1997. SAS-STABLE: Stability analises of balanced an unbalanced data. En: Agronomy Journal. Vol. 89, no. 7; p. 929-932. [ Links ]
Mood, A. M., Graybill, F. A. and Boes, D. C. 1974. Introduction to the theory of statistics. 3 ed. Tokyo, Japan : McGraw Hill. 564 p. [ Links ]
Pauler D. K., Wakefield, J. C. and Kass, R. E. 1999. Bayes factor and approximations for variance component models. En: Journal of the American Statistical Association. Vol. 94, no. 448; p. 1242-1253. [ Links ]
Ripley, B. 1987. Stochastic simulation. New York: John Wiley. 237 p. [ Links ]
Searle, S. R., Casella, G., and Mc-culloch, C. E. 1992. Variance components. New York: John Wiley. 501 p. [ Links ]
Shukla, G. K. 1972. Some statistical aspects of partitioning genotype environment components of variability. En: Heredity. Vol. 29, no. 2; p. 237-245. [ Links ]
Sorensen, D., 1996. Gibbs sampling in quantitative genetics. Tjele, Denmark : National Institute of Animal Science (Internal Report; no. 82) [ Links ]
Sun, L., Hsu, J. S. J., Guttman, I and Leonard, T. 1996. Bayesian methods for variance component models. En: Journal of the American Statistical Association. Vol. 91, no. 10; p. 743-752. [ Links ]
Wolfinger, R. D. and Kass, R. E. 2000. Nonconjugate bayesian analysis of variance component models. En: Biometrics. Vol. 56, no. 3; p. 768-774. [ Links ] ]]>