Área de estudio

El estudio se realizó en bosques mezclados de los ejidos San Diego de Tezains y Salto de Camellones, pertenecientes a la Unidad de Manejo Forestal (Umafor) 1005 “Santiago Papasquiaro y Anexos”. El área de estudio se ubica en la región noroeste del estado de Durango, entre los 24°30’-25°27’ N y 105°01’-106°24’ O. La altura sobre el nivel del mar fluctúa de 800 a 3 103 m. El clima predominante es templado con lluvias en verano. La precipitación media anual varía de 1 000 a 1 200 mm. La temperatura media anual de 5 a 18 °C, la temperatura más fría ocurre en enero (-6 °C) y las más calurosa en mayo (28 °C) (^{García, 1981}; ^{Quiñonez et al., 2018}). Las comunidades forestales más importantes están constituidas por especies de Pinus, Quercus, Juniperus, Cupressus, Pseudotsuga, Arbutus y Alnus. Aproximadamente 90 % del área de estudio corresponde a masas mixtas con estructuras irregulares (^{Quiñonez et al., 2018}).

Origen de los datos

La información dasométrica utilizada para el ajuste de la relación funcional h-d se obtuvo de 44 unidades o sitios permanentes de muestreo y monitoreo del crecimiento y producción que datan del año 2008, y cubren todos los tipos de vegetación, calidades de sitio, densidades y edades en los rodales en producción maderable. Sus dimensiones son de 50 × 50 m y se establecieron bajo un diseño de muestreo sistemático. Para cada individuo con diámetro normal igual o mayor de 7.5 cm, se registró información de especie, diámetro normal (d, cm), altura total del árbol (h, m), diámetro de copa (m) y su posición respecto al centro de la parcela (azimut, grados; distancia, m).

Se procesó una base de datos con un total de 24 especies. La composición estimada fue de 84 % para el género Pinus (P. arizonica Engelm., P. ayacahuite Ehrenb., P. durangensis Martínez, P. herrerae Martínez, P. lumholtzii B. L. Rob. & Fernald, P. teocote Schiede ex Schltdl. & Cham., P. douglasiana Martínez), 8 % para el Quercus (Q. sideroxyla Humb. & Bonpl., Q. arizonica Sarg., Q. mcvaughii Spellenb, Q. durifolia Seemen ex Loes, Q. crassifolia Bonpl., Q. jonesii Trel., Q. rugosa Neé, Q. laeta Liebm.) y 8 % para los taxones de Juniperus, Cupresus, Arbutus y Alnus (J. deppeanna Steud., J. durangensis Martínez, C. lusitánica Mill., A. arizonica (A. Gray) Sarg., A. bicolor M. González & P.D. Sørensen, A. madrensis M. González, A. tesselata P.D. Sørensen, A. xalapensis Kunth., y Alnus firmifolia Mill). Solo se utilizaron los datos de h-d de los taxa del género Pinus. La distribución de los pares de datos se examinó a través de gráficas, para detectar valores atípicos, los cuales fueron eliminados de la base de datos usada para el ajuste de las ecuaciones; en el Cuadro 1 se muestran los principales estadísticos descriptivos de las variables analizadas.

En cada unidad de muestreo (sitio) se estimaron: el número de árboles por hectárea (N, árboles ha^-1), área basal por hectárea (G, m² ha^-1), diámetro medio cuadrático (dg, cm), altura dominante estimada como el promedio de los 25 árboles de mayor diámetro en el sitio (H ₀, m) (^{Assmann, 1970}), diámetro dominante calculado como el promedio de los 25 árboles de mayor diámetro (D ₀ , cm) (^{Assmann, 1970}) e Índice de Hart (%) obtenido mediante la siguiente expresión:

Modelos

Se analizaron y compararon 10 ecuaciones locales y 10 ecuaciones generalizadas para describir la relación h-d (Cuadro 2). El ajuste fue individual por especie y en grupo (todas las del género Pinus), mediante mínimos cuadrados ordinarios no lineales (ONLS) con el programa R (^{R Core Team, 2016}) por medio de la función “nlsLM” del paquete minpack.lm (^{Elzhov et al., 2012}).

La evaluación de la capacidad de ajuste se basó en gráficos de los residuos, comparaciones numéricas de la raíz del error medio cuadrático (RMSE), valores del coeficiente de determinación (R²), del sesgo promedio de los residuales (Em) y el criterio de información Bayesiano de Schwarz (BIC) (^{Schwarz, 1978}) que compara las ecuaciones en términos de su parsimonia (simplicidad). La formulación de estos estadísticos se expresa de la manera siguiente:

