Ajuste de curvas de índice de sitiomediante modelos mixtos para plantaciones de eucalyptus urophylla en Venezuela.

Omar Carrero, Mauricio Jerez, Raúl Macchiavelli, Giampaolo Orlandoni y Jurgen Stock

Omar E. Carrero G. Ingeniero Forestal y M.Sc. en Economía, Universidad de los Andes (ULA), Venezuela. Profesor, ULA, Venezuela. Dirección: Escuela de Ingeniería Forestal, Facultad de Ciencias Forestales y Ambientales. Universidad de Los Andes, Mérida 5101, Venezuela. e-mail: neto@ula.ve.

Mauricio Jerez. Ingeniero Forestal y M.Sc. en Estadística Aplicada y Computación, ULA, Venezuela. Ph.D. en Manejo Ecológico Forestal, Louisiana State University, EEUU. Profesor, ULA, Venezuela. e-mail: jerez@ula.ve

Raúl Macchiavelli. Ingeniero Agrónomo, Universidad Nacional de Córdoba, Argentina. M.Sc. en Biometría, Universidad Nacional de Buenos Aires, Argentina. Ph.D. en Estadística, Pennsylvania State University, EEUU. Profesor, Universidad de Puerto Rico, EEUU. e-mail: macchia@uprm.edu

Giampaolo Orlandoni Merli. Economista, ULA, Venezuela. M.S., Iowa State University, EEUU. Profesor, ULA, Venezuela. e-mail: orlandon@ula.ve

Jurgen Stock. Ingeniero Forestal y M.Sc. en Manejo de Bosques, ULA, Venezuela. Gerente de Investigación, Smurfitkappa Cartón de Venezuela. Venezuela. e-mail: jurgen.stock@smurfitkappa,com.ve

RESUMEN

La clasificación de rodales mediante "curvas de índice de sitio" es útil en el establecimiento y manejo de plantaciones forestales ya que de la calidad del sitio depende en gran parte el crecimiento y rendimiento de las mismas. La mayor parte de los datos forestales son obtenidos longitudinalmente en el mismo sitio y las varianzas a menudo cambian en el tiempo, lo que hace inapropiada la aplicación del modelo general ordinario. Los modelos mixtos son buenas alternativas al permitir modelar estructuras covariantes. En este estudio se utilizaron modelos lineales y no lineales para analizar la relación de la altura mayor con la edad en plantaciones de Eucalyptus urophylla en Venezuela, y hallar la "curva guía" más apropiada para la clasificación de sitios. Los datos consistieron en mediciones repetidas obtenidas a partir de parcelas permanentes. Los modelos ajustados fueron los de Schumacher (lineal), Chapman-Richards, Weibull, Logistica, Gompertz y el Monomolecular (no-lineales). Se emplearon tres parametrizaciones diferentes: 1) original, asumiendo asíntotas aleatorias; 2) reparametrización para que el factor local sea la asíntota; y 3) reparametrización para que el factor local sea el parámetro asociado al dominio de la función. En las dos últimas la mayor altura a edad base tiene asociado un efecto aleatorio. La clasificación de sitio se realizó utilizando el concepto de ‘índice de sitio del sitio’, interpretado como la altura más probable a la edad base, entre todos los rodales hipotéticos que podrían crecer en el siotio. El mejor ajuste correspondió a un modelo no lineal mixto tipo Chapman Richards.

Fitting site index curves using mixed models for eucalyptus urophylla plantations in Venezuela.

SUMMARY

Classification of stands through "site index curves" is useful for establishing and managing forest plantations because growth and yield is highly dependent on site quality. Most forestry data are taken longitudinally at the same sites and the variances often change with time, which makes unsuitable the application of the ordinary general linear model. Mixed models are good alternatives since they permit modeling covariance structures. In this work, linear and nonlinear mixed models were used to analyze the relationship of top height with age for Eucalyptus urophylla plantations in Venezuela, and to find the most appropriate "guide curve" for site classification. Data consisted of repeated measurements on permanent plots. The models fitted were Schumacher (linear), Chapman-Richards, Weibull, Logistic, Gompertz, and Monomolecular (non linear). Three different parametrizations were used: 1) original, assuming random asymptotes; 2) reparametrized so that the local factor is the asymptote; and 3) reparametrized so that the local factor is the parameter associated with the function’s domain. In the last two cases the largest height at the base age has an associated random effect. Site classification was done using the concept of ‘site’s site index’, interpreted as the most probable height at base age in all the hypothetical stands that could grow in that site. The best fitting guide curve was obtained with a mixed-model Chapman-Richards equation.

Ajuste de curvas de ìndice de sitio mediante modelos mistos para plantações de eucalyptus urophylla na Venezuela.

RESUMO

