Área de estudio

El estudio se llevó a cabo en los rodales bajo aprovechamiento forestal localizados en la Umafor 2001, ubicada al noreste del estado de Oaxaca, México, y que abarca una superficie de 195 397 ha. Se localiza a los 17°22’58” N y los 96°22’29” O (Figura 1). La altitud de las áreas de muestreo varía de 1 800 a 3 361 m. La temperatura anual oscila entre -3 °C y 28 °C, con precipitación promedio anual de 3 797 mm (^{Granados-Sánchez, 2009}). El tipo de vegetación predominante corresponde a bosques de pino-encino, con dominancia de Pinus oaxacana, P. douglasiana, P. patula, P. pseudostrobus, P. leiophylla Schiede ex Schltdl. & Cham., P. pringlei Shaw., P. rudis Endl., P. teocote Schiede ex Schltdl. & Cham., Quercus crassifolia Bonpl., Q. laurina Bonpl. y Q. rugosa Née, principalmente.

Figura 1 Localización geográfica del área de estudio. 

Datos

Para ajustar las ecuaciones dinámicas de altura dominante e IS se usó una base de datos procedente de 206 análisis troncales, integrada por 44 árboles de Pinus oaxacana, 45 de P. douglasiana, 58 de P. patula, y 59 de P. pseudostrobus. Los individuos fueron derribados a una altura del tocón que varió de 0.1 m a 0.25 m, de tal manera que la primera sección correspondió a la altura del tocón, la segunda a 1.3 m y posteriormente, cada 2.54 m hasta llegar a la punta del árbol. Los ejemplares se seleccionaron en rodales mezclados y se intentó cubrir las diferentes calidades de estación, así como las categorías de diámetro presentes en las áreas con manejo forestal. Si bien, el concepto de índice de sitio está diseñado para masas coetáneas y puras (^{Daniel et al., 1979}), el carácter de dominante para un árbol implica que estos individuos son poco sensibles al efecto de la densidad, la competencia y la mezcla de especies, ya que su atura se mantiene constante en un intervalo amplio de densidades y grados de mezcla, lo cual refleja a lo largo de su vida el potencial productivo de un sitio en particular, y por ello el concepto se ha hecho extensivo a las masas de pino mezcladas (^{Corral et al., 2004}; ^{Vargas-Larreta, 2010}; ^{Castillo et al., 2013}).

Los estadísticos descriptivos de la base de datos se presentan en el Cuadro 1, se observa que la muestra cubre un intervalo amplio de combinaciones edad-altura dominante, por lo que se considera adecuada para el ajuste de las ecuaciones dinámicas de crecimiento.

Las alturas verdaderas de las secciones se estimaron con el método de ^{Carmean (1972)} modificado por ^{Newberry (1991)}, ya que con este se han tenido buenos resultados en otras investigaciones (^{Dyer y Bailey, 1987}; ^{Fabbio et al., 1994}; ^{Castillo et al., 2013}). El método se basa en dos supuestos: i) entre dos secciones, el árbol crece a un ritmo constante, y ii) el corte se realiza, como promedio, en el centro del crecimiento en altura de un año (^{Dyer y Bailey, 1987}; ^{Amaral et al., 2010}). Las ecuaciones empleadas para calcular la altura verdadera varían según la sección del árbol (ecuaciones 1-3).

Donde:

H ₁ y H ₂ = Alturas de las secciones inferior y superior de la troza

N ₁ y N ₂ = Número de anillos de las secciones inferior y superior de la troza

N ₀ = Edad del árbol; es decir, el número de anillos del tocón (en el tocón N ₀ = N ₁ )

T ₀ = Edad del árbol cuando alcanzó la altura H₁, es decir N ₀ - N ₁

T = Número entero de 1 a N ₁ - N ₂

Ecuaciones dinámicas de índice de sitio

