Introducción

Con la descripción del diámetro a distintas alturas sobre el fuste (ahusamiento: di, cm) se puede detallar la información de un inventario alcuantificar la distribución de productos y ajustar las existencias maderables para un mercado diferenciado (^{Rachid et al., 2014}). En la literatura existen diversas expresiones matemáticas que contribuyen a describir el di. En México se han utilizado modelos de manera satisfactoria en distintas especies y condiciones de crecimiento, los de tipo polinómico de orden n (^{Hernández-Ramos et al., 2018}; ^{Ramírez-Martínez et al., 2018}), exponenciales (^{Pompa et al., 2009}) y segmentados (^{Hernández-Ramos et al., 2017}; ^{Tamarit et al., 2014}). De igual manera, sistemas compatibles de ahusamiento y volumen comercial (^{Cruz-Cobos et al., 2008}; ^{García-Espinoza et al., 2019}) o de razón de volumen (^{Cruz-Cobos et al., 2023}).

Generalmente, el ajuste estadístico se hace mediante mínimos cuadrados ordinarios (MCO) (^{Rachid et al., 2014}), ecuaciones aparentemente no correlacionadas (SUR) (^{Flores et al., 2021}), máxima verosimilitud (FILM) (^{Hernández-Ramos et al., 2017}), mínimos cuadrados generalizados no lineales (MCGNL) (^{Monárrez-González et al., 2024}) y algoritmo de minimización de la suma de cuadrados de Marquardt (^{Niño et al., 2018}). Sin embargo, los modelos de efectos mixtos (MEM) (^{Tamarit et al., 2014}) o la inclusión de variables Dummy, pueden mejorar estadísticamente los resultados (^{Torres et al., 2020}).

En los ajustes de funciones de perfil del fuste o de sistema de ecuaciones simultáneas volumen comercial-ahusamiento se utilizan datos del mismo árbol; entonces, los componentes del error están correlacionados y el ajuste por MCO producirá estimadores teóricos robustos, pero ignorando los errores (^{Cruz-Cobos et al., 2008}; ^{Hernández et al., 2013}). Para solventar problemas de multicolinealidad, la técnica de MEM es una opción viable (^{Pinheiro & Bates, 2000}) en cualquiera de los tres enfoques de su aplicación: (I) Predicción de la variable de interés con mayor certeza mediante una mejora estadística (^{Cruz-Cobos et al., 2008}; ^{Tamarit et al., 2014}); (II) Estimación de los componentes de varianza para ampliar su aplicabilidad o reducir el esfuerzo de muestreo (^{Saygili & Kahriman, 2023}); y (III) La comprensión del fenómeno mediante la explicación de la interacción genotipo-ambiente (^{Balzarini, 2002}; ^{Bandera & Pérez, 2018}).

Los MEM incorporan efectos aleatorios a sus parámetros, lo que influye favorablemente en el término de error (^{Correa & Salazar, 2016}; ^{Pinheiro et al., 2025}), ya que corrigen la estructura de varianzas-covarianzas asociadas con datos lineales o de remediciones de los árboles, y se ha recurrido al uso de esta técnica para modelar los patrones de volumen y crecimiento de árboles, lo cual ha generado mejores resultados, con respecto a MCO (^{Zuur et al., 2009}). En Michoacán, las plantaciones forestales (PF) son una opción para reducir la presión sobre los bosques ante la demanda de madera por la sociedad. Por lo anterior, existe la necesidad de contar con herramientas silvícolas cuantitativas acordes a las condiciones de cada especie para la gestión de los recursos, es así como se planteó el objetivo de ajustar una función de perfil fustal y definir su volumen fustal para árboles de Pinus leiophylla Schiede ex Schltdl. & Cham. en la comunidad indígena Patamban, municipio Tangancícuaro, Michoacán, México.

Materiales y Métodos