A classificação do lotes através do ‘curvas de índice de sitio’ é util para o estabelecimento e manejo de plantações florestais, já que da qualidade do sitio depende o crescimento e rendimento do mesmo. A maior parte dos dados florestais é obtida longitudinalmente no mesmo local e as variações com freqüência mudam no tempo, o que faz inapropriada a aplicação do modelo geral ordinário. Os modelos mistos são boas alternativas ao permitir modelar estruturas covariantes. Neste estudo se utilizaram modelos lineares e não-lineares para ajustar a relação de altura maior com idade em plantações de Eucalyptus urophylla na Venezuela, e encontrar a ‘curva guia’ para a classificação do site. Os dados consistieram de medições repetidas em plantações permanentes. Os modelos ajustados foram os de Schumacher (linear), Chapman-Richards, Weibull, Logistic, Gompertz e o Monomolecular (não-lineares). Empregaram-se três parametrizações diferentes: 1) original, assumindo assíntotas aleatórias; 2) reparametrização para que o fator local seja a assíntota; e 3) reparametrização para que o fator local seja o parâmetro associado ao domínio da função. Nas duas últimas a maior altura na idade base tem associado um efeito aleatório. A classificação de sitio se realizou utilizando o conceito de ‘índice de sitio do local’, interpretado como a altura mais provável na idade base, entre todos os pequenos lotes hipotéticos que poderiam crescer no sitio. O melhor ajuste correspondeu a um modelo não linear misto tipo Chapman Richards.

Palabras clave/ Efectos Aleatorios/ Eucalipto/ Medidas Repetidas/ Modelos No Lineales/ Plantaciones Forestales/

Recibido: 02/09/2007.  Modificado: 05/03/2008.  Aceptado: 07/03/2008.

La clasificación de sitios de acuerdo a su capacidad productiva es de gran importancia para la planificación del establecimiento y manejo de plantaciones forestales ya sea a pequeña o gran escala. Las diferencias en productividad de los sitios determinan aspectos críticos del manejo de una plantación tales como los tratamientos silviculturales, el turno de cosecha, las técnicas de explotación, los productos a obtener y la rentabilidad de la inversión (Scolforo, 1992; Andenmatten y Letourneau, 1998; Gaillard et al., 2000).

Una de las metodologías más usadas para clasificar sitios en plantaciones forestales monoespecíficas, coetáneas, es el ‘índice de sitio’ (SI) definido, de manera general, como la ‘altura mayor’ alcanzada por un rodal a una edad predeterminada, llamada edad base (Carmean, 1975, Clutter et al., 1983). La altura mayor del rodal se define como la altura total promedio de una muestra representativa de los árboles más altos (dominantes y codominantes) del rodal, aunque existen otros criterios, tomando en cuenta su distribución espacial en el rodal. Se considera que esta altura es un buen indicador de la productividad potencial de una especie para un sitio en particular. Esta afirmación se basa en la hipótesis de Eichhorn (Fontes et al., 2003) la cual establece que la producción total de un rodal (volumen en pie + volumen extraído en cortas intermedias) óptimamente ocupado es una función de la altura del rodal. Para diferenciar calidades de sitio para una especie dada se debe definir una familia de curvas de SI a partir de la altura dominante observada para una variedad de sitios tan amplia como sea posible.

Tradicionalmente las curvas de SI se han ajustado usando modelos de regresión lineal y no lineal clásicos que asumen normalidad, igualdad de varianzas e independencia de los errores. Tales supuestos pueden ser violados cuando se analizan datos longitudinales. Los modelos mixtos permiten representar apropiadamente la estructura de varianzas-covarianzas asociadas con datos de medidas repetidas (Davidian y Giltinan, 1995; Gregoire et al., 1995; Verbeke y Molenberghs, 2000; Diggle et al., 2002; Littel et al., 2006). Tales modelos mixtos han sido poco utilizados en el área forestal, a pesar de que permiten modelar la estructura de varianzas-covarianzas y así realizar inferencias más precisas. El objetivo de este trabajo es comparar el ajuste de diferentes modelos estadísticos a datos de altura mayor en función de la edad, considerando efectos fijos y aleatorios (modelos mixtos) a fin de seleccionar los más apropiados para desarrollar curvas de índice de sitio para Eucalyptus urophylla en los Llanos Occidentales de Venezuela, basándose en datos provenientes de parcelas permanentes de crecimiento que son unidades de muestreo en las que se hacen mediciones repetidas para obtener información sobre la dinámica del desarrollo de una masa forestal en condiciones determinadas.

Marco Teórico

Calidad de sitio

