Spanish (pdf)
Article in xml format
Article references
Automatic translation
Send this article by e-mail
Cited by SciELO 
Similars in
SciELO 
uBio 
Permalink
INTRODUCCIÓN
La investigación ha producido numerosos modelos matemáticos para describir la curva de lactación bovina (Beever et al., 1991), la cual típicamente consta de una fase ascendente, un pico de producción y una fase descendente (Pollot, 2000). Dos de los modelos más usados son la función de Wood (1967): y = β xβ2e-β3x y la función de Wilmink (1987): y = β +β eβ3x+ β x. Para ambas funciones, y representa la producción diaria de leche y x son los días de lactación. Los parámetros β1, β2 y β3 de Wood y β1, β2 y β4 de Wilmink están relacionados, respectivamente, con el nivel de producción, el incremento de la producción antes del pico y la subsecuente disminución. El parámetro β de Wilmink está asociado con el tiempo hasta alcanzar el pico de producción y usualmente es reemplazado por una constante (Quinn et al., 2005; Torshizi et al., 2011). Es posible estimar el tiempo en el pico de producción: xmax=β2/β3 y la producción en el pico: ymax=β1(β2/β3)β2e-β2 para Wood y xmax=1/β3log(β4/β2β3) e ymax=β1-(β4/β3)log(β4/β2β3)+(β4/β3) para Wilmink.
Los parámetros β1 son convencionalmente estimados a nivel de hato mediante soluciones iterativas de las ecuaciones normales obtenidas por mínimos cuadrados bajo asunciones de normalidad, homocedasticidad e independencia de errores (Draper y Smith, 1981). En este contexto de modelos de efectos fijos no lineales (MNL), las funciones de Wood y Wilmink han sido extensivamente usadas no solo en el análisis de curvas de lactación de bovinos (Scherchand et al., 1995; Scott et al., 1996), sino también de ovinos (Ángeles et al., 2014), y caprinos (Gipson y Grossman, 1990). En el país, la función de Wood ha sido usada para describir el efecto del número y estación de partos sobre la lactación promedio de vacas en la cuenca lechera de Lima (Rodríguez et al., 2005) y para modelar la curva de lactación de vacas Gyr y cruces F-1 Gyr x Holstein en la región San Martín (Huamán et al., 2018).
El uso de los MNL para la estimación de los parámetros promedio de un gran número de lactaciones no plantea muchas dificultades a las asunciones mencionadas. Sin embargo, la creciente necesidad de definir el manejo alimentario y sanitario y la evaluación genética a nivel de sujetos ha puesto un énfasis en el modelamiento de las curvas individuales de lactación, en la variación de sus parámetros, y en la asociación de esta variación con efectos aleatorios y con efectos sistemáticos internos (raza, edad, número de partos, etc.) y externos (alimentación, estación de parto, etc.) (Scott et al., 1996).
El modelamiento de curvas individuales de lactación plantea un escenario diferente, caracterizado por una alta variabilidad, heterocedasticidad, y errores correlacionados (Macciotta et al., 2005). Adicionalmente, las curvas individuales suelen estar basadas en pocos registros, los cuales frecuentemente entre lactaciones y su relación con factores fijos y aleatorios. Todo lo anterior plantea la necesidad de modelos predictivos con la suficiente flexibilidad como para incorporar heterocedasticidad y correlación de errores, y con la capacidad de describir la variabilidad dentro y entre lactaciones. Los modelos mixtos no lineales (MMNL) parecen satisfacer la mayoría de estos requerimientos (Vonesh y Carter, 1992; Davidian y Giltinan, 1995; Serroyen et al., 2009). A pesar de esta conveniencia, son pocas las referencias en el uso de MMNL para curvas de lactación. Quintero et al. (2007) en Colombia, examinaron siete funciones, entre ellas la de Wood y Wilmink, para curvas de lactación de búfalas, con incorporación de efectos aleatorios solo para los coeficientes lineales, pero sin considerar la heterocedasticidad ni la correlación de errores. Vásquez (2017) probó cuatro funciones (Wood y Wilmink incluidas) usando MNL y MMNL para describir las curvas de lactación de vacas estabuladas del valle de Huaura, Perú.
Los MMNL, también conocidos como modelos jerárquicos no lineales, son un caso particular de los modelos mixtos, donde, en un contexto multinivel, se ajusta el modelo a nivel de individuos, caracterizando la variación dentro de ellos (o dentro de lactaciones, en este caso) a través de una estructura específica de la covarianza individual. A nivel de población, la variación entre lactaciones es representada a través de modelos de regresión no lineal específicos para cada lactación, los cuales pueden incorporar efectos aleatorios y sistemáticos (Lindstrom y Bates, 1990; Davidian y Giltinan, 1995).
Los objetivos del presente estudio fueron (i) evaluar la capacidad de los MMNL para describir curvas de lactación bovina obtenidas a partir de registros del establo de la Estación IVITA El Mantaro (Junín, Perú), (ii) comparar las funciones de Wood y de Wilmink en términos de su bondad de ajuste y en un contexto de variabilidad entre lactaciones, heteroscedasticidad y errores correlacionados y (iii) estimar el efecto del número de partos y la estación del año en la que se produjo el parto sobre los parámetros de las curvas de lactación.
MATERIALES Y MÉTODOS
Población
La población inicial estuvo conformada por 257 lactaciones correspondientes al periodo 2004-2012, provenientes de vacas predominantemente Holstein, Brown Swiss y algunas con trazas de Jersey, de 1 a 3 o más partos, del Establo San Juan de la Estación Experimental IVITA El Mantaro de la Facultad de Medicina Veterinaria de la Universidad Nacional Mayor de San Marcos en la cuenca lechera del Centro, Junín, Perú. Las vacas fueron manejadas al pastoreo, con suplementación alimentaria y confinamiento parcial durante el ordeño.
Tamaño de la muestra
Las lactaciones fueron clasificadas según estación de parto (lluviosa: octubre-marzo; seca: abril-setiembre) y número de partos (1, 2 y 3 o más partos), los cuales son dos de los factores que más afectan las características de una lactación (Rekik et al., 2003; Gloria et al., 2012). Se seleccionaron en forma aleatoria siete lactaciones por cada uno de los seis grupos, conformando una muestra de 42 lactaciones.
Recolección de Datos
Los datos fueron recolectados de los registros del establo en forma de producción diaria de leche (kg) para cada lactación e incorporados en la hoja de cálculo Microsoftâ Excel, como archivos con la extensión csv. Para el análisis de las curvas se incluyeron los registros de los días 5 y 15 después del parto, de allí cada 15 días hasta el día 60 y luego cada 30 días hasta el final de la lactación. Se excluyeron aquellas lactaciones con un primer registro posterior a los 5 días, con menos de 12 registros o aquellas con una extensión mayor a 420 días. Los datos fueron importados por el programa R© 3.5.3 (R Core Team, 2019) mediante la instrucción read.table.
Ajuste de MNL
Para proporcionar una base de comparación se ajustaron las funciones de Wood y Wilmink usando MNL. Los estimados de los parámetros se obtuvieron por mínimos cuadrados mediante el algoritmo iterativo de Gauss-Newton (Draper y Smith, 1981) usando la función nls de R. Los valores iniciales de los parámetros se obtuvieron mediante solución de las formas lineales, log y = log β1 + β log x + β x, para la función de Wood e y = β + β (0.914x) + β x para la función de Wilmink. Para esta última, el coeficiente β fue reemplazado por la constante -0.09, obtenida comparando el ajuste lineal del modelo para diferentes posibles valores β variando entre -0.04 y -0.1 y eligiendo la que producía el menor residuo.
Ajuste de MMNL
El ajuste de los MMNL fue secuencial (Pinheiro y Bates, 2000). Primero, para definir la forma de las matrices de diseño y de covarianza de efectos aleatorios se ajustó el MNL a cada una de las 42 lactaciones por separado por medio de la función nlsList del módulo nlme de R y se examinaron las distribuciones de los intervalos de confianza (IC )0.95 de los parámetros de los modelos. Aquellos parámetros cuyos IC0-95 mostraron superposición completa o casi completa fueron modelados como efectos fijos; de lo contrario fueron modelados como efectos fijos más efectos aleatorios. Para estimar la matriz de covarianza de los efectos aleatorios, se ajustó un MMNL sin covariables, homoscedástico, con errores independientes y con la forma de la matriz de efectos aleatorios sugerida preliminarmente. El ajuste se efectuó por el método de Máxima Verosimilitud usando Linealización Condicional de Primer Orden (Lindstrom y Bates, 1990; Davidian y Giltinan, 1995; Pinheiro y Bates, 2000), con la ayuda de la función nlme de R. El examen de los IC0.95 de los efectos aleatorios y sus correlaciones determinó la forma final de la matriz de covarianza efectos aleatorios.
En segundo lugar, la asunción de homocedasticidad se verificó a través del examen del diagrama de dispersión de los residuos estandarizados con respecto a las respuestas ajustadas (y). La presencia de un patrón sistemático de dispersión sugiere la necesidad de incorporar la heteroce-dasticidad en el MMNL. Esta incorporación se llevó a cabo modelando la varianza como una función de potencia de la media: g(µ θ) = µθ , usando el argumento weights=varPower() de la función nlme, disponible en R. En una tercera etapa, con el fin de flexibilizar la asunción de independencia de errores dentro de lactaciones, se incorporó en el MMNL una matriz de correlación espacial exponencial (Jones y Ackerson, 1990; Pinheiro y Bates, 2000; Diggle et al., 2002), previo examen de los semivariogramas muestrales de los residuos estandarizados.
Los principales parámetros para estimarse con la matriz de correlación espacial son el rango y el efecto «nugget». El rango es el periodo de tiempo a lo largo del cual la semivarianza aumenta o, equiva-lentemente, la correlación entre errores disminuye. El «nugget» expresa la semivarianza de los errores con una separación temporal de cero y es interpretada como aquella variabilidad entre los errores no explicada por su posición (usualmente asignada al error de medición) (Diggle et al., 2002). La matriz de correlación espacial fue incorporada en el MMNL mación de Akaike (AIC) (Akaike, 1974) y bayesiano (BIC) (Schwarz, 1978) y el log verosimilitud. El incremento de la verosimilitud como producto de los sucesivos refinamientos del MMNL fue evaluado mediante la prueba de Razón de Verosimilitud. El nivel de significación empleado a lo largo del estudio fue 0.05. mediante el argumento corr=cor Exp-(form=~x,nugget=TRUE) de la función nlme de R.
Las etapas de la secuencia fueron aplicadas, individual y comparativamente, tanto a la función de Wood como a la de Wilmink. Al final de cada etapa se examinó el efecto de la modificación del MMNL sobre expresiones de bondad de ajuste como el error estándar residual (σ), los criterios de información de Akaike (AIC) (Akaike, 1974) y bayesiano (BIC) (Schwarz, 1978) y el log verosimilitud. El incremento de la verosimilitud como producto de los sucesivos refinamientos del MMNL fue evaluado mediante la prueba de Razón de Verosimilitud. El nivel de significación empleado a lo largo del estudio fue 0.05.
RESULTADOS Y DISCUSIÓN
La Figura 1 muestra las trayectorias que conectan la producción de leche de las observaciones repetidas de cada una de las 42 lactaciones. Es discernible que (i) a pesar de la alta variabilidad, tanto en nivel de producción, como en extensión de la lactación, las trayectorias parecen tener un patrón similar, y colectivamente se distribuyen alrededor de la curva típica descrita previamente; (ii) las lactaciones con alta producción al inicio tienden a mantenerla a lo largo de la trayectoria y viceversa, en una suerte de «seguimiento» o «tracking», característico de datos longitudinales (Diggle et al., 2002); y (iii) la dispersión en la producción de leche tiende a ser menor al inicio de la lactación que al final.
Figura 1.Trayectorias de producción de leche de 42 lactaciones de vacas predominantemente Holstein y Brown Swiss de un establo en Junín, Perú (2004-2012)
El ajuste del MNL a los datos de lactación produjo las ecuacion e s y = 12.55x0.17e-0.0038x para Wood e y = 22.8510.37e-0.09x 0.038x para Wilmink. Todos los parámetros fueron estadísticamente significativos. Los criterios de bondad de ajuste fueron muy similares para Wood y Wilmink: σ 4.27 y 4.24, AIC 3448.57 y 3441.51, BIC 3466.16 y 3459.10 y log Verosimilitud - 1720.29 y -1716.16, respectivamente. La Figura 2 muestra la dispersión de los valores de producción de leche por tiempo y por lactación y su tendencia a aglomerarse alrededor de las curvas de Wood y Wilmink. La casi perfecta concurrencia en el comportamiento de ambas funciones refleja la similaridad en los criterios de bondad de ajuste.
Figura 2 Dispersión de datos de producción de leche de 42 lactaciones de vacas predominantemente Holstein y Brown Swiss de un establo en Junín, Perú (2004-2012) con curvas de los modelos mixtos no lineales de Wood (y = 12.55x0.17e-0.0038x) y de Wilmink (y = 22.85-10.37e-0.09x 0.038x)
El Cuadro 1 resume los resultados de la secuencia de ajuste de los MMNL en términos de los parámetros estimados y de los criterios de bondad de ajuste. El MMNL homocedástico y con errores independientes fue analizado con efectos aleatorios para todos los parámetros, y con una matriz de covarianza de los efectos aleatorios general definida positiva, tal como fue sugerido por el examen de la distribución de los IC0.95 de los parámetros y por los IC0.95 de los componentes de esta matriz. Este ajuste no modificó sustancialmente los estimados de los parámetros βi ni sus errores estándar con respecto al MNL; sin embargo, redujo en un 60% los σ y mejoró los criterios AIC, BIC y log Verosimilitud. Esta mejora fue confirmada por las pruebas de razón de verosimilitud: L.ratio(10,4)=739.94 (p<0.0001) para Wood y L.ratio(10,4)=841.56 (p<0.0001) para Wilmink, las cuales evidencian que la inclusión de los efectos aleatorios incrementaron en forma significativa la verosimilitud de los modelos.
Cuadro 1 Estimados de los parámetros de las funciones de Wood (y = β1xβ2e-β3x) y Wilmink (y = β1 + β2e-0.09x + β3x) e incrementos en la bondad de ajuste como producto de la incorporación consecutiva de heteroscedasticidad y correlación de errores en la secuencia de ajuste de modelos mixtos no lineales para la descripción de curvas de lactación
Los p-valores corresponden al incremento de la verosimilitud del modelo pertinente con respecto al modelo precedente. Errores estándar entre paréntesis
1Error estándar residual
2Criterio de información de Akaike
3 Criterio de información bayesiano
4 Modelo mixto no lineal
El examen de la dispersión de residuos con respecto a las y sugirió el ajuste de un MMNL heterocedástico, el cual produjo θ = 0.372 para la función de Wood y θ = 0.327 para la función de Wilmink, lo cual sugiere que las varianzas residuales no son homogéneas sino que varían en proporción a µ 2(0.372) y a µ 2(0.327), respectivamente. La incorporación de estas funciones de varianza redujo σ en aproximadamente 61%, mejoró los criterios de bondad de ajuste e incrementó significativamente la verosimilitud del modelo para ambas funciones, con respecto al MMNL homocedástico (Cuadro 1). Se considera que funciones de varianza de este tipo no son intuitivas para datos de lactación; donde, para una población heterogénea, la producción de leche es completamente variable desde el inicio de la lactancia. Bajo este mismo argumento, Macciotta et al. (2011) señalan que la evolución de la varianza es más bien ondulante a lo largo de la lactación y Jaffrezic et al. (2000) mencionan que las mayores varianzas residuales se muestran al inicio y al final de la lactación. Desafortunadamente, no se ha encontrado información publicada sobre funciones específicas de varianza para curvas de lactación con la cual referida por otros autores (Quinn et al., 2005; Torshizi et al., 2011). En el contexto de los MMNL, Vásquez (2017) obtuvo que la función de mejor ajuste para lactaciones de un solo parto fue la de Wood. La de Wilmink, en cambio, se comportó mejor en lactaciones de dos a más partos. En términos de parámetros derivados tampoco parecen haber muchas diferencias, con ymax de 20.3 kg para Wood y 21.7 kg para Wilmink. Esta última función, sin embargo, predice un ymax más temprano (37.4 vs 46.9 días). Similar comportamiento ha sido observado por Torshizi et al. (2011) (66 vs 81 días) y por Vásquez (2017) (65 vs 83 días). La semejanza entre ambas funciones, sumada a ciertas ventajas de Wilmink como su estructura aditiva y fácil linealización, justifica el rol de esta función como alternativa a Wood (Macciotta et al., 2005).
Contrariamente a lo esperado, ninguna de las covariables estación o número de partos o sus interacciones tuvo un efecto significativo sobre los parámetros de las funciones, con excepción del efecto «3 partos vs 1 parto» sobre el coeficiente estimado β de Wilmink. Para verificar este resultado, se evaluó la inclusión aislada del efecto de nú número de partos sobre β3 de Wilmink, la cual no incrementó significativamente la verosimilitud del MMNL [L.ratio(11,13)=2.87 (p=0.24)]. La ausencia de efecto de las covariables sugiere que los efectos aleatorios estuvieron exclusivamente representados por la variabilidad individual y otros factores no controlados en este estudio; entre ellos el genético (Gloria et al., 2012), la carga animal (MacDonald et al., 2008) o la extensión de la lactación (Grossman y Koops, 2003)
Dos criterios gráficos de la bondad de ajuste de los MMNL son la comparación de las curvas de Wood y Wilmink por lactación y para efectos fijos (Pinheiro y Bates, 2000) y el diagrama de dispersión de las respuestas observadas con respecto a las respuestas ajustadas (Gauch et al., 2003). La Figura 3 muestra la superposición de las curvas de Wood y de Wilmink, heterocedásticas y con errores correlacionados, tanto para los efec tos fijos como para las predicciones individuales, de una muestra de 6 de las 42 lactaciones. Las curvas predichas para las lactaciones individuales se ajustan a los datos observados y tienden a distribuirse en forma balanceada alrededor de la curva de efectos fijos. En la mayoría de las curvas individuales y particularmente en aquellas que no parecen ajustarse a los datos, existe una suerte de «retracción» (shrinkage) hacia la curva de efectos fijos. El efecto, común en modelos mixtos, es atribuido a la agregación parcial de la información de todas las lactaciones para producir los estimados de cada efecto aleatorio individual (Combes et al., 2013) y presta robustez a la predicción de los parámetros de cada lactación. Esta «retracción» suele ser directamente proporcional al error de predicción e inversamente proporcional al número de observaciones por individuo (Xu et al., 2012).
Figura 3 Curvas de lactación de efectos fijos (línea discontínua) y de efectos aleatorios (Línea continua) para los Modelos Mixtos no Lineales de Wood y Wilmink correspondientes a una muestra de 6 de 42 lactaciones de vacas predominantemente Holstein y Brown Swiss de un establo en Junín, Perú (2004-2012). Los puntos son los datos observados
La aglomeración de los puntos alrededor de la línea de equivalencia (Figura 4) expresa la concordancia entre los valores observados y los ajustados y atestigua la pertinencia del MMNL para ambas funciones.
CONCLUSIONES
Los datos de producción láctea del Establo San Juan de IVITA El Mantaro entre los años 2002 y 20012 describen un sistema altamente variable, de baja productividad y con lactaciones extendidas, característico de un sistema bajo pastoreo.
No parecen existir diferencias considerables en el desempeño de las funciones de Wood y Wilmink, a excepción de pequeñas pero consistentes diferencias en los criterios de bondad de ajuste a favor de esta última, así como la tendencia de la función de Wilmink a predecir picos de producción más tempranos.
No hubo una contribución significativa del número de partos ni de la estación de partos a la variabilidad de los efectos aleatorios para ninguna de las funciones estudiadas.