Se ajustaron tres ecuaciones con el método GADA, las cuales se han documentado en diversos estudios para describir el crecimiento en altura dominante e índice de sitio, esta son: Korf (4), Hossfeld (5) y Bertalanffy-Richards (6) (^{Cieszewski, 2002}; ^{Castillo et al., 2013}; ^{Quiñonez-Barraza et al., 2015}; ^{González et al., 2016}). La expresión de los modelos base, así como su formulación GADA se presentan en el Cuadro 2.

Descripción del método GADA

El desarrollo de cualquier ecuación dinámica con la formulación GADA considera: i) seleccionar una ecuación base e identificar el número de parámetros que serán dependientes de la productividad del sitio; ii) los parámetros seleccionados se expresan como funciones de la calidad de estación definida por la variable X ₀ (variable no observable e independiente que describe la productividad del sitio, como resultado de la suma de factores ecológicos, climáticos, regímenes de manejo, condiciones del suelo y los parámetros nuevos; iii) la ecuación base bidimensional (H=f(t)) seleccionada es expandida a una ecuación tridimensional (H=f(t, X ₀ )) de índice de sitio; y iv) se despeja el valor de X ₀ , a partir de condiciones iniciales del sitio; es decir, de valores de partida de altura dominante - edad (t ₀ , H ₀ ), de forma que el modelo pueda definirse implícitamente y aplicarse en la práctica (H=f(t, t ₀ , H ₀ )) (^{Cieszewski y Bailey, 2000}; ^{Cieszewski, 2001}; ^{Cieszewski, 2002}; ^{Cieszewski, 2003}). Durante el proceso, a menudo se eliminan parámetros redundantes y se obtiene un modelo con un número de parámetros menor o igual que la ecuación base original.

Método y estadísticos de ajuste para la comparación de los modelos

Los parámetros globales y específicos del sitio se estimaron con el método iterativo anidado (nested iterative procedure) (^{Tait et al., 1988}; ^{Cieszewski, 2003}; ^{Krumland y Eng, 2005}), el cual se emplea cuando las bases de datos proceden de análisis troncales o parcelas permanentes. Los modelos se ajustaron con estructura de error autorregresiva de tiempo continuo (CAR_i), para corregir problemas potenciales de autocorrelación del término del error (^{Zimmerman et al., 2001}; ^{Nord-Larsen, 2006}; ^{Crecente-Campo et al., 2009}) y obtener estimadores de los parámetros insesgados y eficientes (^{Parresol y Vissage, 1998}).

Donde:

Hij = Predicción de la altura i utilizando Hj (altura j), ti (edad i), y tj (edad j≠i) como variables predictores

β = Vector de parámetros a estimar

eij = Error

d1=1 para j> 1 y cero cuando j=1

ρ1 = Parámetro autorregresivo de primer orden cuyo valor debe estimarse

tij-tij-1 = Diferencia temporal que separa la observación jth de la observación jth-1 en cada análisis troncal, t_ij > t_ij-1

La determinación de la bondad de ajuste de los modelos se basó en resultados de análisis cuantitativos, como la raíz del error medio cuadrático (REMC), el coeficiente de determinación ajustado (R ² _adj ), el sesgo promedio (e) y gráficos de residuales para comparar los aspectos siguientes: i) comportamiento biológico de las predicciones de los diferentes modelos ajustados; ii) representación de los residuales frente a los valores predichos por el modelo, y iii) representación de los residuales frente a residuales con diferentes retardos para comprobar la corrección de la autocorrelación de los errores, mediante la modelación de la estructura del error (^{Goelz y Burk, 1992}; ^{Sharma et al., 2011}). La presencia de la autocorrelación en los datos utilizados se calculó con el estadístico de Durbin-Watson (^{Durbin y Watson, 1951}). Las expresiones de dichos estadísticos son:

Donde:

H,H^ e H- = Valores observados, predichos y medio en pares de datos de altura dominante - edad

n = Número de observaciones

p = Número de parámetros del modelo

e = Valor residual del modelo ajustado

El ajuste simultáneo de las ecuaciones de crecimiento y de la estructura del error (CAR2) se realizó con el procedimiento MODEL del paquete estadístico SAS/ETS^TM (^{SAS, 2004}), que permite una actualización dinámica de los residuos en el procedimiento de ajuste.

