Introducción

En México, la producción maderable proviene principalmente de bosques naturales, sin embargo, el país tiene déficit para abastecer la cantidad de productos forestales maderables que consume, lo que ha resultado que en la década del 2010 al 2020 se tuviera un déficit de 6.7 miles de millones de pesos anuales en la balanza comercial (^{SEMARNAT 2021}). Las plantaciones forestales tienen importancia creciente para satisfacer las necesidades de madera y subproductos, útiles en todo el mundo. En el país las principales entidades con superficies de plantaciones forestales maderables son: Tabasco (46 435 ha^-1), Veracruz (35 626 ha^-1), Campeche (26 608 ha^-1), Chiapas (16 960 ha^-1), Michoacán (16 076 ha^-1), Puebla (14 245 ha^-1) y Oaxaca (13 047 ha^-1) (^{Conafor 2020}). Debido a la creciente importancia de las plantaciones forestales (^{SEMARNAT 2021}), es necesario generar mecanismos que proporcionen alternativas de evaluación confiable y precisa, tanto para inventarios, como para la estimación de volumen, biomasa y carbono. En este sentido, los modelos predictivos representan esta alternativa, ya que ayudan a reducir costos y a disminuir el error asociado a las mediciones en campo.

Entre los modelos más usados en México se encuentran aquellos que se utilizan para predecir la altura de los árboles a partir del diámetro normal (d, cm) en combinación con variables del rodal, como el diámetro medio cuadrático y la altura dominante (^{Hernández et al. 2015}, ^{Corral et al. 2019}, ^{Camacho et al. 2022}). Estos modelos, comúnmente llamados generalizados, suelen ser útiles, ya que la altura es una variable difícil de medir en campo y susceptible a errores. Las dimensiones de d, medido a 1.30 m del ras de suelo, y la altura total (At, m) son los elementos más importantes para el manejo forestal (^{Ogana et al. 2020}). Ambas variables, se miden rutinariamente en inventarios forestales y son fundamentales para estimar el área basal, volumen total y la evaluación productiva de un sitio o rodal (^{Ogana 2019}).

El diámetro es la variable más precisa y sencilla de medir en campo, ya que se toma directamente del árbol, mientras que, medir la At es costoso y requiere de mayor tiempo, especialmente en rodales densos y terrenos poco accesibles (^{Ozcelik et al. 2018}). Debido al problema asociado con la medición de la At, generalmente sólo se mide una submuestra de árboles (^{Seki y Sakici 2022}). Para completar las mediciones, las alturas de los árboles restantes se estiman a través de modelos de At-d, que se elaboran a partir de numerosos pares de datos At-d obtenidos en campo (^{Kalbi et al. 2018}).

La relación At-d puede ser modelada con ecuaciones locales, aquellas que solo utilizan variables del árbol, como predictoras, como el d, sin embargo, las ecuaciones locales tienen un intervalo de aplicación reducido a nivel de individuos o sitios específicos (^{Hernández-Ramos et al. 2018a}). Mientras que ^{Corral et al. (2019)}, han señalado que modelar las relaciones At-d con una sola variable, por ejemplo, el d, puede no ser adaptable a diferentes dinámicas de rodal y condiciones silvícolas. Por lo tanto, es conveniente desarrollar modelos que incluyan la relación At-d de manera generalizada, que incorporen variables del rodal, para tener en cuenta diferentes variaciones en la estructura del rodal y condiciones del sitio (^{Ogana et al. 2020}), junto con un mayor intervalo de aplicabilidad para un bosque, plantaciones forestales o rodal, donde las áreas presentan distintas condiciones de crecimiento (^{Hernández-Ramos et al. 2015}).

