INDIVIDUAL TREE DIAMETER GROWTH MODELS AND VARIABILITY OF MIXED NOTHOFAGUS SECOND GROWTH FORESTS IN SOUTHERN CHILE By PAULO CESAR MORENO MEYNARD A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PA RTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2017
2017 Paulo Cesar Moreno Meynard
To Elena my daughter Carmen and Neftal, my parents
4 ACKNOWLEDGME NTS I would like to thank the members of my supervisory committee, Drs. Gezan, Escobedo and Cropper, for sharing with me your time and knowledge through my program. Specially, I am very grateful with Dr. Gezan for giving me the opportunity to rea lize one o f my dreams, to study a Master in University of Florida. In addition, I thank him for his time, patience, guidance, energy, and dedication to help me achieve the goals of this research. I also want to thank the help of Dr. Ortega from Universidad Austral d e Chile, with the database on which this studio is based Thanks to Sebastian for sharing your knowledge in programming and to be an excellent coworker. Special thanks to my family, specially to my daughter Elena, who was without he r father alo ng this jo urney; my parents Carmen y Neftal and my sisters for helping me to realize this dream. Also, I would to thank for their support to Andrs and the Chilean community in Gainesville and so many friends, w ho became my family during these last two years
5 TA BLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 LIST OF ABBREVIATIONS ................................ ................................ ........................... 10 ABSTRACT ................................ ................................ ................................ ................... 13 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 15 2 UNDERSTANDING VARIABILITY IN MIXED NOTHOFAGUS SECOND GROW TH FORESTS IN SOUTHERN CHILE: A MULTIVARIATE ANALYSIS ....... 18 Background ................................ ................................ ................................ ............. 18 Materials and Methods ................................ ................................ ............................ 22 Data Description ................................ ................................ ............................... 22 Multivariate Methods ................................ ................................ ........................ 24 Results ................................ ................................ ................................ .................... 26 Discussion ................................ ................................ ................................ .............. 30 3 INDIVIDUAL TREE DIAMETER GROWTH MODELS OF MIXED NOTHOFAGUS SECOND GROWTH FORESTS IN SOUTHERN CHILE .............. 41 Background ................................ ................................ ................................ ............. 41 Materials and Methods ................................ ................................ ............................ 46 Data Description ................................ ................................ ............................... 46 Model Fitting ................................ ................................ ................................ ..... 48 Results ................................ ................................ ................................ .................... 52 Model Fitting ................................ ................................ ................................ ..... 52 Model Projection ................................ ................................ ............................... 56 Selected Model ................................ ................................ ................................ 56 Discussion ................................ ................................ ................................ .............. 57 4 CONCLUSIONS ................................ ................................ ................................ ..... 72 APPENDIX A DOMINANT HEIGHT SITE MODEL ................................ ................................ ....... 75 B SUPPLEMENTARY TABLES ................................ ................................ ................. 76
6 C SUPPLEMENTARY FIGURES ................................ ................................ ............... 77 LIST OF REFERENCES ................................ ................................ ............................... 82 BIOGRAPHICAL SKETCH ................................ ................................ ............................ 88
7 LIST OF TABLES Table page 2 1 Summary statistics of this study dataset for individual and stand variables in each growth zone. ................................ ................................ .............................. 33 2 2 Factor loading for the first five principal componen t (PC) obtained from the tree and stand subsets together. Values in parenthesis correspond to the proportion of the variance explained by each component. ................................ 34 2 3 Factor loading for the first fiv e principal component (PC) obtained from average plot AIDBH and stand subset. Values in parenthesis correspond to the proportion of the variance explained by each component. ........................... 35 3 1 Summary stat istics of individual and stand attributes for fitting database based on PP and TP networks. ................................ ................................ .......... 64 3 2 Parameter estimates of the four selected models for estimation of the logarithm of AIDBH (mm). ................................ ................................ ................... 65 3 3 Goodness of fit statistics of four models for ln(AIDBH) using test data. .............. 65 3 4 Projection goodness of fit statistics for DBH (cm) in a period of 6 and 12 years using projection database. ................................ ................................ ........ 66 3 5 Least square means (standard errors in parentheses) statistics of all pairwise Tukey comparisons between all 11 groups b ased for Specie Zone in CV + SpZone model. ................................ ................................ ................................ ... 67 A 1 Parameter estimates of the dominant height site model by specie and growth zone. ................................ ................................ ................................ ................... 75 B 1 P values with Bonferroni correction obtained from a post hoc test of all six possible pairwise comparisons for growth zones related to two clusters. ........... 76 B 2 Mean diameter breast height g rowth rates (standard deviation in parentheses) by specie and growth zone obtained from the available data (fitting database). ................................ ................................ ................................ 76 B 3 Goodness of fit statistics for factors specie, zone and bot h in the CV regression model using the test data. ................................ ................................ 76
8 LIST OF FIGURES Figure page 2 1 Non Metric Multidimensional Analysis (NMDS) using the tree sub set, identifying growth zones (a) and species (b). ................................ ..................... 36 2 2 Non Metric Multidimensional Analysis (NMDS) using the stand subset, identifying growth zones (a) and dominant species (b). ................................ ...... 37 2 3 Principal coordinates analysis (PCoA) using the tree and stand subsets together, identifying growth zones (a) and species (b). ................................ ...... 38 2 4 F irst two PCA axes with scores and factor loadings obtained using, the tree and stand subsets (a), stand subset together with average plot AIDBH. The two K means clusters are also identified (b). ................................ ...................... 39 2 5 Geographical localization of plots identified according to their two K means clusters in relation to growth zones. ................................ ................................ ... 40 3 1 Plots distribution in the study area (southern Chile) for the permanent plots (PP) and temporary plots (TP) networks. ................................ ........................... 68 3 2 Statistics C p BIC, RSS, R2adj and MSE K for different number of predictors with CV regression selection model using fitting databa se. ................................ 69 3 3 Fit and residuals plots of four selected models using test dataset. CV (a and b), LASSO (c and d), CV + SpZone (e and f), and LASSO + SpZone (g and h). ................................ ................................ ................................ ....................... 70 3 4 RMSE% (a) and BIAS% (b) by DBH class for the four selected models evaluated using the test dataset. ................................ ................................ ........ 71 3 5 Simulation prediction in future DBH (cm) for p rojection periods of 6 years (a) and 12 years (b) based on the CV + SpZone model using projection database. ................................ ................................ ................................ ............ 71 C 1 Scatterplot correlation matrix of predictors (in original units) for the fitti ng database. ................................ ................................ ................................ ............ 77 C 2 Scree plot of PCoA with the cumulative variance of each coordinate based on the continuous variables of the tree and stand subsets. ................................ 78 C 3 Scree plot of PCA for continuous variables of tree and stand subset. ................ 78 C 4 PC1 and PC3 with scores and factor loadings, using tree and stand subset. ..... 79 C 5 PC2 and PC3 with scores and factor loadings, using tree and stand subset. ..... 79
9 C 6 Scree plot of PCA for continuous variables of stand subs et and average plot of AIDBH. ................................ ................................ ................................ ........... 80 C 7 Scree plot of optimal number of cluster using within sum of squares and average silhouette width for stand level information. ................................ .......... 80 C 8 Histograms of expected deltas and observed delta obtaining by MRPP test for two K means cluster using stand level information. ................................ ...... 81
10 LIST OF ABBREVIATIONS A Breast Height Age As Chanc e Corrected Within group Agreement Statistic Ad Dominant Stand Age at Breast Height AIDBH Annual Increment in Diameter at Breast Height BA Stand Basal Area BAL Basal A rea of L arger T rees BALn Basal A rea of L arger T rees of Nothofagus BALr Relative Bas al A rea of L arger trees BAN Basal Area of Nothofagus BIAS% Relative B ias BIC CA Correspondence Analysis CONAF National Forest C orporation of Chile Cp CV C ross V alidation DBH Diameter at Breast H eight DS Dominant Specie s GVIF Generalized Variance Inflaction Factor G&Y Growth and Yield H Total Height Hd Dominant Height INFOR Forest Research Institute of Chile LASSO Least Absolute Shrinkage and Selection Operator ln Natural logarithm
11 MRPP Multiple R esponse Permutation Procedure MSE Mean Square Error MSE K Mean Squared Error averaged over the K folds N Trees per hectare NMDS N on metric Multidimensional S caling PC P rincipal Component PCA P rincipal Component A nalysis PCoA Principal C oordinates A na lysis PP Permanent Plots QD Quadratic mean Diameter R2adj Adjusted r squared R 2 emp Empirical coefficient of correlation RR Ridge Regression RS Relative Spacing Index RSS Residual Sum of S quares RMSE Root Mean Square Error RMSE% Relative Root Mean Square Error RORACO Chilean forest type composed by three Nothofagus species: roble ( Nothofagus obliqua ), raul ( N. alpina ), and c oihue ( N. dombeyi ) SD Standard Deviation SDI Stand Density Index SI Site Index SP Specie SS Sociological Status TP Tem porary Plots
12 VIF Variance Inflaction Factor
13 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science INDIVIDUAL TREE DIAMET ER GROWTH MODELS AND VARIABILITY OF MIXED NOTHOFAGUS SECOND GROWTH FORESTS IN SOUTHERN CHILE By Paulo Cesar Moreno Meynard August 2017 Chair: Salvador A. Gezan Major: Forest Resources and Conservation The second growth forests of Nothofagus obliqua (ro ble) N. alpina (raul) and N. dombeyi (coihue) defined locally as RORACO forest type, is one of the most important natural mixed forests of Chile. S eve ral studies have identified a wide range of factors that influence the phenotypic variability found in th is forests. Currently there is no individu al tree diameter growth model for the full geographical distribution of the RORACO forest type. Th is study, contributed to better understand ing the underlying variability for these forest s through the use of the following multivariate analysis : 1) non metric multidimensional scaling (NMDS); 2) principal coordinates analysis (PCoA) ; and 3) principal component analysis (PCA). In addition, a n individual tree diameter growth model was performed, comparing two selecti on variable procedures : 1) cross validation (CV) ; and 2) Least Absolute Shrinkage and Selection Operator ( LASSO ) regression The data in this study is based on tree and stand variables typically found in available forest inventories Findings indicate that site productivity and growth zones do not explain much of the variability of this sampled population. However, stand development stages, tre e to tree competition, and tree size attributes are critical variables with a high percentage of variance explained ranging from 61% to 67%. In addition, the final
14 diameter growth model used a CV procedure incorporating species and growth zones as a combined factor. A small set of predictors successfully explain ed a large portion of diameter breast height ( DBH ) growth with variables associate d to competition, tree and stand level The goodness of fit statistics for this final model were R 2 emp = 0.56, RMSE% = 44.49% and BIAS% = 1.96% for annual DBH growth prediction, and R 2 emp = 0.98 and 0.97 for DBH projection after 6 and 12 years, respectively. This model constitutes a simple and useful tool to support management and decision making for th ese ecosystems.
15 CHAPTER 1 INTRODUCTION One of the most important forest types in Chile is known locally as RORACO, its name o riginates from their three main Nothofagus species: roble ( Nothofagus obliqua (Mirb.) Oerst ), raul ( N. alpina (Poepp. Endl.) Oerst. ), and c oihue ( N. dombeyi (Mirb.) Oerst. ) These forests are found in t he south of Chile from 36S to 41S latitude approx imately from the co ast to the Andes range, corresponding to 10.8 % of the total native forest area of the country (CONAF 2011 ). In 2015, the saw timber volume from these three Nothofagus species was approximately 45% of the national production (INFOR 2016). The volume growth rates for these forests range between 6 and 10 m 3 ha 1 year 1 (Donoso et al. 2004) corresponding to some of the highest values for southamerican temperate native forests Despite the importance of this resource forest managers and owne rs do not have enough information and tools regarding present conditions and future state under varying silvicultural activities Some information is available locally, but not for the full geographical range of the RORACO forest type In addition, anothe r related problem is the h igh phenotypic variability found in these forests due to different factors, such as large geographical area, varying (or lack of) management and strong anthropic disturbances, together with diverse biophysic al variables ecology and species (Donoso 1993; Chauchard and Sbrancia 2003; Echeverria and Lara 2004; Luebert and Pliscoff 2006). A very important tool to describe and project any stand structure is i ndividual tree growth and yield (G&Y) model wit h the possibility of incorpo rating management and estimate s of timber products (Andreassen and Tomter 2003 ) These model s are
16 constituted by three components (Burkhart and Tome 2012) : 1) individual tree diameter at breast height (DBH) or basal area growth equation ; 2) individual tree height growth equation or, alternatively, a height diameter relationship ; and 3) mortality or survival component Most G&Y models, particularly individual tree growth models, are obtained by fitting multiple linear regressions with many potential predicto rs in their original units or using transformations. These predictors often suffer from multicollinearity due to high correlation be tween them (Welham et al. 2015), but t here are several statistical m ethods for selecting predictors to avoid multicollineari ty such as stepwise or shrinkage methods (James et al. 2013). The data for this research originated from a study from the Universidad Austral de Chile (Orte ga and Gezan 1998) consisting of two plot networks. One with a total of 128 permanent plots, and an other with 60 temporary plots The data were collected using a stratified population, based on a national forest cadastre (CONAF et al. 1999), for the full geographical distribution of the RORACO forest type (> 14 ,000 k m 2 including four growth zones prop osed by Gezan and Moreno ; 1999 ). The variables used in this study are based on tree and stand attributes typically obtained from forest inventories, such as site productivity, tree to tree competition, individual tree and whole stand attributes (Vanclay 19 94) The refore, the aims of this study were: 1) to understand the underlying phenotypic variability for second growth forest s of N. alpina N. obliqua and N dombeyi through unsupervised multivariate methods, and 2) to develop an individual
17 t ree growth model to estimate annual increment in DBH for the same species, usi ng a set of tree and stand level predictors. I will analyze the first aim using s everal m ultivariate analysi s of ordination and clustering (Chapter 2) I will specific ally analyze this obje ctives: 1) to reduce the dimensionality of the data through ordination techniques, such as non metric multidimensional scaling (NMDS), principal coordinates analysis (PCoA) and principal components analysis (PCA ); 2) to identify those variables or factors that explain high amount of data variability; and 3) to generate K means clusters to assess the validity of the growth zones. For the develop ment of an individual t ree growth model (Chapter 3), t wo selection variable procedures were evaluated: 1) model se le ction with cross validation and 2) LASSO regression. The resulting selected models were later extended using growth zone s or specie s as factors, to evaluate if they contribut e to improve the performance of these models, and therefore, better explain dyn amics for this forest type The specific objectives in this chapter were: 1) to compare different multiple linear fitting procedures to obtain the best individual tree growth model; 2) to generalize diameter growth models to incorporate growth zone and spe cie, and 3) to evaluate the final diameter growth model using an independent validation dataset.
18 CHAPTER 2 UNDERSTANDING VARIABILITY IN MIXED NOTHOFAGUS SECOND GROWTH FORESTS IN SOUTHERN CHILE : A MULTIVARIATE ANALY SIS Background The structure an d dynamics of a mixed natural forest are difficult to understand where only certain aspects of their variability can be studied in detail (Landres et al. 1999) Forest variability is often expressed as heterogeneity, defined as any form of environmental v ariation, physical or biotic, occurring in space and/or time (Kuuluvainen et al. 2002). Micro site or local factors influences individual trees, which are the main source of small scale environmental heterogeneity in forest ecosystems. These factors includ e: light quantity and quality due to changes in canopy, interception of precipitation from foliage and branch litter, and root uptake of water and nutrients, among others (Kuuluvainen et al. 2002). Hence, this small scale heterogeneity is often translated into tree to tree competition at many levels. This competition generates a high variability, where complex mixed and natural forests are often more heterogeneous. One of the most relevant forests in the south ern hemisphere are composed by species of the N othofagaceae family, with Nothofagus as their unique genre, which are found in Australia, New Zealand, New Guinea, New Caledonia, Chile, and Argentina. The southamerican Nothofagus species dominate temperate and sub Antarctic forests of Chile and Argentina In Chile, there are nine species including three evergr eens, and six deciduous (Ormazabal and Benoit 1987; Donoso 1993 ). A very important Nothofagus forest types in Chile is called RORACO, its name originates from their three main Nothofagus species: ro ble ( Nothofagus obliqua (Mirb.) Oerst ), raul ( N. alpina (Poepp. Endl.) Oerst. ), and c oihue ( N. dombeyi (Mirb.) Oerst. ) RORACO is the fourth largest forest type in Chile with 10.8% of the total surface
19 (CONAF 2011), with saw timber production reach ing 4 5% of the total wood market (INFOR 2016). The second growth forest stands of RORACO are found in t he south of Chile from 36S to 41S latitude approximately from the co ast to the Andes mountain range. This wide distribution leads to a diverse set of ecosy stems and growing conditions (Donoso 1993 ; Donoso et al. 1993). These forests are commonly m ixed stands, and are characterized by the presence of two or more vertical strata where Nothofagus species are the dominant over story s trata and companion are ofte n shade tolerant (Donoso 1981, 1993). In addition, it is common to find some level of hybridization between Nothofagus species (Donoso 1987). All of the above elements contribute to a high natural variability of this forest type. Several studies have chara cterized the ir different aspects, such as: tree volume, taper equations, growth zones, diameter breast height (DBH) growth models, height diameter relationships mortality, etc. (Cubillos 1987; Donoso et al. 1993 ; Schlatter and Gerding 1995; Ortega and Gez an 1998; Moreno 2001; Chauchard and Sbrancia 2003; Lusk and Ortega 2003; Echeverria and Lara 2004; Salas and Garcia 2006; Donoso and Lusk 2007; Gezan et al. 2007, 2009; Thiers et al. 2008; Esse et al. 2013). RORACO forest type is a complex ecosystem with many origins, where its classification presents diffuse boundaries overlapping with others forests type, and they are present in a continuum vegetation gradients spanning wide geographical ranges in latitude and altitude (Donoso 2003). These forests belong to the bio climate temperate hiperoceanic with a large extension from coastal zones, interior valleys, and m ountainous areas, containing several microclimates (Luebert and Pliscoff 2006).
20 Echeverra and Lara (2004) determined that the variability in diame ter growth in N. alpina and N. obliqua were associated mainly to climate and soils. Another study found that, at a macro region scale, N. obliqua productivity depends mainly on climate, but at a regional scale, edaphic and topographic variables are more re levant (Thiers et al. 2008). Esse et al. (2013) found that the natural spatial distribution of N. dombeyi is explained, in decreasing order, by annual average rainfall, thermal oscillation and soil drainage. Two critical factors that affect the variabilit y of RORACO forests are the presence of pure or mixed forest, and their regeneration strategy (Donoso 1993). Presence of pure or mixed forest is explained with the interaction between autecology of each specie and biophysical variables, where, for example, extreme climate and soils conditions tends to favor pure forest (Donoso 1993). In addition, high and low forests are the result of different regeneration strategies, including sexual (seeds) and asexual propagation (coppice). Forest disturbances are an i mportant factor in the ecology of these Nothofagus species, contributing to an increase in their variability. Large scale disturbances usually promote a vast regeneration event, which promotes even aged strata or cohorts (Luebert and Pliscoff 2006). Howeve r, Donoso (1993) points out differences in the ecology strategies of the Nothofagus species, which range from big disturbances to patch dynamics. Here, N. dombeyi and N. obliqua are more aggressive than N. alpina preferring large disturbances; whereas the latter, due to its shadow tolerance nature, favors patch dynamics. RORACO forest suffer from a high level of anthropic disturbances affecting severely the landscape and their ecology (Luebert and Pliscoff
21 2006). For example, these landscape transformatio ns can influence the development of varieties or even geographic races (Donoso 1993). Individual and stand level g rowth has been always an important variable used to characterize most forest. Donoso et al. (1993) defined four growth zones for N. alpina a nd N. obliqua where the better growth zones corresponded to those with temperate climate. Esse et al. (2013) proposed for N. dombeyi five growth zones, which were define d in decreasing order, by climate, soil, and topography. Also, Gezan and Moreno (1999 ), based mainly on climate, soil, site productivity defined four growth zones and 11 sub zones. Several studies on volume equations, basal area growth, dominant height site, and taper equations have found this zoning useful to model fitting ( Moreno 2001; Gezan et al. 2009). Multivariate analysis is a group of methods that help to describe and/or explore several data and they provide with formal inferences about the population (Everitt and Hothorn 2011). There are two methods of multivariate analysis: supe rvised and unsupervised (James et al. 2013) In supervised methods the goal is predict a desired response from a set of predictor variables, whereas for unsupervised methods the goal is discover relationships among the data from the complete set of variabl es. For example, with the unsupervised methods of ordination or dimensionality reduction it is possible to understand data variability, or by using cluster analysis unknown similarities between objects can be discovered (James et al. 2013). Some multivaria te methods of ordination are: principal component analysis (PCA), principal coordinates analysis (PCoA, also called multidimensional scaling), non metric multidimensional scaling (NMDS), and correspondence analysis (CA) (Everitt and Hothorn 2011). These
22 me thods help reduce the complexity of the data, visualize multidimensional data and develop research hypotheses (Everitt and Hothorn 2011). Cluster analysis allow finding similar subgroups of objects in a dataset, where some of the approaches commonly used a r e: agglomerative hierarchical, K means cluster ing, and model based clustering (Everitt and Hothorn 2011). The main objective of this stu dy is to understand the underlying variability of second growth forest s of N. alpina N. obliqua and N dombeyi thr ough multivariate analysis of ordination (non metric multidimensional scaling, principal coordinates analysis and principal compon ents analysis) and clustering (K means) using a group of tree and stand variables typically obtained from forest inventories. In particular 158 plots (1,1 08 trees with growth data) cover ing the complete geographical distribution of the RORACO forest type. The specific objectives are: 1) to reduce the dimensionality of the data through ordination techniques; 2) to identify those v ariables or factors that explain high amount of data variability; and 3) to generate clusters to assess the validity of the growth zones proposed by Gezan and Moreno (1999). Materials and M ethods Data D escription A total of 1,108 trees of N. alpina N. obl iqua and N dombeyi with diameter growth information were available for this analysis. These originated from a sample of 158 plots from a study of the Universidad Austral de Chile (Ortega and Gezan 1998) consisting on two plot networks. The data were coll ected under a population stratification based on a national forest cadastre (CONAF et al. 1999). Further details of this dataset can be found in Gezan and Moreno (1999 ).
23 This study used the average annual increment in DBH based on the penultimate and antep enultimate growth year s, which was obtained from increment cores ( 864 ) and stem sections (244) at breast height. The last growth year was discarded as p lots were established at different times within the growing season T he dataset includes several tree a nd stand level variables (Table 2 1 ). The tree level attributes were: species, diameter breast height (DBH, cm), total height (H, m), breast height age (A, year), basal area of larger trees (BAL, m 2 ha 1 ), basal area of Nothofa gus of larger trees (BALn m 2 ha 1 ) sociological status (S S defined according to vertical stratification with dominant (1) codominant (2) intermediate (3) or suppressed (4)), relative basal area of larger trees (BALr, with BALr = BAL /BA) and the averag e a nnual increment in DBH (AIDBH, m m year 1 ). The stand level variables considered were: growth zone (Zone), stand basal area (BA, m 2 ha 1 ), density or number of trees per hectare (N, trees ha 1 ), quadratic mean diameter (QD, cm), dominant height (Hd, m), dominant breast height stand age (Ad, years), dominant specie (DS), site index (SI, m), basal area of Nothofagus (BAN, m 2 ha 1 ), stand density index (SDI, trees ha 1 ) and relative spacing index (RS). BAN was calculated as the sum of basal area of N. alpina N. obliqua and N. dombeyi Hd and Ad were calculated using the 100 largest diameter trees per hectare, and DS corresponded to the Nothofagus specie with the largest proporti on of BAN. SI was estimated using available dominant height site models obtained on the same dataset for an index age of 20 years (see appendix A). (Avery and 1.4112 as reported by Gezan et al. (2007) fitted using both the permanent plot ( PP ) and temporary plot ( T P ) networks. Finally, RS was calculated according to: RS = [(10000/N) 0. 5 ]/Hd (Avery and Burkhart 2002; Contreras et al. 2011).
24 Multivariate M ethods To understand the underlying variability of the RORACO forest three methods of ordination were evaluated wi th the available data: 1) non metric multidimensional scaling (NMDS), 2) principal coordinates analysis (PCoA), and 3) principal component analysis (PCA). The goal of ordination methods is to represent the main trends of the data through orthogonal axes (l inearly independent, uncorrelated), reducing the original number of dimensions. The observed trends can be interpreted graphically or be further evaluated with the use of other multivariate methods such as clustering (Borcard et al. 2011). NMDS is probably the most flexible ordination method allowing the use of continuous and categorical variables (Everitt and Hothorn 2011) and unlike PCoA and PCA, it is not an eigenvalue eigenvector based method. NMDS, for the new data representation, keeps the order of t he relationships among the initial set of objects but no t the exact distances (Borcard et al. 2011). This method tolerates missing values, as long as there is enough information to identify the position of each object with respect to others (Borcard et al. 2011). For this study, t he reduction of dimensionality was established using two axes (k = 2) with a modified Gower dissimilarity matrix (Anderson et al. 2006). Two NMDS analyses w ere performed one with the tree subset and another with the stand subset. Based on the analysis, ordination plots were used to represent the relationships among objects (trees or plots), where additional factors such as, growth zones, tree species, and dominant species were used in the identification of objects to help with the visualization of results. Here, object clouds were delimited using the extreme objects positions in relation to each factor level.
25 PCoA, in contrast to NMDS, that orders the data in axes according to their contribution to total variability based on eigenva lues. PCoA does allow the use of any similarity/dissimilarity matrix which represents the relationship between objects or variables For PCoA in this study, both the tree and stand subsets were analyzed together, where the number of coordinates were selec ted by examining a scree plot. Here, the same modified Gower dissimilarity matrix and g raphical output generated for NMDS were also used PCA is the main eigen value eigen vector based method that describe the variation of correlated quantitative variables, in terms of a new set of uncorrelated variables or ax e s As with PCoA, the se new axis or principal components (PC) are derived in decreasing order of importance according to their total variance explained (Everitt et al. 2011). PCA uses a similarity/dissim ilarity matrix calculated based on the variances and covariances of the variables, or the correlations between them (Borcard et al. 2011); thus PCA uses the Euclidean distance to represent the objects. Representation of the objects (scores) and the origina l variables (factor loadings) can be done through a PCA Biplot. In this study, an exploratory PCA was used to identify relevant continuous variables. Here, two analyses were performed, one used both the tree and stand subsets, and a second only stand varia bles with the inclusion of AIDBH plot average. For this method, in order to achieve normality assumptions, some variables were transformed using the inverse, logarithm, square root, or arc sin functions. After transformations, the variables were standardiz ed. The number of PCA axes was determined by using (Quinn and Keough 2002): 1) a latent root criterion by keeping components with eigenvalues >1; 2) through a scree plot; and 3) by having a
26 percent variance explained greater than 90%. In order to interpret th e relationship among objects, distance Biplot were used. K means clustering is a method for dividing a dataset into K different and non overlapping clusters. In this study, t his method was implemented only with the subset of stand variables standardized The determination of the number of clusters was done through scree plots of: 1) within cluster sum of squares; and 2) silhouette method (Everitt and Hothorn 2011). In addition, a Multiple Response Permutation Procedure (MRPP) was performed to test the di fferences between clusters (McCune et al. 2002). Pearson chi squared test of association was used to establish any association between clusters and growth zones as defined by Gezan and Moreno (1999). If any relationship was found, a post hoc test for all p airwise comparisons with Bonferroni corrections of the P values was performed ( MacDonald and Gardner 2000) All statistical analyses were performed using R v. 3.2.2 (R Core Team 2015). For all analyses, where required, a significance level of 5% was used. The R packages used for dissimilarities matrices NMDS PCoA, and MRPP was vegan version 2.4 1 (Oksanen et al. 2013) and for Silhouette plot cluster version 2.0 4 (Maechler 2013) Results The variables of tree and stand subset presented, as expected, a h igh level of correlation between some variables, such as BA BAN, DH RS, BA SDI, BAL BALn, DH QD, A DA, RS QD, and DBH H, with values >0.83 (see Figure C 1 ). NMDS, using categorical and continuous data for the tree subset, con verged with a stress goodness of fit value of 0.09, corresponding to a good representation of original distance matrix For this subset, NMDS ordination with the two selected axes did not identify different groups, only a single group surrounded by isolate d data points
27 (Figure 2 1 ) For this representation, there are no clear patterns for the different growth zones or species (growth zones represent almost exactly the same tree cloud, and N obliqua had minor differences relate d to N. alpina and N. dombeyi ) indicating that these factors do not explain the variability of the two selected axes. Another NMDS was performed only with the stand subset (Figure 2 2 ) resulting in a stress value o f 0.05, w hich represents a good description of the original distance matrix (Everitt and Hothorn 2011) As with the tree subset, the clouds for growth zones were overlapping; however, the dominant specie had a small difference between N. obliqua and N dombeyi whi le N alpina is practically contained within the other two dominant species clouds. The ordination method of PCoA explained ~90% of the cumulative variance using two dimensions (PCoA axes) when both tree and stand variables were evaluated together (Figure C 2 ). Interestingly, trees belonging to a plot were very close in these new coordinates (Figure 2 3 ) ; hence, it appears that the variability associated with the stand variables dominates over the one with individual tree variables. In addition, the results from PCoA confirms the findings from NMDS, with overlapping growth zones, and greater similarities between N. alpina and N. dombeyi The PCA method was also used, based on both subsets, to assist in identification of those variables that better explain the variability of the sample population. T o achieve normality the variables D BH, BA, BALn, and AIDBH were transformed with a square root and BALr was transformed with the arcsine function. The re sults of PCA indicated that four or five principal components (PC) explain the majority of the variability (88% and 92%, respectively) found in these Nothofagus forest s
28 (Figure C 3 ). Factor loadings for these five PC are prese nted in Table 2 2 PC1 with to DBH PC1 is explained mostly by the stand variables BA, QD, Hd, Ad, BAN, SDI, and RS. The other stand variables N and SI were relevant only on the third and fourth component, respectively. In relation to tree variables, it was observed a high contribution in PC1 for variables related to tree size (DBH, H, and A), while tree to tree competition variables (BAL, BALn, and BALr) were important in PC2 with almost identical factor loadings, and AIDBH was also st rongly associated with this component but with opposite direction (i.e., negative factor loading). A Biplot for PC1 and PC2 is presented in Figure 2 4 a that explains 61% of the variability. Other Biplots are shown in Figure C 4 and C 5 Recall, PC1 is mainly representing stand development stages, and PC2 tree to tree competition. The scores associated with each observation show a wide range of tree to tree competitio n at later stand developme nt stages (left side of Figure 2 4 a ), and for early stand development stages, this tree to tree competition is more uniform (right side of Figure 2 4 a ). This Biplot also shows cl early the negative association between AIDBH and tree competition variables; where, as expected, high individual competition means low diameter growth. Surprisingly, productivity (represented by SI) had a negligible effect in explaining the variability of this population. From this Biplot, groups of correlated variables can be easily identified: Hd Ad QD, BA BAN RS, BAL BALn BALr
29 AIDBH, and A N. No clear patterns of the distribution of the trees were detected according to growth zone and species. A second P CA was don e with only the variables from stand subset together with average plot AIDBH. Factor loadings of four components are shown in Table 2 3 The meaning of each PC, as expected, were almost identical to the ones described earlier st PC explained a 67% of the total variance (Figure C 6 ) and they are rep resented by a Biplot in Figure 2 4b Here the plot observations (scores) present a higher variability in stand competitio n at later stand development stages (righ t side of Figure 2 4b ), and vice versa The implementation of K means clustering with the subset of stand variables identified two clusters that were selected using the within cluster s um of squares and silhouette method s (Figure C 7 ) The first cluster was associated with high values of N and RS while the second with high values of BA, BAN, DA, DH, and QD (Figure 2 4b and Tabl e 2 4 ) Hence the main differences among plots were associated to stand development stages where, for one side, there are young forests with a high tree density, and on the other older forests with basal area concentrated in few individuals. A n MRPP test statistic of chance corrected within group agreement (A s ) of 0.16 was obtained for these two clusters, indicating strong differences between them In addition, the statistics of observed delta was 3.51 (P < 0.001) versus an expected delta of 4.17 indicating strong significant differences (P < 0.001) in term s of stand variables (Figure
30 C 8 ). Association between cluster s and growth zones was evaluated with a Pearson chi squared statistics with a value of 8.69 (P = 0.034) The post hoc pairwise testing only established a small association between clusters and the growth zones 1 and 2 (Table B 1 ). T hese two clusters were contrasted in a map with growth zones, where no cl e ar patterns were noted (Figure 2 5 ). Discussion For second growth RORACO forest s high levels of variability were found in the tree and stand subsets, with wide ranges across growth zones in terms of averages and standard dev iation (Table 2 1 ). This variability originates from several factors including biophysics variables ecology and species characteristics, disturbance etc. (Donoso 1993; Chauchard and Sbrancia 2003; Echeverria and Lara 2004; Lue bert and Pliscoff 2006). Stand variability is not only found between provenances but also within stands at the individual level (Donoso 1993). In Chile, t he Andes mountain range vegetation has a marked altitudinal gradient influence d by temperature and ra infall where higher altitude s are characterized by low temperature s and increased rainfall (Luebert and Pliscoff 2006). Together with this altitudinal gradient that in Chile is closely related to longitude, the latitudinal gradient strongly affects the v egetative growth period (Donoso 2003). To simplify the effect of biophysical variables on a forest resource, it is useful to define relatively homogeneous geographical growth zones. The zoning used in this study was defined by combining different sources o f information related to (Gezan and Moreno 1999): climate, soil, vegetation distribution, and site productivity. These biophysical variables have been used in many other zoning studies for Nothofagus ( Novoa and Villaseca 1989; Donoso
31 et al. 1993 ; Schlatter and Gerding 1995; Echeverria and Lara 2004; Thiers et al. 2008 ; Esse et al. 2013 ). Variability found within a stand is related to ecological factors and /or forest dynamics. The results from the unsupervised multivariate methods identified important facto rs at the tree and stand level for the RORACO forest. Both of the PCA presented a high percentage of variance explained related to stand development stages, ranging from 37% to 43%. In addition, two clear cluster s were found, which represent two conditions of stand development stages: 1) young forests with a high tree density; and, 2) adult forests with high basal area concentrated in few large trees. Other authors have reported differences bas ed on stand developmental stage. F or example, Donoso ( 199 3) fou nd, for Nothofagus changes in the diameter distribution, from unimodal to bimodal, as the stand ages, reflecting additional heterogeneity that is present in these forests. This heterogeneity is incremented by the interaction with biophysical variables suc h as climate, altitude, and latitude. The other multivariate methods, NMDS analysis and PCoA evidenced a differentiation among species, but not for growth zones. Lusk and Ortega (2003) found no differences in productivity for comparisons based on the same zoning. The results of NMDS and PCoA at stand level, showed that in the multivariate space N alpina was intermediate between the others two species, being closer to N. dombeyi than N. obliqua Other authors have also found important differences between t hese species. Chauchard and Sbrancia (2003) in relation to shadow tolerant, indicated that N. alpina is the most tolerant, followed by N. dombeyi and then N. obliqua In terms of tree to tree competition, the same authors reported that at older ages, N. d ombeyi stands
32 present higher density than N. obliqua which defines its sociological status at a j uvenile phase and remains in this strata. Donoso (1993) states that N obliqua have a strong natural thinning at 40 50 years of age due to its low ability for tree to tree competition. In the present study, tree to tree competition, as described by the PC2 explained 24% of the total variability for both subsets (tree and stand). This component was mainly described by the variables of BAL, BALr, and BALn. S evera l studies have identify a wide range of factors that influence the variability foun d in the RORACO forest type. This study, contributed to better understand ing the underlying variability for these forest through the use of multivariate analysis (NMDS, PCoA and PCA) based on tree and stand variables typically obtained from forest inventories. The findings indicated that site productivity and growth zones do not explain much the variability of this sampled population. However, stand development stages, tree t o tree competition, and size tree attributes are critical variables to explain this variability. Also, growth zones were identified to have relevant importance at the macro region level. The biophysics variables were not directly included in this study, be ing represented by growth zones so that was a limitation to t ry of explain the phenotypic variability only with tree and stand variables. However, this study was focused to understand the behavior of these variables as first step to development a growth m odel for RORACO forest type.
33 Table 2 1 Summary statistics of this study dataset for individual and stand variables in each growth zone. Zone 1 Zone 2 Zone 3 Zone 4 Variable n mean (SD) range n mean (SD) range n mean (SD) r ange n mean (SD) range DBH 321 15.4 (9.3) 5 60 317 18.7 (11.0) 5 62 194 19.5 (11.5) 5 59 276 17.2 (9.6) 5 45 H 321 14.4 (6.0) 5 35 317 17.2 (7.0) 5 38 194 18.0 (7.3) 4 45 276 15.0 (6.7) 4 39 A 321 29.5 (14.8) 8 104 317 34.0 (14.8) 10 83 194 30. 8 (14.0) 9 65 276 34.2 (12.9) 10 68 BAL 321 20.3 (15.1 ) 0 74 317 24.2 (17.2) 0 83 194 21.6 (15.9) 0 64 276 23.9 (18.9) 0 92 BALn 321 16.0 (12.6) 0 56 317 19.8 (15.2) 0 82 194 17.1 (13.3) 0 55 276 18.3 (15.2) 0 63 SS 321 2.2 (1.0) 1 4 317 2.3 (1. 1) 1 4 194 2.2 (1.0) 1 4 276 2.3 (1.0) 1 4 BALr 321 0.5 (0.3) 0 1 317 0.6 (0.3) 0 1 194 0.5 (0.3) 0 1 276 0.6 (0.3) 0 1 AIDBH 321 3.0 (2.0) 0 11 317 2.9 (2.1) 0 12 194 3.8 (2.2) 0 11 276 2.7 (1.7) 0 10 BA 42 37.9 (15.9) 10 79 46 42.9 (14.8) 16 84 29 43.2 (14.8) 19 80 41 46.9 (19.8) 11 98 N 42 2,821 (1,325) 760 5,600 46 2,337 (1,245) 320 5,000 29 2,226 (1,104) 640 4,640 41 2,794 (1,135) 760 5,560 QD 42 14.3 (5.7) 7 31 46 17.3 (6.6) 7 33 29 17.3 (6.1) 8 32 41 15.2 (4.5) 8 26 Hd 42 19.1 (6.4) 8 30 46 22.3 (21.2) 11 35 29 22.9 (7.2) 11 42 41 19.8 (6.3) 9 37 Ad 42 36.4 (16.9) 14 87 46 42.4 (14.6) 13 77 29 37.2 (13.3) 16 64 41 40.7 (12.0) 16 68 SI 42 11.8 (3.8) 4 20 46 12.5 (4.1) 4 24 29 13.7 (3.4) 7 23 41 10.2 (3.7) 2 18 BAN 42 31.8 (12.4) 10 62 46 37.2 (14.6) 14 83 29 37.4 (14.0) 8 69 41 40.3 (18.3) 9 90 SDI 42 1,054 (344) 410 1,895 46 1,089 (328) 414 1,906 29 1,081 (304) 717 1,830 41 1,245 (438) 406 2,305 RS 42 0.1 (0.0) 0.1 0.2 46 0.1 (0.0) 0.1 0.2 29 0.1 (0.0) 0.1 0 .2 41 0.1 (0.0) 0.1 0.2 Note : n, number of observations; SD, standard deviation; DBH, diameter breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larger trees ( m 2 ha 1 ); BALn, basal area of larger trees for Nothofag us ( m 2 ha 1 ); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year 1 ); BA, stand basal area ( m 2 ha 1 ); N, number of trees (trees ha 1 ); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand a ge (years); SI, site index (m); BAN, basal area of Nothofagus ( m 2 ha 1 ); SDI, stand density index (trees ha 1 ); RS, relative spacing index.
34 Table 2 2 Factor loading for the first five principal component (PC) obtained from the t ree and stand subsets together Values in parenthesis correspond to the proportion of the variance explained by each component. Variable Transformation PC1 (37%) PC2 (24%) PC3 (16%) PC4 (11%) PC5 (4%) DBH square root 0.271 0.366 H 0.304 0.263 0.151 0.220 A 0.307 0.152 0.369 0.144 BAL 0.145 0.452 0.194 BALn square root 0.103 0.438 0.115 BALr arcsin 0.454 0.220 AIDBH square root 0.315 0.199 0.277 0.561 BA square root 0.341 0.127 0.286 0.165 N 0.176 0.509 0.220 QD 0.343 0.279 0.227 Hd 0.354 0.187 0.231 0.213 Ad 0.327 0.148 0.368 SI 0.714 0.151 BAN 0.328 0.109 0.216 0.359 SDI 0.204 0.151 0.494 RS 0.243 0.104 0.336 0.188 0.493 Note : Factor loadings smaller than 0.1 in absolute value are not shown and those larger than 0.3, in absolute value, are in bold. DBH, diameter breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larger trees ( m 2 ha 1 ); BALn, basa l area of larger trees for Nothofagus ( m 2 ha 1 ); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year 1 ); BA, stand basal area ( m 2 ha 1 ); N, number of trees (trees ha 1 ); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (years); SI, site index (m); BAN, basal area of Nothofagus ( m 2 ha 1 ); SDI, stand density index (trees ha 1 ); RS, relative spacing index.
35 Table 2 3 Factor loading for the first five principa l component (PC) obtained from average plot AIDBH and stand subset Values in parenthesis correspond to the proportion of the variance explained by each component. Variable Transformation PC1 (43%) PC2 (24%) PC3 (18%) PC4 (8%) AIDBH square root 0.570 0.555 BA square root 0.429 0.263 0.164 N 0.200 0.559 0.150 QD 0.396 0.324 0.171 Hd 0.391 0.275 0.173 0.317 Ad 0.362 0.150 0.380 SI 0.180 0.668 0.261 BAN 0.406 0.149 0.372 SDI 0.279 0.514 0.101 RS 0.294 0. 308 0.226 0.543 Note : Factor loadings smaller than 0.1 in absolute value are not shown and those larger than 0.3, in absolute value, are in bold. DBH, diameter breast height (cm); H, total height (m); A, breast height age (years); BAL, basal area of larg er trees ( m 2 ha 1 ); BALn, basal area of larger trees for Nothofagus ( m 2 ha 1 ); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year 1 ); BA, stand basal area ( m 2 ha 1 ); N, number of trees (trees ha 1 ); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (years); SI, site index (m); BAN, basal area of Nothofagus ( m 2 ha 1 ); SDI, stand density index (trees ha 1 ); RS, relative spacing index. Table 2 4 Mean values o f stand subset and average plot AIDBH for each K means cluster. Variable Cluster 1 Cluster 2 Total AIDBH 4.3 3.4 3.9 BA 33.0 52.8 42.5 N 3226.3 1845.9 2563.0 QD 11.6 20.7 16.0 Hd 16.1 26.1 20.9 Ad 30.2 49.9 39.6 SI 11.6 12.2 11.9 BAN 27.4 46.0 36. 3 SDI 1029.4 1209.2 1115.8 RS 0.1 0.1 0.1 Note : AIDBH, annual increment in DBH (mm year 1 ); BA, stand basal area ( m 2 ha 1 ); N, number of trees (trees ha 1 ); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (year s); SI, site index (m); BAN, basal area of Nothofagus ( m 2 ha 1 ); SDI, stand density index (trees ha 1 ); RS, relative spacing index.
36 Figure 2 1 Non Metric Multidimensional Analysis (NMDS) using the tree subset, identifying gr owth zones (a) and species (b). a b
37 Figure 2 2 Non Metric Multidimensional Analysis (NMDS) using the stand subset, identifying growth zones (a) and dominant species (b). a b
38 Figure 2 3 Principal coord inates analysis ( PCoA ) using the tree and stand subsets together, identifying growth zones (a) and species (b). a b
39 Figure 2 4 First two PCA axes with scores and factor loadings obtained using, the tree and stand subsets ( a ) stand subset together with average plot AIDBH. The two K means clusters are also identified ( b ). a b
40 Figure 2 5 Geographical localization of plots identified according to their two K means clusters in relation to growth zones.
41 CHAPTER 3 INDIVIDUAL TREE DIAMETER GROWTH MODELS OF MIXED NOTHOFAGUS SECOND GROWTH FORESTS IN SOUTHERN CHILE Background S outhamerican Nothofagus species dominate temperate and sub Antarctic forests of Chile and Argentina from 33 to 56 south latitude. In Chile, there are nine species of this genera including three evergr eens, and six deciduous (Donoso 1987 ; Donoso 1993 ). The second growth forests of r oble ( Nothofagus obliqua (Mirb.) Oerst ), raul ( N. alpina (Poepp. & Endl.) Oer st. ), and c oihue ( N. dombeyi (Mirb.) Oerst. ) known locally as RORACO forest type, is one of the most i mportant natural mixed forests of Chile. According to inventory, there are about 1.4 million hectares of these forests, correspo nding to 10.8 % of the total native forest area of the country (CONAF 2011 ). In 2015, the saw timber from RORACO forest type was approximately 45% of the national timber production (INFOR 2016). The volume growth rates for these forests range between 6 and 10 m 3 ha 1 year 1 (Donoso et al. 2004) corresponding to some of the highest values for southamerican temperate native forests Given its large geographical area, timber market value and attractive growth rates, this resource presents high economic potenti al value within the c hilean forestry sector (Gezan et al. 2007). High variability is found in the RORACO forest type due to different factors, such as large geographical area, biological diversity with many companion species, varying (or lack of) managemen t, anthropic disturbances and site productivity (Chapter 2 ). For example, the se three Nothofagus species have different altitude ranges: N. obliqua is com monly found between 100 and 600 m a.s.l. N. alpina between 600 and 900 m a.s.l., and N. dombeyi is fr equent ly found over 900 m a.s.l. Additionally, i t is possible to find
42 different ecotones as single or m ixed composition forest s (INFOR 2013). Echeverria and Lara ( 2004) found that the main environmental factors that affect growth rates for RORACO forest ar e: longitude, climatic variables (mean annual rainfall, summer humidity index, frost free period) and edaphic variables (sand and silt contents) that account ed for 70.4% of the total variance Definition of growth zones helps to manage this variability, w ith several reported approaches for the RORACO forest type (Donoso et al. 1993; Gezan and Moreno 1999 ; Echeverria and Lara 2004 ); where within a given zone, it is expected to find similar growth patterns silvicultural treatment responses, and compani on sp ecies, among others In order to obtain economic, social and ecological sustainability of these forest, forest managers and owners need information regarding present conditions, and future state under varying silvicultural activities. Forest inventories gi ves critical information of the current stand, and growth and yield (G&Y) models provide with future conditions and help with modeling forest behavior (Shortt and Burkha r t 1996). Currently there is a lack of general G&Y models useful for the geographical range of the RORACO forest type. Some specific studies present models; however, these are limited to small geographical areas and specific conditions (Donoso et al. 1993; Gezan and Moreno 1999; Moreno 2001 ; Echeverria and Lara 2004 ; Salas and Garcia 2006 ; Ugalde 2006; Gezan et al. 2007 ). Birdsey (1990) states that G&Y models are useful tools for inventory updates and decision making regarding silvicultural management, for short or long term targets G&Y models have been used mainly for plantation or even a ges, single species forests as their complexity for mixed or uneven aged forests increases considerably. Some of
43 the challenges for t he construction of mixed forest models include: growth rates that differ among species, stands that include shade tolerant and intolerant species, and requirement for selective harvests with a conti nuous forest cover. All of these, and other factors, increase the uncertainties and rise the complexity and lower prediction accuracy quality of mixed forests G&Y models (Burkhart a nd Tome 2012 ). Numerous G&Y models approaches have been based on : stand initial conditions, species composition, objectives and needs of forest managers, among others (Burkhart and Tome 2012). These approaches range from models that only provide with aggr egate d stand variables (whole stand level) to models with information for each tree (individual tree level) (Burkhart and Tome 2012). Whole stand level are the most common G&Y models, and these can be extended to produce a stand table through the generatio n of a diameter distribution with several size classes (Vanclay 1994) Individual tree models are the highest detail level for modeling growth and yield, and they can be used as feedback for whole stand models ( Burkhart 2003 ; Burkhart and Tome 2012 ; Shortt and Burkhart 1996). Individual tree models ha ve been widely used for different forest ecosystems, including plantations, even aged mixed hardwood stands, uneven aged mixed species, tropical forests, boreal forests, and temperate forests ( Harrison et al. 1 986 ; Murphy and Graney 1998; Huang and Titus 1999 ; Andreassen and Tomter 2003 ; Pukkala et al. 2009; Nunes et al. 2011 ; Ma and Lei 2015 ; Salas et al. 201 6 ). For individual tree models the projection unit is the individual tree, with the distance dependent a nd distance independent model procedures, which differ in the inclusion of the physical coordinates of each tree within a plot. I ndividual tree models offer high flexibility to describe and project any stand structure and allow for the
44 incorporat ion of mor e reasonable process based modules, such as thinning differentiation of growth patterns and, in the case of mixed forest, responses for different species to competition and mortality, for example. However, whole stand models often give more accurately est imate s and projections of the overall stand va riables (Burkhart and Tome 2012 ; Gonzalez Benecke et al. 2012 ), but at a lower level of resolution An individual tree G&Y model consists of three components (Burkhart and Tome 2012) : 1) individual tree diamet er breast height (DBH ) or basal area growth equation; 2) individual tree height growth equation or, alternatively, a height diameter relationship; and 3) mortality or survival component. The DBH or basal area growth is modeled using a function of predictor s, tree and stand attributes including : DBH, total height and age ; dominant age, dominant height or site index (as a measure of site productivity); basal area of larger trees and sociological status (as a measure competition index); stand basal area, numb er of trees per unit area, and quadratic diameter (Huang and Titus 1995; Burkhart and Tome 2012 ; Ma and Lei 2015) For mixed forests and large areas, these growth equations are tuned up or calibrated to allow for different species and/or geographical zones (Pukkala et al. 2009). For mortality/s urvival often logistic models are used that also consider individual tree attributes, together with a measure of competition and stand attributes ( Murphy and Graney 1998 ; Nunes et al. 2011 ; Burkhart and Tome 2012 ; Ma and Lei 2015 ). Most G&Y models, and particularly, individual tree growth one s, are obtained by fitting multiple linear regressions with many potential predictors in their original units or with transformations. These predictors often suffer from multicoll inearity (Welham et al. 2015). There are several statistical methods for selecting predictors including: best
45 subset forward, backward, forward stepwise, backward stepwise, and hybrid selection procedures (James et al. 2013). All of these selection method s have the same logic, which is to fit a multiple linear model with the selected subset of predictors (Hastie et al. 2009). Alternatively, some shrinkage or regularization method has been commonly recommended, where a multiple linear model is fitted with a ll predictors simultaneously incorpora ting a constraint that shrink the model coefficient toward zero. The advantage of these procedures is related to the reduction of the variance of each coefficient estimate with this regularization (Tibshirani 1996; Has tie et al. 2009; Kyung et al. 2010 ; James et al. 2013 ). The two best known of such procedures are: ridge regression (RR) and least absolute shrinkage and selection operator (LASSO) ( Kyung et al. 2010 ; Hastie et al. 2009). RR is similar to least squares, ex cept for the incorporating of a l 2 shrinkage method is related to the bias variance trade off, where for of the coefficient estimates decre ases by inducing an, often small, bias in the estimates ( Hastie et al. 2009; James et al. 2013) The main disadvantage of RR is that it includes all predictors in the final model, which is, among other issues a problem for model interpretation related to t he predictors importance LASSO, in contrast, uses a l 1 shrinkage penalty that makes some coefficient estimates exactly zero, particularly when the tuning parameter is large; thus, LASSO at the same time performs shrinkage and variable selection (Hastie et al. 2009; Kyung et al. 2010 ; James et al. 2013 ). One limitation for LASSO, is that under non Bayesian statistics, due the penalized
46 estimation, it does not provid e P values, confidence intervals or standard e rror of the regression coefficients ( Tibshirani 1996 ; Wu et al. 2009 ; Lockhart et al. 2014) The lack of growth information for the RORACO forestry type and social and economic importance of this resource sugge st the necessity of improvement in growth models for full geographic area of these forests T he main aim of this study is to develop an individual t ree growth model to estimate annual increment in DBH for second growth forest stands of N. alpina N. obliqu a and N dombeyi including a set of tree and stand level predictors to explain diameter growth of these species. The data originates from a study realized by Universidad Austral de Chile (Ortega and Gezan 1998), over a n ecolog ically and geographical ly di verse area (> 14 ,000 k m 2 including four growth zones). The specific study objectives are: 1) to compare different multiple linear fitting procedures to obtain the best individual tree growth model; 2) to generalize diameter growth models to incorporate gr owth zone and specie, and 3) to evaluate the final diameter growth model using an independent validation dataset. Materials and M ethods Data D escription The dataset used to fit the individual tree diameter breast height (DBH) growth mo del is from a study b y the Universidad Austral de Chile (Ortega and Gezan 1998) consisting on two plot networks. One with a total of 128 permanent plots (PP), and another with 6 0 temporary plots (TP) (Figure 3 1 ). The data were collected under a population stratification, based on a national forest cadastre (CONAF et al. 1999), that considered different stand conditions, including: RORACO second growth forest type dominated by their Nothofagus spp. canopy coverage >50%, and with minimum or no sta nd disturbance.
47 Each PP had an area of 500 m 2 but for higher densities (> 4,800 trees ha 1 ) the plot area was reduced to 250 m 2 All trees with 5 cm or more in DBH and that were taller than 2 m in height were measured. T otal height was obtained in a subsa mple of 15 trees per plot and t ree increment cores were extracted from this subsample at breast height to obtain breast height age and past diameter growth The TP were formed by two ci rcular subplots with an area of 125 m 2 each All trees with DBH > 5 cm and taller than 2 m were measured. In each subplot, two trees were felled for complete stem analysis. Further details of these datasets is presented in Gezan and Moreno (1999 ). This study used the annual average DBH growth based on the last years, which o riginated from the increment cores (PP) and tree sections (TP) at breast height of the penultimate and antepenultimate growth year s. The last growth year was discarded as p lots were established at different months within a year, without all completing the current growing season The dataset included the following tree level attributes : species (SP) diameter breast height (DBH, cm), height (H, m), breast height age (A, year), basal area of larger trees (BAL, m 2 ha 1 ), basal area of Nothofagus of larger tre es (BALn m 2 ha 1 ) sociological status (S S defined according to vertical stratification with dominant (1) codominant (2) intermediate (3) or suppressed (4)) and relative basal area of larger trees (BALr, with BALr = BAL /BA) The response variable as indicated earlier, corresponded to the average a nnual increment in DBH (AIDBH, m m year 1 ). Summary statistics of the above att ributes are presented in Table 3 1 The stand level variables considered that characterized each plot were: growth zone (Zone), stand basal area (BA, m 2 ha 1 ), trees per hectare (N, trees ha 1 ), quadratic
48 mean diameter (QD, cm), dominant height (Hd, m), dominant breast height stand age (Ad, years), dominant specie (DS), site index (SI, m), basal area of N othofagus (BAN, m 2 ha 1 ), stand density index (SDI, trees ha 1 ) and relative spacing index (RS). BAN was calculated as the sum of basal area of N. alpina, N. obliqua and N. dombeyi Hd and Ad were calculated using the 100 thickest trees per hectare, and D S corresponded to the Nothofagus specie with the largest proportion of BAN. SI was estimated by using available dominant height site models obtained on the same dataset for an index age of 20 years (see appendix A). QD) (Avery and 1.4112 as reported by Gezan et al. (2007 ) fitted using both the PP and TP networks. Finally, RS was calculated according to: RS = [(10000/N) 0. 5 ]/Hd (Avery and Burkhart 2002; Contreras et al. 2011). Hence, a total of 1,108 trees constituted the final dataset for fitting containing tree and sta nd attributes for a total of 158 plots with a single measurement. The validation or projection database was formed by 20 plots belonging to the PP network also established in 2000, re measured in 2006, and a subset of seven plots was measured a third time in 2012. For 2006 there were a total of 1,568 trees with measurements, and for 2012 this number dropped to 435 trees. Note that the dominant species for these plots are N. ob liqua (18 plots) and N. alpina (2 plots). Model Fitting Multiple linear regression, using several potential predictors, were fitted for estimation of AIDBH. Two selection variable procedures were evaluated: 1) model selection with cross validation (CV), an d 2) LASSO regression. The resulting selected models were later extended using regression with groups to identify which of the
49 factors growth zone or specie contributes to improve the performance of these models, and therefore, explain better the biologi cal process. AIDBH could be influenced by different elements, such as site productivity, tree to tree competition, together with individual tree and whole stand attributes (Andreassen and Tomter 2003 ). The fitted growth model form used in this study was: l n(AIDBH) = f(productivity, competition tree, stand) (3 1) W here, AIDBH is the annual increment in DBH; productivity was represented by Hd and SI; competition predictors were SS, BAL, BALn, BALr, SDI, and RS; tree attributes corresponded to DBH, A, and H; and whole stand attributes included DA, BA, QD, and BAN. Hence, these corresponded to a total of p = 15 predictors. The original predictors and transformations for each of these 15 predictors were considered, these were: natural logarithm, i n verse, inverse square root and quadratic term. Some of these transformations included a constant c = 10, particularly for logarithm transformation. The above transformations were selected as these are the mo st common one reported for individual tree growt h model ( Shortt and Burkhart 1996 ; Murphy and Graney 1998; Huang and Titus 1999 ; Andreassen and Tomter 2003 ; Pukkala et al. 2009; Contreras et al. 2011; Nunes et al. 2011 ; Pukkala et al. 2011; Ma and Lei 2015 ) The model selection procedure with CV regress ion, was setup by using a K fold cross validation scheme that divided the fitting database into K = 10 random groups ( or folds ) of similar size (James et al. 2013) As a first step, a model with all predictors and its transformations (p = 67 ), was evaluate d using the CV procedure, subsequently, a single expression of each predictor was chosen by backward selection and by doing so,
50 a temporary model containing p = 1 5 predictors was available for the next step. Then, a new CV procedure was run with these 15 p re selected predictors to obtain a final candidate model. The best CV model was selected according to the following goodness of fit statistics: a djusted r squared ( R 2 C p information criterion (BIC), residual sum of squares (RSS), and the mean squared error averaged over the K folds (MSE K ) (Lumley and Miller 2009 ; James et al. 2013 ). The second procedure used for model selection and model fit was LASSO regression (James et al. 2013). The was determined by cross v alidation using the fitting database with the mean square error as criteria to obtain the final candidate model. Once the potential predictors were identified for the CV and LASSO regression procedures, the fitting database was then split 50/50 into traini ng and test dataset. The training dataset was used to obtain the model parameters. A model correction factor (Baskerville 1972) was incorporated due to the use of the natural logarithmic of the response AIDBH, with 0 = 0 exp (MSE/2) is the adjusted intercept, and MSE is the mean square error for a given model ob tained based on the training dataset. To evaluate the predictions performance of the models the following goodness of fit statistics were used on the test data ( Huang et al. 2003; Welham et al. 2015): empirical coefficient of correlation (R 2 emp), root mean square error (RMSE), relative root mean square error (RMSE%), bias (BIAS) and relative bias (BIAS%). Their formulae are: (3 2) (3 3) (3 4)
51 (3 5 ) (3 6) W here, a nd are the i th observation and prediction value, respectively; n is the number of observations, is the mean of the observed values. All of the summations go from 1 to n. These statistics were obtained for the complete test dataset, and additionally, th ey were calculated by dividing this dataset into three DBH classes, defined as: 5 15 cm, 15 30 cm, and >30 cm. In order to generalize the candidate diameter growth model from the CV procedure, species and growth zones were incorporated as both independent factors and as interaction among predictors. The species factor was defined with a label for each of the three Nothofagus species. The four growth zones correspond to those defined by Gezan and Moreno (1999) for the same population (see Fi gure 3 1 ). Each factor was evaluated separately or as a combined specie zone factor. The fitting of the models was done with the 50/50 training and test data set, and the goodness of fit statistics described before were also calculated. The sign ificance and selection of each factor, or its interaction with a given predictor, were assess through an F test as reported in the analysis of variance, and model terms were eliminated using a backward selection procedure. Later, the final selected factors were further evaluated by obtaining their least squares means and all pairwise comparisons using Tukey adjustment. For all analyses, a significance level of 5% was used. In the case of the LASSO model candidate dummy variables were defined for each of t he levels of specie zone. A final model was selected by running LASSO regression with a new tuning parameter determined by CV using the fitting database.
52 In order to verify the levels of multicollinearity between the selected predictors in the CV and LASSO regression models variance inflation factor (VIF ) and generalized variance inflation factor (GVIF) (Fox and Monette 1992) were calculated for the models with and without factors, respectively. All predictors with VIF or GVIF greater than 10 were removed from the se four model s, as they were considered highly collinear with other predictors P erformance of the CV and LASSO regression models, with and without groups, was evaluated using an independent validation with the projection database and model simula tions of 6 and 12 years since plot establishment. Here, AIDBH estimates and predictors were updated yearly until projection age. Companion species in this forest type that lack a fitted model, were projected with the same models used for Nothofagus and B ALn values for companion species were assumed to be all equal to the minimum Nothofagus BALn value. Normality of residuals and heterogeneity of variances was assessed and no important departures from these assumptions were noted. Database management and st atistical analyses were all performed using the statistical software R v 3.2.2 (R Core Team 2015) The R packages used for the CV regression was leaps (Lumley and Miller 2009), and for LASSO regression glmnet (Friedman et al. 2016). Results Model Fitting The rates of annual DBH growth for the different species obtained from the fitting database corresponded to averages of 2.6, 3.0 and 3.6 mm for N. alpina N. obliqua and N. dombeyi respectively. These rates vary greatly according to growth zone (see Tabl e
53 B 2 ), where, zone 3 had the highest annual DBH growth average with 3.8 mm and the zone 4 had the lower growth with 2.7 mm. The CV regression procedure of variable selection resulted in five predictors of AIDBH, these were sel ected according to the goodness of fit statistics R2adj C p BIC, RSS, and MSE K that consistently chose five predictors (Figure 3 2 ) All model parameter estimates were significant at P < 0.001, and its estimated par ameters ar e presented in Table 3 2 The resulting expression of the model is: (3 7) W here to are the parameters, e is the erro r term, and the other terms were previously described. The validation of the prediction of CV selection model resulted in a R 2 emp of 0.56 with a RMSE% of 44.16% an d a BIAS% of 1.96% (see Table 3 3 ). Thus, the predictions have moderate quality (large RMSE%), but they have negligible bias (Figure 3 3 ) The model selected by LASSO regression, at the beginning had eight predictors; however, after observing the VIF values, the final LASSO model resulted in seven predictors. by cross validation was 4 In relation to the CV regression model, LASSO incorporated the predictors BALr and DA and the expression of this model is: (3 8)
54 W here to are par ameters to estimate (see Table 3 2 ). This model showed sim ilar goodness of fit statistics as t he CV regression, for the predictions based on the validation data, with R 2 emp = 0.57, RMSE% = 44.16 and BIAS% = 2.94 To extend the diameter growth model to incorporate the factors of specie and growth zone, the CV model was fitted with these factors and their interactions against all predictors. These factors were added individually ( i.e. Zone, SP) and using a combination of both as a unique new factor (SpZone). A backward selection procedure eliminated all of these interactions, resulting in models tha t only describe different intercepts for each of the levels of these factors. The best model corresponded to the combined factor of SpZone (CV + SpZone, see Table B 3 ): (3 9) ij is the intercept parameter for the i th specie in the j th zone, and to are slope parameters. S pZone i j are the dummy variables associated with the ij th group. The s pecie were coded as: 1 ( N. alpina ) 2 ( N. obliqua ), and 3 ( N. dombeyi ) The other model terms were described previously. Note that with the incorporations of S pZone i j the predictor SDI left the model as it was not significant (P = 0.297). The goodness of fit statistics for this extended model are R 2 emp = 0.56, RMSE% = 44.49%, and BIAS% = 1.96%. This model had similar performance than the CV model (Table 3 3 ); however, the factor SpZone resulted highly significant (P < 0.001).
55 Equivalently, the LASSO regression model was extended wi th the factor SpZone (LASSSO + SpZone) using a new tuning parameter ( = 3.9 10 3 ). Here, the final model is: (3 10) Where to are parameters to estimate, and the other terms were described previously. The goodness of fit statistics for this model are R 2 emp = 0.54, RMSE% = 45% and BIAS% = 4.31%, values that are lower in performance with respect to other three previous models evaluated (Table 3 3 ). Note that standard errors and P values for LASSO and LASSO + SpZone are not reported because it is not possible to obtain a reliable unbiased estimate of the standard error ( Kyung et al. 2010; Lockhart et al 2014). The 0 based on Baskerville (1972) correction of CV, LASSO, CV + SpZone, and LASSO + SpZone models were 2.546, 0.160, 2.829, and 1.239, respectively. The evaluation of the four predictive models of AIDBH across different DBH classes (5 15, 15 30, >30 cm) p ro vided some i nteresting differences (Figure 3 4 ). For RMSE%, the CV, LASSO, and CV + SpZone models had similar performances with lower values for the intermediate class (15 30 cm). LASSO + SpZone model had better performance for the other two classes; however, it showed much larger RMSE% in the intermediate class. In the case of BIAS%, the four models resulted in similar performance across all DBH classes. It appears that the best model corresponded to CV + SpZone; however, i t presented larger negative BIAS% for the class >30 cm, and the LASSO regression models had the smallest bias for the same class.
56 Model Projection Model validation usin g the projection dataset found good results for all models with R 2 emp >0.97 (Table 3 4 ). CV and LASSO regression had better goodness of fit statistics, while LASSO + SpZone had the worst performance. A projection of 6 years showed that LASSO regression had a better R 2 emp = 0.99 and RMSE% = 7.26 but with more BIAS% = 0.36 than the CV model, although all models had excellent statistics. CV and LASSO regression for a simulation of 12 years continued with good prediction performance. CV regression presented R 2 emp = 0.97, RMSE% = 9.66 and BIAS% = 1.75, at the same proje ction time, with similar results for LASSO regression with R 2 emp = 0.97, RMSE% = 9.50 and BIAS% = 1.96. Surprisingly, t he extended models with species and growth zone factors were not superior as those without them Selected Model LASSO + SpZone model show ed poor performance for fitting and projection, with a negative bias. In contrast, the CV, LASSO, and CV + SpZone models presented similar performance in terms of goodness of fit statistics. However, LASSO regression had a larger bias and, in addition, inc luded two more predictors than CV. Hence, CV model could be considered the best model, but the analysis of variance for CV + SpZone showed that the factor SpZone is highly significant with P<0.001. Thus, the selected model for this study correspond to CV + SpZone. Based on this mode l the incorporation of species and growth zone confirm the differences reported earlier about DBH growth from the fitting database (see Table B 2 and C hapter 2 ). An asses sment of the selected model corroborates forest ecological processes and shows biological consistency in the parameter estimates (Table 3 2 ) For example, BAL n and S S are competition variables where a higher value indicates a hi gh
57 competition for a given tree, and their model parameters have both negative si gn; hence, with higher competition there are lower rates of growth In the case of the inverse of BALr the estimate d coefficient value is positive thus a high competition pr oduces a smaller growth rate, and vice versa. The most important predictor, DBH has a strong effect on growth, where bigger trees have higher diameter increment because due to its positive sign, which is expected as this effect is c ontrolled by the predict or age. Hence, at a fixed age, larger trees tend to have greater growth rates. Finally, the parameter associated with age is negative, thus for a fixed DBH value, a young tree is capable of more growth than a n older tree For the CV + SpZone model, the le ast squares means showed that there are important differences among the levels of this combined factor (Table 3 5 ) First, these results show similar DBH growth for N. dombeyi across all growth zone s Second, N. obliqua also ha s a similar response in all zones with the exception of zone 3, this is probably a result of been the wettest zone with long growing season Third N. alpina has the largest differences between zones where growth seems to increase with altitude ( e.g. zon e 4 Andes mountain ranges has the largest growth ) These results, justify the inclusion of the combined factor SpZone in the final selected model: CV + SpZone. The validation of the projection with independent data for the CV + SpZone model shows high pr edict ive power for estimating future DBH for a period of 6 years; and even if the projection is 12 years resul ts are very consistent (Figure 3 5 ). Discussion Diameter growth models have been recognized as a crucial component f or projecting individual tree growth ( Vanclay 1994 ; Avery and Burkhart 2002; Burkhart and
58 Tome 2012 ). In this study, a DBH growth distance independent model was development for mixed Nothofagus second growth forest in southern Chile. Individual tree models present good results for short term projections, with a similar level of accuracy respect to whole stand models ( Buckman 1962 ; Burkhart 2003 ). However, these models often have poorer levels of accuracy for long term projections, where whole stand models a re more reliable. Individual tree models are more useful to project changes in the forest structure, since the diameter distribution is directly obtained from the input data. Estimation of timber products, stand structure, differentiation by specie and/or strata, and incorporations of silvicultural methods are also improved with individual tree models (Andreassen and Tomter 2003 ). In relation to the distan ce dependent models, the distance independent models are simpler to implement given that they do not ne ed the physical coordinates of each tree within a plot; thus, they can be used with the information obtained from a typical forest inventory. Salas et al. (2016) mentioned that both procedures do not have significate differences in growth prediction, even in natural forest s ; furthermore spatial dependency is more complex to model. Therefore, it is reasonable to aim at developing a distance independe nt individual tree growth model in this study for Nothofagus forests. The response variable for this study w as the annual average DBH growth obtained over a period of two years that originated from increment cores and tree sections from felled trees. Several crucial elements are associated with the use of this response variable. First, a short time period of two years was chosen in order to have better association of observed tree growth with its present tree and stand condition s, a relationship that get weaker as the periods get longer. This is especially present, on
59 those years with particular conditions ( e.g. climatic, diseases, catastrophic event); however, Andreassen and Tomter (2003) found no effects of the prediction accuracy with different growth period lengths. Second, no evidence of wood compression was found in increment cores when compared to sections (d ata not shown). Finally, the available individual basal area growth was discarded as response as it presented poorer goodness of fit statistics than DBH growth. Individual tree growth models that include predictors of competition, productivity, tree a nd stand level have been used in many reported studies ( Harrison et al. 1986; Murphy and Graney 1998; Huang and Titus 1999; Peng 2000; Andreassen and Tomter 2003 ; Pukkala et al. 2009; Nunes et al. 2011 ; Ma and Lei 2015). Hence, there are potentially many p redictors that might describe the same process, which will suffer from multicollinearity producing inconsistences and unexpected results in the fitted models (Welham et al. 2015). Thus, it is necessary to select a subset of these predictors using some vari able selection method In this study, the fitting database contained, as expected, pairs of predictors with high correlation, such as BA BAN, DH RS, BA SDI, BAL BALn, DH QD, A DA, RS QD, and DBH H, with values >0.83 (Figure C 1 ), and the variable selection methods implemented here were CV and LASSO regression. The predictors selected by CV and LASSO regression corresponded to competition, tree and stand level variables, but they did not included variables associated with produc tivity. In CV regression DBH, A, SS, SDI, and BALn were the selected predictors, while in LASSO regression these were the same predictors with the addition of DA and BALr. The tree level variables selected corresponded to DBH and A, both correspond to a me asure of tree size; where for DBH, growth rates increase with larger diameters;
60 however, age presented a negative parameter, indicating that at a fixed DBH value, older trees present lower growth rates. This agrees with a similar result reported by Cubillo s (1987) in N. alpina In this study, competition variables had a high relevance in explaining individual tree growth, as found in other Nothofagus studies (Cubillos 1987 ; Monserud and Sterba 1996 ; Coomes and Allen 2007 ; Richardson et al. 2011 ; Ivancih et al. 2014). Interestingly, the predictor BALn (basal area of larger trees only for Nothofagus ), was selected instead of BAL (that includes all species), agrees with findings from previous studies that support the theory of additive effects in growth for No thofagus and companion species (Lusk and Ortega 2003; Donoso and Lusk 2007; Donoso and Soto 2016), that translates into independent behavior (and therefore models) for each of these cohorts. In contrast productivity predictors such as SI, were not select ed for modelling DBH growth in these forests. The average DBH growth rates for combinations of specie a nd zone varied (Table B 2 ), with important differences between species and zone; however, extreme cases where found for the c ombinations of N. alpina zone 1, and N. dombeyi zone 3, with an 88% larger growth rate for the latter one ( i.e. 2.15 mm against 4.05 mm). Other studies also have found differences in growth between Nothofagus species and zones (Donoso et al. 1993 ; Sch latter et al. 1994 and 1995 ; Donoso et al. 2004 ; Echeverria and Lara 2004; Gezan et al. 2009 ; Esse et al. 2013). The incorporation of species and/or zone factors to growth model s has been previously reported (Andreassen and Tomter 2003 ; Pukkala 2009). In t he present study, the combined specie zone factor (SpZone with 11 levels), explained more variability than the individual factors (with 7 levels for both specie and zone), indicating that a
61 greater disaggregation of the data is required to model growth rat e. This agrees with the average growth differe nces reported earlier (Table B 2 ). Interestingly, the LASSO + SpZone model selected as predictors for the combined factor SpZone the extreme groups (Table 3 5 ), by combining the central five groups into a single group. The fitting of the four candidate models to predict AIDBH resulted in good overall results, with R 2 emp ranging from 0.54 to 0.57, RMSE% of ~44%, and a BIAS% smaller than 3%. Similar goodness of fit statistics have been reported for AIDBH model with R 2 emp ranging from 0.26 to 0.68 ( Wykoff 1990; Monserud and Sterba 1996; Andreassen and Tomter 2003 ; Nunes et al. 2011 ). AIDBH projections for 6 and 12 years resulted in R 2 emp ranges of 0.23 0.28 an d 0.15 0.24, respectively. This drop in the goodness of fit statistics is been affected by the poor quality of the projection database, that included, for example, at 6 years 20% of the trees with null or negative DBH increments. Th is is likely to be assoc iated with field measurement errors and low growth rates for th ese species. In contrast, DBH projections for 6 and 12 years result in overall R 2 emp values greater than 0.97 for both Nothofagus and companion species. Small diameter trees (DBH < 15 cm) pres ented lower correlations, with values of 0.88 and 0.58 for projections of 6 and 12 years, respectively. As expected, goodness of fit statistics were better for the Nothofagus cohort than the companion species cohort, because the models were built using onl y measurements from Nothofagus T he final selected model in this study is CV + SpZone (Table 3 3 ). The LASSO model was discarded, even that was highly competitive, due to lower prediction performance and it also included some p redictors with high multicollinearity (VIF > 5).
62 LASSO + SpZone had the worst performance with issues for prediction and projection. The CV model presented better goodness of fit statistics than the CV + SpZone model; however, the combined factor SpZone wa s highly significant (P < 0.001). Hence, the inclusion of SpZone is justified given the wide geographical range of this population and its differentia l diameter growth (see Table B 2 ). These results has led to many authors to de fine growth zones (Donoso et al. 1993; Gezan and Moreno 1999 ; Echeverria and Lara 2004 ). Alternatively, it is possible to use the CV model that does not require zone or specie to have reliable prediction of DBH growth. The advantages of CV as a predictor s election procedures is the possibility of providing for a direct estimate of test error (MSE K ) compare to stepwise methods. LASSO is an alternative procedure that regularizes the coefficient estimates, forcing some of them to be exactly equal to zero, perf orming variable selection. Both procedures have been used commonly due their advantages to identify those variables that better predict a desired response (James et al. 2013 ) Similarly, t he importance of incorporate dummy variables (or factors) in a regre ssion is related to the model specificity, where the factors can correct intercepts or slop e from the original regression. Given the high levels of variability present in these mixed natural forests (Chapter 2 ), the above results are promising due to high complexity of these ecosystems, including several strata, species, anthropic disturbances, and stage structure (Donoso 1993) However, future research is warranted to better improve the performance of model s of individual tree growth for this population. S ome of these studies could include : 1) models that do not require age, as this predictor is often difficult to obtain from ty pical forest inventories; 2) increase the additional validation data from a wider
63 geographic range and with addi tional forest condi tions; 3) develop an individual tree growth mod el for companion species; 4) fully implement an individual tree model routine that includes a mortality and a height growth module; and 5) develop a compatible system based on both individual tree and whole s tand models. The results from this study, on the evaluation of individual tree DBH growth model indicated that CV + SpZone model is useful for practical applications; although, this model should be used within its frame of inference: RORACO second growth forest, between 37 and 41.5 south latitude with breast height ages ranging from 20 to 80 years, and with a Nothof agus basal area greater than 60%. It is recommended the use of this proposed model for projections as long as 12 years, where this individua l tree DBH growth model showed to be a robust tool with reasonable predictions. A comparison of two selection predictor procedures to estimate annual increment in DBH was presented using data from a second growth mixed forest conformed by N alpina N obl iqua and N. dombeyi in southern Chile. A selected small set of predictors successfully explain ed a large portion of DBH growth with variables associate competition, tree and stand level components; however, no variables associated with productivity wer e included The final model constitutes a simple and valuable tool to support management decision for this ecosystems. This study, is the first complete effort to obtain an individual tree growth model for DBH useful for the full geographical distribution of the RORACO forest type.
64 Table 3 1 Summary statistics of individual and stand attributes for fitting database based on PP and TP networks. Tree Specie N. alpina N. obliqua N. dombeyi Tree variables n mean SD min max n mean SD min max n mean SD min max DBH 202 16.7 9.8 5.0 42.7 627 17.0 10.0 5.0 62.1 279 19.3 11.3 5.0 60.2 H 202 15.2 6.6 4.2 34.5 627 16.1 7.3 3.5 45.0 279 16.2 6.1 4.2 34.5 A 202 36.7 14.2 11.0 104.0 627 30.6 14.6 8.0 95.0 279 32.5 13.2 9.0 80.0 BAL 202 26.0 16.8 0.0 83.4 627 19.0 14.0 0.0 65.1 279 28.0 20.6 0.0 92.4 BALn 202 20.9 15.0 0.0 82.3 627 15.7 12.3 0.0 54.9 279 20.6 16.6 0.0 82.2 SS 202 2.3 1.1 1.0 4.0 627 2.3 1.0 1.0 4.0 279 2.2 1.0 1.0 4.0 BALr 202 0.6 0.3 0.0 1.0 627 0. 6 0.3 0.0 1.0 279 0.6 0.3 0.0 1.0 AIDBH 202 2.6 1.6 0.2 7.7 627 3.0 2.1 0.2 12.1 279 3.6 2.2 0.1 10.2 Dominant Specie N. alpina N. obliqua N. dombeyi Stand variables n mean SD min max n mean SD min max n mean SD min max BA 24 46.7 13.8 10.8 72 .3 87 35.7 13.0 9.5 71.5 47 55.2 17.8 23.4 98.4 N 24 2,564 1,184 320 5,000 87 2,378 1,119 320 4,640 47 2,900 1,357 880 5,600 QD 24 16.4 4.8 8.5 30.5 87 15.4 6.2 6.8 33.3 47 16.9 5.4 8.4 30.4 Hd 24 21.2 5.7 10.2 35.1 87 20.8 7.0 7.8 42.4 47 20 .8 5.9 9.9 34.1 Ad 24 47.5 10.9 23.0 77.0 87 36.6 14.7 12.7 86.8 47 40.2 13.1 21.3 85.1 SI 24 10.3 3.9 4.1 24.3 87 12.5 3.9 2.0 22.9 47 11.4 3.4 3.9 18.4 BAN 24 37.6 12.4 10.8 59.7 87 30.4 11.5 8.8 57.8 47 46.5 18.7 8.0 89.6 SDI 24 1,212 336 406 1,944 87 956 258 410 1,674 47 1,400 361 640 2,305 RS 24 0.1 0.0 0.1 0.2 87 0.1 0.0 0.1 0.3 47 0.1 0.0 0.1 0.2 Note : n, number of observations; SD, standard deviation; min, minimum; max, maximum; DBH, diameter breast height (cm); H, total h eight (m); A, breast height age (years); BAL, basal area of larger trees ( m 2 ha 1 ); BALn, basal area of larger trees for Nothofagus ( m 2 ha 1 ); SS, sociological status; BALr, relative BAL; AIDBH, annual increment in DBH (mm year 1 ); BA, stand basal area ( m 2 ha 1 ) ; N, number of trees (trees ha 1 ); QD, quadratic diameter (cm); Hd, dominant height (m); Ad, dominant breast height stand age (years); SI, site index (m); BAN, basa l area of Nothofagus ( m 2 ha 1 ); SDI, stand density index (trees ha 1 ); RS, relative spacing index.
65 Table 3 2 Parameter estimates of the four selected models for estimation of the logarithm of AIDBH (mm). Parameter Estimate SE P value VIF* Parameter Estimate VIF* CV (eq. 3 7) LASSO (eq. 3 8) 2.410 0 2.173 1 < 0.001 1 3 3 < 0.001 1.74 3 3.2 4 5 < 0.001 1.22 4 1.4 1 2 < 0.001 3.29 1 5.28 1.138 0 2 < 0.001 2.26 0 6.73 1 2 < 0.001 2.22 1 2.44 2 5.36 0 6.59 CV + SpZone (eq. 3 9) LASSO + SpZone (eq. 3 10) 2.702 0 2.562 1 < 0.001 2.09 0 2.908 0 2.345 1 < 0.001 2.09 2 2.28 3.065 0 2.397 1 < 0.001 2.09 1 2.28 2.538 0 2.120 1 < 0.001 2.09 1 2.28 2.587 0 2.219 1 < 0.001 2.09 1 2.28 2.841 0 2.272 1 < 0.001 2.09 2 2.28 2.678 0 2.329 1 < 0.001 2.09 2 2.28 2.946 0 2.298 1 < 0.001 2.09 3 3.30 2.948 0 2.680 1 < 0.001 2.09 5 2.07 2.941 0 2.353 1 < 0.001 2.09 1 6.21 2.902 0 2.358 1 < 0.001 2.09 0 7.30 6.517 3 2.002 3 0.0012 1.82 1 2.53 9.307 1 6.970 2 < 0.001 3.82 2 5.56 1.175 0 7.797 2 < 0.001 2.61 1 6.91 1.401 1 3.092 2 < 0.001 2.29 Note : SE, Standard error; VIF*, variance inflation factor or generalized variance inflation factor. Table 3 3 Goodness of fit statistics of four models for ln(AIDBH) using test data. Model N R 2 emp RMSE RMSE% BIAS BIAS% CV 551 0.56 1.35 44.16 0.06 1.96 LASSO 551 0.57 1.35 44.16 0.09 2.94 CV + SpZone 551 0.56 1.36 44.49 0.06 1.96 LASSO + SpZone 551 0.54 1.36 45.13 0.13 4.31 Note : SpZone, combination of Specie and Zone; R 2 emp, empirical coefficient of correlation; RMSE, root mean square error; RMSE%, relative root mean square error; BIAS%, relative bias.
66 Table 3 4 Projection goodness of fit statistics for DBH (cm) in a period of 6 and 12 years using projection database. Projection = 6 years CV CV + SpZone n R 2 emp RMSE% BIAS% n R 2 emp RMSE% BIAS% Total 1455 0.98 7.32 0.18 1455 0.98 7.44 0.54 DBH(5 15) 777 0.88 10.06 0.42 777 0.89 9.75 0.00 DBH(15 30) 510 0.88 6.94 0.94 510 0.87 7.27 1.46 DBH(>30) 168 0.97 4.11 0.46 168 0.97 4.19 0.32 Nothofagus 943 0.99 6.33 0.26 943 0.99 6.48 0.31 Companion 512 0.96 10.47 1.46 512 0.96 10.47 1.20 LASSO LASSO + SpZone n R 2 emp RMSE% BIAS% n R 2 emp RMSE% BIAS% Total 1455 0.99 7.26 0.36 1455 0.98 7.32 0.60 DBH(5 15) 777 0.88 9.96 0.52 777 0.89 9.85 0.31 DBH(15 30) 510 0.88 6.94 1.13 510 0.88 7.12 1.51 DBH(>30) 168 0.97 3.96 0.13 168 0.97 4.02 0.03 Nothofagus 943 0.99 6.28 0.05 943 0.99 6.38 0.36 Companion 512 0.96 10.29 1.54 512 0.96 10.38 1.29 Projection = 12 years CV CV + SpZone n R 2 emp RMSE% BIAS% n R 2 emp RMSE% BIAS% Total 389 0.97 9.66 1.75 389 0.97 9.82 2.34 DBH(5 15) 177 0.58 17.07 5.03 177 0.59 16.87 6.02 DBH(15 30) 151 0.83 8.39 2.10 151 0.81 8.86 2.70 DBH(>30) 61 0.89 5.55 1.36 61 0.89 5.60 1.00 Nothofagus 295 0.98 7.93 0.29 295 0.97 8.17 0.53 Companion 94 0.83 18.45 12.49 94 0.83 18.20 11.92 LASSO LASSO + SpZone n R 2 emp RMSE% BIAS% n R 2 emp RMSE% BIAS% Total 389 0.97 9.50 1.96 389 0.97 9.50 2.39 DBH(5 15) 177 0.58 17.07 4.74 177 0. 59 16.87 5.13 DBH(15 30) 151 0.83 8.30 2.42 151 0.82 8.49 2.80 DBH(>30) 61 0.90 5.30 0.79 61 0.90 5.27 0.37 Nothofagus 295 0.98 7.74 0.05 295 0.98 7.83 0.62 Companion 94 0.83 18.45 12.41 94 0.83 18.04 11.76 Note : CV, cross validation procedure; Sp Zone, combination of Specie and Zone; DBH, diameter breast height (cm); R 2 emp, empirical coefficient of correlation; RMSE%, relative root mean square error; BIAS%, relative bias.
67 Table 3 5 Least square means (standard errors i n parentheses) statistics of all pairwise Tukey comparisons between all 11 groups based for Specie Zone in CV + SpZone model. N. alpina N. obliqua N. dombeyi zone 1 0.794 (0.125) abcd 0.652 (0.057) a 1.020 (0.079) cd zone 2 0.996 (0.072) cd 0.697 (0.053) ab 1.027 (0.154) abcd zone 3 0.948 (0.063) bcd 1.034 (0.090) cd zone 4 1.153 (0.086) d 0.784 (0.069) abc 0.960 (0.087) abcd
68 Figure 3 1 Plots distribution in the study area (sou thern C hile ) for the permanent plots (PP) and temporary plots (TP) networks
69 Figure 3 2 Statistics C p BIC, RSS R2adj and MSE K for different number of predictors with CV regression selection model using fitting database.
70 Figure 3 3 Fit and residuals plots of four selected models using test dataset. CV (a and b), LASSO (c and d), CV + SpZone (e and f), and LASSO + SpZone (g and h).
71 Figure 3 4 RMSE% (a) and BIAS% (b ) by DBH class for the four selected models evaluated using the test dataset. Figure 3 5 Simulation prediction in future DBH (cm) for projection periods of 6 years (a) and 12 years (b) based on the CV + SpZone model using pro jection database. a b a b
72 CHAPTER 4 CONCLUSIONS For the second growth mixed RORACO forest type, one of the most important natural forest resources in Chile, currently forest managers and owners do not have enough information and tools to manage this resource ens uring sustainability. This study used data from a total of 128 permanent plots and 60 temporary plots established under a stratification that covers the full geographical range of this forest type in Chile The variability found from the unsupervised multi variate methods identified important variables at the tree and stand level for the RORACO forest. The PCA analysis presented a high percentage of variance explained with components related to stand development stages tree to tree competition, and stocking and stand competition. Also, site productivity variables did not explain much of the variability in this dataset The other ordination methods used, NMDS and PCoA provided evidence of important differentiation among species but not between growth zones The results with stand level variables showed that in the multivariate space ; N alpina was intermediate between the others two species, being closer to N. dombeyi than N. obliqua In addition, two cluster s were enough to separate the sampled populatio n. T he main differences among plots were associated to stand development stages where, there are young forests with a high tree density, and where the other older forests with basal area concentrated in few individuals. The results from the K means clust er analysis also identified that g rowth zones had a relevant impor tance at the macro region level, particularly zones 1 and 2 that are statistically different
73 Supervised methods were successful to model and estimate annual growth in DBH for the combined data where the selected predictors corresponded to the variables obtained from traditional forest inventories. Both of the selection predictor procedures evaluated cross validation (CV) and LASSO regression, showed similar goodness of fit statistics with almost identical variables selected. In both cases, a small set of predictors successfully explain ed a large portion of DBH growth with selected variables associate d to competition, tree and stand level components; however, no predictors were selected as sociated with productivity Th e latter agrees with the unsupervised methods findings where productivity variables had a low contribution to explain the variability of these forests. In contrast, as found with unsupervised methods, tree to tree competition variables had a high contribution to explain variability and individual tree growth. The selection of specific predictors associated only with Nothofagus (such as basal area of larger trees for Nothofagus BALn), instead of those generic predictors that include all species (such as basal area of larger trees, BAL) agrees with others studies that support t he theory of additive effects where Nothofagus growth is independent of the companion species growth Similarly, th e incorporation of species and growt h zones as factors into the models to estimate annual growth in DBH resulted in an increase of the specificity of these regressions. This is reasonable given the wide geographical range of this sampled population and its differential diameter growth by spe cie Also, t he final selected model corresponded to a C V with a factor of specie and zone combined, where all estimated parameters showed reasonable biological
74 interpretation. This model showed moderate goodness of fit statistics for both, the prediction a nd validation datasets. Further developments are required to provide with a complete set of tools to manage this resource properly Some of these include: 1) increasing the additional validation data from a wider geographic range and with additional forest conditions; 2) developing models that do not require age, as this predictor is often difficult to obtain from typical forest inventories; 3) constructing an individual tree growth model for companion species; 4) implementing a full individual tree model r outine that includes mortality and height growth module s ; and 5) creating a compatible system based on both individual tree and whole stand growth models. In summary, the final model obtained from this study, constitutes a simple and valuable tool to suppo rt management decision for these mixed natural forests that present high levels of variability. T his is the first complete effort to develop an individual tree growth model for DBH that is useful for the full geographical distribution of the RORACO forest type in Chile
75 APPENDIX A DOMINANT HEIGHT SITE MODEL The nonlinear models of Chapman Richards are commonly used in others Nothofagus studies (Moreno 2001 ; Salas and Garcia 2006 ; Ugalde 2006). The dominant height site model 1 for this study come from the da ta originated in the project reported by Ortega and Gezan (1998). The model presented below was fitted with the same database, target species, and growth zones of this study (see Table A 1 ). The final selected model has the fol lowing general equation: (A 1) Where, = dominant height (m); = dominant age (years); = site index to 20 years (m); and , are the estimated parameters presented in Table A 1. Table A 1. Parameter estimate s of the dominant height site model by specie and growth zone. Specie Zone a N. alpina 1 40.108798 0.8756364 0.0025575 2 42.716962 0.6814586 0.0235444 4 35.900857 0.8425262 0.0034789 N. obliqua 1 50.781108 0.9007745 0.0055009 2 36.532756 0.7602626 0.0036924 3 52.971650 0.7186384 0.0109518 4 60.381728 0.8868993 0.0103325 N. dombeyi 4 42.829543 0.6484692 0.0173443 all 45.342151 0.6340160 0.0239507 1 The characteristics of the model used, the source of information and procedures ar e described by Gezan and Moreno (2000). Curvas de Sitio Altura Dominante para Renovales de Roble, Raul y Coigue. Documento I nterno Proyecto FONDEF D97I1065 Universidad Austral de Chile.
76 APPENDIX B SUPPLEMENTARY TABLES Table B 1 P values with Bonferroni correction obtained from a post hoc test of all six possible pairwise comparisons for growth zones related to two clusters. Zone 1 Zone 2 Zone 3 Zone 2 0.049 Zone 3 0.631 0.999 Zone 4 0.430 0.999 0.999 Table B 2 Mean diameter breast height growth r ates ( standard deviation in parentheses) by specie and growth zone obtained from the available data ( fitting database ). N. alpina N. obliqua N. dombeyi Total Zone 1 2.15 (1.5) 2.92 (2.0) 3.51 (2.3 ) 3.01 (2.1) Zone 2 2.66 (1.6) 3.06 (2.3) 3.07 (2.5) 2.94 (2.1) Zone 3 3.69 (2.2) 4.05 (2.3) 3.82 (2.2) Zone 4 2.71 (1.7) 2.34 (1.5) 3.34 (1.9) 2.75 (1.7) Total 2.58 (1.6) 3.00 (2.1) 3.56 (2.2) 3.07 (2.1) Table B 3 Goodness of fit statistics for factors specie, zone and both in the CV regression model using the test data. n R 2 emp RMSE RMSE% BIAS BIAS% Zone 551 0.55 1.37 44.82 0.07 2.29 Sp 551 0.56 1.36 44.49 0.07 2.29 Zone + Sp 551 0.55 1.38 45.14 0.07 2.29 SpZone 551 0.56 1.36 44.49 0.06 1.96 Note: Sp, specie; SpZone, combination of Specie and Zone; R 2 emp, empirical coefficient of correlation; RMSE, root mean square error; RMSE%, relative root mean square error; BIAS%, relative bias.
77 APPENDIX C SUPPLEMENTARY FIGURES Figure C 1 Scatterplot correlation m atrix of predictors (in original units) for the fitting database.
78 Figure C 2 Scree plot of PCoA with the cumulative variance of each coordinate based on th e continuous variables of the tree and stand subsets. Figure C 3 Scree plot of PCA for continuous variables of tree and stand subset
79 Figure C 4 PC1 and PC3 with scores and factor loadings, using tree and stand subset Figure C 5 PC2 and PC3 with scores and factor loadings, using tree and stand subset
80 Figure C 6 Scree plot of PCA for continuous variables of stand subset and average plot of AIDB H. Figure C 7 Scree plot of optimal number of cluster using within sum of squares an d average silhouette width for s tand level information.
81 Figure C 8 Histograms of expected deltas and observed delta obtaining by MRPP test for two K means cluster using stand level information.
82 LIST OF REFERENCES Anderson MJ, Ellingsen KE, McArdle BH (2006) Multivariate dispersion as a measure of beta diversity Ecology Letters 9:683 693 Andreassen K, Tomter SM (2003) Basal area growth models for individual trees of Norway spruce, Scots pine, birch and other broadleaves in Norway. For Ecol Manage 180:11 24 Avery TE, Burkhart HE (2002) Forest measurements, 5th edn. McGraw Hill, New York Baskerville GL (1972) Use of lo garithmic regression in the estimation of pl ant biomass. Can J For Res 2:49 53 Birdsey R (1990) Updating methods for forest inventories: an overview. In: LaBau VJ, Cunia T (eds) State of the art methodology of forest inventory: a symposium proceedings. US DA For Serv Gen Tech, pp 429 435 Borcard D, Gillet F, Legendre P (2011) Numerical ecology with R. Springer, New York Buckman RE (1962) Growth and yield of red pine in Minnesota. USDA Forest Service, Technical Bulletin 1272, Washington, DC Burkhart H (20 03) Suggestions for choosing an appropriate level for modeling forest stands. In: Amaro A, Reed D, Soares P (eds ) Modelling forest systems. CAB International, Wallingford, pp 3 10 Burkhart HE, Tome M (2012) Modeling Forest Trees and Stands. Springer, Net herlands Chauchard L, Sbrancia R (2003) Diameter growth models for Nothofagus obliqua Bosque 2 4(3):3 16 CONAF (Corporacin Nacional Forestal, CL), CONAMA (Comisin Nacional del Medio Ambiente, CL), BIRF (Banco Internacional de Reconstruccin y Fomento, USA), Universidad Austral de Chile, Pontificia Universidad Catlica de Chile, Universidad Catlica de Temuco (1999) Proyecto catastro y evaluacin de los recursos vegetacionales nativos de Chile. Ministerio de agricultura, Santiago, Chile CONAF (2011) Cat astro de los recursos vegetacionales nativos de Chile. Monitoreo de cambios y actualizaciones. Periodo 1997 2011. Ministerio de agricultura, Santiago, Chile Contreras M, Affleck D, Chung W (2011) Evaluating tree competition indices as predictors of basa l area increment in western Montana forests. For Ecol Manage 262:1939 1949
83 Coomes DA, Allen RB (2007) Effects of size, competition and altitude on tree growth. Journal of Ecology 95:1084 1097 Cubillos V (1987) Modelos de crecimiento diametral para algun os renovales de Raul. Cienc invest for 1(1):67 76. Donoso C (1981) Ecologa Forestal. Editorial Universitaria, Santiago, Chile Donoso C (1987) Variacin natural en especies de Nothofagus en Chile. Bosque 8(2):85 97 Donoso C (1993) Bosques templados de Chile y Argentina. Variacin, estructura y dinmica. EditoriaI Universitaria, Santiago, Chile Donoso C, Premoli A, Gallo L, Ipinza R (2004) Variacin intraespecfica en las especies arbreas de los bosques templados de Chile y Argentina. Editorial Unive rsitaria, Santiago Chile Donoso PJ, Donoso C, Sandoval V (1993) Proposicin de zonas de crecimiento de renovales de roble ( Nothofagus obliqua) y raul ( Nothofagus alpina) en su rango de distribucin natural. Bosque 14( 2):37 55 Donoso PJ, Lusk CH (2007) Differential effects of emergent Nothofagus dombeyi on growth and basal area of canopy species in an old growth temperate rainforest. J V eg Sci 18(5):675 684 Donoso PJ, Soto D (2016) Does site quality affect the additive basal area phenomenon? Results fro m Chilean old growth temperate rai nforests. Can J For Res 00:1330 1336 Echeverria C, Lara A (2004) Growth patterns of secondary Nothofagus obliqua N. alpina forests in souther n Chile. For Ecol Manage 195:29 43 Esse C, Donoso PJ, Gerding V, Encina Montoya F (2013) Determination of homogeneous edaphoclimatic zones for the secondary forests of Nothofagus dombeyi in central southern Chile. Cien Inv Agr 40(2):351 360 Everitt B, Hothorn T (2011) An introduction to applied multivariate analysis with R. Springer New York Fox J, Monette G (1992) Generalized collinearity diagnostics. JASA 87:178 183 Friedman J, Hastie T, Simon N, Tibshirani R (2016) glmnet: Lasso and elastic net regularized generalized linear models. R package version 2.0 5
84 Gezan SA, Moreno P ( 1999). Establecimiento y medicin de una red de parcelas permanentes para renovales de roble ( Nothofagus obliqua), raul ( Nothofagus alpina) y coige ( Nothofagus dombeyi). Bosque Nativo (Valdivia) 22:9 11 Gezan SA, Moreno P, Ortega A (2009) Taper models f or roble, raul and coige second growth forests in Chile. Bosque 30(2):61 69 Gezan SA, Ortega A, Andenmatten E (2007) Diagramas de manejo de densidad para renovales de roble, raul y coige en Chile. Bosq ue 28(2):97 105 Gonzalez Benecke CA, Gezan SA, Le duc DJ, Martin TA, Cropper Jr. WP, Samuelson LJ (2012) Modeling survival, yield, volume partitioning and their response to thinning for longleaf pine plantations. Forests 3(4):1104 1132 Harrison WC, Burk TE, Beck DE (1986) Individual tree basal area incre ment and total height equations for Appalachian mixed hardwoods after thinning. South J Appl For 10:99 104 Hastie T, Tibshirani R, Friedman J (2009) The elements of statistical learning (2nd edition): data mining, inference, and prediction. Springer Verla g, Berln Huang S, Titus SJ (1995) An i ndividual t ree d iameter i ncrement m odel for w hite s pruce in Alberta. Can J For Res 25:1455 1465 Huang S, Titus SJ (1999) Estimating a system of nonlinear simultaneous individual tree models for white spruce in bore al mixed species stands. Can J For Res 29:1805 1811 Huang S, Yang Y, Wang Y (2003) A critical look at procedures for validating growth and yield models. In: Amaro A, Reed D, Soares P (eds ) Modelling forest systems. CAB Int ernational, Wallingford, pp 271 293 INFOR (2013) Rentabilidad econmica de distintas propuestas silvcolas para renovales de Nothofagus en el centro sur de Chile. Informe Tcnico 193. http://biblioteca.infor.cl/DataFiles/3142 9.pdf INFOR (2016) Bosque nativo 12. http://wef.infor.cl/publicaciones/bn/2016/12/BN201612.pdf Accessed 13 March 2017 Ivancich H, Martnez Pastur G, Lencinas MV, Cellini JM, Pe ri PL (2014) Proposals for Nothofagus antarctica diameter growth estimation: simple vs. global models. J For Sci 60:307 317 James G, Witten D, Hastie T, Tibshirani R (2013) An introduction to statistical learning with applications in R. Springer, New York
85 Kuuluvainen T (2002) Natural variability of forests as a reference for restoring and managing biological diversity in boreal Fennoscandia. Silva Fennica 36(1):97 125 Kyung M, Gill J, Ghosh M, Casells G (2010) Penalized regression, standard errors, and B ayesian lassos. Bayesian Anal. 5(2):369 412 Landres PB, Morgan P, Swanson FJ (1999) Overview of the use of natural variability concepts in managing ecological systems. Ecological Applications 9:1179 1188 Lockhart R, Taylor J, Tibshirani RJ, Tibshirani R (2014) A significance test for the lasso. Ann. Stat. 42(2):413 468 Luebert P, Pliscoff P (2006) Sinopsis bioclimtica y vegetacional de Chile. Editorial Universitaria, Santiago, Chile Lumley T, Miller A (2009) Leaps: regression subset selection. R packag e version 2.9. Lusk CH, Ortega A (2003) Vertical structure and basal area development in second growth Nothofagus stands in Chile. J Appl Ecol 50:639 645 Ma W, Lei X (2015) Nonlinear simultaneous equations for individual tree diameter growth and mortali ty model of natural Mongolian oak forests in northeast China. Forests 6:2261 2280 MacDonald PL, Gardner RC (2000) Type I erro r rate comparisons of post hoc procedures for I J chi squa re tables. Educ Psychol Meas 60(5):735 754 Maechler M (2013) Cluster an alysis e xtended Rousseeuw et al. R CRAN McCune B, Grace JB, Urban, DL (2002). Analysis of ecological communities. MjM Software Design, Gleneden Beach Monserud RA, Sterba H (1996) A basal area increment model for individual trees growing in even and unev en aged forest stands in Austria. For Ecol Manage 80:57 80 Moreno P (2001) Proposicin preliminar de curvas de ndice de sitio para renovales de raul. Tesis Ingeniero Forestal. Universidad de Chile Chile Murphy P, Graney D (1998) Individual tree basal area growth, survival, and total height models for upland hardwoods in the Boston Mountains of Arkansas. South J Appl For 22:184 92 Novoa R, Villaseca C ( 1989 ). Mapa agroclim tico de Chile. Instit uto de Investigaciones Agropecuarias (INIA), Ministerio de Agricultura, Santiago Chile
86 Nunes L, Tome J, Tome M (2011) Prediction of annual tree growth and survival for thinned and unthinned even aged maritime pine stands in Portugal from data with different time measurement intervals. For Ecol Manage 262:1491 14 99 Community ecology package Ormazabal C, Benoit I (1987) Present status of conservation of the genus Nothofagus in Chile. Bosque 8(2):109 120 Ortega A, Gezan S (1998) Cuantificacin de crecimiento y proyeccin de calidad en Nothofagus Bosque 19(1):123 126 Peng C (2000) Growth and yield models for uneven aged stands: past, present and future. For Ecol Manage 132(2):259 279 Pukkala T, Lhde E, Laiho O (2009) Growth and yield models for uneven sized forest stands in F inland. For Ecol Manage 258:207 216 Pukkala T, Lhde E, Laiho O (2011) Using optimization for fitting individual tree growth models for uneven aged s tands. Eur J Forest Res 130:829 839 Quinn GP Keough MJ (2002) Experimental design and data analysis for biologists. Cambridge University Press. Cambridge R Core Team (2015) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austri a. URL https://www.R project.org/ Richardson SJ, Hurst JM, Easdale TA, Wiser SK, Griffiths AD, Allen RB (2011) Diameter growth rates of beech ( Nothofagus ) trees around New Zealand. N Z J For 56(1):3 11 Salas C, Garcia O (2006) Modelling height development of mature Nothofagus obliqua. For Ecol Manage 229:1 6 Salas C, Gregoire T, Craven D, Gilabert H (2016) Forest growth modelling: the state of the art. Bosque 37(1):3 12 Schlatter J, Gerding V (1995) A method o f site classification for forestry production: A Chilean example. Bosque 16(2):13 20 Schlatter J, Gerding V, Adriazola J (1994) Sistema de ordenamiento de la tierra. Herramienta para la planificacin forestal, aplicada a las regiones VII, VIII y IX. S eri e tcnica, F acultad C iencias F orestales, Universidad Austral de Chile. Valdivia, Chile
87 Schlatter J, Gerding V, Huber H (1995) Sistema de ordenamiento de la tierra. Herramienta para la planificacin forestal, aplicada para la X regin. Serie tcnica, F ac ultad C iencias F orestales, Universidad Austral de Chile. Valdivia, Chile Shortt JS, Burkhart HE (1996) A comparison of loblolly pine plantation growth and yield models for inventory updating. South J Appl For 20:15 22 Thiers O, Gerding V Hildebrand E ( 2 008 ) Renovales de Nothofagus obliqua en centro y s ur de Chile: Factores de sitio relevantes para su productividad. In : Bava J, Picco OA, Pildan MB, Lpez P, Orellana I ( eds ) Libro de a ctas de Eco Reuniones: segunda r eunin sobre los Nothofagus en la Patag onia. Esquel, Argentina, pp 255 260 Tibshirani R (1996) Regression shrinkage and selection via the Lasso. J R Statist Soc B 58(1):267 288 Ugalde G (2006) Crecimiento en altura de renovales de lenga ( Nothofagus pumilio (Poepp. Et Endl.) Krasser) en Monte Alto (XII regin) en funcin de la calidad de sitio. Tesis Ingeniero Forestal. Universidad de Chile Chile Vanclay JK (1994) Modelling forest growth and yield: applications to mixed tropical forests. CAB International, Wallingford, UK Welham SJ, Gezan SA Clark SJ, Mead A (2015) Statistical methods in biology: design and analysis of experiments and regression. CRC Press, Boca Raton Wu TT, Chen YF, Hastie T, Sobel E, Lange K (2009) Genome wide association analysis by lasso penalized logistic regre ssion. B ioinformatics 25(6):714 721 Wykoff WR (1990) A basal area increment model for individual conifers in the northern Rock y Mountains. For Sci 36(4):1077 1104
88 BIOGRAPHICAL SKETCH Paulo C. Moreno Meynard has experience in a broad range of aspects related to forest resources. He received a title of Forestry Engineer (as US Bachelor) from the Universidad de Chile, Chile, a certificate in Geographic Information Systems and Remote Sen sing in Agroforestry Management from Universidad Catlica de Temuco Chile, a ce rtificate in Geospatial Analysis from University of Florida, and completed his Master of Science in the School of Forest Resources and Conservation at the University of Florida. After his undergraduate experience, Paulo was work ed as forest researcher at U niversidad Austral de Chile, Valdivia, where he started to work with forest inventories, GIS, biometry and statistics focused in the species N othofagus alpina N obliqua and N. dombeyi After this he worked as regional manager at Forest Research Institu te of Chile in Coyhaique, small city located in Chilean Patagonia. Here, Paulo developed research related to management and genetic of N. pumilio b iometry and g enetic of Pinus ponderosa environmental restoration, and CDM projects Protocol In 2009, Paulo was employed as project engineer at Forestal Mininco SA, in Coyhaique, where he contributed to obtain the FSC certification, additionally he worked in production of native species by cuttings management for Pinus ponderosa Pinus contorta, an d Pseudotsuga menziesii biometry, inventories, and GIS. Paulo in 2014 also was a short consultant for World Bank related to mechanisms to reduce carbon emissions, REDD+