Introducción

La alometría es ampliamente utilizada para modelar la relación entre el cambio proporcional en los atributos de los árboles como: diámetro de copa (dc; m) (^{Cisneros et al., 2019}), altura total (At; m) (^{Guerra-De la Cruz et al., 2019}), volumen fustal (V; m³) (^{Cruz et al., 2016}), biomasa (B; Mg) y carbono (C; kg) (^{Avendaño et al., 2009}), en función del diámetro normal (dn; cm) o del diámetro de tocón (dt; cm) (^{Quiñonez et al., 2012}; ^{García-Cuevas et al., 2017}).

La modelación de las relaciones alométricas entre el dt, dn, dc, At, V, B y C, por lo general, se han ajustado por regresión, usando mínimos cuadrados ordinarios (MCO) (^{García-Cuevas et al., 2017}; ^{Guerra-De la Cruz et al., 2019}) y la generalización de momentos a través de la matriz de covarianza (gmm; por sus siglas en inglés) (^{Hernández et al., 2016}). Sin embargo, estas técnicas de ajuste presentan algunas deficiencias al ser empleadas con datos correlacionados provenientes de la misma unidad de muestreo (sitio o plantación), de registros longitudinales en un mismo individuo o de mediciones repetidas en el tiempo (^{Calama y Montero, 2004}; ^{Tamarit et al., 2014}; ^{Corral et al., 2019}).

El uso de modelos de efectos mixtos (MEM) con la inclusión de covariables se ha convertido en una alternativa de los ajustes tradicionales antes expuestos, porque disminuyen el error y reducen de forma eficiente la variabilidad específica por nivel clasificación (^{Baty et al., 2015}; ^{Correa y Salazar, 2016}; ^{Corral et al., 2019}). Los ajustes con MEM tienen la ventaja de estimar parámetros específicos dentro de los modelos empleados para cada condición evaluada, esto proporciona información útil para identificar patrones estructurados de varianza-covarianza y, mediante la inclusión de covariables, se logra explicar la variabilidad de los efectos aleatorios (^{Pinheiro y Bates, 1998}; ^{García-Espinoza et al., 2018}).

El análisis de las relaciones alométricas en estudios forestales es importante en la gestión de los recursos, ya que estas son una herramienta cuantitativa confiable, que contribuye a disminuir el tiempo y recursos invertidos dentro del levantamiento de la información de campo; estas son fundamentales en la estimación del V utilizado como indicador en la organización de actividades silvícolas y en la evaluación de los aprovechamientos maderables o cortas clandestinas, ya sea en rodales o plantaciones forestales (PF) (^{Prodan et al., 1997}; ^{Castedo et al., 2006}; ^{García-Cuevas et al., 2017}).

Las PF son un elemento para elevar la productividad forestal, almacenar carbono y mitigar los efectos del cambio climático (^{Guariguata, 2009}; ^{Barros e Ipinza, 2011}). En Michoacán, existen alrededor de 16 000 ha de PF establecidas con el género Pinus (^{Comisión Nacional Forestal [Conafor], 2018}); P. pseudostrobus Lindl. y P. montezumae Lamb. son las que representan un potencial de establecimiento importante para extender la superficie y elevar el rendimiento maderable (^{Muñoz et al., 2011}; ^{Muñoz et al., 2015}).

El presente estudio se realizó bajo la hipótesis de que las descripciones de las relaciones alométricas entre dt, dn, dc, altura de fuste limpio (Afl; m), At y V con MEM para PF de Pinus pseudostrobus y Pinus montezumae, establecidas en la Comunidad Indígena de Nuevo San Juan Parangaricutiro (CINSJP), Michoacán, México, se mejoran al incorporar variables aleatorias en los MEM, como la edad y las distintas condiciones de crecimiento (altitud, exposición, pendiente y densidad en las plantaciones) para cada sitio muestreado.

Objetivos

Modelar las relaciones alométricas entre diámetro de tocón, diámetro normal, diámetro de copa, altura de fuste limpio, altura toral y volumen, utilizando modelos de efectos mixtos, para plantaciones forestales de Pinus pseudostrobus y Pinus montezumae establecidas en la Comunidad Indígena de Nuevo San Juan Parangaricutiro (CINSJP), Michoacán, México, con la finalidad de mejorar las estimaciones de las variables dasométricas estudiadas.

