Datos experimentales

Para aplicar e ilustrar la técnica de suavizado enfocado en Psplines, se utilizaron datos reales proporcionados por ^{Lugo (2018)}, correspondientes a su tesis doctoral titulada “Efecto de Bacterias Antagonistas y Extractos Vegetales sobre la Dinámica Poblacional de Ralstonia Solanacearum (Smith) Yabuuchi et al. y la Expresión de Síntomas en Tomate (Solanum Lycopersicum l.)”

Para la obtención de R. solanacearum se planificó un experimento en un diseño completamente aleatorizado, con cinco tratamientos y tres réplicas. La unidad experimental consistió de una maceta de 1500 g de capacidad (una planta por maceta). Los tratamientos se describen en la tabla 1.

La medición del tamaño de la población de R. solanacearum en el suelo, medido en UFC/g de suelo se hizo cada semana, es decir, 20 mediciones de la variable, las cuales se muestran en la tabla 2.

Dinámica poblacional

Para obtener los modelos matemáticos que describen la dinámica poblacional en el suelo de R. solanacearum bajo las condiciones descritas, el autor utilizó las mediciones de UFC/g de suelo obtenidos de cada tratamiento en la fase experimental, graficó y trató de ajustar a uno o más de los modelos matemáticos tomando en cuenta los basamentos teóricos de los modelos, el ajuste gráfico a la tendencia de los datos y la obtención de estimadores iniciales de los parámetros. Trabajó con 32 series de datos confeccionadas en las condiciones de las plantas definidas como total, sanas y enfermas y a su vez, sobre estas condiciones se definieron series de diferentes longitudes: serie completa y en dos fases, creciente y decreciente.

En total se seleccionaron 82 modelos usando el software R-Project y se eligieron los cinco mejores modelos matemáticos considerando tendencias gráficas, pruebas de t y F realizadas por el software R-Project, coeficientes basados en desviación entre estimados y observados y criterios de información. Las series seleccionadas para este trabajo pueden verse mejor en la tabla 3 junto con los nombres de los modelos para cada serie y los valores de los criterios de información de estos tomados de ^{Bandera y Pérez (2018)}.

Análisis de datos

Para el análisis de los datos se hace uso del software R, específicamente la función gam() para lo cual es necesario cargar el paquete mgcv de ^{Wood (2006}). Esta función puede usar p-splines univariantes según lo propuesto por ^{Eilers y Marx (1996}). En realidad, este paquete contiene dos funciones que permiten utilizar P-splines: gam y gamm, la diferencia entre las dos es que la segunda permite elegir el parámetro de suavizado mediante REML, mientras que la primera es similar a la función escrita por ^{Hastie y Tibshirani (1990}), pero permite utilizar splines de rango bajo, además de haber corregido los errores que existían en el cálculo de la varianza de los parámetros y elige el parámetro de suavizado mediante la validación cruzada generalizada (GCV por sus siglas en inglés). En ambos casos se puede imponer un valor arbitrario para el parámetro de suavizado sin que sea elegido por la propia función, elegir la base a utilizar, según las ofrecidas por el paquete y se puede elegir el número de nodos y el orden de la penalización. Los argumentos de la función GAM se muestran en la tabla 4.

El procedimiento para analizar los datos en R se hace de forma similar para las 6 series. Inicialmente se especifica un modelo P spline a través de la función GAM con un término de suavizado de la forma 𝑠 (𝑥, 𝑏𝑠 = "𝑝𝑠", 𝑘=𝑣𝑎??𝑖𝑎𝑏𝑙𝑒, 𝑠𝑝=𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒), el cual especifica una base P-spline con 𝑚=2 por defecto, que significa una base y una penalización de segundo orden. “Las diferencias entre los modelos dependerán de los argumentos variables como el tamaño de base (??) y el parámetro de suavizado (𝑠𝑝), esto con el fin de seleccionar el mejor ajuste p splines para los datos” (^{Burbano et al., 2022}, p. 270).

