INTRODUCCIÓN

Drimys winteri J.R et G. Forster, es una especie arbórea que posee una gran plasticidad ecológica que le permite crecer en la mayor parte de los bosques naturales de Chile, entre los 30° y los 50° S, tanto en la cordillera de la Costa como en la de Los Andes, preferentemente ocupando terrenos húmedos (^{Larjavaara y Muller-Landau 2013}). Drimys winteri es una especie pionera que coloniza e invade terrenos desnudos después de cortas e incendios especialmente en la región sur-austral de Chile formando importantes bosques secundarios (^{Corvalán et al. 1987}a). Específicamente en la isla de Chiloé y la provincia de Llanquihue en el sur de Chile, D. winteri es una de las especies nativas de mayor interés productivo, formando bosques secundarios de un alto potencial de manejo silvícola (^{Corvalán et al. 1987}b, ^{Navarro et al. 1999}).

La relación altura-diámetro de los árboles es muy importante en ecología y manejo de bosques. Esta relación es el resultado de diversos factores que afectan el crecimiento de un árbol, como por ejemplo, el sitio donde crece y la competencia a la cual se ve expuesto. El error y el alto costo de medición de las alturas de los árboles es un tema clásico de la mensura forestal que se resuelve utilizando modelos de regresión empleando el diámetro (d) como variable predictora. Los errores de medición son especialmente comunes en rodales de alta densidad como en los renovales de D. winteri (^{St-Onge et al. 2004}, ^{Larjavaara y Muller-Landau 2013}), por lo que es conveniente desarrollar modelos de altura-diámetro (h-d) a base de mediciones obtenidas en una submuestra de árboles (^{Salas y Real 2013}, ^{Mehtätalo et al. 2015}).

La mayoría de los modelos de altura-diámetro para D. winteri han sido ajustados en forma local y a base de muestras restringidas. Los modelos de altura-diámetro se han clasificado en: “locales” a los que solamente ocupan como variable predictora a la variable diámetro del árbol; y “generalizados” a los que incluyen además alguna(s) variable(s) de rodal (^{Temesgen et al. 2007}, ^{Mehtätalo et al. 2015}). Todos los modelos de h-d para D. winteri reportados en la literatura son locales y se han desarrollado aplicando modelos estadísticos tradicionales, i.e., regresión lineal mediante el método de los mínimos cuadrados (^{Corvalán et al. 1987}b, ^{Navarro et al. 1999}, ^{Reyes et al. 2009}). Los modelos locales de h-d solo son apropiados para predicciones en condiciones de bosque y sitio restringidas. Por el contrario, los modelos generalizados de h-d, permiten su uso en diferentes condiciones de bosque y sitios, lo cual en el caso de D. winteri, que es una especie que crece naturalmente en rodales densos y con una gran amplitud de calidad de sitios, los hace particularmente de interés dado su plasticidad ecológica. Sin embargo, a la fecha no existen estudios al respecto para esta especie.

El efecto de variables de rodal y/o sitio en la modelación de la relación altura-diámetro ha sido deficientemente evaluado desde un punto de vista estadístico. Cuando un modelo generalizado de h-d se ajusta mediante el método de los mínimos cuadrados (e.g., en ^{Staudhammer y LeMay 2000}, ^{Schröder y Álvarez 2001}, ^{Trincado y Leal 2006}), no se considera la estructura jerárquica de los datos (i.e., árboles anidados a unidades de muestreo), induciendo por lo tanto una correlación que debe ser tomada en cuenta al momento de estimar los parámetros de los modelos y sus respectivos estimadores de varianza. Esto trae problemas de pseudoreplicación (^{Hurlbert 1984}), que no han sido siempre bien reconocidos en la literatura. Una alternativa estadísticamente correcta para representar el efecto de las variables de rodal (e.g., densidad, área basal) y del sitio (e.g., índice de sitio, o nivel de crecimiento medio) en un modelo h-d, es ajustarlo en un contexto de efectos mixtos (^{Mehtätalo 2004}, ^{Robinson y Wykoff 2004}, ^{Mehtätalo et al. 2015}). El objetivo del presente estudio es modelar el efecto de variables de rodal y sitio en la alometría altura-diámetro de D. winteri, mediante un enfoque de efectos mixtos. Los modelos que incluyen estas variables pueden ser de gran apoyo para la toma de decisiones silviculturales tales como la selección de sitio y el control de la competencia.

