Introducción

La estimación o cálculo del volumen de árboles individuales es una variable fundamental para la planeación de las actividades silvícolas en el programa de manejo forestal con el que se dictan las pautas silvícolas a seguir, en virtud de que se basa en el crecimiento y rendimiento de cada especie (^{Cruz et al., 2016}). Una de las alternativas más viables para estimar y predecir el volumen del arbolado a partir de los datos de un inventario forestal es el uso de ecuaciones que toman como variables independientes al diámetro normal (Dn) y la altura total (At), las cuales son de fácil y rápida medición en campo (^{Acosta & Carrillo, 2008}; ^{Gómez et al., 2018}).

A través del uso de las variables Dn y At, se han propuesto modelos que estiman el volumen fustal (Vf) para diversas especies de interés maderable en condiciones variables de crecimiento (^{Cruz et al., 2016}; ^{Gómez et al., 2018}; ^{Tschieder et al., 2011}). De igual manera, se han sugerido funciones para modelar la disminución del diámetro a medida que se acerca a la At sobre el fuste de cada árbol (i.e. ahusamiento, d) (^{Maldonado-Ayala & Návar, 2000}; ^{Ramírez-Martínez et al., 2018}; ^{Rentería-Anima & Ramírez-Maldonado, 1998}), estimar el volumen comercial por individuo (Vc) y clasificar la distribución de productos maderables del arbolado en pie de acuerdo a sus dimensiones (^{Hernández et al., 2013}), así como representar la variación de la forma en distintas procedencias en una misma especie (^{Tlaxcala-Méndez et al., 2016}) o la razón de proporción entre el Vc y Vf (^{Hernández-Ramos et al., 2018b}).

En el aprovechamiento maderable en zonas tropicales de México, Quintana Roo es un estado de importancia forestal en términos de extracción maderable, debido a que se registró una tasa de aprovechamiento promedio anual de 45 151 m³ rollo total árbol (rta) entre los años de 1990 a 2017 (^{Secretaría del Medio Ambiente y Recursos Naturales [Semarnat], 2018}), la cual principalmente se obtiene en el centro y sur de la entidad. El estado ocupa el segundo lugar nacional en la extracción maderable de especies tropicales, que junto con Veracruz (1°) y Tamaulipas (3°) representan el 72.4% de la producción nacional (^{Semarnat, 2016}). De acuerdo con un análisis del volumen autorizado entre 2016 y 2020 por la Semarnat, el volumen proveniente de Manilkara zapota (L.) P. Royen contribuye con un promedio de 7%.

Por el volumen potencial para el aprovechamiento maderable dentro de las unidades de manejo forestal (Umafor), 2308 y 2304 en Quintana Roo (24.24 m³ ha^-1), y la diversidad de condiciones de crecimiento que tiene M. zapota, lo cual la hace ser de las especies tropicales más abundantes de las selvas en el estado (^{Gallegos et al., 2012}; ^{Semarnat-Comisión Nacional Forestal [Conafor], 2014}), vale la pena hacer hincapié en el uso de herramientas silvícolas cuantitativas aplicadas en el manejo forestal maderable para el aprovechamiento de los bosques tropicales en el centro y sur de esta entidad. Así, el objetivo de este trabajo fue evaluar y seleccionar modelos de Vf y d para árboles de Manilkara zapota del centro y sur de Quintana Roo, México.

Materiales y métodos

El estudio se realizó en el centro y sur de Quintana Roo, México, entre las coordenadas 17°50´57.48´´ y 20°12´14.72´´N, y 87°10´57.19´´ y 89°16´18.35´´ O (Figura 1). El tipo de vegetación fue selva mediana sub-perennifolia y selva baja caducifolia. Los principales géneros que conforman la selva mediana son Bursera, Lysiloma, Picidia, Coccoloba, Pouteria, Gymnanthes, Lonchocarpus, Manilkara, Metopium y Vitex; mientras que Metopium, Coccoloba, Manilkara, Haematoxylum, Bucida, Crotan, Bursera, Chrysobalanus, Lonchecarpus y Cameraria son parte de las selvas bajas. El clima predominante es de tipo cálido sub-húmedo (Aw), con temperatura media anual de 26 °C, precipitación promedio de 1300 mm y altitud promedio de 15 m (^{Semarnat-Conafor, 2014}).

