Introduction

A preliminary study carried across 10 sites scattered across the Amazon Basin revealed a consistent and uniform relationship between logging intensity (i.e. the percentage of initial biomass loss) and the time to recover intial biomass stocks (Rutishauser et al. 2015). We now aim at understandgin how do timber volume recover after logging and what are the main drivers. Several research questions and a general framework were set up among TmFO partners during a 3 days workshop held in Belem in March 2016.

Research questions and methods

Timber is broadly defined as the volume of trunk of trees >= 50 cm. This volume (cubic meter) was estimated using site-specific volumetric equations. Summing all timber at plot scale provide the timber volume (TV), expressed in m3/ha. TV was further divided into three groups:

The main research questions are:

  1. How long does it take to recover harvested timber volumes (3 levels of analysis, see below)?

    • split into trees DBH > 50 cm present before logging(“old”) and “new” trees DBH > 50 cm recruited post-logging

  2. What are the main drivers of TV recovery post-logging?

Plots selection & descritption

Only plots that were logged (i.e. no control) or were free from major fire disturbances were selected. On 215 plots, 144 were chosen. Harvest intensity ranged from 1.79 and 102.58 m3 ha-1.

Figure 1: TV at each plot and census. Red dot=year of logging & red line = linear regression of which the slope gives the TV recovery rate.

Defining logging intensity & TV recovery

A straightforward definition of logging intensity is to sum all trees harvested within a plot. Plots that experienced post-logging treatments at Paracou, Tapajos and Ecosilva were discarded. For sake of comparison, logging intensity is expressed as a fraction of initial timber volume present before logging (%). Timbre volume (TV) recovery rates (m3 ha-1 y-1) are further computed as linear functions between TV stocks and time since logging.

Figure 2: Volume harvested by initial timber volume per plot

