Caracterización del área de estudio

La información dasométrica que se utilizó fue colectada en predios forestales con programa de manejo para el aprovechamiento maderable de la Unidad de Manejo Forestal (Umafor) 2108, que corresponde a la región de Chignahuapan-Zacatlán, localizada al norte del estado de Puebla; geográficamente se localiza entre los paralelos 20°07’06” y 19°44’18” de latitud norte y entre los meridianos 97°57’18” y 97°38’42” de longitud oeste; se ubica en la Región Hidrológica (RH27) de la cuenca del Río Tecolutla. Esta Umafor está constituida por la Asociación Regional de Silvicultores de Chignahuapan-Zacatlán, A.C., conformada por silvicultores dueños poseedores de terrenos forestales.

Con la actual política pública del gobierno federal mexicano de incrementar la producción y productividad maderable, se le ha identificado como una importante cuenca de abasto forestal, con la incorporación de mejores prácticas y técnicas silvícolas intensivas que contribuyan a fortalecer y a mejorar el manejo forestal.

Tamaño de muestra y base de datos

Se procesó una muestra compuesta por 252 pares de observaciones sobre el número de árboles y el diámetro cuadrático promedio; fue seleccionada con base en un amplio conjunto de sitios de muestreo para inventario maderable, situados en predios forestales de los municipios Aquixtla, Chignahuapan, Ixtacamaxtitlán, Tetela de Ocampo y Zacatlán; los sitios fueron circulares, con una superficie de 1 000 m². Se eligieron aquéllos en los que P. patula conformó rodales coetáneos uniespecíficos y que, preferentemente, presentaban alta densidad, condición necesaria para cumplir el requerimiento de competencia extrema por recursos limitados. Adicionalmente, se tuvo el cuidado de que las observaciones de las variables dasométricas cubrieran todas las condiciones de crecimiento e intervalo de edades de la especie estudiada.

En cada sitio de muestreo, todos y cada uno de los árboles fueron medidos en su diámetro normal, el cual se registró en centímetros, con una cinta diamétrica modelo 283D/5m-CSE, también se contabilizó y registró el número de individuos por sitio. A partir de esta información se generó una base de datos y se derivaron, a nivel de hectárea, las variables dasométricas del rodal: densidad o número de árboles (N) y el diámetro cuadrático (Dq), cuyos estadísticos básicos se presentan en el Cuadro 1.

La base de datos final fue auditada, mediante inspección gráfica, para corroborar que las variables de interés presentaran un comportamiento gráfico lógico, en forma de jota inversa (Figura 1); de esta manera, quedó preparada para el ajuste estadístico por regresión de la función tamaño-densidad basada en el índice de densidad de rodal de Reineke.

Figura 1 Comportamiento gráfico de las variables utilizadas en el análisis por regresión de la función tamaño-densidad de Reineke. 

Análisis por regresión estadística

La función tamaño densidad de Reineke expresada en forma lineal, se ajustó por regresión estadística con el método de mínimos cuadrados ordinarios (MCO) lineales, la estructura matemática linealizada del modelo de Reineke es de la forma:

Donde:

Ln = Función logaritmo natural

ε = Término del error que se distribuye como ε = idd N(0, σ²)

Cuando la función se ajustó utilizando la técnica econométrica de regresión frontera estocástica (RFE), la estructura del modelo resultó:

En este caso, el componente de error se divide en: 1) un componente de error que se asume para dar cuenta de la ineficiencia técnica en los datos (u_i); y 2) un componente de error asociado con la medición de las observaciones individuales (v_i).

v_i se asume como una perturbación simétrica distribuida independientemente de u_i, recoge las variaciones aleatorias en la producción debido a factores como errores aleatorios, errores en la observación y de medición de los datos; se asume que se distribuye en la forma:

u_i es un término asimétrico que recoge la ineficiencia técnica de las observaciones y se asume que se distribuye independientemente de v_i y de los regresores.

Bajo estos supuestos, se deben seleccionar distribuciones estadísticas para u_i, que se distribuyan para un solo lado, como en el caso de la modalidad de RFE seminormal y de la exponencial.