Fuente: Elaboración propia.

Figura 1 Ubicación del área de estudio. 

El muestreo se llevó a cabo en los años 2014 y 2015 en los ejidos Caobas, Chan Santa Cruz, Felipe Carrillo Puerto, Naranjal Poniente, X-Hazil sur y el sitio experimental San Felipe Bacalar del INIFAP, cubriendo la mayor variabilidad posible de condiciones de crecimiento de la especie dentro de las selvas mediana sub-perennifolia y baja caducifolia. De manera aleatoria, se seleccionó una muestra de 186 árboles de M. zapota en rodales de las áreas forestales en donde este taxón tuvo presencia; los ejemplares seleccionados presentaron fuste recto, sin malformaciones, copa completa, libre de plagas o enfermedades y daños mecánicos por vientos, huracanes o incendios.

De cada árbol se registró la información de diámetros del tronco y alturas de distintas secciones sobre el fuste (di y Ai). El diámetro se midió de forma directa con una forcípula graduada en milímetros, iniciando con el diámetro y altura del tocón (d _tocón y A _tocón , a 0.3 m de altura sobre el nivel del suelo), y después a 0.6 m, 0.9 m, 1.30 m y 2.5 m de altura del fuste. El resto de las mediciones sobre el fuste se realizaron con un telerelascopio de Bitterlich® hasta llegar a la At = 0 y d = 0 (i.e. punta del árbol), para obtener 1666 pares de datos de di y Ai.

El cálculo del volumen por troza (Vt) y Vf de cada individuo se realizó a través del método de trozas traslapadas (^{Bailey, 1995}). Con la información completa y depurada a través de eliminar datos atípicos, aberrantes o fuera de la tendencia general de la información, se calculó la estadística descriptiva básica de la muestra, se estimaron el índice de curtosis (IC) y el coeficiente de asimetría (CA) de cada variable; los cuales deben de estar en el intervalo de 3 y -3 para considerar que los datos cumplen con el supuesto de la distribución normal (^{Galán et al., 2016}; ^{Martínez-González et al., 2006}).

Se evaluaron ocho modelos para estimar el Vf en función del Dn y la At (^{Cruz et al., 2016}; ^{Gómez et al, 2018}; ^{Tschieder et al., 2011}), así como dos funciones para modelar el d en función de Dn, At y Ai (^{Ramírez-Martínez et al., 2018}; ^{Rentería et al., 2006}), los cuales han demostrado ser eficientes para estimar estas variables en estudios semejantes (Tabla 1).

Las funciones matemáticas se ajustaron con el programa estadístico RStudio Version 1.1.463® mediante la función de modelos lineales (nl) y no lineales (nls), con la técnica de máxima verosimilitud (ML) (R Core Team, 2016). En los ajustes de ahusamiento se incorporó una variable indicadora en la punta (ϑ = 0.001) para mejorar la convergencia del ajuste (^{Uranga-Valencia et al., 2015}). Para corregir la potencial heterocedasticidad, durante el ajuste de Vf y d se incluyeron las variables de ponderación a los residuales de Dn ² , 1/Dn ² , 1/raíz(Dn ² ), Dn ² At, raíz(Dn ² At), 1/Dn ² At y 1/raíz(Dn ² At); con base en el procedimiento propuesto por ^{Harvey (1976)}, y empleado por ^{Parresol (1993)} y ^{Hernández-Ramos et al. (2017b)} para ajustar funciones de d o generar sistemas compatibles de ahusamiento-Volumen comercial (d-Vc). El cumplimiento del supuesto de normalidad en la frecuencia de los residuales se evaluó de forma gráfica en ambas variables dependientes (^{da Cunha et al., 2009}).

El mejor modelo para estimar tanto Vf como d se seleccionó con base en varios criterios: Se consideró a la variable de ponderación que minimizara la heterocedasticidad; que todos los parámetros fueran significativos (p = 0.05); que presentara los mayores valores en el coeficiente de determinación (ecuación 11: R ² ); los menores valores en la raíz del error cuadrático medio (ecuación 12: RMSE); y los mejores valores de verosimilitud en los criterios de Akaike, Bayesiano y su logaritmo (ecuación 13: AIC; ecuación 14: BIC y LogLik) (^{Findley, 1991}; ^{Martínez-González et al., 2006}; ^{Ramírez-Martínez et al., 2018}; ^{Uranga-Valencia et al., 2015}). Como medida de evaluación de la precisión de los modelos, se calculó el sesgo (ecuación 15: E) y el error medio (ecuación 16: Em) (^{Prodan et al., 1997}; ^{Tlaxcala-Méndez et al., 2016}).