Un modelo generalizado At-d permite a los usuarios y técnicos forestales obtener estimaciones con menor error al estimar la At a partir de mediciones de d (^{Hernández-Ramos 2018b}, ^{Bronisz 2020}, ^{Camacho 2022}). Los modelos At-d son uno de los componentes más importantes en la ciencia y manejo forestal, ya que una de sus funciones principales es la determinación y evaluación del volumen comercializable, crecimiento, rendimiento y condiciones del sitio de los ecosistemas forestales (^{Seki y Sakici 2022}). En México, aunque se han tenido avances importantes en la generación de modelos alométricos, la diversidad de especies y condiciones ambientales incrementan la necesidad de desarrollar nuevos modelos con menos errores. Tal es el caso del estado de Michoacán, en dónde Pinus leiophylla Schl. & Cham. es una de las especies más importantes, económica y ecológicamente, en particular dentro de la Comunidad Indígena (CI) de Patamban y en la región forestal de la Meseta Purépecha, Michoacán, por lo cual el presente trabajo tiene como objetivo determinar un modelo generalizado altura total-diámetro normal para plantaciones forestales de P. leiophylla, en la CI de Patamban, Tangancícuaro, Michoacán, México.

Materiales y métodos

El estudio se realizó en la Meseta Purépecha ubicada entre los 19°45´00´´de Latitud norte y los 102°20´00´´ de Longitud oeste, en la parte centro-noroccidental del estado de Michoacán. Con altitud que oscila entre 1 700 y 3 200 m s. n. m. El clima predominante es templado húmedo y tropical con lluvias en verano. La temperatura media anual varía entre 14-23 °C, con precipitaciones en el intervalo de 600 a 1 100 mm, durante los meses de junio a octubre. El relieve de la región está conformado por pequeños accidentes orográficos constituidos por el Eje Neovolcánico Transversal; aunque predominan algunos relieves de valles y planicies entre los cerros y mesetas. La vegetación arbórea dominante de esta zona está compuesta principalmente por especies de Pinus, Quercus y Abies (^{IPLAEM 2021}, ^{Orduña-Trejo et al. 1999}).

En siete plantaciones forestales, que en total tienen 22 ha, se establecieron 33 sitios cuadrangulares temporales de 20 m x 20 m (400 m²) distribuidos completamente al azar dentro de plantaciones sanas de P. leiophylla. Los árboles que se plantaron provenían de viveros comunitarios y viveros estatales de Conafor. La densidad inicial de plantación fue de 2 500 árboles ha^-1, en las que se aplicaron aclareos al 30% cuando tenían tres años de establecidas. En cada plantación de entre ocho y 28 años, se obtuvieron datos dasométricos de todos los árboles en cada uno de los 33 sitios que incluyeron mediciones del d a 1.30 m con cinta diamétrica de 5 m modelo Forestry Suppliers® y con el Hipsómetro Nikon ® Forestry PRO II con precisión de cm se obtuvo la At. Los datos obtenidos de campo fueron capturados en Microsoft Excel® para integrar la base de datos.

En total se utilizaron para el ajuste estadístico 779 pares de datos At-d provenientes de 25 sitios seleccionados con supervivencia inicial en los primeros tres años mayor o igual que 75% (1 875 árboles ha^-1) y que en promedio tenían una tasa de mortalidad anual del 6.14%, en edad de ocho a 28 años, At de 1.86 m a 19.60 m y d entre 1.5 cm y 36 cm; para el ajuste de los modelos generalizados a través de la función nls del programa R (^{R Core Team 2018}) por mínimos cuadrados no lineales (MCNL). Los modelos seleccionados se tomaron de ^{Hernández et al. (2015)}, ^{Hernández-Ramos et al. (2018b)}, ^{Corral et al. (2019)} y ^{Ogana et al. (2020)} (Tabla 1).

La selección del mejor modelo fue a través de la significancia del valor de los parámetros (a=0.05), los valores del coeficiente de determinación ajustado (R², [8]), la raíz del cuadrado medio del error (RCME, [9]), criterio de información de Akaike y Bayesiano (AIC [10] y BIC [11]); y del Sesgo [12] (^{Bronisz y Mehtätalo 2020}).

Donde, yi, yi^ y y-i son los valores observados, estimados y promedio, respectivamente; 𝑛 es el número total de datos utilizados en el ajuste de los modelos; 𝑝 se refiere al número de parámetros (Tabla 1).