En las ciencias forestales el término ‘sitio’ se utiliza para designar la influencia del ambiente sobre la producción de un bien o servicio del bosque, ya sea madera, forraje o frutos. La calidad de sitio es la capacidad intrínseca que tiene para la producción de un bien y se expresa en términos de la producción de dicho bien (Torres y Magaña, 2001). En plantaciones forestales cuyo propósito es la explotación maderera, la calidad de sitio se define como su potencial para la producción de madera de una especie o un tipo de bosque, considerando que mejores calidades tendrán mayor producción (Clutter et al., 1983). Existen varias técnicas para la clasificación de los sitios. Clutter et al. (1983) las dividen en métodos directos y métodos indirectos; para los primeros es necesario que la especie de interés exista o haya existido en esa localidad y se basa en las características propias de la masa forestal de dicha especie. El segundo grupo relaciona el sitio con características como clima, suelo, topografía, y factores bióticos (Mora y Meza, 2004). Dentro de la primera clasificación, se considera el método conocido como ‘Índice de Sitio’ (SI), donde la altura mayor promedio del rodal es usualmente una medida derivada a partir de la altura total de un número representativo de árboles por unidad de área. La definición más utilizada asume que la altura mayor es la altura promedio de los 100 árboles dominantes y codominantes bien distribuidos espacialmente por ha (Clutter et al., 1983). El método del SI ha sido el más utilizado para determinar la calidad de sitio en plantaciones forestales coetáneas y se basa en los siguientes supuestos: 1) la altura mayor, en un rodal monoespecífico y coetáneo es poco influenciada por la densidad de la plantación; 2) el crecimiento en altura mayor del rodal sigue un patrón predeterminado; y 3) la altura mayor se correlaciona bien con la producción volumétrica (Clutter et al.1983). La edad base (E_b) es un criterio subjetivo que puede depender de la rapidez de crecimiento de la especie y/o el turno de aprovechamiento. Usualmente se escogen edades próximas al turno (Clutter et al., 1983) o cuando se alcanza el máximo crecimiento medio en altura para la especie (Andenmatten y Letourneau, 1998). Este método implica el desarrollo de un conjunto de curvas que representan patrones de crecimiento en altura en función de la edad, cada una con un valor numérico o cuantitativo asignado, de acuerdo a la altura alcanzada a la E_b. Dichos patrones se pueden representar con modelos matemáticos flexibles (Grey, 1989). Para ser confiables como instrumentos de predicción, estos modelos deben poseer las siguientes propiedades: a) dar estimaciones no sesgadas del SI con la misma precisión, independientemente de la edad considerada; b) ser invariantes a la edad base y con altura cero a la edad cero; c) la altura a la E_b debe ser igual a SI; y d) las curvas correspondientes a cada calidad de sitio deberían tener asíntotas estimadas independientemente (Grey, 1989). Este último requisito, aunque se basa en que algunas especies pueden alcanzar diferentes alturas máximas absolutas dependiendo del sitio (p.ej. si los sitios difieren en altitud). Existen tres estrategias o métodos para modelar estas curvas: a) el de la curva guía, b) el de la diferencia algebraica y c) el de predicción de parámetros (Grey, 1989; Clutter et al., 1983).

Los modelos pueden desarrollarse a partir de mediciones de la altura mayor de árboles a diversas edades provenientes de redes de parcelas permanentes, parcelas temporales, o de análisis troncales (Torres y Magaña, 2001; Clutter et al., 1983). El ajuste se hace mediante regresión, donde la altura mayor o una transformación de la misma es la variable dependiente, mientras que la edad (y a veces otras variables) es la independiente. Se selecciona el modelo con mejor ajuste y se usa para generar una familia de curvas que define las categorías de calidad de sitio. Las curvas de la familia pueden ser proporcionales entre sí (anamórficas), poseer diferente forma para cada calidad de sitio pero sin intersectarse (polimórficas disjuntas) o intersectarse (polimórficas no disjuntas; Clutter et al., 1983, Grey, 1989; Avery y Burkhart, 1994). En cualquiera de las formas es posible reparametrizar los modelos (Scolforo, 1992, García, 2005). Adicionalmente, estos modelos pueden tratarse como determinísticos o como estocásticos. El primer caso corresponde a la definición clásica en la que el SI es entendido como la altura observada a la edad base, la cual es una propiedad particular del rodal y es llamada ‘índice de sitio del rodal’. En el segundo caso, se toma en cuenta la variabilidad debida a condiciones climáticas y/o errores de muestreo y medición, que pueden afectar la altura; interpretándose en este caso, el SI, como la altura esperada o más probable a la edad base, considerando todos los rodales hipotéticos que podrían crecer en ese sitio. Según este enfoque, se hablaría del ‘índice de sitio del sitio’ (García, 2005).

Los modelos de índice de sitio incluyen uno o más parámetros a estimar. Puede haber parámetros globales comunes a todos los rodales o parcelas, y parámetros locales específicos para cada rodal o parcela que dependen del sitio en particular (García, 2005). Dichos parámetros pueden estimarse con técnicas como la predicción de parámetros (Clutter et al., 1983), el método de las ecuaciones diferenciales (García, 2005) y el ajuste de modelos lineales y no lineales con efectos mixtos (Hall y Bailey, 2001; Gregoire et al., 1995). Este último es el enfoque adoptado en este trabajo.

Modelos mixtos

En Venezuela y otros países latinoamericanos, los estudios de índice de sitio para especies forestales se han basado en el ajuste de observaciones repetidas provenientes de parcelas permanentes de crecimiento o de análisis troncales, a partir de las cuales se han generado curvas mediante análisis de regresión utilizando como base modelos del tipo

y= ƒ(X,q,e)            (1)

donde y: vector de observaciones (variable dependiente), f: una función lineal o no lineal en los parámetros, X: matriz de diseño (variables independientes), q: vector de parámetros de efectos fijos, y e: vector de errores.