Materiales y métodos

Las PF estudiadas se encuentran dentro de la CINSJP, Michoacán, entre los paralelos 19° 17’ y 19° 30’ de latitud norte y los meridianos 102° 06’ y 102° 17’ de longitud oeste, en un intervalo altitudinal de 2000 m a 3100 m. El relieve del área de estudio es de tipo sierra, con una geología ígnea extrusiva y un suelo dominante de tipo Andosol característico de la subprovincia fisiográfica Neovolcánica Tarasca. El clima es templado húmedo con abundantes lluvias en verano (Cw), precipitación media anual de 1500 mm y temperatura promedio de 15 °C (^{Instituto Nacional de Estadística y Geografía [Inegi], 2009}, ²⁰¹⁷).

Se evaluaron 34 PF con edades entre 7 años y 37 años para P. pseudostrobus y 14 plantaciones de entre 7 años y 32 años de P. montezumae; en dichas PF se ubicaron, de manera sistemática, 196 sitios de muestreo cuadrangulares de 400 m² (20 m × 20 m). Una vez que se determinó la edad (E), a partir de contar y promediar los anillos de crecimiento obtenidos de núcleos de madera en tres árboles en cada sitio de muestreo (588 núcleos en total), se midió dt, dn, dc, Afl y At. Cabe mencionar que las dimensiones de la Afl fueron afectadas por las podas heterogéneas aplicadas en las PF o realizadas en etapas tempranas de crecimiento.

Paralelamente, de los árboles se registró su posición vertical dentro del perímetro (Ps): dominante, codominante, intermedio y suprimido, con base en su dominancia vertical (^{Krafl, 1984}; ^{Chauchard et al., 2001}). Se registraron características físicas del lugar como altitud (Alt; m), exposición; Exp; cenital (C), norte (N), noreste (NE), este (Est), sureste (SE), sur (S), suroeste (SW), oeste (W) y noroeste (NW), y pendiente (Pend; expresada en porcentaje). También, con base en los registros por sitio producto del inventario realizado, se calculó la densidad por hectárea (Na) en cada plantación.

La estimación del volumen fustal de cada árbol de P. montezumae se realizó con la ecuación 1 (^{García-Espinoza et al., 2016}); y para P. pseudostrobus, el cálculo se realizó con la ecuación 2 (^{García Espinoza et al., 2018}).

donde:

V = volumen fustal (m³)

dn = diámetro normal (cm)

At =altura total (m)

Los análisis preliminares revelaron la distribución esta dística y la variabilidad de los parámetros dendrométricos empleados, cuyos valores descriptivos se presentan en la tabla 1.

Para el análisis de la relación entre dn y dt se usó un modelo lineal (^{Hernández et al., 2016}); en el caso de At, Afl y dc en función del dt (^{Domínguez et al., 2006}; ^{Quiñonez et al., 2012}; ^{García-Cuevas et al., 2017}) y de At, Afl y dc en función del dn (^{Domínguez et al., 2006}; ^{Scaranello et al., 2012}; ^{Guerra-De la Cruz et al., 2019}), se utilizó el modelo de Chapman-Richards; y para la estimación de V en función del dt y dn, se usó un modelo alométrico (^{García-Espinoza et al., 2016}; ^{García-Espinoza et al., 2018}; ^{García-Cuevas et al., 2017}) (Tabla 2).

De manera preliminar, se realizó un ajuste de los modelos con MCO para las nueve relaciones presentadas en la tabla 2, con el paquete lm para el modelo lineal [1] (^{Santana y Mateos, 2014}) y con nls para los no lineales [2 y 3] (^{Baty et al., 2015}), ambos en programa R® Project ver 3.6 (^{R Core Team, 2016}). Después, los mismos modelos se ajustaron utilizando MEM (^{Calama y Montero, 2004}; ^{Corral et al., 2019}) con el método de máxima verosimilitud (ML), en virtud de que los datos analizados provienen de una misma unidad experimental (plantación) y se presume la no existencia de independencia de observaciones. Esto garantizaría estimadores insesgados de los parámetros (^{Maindonald, 2008}; ^{Pérez y Correa, 2008}) y en las expresiones [1] y [2] de la ^{tabla 2}; para el V, la estimación se realizó con las ecuaciones propuestas por ^{García-Espinoza et al. (2016)} y ^{García-Espinoza et al. (2018)}.