Donde:

yi, yi^ y >y- = Altura observada, estimada y la media observada

n = Número de observaciones

k= Número de parámetros de la ecuación

ln= Logaritmo natural

Para aquellas ecuaciones locales que comparativamente presentaron mayor ventaja en los estadísticos de bondad de ajuste, se analizó la relación entre los parámetros y las variables del sitio (covariables), con la finalidad de explicar mejor la variabilidad entre la altura total y las características del rodal. Con la inclusión de las covariables en el modelo base también se probaron diferentes combinaciones de las variables del rodal y sus transformaciones (cuadrado, logaritmo, raíz, inverso, etcétera) con el fin de mejorar los ajustes. La relación se analizó mediante la función “cor” de R (^{R Core Team, 2016}) y por análisis gráficos.

Modelo de efectos mixtos

Si bien los pares de datos h-d utilizados en el ajuste de las ecuaciones, se obtuvieron de diferentes unidades de muestreo; se detectó una correlación alta entre las mediciones tomadas dentro de una misma unidad, lo que viola uno de los supuestos fundamentales de la teoría de la regresión. Para tomar en cuenta la falta de independencia entre observaciones, se utilizó la técnica de modelos de efectos mixtos, en la cual la variabilidad entre unidades de muestreo se explica mediante la introducción de un parámetro aleatorio, que se estima de forma simultánea a los parámetros fijos (^{Pinheiro y Bates, 1998}; ^{Calama y Montero, 2004}; ^{Castedo et al., 2006}).

El vector de los parámetros de un modelo mixto no lineal puede ser definido de acuerdo con ^{Pinheiro y Bates, (1998)} de la manera siguiente:

Donde:

yij = La j-ésima altura total (árbol) tomada de i-ésima unidad muestral

xij = La j-ésima medición de la variable predictora (d y variables del sitio) tomada de i-ésima unidad muestral

El vector ϕi se puede dividir en componentes fijos y aleatorios. Los fijos son comunes para toda la población y los aleatorios varían de parcela en parcela con la forma:

Donde:

ϕi = Vector de parámetros rx1 (r es el número total de parámetros en el modelo) específico para la j-ésima unidad muestral

f = Función no lineal

εij = Error aleatorio en forma de vector yi=fϕi,xi+εi

λ= Vector px1 de los parámetros fijos (p el número de parámetros fijos)

b_i= Vector qx1 de los parámetros aleatorios asociados con la i-ésima unidad (q es el número de parámetros aleatorios)

A_i y B_i = matrices de tamaño rxp y rxq para los efectos fijos y aleatorios específicos para la i-ésima unidad, respectivamente

Los supuestos teóricos básicos de los modelos mixtos no lineales asumen que el vector de los residuales (εi) y el vector del efecto aleatorio (bi) tienen una distribución normal con media igual a cero y una matriz de varianzas-covarianzas (R_i y D) que representa la variabilidad entre las distintas unidades de muestreo (^{Littell et al., 2006}).

Estimación de los parámetros aleatorios

Los parámetros aleatorios del modelo se calcularon mediante la aproximación marginal de máxima verosimilitud del mejor predictor lineal empírico insesgado (EBLUP) (^{Vonesh y Chinchilli, 1997}), a partir de la expansión de la serie de Taylor con el algoritmo de ^{Lindstrom y Bates (1990)} implementado en la función “nlme” de R (^{R Core Team, 2016}).

La construcción de la ecuación generalizada con inclusión de efectos mixtos se basó en la selección de la ecuación local que mejor bondad de ajuste presentó para todas las especies estudiadas. Así, se probaron diferentes combinaciones de parámetros fijos y aleatorios en cada ecuación, y se compararon los estadísticos RMSE, R ² , Em y BIC, a fin de determinar cuál o cuáles parámetros fijos deberían ser considerados como mixtos. La ecuación h-d generada de esta forma se comparó con la mejor ecuación generalizada.

Predicción de los parámetros aleatorios

Para la aplicación del modelo de efectos mixtos, se necesita información previa de la variable de respuesta (altura total del árbol) que posteriormente se utiliza para calcular los parámetros aleatorios específicos de esa unidad muestral. Las alturas medidas de uno o más árboles en cada sitio son factibles de usarse para predecir los parámetros aleatorios específicos para esa unidad muestral y ajustar la parte fija del modelo de efectos mixtos (respuesta media de la población) a las condiciones del sitio. Esto se conoce como localización o calibración del modelo de efectos mixtos (^{Sharma et al., 2016}).