El estudio se realizó en la comunidad indígena Patambán, Tangancícuaro, Michoacán, México, ubicada en la Meseta Tarasca que se encuentra en la región fisiográfica del Eje Neovolcánico como parte del sistema de sierra volcánica compleja y pequeños valles montañosos (pendientes 20-65 %) entre 1 700 y 3 500 m de altitud. Esta comunidad se localiza a una altitud de 1 740 m en las coordenadas geográficas 19°53’45.80” LN y 102°12’51.20” LO, en un clima templado subhúmedo con lluvias en verano y temperatura promedio de 16 a 18 °C, así como un suelo de tipo Andosol (^{Instituto Nacional de Estadística y Geografía [INEGI], 2010}).

Se aplicó un muestreo aleatorio de 12 sitios cuadrangulares temporales de 20×20 m (400 m²), en cuatro plantaciones de P. leiophylla, con edades de 8 a 28 años. En cada sitio se recopilaron datos de densidad del arbolado (árboles ha^-1) y se eligieron 36 individuos dominantes en el sitio (tres árboles por sitio), a los cuales se le midió, de forma directa, el diámetro normal a la altura de 1.30 m sobre el suelo (d, m) con una cinta diamétrica marca Forestry Suppliers Inc. ^® modelo P. O. BOX JACKSON, además del diámetro a distintas alturas sobre el fuste (di) a 0.3, 0.7, 1.30 m y a cada 2.5 m de altura iniciando en el diámetro y altura del tocón. De manera indirecta, con el Telerrelascopio Bitterlich ^® modelo SW, hasta llegar a la altura total (At, m), y desde el diámetro del tocón (di=0) se determinaron las dimensiones de diámetro para posteriormente realizar los cálculos de ahusamiento a distintas alturas sobre el fuste de cada individuo (Ai, m).

Se calculó el volumen por sección (V _sección , m³) mediante la fórmula de Smalian (Ecuación 1) y el volumen de la punta (V _punta , m³) con la expresión del cono (Ecuación 2) (^{Niño et al., 2018}). Estos volúmenes (V _sección y V _punta ) se sumaron para obtener el volumen fustal del árbol (Vf, m³).

Donde:

V _sección = Volumen de la troza (m³)

V _punta = Volumen de la punta (m³)

g _n-1 = Área basal del diámetro mayor de la troza (m²)

g _i = Área basal del diámetro menor de la troza (m²)

g _n = Área basal del diámetro mayor de la punta (m²)

L = Longitud de la troza (m)

Con 245 pares de datos de di y Ai se ajustaron ocho expresiones de ahusamiento (^{Pompa et al., 2009}; ^{Ramírez-Martínez et al., 2018}; ^{Torres et al., 2020}) en el programa Rstudio ^® 2024.04.2 versión Build 764, en un primer enfoque a través de mínimos cuadrados no lineales (MCNL) con la función nls (Cuadro 1) (^{Baty et al., 2015}; ^{R Core Team, 2024}). Para evitar problemas de convergencia estadística en los ajustes, se incluyó un valor de delta=0.01 en di para impedir pérdida de observaciones y el valor de cero en At (^{Hernández et al., 2013}).

La selección del modelo fue en función de la significancia del valor de los parámetros (a=0.05), el coeficiente de determinación (R ² , Ecuación 11), la raíz del cuadrado medio del error (RCME, Ecuación 12), criterio de información de Akaike (AIC, Ecuación 13) y Bayesiano (BIC, Ecuación 14), logaritmo natural de la función de verosimilitud ( Fθ, Ecuación 15) y del sesgo (Ecuación 16) (^{Bronisz & Mehtätalo, 2020}; ^{Hernández et al., 2013}).

Donde:

yi, yi^ y y-i = Valores observados, estimados y promedio, respectivamente

n = Número total de datos utilizados en el ajuste de los modelos

p = Número de parámetros (Cuadro 1)

fi = Función de modelo

θ = Probabilidad de un vector de parámetro

ln = Logaritmo natural