Si el valor de la ineficiencia técnica u_i se asume como cero, lo cual es menos probable con valores crecientes de u_i, entonces el modelo idd N⁺(0, σ² _u) se refiere al modelo seminormal.

Si los u_i (i=1..., N) son variables aleatorias no negativas idd N⁺(0, σ² _u), entonces el modelo se conoce como el modelo normal-truncado (en cero) (^{Kumbhakar y Lovell, 2000}; ^{Bi, 2004}; ^{Pretzsch y Biber, 2005}; ^{Zhang et al., 2005}).

Cuando en el ajuste estadístico se usó el enfoque de RFE, se probaron las modalidades de Modelo Semi-Normal (MSN), Modelo Normal Truncado (MNT) y Modelo Normal Exponencial (MNE).

La función tamaño-densidad de Reineke se ajustó por MCO, con el procedimiento Model del paquete estadístico SAS/ETS^® (^{Statistical Analysis System, 2011}) versión 9.3. Cuando se usó RFE, la función se ajustó mediante el procedimiento QLIM del mismo paquete estadístico y versión; en este caso, el algoritmo utiliza la técnica de máxima verosimilitud para estimar la frontera y el parámetro de la ineficiencia técnica de las observaciones (u_i). La técnica de RFE, comparativamente, tiene como ventaja que estima valores extremos (frontera) de un conjunto de datos, en lugar de la media de una función que calcula el método de MCO.

El resultado y la calidad del ajuste por MCO se compararon con los ajustes obtenidos al utilizar RFE, los estadísticos de bondad de ajuste fueron el valor del logaritmo de verosimilitud (logLik), el criterio de información de Akaike (AIC) y el criterio de Schwarz (SchC), a los que, con base en ^{Sakici et al. (2008)} se les aplicó un sistema de calificación y se obtuvo una calificación total para cada modalidad de ajuste.

También, se consideró la significancia de los parámetros y las varianzas de los componentes del error, sobre todo la razón de varianzas de los componentes del error (λ); además de la varianza total (σ) (^{Weiskittel et al., 2009}; ^{Comeau et al., 2010}; ^{Zhang et al., 2013}; ^{Quiñonez-Barraza et al., 2018}). Adicionalmente, para reforzar el criterio de selección, se realizó una comparación gráfica de las líneas de autoaclareo que generó cada una de las técnicas y modalidades de ajuste; principalmente, se observó el comportamiento de la trayectoria de las líneas de autoaclareo sobrepuestas a los datos observados. En el caso del ajuste por MCO, la línea de autoaclareo se determinó de acuerdo con ^{Santiago-García et al (2013)}, en el que se modificó manualmente el valor del intercepto (parámetro α) y se mantuvo fijo el valor de la pendiente (parámetro β), hasta alcanzar un α_máx, que gráficamente se consideró apropiado.

Con los valores de los parámetros del mejor ajuste seleccionado, sobre un gráfico en escala logarítmica, se delimitó la línea de autoaclareo que una hectárea puede soportar sin autoaclareo, y que equivale a un IDR de 100 %, del IDR. Después, con base en la línea de IDR máximo, se estimaron las líneas que definieron cuatro diferentes zonas de crecimiento de Langsaeter, que conforman fajas o bandas de densidades relativas (^{Langsaeter, 1941}; ^{Daniel et al., 1979}; ^{Smith et al., 1997}; ^{Gilmore et al., 2005}) sobre las cuales se deberá gestionar el manejo y crecimiento de los rodales orientados a cumplir el objetivo de producción maderable; de esta manera, queda construido el DMD.

La zona 1 de subutilización se delimitó hasta 25 % con respecto del IDR máximo; la zona 2 de transición se definió entre 25 y 35 %; la zona 3, que corresponde al máximo crecimiento del rodal por hectárea, se localizó entre 35 y 65 %; la zona 4, la de autoaclareo o mortalidad inminente, se ubicó de 65 a 100 %. En la construcción del DMD y para la definición del IDR se utilizó un Dq de referencia de 20 cm.