Table 1: Logging information by plot
plot-ID Logging date TV harvested Initial TV % TV harv. Pre-logging stem density Site
chi-1 2004 27.452179 111.41369 0.2463986 30 ChicoBocao
chi-2 2004 2.983902 75.63877 0.0394494 25 ChicoBocao
chi-3 2004 3.106371 99.07330 0.0313543 34 ChicoBocao
chi-5 2004 18.852582 128.78453 0.1463886 32 ChicoBocao
chi-7 2004 25.717796 70.10773 0.3668325 18 ChicoBocao
cum-2 2007 3.944873 72.65489 0.0542960 24 Cumaru
cum-3 2007 64.944506 166.76745 0.3894315 33 Cumaru
cum-4 2007 22.315112 103.22454 0.2161803 21 Cumaru
cum-6 2007 9.792339 107.24055 0.0913119 28 Cumaru
eco-X1 2005 36.148768 96.63825 0.3740627 26 Ecosilva
eco-X10 2005 102.580911 192.28002 0.5334975 35 Ecosilva
eco-X11 2005 19.563292 130.28788 0.1501544 35 Ecosilva
eco-X12 2005 35.993045 84.53886 0.4257574 23 Ecosilva
eco-X13 2005 20.454547 84.05106 0.2433586 22 Ecosilva
eco-X14 2005 35.554555 127.57503 0.2786952 40 Ecosilva
eco-X15 2005 31.576301 122.10023 0.2586097 28 Ecosilva
eco-X16 2005 25.509817 130.02637 0.1961896 34 Ecosilva
eco-X17 2005 50.843147 118.15151 0.4303216 37 Ecosilva
eco-X18 2005 4.974913 127.86599 0.0389072 29 Ecosilva
eco-X2 2005 16.329594 88.17447 0.1851964 20 Ecosilva
eco-X3 2005 44.586666 115.71796 0.3853046 32 Ecosilva
eco-X4 2005 31.549139 103.53018 0.3047337 29 Ecosilva
eco-X5 2005 42.594313 116.00002 0.3671923 32 Ecosilva
eco-X6 2005 21.077340 101.11611 0.2084469 28 Ecosilva
eco-X7 2005 33.040448 88.54397 0.3731530 26 Ecosilva
eco-X8 2005 9.963002 56.98787 0.1748267 18 Ecosilva
eco-X9 2005 30.893274 130.59006 0.2365668 35 Ecosilva
ira-10 2004 28.936614 51.38456 0.5631383 12 Iracema
ira-11 2004 6.013837 37.93177 0.1585435 14 Iracema
ira-13 2004 12.229531 61.07040 0.2002530 18 Iracema
ira-14 2004 10.649455 51.44328 0.2070135 15 Iracema
ira-15 2004 8.430066 46.65775 0.1806788 16 Iracema
ira-16 2004 11.268172 59.73868 0.1886244 16 Iracema
ira-17 2004 16.341662 58.47762 0.2794516 18 Iracema
ira-18 2004 4.104867 71.82236 0.0571530 22 Iracema
ira-19 2004 3.523023 93.41676 0.0377130 26 Iracema
ira-2 2004 6.015316 60.50047 0.0994259 17 Iracema
ira-20 2004 7.659879 39.05037 0.1961538 12 Iracema
ira-22 2004 6.162740 103.89010 0.0593198 25 Iracema
ira-25 2004 10.630246 43.96495 0.2417891 12 Iracema
ira-26 2004 4.656273 65.86486 0.0706943 21 Iracema
ira-27 2004 11.130617 80.10165 0.1389561 19 Iracema
ira-29 2004 5.694724 65.03477 0.0875643 18 Iracema
ira-3 2004 5.340486 62.06409 0.0860479 17 Iracema
ira-30 2004 10.271264 71.43515 0.1437845 20 Iracema
ira-4 2004 9.299168 77.88491 0.1193963 16 Iracema
ira-5 2004 19.461608 114.13859 0.1705086 24 Iracema
ira-6 2004 6.420440 72.65411 0.0883700 16 Iracema
ira-7 2004 3.024427 51.22692 0.0590398 17 Iracema
itab-1 1996 15.195549 106.98033 0.1420406 38 Itacoatiara B
itab-10 1996 23.531854 96.26492 0.2444489 35 Itacoatiara B
itab-11 1996 30.867834 96.02404 0.3214594 33 Itacoatiara B
itab-12 1996 34.592513 104.53683 0.3309122 35 Itacoatiara B
itab-13 1996 44.692951 148.40279 0.3011598 41 Itacoatiara B
itab-14 1996 4.467986 80.63414 0.0554106 22 Itacoatiara B
itab-2 1996 31.052926 133.28167 0.2329872 45 Itacoatiara B
itab-3 1996 20.819199 92.70594 0.2245724 26 Itacoatiara B
itab-5 1996 34.493212 122.39915 0.2818092 31 Itacoatiara B
itab-6 1996 34.329086 94.16960 0.3645453 32 Itacoatiara B
itab-7 1996 5.593218 106.00656 0.0527629 33 Itacoatiara B
itab-8 1996 33.359417 104.39135 0.3195611 38 Itacoatiara B
itab-9 1996 5.266368 46.27985 0.1137940 17 Itacoatiara B
itac-1 1997 13.400004 82.10561 0.1632045 28 Itacoatiara C
itac-10 1997 8.611107 59.88961 0.1437830 20 Itacoatiara C
itac-11 1997 2.785741 80.25420 0.0347115 24 Itacoatiara C
itac-13 1997 6.846188 66.05611 0.1036420 12 Itacoatiara C
itac-2 1997 10.193673 65.57412 0.1554527 21 Itacoatiara C
itac-4 1997 14.197146 88.92517 0.1596527 22 Itacoatiara C
itac-6 1997 14.475289 131.08110 0.1104300 32 Itacoatiara C
itac-7 1997 1.976824 70.40800 0.0280767 24 Itacoatiara C
itac-8 1997 32.125217 111.15498 0.2890129 32 Itacoatiara C
itac-9 1997 8.252058 27.96598 0.2950749 13 Itacoatiara C
itad-1 1998 21.173490 127.58390 0.1659574 39 Itacoatiara D
itad-11 1998 4.620159 139.41658 0.0331392 35 Itacoatiara D
itad-12 1998 30.779182 114.44401 0.2689453 36 Itacoatiara D
itad-13 1998 8.117371 87.49357 0.0927768 28 Itacoatiara D
itad-14 1998 43.444763 176.94851 0.2455221 50 Itacoatiara D
itad-2 1998 8.924915 133.57555 0.0668155 38 Itacoatiara D
itad-3 1998 15.236311 120.49003 0.1264529 39 Itacoatiara D
itad-4 1998 34.224032 143.27283 0.2388731 36 Itacoatiara D
itad-5 1998 25.555566 152.81034 0.1672372 40 Itacoatiara D
itad-6 1998 11.599127 119.29894 0.0972274 37 Itacoatiara D
itad-7 1998 23.267348 167.74809 0.1387041 38 Itacoatiara D
itad-8 1998 7.597272 105.23348 0.0721944 27 Itacoatiara D
itad-9 1998 14.748177 101.06605 0.1459261 32 Itacoatiara D
jar-101 1986 28.006961 139.44948 0.2008395 32 Jari
jar-102 1986 69.673088 168.34552 0.4138696 36 Jari
jar-103 1986 15.050889 100.46717 0.1498090 32 Jari
jar-104 1986 31.542695 110.16755 0.2863157 30 Jari
jar-105 1986 37.880698 99.99041 0.3788433 30 Jari
jar-106 1986 13.765716 113.25966 0.1215412 29 Jari
jar-107 1986 18.972694 92.18202 0.2058177 28 Jari
jar-108 1986 26.817577 108.20638 0.2478373 34 Jari
jar-109 1986 5.729426 108.40359 0.0528527 29 Jari
jar-110 1986 25.497145 151.32935 0.1684878 46 Jari
jar-111 1986 27.336089 115.56600 0.2365409 35 Jari
jar-112 1986 20.147254 148.84919 0.1353535 29 Jari
jar-201 1986 9.886241 142.54580 0.0693548 45 Jari
jar-202 1986 33.850733 113.39700 0.2985152 33 Jari
jar-203 1986 19.319679 126.29126 0.1529772 41 Jari
jar-204 1986 25.729300 138.55720 0.1856944 34 Jari
jar-205 1986 38.942339 152.87001 0.2547415 35 Jari
jar-206 1986 7.779658 79.73304 0.0975713 28 Jari
jar-207 1986 21.387511 75.89608 0.2817999 21 Jari
jar-208 1986 35.918791 120.03841 0.2992275 39 Jari
jar-209 1986 45.005265 138.94824 0.3238995 34 Jari
jar-210 1986 10.810684 141.91705 0.0761761 42 Jari
jar-211 1986 14.276482 93.66656 0.1524181 29 Jari
jar-212 1986 36.494251 108.35283 0.3368094 31 Jari
jar-301 1986 30.243143 111.40384 0.2714731 35 Jari
jar-302 1986 26.665473 91.86113 0.2902803 30 Jari
jar-303 1986 25.673819 74.70678 0.3436612 22 Jari
jar-304 1986 41.564297 99.58306 0.4173832 25 Jari
jar-305 1986 8.975882 92.60434 0.0969272 27 Jari
jar-306 1986 12.502798 122.32393 0.1022106 37 Jari
jar-307 1986 9.696712 96.80833 0.1001640 27 Jari
jar-308 1986 22.369825 116.27704 0.1923839 35 Jari
jar-309 1986 22.001872 103.65395 0.2122627 29 Jari
jar-310 1986 40.980410 111.24313 0.3683860 32 Jari
jar-311 1986 21.680037 114.73715 0.1889539 34 Jari
jar-312 1986 25.373064 90.76530 0.2795459 27 Jari
pag-CL 1993 14.687860 50.91820 0.2884600 362 Paragominas
pag-RIL 1993 12.337652 50.00440 0.2467313 346 Paragominas
par-2 1987 31.150818 83.23657 0.3742444 175 Paracou
par-7 1987 27.064622 80.52307 0.3361102 183 Paracou
par-9 1987 26.305911 65.21304 0.4033842 143 Paracou
pet-4 2003 18.023931 83.16657 0.2167209 26 Peteco
pet-5 2003 18.152080 125.79247 0.1443018 34 Peteco
pet-6 2003 45.052938 103.45480 0.4354843 27 Peteco
pet-7 2003 22.572429 80.08924 0.2818410 20 Peteco
pet-9 2003 4.782163 116.77622 0.0409515 31 Peteco
tab-1 2001 14.083986 66.60211 0.2114646 21 Tabocal
tab-2 2001 4.886323 58.89183 0.0829712 22 Tabocal
tab-3 2001 3.793691 52.76877 0.0718927 15 Tabocal
tab-4 2001 1.788432 37.68834 0.0474532 16 Tabocal
tab-5 2001 10.018218 38.97698 0.2570291 9 Tabocal
tab-6 2001 5.925447 42.93508 0.1380095 12 Tabocal
tap-2 1983 31.986567 80.49897 0.3973537 24 Tapajos
tap-3 1983 26.451640 124.80372 0.2119459 40 Tapajos
tap-4 1983 43.110006 153.39405 0.2810409 34 Tapajos
tap-5 1983 44.636167 119.47554 0.3736009 35 Tapajos
tap-6 1983 27.160839 122.27577 0.2221277 36 Tapajos
tap-7 1983 15.870821 52.17835 0.3041649 17 Tapajos
tap-8 1983 56.564566 106.45033 0.5313705 29 Tapajos