La ecuación 3 (^{Correa y Salazar, 2016}), mostrada a continuación, define la estructura del modelo lineal mixto y las expresiones 4 y 5 detallan sus componentes.

donde:

𝑌 _ij = respuesta para la i-ésima observación de la variable dependiente; i = 1, 2,…, m; j = 1, 2, …, n _i

𝑋 _ij = respuesta para la i-ésima observación de la variable independiente; i = 1, 2,…, m; j = 1, 2, …, n _i

𝑍 _i = matriz de orden o clasificación (n _i × q) de las covariables conocidas

𝛽_𝑖 = vector de q-dimensional de los coeficientes de regression

𝐾_𝑖 = matriz q × p de las covariables conocidas

𝛽 = vector p-dimensional de parámetros de regresión desconocido

D = matriz orden q × q de varianza y covarianza

Ɛ_i= error producto de la sumatoria de la matriz de varianza y covarianzas (∑i) ; se asume N(O,∑i)

En los modelos no lineales, la estructura estuvo representada por las ecuaciones 6 y 7 (^{Baty et al., 2015}; ^{Corral et al., 2019}).

donde:

Y_𝑖 = i-ésima observación dependiente, tomada de la j-ésima unidad de clasificación (covariable)

𝑋_𝑖 = i-ésima observación independiente, tomada de la j-ésima unidad de clasificación (covariable)

𝑓 = función no lineal

𝜀_𝑖 = error de la i-ésima observación dependiente, tomada de la j-ésima unidad de clasificación (covariable)

El vector 𝜃_𝑖 puede ser dividido para los parámetros fijos (𝐴_𝑖) o aleatorios (𝐵_𝑖).

donde:

𝜃 _𝑖j = vector de parámetros r×1; específico del j-ésimo nivel de clasificación por covariable

r = número de parámetros en el modelo

𝐴_𝑖 = matriz de tamaño r×p para los efectos fijos, específicos en cada nivel de clasificación

𝐵_𝑖 = matriz de tamaño r×q para los efectos aleatorios, específicos en cada nivel de clasificación

𝛾 = vector p×1 de los parámetros fijos (covariables), siendo p el número de parámetros fijos

𝑏_𝑖 = vector q×1 de los parámetros aleatorios (covariables), siendo q el número de parámetros aleatorios

El ajuste con MEM se realizó con el paquete nlme y la técnica ML del programa R Project ver 3.6 (^{R Core Team, 2016}), en donde los efectos fijos son utilizados para explicar la variabilidad muestral, mientras que los efectos aleatorios agrupan la variabilidad por nivel de organización (^{Robinson, 2010}; ^{Correa y Salazar, 2016}). La elección de los efectos se realizó bajo el supuesto de tener un valor de agrupación aleatorio muestreado a través del sitio de un intervalo de posibles niveles para cada covariable (Alt, Exp o Pend), de esta manera se margina el valor específico para la condición evaluada (^{Robinson, 2016}).

La inclusión del efecto aleatorio (parámetro a _ij ) dentro de cada parámetro fijo (a ₀ , a ₁ y a ₂ ) (Tabla 2) se realizó de forma individual en cada uno de los modelos [𝑦𝑦=(𝑎₀+𝑎_0𝑗)+𝑎_1𝑥; 𝑦=𝑎₀+(𝑎₁+𝑎_1𝑗)𝑥; 𝑦=(𝑎₀+𝑎_0𝑗)+(𝑎₁+𝑎_1𝑗)𝑥 ] para tratar de englobar y explicar la variabilidad por cada covariable; de esta manera, se estimaron simultáneamente los valores de los parámetros (^{Pinheiro y Bates, 1998}; ^{Calama y Montero, 2004}; ^{Castedo et al., 2006}). Las covariables de E y Na se incluyeron en los MEM como variables continuas al considerar valores entre tres intervalos definidos (Na: 151 - 250 ha^-1, 251 - 350 ha^-1 hasta 1251 - 1350 ha^-1), mientras que, las covariables de Ps, Alt, Exp y Pend fueron empleadas como variables categóricas de clasificación (Exp: N, NE, NW).

