Introducción

Metopium brownei (Jacq.) Urb. (Anacardiaceae) o chechén negro es un árbol de hasta 25 m de altura total (at, m) y diámetro normal (d, cm) de 60 cm. Se distribuye en la vertiente del Golfo de México, desde el sur de Veracruz, Tabasco, la Península de Yucatán y norte de Chiapas. Es abundante en el estrato medio y superior de las selvas medianas subperennifolias y subcaducifolias, con menor abundancia en las selvas bajas; tanto en suelos someros de buen drenaje superficial, como en suelos profundos que soportan inundaciones periódicas (^{Pennington y Sarukhán, 2005}). Por el atractivo veteado de su madera dado por su estructura anatómica y la combinación de colores del duramen y la albura, se aprovecha comercialmente para la fabricación de chapa, duelas, pisos, lambrines, muebles y artesanías con un alto valor en el mercado (^{Tamarit y López, 2007}). Es un taxón de interés para la industria forestal (^{Gallegos et al., 2012}; ^{Silva y Fuentes, 2012}).

El estado de Quintana Roo ocupa el segundo lugar nacional en producción maderable de especies tropicales provenientes de selvas, con un aprovechamiento de 47 139 m³ de madera en rollo de especies duras y blandas tropicales (12.10 % de la producción nacional) y de 4 807 m³ de taxones preciosos (15.18 % de la producción nacional), con una derrama económica importante (^{Semarnat, 2016}).

En la realización de inventarios forestales maderables, el d y la at son las dos variables más importantes (^{Liu et al., 2017}), pero medir la at de todos los árboles es muy costoso en términos de tiempo y recursos. Además, es la variable en la cual se incurre en los mayores errores durante su medición indirecta (^{Hunter et al., 2013}). Por el contrario, la medición del d es relativamente sencilla, precisa y menos costosa (^{Larjavaara y Muller-Landau, 2013}).

La aplicación del manejo forestal requiere del conocimiento del d y la at a nivel de árbol y rodal para estimar los volúmenes existentes en las masas forestales ( ^{Pompa-García et al., 2011}). Lo anterior se logra con el empleo de modelos de crecimiento o alométricos que a partir del d y at estiman el volumen (^{Liu et al., 2017}). Por ello, la construcción y el uso de modelos que describen las relaciones funcionales y la dinámica de crecimiento a lo largo del tiempo son de gran utilidad (^{Regalado et al., 2008}). Un modelo alométrico se define como una ecuación que describe cuantitativamente una relación que permite predecir una variable del árbol (at, volumen o biomasa) en función de otra, por ejemplo, el d (^{Picard et al., 2012}).

Las relaciones entre at y d, o alguna otra variable como el diámetro de tocón (dt), obedecen a una regla de proporcionalidad que es la misma para los árboles que crecen en condiciones semejantes (^{Bohlman y O’Brien, 2006}). Esto es un principio básico de la alometría que permite predecir variables difíciles de medir en función de otras de fácil medición (^{Huxley, 1924}).

La generación de ecuaciones confiables para algunas variables del árbol es compleja debido a la variabilidad que presentan y a las diferentes condiciones en que se desarrollan (^{Quiñonez et al., 2012}).

Desde hace algunos años, el ajuste de modelos de efectos mixtos (MEM) se ha convertido en una opción viable para este tipo de estudios (^{Corral-Rivas et al., 2014}), los cuales reemplazan a las técnicas clásicas de regresión. La inclusión de parámetros aleatorios y específicos por unidad de muestreo, posibilitan la modelación de la variabilidad de un fenómeno dado como el crecimiento en altura, entre las diferentes localizaciones o factores de agrupación de la variable de interés, y con ello, la varianza del error se reduce notablemente (^{Rijal et al., 2012}; ^{Seoane, 2014}). La técnica de MEM incluye parámetros fijos comunes a toda la población y otros aleatorios específicos por niveles de agrupación (^{Chenge, 2021}; ^{Raptis et al., 2021}). Al utilizar MEM se obtienen estimadores más eficientes, precisos y confiables que los parámetros fijos del modelo ajustados por MCNL, por lo tanto, se puede reducir la fuente de variación de la variable de interés (^{Carrero et al., 2008}).