Regional list of commercial species

All site leaders were asked to provide a list of the main commercial species (compiled below). In total, there is 265 species considered.

Table 2: Number of commercial species and proportion of timber tree population per site
label Nsp prop
ChicoBocao 9 12.9
Cumaru 9 15.3
Ecosilva 28 30.4
Iracema 11 10.3
Itacoatiara B 31 32.3
Itacoatiara C 14 18.7
Itacoatiara D 32 32.7
Jari 25 15.6
Paragominas 28 24.8
Paracou 33 31.4
Peteco 15 19.5
Tabocal 4 8.0
Tapajos 20 22.2 
## Warning in matrix(unique(list$nameV[order(list$nameV)]), ncol = 3): data
## length [265] is not a sub-multiple or multiple of the number of rows [89]
Table 3: List of commercial species across all sites (regional list)
Scientific name Scientific name Scientific name
Agonandra brasiliensis Enterolobium schomburkii Ocotea tomentella
Albizia pedicellaris Eperua falcata Ormosia coutinhoi
Alchornea triplinervia Eperua grandiflora Paramachaerium ormosioides
Alexa grandiflora Eperua grandifolia Parinari campestris
Amanoa guianensis Eperua rubiginosa Parkia decussata
Anacardium giganteum Eriotheca longipedicellata Parkia gigantocarpa
Anacardium parvifolium Erisma uncinatum Parkia multijuga
Andira laurifolia Eschweilera coriacea Parkia paraensis
Andira parviflora Eschweilera odorata Parkia pendula
Andira unifolialata Eschweilera ovata Parkia ulei
Aniba burchellii Eschweilera paniculata Parkia velutina
Aniba canelilla Eschweilera sagotiana Peltogyne catingae
Antonia ovata Euxylophora paraensis Persea laevigata
Apeiba albiflora Ficus nymphaeifolia Phyllocarpus riedellii
Apeiba petoumo Ficus piresiana Piptadenia suaveolens
Apuleia leiocarpa Glycydendron amazonicum Pithecellobium incuriale
Apuleia leocarpa Goupia glabra Pithecellobium racemosum
Apuleia molaris Guazuma ulmifolia Platonia insignis
Astronium gracile Handroanthus impetiginosus Platymiscium filipes
Astronium lecointei Handroanthus serratifolia Poelcilanthe effusa
Astronium leicoitei Handroanthus serratifolius Pouteria anomala
Bagassa guianensis Humiria balsamifera Pouteria bilocularis
Balizia pedicellaris Hymenaea courbaril Pouteria engleri
Bertholletia excelsa Hymenaea oblongifolia Pouteria eugeniifolia
Bocoa prouacensis Hymenaea parvifolia Pouteria flavilatex
Bombacopsis nervosa Hymeneaea courbaril Pouteria guianensis
Bombax globosum Hymeneaea oblongifolia Pouteria laevigata
Bowdichia nitida Hymenolobium excelsum Pouteria opposita
Brosimum acutifolium Hymenolobium heterocarpum Pouteria oppositifolia
Brosimum guianense Hymenolobium modestum Pouteria platyphylla
Brosimum lactescens Hymenolobium nitidum Pouteria rodriguesiana
Brosimum lanciferum Hymenolobium petraeum Pradosia cochlearia
Brosimum parinarioides Hymenolobium pulcherrimum Protium paniculatum
Brosimum potabile Hymenolobium sericeum Protium puncticulatum
Brosimum rubescens Ilex inundata Pseudopiptadenia psilostachya
Brosimum uleanum Inga alba Pseudopiptadenia suaveolens
Brosimum utile Iryanthera crassifolia Pterocarpus officinalis
Buchenavia capitata Iryanthera grandis Qualea paraensis
Buchenavia guianensis Iryanthera paraensis Qualea rosea
Buchenavia nitidissima Jacaranda copaia Qualea tesmannii
Buchenavia parvifolia Laetia procera Roupala montana
Buchenavia viridiflora Lecythis lurida Ruizterania albiflora
Byrsonima laevigata Lecythis pisonis Sacoglottis guianensis
Calophyllum brasiliense Lecythis poiteaui Schefflera morototoni
Candolleodendron brachystachyum Lecythis prancei Scleronema micranthum
Capirona huberiana Lecythis zabucajo Sextonia rubra
Carapa guianensis Licania alba Simarouba amara
Carapa guiansensis Licania ovalifolia Sterculia pruriens
Carapa surinamensis Licaria aritu Stryphnodendron adstringens
Cariniana micrantha Licaria brasiliensis Stryphnodendron pulcherrimum
Caryocar glabrum Licaria cannella Swartzia grandifolia
Caryocar microcarpum Licaria crassifolia Swartzia panacoco
Caryocar villosum Licaria rigida Swartzia polyphylla
Cedrela fissilis Luehea grandiflora Symphonia globulifera
Cedrela odorata Lueheopsis rugosa Tabebuia impetiginosa
Ceiba pentandra Manilkara bidentata Tabebuia insignis
Chaetocarpus schomburgkianus Manilkara cavalcantei Tabebuia serratifolia
Chimarrhis turbinata Manilkara huberi Tachigali goeldianum
Chrysophyllum lucentifolium Manilkara paraensis Tachigali melinonii
Chrysophyllum pachycarpa Manilkara surinamensis Tachigali myrmecophyla
Chrysophyllum pomiferum Mezilaurus duckei Tachigali myrmecophylla
Chrysophyllum prieurii Mezilaurus itauba Tachigali paraensis
Chrysophyllum sanguinolentum Mezilaurus sinandra Tapura capitulifera
Chrysophyllum venezuelanense Mezilaurus synandra Tetragastris altissima
Clarisia racemosa Micrandropsis scleroxylon Torresea acreana
Copaifera multijuga Micropholis egensis Trattinnickia rhoifolia
Cordia bicolor Micropholis guyanensis Vantanea guianensis
Cordia goeldiana Micropholis melioniana Vantanea parviflora
Couma guianensis Micropholis venulosa Vatairea eritrocarpa
Couratari guianensis Minquartia guianensis Vatairea guianensis
Couratari multiflora Monopteryx inpae Vatairea paraensis
Couratari oblongifolia Moronobea coccinea Vatairea sericea
Couratari stellata Myroxylon balsamo Virola duckei
Dendrobangia boliviana Nectandra cissiflora Virola kwatae
Dialium guianense Nectandra micranthera Virola michelii
Diatenopteryx sorbifolia Nectandra purusensis Virola michellii
Dicorynia guianensis Ocotea acutangula Virola mickelii
Dinizia excelsa Ocotea caudata Virola multinervia
Diplotropis purpurea Ocotea cernua Virola surinamensis
Diplotropis racemosa Ocotea costulata Vochysia guianensis
Diplotropis triloba Ocotea fragantissima Vochysia maxima
Dipteryx odorata Ocotea fragrantissima Vochysia neyratii
Dipteryx punctata Ocotea glomerata Vochysia surinamensis
Ecclinusa guianensis Ocotea longifolia Vochysia tomentosa
Endlicheria bracteata Ocotea neesiana Vouacapoua americana
Endopleura uchi Ocotea opifera Zollernia paraensis
Enterolobium maximum Ocotea petalanthera Zygia racemosa
Enterolobium oldemanii Ocotea puberula Agonandra brasiliensis
Enterolobium schomburgkii Ocotea rubra Albizia pedicellaris