Una vez elegida la estructura generalizada At-d a través de mínimos cuadrados no lineales (MCNL), se reajustó la expresión a través de la función nlme con el enfoque de modelo de efectos mixtos (MEM) con la finalidad de estimar la variación dentro y entre sitios de muestreo (^{Corral et al. 2019}, ^{Ogana et al. 2020}). Este último enfoque es más verosímil, tiene una mayor explicación de la muestra y un menor error de estimación en contraste con MCNL, ya que agrupa la variabilidad por nivel o factor (^{Bronisz y Mehtätalo 2020}, ^{Pinheiro y Bates 2023}). La estimación de los parámetros se llevó a cabo a través del método FOCE (First Conditional Expectation, por sus siglas en inglés) (^{Lindstrom y Bates 1990}). La inclusión del efecto aleatorio se realizó de forma aditiva (+a_i) e individual en cada uno de los parámetros de la forma siguiente al tomar de referencia el modelo generalizado de Wang y Tang (6):

Dónde la forma matricial es equivalente a: Aij=xij;θi+εij . Atij y xij se refiere al valor de la altura y las variables de diámetro ( dij ) o inherentes al rodal (e. i. diámetro cuadrático, Dqij ) del j-ésimo individuo de cada sitio (i); 𝜃 𝑖𝑗 es igual a Aiγ+βibi donde Ai y βi es la matriz de dimensión r×p y r×q para el número de parámetros ( r ) fijos ( p ) y aleatorios ( q ) del modelo por sitio, respectivamente; mientras que, γ y bi son los vectores de cada uno; y εij es el error aleatorio asociado al ajuste.

Debido a que es común encontrar heterogeneidad de varianzas en los datos (^{Lindstrom y Bates 1990}, ^{Pinheiro y Bates 2013}), se incluyó en la modelación una estructura de potencia de la covariable o valores predichos: Varε1=σ2∙di2∙ϑ [13], donde 𝑑 𝑖 se refiere al valor de diámetro de cada individuo (i) y ϑ es el valor de potencia en la corrección de la heterogeneidad de varianzas en los residuos (^{Garibaldi et al. 2019}, ^{Pinheiro y Bates 2023}). La selección de la mejor expresión se realizó de la misma forma que con MCNL al tomar la significancia de los parámetros y de referencia los valores del R², AIC, BIC, RCME y Sesgo como indicadores de bondad de ajuste y de robustez del modelo, los cuales han sido utilizados en estudios semejantes (^{Corral et al. 2019}, ^{Camacho et al. 2022}). Se realizó un contraste entre la distribución de los residuales de los ajustes de MCNL y MEM para visualizar la corrección de heterogeneidad de varianzas propuesta, además de forma gráfica se verificó la frecuencia de los residuos los cuales deben de tener una forma de campana de Gauss para asumir normalidad en el ajuste.

La calibración o localización del modelo de efectos mixtos se realizó al utilizar sólo ocho sitios de muestreo con supervivencia inicial a los tres años de establecidas menor que 75% (125 pares de datos At-d) de acuerdo con la metodología propuesta por ^{Vonesh y Chinchilli (1997)} y aplicada por ^{Corral et al. (2019)} y ^{Ogana et al. (2020)} de tres formas: i) utilizar el valor del d mínimo de la parcela, ii) emplear los dos diámetros más pequeños del sitio, iii) tomar de referencia los dos valores extremos de la unidad de muestreo, y iv) el diámetro mínimo y la dimensión de los dos árboles más gruesos de la parcela. Para este procedimiento se seleccionaron seis parcelas independientes al ajuste, las cuales fueron evaluados sus resultados de acuerdo con los valores de la RCME, Sesgo y error porcentual absoluto medio (E%AM: 100∙∑yi-yi^yin , [14]).

En este paso se utilizó el algoritmo Lindstrom-Bates a través del método FOCE, que tiene la siguiente estructura para determinar el valor del parámetro aleatorio por sitio de calibración en los parámetros fijos del modelo:

Dónde: D^ , es la matríz de varianza-covarianza q×q de los parámetros aleatorios incluidos en el modelo ( q ), Ri^ es la matriz de varianza-covarianza 𝑚 𝑗 × 𝑚 𝑗 del termino del error y Z^i es la matriz de las derivadas parciales m×q de los parámetros aleatorios evaluados en b^=0 ; y ei^ se refiere al vector de residuos m×1 donde los componentes se obtienen por la diferencia entre el valor observado y predicho.