El uso de funciones no lineales tiene la ventaja de que pueden ser derivadas de consideraciones teóricas y sus parámetros pueden reflejar aspectos de interés para el modelador o usuario, pues permiten parametrizar mejor las distintas etapas del crecimiento (Schabenberger y Pierce, 2002). En ambos casos, se asume que e se distribuye normal e independientemente con media de 0 y matriz de varianzas constantes S (Draper y Smith, 1981; Kuehl, 2001; Balzarini et al., 2005). Las aplicaciones de estos modelos con supuestos simplistas para la componente de error, no tienen en cuenta que en el caso de medidas longitudinales (repetidas en tiempo sobre una misma unidad experimental), la estructura de la matriz de varianzas y covarianzas de los errores no se corresponde con tales asunciones. El hecho de realizar mediciones repetidas sobre el mismo sujeto implica que no es posible aleatorizar el factor tiempo, por lo que las medidas tomadas sobre un mismo individuo están autocorrelacionadas y, por tanto, no se cumple el supuesto de independencia de los errores (Schabenberger y Pierce, 2002). Asimismo, los modelos de crecimiento en altura se caracterizan por un incremento de las varianzas de los errores con la edad (heterocedasticidad). Estos problemas conllevan deficiencias en la precisión y capacidad de predicción de los modelos ajustados bajo los supuestos clásicos.

Los modelos lineales y no lineales mixtos surgen de incorporar efectos aleatorios, diferentes de los asociados con el término de error, por lo que podrían estar representados por

y= ƒ(X,q,b,e)    (2)

donde y, f, X, q y e son como en la Ec. 1, y b: vector de efectos aleatorios que modela la heterogeneidad entre sujetos.

A diferencia de los modelos lineales y no lineales clásicos, los modelos de efectos mixtos permiten modelar la estructura de los errores en datos provenientes de mediciones longitudinales. La modelización considera al efecto local (q) como una variable aleatoria que varía dentro de poblaciones de rodales (efectos aleatorios de rodal). Se asume que q tiene una cierta distribución de probabilidad en la población y que el valor observado del parámetro en la muestra proviene de esa distribución; sin embargo en ocasiones, cuando tiene una varianza grande, el valor del parámetro puede estar muy alejado (García, 2005). Varios problemas de estimación del crecimiento de árboles y rodales han sido tratados empleando modelos mixtos (Lappi y Bailey, 1988; Hall y Bailey 2001; Durán et al., 2002; Hall, 2004; Calegario et al., 2005; Gregoire et al., 1995). La mayor ventaja de los modelos mixtos es la generalidad en la inferencia, así como la posibilidad de modelar la correlación entre observaciones del mismo rodal (Calegario et al., 2005). La estimación de parámetros en estos modelos se hace por métodos de verosimilitud.

Criterios de selección de modelos

Los modelos mixtos se ajustan por máxima verosimilitud (MV) ó máxima verosimilitud restringida (MVR), por lo que la comparación de los modelos se puede realizar con la prueba de razón de máxima verosimilitud (PRMV) definida como

PRMV = -2log Verosimilitud (MR)

           -2log Verosimilitud (MC)