Por la importancia ecológica y económica de M. brownei, es fundamental que este taxón se estudie con detalle para generar herramientas biométricas que se deriven del uso de técnicas modernas de modelización que garanticen la estimación, con la más alta precisión posible, de variables de interés como la altura total, por lo que bajo la hipótesis de que al incluir la técnica MEM y covariables de agrupación dentro de un modelo at-d, se contribuirá a mejorar la calidad de ajuste estadístico, y por tanto, la capacidad predictiva. Los objetivos del presente estudio fueron: 1) Evaluar el ajuste de modelos alométricos y de crecimiento bajo la técnica de MCNL para predecir la at en función del d en árboles de M. brownei para seleccionar al de mejor calidad de ajuste, y 2) Ajustar el modelo base seleccionado mediante MEM, bajo la combinación de dos factores de agrupación para que se incremente la capacidad predictiva de la relación at-d para árboles del taxón en Quintana Roo, México.

Materiales y Métodos

El estudio se realizó en el estado de Quintana Roo, México, en selvas medianas subperennifolias, medianas subcaducifolias, bajas subperennifolias y subcaducifolias que se desarrollan en lomeríos y en la llanura rocosa. El clima predominante es cálido subhúmedo (Aw), con temperatura media anual de 26 °C y precipitación media de 1 300 mm (^{Pennington y Sarukhán, 2005}; ^{Conafor, 2014}).

Se integró una base de datos conformada por 2 794 pares de datos de at-d de M. brownei, la cual se obtuvo de 172 conglomerados del Inventario Nacional Forestal y de Suelos 2004-2009 (^{Conafor, 2012}) (Figura 1).

Figura 1 Distribución espacial de los conglomerados por municipio con presencia de Metopium brownei (Jacq.) Urb. en Quintana Roo, México. 

Para efectos de aplicación práctica del modelo at-d por generar, la información se agrupó al combinar los niveles de clasificación de tipo de vegetación (Veg): selva mediada subcaducifolia (SMS), selva mediana subperennifolia (SMQ), selva baja subperennifolia (SBQ) y selva baja caducifolia (SBC), además de las agrupaciones de localidad como municipio (Mpio): Bacalar (Bac), Benito Juárez (BJ), Felipe Carrillo Puerto (FCP), Isla Mujeres (IM), José María Morelos (JMM), Lázaro Cárdenas (LC), Othón Pompeyo Blanco (OPB) y Solidaridad (Sol), ya que el aprovechamiento en Felipe Carrillo Puerto de esta especie para la elaboración de durmientes durante varias décadas eliminó los mejores individuos y probablemente afectó la morfometría de los árboles residuales.

El d con corteza se midió a 1.30 m del suelo, con cinta diamétrica graduada de cinco metros marca Qualitäs-bandmaβ de Forestry Suppliers Inc., graduada en centímetros y milímetros. La at se obtuvo desde la base del árbol hasta la punta de la copa con clinómetro Suunto^® graduado en grados y porcentajes.

La base de datos se dividió de manera aleatoria en 60 %, que se utilizó para realizar el ajuste estadístico de los modelos at-d y 40 % para validar el modelo seleccionado. Seis modelos no lineales se ajustaron, así como uno alométrico de tipo potencial (^{Burkhart y Tomé, 2012}; ^{Panik, 2014}; ^{Tuan et al., 2019}; ^{Hernández-Ramos et al., 2020}) (Cuadro 1). En una primera etapa, con el fin de obtener una ecuación base, el ajuste se realizó mediante el enfoque de mínimos cuadrados no lineales (MCNL) en el paquete estadístico R con el comando nls (^{R Core Team, 2016}).

La evaluación de bondad de ajuste y selección del modelo se realizó considerando la significancia de los parámetros (α≤0.05) y los mejores valores en los estadísticos coeficiente de determinación (R ² ), criterios de información de Akaike y Bayesiano (AIC y BIC) y logaritmo de la verosimilitud (loglik) (^{Corral-Rivas et al., 2014}; ^{Guerra-De la Cruz et al., 2019}). Después de seleccionar el mejor modelo, se ejecutó el MEM bajo la combinación de los factores de agrupación de Veg-Mpio.