Una vez seleccionada la expresión base, esta se ajustó con la técnica de MEM con la función nlme por máxima verosimilitud (^{Pinheiro & Bates, 2000}; ^{Pinheiro et al., 2025}; ^{R Core Team, 2024}) y el método de expansión de Primer Orden (FO, First-Order por sus siglas en inglés) (^{Fu et al., 2014}; ^{Yang & Huang, 2013}), en el cual, el nivel de agrupación por individuo se incluyó de forma aditiva (+ui) al valor de uno o la combinación de varios de sus parámetros fijos (βi) (^{Zuur et al., 2009}). Por ejemplo, para la ecuación de Clutter: di=β1dβ2At-Ai(β3+ui)Atβ4+ε (Ecuación 17), y en donde el parámetro seleccionado para la inclusión de los efectos aleatorios es β3 deberá de incluir los mejores valores de los estadísticos de AIC, BIC y logLik, que la RCME y Sesgo sean los más bajos, y que todos sus parámetros fuesen significativos (a=0.05).

La corrección por heterocedasticidad se realizó a través de la expresión varExp: varεij=exp2∙δ1∙vi (Ecuación 18), donde varεij es la función de varianza evaluada en la covariable de varianza de los residuales del predictor (vi); mientras que δ1 y δ2 se refieren a los coeficientes de la función de varianza, los cuales serán específicos para cada nivel ( δi) (^{Bronisz & Mehtätalo, 2020}; ^{Pinheiro & Bates, 2000}). Los problemas de medidas repetidas en un mismo individuo se abordaron a través de una estructura de correlación de media móvil de autocorrelación (ARMA) de orden (p, q), donde p y q son los números enteros no negativos que especifican respectivamente el orden autorregresivo y de promedio móvil de la estructura ARMA, en cuyo caso ambos tienen el valor predeterminado de 0 (^{Pinheiro et al., 2025}; ^{Zuur et al., 2009}).

El cumplimiento de los supuestos de regresión se llevó a cabo mediante las pruebas gráficas de normalidad en la frecuencia de los residuos, distribución homocedástica de los residuales y la autocorrelación de los errores (^{Martínez-González et al., 2014}). De igual manera, se verificó la distribución de los valores aleatorios de los parámetros por árbol obtenida en el ajuste de MEM.

Se determinó la ecuación de Vf implícita de todas las expresiones de di (Cuadro 1) a través del valor de sus parámetros, para posteriormente integrar el valor del d a lo largo del fuste con respecto a Ai como un sólido en revolución; para el Vt se utilizó la Expresión 19 (^{Flores et al., 2021}; ^{Hernández-Ramos et al., 2018}; ^{Pompa et al., 2009}; ^{Ramírez-Martínez et al., 2018}):

Donde:

Vt = Volumen fustal

k = Constante volumétrica

δ = Diferencial sobre la altura del fuste

A ₁ y A ₂ = Alturas que definen la integral y que pueden variar de A ₁ =0 hasta A ₂ =At

d = Diámetro normal (m)

Las estimaciones de Vf con los modelos implícitos se contrastaron a través de un análisis de los valores observados como medias independientes con una prueba de t a un p=0.05 (^{Infante & Zarate, 2012}). Las hipótesis fueron: H₀: no existe diferencia entre las estimaciones (µ ₁ =µ ₂ ), y H_a: el valor real de la media poblacional (Vf) es distinto del valor que establece la H₀ (µ ₁ ≠µ ₂ ).

Resultados

El ajuste estadístico mostró parámetros no significativos (a<0.05) en las expresiones de Cielito 1, Rentería y Rustagi-Loveless; las ecuaciones de Amidon y Newnham solo explicaron 42.6 y 78.9 % de la variabilidad muestral, respectivamente; mientras que Clutter, González y Forslund tuvieron un valor de R ² mayor del 0.934. La expresión de Clutter presentó los menores valores en AIC, BIC y logLik, así como las más pequeñas desviaciones de RCME y Sesgo; por lo cual fue seleccionada como la ecuación base para desarrollar el modelo compatible de ahusamiento y volumen (Cuadro 2).