donde MR: modelo reducido y MC: modelo completo. El valor calculado de PMV se compara con una c²_{(kmc-kmr, α)}, donde kmc: número de parámetros del modelo completo, kmr: número de parámetros del modelo reducido, y α: nivel de significación. Si PMV >c²_{(kmc-kmr, α)} entonces el modelo correcto es el completo, en caso contrario debe seleccionarse el reducido (Schabenberger y Pierce, 2002; Balzarini et al., 2005; Calegario et al., 2005). Otros criterios utilizados son el criterio de información de Akaike [AIC= -2log (MV) + 2(q+p)], el criterio de información de Akaike corregido [AICC= -2log (MV) +2n(p+q+1)/(n-(p+q)-2] y el criterio de información bayesiano de Schwartz (BIC= -2log (MV) + 2(p+q) log (n), donde MV: máxima verosimilitud, q: número de parámetros de la estructura de covarianza, p: rango de la matriz de diseño X, y n: número de observaciones. En todos estos criterios, a menor valor, mejor será el ajuste (Calegario et al., 2005).

Materiales y Métodos

Datos

Los datos para ajustar las curvas de índice de sitio provienen de mediciones de parcelas permanentes localizadas en plantaciones de Eucalyptus urophylla (6800ha) de 4-12 años de edad, establecidas por Smurfit Cartón de Venezuela. Dichas parcelas abarcan una variedad de condiciones de sitio en los estados Portuguesa, Lara y Cojedes. La temperatura promedio anual varía entre 24 y 28,5°C, y la precipitación promedio anual es 1330-1715mm. Los suelos de la región fueron formados en el cuaternario por sedimentación de los diferentes ríos que la cruzan. Las plantaciones crecen en variados tipos de suelos: Alfisoles, Ultisoles, Inceptisoles, Entisoles, Molisoles y Vertisoles (Henri, 2001). Los mejores rendimientos se obtienen en zonas con suelos profundos, con texturas areno-francosa hasta franco-arcillo-limosa. En esas localidades es posible encontrar Dystropepts en combinación con Tropaquepts, con pendientes de 3-8% y Dystropepts con Haplustults, con textura pesada y pendiente de 8-20%. Los suelos presentan buen drenaje y pH moderadamente acido, variables estrechamente relacionadas con el relieve. Los peores rendimientos para E. urophilla se encuentran en localidades donde predominan Tropaqualfs, Haplustalfs y Ustropepts. Estos suelos por lo general tienen texturas pesadas, presentando en algunos casos pedregosidad, mal drenaje, pH ácido y se observa con frecuencia moteados grises (Acosta, 2004).

Todas las parcelas se ubicaron en plantaciones establecidas utilizando las mismas técnicas de preparación de sitio e implantación. El distanciamiento inicial fue de 3×3m (1111 árboles/ha) para todas las plantaciones, sin embargo hubo diferencias en la fuente de la semilla (Tabla I).

Las parcelas tienen una superficie de 495m² (55 árboles/parcela) y se midieron anualmente a partir del segundo año del establecimiento de la plantación; sin embargo, por lo general no se pudo medir en el mismo mes, por lo que el intervalo de medición no es exactamente de un año. En cada medición se evaluó todos los árboles vivos presentes. A cada individuo se le midió la circunferencia a la altura de pecho (mm), la altura total (m) y la presencia de bifurcación por debajo de 1,30m. Se obtuvieron 4-10 mediciones anuales repetidas, dependiendo de la edad de la plantación en la que se registraba la observación (2-12 años). Por tanto, el conjunto de datos para este estudio es de tipo longitudinal con observaciones irregularmente espaciadas en el tiempo y desbalanceado, ya que el número de observaciones difirió entre unidades experimentales (Figura 1).

Selección de modelos

Se seleccionó un conjunto de modelos lineales y no lineales para ajustar los pares de observaciones de altura mayor (m) y la edad (años) para todas las parcelas a fin de generar una ‘curva guía’ (Tabla II). Dichos modelos han sido utilizados para el desarrollo de curvas de índice de sitio y modelos de crecimiento en plantaciones forestales de diferentes especies (Torres y Magaña, 2001, Pagel et al., 2006). La ventaja de estos modelos es la posibilidad de darle una interpretación biológica a sus parámetros (Palahí et al., 2004; García, 2005; Pagel et al., 2006). Para todos los modelos, el parámetro ‘a’ representa la asíntota o valor máximo que puede alcanzar la curva; ‘b’ corresponde a un parámetro global a todos los rodales relacionado con la altura alcanzada a la edad base y ‘c’ es un parámetro local específico de cada rodal y que podría representar alguna condición particular del sitio (García, 2005). El modelo de Schumacher se ajustó por mínimos cuadrados ordinarios usando el procedimiento GLM del software SAS v.8.1 (SAS, 2001) como un modelo lineal de efectos fijos, siendo la variable dependiente el logaritmo natural de la altura dominante en cada parcela. Esta transformación es comúnmente usada ya que podría eliminar problemas de normalidad y heterogeneidad de varianzas en la distribución de los términos de error. Se realizó un análisis de residuos a fin de controlar el cumplimiento de los supuestos del modelo lineal en este conjunto de datos. Posteriormente, este mismo modelo se ajustó como un modelo lineal mixto, asumiendo diferentes estructuras para la matriz de varianzas-covarianzas de los términos de error, a fin de tomar en cuenta las desviaciones de los supuestos del modelo lineal general clásico, particularmente la correlación temporal. Se probaron ocho estructuras de la matriz de varianzas-covarianzas: 1) sin estructura (UN), 2) componentes de varianza (VC), 3) autoregresiva de primer orden (AR(1)), 4) simetría compuesta (CS), 5) Toeplitz (TOEP), 6) autoregresiva de primer orden heterocedástica (ARH(1)), 7) simetría compuesta heterocedástica (CSH) y 8) Toeplitz heterocedástica (TOEPH).

Los modelos no lineales considerando solo efectos fijos (Tabla II) se ajustaron mediante el procedimiento NLIN de SAS con el método de máxima verosimilitud, usando el algoritmo de Newton-Raphson. Los modelos no lineales de efectos mixtos se ajustaron por el método de máxima verosimilitud con el procedimiento NLMIXED de SAS. Se ajustaron tres grupos de modelos mixtos no lineales, dependiendo de si el modelo fue reparametrizado o no, y del parámetro al que fue asociado un efecto aleatorio (parámetro local). De tal manera se tiene: a) modelos no lineales (fórmula original) con asíntota aleatoria (a+u); b) modelos no lineales reparametrizados en función de la asíntota ‘a’ y asumiendo que el efecto aleatorio se asocia con la altura mayor a la edad base (S+u), y c) modelos no lineales reparametrizados en función de ‘b’ y suponiendo que el efecto aleatorio se asocia con la altura mayor a la edad base (S+u). La edad base (E_b) seleccionada fue de 6,5 años, que corresponde al turno de cosecha para esta especie. Los modelos reparametrizados dan origen a familias de curvas si se asignan valores fijos a S (p.ej. 15, 20, 25 y 30m). El primer grupo de modelos no lineales generará una familia de curvas anamórficas en la estructura de medias mientras que el segundo (reparametrizados), genera una familia de curvas polimórficas. Para ajustar modelos no lineales fue necesario establecer valores iniciales para los parámetros a estimar. Buenos valores iniciales fueron obtenidos de acuerdo a los procedimientos sugeridos por Schabenberger y Pierce (2002). Para identificar el mejor modelo se utilizaron los criterios estadísticos -2 Log Verosimilitud, AIC, AICC y BIC. Dado que la calidad del ajuste no refleja necesariamente la calidad de las predicciones futuras, los modelos deberían ser validados usando un conjunto de datos independientes (Kozak y Kozak, 2003; Quevedo et al., 2004) o bien ajustando los modelos a una serie más corta de observaciones (datos de calibración) y dejando la ultima porción de la serie observada como datos de validación. Debido a que la longitud de las series disponible es reducida, se optó por usar todos los datos disponibles para ajustar los modelos y esperar hasta tener nuevos datos a fin de llevar a cabo una validación más estricta. Se realizó una inspección del comportamiento de los residuos (valores observados - valores predichos) de cada modelo ajustado. Además de los criterios estadísticos, se tomó en cuenta si la curva generada cumplió con los requerimientos mencionados en la introducción para ser un buen modelo de índice de sitio. Con la finalidad de ilustrar la aplicabilidad de estos modelos para la clasificación de sitios, con el modelo de mejor ajuste se generaron familias de curvas anamórficas y polimórficas para delimitar las distintas calidades de sitio. (Avery y Burkhart, 1994; Torres y Magaña, 2001) sustituyendo en el modelo resultante los valores de S_b y E_b.

