Ideas destacadas:

Introducción

La estimación del potencial productivo de sitios forestales es importante para el manejo forestal porque ayuda a determinar la cosecha anual, establecer el turno y seleccionar las especies más adecuadas para un área (^{Burkhart & Tomé, 2012}). La productividad forestal como concepto biológico complejo se estima indirectamente por medio de modelos de crecimiento (^{Tamarit-Urias et al., 2014}). La forma habitual para referenciar la productividad es utilizar la relación altura dominante-edad mediante el índice de sitio (IS) que se define como la altura promedio que alcanzan los árboles dominantes de un rodal a una edad base determinada (^{Quiñonez-Barraza et al., 2015}). Con las ecuaciones de IS se construyen curvas que muestran el patrón de crecimiento en altura de los árboles de una especie en una ecorregión (^{Akbas & Senyurt, 2018}).

Las curvas de crecimiento que las ecuaciones dinámicas generan deben cumplir con los atributos siguientes: realismo biológico, crecimiento sigmoide, punto de inflexión, polimorfismo, asíntotas variables, invarianza de la edad base y del camino de simulación, y la altura a la edad base debe ser igual al IS (^{Guerra-Hernández et al., 2021}). La invarianza de la edad base se ha formalizado mediante el enfoque de diferencia algebraica, donde un solo parámetro del modelo base se expresa en función del sitio; sin embargo, este enfoque tiene como principal limitante que las ecuaciones derivadas generan únicamente curvas anamórficas o polimórficas, con lo cual se restringe la posibilidad de describir una variedad mayor de condiciones de crecimiento comunes en suelos forestales (Cieszewski, 2021).

En México, pese a las limitaciones de la diferencia algebraica, en las últimas décadas del siglo pasado y aún en la primera del presente, los manejadores forestales utilizaron este enfoque para desarrollar ecuaciones de IS para especies del género Pinus. Por tanto, es necesario actualizar esos modelos mediante enfoques y técnicas de modelación más flexibles. Con ello, de acuerdo con ^{Seki y Sakici (2017}) y ^{Özçelik, Cao, Gómez-García, Crecente-Campo, y Eler (2019}), se garantizará mayor fiabilidad y precisión en las estimaciones del crecimiento de la altura dominante e IS. Para lograr este propósito, se cuenta con el enfoque de diferencia algebraica generalizada (DAG). En este enfoque, el modelo base se expande para que dos de sus parámetros se expresen en función del sitio mediante una variable teórica X, la cual se asume que engloba los factores de la productividad del sitio. La DAG permite derivar ecuaciones dinámicas que reproducen curvas de IS polimórficas con múltiples asíntotas de mayor realismo biológico (Cieszewski, 2021; ^{Socha, Tyminska-Czabańska, Grabska, & Orzel, 2020}; ^{Sprengel, Spiecker, & Wu, 2022}).

Pinus chiapensis (Martínez) Andresen es una conífera subtropical endémica de México y Guatemala que forma parte del bosque mesófilo de montaña en Puebla, Veracruz, Guerrero, Chiapas y Oaxaca. La especie destaca por su crecimiento rápido y potencial maderable alto; se distribuye en forma natural en altitudes de 600 a 2 200 m, con una precipitación anual promedio de 1 500 a 2 000 mm temperaturas de 20 a 22 °C y pendientes pronunciadas con suelos ácidos (^{Perry, 2009}). Esta especie desempeña funciones ecológicas y económicas clave en el neotrópico, es un recurso genético valioso (^{Del Castillo, Trujillo, & Saénz-Romero, 2009}) y está catalogada con estatus de protección especial (^{Chávez-Pascual et al., 2013}). La generación de una nueva función de IS mediante DAG podrá ser de capacidad predictiva mayor que los modelos polimórficos previos derivados de diferencia algebraica.

El objetivo de este trabajo fue desarrollar una ecuación dinámica de IS mediante el enfoque de diferencia algebraica generalizada como una herramienta silvícola para describir el patrón de crecimiento en altura dominante y clasificar la productividad de rodales naturales de P. chiapensis de Puebla y Veracruz, México.

Materiales y métodos

Resultados y discusión