Al incluir el efecto aleatorio de manera aditiva por individuo, como un nivel de agrupación en la expresión de Clutter, resultó que al considerar el efecto en dos o más parámetros, los valores no fueron significativos (a<0.05); mientras que, al utilizar los parámetros β1, β2 y β4 los valores de R ² , AIC, BIC, logLik, RCME y Sesgo no mejoraron las estimaciones en el di; en contraste de cuando se utilizó β3, el cual está relacionado con la diferencia entre At y Ai de cada árbol. Además, el ajuste por MEM mejoró en 2.28 % los valores de verosimilitud, con respecto al ajuste a través de MCNL; mientras que el Sesgo se redujo alrededor de 56.8 % (Cuadro 3).

Al evaluar de forma gráfica los supuestos de regresión, se observó una distribución en la frecuencia de los residuales con forma de campana de Gauss (Figura 1A ), y residuos homocedásticos alrededor de cero (Figura 1B), lo cual indicó el cumplimiento de normalidad y la falta de heterocedasticidad en los residuos.

A = Normalidad; B = Homocedasticidad.

Figura 1 Pruebas gráficas del modelo de ahusamiento Clutter ajustado por modelos de efectos mixtos. 

Los valores de autocorrelación parcial fueron menores a 0.2 (Figura 2A) y señalaron el cumplimiento del supuesto de regresión en los errores obtenidos del ajuste. Además, se graficó el valor del efecto aleatorio por individuo dentro del parámetro fijo seleccionado (β3), lo que demostró la variabilidad del di en la muestra utilizada (Figura 2B).

A = Autocorrelación de los residuales; B = Distribución de los valores aleatorios de los parámetros árbol por árbol.

Figura 2 Pruebas gráficas del modelo de ahusamiento Clutter ajustado por modelos de efectos mixtos. 

Con el objetivo de ampliar la aplicabilidad del ajuste presentado, se determinaron los valores de la matriz de varianza-covarianza (vcov) mediante los MEM (Cuadro 4), con la idea de una posible calibración con datos independientes a la muestra del parámetro relacionado a la diferencia entre At y Ai (β3), cuando se trata de evaluar la variabilidad del di en individuos con características distintas a la información empleada para el ajuste o fuera del área geográfica donde se desarrollan las PF de la especie en estudio.

Con la finalidad de identificar los modelos implícitos de Vf de cada expresión de di, se obtuvieron formas matemáticas para estimar el Vf de tipo Schumacher-Hall para los modelos de Clutter, Cielito 1 y Rentería; González, Amidon y Newnham en su estructura incluyen la forma de variable combinada (d2×At), pero en lo global la expresión matemática implícita es poco aplicable; los modelos Forslund y Rustagi-Loveless no tienen una solución matemática explícita, por lo cual no se incluyen en el Cuadro 5.

Al contrastar las estimaciones con estos modelos implícitos y el Vf observado se tuvieron diferencias significativas (DMS=0.1313, F=84.36; error=0.0367, p<0.0001). De acuerdo con Tukey (a=0.05), las estimaciones formaron dos grupos: uno en el cual los modelos implícitos de Clutter, Renteria, Amidon y Newnham, con respecto de los datos observados son iguales, y otro con el de Gonzáles. Sin embargo, los parámetros de la expresión de Schumacher-Hall (Ecuación 17) derivada de la función di de Clutter (Ecuación 3) son compatibles y es posible llegar a una expresión de volumen comercial (Vc), si así se desea.

La Figura 3 muestra una tendencia lineal entre las estimaciones de ahusamiento realizadas con la ecuación de Clutter versus los valores observados y una descripción acorde con la forma de los árboles, por lo que se confirma la explicación de la tendencia de la información al describir el perfil de los árboles establecidos en plantaciones de Pinus leiophylla.