Resultados y Discusión

Los valores obtenidos con cada modelo para los criterios de bondad de ajuste se presentan en la Tabla III. Debe tenerse en cuenta que los resultados obtenidos para los modelos lineales no pueden compararse con los de los modelos no lineales, ya que las variables dependientes no son las mismas.

Entre los modelos lineales el peor ajuste se obtuvo para la variante de efectos fijos. El análisis de residuos para este modelo demostró que éstos siguen una distribución normal (prueba Kolmogorov-Smirnov (p>0,15)). Un gráfico de valores observados vs predichos no mostró un patrón evidente en la dispersión de los residuos, por lo que no se tuvo evidencia de heteroscedasticidad. Sin embargo, se observó una fuerte autocorrelación entre los residuos (Durbin Watson= 0,646; coeficiente de autocorrelación de primer orden de 0,672), lo que justificó el modelado de la matriz de varianzas y covarianzas mediante un modelo mixto de covarianza residual. El mejor ajuste de entre las variantes lineales lo presentó el modelo con matriz autoregresiva de primer orden (AR1) que obtuvo el menor valor para los criterios AIC, AICC y BIC. La desventaja de usar modelos linearizados es que los valores predichos (en este caso la transformación lnH) deben retransformarse para obtener el valor de las alturas que puede ser interpretado de manera directa. Según Arabatzis y Burkhart (1992) el modelo retransformado tiene propiedades distribucionales diferentes a los del modelo linealizado ajustado, por lo que usualmente debe usarse un factor de corrección.

En los modelos no lineales, los mixtos, todos los casos tuvieron un ajuste considerablemente mejor que los de efectos fijos (Tabla III). Entre los de efectos fijos, los mejores modelos fueron el monomolecular, Weibull y Chapman-Richards. Entre los mixtos el mejor ajuste fue para los modelos sin reparametrizar con asíntota aleatoria y los reparametrizados en función de ‘a’ siendo en ambos casos los ajustes muy similares para el monomolecular, Weibull y Chapman-Richards (Figura 2a). Para los modelos mixtos reparametrizados en función de ‘b’, el modelo de Chapman-Richards fue mejor que los restantes. En este caso, no se logró encontrar una solución para el modelo de Weibull, lo cual no fue sorpresivo, ya que para reparametrizaciones particulares de muchos modelos no lineales no es posible encontrar una solución (García, 2005).

Los modelos reparametrizados deberían presentar ventajas en la confiabilidad del ajuste, ya que este se hace considerando a S como parámetro cuya media y varianza son conocidas a priori ya que las alturas a la E_b han sido realmente medidas y pueden ser usadas como parámetros iniciales, mientras que los modelos sin reparametrizar utilizan una estimación del valor asintótico de las curvas (altura máxima que potencialmente pueden alcanzar los árboles) que no es conocido en la práctica. Adicionalmente al criterio de bondad de ajuste y la posibilidad de ajustar modelos reparametrizados, la función de Chapman-Richards cumple con las propiedades deseables para un modelo de índice de sitio, y presenta un comportamiento razonable que permite cierta flexibilidad para la extrapolación mas allá de las edades de medición, por lo que este modelo fue considerado como el mas apropiado entre los probados para determinar índices de sitio para la especie y localidades estudiadas. Las curvas ajustadas para las cuatro variantes del modelo Champan Richards sobrepuestas a los datos originales muestran que las curvas se comportan de manera muy similar hasta E_b, posteriormente el modelo reparametrizado en función de ‘b’ tiende a una asíntota más elevada que los restantes (Figura 2b). Al examinar los estadísticos para las diferentes variantes del modelo, se observa valores altamente significativos para los parámetros estimados, así como intervalos de confianza estrechos (Tabla IV). El análisis de la matriz de correlación entre los parámetros mostró que los tres primeros modelos presentaron una alta correlación entre los parámetros estimados, lo cual no es deseable ya que podría ser un indicativo de sobreparametrización. Sin embargo, el modelo reparametrizado en función de ‘b’ presentó valores muy bajos de correlación y por tanto no estaría sobreparametrizado (Tabla V). Las variantes ajustadas permiten generar familias de curvas anamórficas con asíntotas diferentes y polimórficas con asíntota común. Como ejemplo, se presenta familias de curvas generadas por el modelo mixto reparametrizado en función de ‘a’ (Figura 3a) y en función de ‘b’ (Figura 3b). Es posible reparametrizar el modelo de forma tal que las asíntotas puedan ser estimadas independientemente (Dieguez-Aranda et al., 2006). Sin embargo, los datos utilizados no dan información sobre la altura máxima que podría alcanzar la especie en cada calidad de sitio. Además de presentar mejores propiedades estadísticas que los modelos de efectos fijos, en los modelos mixtos se puede interpretar la altura observada para un rodal a una edad dada como un valor que tendría una cierta probabilidad de pertenecer a un sitio con características ambientales particulares, representado por una curva promedio de altura para ese sitio en particular (índice de sitio del sitio).