Exploratory analysis

Relationship between timber recovery and logging intensity

Here, we check whether recovery rates (Figure 3) or the amount of timber cummulated over the study period (Figure 4) are somehow related to logging intensity.

Figure 3: TV recovery rates vs logging intenstiy (% of initial volume removed). A. site-specific list, B. regional list, C. all timber.

Figure 4: Cumulated TV vs logging intenstiy (% of initial volume removed) for site-specific (A), regional (B) and all (C) timber species.

Time to recover TV harvested

Recovery times (year) are obtained by dividing the volume of timber harvested in each plot by the post-logging TV recovery rates. Effects of different biometric response variables on recovery times were further explored through linear mixed effects models.

Figure 5: Range of TV recovery rates accounting for site-specific (A), regional (B) and all (C) timber species

As for the biomass study, we used linear mixed-effects models, where sites are regarded as “random” effects. Plots were weighted in function of their size and lenght between initial and final censuses. The best models are found through an exhaustive screening and ranking using Bayesian Information Criterion (BIC) (package lmerTest). To reduce residual heteroscedasticity, recovery time was log-transformed.

Bioclimatic variables were extracted from WorldClim: Mean Temperature of Warmest Quarter(TWQ), Mean Temperature of Coldest Quarter (TCQ), Annual Precipitation (Precip),Precipitation of Wettest Month (PWeM), Precipitation of Driest Month (PDM), Pre-cipitation Seasonality (Coefficient of Variation) (Season), Precipitation of Wettest Quarter (PWeQ), Precipitation of Driest Quarter (PDQ). Soil properties were extracted from the Harmonized World Soil raster at a resolution of 30 arc-seconds. Information on top soil (0-30 cm) was extracted at each site: texture, drainage, available water content (range), clay,silt and sand content (%), cation-exchange capacity (CEC,cmol.kg-1) and bulk density (bulk.dens,kg.dm-3).