A = Para árboles de Pinus leiophylla Schiede ex Schltdl. & Cham.; B = Descripción del ahusamiento versus estimaciones de di a diferentes Ai con la ecuación de Clutter.

Figura 3 Comparación del diámetro observado (di) a distintas secciones del fuste (Ai) contra el predicho (ahusamiento). 

Discusión

Aun cuando la presencia de varianza no constante en la distribución de residuales (heterocedasticidad) es común al relacionar variables biológicas, no invalida los estimadores obtenidos por MCNL, ya que continúan siendo insesgados y dejan de tener una varianza mínima (^{Hernández et al., 2013}); el enfoque de MEM aplicado en el presente estudio fue viable para la corrección de ese problema al modelar el di. En este sentido, ^{Cruz-Cobos et al. (2008)} al utilizar los MEM, flexibilizaron un modelo lineal polinómico e incrementaron la capacidad predictiva de la expresión di de Cielito 1 para árboles de Pinus cooperi C. E. Blanco en Durango, México.

^{Tamarit et al. (2014)} utilizaron los MEM e incluyeron el efecto aleatorio dentro de un sistema compatible de di-Vc, con ello lograron una mejora estadística con respecto a MCNL, porque se propició un control de la variabilidad específica por individuo; además, se mejoró significativamente el valor del error estándar de los parámetros del modelo y se redujo el Sesgo individual por estimación debido a que la técnica MEM agrupa la variabilidad de forma específica por nivel (individuo), lo que incrementa la verosimilitud en el ajuste estadístico del valor de cada parámetro y se reducen las desviaciones (^{Pinheiro et al., 2025}; ^{Zuur et al., 2009}).

La corrección de autocorrelación de los errores aplicada en este estudio, al incluir una estructura de tipo ARMA (p=1, q=0), muestra resultados satisfactorios, semejantes a los obtenidos por ^{Hernández-Ramos et al. (2018)} y ^{Flores et al. (2021)} al modelar el di de árboles clonales de Eucalyptus urophylla S. T. Blake en Huimanguillo, Tabasco, México, y de Pinus pseudostrobus Lindl. en el ejido Corona del Rosal, Nuevo León, México, respectivamente. En ambos casos se aplicó un retardo en los residuales (lar1) con una estructura CAR (1). El mismo efecto favorable lo obtuvieron ^{Torres et al. (2020)} para tres procedencias de Pinus caribaea Morelet var. hondurensis (Sénécl.) W. H. Barrett & Golfari y Pinus elliottii Engelm. var. elliottii Engelm. en PF establecidas en las Choapas, Veracruz, México, cuando consideraron una estructura autorregresiva de Primer Orden (AR1).

La robustez del ajuste estadístico tal y como lo citan ^{Zhang et al. (2021)} y demuestran ^{Shin et al. (2022)} al comparar los estadísticos clásicos documentados en estudios de ahusamiento (R ² , RCME, AIC, BIC y logLik), proporcionan la oportunidad de realizar estimaciones confiables de ahusamiento. Además, complementariamente, se puede integrar algebraicamente la ecuación obtenida para dar como resultado su ecuación de Vf implícita o de volumen comercial (Vc) que describe la distribución de productos siguiendo el procedimiento descrito por ^{McTague y Weiskittel (2021)} para expresiones matemáticas de di, como el modelo de Clutter.

Autores como ^{Niño et al. (2018)}, al derivar el modelo de volumen comercial compatible para ajustarlo de manera simultánea, logran resultados satisfactorios, ya que obtienen, justamente, expresiones de Vf y Vc, además de modelar el ahusamiento. Esta estrategia también la aplican ^{Cruz-Cobos et al. (2008)}, ^{Hernández-Ramos et al. (2018)} y ^{Torres et al. (2020)} para modelar el di con distintas técnicas de ajuste y enfoques con lo que generan, en todos los casos, una expresión de Vf aplicable a la muestra de árboles analizada para taxones y regiones específicas.