MÉTODOS

Datos

Se establecieron 122 unidades de muestreo (i.e., parcelas) en bosques no-intervenidos de D. winteri en las provincias de Chiloé y Llanquihue en el sur de Chile (figura 1). Dichas unidades fueron establecidas entre los 41° y 43° S por ^{Corvalán et al. (1987}a) y ^{Donoso et al. (2007}). La superficie de las parcelas varía entre 100 y 400 m². En cada unidad de muestreo se midieron variables dendrométricas a todos los árboles con un diámetro a la altura del pecho (d) mayor a 5 cm, tales como: d, especie, posición sociológica, y otras variables. Además, a una submuestra aleatoria de árboles por parcela, se les midió la altura total (h) con hipsómetro. El tamaño de la submuestra de altura por parcela fue variable, oscilando entre 3 y 43 árboles. Además, a una submuestra de tres o cinco árboles dominantes (que pertenecen al dosel superior) de cada parcela se les determinó la edad, mediante análisis fustal (^{Corvalán et al. 1987}a) o extracción de un tarugo de incremento (^{Donoso et al. 2007}), respectivamente. La dispersión entre la altura y el diámetro muestra un amplio rango de variación de la especie en el área de estudio, demostrando la plasticidad ecológica de la especie (figura 2).

Figura 1 Área de estudio y distribución de las unidades de muestreo (i.e., parcelas) en bosques de Drimys winteri en las provincias de Llanquihue y Chiloé, sur de Chile. Nótese que debido a la escala del mapa, existen parcelas que están traslapadas 

Figura 2 Dispersión entre la altura y el diámetro para árboles de Drimys winteri en las provincias de Llanquihue y Chiloé, en el sur de Chile 

Con la información medida en cada unidad muestral, se calcularon variables a nivel agregado o de estado de rodal, como la densidad (N, número de árboles), área basal (G), diámetro del árbol de área basal media ( ^_dg ), según lo propuesto por ^{van Laar y Akça (2007}). La altura dominante ( ^_Hd ) fue estimada mediante el algoritmo de ^{García y Batho (2005}), el cual toma en consideración la superficie de las parcelas, lo cual afecta su estimación. A base de las mediciones de los tarugos de incremento de los árboles dominantes se determinó la edad dominante de la unidad de muestreo ( ^_Ed ), como el promedio de las edades de los árboles dominantes. Finalmente, se obtuvo el incremento medio anual en altura dominante ( ^_imaHd ), que es el cuociente entre ^_Hd y ^_Ed , como una aproximación a la productividad del sitio (^{van Laar y Akça 2007}). El amplio rango de las variables a nivel individual y de parcela (cuadro 1), así como también su gran variabilidad demuestra lo exhaustivo del muestreo en terreno.

Cuadro 1 Estadística descriptiva para las principales variables a nivel de árbol y de rodal. Las variables de árbol son: diámetro (d) y altura total (h), mientras que las de rodal son densidad (N), área basal (G), diámetro del árbol de área basal media ( ^_dg ), edad dominante ( ^_Ed ), e incremento medio anual en altura dominante ( ^_imaHd ). Descriptive statistics of tree- and stand-level variables. Tree-level variables are: diameter (d) and total height (h); meanwhile stand-level variables are: tree density (N), basal area (G), diameter of the mean basal area tree ( ^_dg ), dominant age ( ^_Ed ), and mean annual increment of dominant height ( ^_imaHd ). 

Modelación

Primero se seleccionó una ecuación base de altura-diámetro, y en segunda instancia, se empleó dicha ecuación para ajustarlo mediante modelos de efectos mixtos. Para seleccionar la ecuación base, se ajustaron una serie de modelos no-lineales, dentro de los cuales destacan el de ^{Richards (1959}) y el de Michaelis-Menten (^{Bates y Watts 1980}, ^{Mehtätalo et al. 2015}) mediante el método de mínimos cuadrados no-lineales. Luego de una evaluación predictiva y de las características biológicas de estos modelos, i.e., que al menos uno, o la mayoría, de sus parámetros tenga una interpretación directa dimensional (e.g., la asíntota, ver anexo), se seleccionó el modelo de Michaelis-Menten, que se expresa como