Conclusiones

Los modelos lineales mixtos probados tuvieron un mejor ajuste que el modelo lineal de efectos fijos. El modelo que mejor se comportó tuvo una estructura de residuos autoregresiva de primer orden heterocedástica (ARH1). Sin embargo los modelos lineales probados presentan dificultad en el momento de su interpretación, ya que la variable dependiente está transformada. Por ello, parece más razonable usar modelos no lineales, de los cuales el de Chapman Richards fue el más prometedor dadas sus buenas propiedades estadísticas, flexibilidad e interpretabilidad biológica. Cuando este modelo matemático fue tratado como un modelo mixto reparametrizado en función del índice de sitio se obtuvo el mejor ajuste entre todas las opciones utilizadas. Desde el punto de vista práctico este enfoque permite desarrollar modelos predictivos más precisos y confiables, particularmente si se desea determinar la calidad de sitio en la que se desarrollan plantaciones jóvenes a fin de predecir su rendimiento al momento de la cosecha.

REFERENCIAS

1. Acosta B (2004) Relación suelo-productividad en plantaciones de Eucalyptus urophylla y E: uroXgrandis en la División Forestal de Smurfit Cartón de Venezuela, Portuguesa. Universidad de Los Andes. Mérida, Venezuela. Mimeo. 98 pp.        [ Links ]

2. Arabatzis A, Burkhart H (1992) An evaluation of sampling methods and model forms for estimating height-diameter relationships in loblolly pine plantations. Forest Sci. 38: 192-198.        [ Links ]

3. Andenmatten E, Letourneau F (1998) Curvas de índice de sitio y crecimiento en altura para pino oregón (Pseudotsuga menziesi Mirb. Franco) de aplicación en la región andino patagónica de las provincias de Río Negro y Chubut, Argentina. Rev. Fac. Agron. La Plata. 103: 69-75.        [ Links ]

4. Balzarini M, Macchiavelli R, Casanoves F (2005) Aplicaciones de modelos mixtos en agricultura y forestería. Notas de Clase. Curso Internacional Aplicaciones de Modelos Mixtos en Agricultura y Foresteria. CATIE. Turrialba, Costa Rica. Mimeo. 189 pp.        [ Links ]

5. Avery T, Burkhart H (1994) Forest measurement. 4ª ed. Mcgraw-Hill. Nueva York, EEUU. 394 pp.        [ Links ]

6. Batista J, Couto H (1986) Escolha de modelos matemáticos para a construção de curvas de índice de sítio para florestas implantadas de Eucalyptus sp no estado de São Paulo. IPEF 32: 33-42.        [ Links ]

7. Calegario N, Maestri R, Leal C, Daniels R (2005) Estimativa do crescimento de povoamentos de Eucalyptus baseada na teoria dos modelos não lineares em multiníveis de efeito misto. Ciência Florestal 15: 285-292.        [ Links ]

8. Carmean W (1975) Forest site quality evaluation in the US. Adv. Agron. 27:209-269.        [ Links ]

9. Clutter J, Forston J, Pienaar L, Brister G, Bailey R (1983) Timber management: a quantitative approach. Wiley. Nueva York, EEUU. 333 pp.        [ Links ]

10. Davidian M, Giltinan D (1995) Nonlinear Models for Repeated Measurement Data. Chapman and Hall/CRC. Boca Raton, Florida. 360 pp.        [ Links ]

11. Diéguez U, Burkhart H, Amateis R (2006) Dynamic site model for Loblolly Pine (Pinus taeda L.) plantations in the United States. Forest Sci. 52: 262-272.        [ Links ]

12. Diggle P, Heagerty p, Liang K, Zeger S (2002) Analysis of Longitudinal Data. 2^a ed. Oxford University Press. Nueva York, EEUU. 396 pp.        [ Links ]

13. Draper N, Smith H (1981) Applied Regression Analysis. 2^a ed. Wiley. Nueva York, EEUU. 736 pp.        [ Links ]