Es por ello que primero, se especifica un modelo para la serie en cuestión variando el tamaño de base 𝑘 y se hace una selección de los tres mejores modelos, dependiendo de los criterios de información, métodos de selección del parámetro de suavizado e índices, los cuales deben ser mínimos o máximos según sea el caso. Después de obtener el mejor ajuste P spline, mediante un gráfico de dispersión, se verifica el efecto del suavizado comparando este ajuste versus el modelo proporcionado por ^{Lugo (2018)} como mejor para la descripción de la dinámica poblacional de R. solanacearum en la serie especificada y se muestra la fórmula del modelo aproximado para el mejor o los mejores ajuste P spline. Luego se cotejan las curvas ajustadas para los distintos valores de k.

Segundo, al obtener el mejor ajuste en la parte anterior, se procede a probar el mismo bajo diferentes valores de 𝜆, con el fin de obtener un valor de este parámetro que optimice el ajuste. Para realizar la variación de λ, se toma como mínimo y máximo, aquel valor donde el AIC tienda a aumentar. Para esto, se consideran los mismos ítems que en la parte anterior. Seguidamente, sobre este ajuste se contrasta gráficamente el efecto de distintos valores del parámetro de suavizado. Los ítems que resumen la información arrojada por los tres mejores modelos seleccionados para escoger el modelo final se describen en la tabla 5.

Resultados de serie TOTAL

Se presenta el modelo planteado por ^{Lugo (2018}) para describir la dinámica poblacional de Ralstonia solanacearum en el suelo para la serie TOTAL:

En la tabla 6, se resumen los resultados obtenidos donde se puede observar que en los tres ajustes P splines, el término de suavizado (𝑠(𝑡)), resulta significativo. El modelo fitT.7 arrojó el menor GCV y BIC, sin embargo, el ajuste denotado por el modelo fitT.9, presenta como favorables, tres de los cinco criterios considerados para la selección del mejor, a saber, el menor AIC y el mayor logLik y R²adj, por lo cual se considera como el modelo con mejor ajuste y el modelo aproximado es el siguiente:

fitT.9<-gam (y ~ s(t, bs="ps",k=9))

Se representa las variantes del ajuste P spline. Note que los ajustes no son tan distintos. En azul, el mejor ajuste. No obstante, al comparar con modelo total, se evidencia que este presenta un menor AIC y un menor BIC en comparación con el mejor ajuste P spline (Figura 1).

Elaboración: El autor.

Figura 1. Representación gráfica del efecto del suavizado con IC del 95% y 3 valores de k para la serie TOTAL. 

En la figura 2, se contrasta el mejor ajuste P spline seleccionado versus modelo total y puede verse que los modelos no difieren mucho, sin embargo, se nota como el ajuste P spline genera una ondulación en la segunda mitad de los datos.

Elaboración: El autor.

Figura 2. Mejor ajuste P spline vs. modelo total. 

Ahora bien, en el modelo seleccionado como mejor, fitT.9, se prueban cinco valores distintos para λ entre 0,01 y 10, para un valor fuera de este rango el valor del AIC tiende al aumento. En la tabla 7 se resume la información arrojada. Tome en cuenta que el modelo fitT.9c contiene un valor de 𝜆 estimado por la función GAM, este tiende a ser pequeño, lo que indica que se está frente a un ajuste P splines.

Al hacer las comparaciones, se encuentra que los modelos fitT.9a y fitT.9c presentan igual número de criterios favorables para optimizar el ajuste. Los modelos aproximados son:

fitT.9a<-gam (y ~ s(t, bs="ps",k=9,sp=0,01)) y fitT.9c<-gam(y ~ (t,bs="ps",k=9,sp=0,36))

Cuando la penalidad es más débil se obtiene una curva más ondulada, caso contrario ocurre para un mayor 𝜆 que produce curvas más suaves. Esto puede verse mejor en la figura 3. A la izquierda y en azul, el menor valor de 𝜆 genera una curva con más picos, en el centro, el modelo con 𝜆 estimado por la función GAM y a la derecha un 𝜆 mayor.