Ecuaciones GADA generadas, sin considerar la estructura del error

Los parámetros estimados, los errores estándar y los estadísticos de ajuste de las ecuaciones dinámicas analizadas para Pinus oaxacana, P. douglasiana, P. patula y P. pseudostrobus, se obtuvieron para todos los pares de altura dominante - edad, con el ajuste de las ecuaciones 4, 5 y 6, sin considerar la modelación de la estructura del error (Cuadro 3). En los resultados se observan problemas de autocorrelación de los residuos, debido a las mediciones de altura en un mismo árbol (Figura 4).

Figura 2 Sobreposición de las curvas de IS generadas con los modelos de Korf (línea discontinua), Hosffeld (línea punteada) y Bertalanffy-Richards (línea continua). 

Las curvas se generaron con el modelo dinámico de Bertalanffy-Richards sin (línea punteada) y con (línea continua) corrección de la autocorrelación mediante una estructura autorregresiva CAR2.

Figura 3 Curvas de índice de sitio a la edad de referencia de 40 años para Pinus oaxacana Mirov (Po), P. douglasiana Martínez (Pd), P. patula Schiede ex Schltdl. & Cham. (Pp) y P. pseudostrobus Lindl. (Pps).  

Figura 4 Gráficas de la Función de Autocorrelación Parcial de la ecuación de Bertalanffy-Richards (Ecuación 6): a) ajuste sin corrección de la autocorrelación, b) al considerar un modelo CAR1, y c) con un modelo CAR2. 

Los coeficientes de determinación ajustados (R ² _adj ) por el número de parámetros en los modelos de Korf, Hossfeld y Bertalanffy-Richards para las diferentes especies muestran que es posible explicar más de 96 % de la varianza observada en el crecimiento en altura dominante, en función de la edad, con estimaciones de parámetros significativos para la mayoría de las combinaciones de ecuaciones y especies (p<0.01). Los errores medios en los tres modelos oscilan entre 1.18 m (Pd) y 1.85 m (Pp) (Cuadro 3).

En general, las tres ecuaciones ajustadas muestran estadísticos de ajuste satisfactorios para predecir el crecimiento en altura dominante en función de la edad, sin expandir los términos del error para corregir la autocorrelación. Por otro lado, el análisis gráfico muestra un comportamiento biológico realista (predicciones similares a los datos observados). Sin embargo, en Pinus oaxacana el parámetro b ₁ fue no significativo en los modelos Korf (4) y Hossfeld (5) por lo que se acepta la hipótesis nula de que los estimadores de los parámetros son iguales a cero; escenario que también ocurrió para el modelo Bertalanffy-Richards (6) en el parámetro b ₂ . El parámetro b ₂ no fue significativo en el modelo de Hossfeld (5) para Pinus patula Schiede ex Schltdl. & Cham., situación que se presentó para el mismo parámetro en el modelo Bertalanffy-Richards (6) para Pinus pseudostrobus, todos a un nivel de significancia de 5 %.

Los resultados indican que las diferencias entre los estadísticos de bondad de ajuste fueron mínimas entre los modelos evaluados, por lo que la selección del más apropiado se apoyó en el análisis gráfico de sobreponer las familias de curvas de crecimiento de altura dominante a las trayectorias de crecimiento de los datos observados; se consideraron categorías de IS de 15 m a 35 m para Pinus oaxacana, de 10 m a 30 m para Pinus douglasiana, con intervalos de 5 m; y de 14 m a 38 m, para Pinus patula y Pinus pseudostrobus, con intervalos de 6 m a una edad de referencia de 40 años para todos los casos (Figura 2). Lo anterior, con base en que diferentes modelos pueden presentar los mismos estadísticos de bondad de ajuste o de comparación, pero una respuesta distinta a las trayectorias de los datos observados; así, unos subestiman en las primeras edades y sobreestiman en las avanzadas, o viceversa (^{Diéguez-Aranda et al., 2006}).