Calidad de ajuste y determinación de la línea de autoaclareo

El resultado del ajuste estadístico de la función tamaño-densidad de Reineke al utilizar la técnica de mínimos cuadrados ordinarios (MCO) lineales y la de regresión frontera estocástica (RFE), en las modalidades de Modelo Semi-Normal (MSN), Modelo Normal Truncado (MNT) y Modelo Normal Exponencial (MNE) evidencia que, en todos los casos, los parámetros de la función son significativos (Cuadro 2).

Sin embargo, para definir la línea de autoaclareo con la técnica de MCO es necesaria la intervención manual posterior al ajuste, debido a que la línea ajustada corresponde al promedio de la nube de datos observados y no al límite superior absoluto; tal situación ocasiona que la determinación de la referida línea en la frontera superior de los datos sea ineficiente, ya que para modificar el valor del intercepto se utiliza alguna metodología que considera diferentes criterios, como que coincida el número de árboles determinado por el IDR máximo (^{Comeau et al., 2010}); o bien aumentar el valor del intercepto en 1.96 desviaciones estándar del error del modelo (^{Gezan et al., 2007}).

Estos enfoques son prácticos, pero resultan subjetivos para estimar la línea de auto-aclareo (^{Santiago-García et al., 2013}; ^{Quiñonez-Barraza et al., 2018}). Además, comparativamente, la técnica de MCO es sensible al conjunto de observaciones seleccionado a máxima densidad, ya que al ajustar el modelo se puede producir una línea de autoaclareo con pendiente inapropiada (^{Zhang et al., 2005}). Por esas razones, se decidió seleccionar la mejor calidad de ajuste entre las diferentes modalidades probadas de la técnica de RFE, lo que garantiza que se obtenga una línea del máximo absoluto en forma directa e inmediata, y que técnicamente, es el límite correcto (^{Bi, 2004}; ^{Comeau et al., 2010}).

La evaluación comparativa y un sistema de calificación implementado para seleccionar la mejor calidad de ajuste entre las diferentes modalidades al usar RFE revelan que la modalidad de MSN fue superior, ya que presentó los valores más bajos en los estadísticos de ajuste que corresponden al criterio de información de Akaike (AIC), y al criterio de información de Schwarz (SchC); así como con los segundos mejores valores del logaritmo de verosimilitud (logLik), de la varianza total del error (σ) y de la razón de varianzas de los componentes del error (λ), aunados a la alta significancia de los parámetros (cuadro 2 y 3). Adicionalmente, presentó la menor varianza del componente del error de medición (σ_v), y el segundo menor valor en la varianza del componente del error relativo al término de ineficiencia técnica de los datos observados (σ_u).

Este resultado se ratificó mediante una comparación e inspección visual grafica de la línea de autoaclareo que cada ajuste genera (Figura 2), lo cual confirmó que el ajuste por RFE-MSN tuvo un mejor patrón, ya que la línea se ajusta mejor al límite superior de los datos observados, que produjo la mejor definición de la frontera máxima; y por tanto, de la línea de autoaclareo.

Figura 2 Comportamiento gráfico de la línea de autoaclareo generada a partir de la función tamaño-densidad de Reineke por las técnicas de ajuste evaluadas. 

Al respecto, ^{Zhang et al. (2005)} señalan que esta modalidad de ajuste es el método de regresión más apropiado para estimar los coeficientes de la función de Reineke y la posterior línea de autoaclareo, ya que es posible utilizar la totalidad de la información N-Dq disponible de parcelas o sitios que se localizan cerca de un límite superior. ^{Burkhart y Tomé (2012)} refieren que con ello se elimina la necesidad de construir, de manera subjetiva, la base de datos para estimar los valores extremos de la función tamaño-densidad; con ese propósito se ha utilizado en diferentes especies de coníferas y de latifoliadas (^{Weiskittel et al., 2009}; ^{Comeau et al., 2010}; ^{Reyes-Hernandez et al., 2013}; ^{Zhang et al., 2013}; ^{Lopes et al., 2016}; ^{Martínez, 2017}).

