INTRODUCCIÓN
En el manejo forestal sustentable es necesario el conocimiento ágil y preciso de los ecosistemas; el avance tecnológico, la disponibilidad de programas estadísticos y técnicas biométricas potentes, permite obtener información certera e implementar herramientas de apoyo para la toma de decisiones (Santiago-García et al. 2020a). La estimación del crecimiento y rendimiento de las masas forestales mediante la modelación matemática, constituye una base científica que proporciona información imprescindible para el manejo forestal sustentable (Salas y Real 2013, Hernández-Cuevas et al. 2018).
La estimación del crecimiento de las masas forestales requiere evaluar el potencial de producción del espacio físico, conocido como calidad de sitio, siendo para este procedimiento, el método de índice de sitio el más utilizado, en el que la altura dominante es el indicador de la productividad forestal por excelencia (Clutter et al. 1983, Scolforo et al. 2013). El índice de sitio se define como la altura dominante que puede alcanzar el rodal a una edad determinada, reflejo de los factores bióticos y abióticos del área forestal (Gadow et al. 2012, Özçelik et al. 2018, Seppänen y Mäkinen 2020, Socha et al. 2021). De acuerdo con Carrero et al. (2008) y García et al. (2021), este método permite el análisis gráfico mediante el desarrollo de un conjunto de curvas que expresan el crecimiento promedio en altura dominante, proyectadas mediante ecuaciones no lineales, lo que permite calificar la productividad del sitio forestal. Las curvas de índice de sitio poseen ciertas propiedades, entre las que destacan: crecimiento sigmoidal con un punto de inflexión, una asíntota horizontal a edades avanzadas, comportamiento lógico, son invariantes con respecto a la edad de referencia y al camino de simulación (Bailey y Clutter 1974, Cieszewski y Bailey 2000).
En el campo forestal existen diversos estudios del ajuste de modelos de altura dominante e índice de sitio con el fin de aproximar de forma numérica la productividad maderable del área forestal (Vargas-Larreta et al. 2013, Scolforo et al. 2013, Nava-Nava et al. 2020), no obstante, previo a la aplicación de estos, es conveniente comprobar su validez, debido a que un modelo es útil siempre y cuando sea realista, simplifique la realidad y apoye a los intereses del silvicultor (Snee 1977). El proceso de validación de modelos, debe ser objetivo, porque podría suceder que el modelo elegido no sea el adecuado. Por lo anterior, Tedeschi (2006) sugiere tomar en cuenta tres aspectos para validar un modelo, el primero, todos los modelos presentan errores, pero unos describen o calculan mejor que otros, segundo, no aficionarse con un modelo y la exclusión de otros, tercero, analizar y acreditar el ajuste del modelo a los datos, y con esto, generar confianza en el modelo actual o permitir la selección de modelos alternativos.
Snee (1977) establece que existen 4 procedimientos para comprobar la validez de modelos de regresión: a) comprobación de las predicciones del modelo y coeficientes, b) recolección de datos nuevos, para su cotejo con las predicciones del modelo, c) comparación de resultados con modelos teóricos y datos simulados, d) partición de los datos disponibles (validación cruzada), dirigiendo una parte al ajuste y la otra a la validación del modelo. De acuerdo con Arlot y Celisse (2010), la validación tiene como objetivo seleccionar el modelo más apropiado para predecir, con parámetros adecuados y, cuando no es posible, recopilar datos independientes del ajuste. Esta técnica es aplicada por su eficiencia y depende principalmente de la partición de datos porque esto influye en los estadísticos de bondad de ajuste (Zhang y Yang 2015).
En México, Ixtlán de Juárez, Oaxaca, es una comunidad modelo del manejo forestal comunitario, con más de 30 años de experiencia, donde la principal actividad económica es el aprovechamiento maderable de los rodales con diversas especies, principalmente de pino (Santiago-García et al. 2022). Por tanto, para la toma de decisiones de manejo forestal, se requiere de herramientas silvícolas cuantitativas que permitan determinar el potencial productivo de sus bosques. Entre los estudios que se han realizado en el área se puede mencionar a Hernández-Cuevas et al. (2018), quienes ajustaron modelos de altura dominante e índice de sitio para Pinus ayacahuite Ehrenb. ex Schltdl., con datos de análisis troncal. Por su parte, Castillo-López et al. (2018) desarrollaron ecuaciones dinámicas a nivel regional para Pinus oaxacana Mirov, P. patula Schiede ex Schltdl. & Cham., P. douglasiana Martínez y P. pseudostrobus Lindl., mediante la misma técnica de análisis, en tanto que, Nava-Nava et al. (2020) efectuaron ajustes de ecuaciones dinámicas para Pinus patula de Ixtlán de Juárez, con una base principalmente de parcelas permanentes. Pinus oaxacana es una de las especies con mayor importancia comercial en los boques comunales de Ixtlán de Juárez, por tanto, es necesario conocer la productividad maderable de los rodales. El objetivo del estudio fue ajustar y validar modelos de crecimiento en altura dominante e índices de sitio. La hipótesis fue que modelos de polimorfismo asintótico al presentar mayor flexibilidad, generan curvas que proyectan adecuadamente el crecimiento y productividad de P. oaxacana, dada la variabilidad natural de la especie.
MÉTODOS
Área de estudio
El estudio se realizó en el bosque natural bajo manejo de Ixtlán de Juárez (figura 1), Oaxaca, México, localizado entre las coordenadas geográficas 17° 18’ 16” - 17° 30’ 00” N y 96° 21’ 29” - 96° 31’ 38” O, en un intervalo de altitudes que varían entre 2.200 y 3.100 m s.n.m. El área que ocupa el bosque de pino-encino presenta un clima templado húmedo con lluvias en verano, siendo este tipo de vegetación el que ocupa el primer lugar en extensión, donde se realiza principalmente aprovechamiento maderable de Pinus patula, P. ayacahuite, P. douglasiana, P. pseudostrobus, P. teocote Schltdl. & Cham., P. leiophylla Schiede & Deppe, P. rudis Endl., y P. oaxacana (Santiago-García et al. 2020b).
Datos dasométricos
En el análisis se utilizaron 44 parcelas permanentes de investigación silvícola de 400 m2 cada una, con cuatro inventarios anuales (2015-2018), así como de 58 parcelas temporales establecidas en 2005 y proyectadas al año 2006 con un modelo anamórfico de altura dominante para P. oaxacana reportado en el programa de manejo forestal (ciclo de corta 2015-2024) de Ixtlán de Juárez, Oaxaca (STF 2015), así, se crearon 58 pares de datos adicionales, de modo que la trayectoria de dichas parcelas hizo posible integrar la información para el ajuste de los modelos dinámicos; lo cual ayudó a enriquecer la base de datos. Lo anterior, se realizó de acuerdo con Bates y Granger (1969), Morris (1974) y Bunn (1989), porque para complementar la base de datos se usaron pronósticos del modelo previo, de modo que se empleó toda la evidencia disponible y fuentes de información sólidas.
Las parcelas están localizadas en rodales puros y coetáneos de P. oaxacana. En cada parcela fue medido el diámetro normal (Dn, cm) de todos los individuos y la altura total (H, m) de al menos 8 árboles por sitio, de los cuales 4 eran dominantes, estos últimos describen exclusivamente a los árboles que tienen éxito en la lucha por la luz y por tal razón, representan el crecimiento máximo en altura de los rodales (Bengoa 1999). Con el fin de captar la mayor parte de las condiciones en las que se desarrolla y domina P. oaxacana, las parcelas permanentes y temporales presentaron diferentes densidades y calidades de sitio en las que se encontraba asociado con especies de pinos y latifoliadas, así, se obtuvo una distribución de datos con un intervalo amplio de edades y alturas dominantes (cuadro 1).
Modelos ajustados
Para la modelación del crecimiento en altura dominante de P. oaxacana se ajustaron las ecuaciones no lineales de Korf, Hossfeld IV y Chapman-Richards en sus expresiones ADA y GADA (por sus siglas en inglés) (Krumland y Eng 2005, Santiago-García et al. 2015, Hernández-Cuevas et al. 2018, Castillo-López et al. 2018). De acuerdo con Bailey y Clutter (1974) la generación de ecuaciones ADA (cuadro 2) implica plantear el modelo base en las condiciones inicial y futura, de manera que de la condición inicial se despeja el parámetro específico del sitio y se sustituye la solución en la condición futura. Santiago-García et al. (2015), establecen que las ecuaciones tipo ADA se constituyen por una ecuación de predicción y una ecuación de proyección, y que estas últimas permiten calificar la calidad de sitio en los rodales, a partir de una edad base, altura dominante y la edad del rodal. Por su parte, la derivación de ecuaciones GADA (cuadro 3) implica seleccionar una ecuación base en la que se identifican los parámetros dependientes de la calidad de sitio. Los parámetros seleccionados se expresan como funciones de la calidad de sitio definida por la variable X0, la ecuación base bidimensional se expande a una ecuación explicita tridimensional, de donde la variable X0 se despeja a partir de condiciones iniciales en edad y altura, lo que permite obtener ecuaciones dinámicas que proveen curvas polimórficas con múltiples asíntotas (Cieszewski 2002, Krumland y Eng 2005, Quiñonez-Barraza et al. 2015).
AD1: es la altura dominante (m) en la condición inicial (E1, años), AD2: altura dominante (m) correspondiente a una condición futura (E2, años), ln: logaritmo natural, exp: función exponencial y bi: parámetros a estimar.
Ajuste y validación cruzada
En el ajuste y validación de las ecuaciones dinámicas ADA y GADA, se realizó la aleatorización de las parcelas con el paquete “caTools” de R-Project® (R Core Team 2017), la base de datos fue particionada en 4 conjuntos con proporciones del 50 %, 60 %, 70 % y 80 %. Estos conjuntos se utilizaron para el ajuste de las ecuaciones dinámicas mediante el método iterativo anidado (nested iterative procedure) (Vargas-Larreta et al. 2010, Castillo-López et al. 2018, Nava-Nava et al. 2020) bajo la técnica de mínimos cuadrados no lineales (NLS). Mientras que las proporciones de datos restantes del 50 %, 40 %, 30 % y 20 % fueron empleadas para evaluar estadísticamente la capacidad predictiva de las ecuaciones ajustadas. Los ajustes se implementaron con el procedimiento MODEL del paquete estadístico SAS/ETS ® 9.4 (SAS Institute Inc. 2017).
Indicadores estadísticos
Para evaluar la bondad de ajuste de las ecuaciones se emplearon los estadísticos siguientes (cuadro 4): suma de cuadrados del error (SCE) [1], la cual representa la variación que no se puede explicar por la ecuación; el cuadrado medio del error (CME) [2], que refleja la varianza del modelo; la raíz del cuadrado medio del error (RCME) [3], que estima el error promedio en términos de la variable dependiente; el coeficiente de determinación ajustado (R 2 adj) [4], que expresa el porcentaje de la variabilidad total explicada por el modelo, tomando en cuenta el número de parámetros (Hernández-Cuevas et al. 2018, Santiago-García et al. 2020b). En la validación de los modelos de altura dominante se utilizó SCE [1], el error medio cuadrático EMC [5] que analiza la precisión de las estimaciones del modelo, el Ē, el cual estima el error promedio absoluto y finalmente, el IE que expresa en términos relativos la eficacia del modelo para estimar determinada variable (Vanclay 1994). La homogeneidad de varianzas en la distribución de los residuales se corroboró mediante la prueba de White (White 1980), en tanto que, la autocorrelación de los residuales se evaluó a través del estadístico de Durbin-Watson (Castillo-López et al. 2018, Nava-Nava et al. 2020), además, estos análisis se complementaron con gráficos de distribución de los residuales.
La comparación de los estadísticos de bondad de ajuste, el comportamiento gráfico y principalmente los valores de SCE, EMC, Ē e IE obtendidos en la validación de las ecuaciones ensayadas con las diferentes proporciones de datos, ayudaron a seleccionar los modelos que mejor predicen el crecimiento en altura dominante. Conjuntamente con el análisis numérico fue necesario analizar el comportamiento gráfico de cada modelo, debido que en algunas ocasiones se sobreestima a edades tempranas y en otras se subestima a edades avanzadas o viceversa (Diéguez-Aranda et al. 2006).
RESULTADOS
Ajuste de ecuaciones ADA
Mediante el uso de los cuatro conjuntos de datos con diferentes proporciones (50 %, 60 %, 70 % y 80 %) se obtuvieron ajustes satisfactorios, debido a que arrojaron valores altos en R 2adj, al explicar más del 99,9 % de la varianza total observada en altura dominante, además de que el estadístico de Durbin-Watson osciló entre 1,87 - 2,16 (figura 2), lo que indica que no existe autocorrelación en la distribución de los residuales de los modelos, cumpliendo así con el supuesto de independencia. De acuerdo con la prueba de White, la hipótesis nula no es rechazada (P > 0,05) por lo que en todos los casos se cumple con el supuesto de homogeneidad de varianzas. Respecto a la estimación de los parámetros, estos fueron significativos (P < 0,05) y con valores bajos en el error estándar, con excepción del parámetro b 0 de los modelos KA-b 1 -50, KA-b 1 -80, KA-b2-70 y KA-b2-80 los cuales obtuvieron errores estándar mayores a 2 (cuadro 5).
Ajuste de ecuaciones GADA
Las ecuaciones GADA presentaron ajustes adecuados debido a que se obtuvieron valores altos en R 2 adj, al explicar más del 99,8 % de la varianza total observada en la altura dominante y resultados cercanos a 2 en el estadístico de Durbin-Watson (figura 3). Asimismo, de acuerdo con la prueba de White (P > 0,05) todos los modelos presentaron homogeneidad de varianzas en la distribución de los residuales. Respecto a los estimadores de los parámetros, estos fueron significativos (P < 0,05) y con errores estándar bajos (con excepción de los parámetros b 1 de las ecuaciones HG-50 y HG-70, con errores altos) (cuadro 6).
Validación de ecuaciones ADA y GADA
La bondad de ajuste y el comportamiento gráfico de las ecuaciones ensayadas fueron aspectos importantes al seleccionar el modelo mejor para estimar la altura dominante y el índice de sitio. Sin embargo, fue indispensable conocer el error de predicción, por lo cual se realizó la validación de cada modelo ajustado, mediante la comparación estadística de las alturas estimadas por los modelos ajustados ADA y GADA y las alturas observadas en campo (bases de datos con proporciones del 50 %, 40 %, 30 % y 20 %). En las figuras 4 y 5 se muestran los estadísticos de la validación (SCE, EMC, Ē e IE), los cuales permitieron evaluar la capacidad predictiva de los modelos ensayados.
Curvas de índice de sitio
Para la construcción de las curvas de crecimiento en altura dominante, se calculó la curva promedio de cada modelo, y las alturas dominantes (AD2) al emplear los parámetros estimados de cada modelo ajustado, la edad de proyección (E2), la edad de referencia de 40 años (E1) y el valor del índice de sitio (AD1), posteriormente se trazaron tres curvas equidistantes por arriba y debajo de la curva promedio. En la figura 6 se presentan las curvas de índice de sitio generadas con los modelos mejores de cada metodología.
Graficas de residuales
Los valores de los residuales de los modelos seleccionados se encontraron en un intervalo de -0,9 a 1,8, además presentaron un patrón aleatorio por abajo y arriba del valor cero de referencia y no existió una tendencia definida, lo que corrobora la homogeneidad de varianzas de los errores (figura 7).
En el análisis de los residuales, para detectar problemas de autocorrelación, además del estadístico de Durbin-Watson, se graficaron los residuales frente a los residuales retardados, una (Lag1), dos (Lag2) y tres (Lag3) observaciones anteriores (figura 8). De esta manera se corroboró que los residuales se distribuyen de forma independiente. Los gráficos correspondientes al primer retardo (figura 8A1, 8B1, 8C1 y 8D1) en todos los modelos, parecen representar mejor el patrón aleatorio de los residuales alrededor de las líneas cero de referencia, en tanto que, con 2 y 3 retardos exhiben una tendencia ligera, pero sin problemas de autocorrelación.
DISCUSIÓN
En esta investigación, al usar de forma pragmática la combinación tanto de dos fuentes de información (parcelas permanentes y temporales), como de los pronósticos generados con el modelo previo, y el método de ajuste iterativo, resultó en una herramienta potente para la optimización de la capacidad predictiva de los modelos de altura dominante.
Bates y Granger (1969), Morris (1974) y Bunn (1989), sugieren que el uso de pronósticos de modelos previos como variables independientes, equivale a utilizar toda la evidencia disponible y fuentes de información sólidas, para mejorar los pronósticos de los modelos nuevos. Asimismo, en la presente investigación se aplicó la suposición fundamental de que no se puede identificar exactamente el verdadero proceso, pero diferentes modelos pueden desempeñar un papel complementario en la generación de datos (Terui y Van Dijk 2002), como en el caso de la altura dominante de los rodales. La estrategia fue similar a la implementada por Nava-Nava et al. (2022), quienes señalan que la altura de los árboles no se puede conocer con certeza, incluso con mediciones de campo convencionales, por lo que la estimación con modelos permite una fuente de datos adecuada. Por tanto, los resultados obtenidos en la presente investigación corroboran la robustez derivada de la combinación empírica de las fuentes de información.
Los estadísticos de bondad de ajuste obtenidos, coinciden con los generados en otros estudios para especies de coníferas (González-Méndez et al. 2016, Socha y Tymińska-Czabańska 2019, Nava-Nava et al. 2020, Guzmán-Santiago et al. 2021). Los coeficientes de determinación (R 2 adj) de los modelos ADA y GADA (figuras 2 y 3) muestran que es posible explicar más del 99 % de la variabilidad en altura dominante para Pinus oaxacana, valor similar al obtenido por Castillo-López et al. (2018) para la especie en Oaxaca, México.
En general todos los modelos del presente estudio mostraron estadísticos de ajuste satisfactorios y parámetros significativos, sin embargo, el modelo de Hossfeld GADA fue descartado para su elección debido a que presentaba errores estándar altos de hasta 3.822,0 en alguno de sus parámetros (cuadro 6). Este error alto también ha sido reportado por Castillo-López et al. (2013), Galindo-Soto et al. (2017) y Nava-Nava et al. (2020).
Los modelos seleccionados (figura 6) presentaron valores similares para los estadísticos SCE, EMC y Ē en la validación (figura 4 y 5). Para los modelos ADA se obtuvieron intervalos de 5,4 - 22,5, 0,1 - 0,2 y 0,001 - 0,2, mientras que en GADA fueron de 5,5 - 21,9, 0,1 - 0,2 y 0,03 - 0,1, para cada indicador respectivamente. Sin embargo, los modelos pueden presentar los mismos estadísticos, pero proyecciones gráficas diferentes (Vargas-Larreta et al. 2010). Finalmente, mediante el análisis gráfico se seleccionaron los modelos GADA, porque las curvas presentaron una flexibilidad mayor (figura 6C y 6D). La flexibilidad de estas curvas se debe a que más de un parámetro depende de la calidad de estación, son polimórficas y con múltiples asíntotas (Cieszewski 2002) y se adaptan fácilmente a la modificación de la edad de referencia (Cieszewski y Bailey 2000).
El modelo GADA de Korf por su parte, ha sido ensayado por Vargas-Larreta et al. (2013) y Castillo-López et al. (2018) para modelar el crecimiento en altura dominante de diferentes especies de pinos en México. Asimismo, en otras partes del mundo también se ha utilizado, por ejemplo, Martín-Benito et al. (2007), estudiaron la variabilidad interregional de los modelos de crecimiento en altura dominante e índice de sitio de Pinus nigra Arn. en la Península Ibérica, seleccionando el modelo de Korf, al obtener resultados mejores, tanto en el ajuste como en la validación.
El modelo GADA de Chapman-Richards ha sido ampliamente utilizado para modelar el crecimiento en altura dominante, por ejemplo, Diéguez-Aranda et al. (2006) construyeron curvas de índice de sitio para Pinus taeda L. en plantaciones de Estados Unidos; Vargas-Larreta et al. (2010) ajustaron modelos dinámicos de índice de sitio para Pinus cooperi Blanco en Durango, México; Tamarit-Urias et al. (2014) ajustaron tres modelos para Tectona grandis L., en el valle de Edzná, Campeche, México, incluyendo el modelo de Chapman-Richards; por su parte, Quiñonez-Barraza et al. (2015), a partir de este modelo, ajustaron ecuaciones de índice de sitio con polimorfismo complejo para masas forestales de Durango, México. En tanto que, Hernández-Cuevas et al. (2018), ajustaron modelos de índice de sitio para la clasificación de la productividad de rodales de P. ayacahuite, empleando entre otros, al modelo de Chapman-Richards. Más recientemente, Nava-Nava et al. (2020) proponen dicho modelo en su formulación GADA para la evaluación de la productividad forestal de P. patula, en la región de estudio.
Los modelos mejores fueron Korf y Chapman-Richards en su formulación GADA ajustados con las proporciones del 80 % y 70 % de la base de datos, respectivamente, los cuales presentaron un valor para RCME de 0,3665 m y 0,3720 m, respectivamente, con valores altos en el coeficiente de determinación ajustado (R2 adj) (figura 3). Otro aspecto crucial en la selección de los modelos mejores fue el comportamiento gráfico, los modelos elegidos agruparon adecuadamente el intervalo de dispersión de los datos observados (figura 6), con asíntotas que oscilan entre los 30 y 40 m, resultados similares fueron obtenidos por Nava-Nava et al. (2020) para P. patula en el área de estudio y por Castillo López et al. (2018) para la especie en la región.
De acuerdo con Cieszewski (2002) y Vargas-Larreta et al. (2010), las ecuaciones GADA generan curvas polimórficas y a la vez anamórficas, así como invariantes respecto a la edad referencia y con respecto al camino de simulación, presentan gran flexibilidad y a partir de un modelo base se pueden derivar otros, permitiendo que más de un parámetro dependa de la calidad de estación. Por estas ventajas, además del comportamiento gráfico, y los estadísticos de ajuste y sobre todo de validación, se optó por seleccionar la ecuación GADA de Chapman-Richards ajustada con el 70 % de los datos (CG-70).
La ecuación de Chapman-Richards representa de forma adecuada los patrones de crecimiento en altura dominante y esta eficiencia se debe a la flexibilidad del modelo y a la representación adecuada de las trayectorias de las curvas sobre los datos observados (Vargas-Larreta et al. 2010, Nava-Nava et al. 2020). En la modelación de la altura dominante para P. oaxacana y P. patula, Castillo-López et al. (2018) y Nava-Nava et al. (2020), respectivamente, seleccionaron la ecuación de Chapman-Richards.
En la región de la Sierra Norte de Oaxaca, México, existen modelos de índice de sitio para P. oaxacana, como el ajustado por Castillo-López et al. (2018), quienes utilizaron información proveniente de análisis troncales, pero a una escala gruesa de análisis, o el modelo anamórfico para P. oaxacana reportado en el programa de manejo de Ixtlán de Juárez, Oaxaca (STF, 2015) y construido con datos de corte transversal, sin embargo, se sabe que este tipo de modelos (anamórficos) tiene ciertas limitaciones en cuanto a las estimaciones, dado que establecen el mismo turno forestal en los rodales, independientemente de su productividad, por lo que no se puede optimizar el potencial real del sitio (Clutter et al. 1983). Por tanto, una ecuación a nivel predial y con datos de tipo longitudinal puede abarcar condiciones de sitio más específicas y apegarse con mayor fidelidad a la dinámica de las masas forestales, de modo que pueda combinarse con herramientas para el manejo de la densidad y modelos de crecimiento y rendimiento, lo que permitirá aplicar el régimen de manejo más apropiado para optimizar la productividad forestal (Santiago-García et al. 2013). Tratándose de una comunidad modelo a nivel mundial, en el manejo sustentable de los recursos forestales, como es el caso de la comunidad de Ixtlán de Juárez, Oaxaca, contar con herramientas propias para la toma de decisiones es una necesidad apremiante (Bray et al. 2005, Klooster et al. 2005, Torres-Rojo et al. 2016, Santiago-García et al. 2022).
CONCLUSIONES
El análisis gráfico, los estadísticos de bondad de ajuste y principalmente la validación de los modelos permitió seleccionar las ecuaciones en diferencia algebraica (ADA): Korf anamórfico b 0 (70 %) y Korf polimórfico b 1 (80 %) y las ecuaciones en diferencia algebraica generalizada (GADA): Chapman-Richards (70 %) y Korf (80 %). De estos, el modelo GADA de Chapman-Richards presentó capacidad predictiva mayor y fue elegido como el mejor para la estimación de la altura dominante y la evaluación del potencial productivo de los rodales de P. oaxacana, en bosques de Ixtlán de Juárez, Oaxaca, México, constituyendo así una herramienta importante para la toma de decisiones de manejo forestal.