[1]

donde: h 𝑖 y 𝑑 𝑖 es la altura total y diámetro a los 1,3 m sobre el fuste para la i-ésima observación, respectivamente; α y β son parámetros; y 𝜀 𝑖 es el i-ésimo error aleatorio del modelo con distribución Gaussiana con media 0 y varianza 𝜎 𝜀 2 .

Este modelo ha sido ampliamente empleado para estudios de cinética enzimática en reacciones químicas y en estudios del crecimiento de plantas (^{Kobe 1999}). Una de las particularidades de la ecuación [1] es que posee una interpretación biológica de sus parámetros, donde α representa a la asíntota (i.e., máxima altura) y β al diámetro donde se alcanza la mitad de la asíntota. En el anexo se detalla una demostración matemática al respecto.

En una segunda etapa, y empleando el modelo de Michaelis-Menten como base se ajustaron diferentes modelos de efectos mixtos, en donde se agregaron efectos aleatorios a cada uno de los parámetros-coeficientes de ecuación [1]. Los modelos de efectos mixtos permiten representar la estructura jerárquica de los datos (^{Pinheiro y Bates 2000}, ^{Trincado et al. 2007}, Salas et al. 2007), ya que las alturas individuales son medidas dentro de unidades de muestreo, por lo tanto, existen datos anidados. Se evaluaron diferentes variantes de modelos de efectos mixtos, al agregar efectos aleatorios a uno y/o a todos o diferentes combinaciones de estos. El ajuste se realizó mediante máxima verosimilitud empleando el paquete nlme (^{Pinheiro et al. 2017}) de R (^{R Core Team 2017}). Luego de comparar los estadísticos de bondad de ajuste entre las diferentes variantes de modelos de efectos mixtos (cuadro 2), tales como AIC, BIC y test de razón de verosimilitud (^{Pinheiro y Bates 2000}), se seleccionó el siguiente modelo de efectos mixtos como la mejor variante

[2]

donde: ℎ 𝑖𝑗 y 𝑑 𝑖𝑗 es altura y diámetro para el i-ésimo árbol en la j-ésima parcela, respectivamente; α y β son parámetros; 𝑎 𝑗 y 𝑏 𝑗 son efectos aleatorios a nivel de parcela y 𝜀 𝑖𝑗 es el error aleatorio del modelo.

El modelo [2] se ajustó además considerando una varianza de los residuales diferente para cada clase de la variable incremento medio anual (i.e., cuatro clases definidas por los cuartiles de la variable, ima _Hd ) y fuente de datos, tal como lo aplicaron ^{Salas et al. (2008}).

Finalmente, se modelaron los efectos aleatorios en función de variables de estado de rodal que representan la densidad de los bosques (e.g., N, G, d _g ) y la calidad del sitio (e.g., ima _Hd ), debido a su potencial efecto en la alometría altura-diámetro. Para tales efectos, se obtuvieron parámetros del modelo de altura-diámetro a nivel de parcela como sigue,

[3]

[4]

donde: 𝛼 y 𝛽 son parámetros estimados; 𝑎 𝑗 y 𝑏 𝑗 son efectos aleatorios predichos por el modelo de efectos mixtos para la j-ésima parcela.

Se evaluaron diferentes modelos para predecir 𝛼 𝑗 y 𝛽 𝑗 de la forma general,

[5]

[6]

donde: 𝑓( ) es un modelo lineal o no-lineal; 𝜃 y ∅ son vectores de parámetros a ser estimados; 𝑋 𝛼 y 𝑋 𝛽 son matrices de variables predictoras para el modelo [5] y [6], respectivamente; y 𝜀 𝛼 𝑗 y 𝜀 𝛽 𝑗 son vectores de errores aleatorios para el modelo de 𝛼 𝑗 y 𝛽 𝑗 , respectivamente.