La prueba de Wald y de razón de verosimilitud demostraron que el valor de la pendiente (β) del modelo tamaño-densidad, estimado mediante RFE-MSN, es estadísticamente diferente (con α=0.05) del valor teórico propuesto por ^{Reineke (1933)}, de -1.605. Resultado que coincide con lo registrado por ^{Reyes-Hernandez et al. (2013)}, quienes usaron la misma técnica y modalidad de ajuste en rodales de Picea glauca (Moench) Voss de la región de Alberta, Canadá, cuyo valor de la pendiente fue de -1.96. Sin embargo, es diferente a lo determinado por ^{Santiago-García et al. (2013)} para P. patula en la región de Zacualtipán, donde el intervalo de confianza al 95 % para el valor de la pendiente de -1.565 ±0.208 sí incluyó al valor teórico. En ese sentido, el intervalo al 95 % de confianza aquí determinado, para el parámetro en cuestión, fue de -1.998 a -2.341, el cual también es diferente al registrado por ^{Quiñonez-Barraza et al. (2018)} de -1.714 a -1.910, quienes usaron la misma técnica y modalidad para un grupo de especies del género Pinus en Durango, México.

En este contexto, se reafirma el argumento de que la pendiente (β) no siempre está cerca del valor teórico y que puede diferir significativamente entre taxones (^{Comeau et al., 2010}; ^{Santiago-García et al., 2013}), lo cual se debe, en parte, a que distintas poblaciones presentan diferentes tasas de mortalidad, en función de su densidad, hábitos de crecimiento, factores de la productividad del sitio e incluso la edad de la masa (^{Bi et al., 2000}; ^{Weiskittel et al., 2009}; ^{Reyes-Hernandez et al., 2013}). Lo anterior, conduce a la necesidad de desarrollar una alometría particular para cada especie de interés y con ello, evitar que se presenten errores al estimar y controlar la densidad (^{Pretzsch y Biber, 2005}).

El IDR máximo estimado con el ajuste de RFE-MSN, cuando se usó un Dq de referencia de 20 cm, fue de 1 078 árboles, el cual difiere del documentado por ^{Santiago-García et al. (2013)} de 1 662 para la misma especie en la región de Zacualtipán, Hidalgo; en donde determinaron que la mejor técnica de ajuste fue la de RFE en su modalidad de MNT.

El IDR máximo es menor al que se estima en el estudio de ^{Monroy (1997)} para la misma especie en Huayacocotla, Veracruz, de 1 215 árboles con el Dq de referencia de 20 cm, en la que la función fue ajustada en su forma no lineal por MCO, con información de 42 sitios de muestreo temporales. También, el IDR máximo es diferente al consignado por ^{Quiñonez-Barraza et al. (2018)} de 1 603 para el mismo valor del Dq de referencia, en taxa de Pinus en Durango, México, mediante la técnica de RFE en la modalidad de MNT como la mejor para definir la línea de autoaclareo.

Construcción y uso del diagrama para manejar la densidad

Con base en la línea de máxima densidad, que en este caso correspondió a la que se derivó del ajuste bajo el enfoque de RFE-MSN, se construyó el diagrama para el manejo de la densidad de rodales coetáneos y uniespecíficos de P. patula. El número de árboles estimados, para los valores límites, en porcentaje del IDR, que define cada línea que delimita las zonas de crecimiento fueron 270, 380, 700 y 1 078 árboles ha^-1, para 25 %, 35 %, 65 % y 100 %, respectivamente; que a su vez, corresponden a las líneas de crecimiento libre (1), de transición (2), de máximo crecimiento por hectárea (3) y de autoaclareo (4) (Figura 3).

Figura 3 DMD basado en el IDR de Reineke para Pinus patula Schiede ex Schltdl. & Cham. y ejemplo de prescripción de un programa de aclareos. 