El ajuste de las ecuaciones dinámicas indicó que todos los parámetros globales y los de autocorrelación (ρ _i ) fueron significativos al nivel significancia del 5 % (Cuadro 3). Los estadísticos de bondad de ajuste, mostrados en el Cuadro 4, revelan que el Radj2, la RCME y el AIC son próximos entre sí. Todos los modelos explicaron la varianza total de la altura dominante en función de la edad en un nivel mayor de 98.08 %. El sistema de calificación cuantitativo evidenció que la ecuación con los mejores estadísticos fue la M4 (modelo de Hossfeld), la segunda mejor fue la M2 seguida de M3, M6, M5 y al final M1. En términos del sesgo destacan las ecuaciones M1, M5 y M6, ya que comparativamente presentaron los valores más bajos. Con relación al estadístico DW, los valores fluctuaron de 1.61 a 1.92, los cuales son próximos a 2, lo que de acuerdo con ^{Trim et al. (2020}), es evidencia de que el problema de autocorrelación se redujo al mínimo. El orden de la estructura autorregresiva AR(ρ) que se indica para cada ecuación, correspondió a aquella en donde se logró alta significancia en todos los parámetros, convergencia plena con un nivel de exigencia de 0.0001 y el mejor comportamiento gráfico de las curvas de IS.

Los residuales estandarizados contra los predichos de las ecuaciones mantuvieron una distribución con tendencia aleatoria en todos los casos (Figura 1), lo que sugiere que la heterocedasticidad se redujo al mínimo para lograr un patrón homocedástico; por lo tanto, se infiere que las ecuaciones evaluadas se ajustan de forma aceptable a los datos observados. Bajo este escenario, ambas correcciones permiten que los estimadores obtenidos en el ajuste sean de varianza mínima, insesgados y más eficientes, lo que mejora sustancialmente las propiedades estadísticas y predictivas de las ecuaciones. De este modo, es posible estimar intervalos de confianza realistas y efectuar los contrastes estadísticos de significación sobre los parámetros con propiedad.

Figura 1 Comportamiento gráfico de los residuales estandarizados con respecto a los valores predichos de la altura dominante de Pinus chiapensis para las ecuaciones evaluadas. 

La Figura 2 muestra las tendencias de crecimiento en altura dominante que las ecuaciones dinámicas reproducen para valores de IS de 22, 27, 32, 37 y 42 m (Eb = 50 años), sobrepuestas a las trayectorias de crecimiento observadas. Una inspección detallada muestra que las funciones M2 y M3 tienden a sobrestimar las alturas dominantes en la mejor calidad de sitio en edades superiores a 60 años; antes de esta edad, la tendencia sugiere ligera subestimación. Por su parte, las ecuaciones M4 y M5 tienden a ligeras subestimaciones en las mejores calidades de IS para todo el intervalo de edad. La ecuación M1 exhibió el comportamiento más deficiente, ya que subestima en todos los IS a edades tempranas. Esta tendencia de subestimación se mantiene para los mejores sitios a edades maduras y longevas.

En contraparte, la ecuación dinámica M6 (modelo base de Levakovic II) mostró mejor realismo biológico, porque las tendencias de crecimiento que describe se sobreponen de mejor manera a las trayectorias observadas en todos los niveles de IS y en todo el intervalo de edad analizado. Por esta razón, la ecuación M6 se seleccionó para describir el patrón de crecimiento en altura dominante y para clasificar la productividad de rodales naturales de P. chiapensis en la zona de estudio. Esta selección ratifica lo señalado por ^{Özçelik et al. (2019}), en el sentido de que las funciones con los mejores estadísticos de ajuste no siempre reproducen de mejor manera las tendencias observadas. En este caso, estadísticamente, la ecuación M6 se posicionó en el antepenúltimo lugar; sin embargo, dado el realismo biológico que recrea, es claro que tiene la capacidad predictiva más alta. Esta elección más que un ejercicio puro de inferencia estadística fue producto de un compromiso y balance entre desempeño estadístico y razones biológicas. La ecuación M6 tiene parsimonia con solo dos parámetros, suficientes para revelar las verdaderas tendencias de crecimiento en las alturas dominantes observadas.

Figura 2 Comportamiento de las curvas de índice de sitio (IS) reproducidas por las ecuaciones dinámicas evaluadas sobrepuestas a las trayectorias de crecimiento observadas en Pinus chiapensis. 

Con los valores de los parámetros de la ecuación dinámica M6 y una Eb = 50 años, se construyeron curvas de IS para tres niveles de productividad etiquetadas como baja, media y alta, cuyos valores medios corresponden a los IS de 25, 32 y 39 m, respectivamente (Figura 3a). Estas curvas junto con sus límites respectivos, en particular el inferior y superior de las clases baja y alta, cubrieron la amplitud y las tendencias de todas las series de crecimiento en el intervalo de edad analizado. Por esta razón, las clases de IS establecidas se consideraron suficientes para la clasificación de la productividad de los rodales naturales de P. chiapensis en la zona de estudio.