donde yi, yi^ y y- son el volumen fustal (Vf) y ahusamiento (d), observado, estimado y promedio, respectivamente; n es el número de observaciones; k es el número de parámetros de cada ecuación; log es el logaritmo natural; y Lik es la verosimilitud.

Después de seleccionar los mejores modelos, se verificó su desviación con respecto a los datos observados mediante un gráfico de dispersión y una regresión lineal, comparando el valor de R ² obtenido entre ellos (^{Gómez et al., 2018}; ^{Medina-Peralta et al., 2010}; ^{Tedeschi, 2006}).

Para la ecuación de ahusamiento seleccionada, con la finalidad de comparar los resultados en la estimación del volumen, se derivó la ecuación de Vf implícita mediante la integración del diámetro a distintas secciones sobre el fuste (d) con respecto a Ai como un sólido en revolución mediante la expresión (17):

donde k = pi/40000 es la constante volumétrica, con pi = 3.1416, y δAi es el diferencial sobre la altura del fuste, mientras que A ₁ y A ₂ son las alturas que definen la integral y que pueden variar de A ₁ = 0 hasta A ₂ = At (altura total).

Se realizó la comparación entre las estimaciones obtenidas al aplicar la ecuación de Vf, seleccionada como la mejor, y la expresión implícita de la función de ahusamiento (AhusModVol) para verificar si el modelo implícito en el ahusamiento puede ser utilizado de forma confiable en la estimación de Vf o tiene menores desviaciones que el ajuste de un modelo tradicional (Tabla 1, ecuaciones 1-8).

Resultados

La estadística descriptiva de la muestra indicó un intervalo de aplicabilidad de los modelos en las categorías diamétricas de 10 cm a 85 cm, con altura total promedio de 14.5 m; el Dn presentó los valores más altos en el error típico, desviación estándar y la varianza. Además, la distribución se consideró simétrica en todas las variables, mientras que la curtosis resultó negativa para Dn y At, y positiva en Vf, sin problemas de desviaciones (Tabla 2).

En el ajuste estadístico de los modelos para volumen fustal y ahusamiento, se determinó que la variable de ponderación que presentó la mejor corrección de la heterocedasticidad fue 1/raíz(Dn ² ); ya que los resultados indican una mejor distribución de residuales, además de parámetros significativos en la mayoría de los modelos (p = 0.05), con excepción de a ₀ y b ₃ en las ecuaciones 7 y 9, respectivamente; mismas que no fueron utilizadas para el análisis posterior. Los valores del factor de forma fueron de 0.537 a 0.503 (a ₀ : ecuación 2; a ₁ : ecuación 5) (Tabla 3), superiores al promedio de la muestra y al obtenido por la expresión que se genera al integrar la expresión de ahusamiento (10) que fue de 0.484.

Con base en los criterios de bondad de ajuste, se determinó que las expresiones de Schumacher-Hall (ecuaciones 1) y Cielito 3 (ecuación 10) son las más adecuadas para estimar el Vf y d, respectivamente; ya que presentaron los mejores valores en R ² , criterios de verosimilitud (AIC, BIC y LogLik) y RMSE (Tabla 4). El bajo valor del sesgo (E: 0.00000492 m³) y del error medio (Em: 0.00000871 m³) para Vf ratificó que el modelo 1 es el más adecuado. La integración matemática del modelo de ahusamiento 10 generó un modelo de volumen de factor de forma constante: Vf=π40000∙0.4841994917∙(Dn2At)1, en el cual se observa que hay sobreestimación del Vf (Figura 2), debido a que tanto su sesgo promedio como su error medio son altos (E = -0.01246861 m³ y Em = -0.02206943 m³). Por lo tanto, solo se debe usar como indicador del factor de forma de la especie y no para predecir el volumen fustal de ejemplares de M. zapota. Esto ratifica que para la estimación confiable del Vf se debe de utilizar la expresión de Schumacher-Hall.