El ajuste de modelos lineales o no lineales con MEM asume que los residuales siguen una distribución normal con media igual a cero y una matriz de varianzascovarianzas capaz de englobar la variabilidad de la muestra utilizada por nivel de clasificación (Ecuación 8) (^{Littell et al., 2006}; ^{Corral et al., 2019}).

donde:

𝑅𝜇 = matriz diagonal cuya 𝑖𝑖 es un elemento diagonal 𝑉ar(𝜇) en función de varianza para la observación 𝑖

Rμ0.5 = matriz diagonal de raíces cuadradas de los elementos correspondientes de 𝑅𝜇

𝐴 = matriz de parámetros de escala

La normalidad y homocedasticidad se analizaron con la prueba de Shapiro-Wilk (SW), la cual indica que valores cercanos a la unidad señalan normalidad en la distribución, y gráficamente la distribución de los residuales no debe tener una tendencia más amplia a medida que incrementan las dimensiones de las variables (^{Martínez-González et al., 2006}; ^{SAS Institute Inc., 2014}).

La bondad de ajuste de los dos enfoques (MCO y MEM) se evaluó a través de la comparación numérica del valor del coeficiente de determinación (R ² ), raíz del cuadrado medio del error (RCME) y criterio de información de Akaike (AIC) (^{Kozak y Smith, 1993}; ^{Prodan et al., 1997}; ^{Lencinas y Mohr-Bell, 2007}). Para seleccionar el mejor modelo, se generó un sistema de calificación (SC), en el cual se jerarquizaron los valores de los estadísticos de ajuste: valores de 1 a n en función de importancia (1 al mejor valor del estadístico y n al peor valor), para después realizar una sumatoria entre las calificaciones y obtener como mejor modelo aquel con el menor valor global (^{Sakici et al., 2008}; ^{Tamarit et al., 2014}; ^{Ordóñez-Prado et al., 2023}).

Posteriormente, se evaluó la ganancia estadística porcentual de ambos enfoques de ajuste y se determinó la precisión de las estimaciones de los modelos seleccionados como los mejores a través del error promedio (Ep) y error promedio relativo (𝐸𝑅). Además, se graficaron las tendencias alométricas particulares entre las relaciones por covariable al emplear la condición de agrupación específica para cada caso.

Resultados

Las PF de P. pseudostrobus y P. montezumae se ubicaron a una altitud en el límite inferior de 2000 m y 2200 m s.n.m., mientras que, el límite superior fue a 3000 m s.n.m. Además, las plantaciones se encuentran en todas las exposiciones topográficas y pendientes que oscilan entre 0% y 25%. La densidad de árboles fue entre 150 ha^-1 y 1400 ha^-1 en las PF de P. pseudostrobus, mientras que, en las de P. montezumae lo hizo entre 140 ha^-1 y 800 ha^-1.

Al ajustar los modelos de la tabla 2 para todas las relaciones alométricas, por MCO (18 relaciones) e MEM (156 relaciones) incluir las covariables (E, Ps, Alt, Exp, Pend y Na) en cada uno de los parámetros y la combinación entre ellos como se muestra en las ecuaciones 9 a 14, se presentaron, en algunos casos, resultados de ajuste con parámetros no significativos (p > 0.05), incluso, con errores estándar (Ee) y valores de t muy grandes con respecto a los demás ajustes en la misma variable y una covariable específica incluida en el modelo, o bien, una falta de convergencia de los parámetros, por lo que se eliminaron de la tabla de calificación final. Al comparar los resultados entre MCO y MEM, se observó que en el segundo enfoque incrementaron los valores del R² y se redujeron la RCME, el AIC, y el Ee en los parámetros; igualmente, como mejora la explicación de la variabilidad (Tabla 3).