14. Duran M, Donoso P, Quintero M (2002) Análisis del crecimiento en altura de árboles bajo un sistema de selección mediante cortas cíclicas usando modelos lineales mixtos. Rev. Forest. Ven 46: 5-62.        [ Links ]

15. Fontes L, Tome M, Coelho M, Wright H, Luis JS, Savill P (2003) Modelling dominant height growth of Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) in Portugal. Forestry 76: 509-523.        [ Links ]

16. Gaillard C, Pece M, Juárez M, Ríos N (2000) Estimación del crecimiento y producción de una plantación de paraíso gigante (Melia azedarach var. gigantea) sin riego en la provincia de Santiago del Estero, Argentina. Quebracho 9: 127-140.        [ Links ]

17. García O (2005) Site index: concepts and methods. University of British Columbia. web.unbc.ca/~garcia/publ/SiteSDEj2.pdf (Cons. 24/12/2005).        [ Links ]

18. Gregoire T, Schabenberger O, Barrett J (1995) Linear modelling of irregularly spaced, unbalanced, longitudinal data from permanent-plot measurements. Can. J. For. Res. 25: 137-156.        [ Links ]

19. Grey D (1989) Site index: A review. South Afr. For. J. 148: 28-32.        [ Links ]

20. Hall D, Bailey R (2001) Modeling and prediction of forest growth variables based on multilevel nonlinear mixed models. For. Sci. 47: 311-321.        [ Links ]

21. Hal D (2004) Multivariate multilevel nonlinear mixed effects models for timber yield predictions. Biometrics 60: 16-24.        [ Links ]

22. Henri C (2001) Soil site productivity of Gmelina arborea, Eucalyptus urophylla and Eucalyptus grandis forest plantations in western Venezuela. For. Ecol. Manag. 144: 255-264.        [ Links ]

23. Kozak A, Kozak R (2003) Does cross validation provide additional information in the evaluation of regression models?. Can. J. For. Res. 33: 976-987.        [ Links ]

24. Kuehl R (2001) Diseño de Experimentos: Principios Estadísticos de Diseño y Análisis de Investigación. 2^da ed. Thomson. México. 666 pp.        [ Links ]

25. Lappi J, Bailey R (1988) A height prediction model with random stand and tree parameters: an alternative to traditional site index methods. Forest .Sci. 34: 907-927.        [ Links ]

26. Littell R, Milliken G, Stroup W, Wolfinger R, Schabenberger O (2006) SAS for Mixed Models. 2^a ed. SAS Institute Inc. Cary, NC, EEUU. 795 pp.        [ Links ]

27. Mora F, Meza V (2004). Comparación del Crecimiento en Altura de la Teca (Tectona grandis) en Costa Rica con otros Trabajos Previos y con otras Regiones del Mundo. Seminario y Grupo de Discusión Virtual en Teca (Tectona grandis L.). Heredia, Costa Rica. [en línea]. <www.una.ac.cr/inis/teca/> (Cons. 24/12/2005).        [ Links ]

28. Ortega A, Montero G (1988) Evaluación de la calidad de las estaciones forestales. Revisión bibliográfica. Ecología 2: 155-184.        [ Links ]

29. Pagel E, Muller I, Guimaraes C, Schneider P, (2006) Ajuste e seleção de modelos tradicionais para série temporal de dados de altura de árvores. Ciência Florestal 16: 177-199.        [ Links ]

30. Palahí M, Tome M, Pukkala T, Trasobares AG (2004) Site index model for Pinus Sylvestris in North-East Spain. For.Ecol. Manag. 187: 35-47.        [ Links ]

31. Quevedo A, Jerez M, Moret A (2004) Selección de modelos para predecir distribuciones diamétricas en plantaciones de teca a partir de variables del rodal con la función Weibull. Rev. For. Ven. 48: 68-74.        [ Links ]

32. Schabenberger O, Pierce F (2002) Contemporary Statistical Models for the Plant and Soil Sciences. CRC Press. EEUU. 737 pp.        [ Links ]

33. Scolforo R (1992) Curvas de índice de sítio para Pinus caribaea var: Hondurensis. IPEF 45: 40-47.        [ Links ]

34. Smurfitkappa (s/f). Historia. Valencia, Venezuela. www.smurfitkappa.com.ve/ (Cons. 29/01/2007).        [ Links ]

35. Sullivan A, Clutter J (1972) A simultaneous growth and yield model for loblolly pine. For. Sci. 18: 76-86.        [ Links ]

36. Torres J (2001) Curvas de índice de sitio de forma y escala variables en la investigación forestal. Agrociencia 35: 87-98.        [ Links ]

37. Torres JM, Magaña O (2001) Evaluación de Plantaciones Forestales. Limusa. México. 472 pp.        [ Links ]

38. SAS (2001) SAS/STAT User’s Guide, v.8.1, SAS On-Line Documentation. SAS Institute Inc. Cary, NC, EEUU.        [ Links ]

39. Verbeke G, Molenberghs G (2000) Linear Mixed Models for Longitudinal Data. Springer. Nueva York , EEUU. 608 pp.        [ Links ]