Clay content was well correlated (cor = -0.98) with the firs PCA axis and bulk density with the second PCA axis (cor = 0.81). Only both variables were further used as soil proxy. Precipitation of the wettest quartile (PWQ) was well correlated (cor = 0.98) with the firs PCA axis and Precipitation of the driest quartile (PDQ) with the second PCA axis (cor = 0.88), and both were used as bioclimatic proxy. Initial forest structures, such as initial TV and forest development stage (i.e. difference between plot-TV and site averaged TV), and logging intensity were included, as explanatory variables.

A. Site-specific list of commercial species

Here, only the commercial species logged at a site are accounted for. Recovery rates and times depend on a small pool of species.

# Variable selection recovery time
require(glmulti)
require(lme4)
require(MuMIn)

GLM1 <- glmulti(log(rec1) ~ ratioV  + vol.ini  + vol.rel    + bulk.dens   + clay  + PWQ + PDQ -1 , data=REL[!is.na(REL$rec1),],  random="+ (1 | site)", fitfunc =lmer.glmulti,report=F, level = 1,crit="bic",method="g",plotty=F) 
## TASK: Genetic algorithm in the candidate set.
## Initialization...
## Algorithm started...
## Improvements in best and average IC have bebingo en below the specified goals.
## Algorithm is declared to have converged.
## Completed.
print(GLM1)
## glmulti.analysis
## Method: g / Fitting: lmer.glmulti / IC used: bic
## Level: 1 / Marginality: FALSE
## From 100 models:
## Best IC: 238.710602396773
## Best model:
## [1] "log(rec1) ~ 1 + ratioV"
## Evidence weight: 0.301462205597357
## Worst IC: 267.132871897282
## 2 models within 2 IC units.
## 23 models to reach 95% of evidence weight.
## Convergence after 200 generations.
## Time elapsed: 7.40018391609192 minutes.
tmp<- weightable(GLM1)
tmp <- tmp[tmp$bic <= min(tmp$bic) + 5,]
tmp
##                                          model      bic    weights
## 1                       log(rec1) ~ 1 + ratioV 238.7106 0.30146221
## 2             log(rec1) ~ 1 + ratioV + vol.ini 239.6410 0.18931703
## 3             log(rec1) ~ 1 + ratioV + vol.rel 241.1553 0.08879106
## 4                 log(rec1) ~ 1 + ratioV + PDQ 242.6081 0.04294442
## 5 log(rec1) ~ 1 + ratioV + vol.ini + bulk.dens 242.6410 0.04224382
## 6           log(rec1) ~ 1 + ratioV + bulk.dens 242.8867 0.03736061
## 7                 log(rec1) ~ 1 + ratioV + PWQ 243.0724 0.03404724
## 8                log(rec1) ~ 1 + ratioV + clay 243.0800 0.03391878
## 9       log(rec1) ~ 1 + ratioV + vol.ini + PDQ 243.6436 0.02558865
plot(GLM1,type="s")