Resultados

En la estadística descriptiva de la muestra utilizada para el ajuste estadístico de los modelos At-d generalizados no se observan problemas de curtosis y asimetría en ninguna de las variables utilizadas (valor <3.0 unidades), pero es evidente la variabilidad muestral alta en las variables de d, dp, N y Dq, así como valores contrastantes entre la At y Ap (Tabla 2). Además, la relación alométrica promedio de la especie indica que, por cada cm de incremento en las dimensiones del d, la At aumentará en promedio 0.46 m (índice de esbeltez).

En el ajuste por MCNL, los modelos de la tabla 1 SAI (1), Mirkovich (2) y Cox III (5) presentaron parámetros no significativos por lo que fueron descartados para su posterior análisis, resultado probablemente a que su estructura matemática no alcanza a modelar la variabilidad muestral dada por las variables explicativas utilizadas (d, Dq y N). Esto contrasta con las dos expresiones que utilizan la Ad como variable inherente a la población para explicar la At del sitio; Michaelis-Menten (4) y Wang y Tang (6), donde esta última es la que mejor indicadores de ajuste tiene y mayor robustez estadística muestra (R², AIC, BIC, RCME y Sesgo) por lo que fue seleccionada como la mejor (Tabla 3).

En un segundo enfoque de ajuste del modelo de Wang y Tang, se observa que en los MEM a medida que se tiene una mejor corrección de la heteroscedasticidad en la distribución de los residuales; valor más cercano a la unidad en el parámetro relacionado a la función de varianza (ω), los valores de R ² incrementan, se hace más verosímil el modelo (AIC y BIC) y disminuyen las desviaciones (RCME, Sesgo y Residual) con respecto al enfoque de MCNL (Tabla 4). Por tal motivo, se eligió la estructura de MEM con la inclusión del efecto aditivo en el parámetro a₀ como la mejor.

En la Figura 1 se muestra la mejora en la distribución de lo residuales entre el ajuste por MCNL (a) y al incluir la estructura de tipo potencia en los MEM (c) en donde se puede observar que este procedimiento elimina la heteroscedasticidad común de este tipo de ajustes At-d. La prueba gráfica en la frecuencia de los residuos señala el cumplimiento de la normalidad (Figura 1b), mientras que, el valor aleatorio del parámetro 𝑎 2 muestras variabilidad entre sitios (Figura 1d).

Al combinar la derivada parcial del modelo de Wang y Tang (16: 0=-a0∙Ada1∙e-a2d∙a1∙e-a2d∙lnAdd ) con los valores de la matriz de varianza-covarianza (Tabla 4) y aplicar las cuatro propuestas de calibración, se obtiene que al emplear el diámetro mínimo y la dimensión de los dos árboles más gruesos de la parcela es donde se obtienen la mayor precisión de las estimaciones al emplear los MEM con una muestra independiente del ajuste (Tabla 5). En la Figura 2 también se observa que al agrupar las desviaciones por categoría diamétrica de las distintas calibraciones se obtuvo que el sesgo por estimación es más consistente que las demás formas de localización de los parámetros.

Figura 1 Distribución de residuales del modelo de Wang y Tang obtenidos por mínimos cuadrados no lineales (a: MCNL) y modelos de efectos mixtos (b: MEM); prueba de normalidad (c) y distribución de los valores aleatorios del parámetro a₂ (d). 

Figura 2 Desviaciones por categoría diamétrica (CD) por método de calibración. 

Discusión

La variabilidad de los datos se atribuye a la edad múltiple de las plantaciones y calidad de estación en donde se establecieron, la densidad inicial y actual, ya que por hábitos de crecimiento la especie es tolerante a la sombra y no se presenta la muerte por competencia de manera temprana, además de que no todas las masas fueron intervenidas con un aclareo de 30% de forma homogénea y en el tiempo requerido (Turno técnico en altura: Incremento corriente anual = Incremento medio anual) por el índice de sitios diferenciado de cada una de ellas.