Para las relaciones entre dn-dt, At-dt, Afl-dt, dc-dt, At-dn, Afldn y dc-dn, en Pinus pseudostrobus, se notó una mejora estadística al utilizar MEM, que en promedio (valores de R², RCME y AIC, para cada relación; Tabla 4), fue de 1.0%, 2.4%, 7.8%, 5.0%, 4.2%, 9.5% y 4.5%, respectivamente; mientras que para, P. montezumae, los resultados mostraron un incremento de 0.7%, 4.0%, 19.6%, 19.7%, 3.8%, 20.6% y 17.5%, en el mismo orden. El error promedio, en todos los casos, estuvo por debajo de las 1.5 unidades, en tanto que, el error relativo en particular para estimar la Afl fue alto, probablemente debido a la heterogeneidad en la forma y tiempo de aplicación de podas en los individuos de ambas especies (Tabla 4).

La inclusión del efecto aleatorio en la relación alométrica dn-dt para P. pseudostrobus indicó que la edad influye directamente en la ordenada al origen (a ₀ ) o constante del modelo lineal (Tabla 3). Para los casos At-dt, Afl-dt, dc-dt, At-dn, Afl-dn y dc-dn, el parámetro referido al máximo crecimiento asintótico de cada variable en el modelo de Chapman-Richards (Tabla 2) fue el de mayor influencia ( a0 ), en el cual la Alt, Exp o Na son los factores que determinan en mayor medida las dimensiones y dinámica de las variables (Tabla 3). Los parámetros globales de cada relación alométrica fueron significativos (p < 0.05), tal como se observa en la tabla 5.

Para las relaciones alométricas de P. montezumae, las covariables muestran una mayor influencia sobre el parámetro relacionado con la pendiente ( a1 ) dentro del modelo lineal. Mientras que, en Chapman-Richards (Tabla 2), de acuerdo con la covariable afecta tanto a la asíntota del modelo ( a0 ) como a las distintas velocidades de cambio en las dimensiones del arbolado ( a1 ). En este caso, la densidad (Na) es la característica que mayor influencia ejerce en las relaciones alométricas, en consecuencia, sobre crecimiento y rendimiento de la especie (Tabla 3). Al ajustar las relaciones alométricas At-dt, Afl-dt, At-dn y Afl-dn, el parámetro a ₂ toma el valor de 1 y el modelo de ChapmanRichards cambia a un modelo Monomolecular (Tabla 5), cuya estructura es:

En la estimación del V en función del dt y dn para las dos especies, se observó que no existe mejora en las estimaciones al incluir dentro de la expresión [3] los efectos con MEM, ya que el valor de la R ² disminuye, los valores de RCME y AIC aumentan, y se incrementan los errores expresados por el error promedio y ER. Por esta situación y aunado a que los volúmenes se estimaron a través de ecuaciones ya existentes, se decidió proponer únicamente los parámetros del modelo ajustado por MCO (Tabla 5) y no realizar correcciones por distribución heterocedástica de los residuales, como se observa en el material complementario: tabla S1 (j; r) para P. pseudostrobus y tabla S2 (j; r) para P. montezumae.

La prueba de Shapiro-Wilks para cada relación alométrica indica normalidad, ya que en todos los casos los valores fueron mayores a 0.91; además, la distribución en la frecuencia de los residuales tiene una forma de campana de Gauss (Material complementario). Para ambas especies en las relaciones dn-dt, At-dt, At-dn y dc-dn, se observó homocedasticidad con respecto a los datos predichos, mientras que para las demás relaciones alométricas se tiene un aumento del residual a medida que incrementa la dimensión de la variable dependiente, sin que esta repercuta en las estimaciones al aplicar la expresión propuesta (Material complementario).

Al utilizar, por especie, el valor de los parámetros globales en las expresiones obtenidas (Tabla 5) y contrastar las estimaciones de dn, At, Afl, dc y V (Fig. 1), se observó que la relación dn-dt para las dos especies gráficamente fue similar (Fig. 1a). En las relaciones de At, dc y V en función del dt y dn, respectivamente, P. pseudostrobus fue superior cuantitativamente en las proporciones alométricas a P. montezumae. Esta última especie refleja un menor desarrollo y crecimiento en PF. De manera consistente, los modelos de dt fueron más conservadores que los que emplean el dn (Fig. 1 b, d, c y e).

