Introducción
En algunas investigaciones se registran entre las variables de interés el recuento de eventos, es decir el conteo de un evento de interés que ocurre en un tiempo o espacio dado; en algunas disciplinas estas variables se analizan utilizando el Modelo Lineal General (ML) con el propósito de expresar la relación entre la variable de respuesta llamada independiente y una o más variables explicativas. Este modelo asume linealidad, normalidad, homogeneidad en las varianzas e independencia en los errores aleatorios, (Draper y Smith 1981). Con el objeto de lograr que estos supuestos se cumplan se recurre a transformaciones de los datos para lograr la homogeneidad de varianzas (Kuehl 2000; Melo Martínez et al. 2007). Cuando no se logra la normalidad, se acude a la estadística no paramétrica para realizar el análisis (Zar 1996; Lehmann 2006). La dificultad con las transformaciones se da en términos de la interpretación de los resultados los cuales no corresponden a la escala inicial de medición de la variable de respuesta.
Los Modelos Lineales Generalizados (MLG), son una extensión del ML, para modelar variables que pertenezcan a la familia de distribución de probabilidad exponencial (Nelder y Wedderburn 1972); esta familia contiene las distribuciones: normal, binomial, Poisson y binomial negativa, entre otras. El MLG propuesto, en el presente estudio, está conformado por: i) una componente aleatoria dada por la variable de respuesta 𝑌𝑖 y su distribución de probabilidad; ii) una componente sistemática conformada por las variables explicativas 𝑋𝑖 expresadas como una función lineal conocida como el predictor lineal 𝜂𝑖; iii) una función de enlace denotada por g (.) que relaciona el valor esperado de la variable de respuesta μi con las variables explicativas 𝑋𝑖 de la siguiente manera: 𝑔(𝜇𝑖)=𝛽0 + 𝛽1𝑥𝑖1 + …‚ + 𝛽𝑝𝑥𝑖𝑝 𝑝𝑎𝑟𝑎 𝑖= 1‚ 2,…‚ 𝑛; las p variables del lado derecho de la expresión corresponden a las variables explicativas y pueden ser factores de estudio, covariables e interacciones entre estas. Los parámetros 𝛽0, 𝛽1, … + 𝛽𝑝 son estimados utilizando máxima verosimilitud. La diferencia entre el ML y el MLG está en que este último incluye una función de enlace g (.) que relaciona la componente aleatoria y la componente sistemática de tal manera que 𝑔(𝜇𝑖)=𝛽0 + 𝛽1𝑥𝑖1 + …‚ + 𝛽𝑝𝑥𝑖𝑝 y no como en el ML que 𝜇𝑖 = 𝛽0 + 𝛽1𝑥𝑖1 + …‚ + 𝛽𝑝𝑥𝑖𝑝 ; de esta forma, la función de enlace lo que garantiza es que los valores estimados tomen valores en la misma escala de los esperados de la variable de respuesta, lo que no ocurre en el ML, en el cual los valores estimados pueden tener valores para 𝜇𝑖<0.
Este artículo tiene como propósito describir y aplicar el Modelo Lineal Generalizado (MLG) para el análisis de datos de un ensayo que se realizó en el Campo Experimental El Palmar de la Vizcaína (CEPV) de Cenipalma, dirigido a validar el control de Leptopharsa gibbicarina (Froeschner, 1976) (Hemiptera: Tingidae) a través del uso del hongo entomopatógeno Purpureocillium lilacinum (Thom.).
Materiales y métodos
El MLG, propuesto por Nelder y Wedderburn (1972) es una extensión del ML, en el cual la variable de respuesta puede tener una distribución normal, binomial, Poisson o binomial negativa entre otras. El MLG está conformado por tres componentes: aleatoria, sistemática y una función de enlace; la componente aleatoria se refiere a los valores observados de la variable de respuesta con función de distribución f (y) la cual debe pertenecer a la familia de funciones exponenciales. La componente sistemática θi = 𝑥𝑖𝛽 corresponde a la variabilidad de la variable de respuesta Y que es explicada por X a partir de una combinación lineal 𝑥𝑖𝛽, conocida como el predictor o conjunto de parámetros y variables explicativas; la función de enlace relaciona las componentes aleatoria y sistemática, y se expresa como 𝑔(𝜇𝑖). La importancia de esta función es que permite transformar los valores esperados a la misma escala de los del predictor, es decir 𝑔(𝜇𝑖)= 𝑥𝑖𝛽; por ejemplo, en el caso particular de la distribución de Poisson, se espera que 𝜇𝑖 tome valores 1, 2, 3,..., al igual que los estimados, sin embargo estos últimos pueden tomar valores negativos cuando no se considera la función de enlace como en el caso del ML.
Si la variable de respuesta de interés es cuantitativa y discreta y se registra como el conteo de un evento en un espacio o tiempo dados, puede hacerse uso de la distribución Poisson (Agresti 2002). Esta distribución se expresa de la forma
El uso del MLG sigue un conjunto de pasos al igual que el ML para el análisis de los datos: i) identificación y selección del modelo: se evalúan los supuestos de la componente aleatoria, se establece la componente sistemática de acuerdo con la forma como se observaron o registraron los datos y se determina la función de enlace; se debe considerar la significancia o no de cada uno de los términos del modelo y la interpretación de los mismos (Lindsey 1995). Los modelos propuestos deben evaluarse en términos de la magnitud de la discrepancia entre los datos observados y los valores esperados, la cual se formula acorde a la distribución de la componente aleatoria. ii) Estimación de parámetros: se hace a partir de métodos iterativos para ecuaciones no lineales como las generadas por este tipo de modelos (Nelder y Wedderburn 1972), a partir de la resolución de estas ecuaciones se determina el máximo valor (Dobson y Barnett 2018); iii) Evaluación del modelo: esta se centra en la propuesta de las componentes del MLG: distribución de la componente aleatoria, especificación de las variables explicativas o componente sistemática y selección de la función de enlace asociada a la componente aleatoria (Agresti 2002); una especificación incorrecta de la componente sistemática genera errores que pueden generar violación de supuestos en la distribución de la variable de respuesta, por ejemplo, en el caso de la distribución de Poisson la relación de igualdad entre la media y la varianza, conocida como equidispersión (Agresti 2002; Morales y López 2009; Torres Blanco 2011). Cuando este supuesto no se cumple, la media de la variable de respuesta puede estar cambiando y esta variación debe ser modelada. Si esto ocurre y, la variable de respuesta de interés Y ~ P (λ), entonces el parámetro λ sigue una distribución Gamma, es decir 𝜆~𝐺 (𝑘‚ 𝜇), y su expresión es de la forma
esta distribución tiene 𝐸(𝜆) = 𝜇 ?? 𝑉(𝜆) = 𝜇2/𝑘 (Agresti 2002); la mezcla de estas dos distribuciones (Poisson y Gamma) genera la distribución Binomial Negativa (BN), que es la de interés, dado que permite modelar conteos en casos de sobredispersión (Agresti 2002). La distribución BN se expresa de la siguiente forma:
su esperanza y varianza están dadas por
Para mostrar el uso del MLG, se utilizaron datos generados en un ensayo, llevado a cabo en Cenipalma, para validar el control del hemíptero Leptopharsa gibbicarina, responsable de ocasionar lesiones en la hoja de la palma de aceite, al facilitar el desarrollo del hongo Pestalotiopsis palmarum. Esta especie generanecrosis, destruye la lámina foliar y disminuye así la capacidad fotosintética de la planta. Esta enfermedad se denomina Pestalotiopsis y afecta la producción. En Cenipalma se encontró que la aplicación del hongo entomopatógeno Purpureocillium lilacinum (CPPl0601) en dosis de 1x1013 conidias/ha, controla el estado adulto de L. gibbicarina (Barrios et al. 2016).
Con el propósito de aplicar el MLG, se utilizaron los datos registrados en un ensayo en el Campo Experimental El Palmar de La Vizcaína, realizado en un lote de 3,9 Has, ubicado en la vereda Peroles, de Barrancabermeja, Santander, Colombia (6º59’32”N 73º43’04”O), con presencia de pestalotiopsis y de L. gibbicarina. El lote se dividió en dos áreas de similar tamaño y se asignó aleatoriamente a una de las áreas la aplicación de P. lilacinum y a la otra no se le hizo aplicación, siendo el testigo. La aplicación del entomopatógeno se realizó en abril de 2016 en el híbrido OxG cultivar Coarí x LaMé, año de siembra 2012, de acuerdo con los criterios y especificaciones indicados en Barrios et al. (2016). En cada una de las áreas se hizo el conteo inicial de adultos de L. gibbicarina en la hoja 25 (según filotaxia de la palma) utilizando muestreo sistemático descrito por Särndal et al. (2003); se seleccionó aleatoriamente una palma como punto de partida y, seguidamente, cada dos palmas se realizó el conteo de insectos hasta recorrer todas las palmas de cada área; además se aplicó el entomopatógeno en el área seleccionada. A los 30 días después de la aplicación (30 dda) del hongo, se realizó el conteo de insectos para las dos áreas, en la misma palma y hoja evaluada al inicio. Los datos se analizaron con el uso del software estadístico SAS Institute Inc (2008).
Resultados
Descripción de los datos. En la Tabla 1 se registran estadísticas descriptivas para el tratamiento control y para el tratamiento donde se aplicó Purpureocillium lilacinum; el conteo de L. gibbicarina para los dos tratamientos antes de aplicar el hongo presenta diferencias considerables en sus promedios y mayor variabilidad en el tratamiento de validación, situación que se presentó por la dispersión del insecto y que se encontró después de aleatorizar el área de aplicación de los tratamientos. La evaluación a los 30 días después de aplicar el hongo presenta menores diferencias en términos de promedio.
Control | Validación | |||||
---|---|---|---|---|---|---|
Variable | n | Media | Varianza | N | Media | Varianza |
dda0 | 13 | 25,23 | 745,69 | 14 | 89,21 | 6460,64 |
dda30 | 13 | 69,15 | 11139,42 | 14 | 87,78 | 4044,64 |
dda: días después de la aplicación.
Análisis de los datos utilizando el MLG. Para el análisis se utilizó la distribución de Poisson con un factor de estudio (tratamientos) y se consideró la lectura antes de aplicar los tratamientos como una covariable (Arango Londoño et al. 2009). Se utilizó como función de enlace 𝑔 (𝜇) =𝑙𝑜𝑔 (𝜇𝑖). En la Tabla 2 se muestran resultados del análisis que permiten verificar la bondad de ajuste del modelo propuesto.
Criterios | GL | Valor | Valor/GL |
---|---|---|---|
Desviación | 50 | 3236,85 | 64,74 |
Desviación escalada | 50 | 3236,85 | 64,74 |
χ2 de Pearson | 50 | 3815,21 | 76,30 |
Log de verosimilitud | 12296,67 | ||
Log de verosimilitud total | - 1761,72 | ||
AIC (más pequeño mejor) | 3531,43 | ||
AICC (más pequeño mejor) | 3532,25 | ||
BIC (más pequeño mejor) | 3539,39 |
A partir de la desviación se evalúa primero la equidispersión del modelo, es decir, si la media y la varianza del modelo son iguales como debe ocurrir en el modelo de Poisson. El valor de la desviación para este caso fue de 64,41 que es mucho mayor de 1,0 indicando que esta afirmación no se cumple (McCullagh y Nelder 1989), es decir el modelo presenta sobredispersión, la varianza es mayor que la media para los datos analizados, situación que se presenta con cierta frecuencia en la distribución de Poisson (Hinde y Demétrio 1998; Agresti 2002).
Una alternativa para el análisis de los datos considerada en este trabajo fue analizarlos a través de la regresión binomial negativa que asume la varianza como una función de la media (Navarro et al. 2001; Morales y López 2009). La bondad de ajuste de este modelo se muestra en la Tabla 3.
Criterios | GL | Valor | Valor/GL |
---|---|---|---|
Desviación | 50 | 62,65 | 1,25 |
Desviación escalada | 50 | 62,65 | 1,25 |
χ2 de Pearson | 50 | 53,49 | 1,07 |
Log de verosimilitud | 13781,34 | ||
Log de verosimilitud total | - 277,05 | ||
AIC (más pequeño mejor) | 564,09 | ||
AICC (más pequeño mejor) | 565,34 | ||
BIC (más pequeño mejor) | 574,04 |
Los resultados muestran un efecto importante de los tratamientos con P = 0,0019, indicando diferencias significativas entre los dos tratamientos; no se observa efecto de la covariable con P = 0,9674 (Tabla 4). La estimación del parámetro de dispersión es de 1,084 con un intervalo de confianza del 95 % de (0,764; 1,537); este intervalo no contiene el cero, luego el parámetro > 0, indicando que el modelo Binomial Negativo es más apropiado que el de Poisson para modelar los datos. El estimador para el tratamiento 1 (aplicación del hongo) fue de - 1,263 e indica la diferencia esperada entre este tratamiento y el tratamiento control, así 𝑒−1,263= 0‚2828 indica que la tasa de incidencia del tratamiento 1 (aplicación del hongo) es 0,2828 veces mayor que el testigo.
Estimación de parámetros de máxima verosimilitud | ||||||||
---|---|---|---|---|---|---|---|---|
Parámetro | GL | Estimador | Error estándar | Wald 95 % Límites de confianza | Wald χ2 | Pr > χ2 | ||
Intercepto | 1 | 4,491 | 0,279 | 3,943 | 5,039 | 257,85 | < 0001 | |
Covariable | 1 | - 0,016 | 0,395 | - 0,791 | 0,759 | 0,00 | 0,9674 | |
Tratamiento | 1 | 1 | - 1,263 | 0,406 | - 2,058 | - 0,468 | 9,69 | 0,0019 |
Covariab*trat | 1 | 1 | 1,031 | 0,572 | - 0,090 | 2,152 | 3,25 | 0,0715 |
Dispersión | 1 | 1,084 | 0,193 | 0,764 | 1,537 |
Discusión
En el análisis de regresión de Poisson antes de interpretar y concluir acerca de la hipótesis de interés (efecto de la aplicación el hongo entomopatógeno), es necesario verificar si el modelo propuesto es el indicado para la modelación de los datos; uno de los supuestos para este modelo es la equidispersión, es decir la media y la varianza deben ser iguales; la prueba muestra que este supuesto no se cumple indicando sobredispersión (Var (Y) > E (Y), Tabla 2); este resultado se esperaba dada la descripción del recuento de insectos antes y después de la aplicación del tratamiento la cual muestra relaciones varianza-promedio distantes de ser iguales para las dos mediciones (Tabla 1). La sobredispersión puede presentarse por alta variabilidad entre las unidades observadas y no considerar una estrategia de agrupación de estas con condiciones de variación similares, aunque esta situación se presenta con cierta frecuencia en el caso de la distribución de Poisson (Morales y López 2009). Una alternativa para solucionar esta situación es hacer uso de la distribución binomial negativa, la cual asume que la varianza presenta un valor superior a la media y se utilizó en este trabajo para completar el análisis de los datos utilizando MLG.
La Tabla 3 muestra el análisis de los datos utilizando la regresión binomial negativa; el estadístico chi cuadrado de Pearson no evidencia rechazar la hipótesis del modelo binomial negativo para la modelación de los datos, lo que permite su uso para el análisis estos; esta distribución ha sido utilizada en experimentos que impliquen conteo de insectos (Silva et al. 2016). La comparación de los dos tratamientos a través del estadístico χ2 (Tabla 4), muestra diferencia entre los dos tratamientos y el coeficiente negativo para el tratamiento aplicación de P. lilacinum, indica que la aplicación de este hongo entomopatógeno genera una reducción en la población de L. gibbicarina, resultado acorde con lo encontrado por Barrios et al. (2016), quienes evaluaron la misma cepa de P. lilacinum en ecosistemas palmeros para el control de L. gibbicarina.
No se observó efecto de la covariable (con un α = 0,05) para las condiciones del lote donde se realizó la evaluación; estos aspectos admiten una recomendación confiable del hongo a nivel comercial. Sin embargo, debido a las diferentes condiciones climáticas presentes dentro de las zonas palmeras en Colombia y teniendo en cuenta la influencia que tienen factores climáticos como la temperatura, humedad relativa y radicación UV en la eficacia de los hongos entomopatógenos (Pucheta Díaz et al. 2006; Carrillo-Rayas y Blanco-Labra 2009), se sugiere realizar otros ensayos de validación bajo diferentes condiciones ambientales. Además, se recomienda: i) verificar antes de iniciar el ensayo la distribución espacial del insecto; ii) contar con repeticiones para cada uno de los tratamientos, que permita hacer inferencia acerca de la bondad del hongo para el control del insecto y finalmente iii) evaluar la posible migración de insectos de un lote al otro dada la cercanía de los mismos en el ensayo realizado.