Estos modelos (ecuaciones [5] y [6]) de efectos aleatorios se denominaron sub-modelos, y se ajustaron mediante mínimos cuadrados, evaluando diferentes potenciales variables predictoras que puedan representar el efecto de la competencia y de la calidad del sitio de los bosques donde se desarrollan los árboles. La selección de los mejores sub-modelos, contempló evaluar el efecto de su predicción en el comportamiento del modelo de altura-diámetro. Para esto se desarrollaron análisis gráficos de comportamiento bajo diferentes niveles de las variables predictoras. Para todos los modelos y sub-modelos ajustados, se calcularon estadísticos de capacidades predictivas como la raíz cuadrada media de las desviaciones, la diferencia agregada, y diferencia absoluta (^{Salas et al. 2010}).

RESULTADOS

La dispersión de la relación altura-diámetro muestra la alta variabilidad de D. winteri creciendo en condiciones naturales, así como demuestra la plasticidad ecológica de la especie o la variabilidad genética (figura 2). La mejor variante de modelos de efectos mixtos fue la que considera agregar efectos aleatorios a los dos parámetros de la ecuación de Michaelis-Menten (cuadro 2), es decir un modelo completo (ecuación [2]). Tal como lo sugieren ^{Pinheiro y Bates (2000}), se evaluó la sobre-parametrización del modelo de efectos mixtos al analizar la relación entre los efectos aleatorios, no encontrándose evidencia de esta.

Cuadro 2 Estadísticos de bondad de ajuste basados en máxima verosimilitud para tres variantes del modelo de efectos mixtos de la ecuación altura-diámetro de Michaelis-Menten (ecuación [1]). La variante α es cuando los efectos aleatorios son agregados a dicho parámetro del modelo base, y así sucesivamente para las otras dos variantes. Nótese que la variante α y β es representada en la ecuación [2]. AIC es el criterio de información de Akaike y BIC es el criterio de información Bayesiano. Maximum likelihood based statistics of the height-diameter Michaelis-Menten (equation [1]), for three mixed-effects model variants. The α variant add random effects into that parameter, and so on for the other variants. Notice that the α and β variant is represented by equation [2]. AIC is the Akaike’s information criterion and BIC is the Bayesian information criterion. 

La estructura del modelo de efectos mixtos propuesto, implica un efecto de las características de los bosques donde crecen los árboles en la asíntota del modelo (α en ecuación [1]) y en el parámetro que determina el diámetro cuando se alcanza la mitad de la asíntota (β en ecuación [1]). Así también el modelo propuesto toma en consideración la estructura jerárquica de los datos, por lo tanto, estimando eficientemente las varianzas de los parámetros estimados del modelo, lo cual es clave para la correcta inferencia estadística.

Luego de analizados una gran variedad de sub-modelos para 𝛼 𝑗 y 𝛽 𝑗 se seleccionó una misma relación funcional para ambos parámetros, la cual fue un modelo lineal del siguiente tipo

[7]

[8]

donde: 𝑖𝑚𝑎 𝐻𝑑 𝑗 , 𝑑 𝑔 𝑗 y 𝑁 𝑗 son el incremento medio anual en altura dominante, diámetro del árbol de área basal media, y densidad para la j-ésima parcela, respectivamente; α₀, ..., α₄ y β₀, ..., β₄ son parámetros, y 𝜀 𝛼𝑗 y 𝜀 𝛽𝑗 son errores aleatorios de los sub-modelos de 𝛼 𝑗 y 𝛽 𝑗 , respectivamente.

El ajuste de estos sub-modelos fueron estadísticamente significativos (P < 0,05), y con errores medios de 2,38 (m) y 1,56 (cm) para 𝛼 𝑗 y 𝛽 𝑗 , respectivamente (cuadro 3).

Cuadro 3 Parámetros estimados del modelo de altura-diámetro de efectos mixtos (ecuación [2]) y sub-modelos para predecir α _j (ecuación [7]) y β _j (ecuación [8]). Estimated parameters of the mixed-effects height-diameter model (equation [2]) and sub-models α _j (equation [3]) and β _j (equation [4]). 