Fuente: Elaboración propia.

Figura 2 Errores medios y sesgos de las predicciones de los modelos de volumen fustal para Manilkara zapota en Quintana Roo, México. 

El análisis grafico de los supuestos de homocedasticidad y normalidad de los modelos seleccionados mostraron una adecuada corrección de la distribución de los residuales al incluir la ponderación de 1/raíz(Dn ² ); cuando dicha corrección no se realiza, se observa que los errores aumentan a medida que la variable dependiente es mayor (Figuras 3a y 3c), mientras que, al incluir la ponderación, se evidenció una distribución homocedástica de los residuales (Figuras 2b y 2d). La frecuencia de los residuales expresados en porcentajes presentó una tendencia normal (Figuras 4a y 4b), lo que evidencia el cumplimiento de los supuestos de regresión de normalidad y homocedasticidad.

Fuente: Elaboración propia.

Figura 3 Distribución de residuales sin y con corrección por heterocedasticidad para los modelos de volumen fustal (Vf) (a y b) y ahusamiento (c y d), respectivamente, para Manilkara zapota en Quintana Roo, México. 

Fuente: Elaboración propia.

Figura 4 Pruebas gráficas de normalidad de residuales para los modelos de volumen fustal (Vf) (a) y ahusamiento (d) (b) para Manilkara zapota en Quintana Roo, México. 

El contraste de las estimaciones obtenidas con los modelos de Schumacher-Hall (Vf) y Cielito 3 (d) con los datos observados reprodujo en ambos casos una tendencia de línea recta con desviación global de 6.28% para el Vf (Figura 5a) y de 3.33% en el d (Figura 5b). La forma de la dispersión de las estimaciones para el d con respecto a la altura en distintas secciones sobre el fuste coincidió en los perfiles fustales observados, con tendencia a la reducción del diámetro a medida que se aproxima a la At, sin presentar problemas de convergencia en la punta del árbol (diámetro cero); sin embargo, se observa una mayor variabilidad de forma en la sección cercana al diámetro del tocón (Figura 5c). La tendencia de la estimación en el Vf concuerda con la distribución de los datos observados (Figura 5d).

Fuente: Elaboración propia.

Figura 5 Datos observados vs. predichos de los modelos 1 y 10; volumen fustal (a), ahusamiento (b), tendencias de diámetros a distintas secciones sobre el fuste (c) y al usar el modelo de Schumacher-Hall (d) para árboles de Manilkara zapota en Quintana Roo, México. 

De manera práctica y como una forma de identificar la dinámica de incremento y cambio en el volumen fustal (m³) a medida que los individuos de Manilkara zapota aumentan sus dimensiones y, por consiguiente, se hacen más longevos, se presenta una tabla de volumen de doble entrada, considerando las distintas dimensiones de los árboles evaluados (Tabla 5).

La Tabla 6 muestra la disminución del ahusamiento en los individuos a medida que sus dimensiones se acercan a la altura total y diámetro cero. De igual forma, se observa la variabilidad del diámetro a distintas secciones, propiciada por las diferentes alturas totales en una misma categoría de diámetro normal para árboles de M. zapota, por lo cual su proyección en la distribución de productos maderables será distintinta en cada caso.

Discusión

En modelos de volumen fustal y ahusamiento es importante modelar la varianza del error durante el ajuste utilizando diferentes factores de ponderación, ya que esto permite corregir problemas de heterocedasticidad; así, se pueden realizar mejores estimaciones con el modelo ajustado, dado que los valores de los parámetros que se obtienen tienen varianza mínima y sus estimadores son insesgados y eficientes (^{Harvey, 1976}; ^{Parresol, 1993}).