Cuando se dispone de una submuestra de m_j alturas de árboles, estas pueden utilizarse para predecir los parámetros aleatorios del vector b _j , con la ecuación siguiente (^{Vonesh y Chinchilli, 1997}):

Donde:

D^ = Matriz q × q de varianzas-covarianzas asociada a los parámetros aleatorios (q es el número de parámetros aleatorios incluidos en el modelo), común para todas las unidades y estimado en el ajuste general del modelo

R^i = Matriz m_j × m_j de varianzas-covarianzas del término del error del modelo

ε^i = Vector de los residuales m × 1, cuyos componentes se obtienen por la diferencia entre la altura observada de cada árbol y el valor predicho utilizando el modelo solo con los parámetros fijos

Zi = Matriz m × q de las derivadas parciales de los parámetros aleatorios evaluados en b^i=0

Para determinar el tamaño óptimo de la submuestra de árboles por medir dentro de cada unidad muestral y obtener la respuesta de la calibración (RC), se analizaron las siguientes opciones:

Las seis alternativas descritas se evaluaron en términos del RMSE, R ² y Em y se compararon con los estadísticos obtenidos mediante las ecuaciones ajustadas con el procedimiento ONLS y NLME (i.e. cuando se incluye a todos los árboles de la unidad).

Las predicciones de las alturas de los árboles en las diferentes opciones de calibración evaluadas, se calcularon con el método para estimar los parámetros aleatorios a través de la expresión siguiente:

Donde: h^ij = Altura total estimada (m)

ϕ^i=Aiλ^ λ^ = Matrices para los efectos fijos y aleatorios específicos para la i-ésima unidad muestral, respectivamente

xij = Vector de las variables predictoras

Ecuaciones locales y generalizadas

Para seleccionar la ecuación local y generalizada que mejor bondad de ajuste presentó de manera individual para cada taxon, en primera instancia se ajustaron todas las ecuaciones del Cuadro 2 por especie. Durante esta etapa, se identificaron problemas de convergencia y parámetros no significativos para algunas de ellas (Pinus ayacahuite, P. herrerae, P. lumholtzii y P. douglasiana) con las ecuaciones (10), (13) y (17); por lo tanto, se decidió realizar los análisis posteriores con la base de datos conjunta para todas los taxa.

Al comparar los valores de los estadísticos de bondad de ajuste del grupo de ecuaciones locales (Cuadro 3), las ecuaciones (3), (8), (9) y (10) resultaron ser las más precisas y con el menor sesgo. Si bien, proporcionaron valores muy similares en los estadísticos de comparación utilizados, la Ecuación (3) tuvo parámetros no significativos en el ajuste, aun cuando el valor del BIC fue ligeramente menor al resto de las ecuaciones. Finalmente, los análisis gráficos del sesgo en la estimación de la altura por clases diamétricas indicaron que la Ecuación (8) evidenció el mejor comportamiento en todo el intervalo de datos, por lo que fue seleccionada como el modelo base para estudiar la inclusión de variables de rodal en su estructura (Figura 1a).

Figura 1 Comportamiento del sesgo por clase diamétrica de las ecuaciones locales y generalizadas con mejor bondad de ajuste.  

^{Huang et al. (2000)} y ^{Peng et al. (2001)} señalan que la Ecuación (8) de Chapman-Richards es quizás la más flexible y versátil para modelar la relación h-d en rodales forestales. Al respecto, ^{Uzoh (2017)} demostró es estadísticamente la más robusta (i.e. menor sesgo) de entre 12 distintas opciones para estimar la altura de Pinus ponderosa Douglas ex Lawson en un modelo de efectos mixtos en el oeste de los Estados Unidos de América.

Al comparar los valores de los estadísticos de bondad de ajuste del grupo de ecuaciones generalizadas (Cuadro 3), se advierte que las ecuaciones (11) y (13) son las más precisas, al considerar los estadísticos R ² y RMSE. Sin embargo, la Ecuación (13) exhibió parámetros no significativos, aun cuando resulta el valor más bajo del BIC. Por lo tanto, la ecuación generalizada de ^{Sharma y Parton (2007)} (11) presentó mejores ventajas en precisión y parsimonia (BIC) del total evaluadas; además, los análisis gráficos del comportamiento del sesgo al estimar las alturas por clases diamétricas tuvieron un comportamiento insesgado (Figura 1b). Esta ecuación es bastante flexible para describir la relación h-d, y ha sido ampliamente usada en estudios similares (^{Vargas et al., 2009}; ^{Crecente et al., 2010}; ^{Corral et al., 2014}).