summary(GLM1@objects[[1]])
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(rec1) ~ 1 + ratioV + (1 | site)
##    Data: data
## Weights: dat$wght
## 
##      AIC      BIC   logLik deviance df.resid 
##    229.2    238.7   -110.6    221.2       76 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3761 -0.5951 -0.0208  0.5536  3.4694 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.2040   0.4517  
##  Residual             0.8239   0.9077  
## Number of obs: 80, groups:  site, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   3.6731     0.2950  12.452
## ratioV        4.9019     0.9622   5.094
## 
## Correlation of Fixed Effects:
##        (Intr)
## ratioV -0.769
r.squaredGLMM(GLM1@objects[[1]])
##       R2m       R2c 
## 0.2471778 0.3966079
Figure 6: Recovery time (log scale) of timber volume of site-specific timber species in function of harvested TV (% of initial TV).

B. Regional list of commercial species

Here, timber species harvested across all study sites were used to compute TV.

# Variable selection recovery time
GLM2 <- glmulti(log(rec2) ~ ratioV  + vol.ini  + vol.rel    + bulk.dens   + clay  + PWQ + PDQ -1 , data=REL[!is.na(REL$rec2),],  random="+ (1 | site)", fitfunc =lmer.glmulti,report=F, level = 1,crit="bic",method="g",plotty=F)
## TASK: Genetic algorithm in the candidate set.
## Initialization...
## Algorithm started...
## Improvements in best and average IC have bebingo en below the specified goals.
## Algorithm is declared to have converged.
## Completed.
print(GLM2)
## glmulti.analysis
## Method: g / Fitting: lmer.glmulti / IC used: bic
## Level: 1 / Marginality: FALSE
## From 100 models:
## Best IC: 228.765404568569
## Best model:
## [1] "log(rec2) ~ 1 + ratioV + vol.ini"
## Evidence weight: 0.236906670596473
## Worst IC: 246.785994052981
## 2 models within 2 IC units.
## 34 models to reach 95% of evidence weight.
## Convergence after 210 generations.
## Time elapsed: 8.64338612556458 minutes.
tmp<- weightable(GLM2)
tmp <- tmp[tmp$bic <= min(tmp$bic) + 5,]
tmp
##                                                  model      bic    weights
## 1                     log(rec2) ~ 1 + ratioV + vol.ini 228.7654 0.23690667
## 2  log(rec2) ~ 1 + ratioV + vol.ini + clay + PWQ + PDQ 230.6550 0.09209887
## 3        log(rec2) ~ 1 + ratioV + vol.ini + clay + PWQ 231.0558 0.07537619
## 4                     log(rec2) ~ 1 + ratioV + vol.rel 231.7110 0.05431875
## 5        log(rec2) ~ 1 + ratioV + vol.rel + clay + PWQ 231.7704 0.05273020
## 6               log(rec2) ~ 1 + ratioV + vol.ini + PWQ 232.1905 0.04274027
## 7              log(rec2) ~ 1 + ratioV + vol.ini + clay 232.2192 0.04213036
## 8               log(rec2) ~ 1 + ratioV + vol.ini + PDQ 232.6493 0.03397819
## 9           log(rec2) ~ 1 + ratioV + vol.ini + vol.rel 233.0768 0.02743900
## 10 log(rec2) ~ 1 + ratioV + vol.rel + clay + PWQ + PDQ 233.0779 0.02742364
## 11        log(rec2) ~ 1 + ratioV + vol.ini + bulk.dens 233.0901 0.02725665
## 12             log(rec2) ~ 1 + ratioV + vol.rel + clay 233.1521 0.02642525
## 13                              log(rec2) ~ 1 + ratioV 233.4267 0.02303471
plot(GLM2,type="s")