Los resultados de las ponderaciones realizadas para corregir la heterocedasticidad en las ecuaciones de volumen fustal y ahusamiento son similares a lo realizado por ^{Tschieder et al. (2011)}, quienes utilizaron factores de ponderación en forma de valor recíproco, dados por el Dn y la variable combinada (Dn ² At) al ajustar ecuaciones de volumen total en Populus deltoides Marshall., en donde se obtuvo un menor valor al cuadrado medio del error y una distribución de la varianza homocedástica. También, los resultados son similares a lo reportado por ^{Cruz et al. (2016)}, quienes usaron como factor de ponderación a la variable 1/Dn ^k al ajustar un modelo de volumen total con corteza de Arbutus spp., en donde se eliminó el patrón de distribución heterocedástica de la varianza. De igual forma, los resultados obtenidos concuerdan con lo reportado por ^{Hernández-Ramos et al. (2018a}, ^b), quienes evaluaron diferentes factores de ponderación al ajustar sistemas compatibles de ahusamiento-volumen comercial (d-Vc), así como ecuaciones de volumen fustal (Vf) y de razón de volumen (r) para Swietenia macrophylla King, donde obtuvieron mejores indicadores de bondad de ajuste que cuando no se usaron tales factores. Los contrastes ratifican que con la corrección implementada se mejoran sustancialmente las propiedades estadísticas y predictivas de las ecuaciones.

La explicación de la variabilidad muestral en Vf de la ecuación 1, expresada a través del R ² (93.55%), muestra que la función utilizada es adecuada porque el valor de este estadístico es superior al referido por ^{Alder (1980)} para ser considerado un buen ajuste (>80%). Además, la alta calidad de ajuste obtenido en este trabajo es parecida a la que obtuvo ^{Hernández-Ramos et al. (2017a)}, quienes usaron las ecuaciones 1 y 2 para estimar el Vf en tres clones de Eucalyptus urophylla Blake y obtuvieron valores en el R ² superiores a 0.96, y es similar al de ^{Gómez et al. (2018)} al ajustar modelos de Vf para Hevea brasiliensis Müell Arg. y determinar un valor de R ² = 0.90 con el modelo 1.

Los resultados del ajuste en el ahusamiento, los cuales muestran una explicación del 95.18% de la muestra y valores promedio en la RMSE de 4.3 cm, son similares a los reportados por ^{Hernández-Ramos et al. (2016)} y ^{Ramírez-Martínez et al. (2018)} al modelar el ahusamiento de árboles de S. macrophylla y Pinus ayacahuite Ehren, al utilizar los modelos de ^{Biging (1984)} y Cielito 3 (^{Rentería-Anima & Ramírez-Maldonado, 1998}), respectivamente. Los autores obtuvieron valores de R ² promedio de 96% y desviaciones de 3.2 cm. Se puede inferir que las ecuaciones 1 y 10 tienen gran flexibilidad para explicar la variabilidad muestral en distintas especies con condiciones o hábitats de crecimiento distintos.

La precisión de las estimaciones para volumen fustal, E: 0.00000492 m³, fue mejor que los valores reportados por ^{Montalvo et al. (2001)}, quienes usaron la ecuación 1 para Pinus occidentalis Sw. en República Dominicana (E: 0.20 m³), y fue mejor que los valores reportados por ^{Henry et al. (2004)} al utilizar expresiones de volumen de tipo logarítmica para estimar el Vf de tres especies de coníferas en Cuba (E: 0.0015 m³). La ecuación 1 se ha reportado como superior con respecto a otras expresiones no lineales y lineales, lo cual puede atribuirse a su estructura matemática robusta y parsimoniosa, lo que le permite asimilar de mejor manera la variabilidad de los datos empleados.

El sesgo obtenido para el ahusamiento (E: 0.01194567 cm) fue menor al reportado por ^{Maldonado-Ayala & Návar (2000)}, quienes usaron los modelos de ^{Biging (1984)} y de ^{Newnham (1992)} en cinco pináceas de Durango, México (E: -0.11 cm y -0.42 cm, respectivamente); sin embargo, el sesgo es mayor al reportado por ^{Hernández-Ramos et al. (2018a)}, quienes con un modelo polinómico de quinto orden aplicado en clones de E. urophylla en Huimanguillo, Tabasco, México, determinaron un valor de E = 0.0017 cm. Esta situación probablemente se debe a que en este trabajo se estudió una sola especie, mientras que en las investigaciones referidas se analizó un conjunto de especies y clones; estos últimos presentan menor variabilidad genética que se manifiesta en fenotipos uniformes, lo que redunda y contribuye a que el sesgo sea menor.