Figura 1 Estimación de variables alométricas para las especies de Pinus pseudostrobus y Pinus montezumae establecidas en plantaciones forestales en Michoacán, México. 

Al graficar las tendencias de las relaciones alométricas ajustadas por MEM, dn-dt no mostró diferencias evidentes por nivel de agrupación para ninguna de las especies, pero sí una mejora estadística, por lo cual no se incluyen en las figuras 2 y 3. Para las relaciones Afl-dt, dc-dt, At-dn y dc-dn, P. pseudostrobus presenta un mejor desarrollo entre las altitudes de 2400 m y 2600 m pero menor a los 3000 m (Fig. 2 b, c, d y g). En tanto que la At en función del dt muestra mayores dimensiones en las exposiciones orientadas al sur (Fig. 4a) y la Afl en función del dn tiene menores dimensiones a partir de la densidad de 600 árboles ha^-1 en adelante (Fig. 4f). Este último resultado deberá de ser tomado con precaución, ya que la aplicación de podas y manejo no fueron consistentes a lo largo del desarrollo de las PF.

Figura 2 Tendencias entre las relaciones alométricas por covariable de Pinus pseudostrobus en plantaciones forestales en Michoacán, México. 

Figura 3 Tendencias entre las relaciones alométricas por covariable de Pinus montezumae en plantaciones forestales. 

Para las relaciones alométricas en P. montezumae, la variable que mayores dimensiones propició para la Afl-dt, dc-dt y Afldn (Fig. 3b, c y e), fue la densidad de 700 ha^-1 árboles. La densidad mencionada junto con la densidad de 800 ha^-1 árboles, son las que menores alturas registraron en la relación At-dt (Fig. 3a); dicho resultado posiblemente se debe a la apertura de espacios de crecimiento debido al manejo silvícola de las PF. No obstante, al analizar los otros componentes de agrupación de las plantaciones ocurrió lo contrario. Por ejemplo, al evaluar las tendencias de la At-dn, la exposición topográfica sureste es la que mayores dimensiones registró (Fig. 3d), en tanto que, la relación alométrica entre el dc-dn mostró un desarrollo foliar mayor en el intervalo de altitud de 2700 m a 2900 m (Fig. 3f).

Discusión

El enfoque MEM permitió identificar relaciones alométricas significativas que el método MCO no detectó, mejorando así la comprensión de la dinámica de crecimiento de las especies de pino estudiadas. Estos resultados destacan la relevancia de incluir factores aleatorios en el modelado de relaciones alométricas, ya que permiten representar de manera más precisa la complejidad y heterogeneidad de los ecosistemas forestales; además, al agrupar la variabilidad específica por nivel de clasificación, mejora la explicación de las relaciones alométricas descritas en esta investigación, tal y como lo han discutido ^{Calama y Montero (2004)}, ^{Baty et al. (2015)}, ^{Correa y Salazar (2016)} y ^{Corral et al. (2019)} en distintos enfoques de aplicación de los MEM. También, los ajustes estadísticos mostraron parámetros específicos para cadanivel de agrupación, los cuales dan origen a tendencias particulares de crecimiento en cada especie (Fig. 2 y 3; ^{Pinheiro y Bates, 1998}; ^{Castedo et al., 2006}; ^{Tamarit et al., 2014}; ^{Corral et al., 2019}), por lo que puede aplicarse un manejo diferenciado. Por ejemplo, el análisis con MEM permitió deducir que el índice de esbeltez (IE: At dn ^-1 ), al considerar el dt de P. pseudostrobus, es influenciado directamente por la exposición del sitio: la orientación sur es la que propicia las mejores condiciones de desarrollo para la especie. Este análisis también evidenció que el dc está limitado por la altitud. Para P. montezumae, la covariable con mayor influencia en At, Afl y dc en función del dt y Afl-dn, es la densidad; mientras que la velocidad de cambio de las proporciones en las relaciones de ambas especies entre el dn-dt disminuye con respecto a la edad o la altitud.