En la figura 4 se evidencia de mejor forma como el valor de 𝜆 controla la suavidad de la curva. Note la ondulación de la curva cuando 𝜆=0.01 y la tendencia lineal que toma cuando 𝜆=10.

Elaboración: El autor

Figura 3. Mejor ajuste P spline para 3 valores distintos de sp en serie TOTAL. 

Elaboración: El autor.

Figura 4. Efecto de sp sobre el mejor ajuste P spline en serie TOTAL. 

Resultados de Serie T1

El modelo planteado por ^{Lugo (2018)} para la serie T1, es el siguiente:

La tabla 8, contiene los resultados para esta serie, donde se evidencia que los tres modelos presentan un suavizado, s(t), significativo, no obstante, el modelo fitT1.19, arrojó los menores valores para CGV, AIC, BIC y los mayores valores para logLik y R²adj en comparación con los modelos restantes y modeloT1, por lo que es considerado como el mejor modelo para el ajuste de esta serie y el modelo aproximado es el siguiente:

FitT1.19<-gam (y ~ s(t,bs="ps",k=19))

Visualice en la figura 5 como cambia la curva para distintos valores de k. Es fácil notar en el centro y en azul como el ajuste presenta una ondulación mayor en la curva en comparación con los demás que producen curvas más suaves.

Elaboración: El autor

Figura 5. Representación gráfica del efecto del suavizado con IC del 95% y 3 valores de k para la serie T1. 

En la figura 6, se puede apreciar la diferencia entre las curvas que producen el mejor ajuste P splines y modeloT1.

Elaboración: El autor

Figura 6. Mejor ajuste P spline vs. modeloT1. 

En la tabla 9 se encuentran los resultados para el control del suavizado del ajuste, con un rango de variación de 𝜆 entre 0,00001 y 10. El modelo fitT1.19a presenta los menores valores de los criterios de selección, sin embargo, este suavizado no es significativo. El modelo fitT1.19c arrojó el menor GCV, pero el modelo fitT1.19b generó el menor AIC y BIC y los mayores valores de logLik y R²adj, para el mejor ajuste de la serie T1, aunque la incertidumbre aumenta en los extremos. El modelo aproximado anterior es:

fitT.19b<-gam (y ~ s(t,bs="ps",k=19,sp=0,00001))

Elaboración: El autor

Figura 7. Mejor ajuste P spline para 3 valores distintos de sp en serie T1. 

La figura 8, permite comparar las curvas generadas por los ajustes P splines con los valores de 𝜆 fijados arbitrariamente y el estimado por la función GAM.

Elaboración: El autor

Figura 8. Efecto de sp sobre el mejor ajuste P spline en serie T1. 

Resultados de Serie T2

El modelo sugerido por ^{Lugo (2018)} para la serie T2, es el siguiente:

El modelo anterior, no describe de forma correcta el patrón de datos correspondiente, lo cual puede deberse a un error de transcripción. Por lo anterior, el modelo aproximado se muestra a continuación:

La tabla 10, recoge los resultados arrojados por los modelos P splines planteados para esta serie. Se puede ver que el modeloT1 presenta el menor BIC mientras que el modelo fitT2.8 generó el menor GCV y AIC y el mayor logLik y R²adj, así se considera a este modelo como la mejor aproximación para la serie T2 y su modelo aproximado es el siguiente: fitT2.8<-gam (y ~ s(t, bs="ps",k=8))

La figura 9, permite contrastar los ajustes P splines para distintos valores de k. En el centro y en azul la mejor aproximación. Note la ondulación que presenta la curva en el extremo derecho y que los dos ajustes restantes no generan.

Elaboración: El autor

Figura 9. Representación gráfica del efecto del suavizado con IC del 95% y 3 valores de k para la serie T2. 