Estos sub-modelos implican que la competencia y la calidad del sitio donde crecen los árboles afectan en forma simultánea a la relación alométrica altura-diámetro. La competencia es representada, en el modelo propuesto, por la densidad y el diámetro medio de los bosques, mientras que la calidad del sitio, por el incremento medio anual en altura dominante. Según los análisis de comportamiento, mayores densidades producen relaciones altura-diámetro de mayor magnitud, es decir, para un mismo diámetro de referencia, un árbol puede alcanzar una altura mayor (figura 3).

Figura 3 Comportamiento del modelo propuesto de altura-diámetro para Drimys winteri según diferentes niveles de competencia (representado por la densidad, N) y productividad de sitio (representado por el incremento medio anual en altura dominante, imaHd). Para las tres figuras se empleó un diámetro del árbol de área basal media (dg) igual a 15 cm. La curva “poblacional” se refiere al valor esperado del modelo sin considerar los efectos de la competencia y sitio 

Los estadísticos de validación del modelo propuesto, es decir: predicción de α y β de la ecuación [1] al emplear los valores respectivos para las variables de rodal y sitio de la j-ésima parcela, y luego calculó del valor predicho para el i-ésimo árbol dentro de la j-ésima parcela, muestran una disminución en relación al ajuste del modelo base mediante mínimos cuadrados. En particular la raíz cuadrada media de las diferencias (RMSD) disminuye un 6 % (cuadro 4).

Cuadro 4 Estadísticos de capacidades predictivas para modelo base (ecuación [1]) y modelo de efectos mixtos (ecuación [2]), con efectos aleatorios predichos según los sub-modelos αj (ecuación [7]) y βj (ecuación [8]). RMSD representa la raíz cuadrada media de las diferencias, DA a la diferencia agregada (o error medio), y DAA a la diferencia media absoluta. Statistics of prediction capabilities for the base model (equation [1]) and the mixed-effects model (equation [2]), with predicted random-effects from submodels αj (equation [3]) and βj (equation [4]). RMSD is the root mean square differences; DA is the aggregated difference (or average error); and DAA the aggregated of the absolute value differences.  

DISCUSIÓN

Las alturas y diámetros reportadas en esta investigación, se encuentran dentro del rango esperado para el área de estudio, con alturas máximas de 25 m y diámetros inferiores a 100 cm (^{Donoso et al. 2006}). Si bien los valores máximos de estas variables son menores a los reportados para otras especies del bosque nativo chileno, por ejemplo, para especies del género Nothofagus (^{Salas y Real 2013}), son valores típicos para D. winteri en bosques de segundo crecimiento (^{Corvalán et al. 1987}b, ^{Navarro et al. 1999}, ^{Soto y Donoso 2006}). La gran variabilidad de las alturas que D. winteri puede alcanzar para un mismo diámetro demuestra la plasticidad de la especie para crecer en diferentes condiciones de sitio y de competencia.

En general, para representar el efecto de variables de rodal y/o de sitio en modelos de altura-diámetro, se han ajustado modelos mediante el método de mínimos cuadrados (^{Staudhammer y LeMay 2000}, ^{Trincado y Leal 2006}), lo que tiende a sobreestimar el efecto o la importancia estadística de las variables de rodal debido a la pseudoreplicación (^{Hurlbert 1984}). En cambio, en el presente estudio, se analizó el efecto de variables de rodal mediante modelos de efectos mixtos, lo cual permite asignar eficientemente el efecto de dichas variables en la jerarquía correspondiente. Hubo una mejora en el ajuste empleando efectos mixtos, donde los estadísticos de bondad y de validación disminuyeron en comparación al ajuste por mínimos cuadrados. Así también, al definir diferentes escenarios y condiciones para las variables predictoras (figura 2), se logró cubrir el rango de variabilidad de la alometría altura-diámetro.