El enfoque de DAG para generar la ecuación M6 en combinación con la técnica de variables Dummy para ajustarla, aunado a las correcciones de autocorrelación y heterocedasticidad, logran una ecuación dinámica robusta para realizar predicciones fiables. Este aspecto se aprecia mejor al contrastar las curvas de IS con las que reproduce el modelo de IS propuesto por ^{Rodríguez-Acosta y Arteaga-Martínez (2005}) (Figura 3b), que corresponde a un modelo polimórfico, cuyo modelo base es el de Chapman-Richards. Se evidencia que para una Eb = 50 años y los valores de IS referidos, el modelo polimórfico no describe adecuadamente el crecimiento en altura dominante de P. chiapensis. La mayor deficiencia se aprecia para la clase de IS alta, donde las tendencias de crecimiento son ilógicas porque las curvas colapsan, debido a que el valor de la asíntota horizontal del modelo es menor que la altura real que los datos observados sugieren. Las diferencias pueden explicarse desde el punto de vista matemático-estadístico, porque para generar el modelo polimórfico usaron el enfoque de diferencias algebraicas, en el que solo el parámetro de forma fue dependiente del sitio y, además, no se corrigieron la autocorrelación ni la heterocedasticidad.

Las alturas dominantes de 25, 32 y 39 m estimadas por la ecuación nueva, que equivalen a los valores de IS (Eb = 50 años), son ligeramente superiores a 24.2, 29.7 y 35.5 m (Eb = 50 años), alturas determinadas mediante el modelo propuesto por ^{Sánchez y Del Castillo (2001}). Estos autores utilizaron datos de sitios temporales y ajustaron modelos lineales por condición de crecimiento en la estimación del IS para rodales naturales coetáneos de P. chiapensis en la Sierra Norte de Oaxaca, México. Esta discrepancia sugiere que las condiciones de crecimiento dadas por condiciones físicas y ambientales en los sitios de Veracruz y Puebla son comparativamente mejores. En este sentido, Del Castillo et al. (2009) refieren que árboles individuales longevos de esta especie, en las mejores calidades de estación, pueden alcanzar alturas de hasta 50 m, tal como se evidencia en las observaciones del presente estudio para la clase alta de IS.

Figura 3 Comparativo gráfico de las curvas de índice de sitio (IS) generadas por (a) la nueva ecuación dinámica y (b) el modelo de ^{Rodríguez-Acosta y Arteaga-Martínez (2005)}. 

Un análisis de la función polimórfica derivada mediante diferencias algebraicas con el modelo base de Hossfeld IV, desarrollado por F^{ierros-Mateo, De los Santos-Posadas, Fierros-González, y Cruz-Cobos (2017}) para la misma especie, mostró que con una Eb = 8 años y clases de IS de 6, 10 y 14 m, se pudo modelar el crecimiento en altura dominante de los datos observados en forma aceptable hasta los 20 años; posteriormente, tiende a un techo asintótico máximo de 33 m, el cual es menor que las alturas máximas observadas. Esta disparidad, aun cuando la zona de estudio fue similar (UMAFOR 2103), puede explicarse porque en dicho estudio se procesó información de plantaciones comerciales jóvenes con un intervalo de edad de 1 a 8 años.

Con valores de IS de 25, 32 y 39 m y Eb = 50 años, la ecuación dinámica en DAG del modelo base de Chapman-Richards [a=exp⁡() y c=(β2+β3)/], reportada en la plataforma SiBiFor por ^{Vargas-Larreta et al. (2017}) para la misma especie en la UMAFOR 1203 del estado de Guerrero, tiende a realizar predicciones superiores a las del modelo M6, sobre todo en la calidad de IS media y alta; en tanto que, con el IS bajo subestima hasta los 20 años y después sobrestima. La sobrestimación para los tres IS se mantiene hasta los 50 años, después tiende a techos asintóticos de 26, 33 y 40 m para los IS bajo, medio y alto, respectivamente; tales valores son inferiores a las trayectorias máximas observadas. En este sentido, aun cuando en el presente trabajo como en el de Vargas-Larreta et al. (2017) se utilizaron ejemplares con régimen de regeneración natural y la técnica de análisis troncal, el comportamiento sugiere que las calidades de estación en el área de estudio son mejores.