summary(GLM2@objects[[1]])
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(rec2) ~ 1 + ratioV + vol.ini + (1 | site)
##    Data: data
## Weights: dat$wght
## 
##      AIC      BIC   logLik deviance df.resid 
##    217.0    228.8   -103.5    207.0       72 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4479 -0.6432 -0.2955  0.6052  3.1554 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.06645  0.2578  
##  Residual             0.81332  0.9018  
## Number of obs: 77, groups:  site, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 1.850072   0.430916   4.293
## ratioV      3.356719   0.894955   3.751
## vol.ini     0.011789   0.003809   3.096
## 
## Correlation of Fixed Effects:
##         (Intr) ratioV
## ratioV  -0.424       
## vol.ini -0.805 -0.085
r.squaredGLMM(GLM2@objects[[1]])
##       R2m       R2c 
## 0.2678426 0.3231444

Figure 7: Recovery time (log scale) of timber volume of regional timber species in function of harvested TV (% of initial TV).

C. All timber species

Here, all trees DBH > 50 cm were lumped together to compute TV.

## TASK: Genetic algorithm in the candidate set.
## Initialization...
## Algorithm started...
## Improvements in best and average IC have bebingo en below the specified goals.
## Algorithm is declared to have converged.
## Completed.
## glmulti.analysis
## Method: g / Fitting: lmer.glmulti / IC used: bic
## Level: 1 / Marginality: FALSE
## From 100 models:
## Best IC: 219.746546823389
## Best model:
## [1] "log(rec3) ~ 1 + ratioV + vol.ini"
## Evidence weight: 0.213971699976682
## Worst IC: 274.12893462338
## 3 models within 2 IC units.
## 34 models to reach 95% of evidence weight.
## Convergence after 250 generations.
## Time elapsed: 8.28837203979492 minutes.
##                                                  model      bic    weights
## 1                     log(rec3) ~ 1 + ratioV + vol.ini 219.7465 0.21397170
## 2              log(rec3) ~ 1 + ratioV + vol.ini + clay 220.8099 0.12573473
## 3        log(rec3) ~ 1 + ratioV + vol.ini + clay + PWQ 221.0795 0.10987696
## 4                        log(rec3) ~ 1 + ratioV + clay 222.7325 0.04807929
## 5                  log(rec3) ~ 1 + ratioV + clay + PWQ 222.9611 0.04288783
## 6                               log(rec3) ~ 1 + ratioV 223.1801 0.03843829
## 7           log(rec3) ~ 1 + ratioV + vol.ini + vol.rel 223.5488 0.03196685
## 8               log(rec3) ~ 1 + ratioV + vol.ini + PDQ 223.6631 0.03019231
## 9               log(rec3) ~ 1 + ratioV + vol.ini + PWQ 223.9445 0.02622927
## 10        log(rec3) ~ 1 + ratioV + vol.ini + bulk.dens 224.0114 0.02536576
## 11             log(rec3) ~ 1 + ratioV + vol.rel + clay 224.0889 0.02440180
## 12 log(rec3) ~ 1 + ratioV + vol.ini + bulk.dens + clay 224.4065 0.02081890
## 13       log(rec3) ~ 1 + ratioV + vol.ini + clay + PDQ 224.4278 0.02059829
## 14       log(rec3) ~ 1 + ratioV + vol.rel + clay + PWQ 224.4970 0.01989806

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(rec3) ~ 1 + ratioV + vol.ini + (1 | site)
##    Data: data
## Weights: dat$wght
## 
##      AIC      BIC   logLik deviance df.resid 
##    208.2    219.7    -99.1    198.2       69 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4558 -0.6995 -0.0726  0.4825  2.2454 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.02137  0.1462  
##  Residual             0.83386  0.9132  
## Number of obs: 74, groups:  site, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.256994   0.424342   0.606
## ratioV      8.121137   0.930693   8.726
## vol.ini     0.011520   0.003881   2.968
## 
## Correlation of Fixed Effects:
##         (Intr) ratioV
## ratioV  -0.446       
## vol.ini -0.810 -0.082
##       R2m       R2c 
## 0.5621281 0.5730705

Figure 8: Recovery time (log scale) of timber volume of all species in function of harvested TV (% of initial TV).