La estructura matemática en forma matricial del modelo por ajustar bajo MEM, se basó en ^{Bronisz y Mehtätalo (2020)} y ^{Chenge et al. (2021)}, la cual para un modelo no lineal se expresó como:

Donde:

f = Función no lineal

Yij y Xij = i-ésima observación dependiente e independiente, respectivamente procedentes de la i-ésima unidad o nivel de clasificación

θij = Vector de parámetros r×1. r = Número de parámetros en el modelo y específico del j-ésimo grupo de clasificación, este vector puede dividirse para los parámetros fijos y aleatorios definidos como:

Donde:

Ai y Bi = Matrices de tamaño r×p y r×q para los efectos fijos y aleatorios, respectivamente, de cada nivel de clasificación

γ y bi = Vector p×1 y q×1 de los parámetros fijos y aleatorios. p×1 = Número de parámetros fijos. q×1 = Parámetros aleatorios

La inclusión de los efectos aleatorios (βi) se realizó de manera individual en cada uno de sus parámetros de los modelos de forma separada; por ejemplo:

Además, en el ajuste de MEM se incluyó una estructura de tipo varPower que representa una estructura de función de varianza de tipo potencia en la cual la distribución de los residuos se corrige por cada factor o nivel de agrupación. En este caso, los ajustes se realizaron mediante máxima verosimilitud en el paquete estadístico R con el comando nlme (^{Pinheiro y Bates, 2000}; ^{R Core Team, 2016}).

La selección del mejor modelo se hizo con los mismos criterios que los modelos de efectos fijos con MCNL. Asimismo, se verificaron los supuestos de regresión de la ecuación seleccionada mediante la comprobación gráfica de los supuestos de normalidad de los datos con la distribución de frecuencias y homocedasticidad de los residuales (^{R Core Team, 2016}). Las desviaciones globales del modelo seleccionado con MEM se verificaron con la raíz del cuadrado medio del error (RCME), mientras que la capacidad de predicción se confirmó a través del sesgo promedio de los residuales (E-, m) y la diferencia agregada expresada en porcentaje (DA %) (^{Lencinas y Mohr-Bell, 2007}; ^{Hernández-Ramos et al., 2020}).

El proceso de validación se realizó con dos muestras independientes: (i) con 40 % de la muestra separada de la base de datos original, y (ii) con 1 861 pares de datos at-d obtenidos en predios forestales en la zona de estudio. En ambos casos, se aplicó una prueba de t para comparar medias independientes (α=0.01) (^{Infante y Zarate, 2012}).

Resultados

El valor promedio para el d fue 17.88 cm, con valores mínimos y máximos de 7 y 80 cm, respectivamente; el coeficiente de variación fue de 49.36 % y presentaron una distribución leptocúrtica y asimétrica positiva. La at promedio fue 10.88 m, con valores extremos de 6 y 22 m, un coeficiente de variación de 27.58 %, y una distribución leptocúrtica con tendencia a ser simétrica (Cuadro 2).

En el ajuste de las ecuaciones locales de efectos fijos, se obtuvieron parámetros significativamente diferentes de cero a un nivel de significancia de 5 % (p≤0.05) en todos los casos. Los valores de R ² explicaron entre 39 % y 45 % de la variabilidad de los datos; los valores en los criterios de información Akaike y Bayesiano variaron entre 7 688 y 7 901, y el logaritmo de verosimilitud promedio fue de -3 857 (Cuadro 3).

Los mejores ajustes se obtuvieron para el modelo 4 que corresponde al de Chapman-Richards (Cuadro 3 y Figura 2). Al analizar las tendencias de las estimaciones, se observó que este modelo reproduce una trayectoria promedio similar a la de los datos observados (Figura 2) y presentó los valores más bajos del AIC y BIC, además del más alto del loglik, por lo tanto, fue seleccionado para ajustarlo mediante MEM y analizar el efecto de la inclusión de variables de clasificación. El modelo potencial sobreestima las dimensiones de altura a partir de las categorías mayores de 35 cm y los modelos de Strand, Hossfeld I Modificado y Verhulst-logístico subestiman la altura en las mismas categorías, por lo que fueron desestimados (Figura 2).