La relación altura-diámetro de D. winteri se ve influenciada por dos grandes factores: la competencia y la calidad de sitio. El efecto de la competencia ha sido representado mayormente a nivel agregado por el número de árboles (^{Parresol 1992}, ^{Zeide y Vanderschaaf 2002}, ^{Soares y Tomé 2002}) y el área basal (^{Parresol 1992}, ^{Staudhammer y LeMay 2000}), y se basa en el hecho que los árboles que crecen en rodales más densos tienden a desarrollarse más en altura que en diámetro, para así competir de mejor forma por la luz. Por lo tanto, esto altera la alometría altura-diámetro, y el modelo propuesto (que depende del número de árboles) así lo representa. El efecto del tamaño medio de los árboles en un bosque, representado por dg, tiene que analizarse en relación al diámetro del árbol a ser evaluado. Esto es, si el árbol tiene un diámetro menor al dg del bosque, representaría a un árbol que está sometido a mayor competencia que un árbol de área basal media. Por lo tanto, debería considerarse su efecto en el presente modelo, como un tipo de índice de competencia a nivel individual (^{van Laar y Akça 2007}). En los estudios de ^{Mehtätalo (2004}), ^{Trincado y Leal (2006}) y ^{Mehtätalo et al. (2015)} se han empleado con éxito también el dg como variable de rodal en modelos de h-d. Considerando todo lo anterior, el modelo de h-d propuesto considera el efecto de la competencia a nivel agregado y a nivel individual, con variables de fácil obtención.

El segundo factor, la calidad del sitio, en modelos h-d ha sido representada casi exclusivamente por la altura dominante (Hd) del bosque, variable de rodal que se debe calcular (o estimar) a base de mediciones (o estimaciones) de alturas individuales que normalmente provienen de un modelo h-d, por lo tanto, su uso es circular en el contexto del presente estudio. Es por esto que en el presente estudio deliberadamente se excluyó emplear la variable Hd. Como indicador indirecto de la productividad de un sitio, se empleó al imaHd, que incluso puede ser categorizado en sitios buenos, regulares y malos, para su uso más práctico. Debido a la mayor productividad del sitio, producto de mejores condiciones edáficas, topográficas y climáticas, existe una mayor disponibilidad de recursos para el crecimiento de un árbol, lo cual se traduce en un mayor desarrollo en altura (i.e., efecto positivo). Es importante destacar que el efecto combinado de ambos factores (i.e., competencia y sitio), implica que la competencia produce árboles de mayor altura en sitios de mayor productividad (^{Pretzsch 2009}).

Finalmente, los resultados acá presentados proveen modelos de altura para diferentes sitios cubiertos por bosques secundarios dominados por D. winteri, y adicionalmente consideran el efecto de la densidad sobre el crecimiento en altura dentro de sitios de la misma calidad. Este es el primer estudio en Chile que evidencia que la altura de los árboles dominantes es influida por la densidad, con una relación negativa como es de esperar, ya que a menor densidad los árboles pueden asignar menos carbohidratos al crecimiento en altura en favor del desarrollo de las copas y los fustes (^{Kramer y Koslowski 1960}, ^{Lambers et al. 2008}). Evidencias de este tipo han sido reportadas en otras especies, como por ejemplo en Pseudotsuga menziesii (Mirbel) Franco, en Estados Unidos (^{Salas et al. 2008}). Los modelos acá provistos en consecuencia constituyen una herramienta de gran valor para mejorar las estimaciones de crecimiento y productividad de los bosques secundarios dominados por D. winteri en Chile, los cuales constituyen una reserva de bosques de alto potencial para el manejo y para el desarrollo forestal especialmente entre los 41° y 43° S (^{Navarro et al. 1999}, ^{Navarro y Cabello 2018}).

CONCLUSIONES

La modelación de la alometría altura-diámetro mediante un enfoque de modelos no-lineales de efectos mixtos, permitió representar satisfactoriamente su comportamiento bajo diferentes factores que la afectan. Estos factores son la competencia y la calidad del sitio, los cuales fueron representados por variables predictoras, cuya significancia estadística fue correctamente estimada, dado que el modelo toma en consideración la estructura jerárquica de los datos. Además, el modelo propuesto posee una interpretación ecológica coherente con el conocimiento actual del crecimiento de árboles y, ofrece mejores capacidades predictivas que un modelo de altura-diámetro local.