Los mejores en términos de precisión y parsimonia resultan de las ecuaciones generalizadas (Cuadro 3), debido a que incluyen características del rodal como variables predictoras, con lo que disminuye la variabilidad dentro de cada unidad muestral, y, en general, son ampliamente recomendadas para describir la relación h-d (^{López et al., 2003}; ^{Barrio et al., 2004}; ^{Canga et al., 2007}; ^{Vargas et al., 2009}; ^{Crecente et al., 2010}; ^{Corral et al., 2014}; ^{Uzoh, 2017}).

El análisis de correlación entre variables del rodal con los parámetros de la Ecuación (8) (Cuadro 4) indicó que solo tres de ellas (H ₀, N y G) se correlacionaron de manera significativa para mejorar su precisión, de lo que resultó la ecuación generalizada con la expresión siguiente:

Donde:

b₀ - b₄ = Parámetros de la ecuación; el resto de las variables previamente definidas

Donde:

H ₀ = Altura dominante (m)

D ₀ = Diámetro dominante (cm)

N = Densidad (número de árboles ha^-1)

G = Área basal total (m² ha^-1)

d_g = Diámetro medio cuadrático (cm)

IH = Índice de Hart (%)

b_i = Parámetros estimados

Con la inclusión de la altura dominante, el área basal y la densidad del rodal, además del diámetro normal en la Ecuación (28), la RMSE se redujo 18 % con respecto al modelo base de Chapman-Richards (Ecuación 8) (de 2.95 m a 2.43 m) y, por consecuencia, el valor de la R ² se incrementó 12 % (de 0.731 a 0.818). El valor del BIC se redujo 8 % (de 20 200 a 18 651) dado que en su expresión considera la suma total de los cuadrados residuales.

Por otra parte, al comparar la precisión y parsimonia de la Ecuación (28) respecto a la Ecuación (11) desarrollada por ^{Sharma y Parton (2007)}, esta también tuvo una mejoría marginal, ya que la RMSE disminuyó en un 0.4 % y el valor de la R ² aumentó 0.16 %, en tanto que BIC tuvo una reducción de 1.5 %.

Con base en los estadísticos de bondad de ajuste y en el comportamiento gráfico del sesgo por clase diamétrica (Figura 1b), la Ecuación (28) se consideró el modelo base para ajustarla mediante la técnica de efectos mixtos en la relación h-d de las siete especies estudiadas. Dicha expresión incluye a la altura dominante como un indicador de la productividad, al área basal como una medida simple y objetiva de la ocupación de los árboles y que en general es mejor que el número de árboles por hectárea, ya que combina este con su tamaño (^{Huang et al., 2009}), y el número de árboles genera predicciones sujetas a cambios instantáneos de acuerdo con los aclareos o intervenciones silvícolas a los que son sometidas las masas (^{Crecente et al., 2010}). Esas variables las han utilizado en la modelación de la relación h-d diversos autores, para una diversidad amplia de taxones (^{López et al., 2003}; ^{Calama y Montero, 2004}; ^{Canga et al., 2007}; ^{Crecente et al., 2010}; ^{Corral et al. 2014}; ^{Sharma y Breidenbach, 2015}).

Modelo de efectos mixtos

En el ajuste de la Ecuación (28), se observó que los parámetros con mayor variabilidad fueron los que definen a la asíntota y la forma de la curva h-d (Cuadro 5). Por lo tanto, en un primer paso se ajustó con los parámetros que definen la asíntota (b₀ , b₁ ) y la forma (b₃ , b₄ ) como aleatorios; sin embargo, al no obtener convergencia, se probaron varias combinaciones con uno y dos parámetros aleatorios asociados a los parámetros fijos hasta lograr la convergencia.

La mejor combinación resultó al asociar el parámetro de forma (b₃ ) con uno aleatorio u_i de forma aditiva (b₃ +u_i ). Resultados similares consignan ^{Sharma y Parton (2007)}, ^{Vargas et al. (2009)} y ^{Corral et al. (2014)}, quienes consideraron solo un parámetro aleatorio dentro del modelo mixto.