A partir del hecho de que el manejo de la densidad es el medio más eficiente de que dispone el silvicultor para alcanzar el objetivo de producción maderable de un rodal, y bajo el supuesto de que la densidad inicial es el punto de partida para prescribir un régimen de aclareos, y que después de practicar un aclareo no existe mortalidad por competencia; por lo que la densidad se mantiene constante hasta un siguiente aclareo.

Como ejemplo de la utilidad práctica y operativa del DMD, se ilustra la forma en que se genera un programa de aclareos como tratamientos silvícolas (cortas intermedias). Se parte de una densidad inicial (N) de 1 050 árboles ha^-1 y se asume un valor teórico de Dq=1 cm, y que para cualquier densidad de un rodal, existe un tamaño máximo promedio en Dq o su equivalente en área basal que puede ser alcanzado, por lo que antes de que se llegue a la zona 4 de mortalidad inminente, se aplica un primer aclareo, en donde se procura que la densidad residual se ubique en el límite inferior de la zona de crecimiento 3 y, en este caso, se alcance un Dq=16 cm.

Desde ese punto, se deja crecer la masa hasta que su crecimiento en Dq, nuevamente, se presente en la vecindad de la línea de autoaclareo, y se requiera la práctica de un segundo aclareo (Figura 3).

Bajo los mismos criterios pueden realizarse aclareos subsecuentes hasta tener el crecimiento máximo posible en Dq, y lograr la densidad óptima deseada. Para maximizar la producción, la densidad de la masa debe mantenerse durante la mayor parte del turno entre la zona 3 y el comienzo de la zona 4.

En el Cuadro 3, para el ejemplo planteado, se presentan los respectivos valores de densidad y su equivalente en área basal por remover y el residual en cada aclareo, además de la intensidad de corta correspondiente. En este caso, en la cosecha final se tendrían 150 árboles con un Dq promedio de 40 cm, que equivaldría a un área basal de 18.8 m² ha^-1.

De manera similar, con el DMD es posible prescribir un programa de aclareo para un rodal en determinada etapa de desarrollo; para ello, con base en un sistema de crecimiento y rendimiento maderable, o bien mediante sitios de muestreo para inventario maderable, se determina su densidad a nivel de hectárea y su respectivo Dq, procediéndose a determinar el IDR en que está el rodal en cuestión; o sobre el DMD se ubica el punto (N, Dq) en función de la zona de crecimiento en que este se localice, y se define la necesidad o no de aplicar un aclareo. Así, conforme a los criterios antes señalados, se analizan diferentes regímenes alternativos de manejo de la densidad y finalmente, se define el más apropiado con base en los objetivos de manejo y producción.

Con relación al primer aclareo, a pesar de que las dimensiones del producto que arroja no son comerciales, debe pensarse que el costo por realizarlo es en realidad una inversión, cuyo retorno y ganancia se reflejará en la cosecha final. Respecto al último aclareo, debe ser selectivo y efectuarse con más cuidado y conocimiento técnico. Hay que tener presente, que a menor calidad de sitio, mayor tiempo le tomará a un rodal alcanzar la línea de autoaclareo.

Para un determinado rodal, se espera que el incremento en volumen del fuste sea cercano al máximo, cuando la densidad relativa esté en la zona (3) (^{Comeau et al., 2010}), condición en la que los árboles aprovechan en forma más efectiva el espacio de crecimiento, sin presentar mortalidad por competencia. Para maximizar crecimiento e incremento en volumen, se requiere manejar la densidad, a través de la aplicación de aclareos, lo que de acuerdo con ^{Pérez y Kanninen (2005)} lleva a estimar un volumen inferior por unidad de superficie, pero con árboles de diámetros más grandes y mejor conformados a edades más jóvenes, lo que conlleva a incrementar su valor.

La frecuencia e intensidad de los aclareos es un compromiso entre las partes silvícola y financiera, esto garantizará la obtención de árboles con las dimensiones y calidad requeridas, de acuerdo al uso industrial, en el mínimo tiempo posible. En la aplicación práctica del DMD, es fundamental el criterio del técnico que realiza el manejo para identificar los árboles que deberán ser removidos con base en los objetivos comerciales establecidos.