La tendencia observada sobre la disminución del d, a medida que se acerca a la At sobre el fuste de las estimaciones de ahusamiento, es lógica y similar a la indicada por ^{Hernández-Ramos et al. (2017b)} y ^{Hernández et al. (2018)}, quienes refieren que el incumplimiento de esta tendencia lógica puede ocasionar predicciones con valores negativos o sesgos considerables al estimar cualquier di y proyectar una distribución de productos.

Los valores de los factores de forma obtenidos (0.537 a 0.503, a ₀ para la ecuación 2; a ₁ para la ecuación 5) para estimar Vf con las ocho ecuaciones de volumen ajustadas (tabla 1) y la expresión de Vf derivada del modelo de ahusamiento (ecuación 10: Vf=π40000∙0.4841994917∙(dn2At)1) de 0.484 fueron menores con respecto a los determinados por ^{De los Santos (1976)} (0.768 para volumen de fuste limpio) y por ^{Patiño et al. (1994)} (0.553 para Vf).

Estas diferencias probablemente se deben a que ^{De los Santos (1976)} construyó ecuaciones para volumen de fuste limpio por grupo de especies, denominadas como duras por la densidad de su madera, y en el cual se incluye a M. zapota (0.86 g/cm²) (^{Ordóñez et al., 2015}), mientras que los factores de forma referidos por ^{Patiño et al. (1994)} corresponden a especies que crecen en el estado de Campeche, en donde las condiciones de crecimiento son distintas en tipo de clima, suelo y vegetación con respecto a los bosques tropicales del centro y sur de Quintana Roo. Sin embargo, las ecuaciones propuestas por ^{De los Santos (1976)} y ^{Patiño et al. (1994)} han sido utilizadas con frecuencia en los programas de manejo de Quintana Roo, situación que refuerza y confirma la necesidad de desarrollar y disponer de una ecuación de volumen específica para este taxón en la entidad.

Al contrastar los factores de forma determinados con respecto a las expresiones propuestas por ^{Vargas-Larreta et al. (2017)} para la misma especie en el estado a nivel de Umafor, se encontró semejanza con el reportado para la Umafor 2301 (0.494), mientras que fue superior a los sugeridos para las Umafor 2302 (0.437), 2303 (0.435) y 2304 (0.407). Estas diferencias pueden ser atribuibles a las distintas condiciones medio ambientales de cada Umafor, en las que crece la especie dentro de Quintana Roo, y a la variabilidad de la muestra utilizada para su construcción.

Las diferencias encontradas entre los resultados de este trabajo con los factores de forma implícitos en cada expresión de volumen empleada en M. zapota por los técnicos forestales en Quintana Roo reportadas en estudios previos (^{De los Santos, 1976}; ^{Patiño et al., 1994}; ^{Vargas-Larreta et al., 2017}) permiten sugerir que, antes de utilizarse, se debe de considerar cada ecuación o factor de forma analizando diversos factores tales como: el objetivo con el que fue generado cada modelo; la amplitud de aplicabilidad de la función en términos del intervalo de clases diamétricas o de altura; la calidad de la información con la que fue construida; las condiciones climáticas, escala y cobertura de la región (estatal o regional: Umafor); las especies o grupos de especie para la que fue propuesta; así como las dimensiones del arbolado en el cual se utilizará. Lo anterior permitirá minimizar errores de estimación y proyecciones de rendimiento maderable al elaborar planes de gestión forestal para la especie en el centro y sur de Quintana Roo.

Conclusión

El modelo de volumen fustal (Vf), que corresponde al de Schumacher-Hall, y el de ahusamiento (d), que corresponde al de Cielito 3, presentaron la mejorar calidad de ajuste y fueron estadísticamente robustos y precisos para ser aplicados de forma confiable en árboles de Manilkara zapota que crecen en los bosques tropicales del centro y sur del estado de Quintana Roo. La aplicación práctica del modelo de volumen es precisamente una tabla de volumen de doble entrada, o bien, su uso directo para estimar inventarios maderables a partir de información de sitios de muestreo. Los factores de forma determinados para este taxón, sobre todo el derivado del modelo de ahusamiento, son biológicamente más realistas y apegados a la realidad observada en la muestra analizada, por lo que pueden sustituir a los reportados el siglo pasado para esta especie, lo que evitará caer en sobrestimaciones al estimar el Vf y al proyectar distribuciones de productos maderables.

Conflicto de intereses

Los autores declaramos no tener conflicto de intereses.