Las expresiones 1 y 2, cuyo componente X ₀ está indicado en el Cuadro 2, derivadas de la ecuación dinámica M6, son útiles para construir curvas de IMA e ICA por clase de IS. En tanto que las expresiones 3 y 4, derivadas de la misma ecuación, sirven para determinar el incremento medio anual máximo (IMA_máx) y el incremento corriente anual máximo (ICA_máx) por nivel de productividad.

IMA=exp(X0)tβ2+tβ1/0t …………….. (5)

ICA=β1 β2tβ2+tβ1/0exp(0)t(β2+t)0…………….. (6)

IMAmáx=β2β1-00…………….. (7)

ICAmáx=12β2β1-00…………….. (8)

Las determinaciones de ICA, IMA y los incrementos máximos evidencian que la ecuación M6 satisface consistentemente el supuesto lógico de que la edad a la que se alcanzan tales incrementos decrece conforme el IS aumenta (Cuadro 5). También se cumple que los incrementos son mayores a medida que el valor del IS es mayor. De acuerdo con ^{Seki y Sakici (2017}) y ^{Socha et al. (2020}), se ratifica que la ecuación dinámica desarrollada cumple a cabalidad con las leyes de crecimiento biológico.

El hecho de que se cumplan estos supuestos relevantes coincide con lo referido por ^{Ercanli, Kahriman, y Yavuz (2014}) y por ^{Seki y Sakici (2017}), en el sentido de que las ecuaciones dinámicas derivadas mediante DAG poseen mayor flexibilidad para modelar una variedad amplia de tendencias de crecimiento. Tales ecuaciones cumplen con todos los atributos deseables que se traducen en consistencia y realismo biológico, por lo que deben de ser utilizadas para modelar el crecimiento de bosques naturales y artificiales. En este contexto, se considera que, para el desarrollo de nuevas ecuaciones dinámicas de IS en México, se debe priorizar el uso del enfoque de DAG y que sean dos los parámetros específicos del sitio en cada modelo de crecimiento base que se utilice. Además, se debe usar algún método de ajuste que sea invariante de la edad base, tal como el de variables Dummy (^{Sghaier et al., 2012}), el procedimiento iterativo anidado (^{Sharma, Brunner, Eid, & Øyen, 2011}) o bien la técnica de modelo de efectos mixtos (^{Tamarit-Urias et al., 2014}). También se debe corregir la autocorrelación y la heterocedasticidad para obtener modelos más robustos que cumplan con los supuestos teóricos del análisis de regresión. Esto permitirá estar al nivel de la tendencia internacional ya que, como refieren ^{Özçelik et al. (2019}), es común que cuando se disponga de datos provenientes de análisis troncales o de parcelas permanentes con remediciones se usen los refinamientos matemáticos y estadísticos referidos.

Mediante el contraste de las edades en que ocurren los turnos técnicos (IMA = ICA) en la altura dominante, por clase de IS de la ecuación dinámica M6, se infiere que P. chiapensis tiende a tasas constantes de crecimiento en la zona de estudio. Lo anterior debido a que las diferencias en la ocurrencia son pequeñas, efecto que conduce a que las curvas de IS pudieran considerarse de tipo cuasianamórficas. Esta tendencia de crecimiento también fue observada por ^{Calegario, Daniels, Maestri, y Neiva (2005}) en clones del género Eucalyptus en los estados Bahía y Espírito de Brasil.

Conclusiones

La ecuación dinámica desarrollada de índice de sitio, derivada del modelo de crecimiento de Levakovic II, tiene base biológica y rigor matemático-estadístico sólido; además, genera curvas polimórficas con asíntotas variables que se apegan con fidelidad a la realidad del crecimiento observado. Asimismo, la ecuación cumple con atributos deseables lógicos entre los que destacan realismo biológico, invarianza con respecto a la edad base y de ruta, flexibilidad, parsimonia y consistencia. Se recomienda su uso como una herramienta actual y fiable para la gestión silvícola en términos de la predicción del crecimiento en altura dominante y determinación del potencial productivo de rodales naturales de Pinus chiapensis en Puebla y Veracruz, México. La ecuación podrá utilizarse como componente básico de sistemas de crecimiento y rendimiento maderable. Esta ecuación nueva presentó mejor desempeño dado por una capacidad predictiva superior, comparada con un modelo polimórfico de índice de sitio previo para la misma especie y región.