Figura 2 Distribución de datos observados y estimados por las ecuaciones de efectos fijos para los modelos at-d para Metopium brownei (Jacq.) Urb. en Quintana Roo, México. 

Cuando en el modelo se incorporó el componente aleatorio por el nivel de agrupación (Veg-Mpio), se observó que algunos de los parámetros no todos sus valores fueron significativos. Por ejemplo:

Donde:

βi: Inclusión del efecto aleatorio

En el Cuadro 4 solo se presentan los efectos que fueron mejores y significativos. Se determinó que al incluir el efecto mixto en el parámetro asintótico (β0+βi), se obtuvo el mejor comportamiento gráfico de los residuales, indicativo de que la modelación de la varianza es adecuada (Figura 3).

Figura 3 Distribución de residuales estandarizados con respecto a los valores predichos en los ajustes por la técnica de MCNL (a) y por la técnica de MEM (b) con el modelo de Chapman-Richards para Metopium brownei (Jacq.) Urb. en Quintana Roo, México. 

La inclusión de las covariables combinadas (Veg-Mpio) en el modelo, condujo a una ganancia estadística porque el valor de R ² mejoró 11.11 % con respecto al valor del modelo ajustado con MCNL. Los indicadores de verosimilitud (AIC, BIC y logLik), también mejoraron en promedio 2.5 % (Cuadro 4). Al contrastar la distribución de los residuales entre el ajuste con MCNL (Figura 3a) y la técnica MEM (Figura 3b), se obtuvo que esta última proporciona un patrón homocedástico con una distribución que se centra sobre el cero. La prueba de normalidad de Shapiro-Wilk mostró valores superiores a 0.91, lo que evidencia el cumplimiento de este supuesto de regresión.

Al evaluar la capacidad de ajuste por MEM, se determinó que el valor de la RCME fue de 2.36 m, el sesgo promedio de 0.003 y la DA fue de 0.0002 %. La validación por comparación estadística de medias en muestras independientes al 99 % de confiabilidad (α=0.01) evidenció que no existen diferencias significativas entre ellas. Para la primera muestra independiente se obtuvo t=0.582 y valor-p=0.5604, mientras que para la segunda muestra de validación, se produjo un valor de t=1.623 y valor-p=0.102, aspectos que pueden constatarse de manera gráfica en la Figura 4 (a y b).

Figura 4 Comparación de medias del proceso de validación del modelo mediante muestras independientes versus los datos estimados: (a) muestra compuesta por 40 % de las observaciones de la base original, y (b) muestra independiente de inventarios forestales. 

Debido a la robustez estadística y mejora de los MEM con respecto al ajuste por MCNL (Cuadros 3 y 4), así como su validación mediante muestras independientes (Figura 4), se obtuvieron los parámetros aleatorios específicos para cada nivel combinado de agrupación (Veg-Mpio) que podrán emplearse para cada condición específica de crecimiento (Cuadro 5).

Además, al contrastar tendencias extremas de at-d, se observó una diferencia promedio de 3.20 m en la altura, y en las categorías diamétricas de 35 y 40 cm (Figura 5).

Figura 5 Tendencia entre at-d en algunas condiciones de referencia específicas para Metopium brownei (Jacq.) Urb. en Quintana Roo, México. 

Discusión

La ganancia en la calidad de ajuste del modelo 6 obtenida bajo la técnica de MEM, se debe a que de acuerdo con ^{García y Rapelli (2011)}, al agrupar la información por alguna covariable, se reducen las desviaciones con respecto a los datos observados y el error de estimación, ya que se corrige la estructura de varianzas-covarianzas (^{Littell et al., 2006}) porque con esta técnica de análisis, la varianza es específica para cada nivel de agrupación (^{Rijal et al., 2012}; ^{Seoane, 2014}), y con ello se logra incrementar el R ² y mejorar los valores de otros estadísticos de bondad de ajuste (AIC, BIC y loglik).

La inclusión de un solo efecto aleatorio específico en los niveles de agrupación (Veg-Mpio) ofreció la más alta calidad predictiva cuando se asoció al parámetro asintótico, esto se explica porque es el atributo más variable y menos lineal (^{Tamarit-Urias et al., 2014}). Esa situación condujo a obtener estimadores más eficientes, precisos y confiables de los parámetros fijos del modelo y predecir parámetros aleatorios de cada unidad experimental que reflejan apropiadamente el patrón de desviación en relación a la media (^{De los Santos-Posadas et al., 2006}; ^{García-Espinoza et al., 2019}).