Para verificar la significancia del suavizado P splines en los datos, observe la figura 10, donde se contrastan las curvas generadas por el mejor ajuste P spline y modeloT2.

Elaboración: El autor

Figura 10. Mejor ajuste P spline vs. modeloT2. 

Para obtener el valor de 𝜆 que optimice el ajuste P spline, se varia en un rango entre 0.001 y 10. Sin embargo, esta variación no causo un efecto significativo, ya que los criterios de selección favorables son proporcionados por el modelo con el valor de 𝜆 estimado por la función GAM. El modelo fitT2.8a arrojó el mayor loglik, mientras que el modelo fitT2.8c, con 𝜆 estimado, presenta el menor GCV, AIC, BIC y el mayor R²adj. Estos resultados se resumen en la tabla 11. El modelo aproximado para este ajuste es:

fitT2.8c<-gam (y ~ s (t, bs="ps",k=8,sp=0,178))

La figura 11, permite comparar los ajustes generados por las variaciones en 𝜆. Las curvas ajustadas no difieren una de la otra en mayor proporción, excepto la curva con un 𝜆=10 que tiende a la linealidad. Aprecie esto de mejor forma en la figura 12.

Elaboración: El autor

Figura 11. Mejor ajuste P spline para 3 valores distintos de sp en serie T2. 

Elaboración: El autor

Figura 12. Efecto de sp sobre el mejor ajuste P spline en serie T2. 

Resultados de Serie T3

El modelo a partir de los datos sugerido por ^{Lugo (2018)} para la serie T3 es:

En la tabla 12, se plasman los resultados obtenidos para los ajustes P splines. El modeloT3 obtuvo el menor BIC mientras que el modelo fitT3.9 con arrojó el menor GCV y AIC y el mayor logLik y R²adj. Así, el mejor modelo para ajustar esta serie queda denotado por fitT3.9. el cual tiene la siguiente forma: fitT3.9<-gam (y ~ s(t,bs="ps",k=9))

La figura 13, muestra el efecto del suavizado en los datos para los distintos valores de k utilizados para la serie T3. A la izquierda y en azul, el modelo seleccionado como mejor ajuste P spline. Presenta una mayor ondulación. Se aprecia, el efecto del suavizado en las curvas.

Elaboración: El autor

Figura 13. Representación gráfica del efecto del suavizado con IC del 95% y 3 valores de k para la serie T3. 

En la figura 14, se pueden cotejar las curvas generadas por el mejor ajuste P splines en la serie T3 y modeloT3. El modelo P spline genera una curva menos lineal.

Elaboración: El autor

Figura 14. Mejor ajuste P spline vs. modeloT3. 

La tabla 13, resume los resultados arrojados para la selección del valor óptimo de λ sobre el mejor ajuste P spline para la serie T3. El modelo fitT3.9a generó el mayor logLik y R2adj, este último valor compartido con el modelo fitT3.9b que a su vez obtuvo el menor AIC. Así mismo, el modelo fitT3.9c arrojó el menor GCV y BIC. Cada uno de los modelos tienen dos criterios de selección a su favor, sin embargo, se considera más relevante el valor del AIC, así se toma como mejor ajuste al modelo denotado por fitT3.9b y su modelo aproximado es: fitT3.9b<-gam (y ~ s(t,bs="ps",k=9,sp=0,0029))

En la figura 15, pueden compararse estos ajuste y a su vez notarse como se comportan las curvas ajustadas cuando 𝜆 tiende a aumentar, el ajuste se alisa más.

Elaboración: El autor

Figura 15. Mejor ajuste P spline para 3 valores distintos de sp en serie T3 

No se aprecia mucha diferencia entre las curvas ajustadas de los modelos fitT3.9b y fitT3.9c. Observe la figura 16 para identificar que el ajuste en las dos curvas es casi idéntico.

Elaboración: El autor

Figura 16. Efecto de sp sobre el mejor ajuste P spline en serie T3.