Estas variaciones alométricas y las distintas tendencias de desarrollo por covariable se pueden atribuir a las múltiples condiciones que generan en conjunto la altitud, el clima y la topografía, tal y como lo documentan ^{Rehfeldt (1988)}, ^{Sundqvist et al. (2013)} y ^{Pyatt et al. (2016)}, al mencionar que las características de cada sitio y la dinámica de interacción entre los individuos del rodal influyen significativamente en las tasas de crecimiento, incremento y rendimientos de las comunidades vegetales. Además, las tendencias en las relaciones alometrícas entre las variables de interés pueden variar como respuesta a las alteraciones climáticas y de altitud (^{Sundqvist et al., 2013}; ^{Sáenz-Romero et al., 2016}).

Como un ejemplo de lo expuesto en el párrafo anterior, ^{Körner y Paulsen (2004)}, ^{Smith et al. (2009}) y ^{Ayma-Romay et al. (2015)} señalan que la variación de temperatura causada por la altitud y condiciones topográficas (isotermas) delimitan la distribución y el desarrollo de las especies, y que condiciones microclimáticas propiciadas por la exposición del terreno a vientos o radiación solar tienen una influencia importante en sobrevivencia, desarrollo y crecimiento de plántulas. Estos resultados se ven reflejados en esta investigación, especialmente en P. pseudostrobus, donde las principales covariables incorporadas al MEM, que mejoraron los modelos alométricos (At-dt, Afl-dt, dc-dt, At-dn, y dc-dn), fueron la altitud y la exposición; sin embargo, para P. montezumae solo las relaciones dn -dt y At-dn mejoraron incorporando las covariables altitud y exposición, respectivamente.

En cambio, para P. montezumae, la covariable involucrada en los MEM (At-dt, Afl-dt, dc-dt, Afl-dn y dc-dn) que incrementa la ganancia estadística fue la densidad, lo cual concuerda con lo encontrado por ^{Anthelme y Dangles (2012)}, quienes señalan que las características bióticas y abióticas son determinantes en las interacciones plantaplanta, y las limitantes ambientales generan respuestas diversas de crecimiento y de nicho en los organismos. En cambio, para P. pseudostrobus, solo la relación At-dn, incorporando covariable densidad tuvo un incremento al usar MEM.

El modelo de Chapman-Richards (Tabla 2), ajustado con MEM, fue el más eficiente para la mayoría de las relaciones alométricas. En este sentido y con fundamento en los estudios sobre los modelos de crecimiento publicados por ^{Kiviste et al. (2002)}, ^{Chiappa et al. (2009)} y ^{Gutiérrez y Sánchez (2017)}, se puede observar, en P. montezumae, que entre dc-dt el parámetro de efecto aleatorio fue el referente a la velocidad de cambio ( a1 ), mientras que en las relaciones restantes para ambas especies (A, Alf y dc - dn y dt) fue el parámetro asintótico, que es un indicativo de la mayor dimensión que se puede alcanzar en cada relación alométrica ( a0 ).

Lo anterior es evidencia de que la elección del lugar de plantación es crucial para el crecimiento y rendimiento de las especies, ya que las condiciones específicas de microclima que propician altitud, exposición y pendiente en conjunto determinan la dinámica en las relaciones alométricas en las PF de P. pseudostrobus y P. montezumae. Esta información puede ser pauta para el establecimiento de nuevas plantaciones, en las condiciones que incrementen las posibilidades de éxito y mejoren su rendimiento maderable, mediante la consideración de los gradientes altitudinales (^{Smith et al., 2009}), las características particulares de cada lugar (^{Anthelme y Dangles, 2012}), o las condiciones climáticas cambiantes (^{Ayma-Romay et al., 2015}; ^{Sáenz-Romero et al., 2016}), factores que en conjunto hacen posible el hábitat óptimo de cada especie.