El uso de modelos generalizados de At-d contribuye a incluir las distintas dinámicas de crecimiento de los rodales dentro de la modelación o la respuesta diferenciada de los tratamientos silvícolas aplicados, lo que no pueden ser representados dentro de los modelos locales At-d (^{Corral et al. 2019}). La inclusión de la Ad dentro de la estructura mejora los resultados y concuerda con lo mencionado por ^{Uzoh (2017)} donde esta variable es un indicador indirecto de la calidad de estación del sitio y es la menos afectada por la densidad; o bien puede ser sustituida por las etiquetas de índice de sitio al combinar la información con la edad de cada plantación (^{Mensah et al. 2018}).

La selección del mejor modelo difiere con lo señalado por ^{Corral et al. (2019)} donde mencionan que solo incluir la Ad en la relación At-d no es suficiente para describir la variabilidad dada por las condiciones del sitio, pero en este trabajo combinar la inclusión del efecto aleatorio y la calibración se consideran suficientes para cubrir las diferencias en la relación alométrica At-d para los árboles establecidos como plantaciones forestales de Pinus leiophylla. Además, contrasta con lo referido por ^{Bronisz y Mehtätalo (2020)} los cuales señalan que se deberá de priorizar la inclusión de las variables inherentes al rodal que no tengan un costo adicional de tiempo y esfuerzo dentro de los modelos generalizados At-d como el dp, Ab, Dq o la N. Sin embargo, en este trabajo los modelos que incluían estas variables mostraban un 5% de menor explicación de la variabilidad muestral y valores de RCME mayores en un 26.6% con respecto al modelo de Wang y Tang que incluye la Ad, además, en México, medir la Ad es una variable común, lo que justifica su inclusión en el modelo, ya que no implicaría más esfuerzo del que ya se hace dentro de un inventario forestal tradicional.

Los MEM mostraron una mejora tanto en la explicación de la variabilidad muestral (R²: 2.4%), estadísticos de ajuste (AIC y BIC de 24.5% en promedio) y en la precisión de las estimaciones con respecto al enfoque de MCNL (RCME y Sesgo: >14%), por ser una alternativa viable para representar satisfactoriamente las diferencias en la relación alométrica dadas por las distintas condiciones de crecimiento (^{Pinheiro y Bates 2013}, ^{Bronisz y Mehtätalo 2020}). Además, este enfoque ayuda a superar el problema de dependencia de los datos, debido a las estructuras jerárquicas intrínsecas de los inventarios forestales.

La elección de la muestra empleada para la calibración del modelo concuerda con lo señalado por ^{Bronisz y Mehtätalo (2020)} los cuales al evaluar empíricamente las estrategias de selección indican que los datos extremos del sitio son donde se obtienen los menores errores residuales y es más eficiente en la localización del parámetro aleatorio. Además, este enfoque permitió ajustar satisfactoriamente el valor de At estimado por unidad de muestreo, lo que permite ampliar la aplicabilidad de la expresión matemática seleccionada tal y como lo reportan ^{Burkhart y Tomé (2012)}, ^{Uzoh (2017)} y ^{Camacho et al. (2022)} al emplear este enfoque de localización en distintos rodales de especies forestales.

Conclusiones

El modelo generalizado altura total (At)-diámetro normal (d) de Wang y Tang que incluye como variable la altura dominante del sitio (Ad) fue eficiente para representar la variabilidad muestral de las plantaciones forestales comerciales de Pinus leiophylla Schl. & Cham. establecidas en la Comunidad estudiada. El enfoque de modelos de efectos mixtos (MEM) mostró mejor precisión en las estimaciones de la At mediante la reducción del sesgo implícito en la aplicación de modelo generalizado At-d de Wang y Tang. También fue eficiente en agrupar la variabilidad muestral por nivel de clasificación (Sitio) en la información empleada para el ajuste con respecto a mínimos cuadrados no lineales. La corrección de heterogeneidad de varianzas en los datos al incluir una estructura de potencia de la covariable en la modelación volvió más homocedástica los residuales del modelo propuesto. La calibración mediante la inclusión de los valores del diámetro mínimo y la dimensión de los dos árboles más gruesos de la parcela señaló que ajusta el valor del parámetro aleatorio y es donde se obtiene la mayor precisión de las estimaciones al emplear los MEM de una muestra independiente del ajuste. Esto permite ampliar su aplicabilidad del ajuste propuesto.