El modelo final generalizado con efectos mixtos se expresa de la manera siguiente:

Donde:

b₀ - b₄ = Parámetros fijos (comunes a todas las unidades muestrales)

uj~N0,σu = Parámetro aleatorio (específico para cada unidad muestral) con una distribución normal con media igual a cero y varianza debida al efecto aleatorio (σu2)

hij y εij = Altura total y error estimados por el modelo para la i-ésima observación (árbol) de la j-ésima unidad muestral.

El modelo de efectos mixtos (Ecuación 29) registró una mejora en precisión de 6.6 % en el valor de la RMSE, con respecto al modelo base (Ecuación 28); además, con la inclusión del parámetro aleatorio entre las unidades muestrales, se logró un incremento de 2.8 % sobre el valor de R² , al explicar la variabilidad observada. El valor del sesgo promedio en ambas ecuaciones fue similar; pero, el BIC, que evalúa la simplicidad del modelo, se redujo 2.9 % (Cuadro 5).

La Figura 2a ilustra la dispersión de los residuos frente a las alturas predichas con la Ecuación (29); se observa una varianza homogénea en todo el intervalo de datos, además de que no existe un patrón sistemático en la variabilidad de los residuos. Por lo tanto, el efecto aleatorio ayudó a eliminar la heterogeneidad en la varianza (^{Fang y Bailey, 2001}). El respectivo gráfico de residuos frente al primer retraso (Figura 1b), no evidencia una tendencia que muestre autocorrelación, lo que indica que se le removió cuando el efecto aleatorio entre unidades muestrales fue modelado al expandir el parámetro b ₃.

Figura 2 Residuales frente a las alturas predichas (a), residuales frente al primer retraso (b) estimados con la Ecuación 29. 

Calibración del modelo

Los resultados de las alternativas de calibración para mejorar las predicciones se muestran en el Cuadro 6. Los valores más altos de la RMSE se obtuvieron al incluir una submuestra de uno a tres árboles con categorías diamétricas extremas (mínimo y máximo). En contraste, con la inclusión de individuos pertenecientes a categorías diamétricas centrales a la distribución, se identificó una ligera reducción en la RMSE. Pese a que no existe una metodología para la definición del tamaño de muestra más apropiado en la calibración de un modelo de efectos mixtos, se advierte que la alternativa RC6 resultó con más precisión, cuando se mide la altura total de una submuestra de tres árboles seleccionados al azar con diámetros cercanos al segundo cuartil de la estructura diamétrica del sitio, dado que presentó una reducción de 4.6 % en la RMSE respecto al valor estimado durante el ajuste de la Ecuación (28) por ONLS. Este valor fue muy cercano al generado por el ajuste del modelo de efectos mixtos (medición de la altura total de todos los árboles) de la Ecuación (29) por el procedimiento NLME.

El modelo de efectos mixtos registró un incremento en la precisión, debido a que incluye la variabilidad que existe en la altura total de los árboles dentro y entre parcelas (^{Calama y Montero, 2004}). En virtud de que los bosques mezclados tienen una variabilidad mayor en la altura total de los árboles, con respecto a aquellos que provienen de rodales con una o dos especies, se estima que las predicciones del modelo de efectos mixtos desarrollado son aceptables, ya que explican 82 % de la variabilidad de los datos. En experiencias con datos correspondientes a masas regulares se ha logrado una capacidad de predicción cercana a 90 % (^{Calama y Montero, 2004}; ^{Diéguez et al., 2005}; ^{López et al., 2012}).

La calibración de los modelos de efectos mixtos es una ventaja para los manejadores forestales, ya que a partir de pocos árboles es posible determinar la altura de todos los individuos presentes en el rodal. ^{Calama y Montero (2004)} proponen cuatro árboles seleccionados aleatoriamente; ^{Corral et al. (2014)} definen que tres árboles cercanos a 10 % del diámetro promedio son suficientes; mientras que ^{Castedo et al. (2006)} y ^{Crecente et al. (2010)} proponen los tres árboles más pequeños de la estructura diamétrica del rodal. A medida que sea mayor el número de árboles seleccionados para calibrar un modelo de efectos mixtos, por lo general, se reduce la RMSE y mejora su capacidad predictiva; sin embargo, lo anterior no siempre se justifica debido a que se incrementa el costo de muestreo, por lo que la calibración con una muestra pequeña de árboles es la principal utilidad de los modelos de efectos mixtos para su aplicación en el campo forestal (^{Castedo et al., 2006}).