Los hallazgos de este estudio tienen implicaciones clave para el manejo silvícola, ya que demuestran que factores como la exposición y la densidad de árboles influyen significativamente en las relaciones alométricas de P. pseudostrobus y P. montezumae. Estudios previos han señalado que la apertura del dosel y la competencia intraespecífica afectan directamente el crecimiento y la forma de los árboles en plantaciones forestales (^{Pretzsch, 2018}). La manipulación del espaciamiento y la densidad de plantación puede mejorar la eficiencia en la producción de madera, optimizando características como el índice de esbeltez y la altura del fuste limpio (^{Ducey y Knapp, 2010}).

En la práctica, se recomienda implementar estrategias como podas y aclareos selectivos para mejorar el rendimiento y la calidad de los individuos, maximizando así la producción forestal sostenible (^{Weiskittel et al., 2011}). Además, los resultados obtenidos respaldan la hipótesis planteada, ya que la inclusión de variables aleatorias en MEM mejoró la estimación de las relaciones alométricas en Pinus pseudostrobus y P. montezumae.

Se identificó que la exposición influye en la relación de las variables alométricas (IE) de P. pseudostrobus, mientras que, la densidad arbórea tiene una interacción importante en la respuesta de las variables dasométricas en P. montezumae, lo que concuerda con estudios previos que evidencian la influencia del manejo silvícola y la competencia intraespecífica en la arquitectura y crecimiento de los árboles (Pretzsch et al., 2015; ^{Weiskittel et al., 2011}). Estos resultados refuerzan la importancia de considerar factores estructurales y ambientales para mejorar la modelación de parámetros dasométricos en PF.

El estudio presenta fortalezas metodológicas notables, como el uso de MEM y pruebas estadísticas robustas, lo que permitió captar variabilidad no explicada por modelos tradicionales. No obstante, se identificaron limitaciones, como la falta de convergencia de algunos parámetros y la variabilidad en la altura del fuste limpio, posiblemente atribuida a diferencias en la poda. Para optimizar la estimación del volumen fustal, se recomienda evaluar modelos heterocedásticos y considerar interacciones genéticas y estructurales en la modelación (^{Zhao et al., 2019}; ^{Ducey y Knapp, 2018}). Asimismo, se sugiere ajustar estrategias de manejo silvícola, como la regulación del espaciamiento y la poda, para mejorar la homogeneidad del crecimiento y la productividad de las plantaciones.

Conclusiones

El estudio confirma que la inclusión de covariables en el ajuste de los modelos de efectos mixtos (MEM) conlleva a una mayor robustez estadística con respecto a los modelos ajustados por el método de mínimos cuadrados ordinarios. A través de los MEM, se confirma, en gran medida, la hipótesis planteada, ya que, al incorporar covariables como la altitud, la exposición topográfica y la densidad de plantación, incrementó el coeficiente de determinación (R²): en el caso de P. pseudostrobus, la mejora en la calidad del ajuste osciló entre 1.0% y 9.5%, mientras que en P. montezumae los incrementos fueron aún mayores, alcanzando hasta 20.6% en la relación alométrica de Afl-dn. Asimismo, se redujeron tanto la raíz del cuadrado medio del error (RCME) como el criterio de información de Akaike (AIC).

Además, se demostró la influencia cuantitativa de las covariables en la ordenada al origen del modelo lineal, así como en la máxima dimensión (asíntota: a ₀ ) y velocidad de crecimiento (a ₁ ) del modelo de Chapman-Richards, la cual es propiciada por la altitud y topografía (variables aleatorias) en donde se desarrollan las plantaciones forestales (PF).

La diferencia de edad y densidad entre las plantaciones, así como de las condiciones ambientales de cada lugar de establecimiento propiciadas por la altitud, exposición y pendiente, tienen una influencia directa en las relaciones alométricas entre las variables de diámetro de tocón, diámetro normal, diámetro de copa, altura de fuste limpio, altura total y, por consecuencia, en volumen (rendimiento maderable) en las PF de estas dos especies.

Entre las fortalezas del estudio, destaca la utilización de un enfoque metodológico robusto que permite una mayor precisión en la predicción de variables dasométricas mediante la integración de factores ambientales y estructurales en los modelos. Si bien los hallazgos sugieren que las modificaciones en el enfoque utilizado contribuyen a una mayor precisión en las estimaciones, aún existen factores que requieren mayor control y análisis como la genética, la sanidad y la competencia intraespecífica, entre otras.