Por tanto, la expresión de Chapman-Richards ajustada bajo el enfoque de MEM es confiable para realizar estimaciones de at en función del d en árboles de M. brownei en los bosques naturales de Quintana Roo, México, tal como lo hizo ^{Arias (2004)} quien obtuvo un R ² de 0.64 para seis especies maderables de Costa Rica y ^{Castillo-Gallegos et al. (2018)} que citan valores de R ² =0.49 para Pinus chiapensis (Martínez) Andresen en plantaciones forestales de Tlapacoyan, Veracruz. En ambos estudios se utiliza la expresión de Chapman-Richards.

El modelo de Chapman-Richards también lo utilizaron ^{Saunders y Wagner (2008)} con éxito en una relación at-d bajo la técnica MEM e inclusión de covariables para nueve especies de árboles del noreste de los Estados Unidos de América. Dichos autores señalan que los modelos sigmoides son biológicamente más apropiados para estudiar las relaciones at-d. ^{Rijal et al. (2012)} y ^{Raptis et al. (2021)} también indican que, con la inclusión de covariables del rodal como predictores, se mejora la precisión de las estimaciones del modelo de Chapman-Richards al ajustarlo mediante MEM.

Es notable el efecto favorable que se logra cuando se usa la técnica de MEM en relaciones at-d. Así, ^{Hernández-Ramos et al. (2020)} obtuvieron un R ² de 0.46 y un incremento del valor referido a este estadístico a 0.70 con el modelo de Hossfeld IV bajo los enfoques de efectos fijos y aleatorios, respectivamente para Lysiloma latisiliquum (L.) Benth. en Quintana Roo, mientras que ^{García et al. (2017)} determinaron modelos locales de efectos fijos con R ² de 0.92 a 0.97 para predecir la at en función del d en ocho especies tropicales de Quintana Roo, México.

^{Hernández-Ramos et al. (2019)} ajustaron modelos at-d para árboles de Bucida buceras L. (pukté) bajo el enfoque de efectos mixtos. En ese estudio, el control de la variabilidad, al igual que para M. brownei, se realizó por agrupaciones al emplear el conglomerado, y obtuvieron una mejora estadística importante con el enfoque de MEM respecto al ajuste de parámetros fijos. Sin embargo, los autores referidos registraron un sesgo de -0.46 m, el cual es superior a lo estimado para M. brownei de 0.003 m, por ello, la ecuación propuesta es confiable para la estimación de la at, en función del d para esta especie en los bosques tropicales de Quintana Roo.

Al emplear el mejor modelo ajustado por MCNL se observa una tendencia de crecimiento constante (Figura 5), mientras que al aplicar los resultados obtenidos con el ajuste de MEM esta tendencia tiende a una asíntota con una diferencia entre ambos de 4.16 m, por lo que este enfoque contribuye a disminuir los errores de sobrestimación en la at utilizada en la elaboración de planes de manejo forestal para M. brownei, especie de interés comercial.

Conclusiones

La selección de un modelo base y su posterior ajuste mediante el uso combinado de la técnica de MEM con factores de agrupación dados por covariables, es un procedimiento eficiente para generar expresiones que permitan estimar la altura total en función del diámetro normal en árboles de Metopium brownei en Quintana Roo, México, por lo que puede adoptarse como una excelente estrategia para generar este tipo de herramientas biométricas.

El modelo de Chapman-Richards parametrizado a una relación at-d y ajustado mediante MEM con el factor combinado entre los tipos de vegetación y el municipio en el parámetro asintótico (β ₀ ), es el de mejor calidad predictiva porque aumenta la precisión de las predicciones, por lo tanto, se recomienda su uso en la realización de inventarios maderables de M. brownei en la región de estudio.

La aplicación de las ecuaciones obtenidas permitirá completar las bases de datos de inventarios forestales en las que por cuestiones logísticas de tiempo y económicas, no es posible la medición de las alturas de todos los árboles de la muestra.