El modelo de ahusamiento de Clutter es robusto estadísticamente (^{Shin et al., 2022}; ^{Zhang et al., 2021}), por lo que es confiable para describir el perfil fustal de los individuos en las PF de Pinus leiophylla. ^{Ramírez-Martínez et al. (2018)} al modelar el di mediante la expresión de Cielito 2, también determinaron a la expresión de Schumacher-Hall como adecuada para predecir el Vf en árboles de Pinus ayacahuite C. Ehrenb. ex Schltdl. en bosques de Ixtlán de Juárez, Oaxaca, México. Además, al usar expresiones compatibles para estimar el di y el Vf se reducen los sesgos de estimación de forma operativa, como lo describen ^{Pompa et al. (2009)} con la expresión de Biging y la expresión implícita de tipo coeficiente mórfico constante para Pinus arizonica Engelm., en el suroeste de Chihuahua, México.

El modelo para el Vf de Schumacher-Hall derivado de la expresión de ahusamiento de Clutter concuerda con la expresión de ^{Cruz-Cobos et al. (2008)}, pero proveniente del modelo de di de Cielito 1. Los resultados difieren de lo señalado por ^{Torres et al. (2020)}, quienes presentan una forma matemática compleja, poco aplicable en la práctica, la cual contiene en su estructura una expresión de variable combinada, semejante a los modelos de Vf obtenidos en este trabajo para las ecuaciones de Amidon (Ecuación 6) y Newnham (Ecuación 9) (Cuadro 3); también contrasta con lo indicado por ^{Pompa et al. (2009)}, cuyo modelo para Vf es de tipo coeficiente mórfico, el cual es compatible con la expresión de Biging.

En general, el ajuste bajo el enfoque de MEM presentado en este trabajo supera los problemas de heterocedasticidad y multicolinealidad clásicos en datos biológicos y longitudinales (^{Correa & Salazar, 2016}; ^{Pinheiro & Bates, 2000}; ^{Pinheiro et al., 2025}; ^{Zuur et al., 2009}). Además, reduce las desviaciones en la estimación de la variable de interés con mayor certeza mediante una mejora estadística (^{Cruz-Cobos et al., 2008}; ^{Tamarit et al., 2014}). Sin embargo, están aún pendientes por explorar otras estrategias y enfoques de análisis, tal como los modelos basados en la biología de la forma y el crecimiento del tallo, o bien en cuantificar el efecto de la ubicación geográfica de las especies (^{McTague & Weiskittel, 2021}).

Los valores de la matriz de varianza-covarianza registrados dan la pauta para realizar una calibración posterior para expandir su aplicabilidad, como lo proponen ^{Çakir y Kahriman (2018)} y ^{Saygili y Kahriman (2023)} con base en la metodología propuesta por ^{Şenyurt et al. (2017)}, ^{Yang et al. (2009)} y ^{Zhang et al. (2021)}, de modo que pueda usarse para el modelo de ahusamiento de Clutter en el parámetro β3 relacionado a la diferencia entre At y Ai.

Conclusiones

La estrategia de ajuste de MEM aplicada para modelar el ahusamiento de árboles de Pinus leiophylla es satisfactoria para corregir los problemas de heterocedasticidad y de autocorrelación de los errores. También se registra una mejora estadística y una reducción de sesgo individual en las estimaciones, en contraste con el ajuste por Mínimos cuadrados no lineales (MCNL).

El modelo de ahusamiento de Clutter y su respectiva expresión de volumen fustal de tipo Schumacher-Hall son compatibles en sus parámetros y confiables para describir el perfil fustal de árboles de Pinus leiophylla, por lo que pueden utilizarse con certidumbre en la proyección de existencias maderables para un mercado diferenciado de las plantaciones forestales evaluadas de P. leiophylla.