Se seleccionó el modelo de Bertalanffy-Richards (6) como el más adecuado para describir las trayectorias de crecimiento de los árboles que integraron la base de datos, y para comparar las estimaciones de crecimiento frente al mismo modelo, a partir de la modelación de la estructura del error, con una estructura autocorregresiva de errores continua tipo CAR1 y CAR2, en virtud de que este modelo describió ligeramente mejor las tendencias individuales de crecimiento en altura dominante, principalmente, a edades superiores a 50 años (Figura 2).

La ecuación de Bertalanffy-Richards, en la que se considera una estructura CAR2 fue capaz de modelar las tendencias de crecimiento de los datos experimentales para las especies estudiadas, y describir patrones de variación en el crecimiento de las curvas construidas, así como potencialidades de sitio diferentes. Este resultado es consistente con otros estudios en los que también se recomienda dicho modelo para la construcción de curvas de índice de sitio de otras especies de interés comercial (^{Kiviste et al., 2002}; ^{Rodríguez-Carrillo et al., 2015}; ^{Castillo et al., 2013}).

Modelación de la estructura del error

En el Cuadro 4 se muestran los estadísticos del ajuste del modelo de Bertalanffy-Richards, en el que se integra un modelo autorregresivo continuo de segundo orden (CAR2), el cual corrigió de manera adecuada la autocorrelación de los errores de la base de datos (valor de DW cercano a 2).

Los errores promedio de la REMC del modelo 6, con la estructura del error oscilaron de 0.79 m (Pd) a 1.18 m (Pp) (Cuadro 4), y explicaron más de 98 % de la varianza del crecimiento en altura dominante. El modelo autorregresivo continuo de segundo orden (CAR2) garantizó la independencia de los residuales en el proceso de ajuste.

En la Figura 3 se observan las curvas de IS obtenidas con la ecuación de Bertalanffy-Richards sin considerar la estructura del error (línea punteada), así como aquéllas en las que se tomó en cuenta la estructura del error CAR2 (línea continua) a una edad de referencia de 40 años, sobrepuestas a los datos experimentales que se utilizaron en su ajuste. Las curvas de IS cubren las diferentes edades y describen la relación entre crecimiento en altura dominante y edad con un comportamiento lógico y biológico (polimorfismo complejo) a edades mayores a los 50 años (^{Cieszewski, 2003}; ^{Álvarez-González et al., 2010}). Se aprecia, gráficamente, que las especies estudiadas presentan distintos patrones de crecimiento en altura dominante, lo cual indica que, el manejo de rodales con dominancia de las cuatro especies de pino requiere de esquemas silvícolas distintos (^{Attis et al., 2015}).

En los gráficos de la función de autocorrelación parcial (PACF), obtenida del ajuste de la ecuación Bertalanffy-Richards con 10 retardos (Lags) de los residuales para cada árbol (Figura 4), se observa que se pasa de un correlograma con autocorrelación positiva de segundo orden, a un correlograma sin evidencia de autocorrelación, estimada con el ajuste de la ecuación dinámica con una estructura autocorregresiva de errores continua tipo CAR2; resultados que son similares a los registrados por ^{Vargas-Larreta et al. (2010)}, ^{Vargas-Larreta. (2013)} y ^{Quiñonez-Barraza et al. (2015)}, para especies forestales en diferentes regiones de México.

Al corregir la autocorrelación, se determina el patrón aleatorio de los residuales alrededor de la línea de cero, con varianza homogénea y sin detectarse una tendencia definida (Figura 5). También, se observa la evolución del sesgo y RCME, y se nota que a partir de la categoría de edad mayor a 50 años, se presenta un aumento considerable en la dispersión y alturas predichas debido, principalmente, a la falta de pares de datos de altura-edad de esas clases de edad.

Figura 5 Residuales vs alturas predichas y evolución del sesgo y raíz del cuadrado medio del error (RCME) por categoría de edad del modelo de Bertalanffy-Richards, con el modelo autorregresivo CAR(2).