In sylvainschmitt/Paracou: Paracou data set analysis
Abstract
Keywords
Beta-diversity, Selective logging, Tropical evergreen forest
Introduction
Tropical logging accounts for one eighth of global timber production [@blaser11]. It represents a major threat to biodiversity [@gibson11]. Most tropical timber production originates from selective logging, the targeted harvesting of timber from commercial species in a single cutting cycle [@martin15]. Others sylvicultural practices, such as tree removal by poisoning, are sometimes also involved in selective logging.
The short-term negative impact of selective logging on biodiversity is well known [@deAvila15; @carreno12; @gibson11; @martin15]. @deAvila15 found that additional sylvicultural practices and logging damages had stronger negative effect on diversity than logging itself. @carreno12 found no direct effect of medium intensity logging on species diversity but a community drift towards 'fast' [@reich14] acquisitive species. However, even the short term impact on biodiversity is unclear, as sometimes short term increases in diversity can occur at low sylvicultural intensities [@martin15]. Much less is known about the long term impact [@osazuwa15a; @osazuwa15b].
There are a number of ways diversity can be measured. The first ways is to count the number of species, known as species richness. This measurement highlight the richness but not the evenness of the abundance distribution of species within communities. Second, the Shannon and Simpson indices takes more evenness into account, particularly the latter [@shannon49]. Finally, @tsallis88 combined these measurements with the help of entropy, i.e. the amount of uncertainty calculated from the frequency distribution of a community [@marcon14b]. @jost07 helped to clarify the concept of partitioning diversity $\gamma$, the total diversity of the metacommunity into $\beta$ and $\alpha$ diversities. But the debate is ongoing, in spite of an attempt to resolve it by @chao12.
$\alpha$-diversity, the local diversity, has been highly documented in tropical forests [@valencia94]. Whereas $\beta$-diversity, how species composition changes with distances and among communities, has been less investigated [@condit02]. Still, $\beta$-diversity is a central concept about diversity organisation inside tropical forests. $\beta$-diversity can result from species adaptations to climates or substrates [@paoli06], from limited dispersal [@qian09], or from other historical effects such as logging. Consequently $\beta$-diversity seems at least as important as $\alpha$-diversity for conservation and forest management.
In this paper we designed communities from individual cohorts with similar history. And we used $\beta$-diversity indices between those communities to study and compare turnover of logged and unlogged forests along time.
Material and methods
Study site
The study was carried out in Paracou research station ($5^\circ 18'N$ and $52^\circ 53'W$) in French Guiana. The site has a tropical wet climate, with an average annual rainfall of 2920mm and temperature of $26^\circ C$ marked with a 4 month dry season and a 8 month rainy season interrupted by a short drier period. The forest is a lowland tropical rain forest with a dominance of Fabaceae, Chrysobalanaceae, Lecythidaceae and Sapotaceae.
In 1984, twelves plots of 6.25ha were established to study the effects of different intensities of selective logging. From 1986 to 1988, four logging treatments (T0, T1, T2 and T3) were applied to the plots according to randomized block design with three replicates [@baraloto12]:
- T0 plots are the control plots were no selective logging or other sylvicultural practices were applied.
- T1 plots underwent selective felling for timber of 50-60cm diameter at breast height (dbh) trees from 58 commercial species at an intensity of about 10 trees/ha. Selective logging in T1 plots resulted in 12-33% above ground biomass (AGB) loss and $41.7 \pm 4.0 ~ ( mean \pm standard~deviation)$ species logged.
- T2 plots underwent selective felling for timber of about 30 trees/ha from same commercial species. Additionally, non-commercial species (NCS) over 50cm dbh and commercial species with major defects were eliminated by poison-girdling. Selective logging in T2 plots resulted in 33-56% AGB loss and $41.7 \pm 6.7$ species logged.
- T3 plots underwent selective felling for timber from same commercial species plus logging of NCS of 40-50cm dbh for fuelwood. Additionally, NCS over 50cm dbh and commercial species with major defects were eliminated by poison-girdling. Selective logging in T3 plots resulted in 35-56% AGB loss and $53.0 \pm 3.0$ species logged.
Floral composition
Before logging, each tree was identified by its vernacular name. Vernacular names can encompass different species and sometimes different genera. Since 1990, survivors and recruited trees have begun to be identified by species name. Consequently, trees removed in 1987 or trees that died before they could be categorized were not identified to the species level. But use of vernacular names has been shown to only induce a small discrepancy (Mirabel, personal communications). Our analysis will then use vernacular names. Censuses included 263 vernacular names.
Analysis
Time sampling
Inventories of all plots were conducted annually from 1984 to 1994, every two years from 1995 to 2012, and annually from 2012 on trees above 10 cm gbh. For our analysis we generate recruit figures for missing years by replicating that of the previous year (n-1) in order to simulate turnover measurements.
Community design
To study floral turnovers we divided the plots into communities with similar histories. We used the logging year as the start date for both logged and unlogged plots. Trees which attain 10cm dbh after this date were defined as recruits (R) whereas trees which had 10cm dbh before this date were defined as survivors (S), both in logged and unlogged plots. In each census, recruits and survivors were divided in two communities: alive trees (A) and dead trees (D). We thus obtain 4 communities: Alive Recruits (AR), Alive Survivors (AS), Dead Recruits (DR), and Dead Survivors (DS) for each census in both logged and unlogged plots (see $figure~1$).
We used this design to examine the evolution of floral turnover between the start date and alive recruits from a given census. The floral composition at start date comprised alive and dead survivors (AS and DS). Thus to study floral turnover between start date and alive recruits from a given census, we need to compare alive recruits (AR) to survivors (AS and DS).
Diversity indices
Diversity of communities can be partitioned into different level of diversity. We chose to use $\beta$-diversity, as extra information is obtained by knowing each community’s composition, compared to the mean community composition ($\alpha$-diversity). We more specifically used the Shannon diversity index. Shannon index, is an order 1 diversity index, giving a trade-off between the measurement of richness and evenness in floral composition. Shannon index of $\beta$-diversity was calculated with Tsallis generalized entropy formula [@marcon14b], and its exponential which gives the equivalent number of communities, with the following equations :
$$^1 H_{\beta} = \sum_{i}{w_{i} \sum_{s}{{p_{si} ~ ln \frac{p_{si}}{p_{s}} } }} ~~ and ~~ ^1 D_{\beta} = e^{ ^1 H_{\beta}} $$
with $^1 H_{\beta}$ the $\beta$-entropy and $^1 D_{\beta}$ the $\beta$-diversity of order 1 (Shannon) for $i$ communities of weights $w_{i}$ containing $s$ species of probability $p_{si}$ in the community $i$ and probability $p_{s}$ in the metacommunity.
Bias was reduced using the ‘best’ transformation following recommendations from @marcon14a.
All analyses were conducted using the R 3.3.0 software environment [@R] with entropart package [@entropart].
Results
Number of stems in alive recruits (AR) community from logged plots went from $2073 \pm 307$ in 1988 to $4132 \pm 301$ in 2015, representing $132 \pm 9$ and $172 \pm 6$ vernacular names respectively. Whereas, number of stems in survivors (S) from logged plots was $1780 \pm 279$ representing $113 \pm 8$ vernacular names
$\beta$-diversity of Tsallis entropy of order 1 ($\beta {treatment}$) between survivors (S) and alive recruits (AR) communities for all treatments is represented by $figure~1$. This result show that $\beta {T0}$ (unlogged) is almost stable in time with a mean of $1.08 \pm 0.018$. Consequently, a slow turnover exists on unlogged plots. There is an effect of logging represented by $\beta {logged}$. $\beta$-diversity from logged plots ($\beta {logged}$) are 4.4 times higher than $\beta {T0}$ (if we consider $\beta = 0$ as a null turnover) for the whole period. Moreover, $\beta$-diversity increases with treatment intensity ($\beta {logged} < \beta {logged} < \beta {logged}$). As unlogged plots $\beta$-diversity, the $\beta$-diversities are stable in time for the two lowest intensity treatment plots (T1 and T2). For the most intense treatment (T3), $\beta$-diversity has a tendency to increase between 1990 and 2003.
$\beta$-diversity rates remaining stable in time, log-abundance distributions should be similar from one year to the other. Consequently, we chose to look 2015 log-abundance figures. The log-abundance figures of the two community (AR and S) are relatively similar between all logged plots in 2015 (see $figure~2$). Survivors (S) distribution is always lower than the alive recruits' one. Log-abundance of logged plots are inverted between survivors (S) and alive recruits (AR) compare to unlogged plots. However, first-rank species abundance of survivors (S) in logged plots are bigger than the one of alive recruits (AR).
We can observe a pike in 1988 for unlogged plots. It is due to a smaller headcount of tree and species in the alive recruits (AR) community compare to survivors (S). It is shown by the low ratio (0.1) between the number of species of AR and S in 1988 (see $supplementary~materials~1$). Contrary to unlogged plots, the selective logging liberated space allowing many tree species to be recruited. So, alive recruits (AR) headcount is sufficient for the logged plots and they did not present this problem. The ratio for unlogged plots $supplementary~materials~2$).
We did not find clear differences in $\beta$-diversity between the four treatment for order 0 entropy (see $supplementary~materials~2$). Consequently, previous results concern mainly evenness (not present at entropy of order 0) and not richness. Still, $\beta$-diversity of order 0 is upper than 1 but mingled in time.
Additionally, if we consider more evenness by using $\beta$-diversity of order 2, results remains the same than order 1 concerning the comparison between logged and unlogged plots. But the differences between the three intensities of logging diminish (see $supplementary~materials~3$).
If we consider more evenness by using $\beta$-diversity at order 2, the results is still the same than order 1 for the comparison between logged and T0, but the differences between the three levels of log treatments disappear (see $supplementary~materials~3$).
rm(list = ls()) ; invisible(gc())
library(knitr)
library(parallel)
cores <- detectCores() - 1
library(Paracou)
opts_chunk$set( echo = F, message = F, warning = F, include = F, fig.height = 8, fig.width = 8, cache = T, cache.lazy = F)
tree <- data.table::fread('~/Documents/BIOGET/Projet/Paracou/inst/extdata/paracou_p1_15.csv', data.table = F)
# Extracting interest data
tree <- tree[c('n_parcelle', 'i_arbre', 'nomPilote', 'campagne', 'Genre', 'espece', 'code_vivant')]
names(tree) <- c('plot', 'id', 'vern', 'census', 'genus', 'species', 'alive')
# Treatment
tree$treatment <- NA
tree$treatment[which(tree$plot %in% c(1, 6, 11))] <- 0 # Unlogged
tree$treatment[which(tree$plot %in% c(2, 7, 9))] <- 1 # T1
tree$treatment[which(tree$plot %in% c(3, 5, 10))] <- 2 # T2
tree$treatment[which(tree$plot %in% c(4, 8, 12))] <- 3 # T3
bd <- tree[which(!duplicated(tree$id)), c('id', 'treatment', 'plot', 'vern', 'genus', 'species')]
bd <- merge(bd, tree[which(tree$alive == 'FAUX'), c('id', 'census')], all.x = T)
names(bd)[7] <- 'death'
# Missing years
bad_years <- c(1996, 1998, 2000, 2002, 2004, 2006)
age <- table(tree[-which(tree$alive == 'FAUX'),'id'], tree[-which(tree$alive == 'FAUX'),'census'])
age <- as.data.frame.matrix(age)
age[as.character(bad_years)] <- age[as.character(bad_years - 1)]
age <- rowSums(age) - 1
age <- data.frame(id = names(age), age = age)
bd <- merge(bd, age)
sel <- which(is.na(bd$death))
bd$birth <- NA
bd$birth[-sel] <- bd$death[-sel] - bd$age[-sel]
bd$birth[sel] <- 2015 - bd$age[sel]
bd <- bd[which(!duplicated(bd$id)),]
row.names(bd) <- bd$id
library(entropart)
years <- c(1988:2015)
q = 1 # Tsallis entropy order
diversity <- sapply(years, function(y){
# Attributing communities
bd_y <- bd
bd_y$com <- com(tree, bd_y, y, 1987)
if(length(which(is.na(bd_y$com))) > 0){
bd_y <- bd_y[-which(is.na(bd_y$com)),]
}
# Building lists of treatment and repetition
Coms <- list(T0 = subset(bd_y, treatment == 0),
T1 = subset(bd_y, treatment == 1),
T2 = subset(bd_y, treatment == 2),
T3 = subset(bd_y, treatment == 3))
Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1),
P6 = subset(Coms$T0, plot == 6),
P11 = subset(Coms$T0, plot == 11))
Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2),
P7 = subset(Coms$T1, plot == 7),
P9 = subset(Coms$T1, plot == 9))
Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3),
P5 = subset(Coms$T2, plot == 5),
P10 = subset(Coms$T2, plot == 10))
Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4),
P8 = subset(Coms$T3, plot == 8),
P12 = subset(Coms$T3, plot == 12))
# Extracting searched beta-diversities
## T0.DS.AR
T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)})
T0.mc <- lapply(T0.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T0.S.AR <- lapply(T0.mc, function(x){BetaDiversity(x, q = q)$Total})
## T1.DS.AR
T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T1.df <- lapply(T1.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T1.mc <- lapply(T1.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T1.S.AR <- lapply(T1.mc, function(x){BetaDiversity(x, q = q)$Total})
## T2.DS.AR
T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T2.df <- lapply(T2.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T2.mc <- lapply(T2.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T2.S.AR <- lapply(T2.mc, function(x){BetaDiversity(x, q = q)$Total})
## T3.DS.AR
T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T3.df <- lapply(T3.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T3.mc <- lapply(T3.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T3.S.AR <- lapply(T3.mc, function(x){BetaDiversity(x, q = q)$Total})
return(list(
T0.S.AR = T0.S.AR,
T1.S.AR = T1.S.AR,
T2.S.AR = T2.S.AR,
T3.S.AR = T3.S.AR
))
}, simplify = F)
names(diversity) <- years
diversity <- lapply(diversity, `[`, names(diversity[[1]]))
diversity <- apply(do.call(rbind, diversity), 2, as.list)
diversity <- lapply(diversity, function(x){do.call('rbind', lapply(x, rbind))})
diversity <- lapply(diversity, function(x){row.names(x) <- years ; return(x)})
diversity <- lapply(diversity, data.frame)
library(ggplot2)
# data
diversity <- lapply(diversity, function(x){
data.frame(apply(x, 2, as.numeric), row.names = years)
})
diversity.stats <- lapply(diversity, function(x){
data.frame(year = years,
mean = apply(x, 1, mean, na.rm = T),
se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))}))
})
diversity.stats <- do.call('rbind', diversity.stats)
diversity.stats$coms <- as.factor(rep(names(diversity), each = length(years)))
# plot
ggplot(diversity.stats, aes(x = year, y = mean, colour = coms)) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se),
colour = "black", width=.1) +
geom_line() +
geom_point(size=2, shape=21, fill="white") +
xlab('year') +
ylab('Beta-diversity of order 1 (Shannon index)') +
scale_colour_hue(name="Beta-diversity index",
breaks=c("T0.S.AR",
"T1.S.AR",
"T2.S.AR",
"T3.S.AR"),
labels=c("AR and S in T0",
"AR and S in T1",
"AR and S in T2",
"AR and S in T3"),
l=40) +
expand_limits(y=1) +
theme_bw() +
theme(legend.justification = c(1,0),
legend.position = c(1,0))
Figure 2: Beta-diversity of order 1 (Shannon index) evolution after logging. Beta-diversity was measured with Shannon beta-diversity from Tsallis entropy at order 2 between alive recruits (AR) and survivors (S) in the four treatments: T0 the control, T1, T2, and T3 represented by red, green, light blue and dark blue lines respectively. Vertical bars represents the standard errors.
# Attributing communities
bd_y <- bd
bd_y$com <- com(tree, bd_y, 2015, 1987)
if(length(which(is.na(bd_y$com))) > 0){
bd_y <- bd_y[-which(is.na(bd_y$com)),]
}
# Building lists of treatment and repetition
Coms <- list(T0 = subset(bd_y, treatment == 0),
T1 = subset(bd_y, treatment == 1),
T2 = subset(bd_y, treatment == 2),
T3 = subset(bd_y, treatment == 3))
Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1),
P6 = subset(Coms$T0, plot == 6),
P11 = subset(Coms$T0, plot == 11))
Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2),
P7 = subset(Coms$T1, plot == 7),
P9 = subset(Coms$T1, plot == 9))
Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3),
P5 = subset(Coms$T2, plot == 5),
P10 = subset(Coms$T2, plot == 10))
Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4),
P8 = subset(Coms$T3, plot == 8),
P12 = subset(Coms$T3, plot == 12))
# Data
Coms <- lapply(Coms, function(Ti){
lapply(Ti, function(x) {
with(x, as.data.frame.matrix(table(vern, com)))
})
})
Coms <- lapply(Coms, function(Ti){
lapply(Ti, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x[c('AR', 'S')])
})
})
Coms <- lapply(Coms, function(Ti){
lapply(Ti, function(x) {
data.frame(AR = log(sort(x$AR, decreasing = T)),
S = log(sort(x$S, decreasing = T)),
rank = 1:length(x$AR))
})
})
Coms.df <- lapply(Coms, function(x){do.call('rbind',x)})
reps <- lapply(Coms, function(x){
unlist(lapply(x, function(y){
length(y[,1])
}))
})
Coms.df <- mapply(function(x, y){
x$plot <- as.factor(rep(1:3, y))
return(x)
}, x = Coms.df, y = reps, SIMPLIFY = F)
Coms.df <- lapply(Coms.df, function(x){
AR <- cbind(x[-2], 'AR')
names(AR) <- c('density', 'rank', 'plot', 'Com')
S <- cbind(x[-1], 'S')
names(S) <- names(AR)
return(rbind(AR, S))
})
Coms.df <- lapply(Coms.df, function(x){
x[which(x$density > 0),]
})
# plot
p1 <- ggplot(Coms.df$T0, aes(x = rank, y = density, colour=Com)) +
geom_point(size=2, fill = "white", aes(shape=plot)) +
xlab('Rank') +
ylab('Log of abundance') +
scale_colour_hue(name="Community",
breaks=c("AR",
"S"),
labels=c("Alive recruits",
"Survivors"),
l=40) +
scale_shape(solid=F,
name="Plot",
labels=c("Plot 1",
"Plot 6",
"Plot 11")
) +
theme_bw() +
theme(legend.justification = c(1,1),
legend.position = c(1,1))
p2 <- ggplot(Coms.df$T1, aes(x = rank, y = density, colour=Com)) +
geom_point(size=2, fill = "white", aes(shape=plot)) +
xlab('Rank') +
ylab('Log of abundance') +
scale_colour_hue(name="Community",
breaks=c("AR",
"S"),
labels=c("Alive recruits",
"Survivors"),
l=40) +
scale_shape(solid=F,
name="Plot",
labels=c("Plot 2",
"Plot 7",
"Plot 9")
) +
theme_bw() +
theme(legend.justification = c(1,1),
legend.position = c(1,1))
p3 <- ggplot(Coms.df$T2, aes(x = rank, y = density, colour=Com)) +
geom_point(size=2, fill = "white", aes(shape=plot)) +
xlab('Rank') +
ylab('Log of abundance') +
scale_colour_hue(name="Community",
breaks=c("AR",
"S"),
labels=c("Alive recruits",
"Survivors"),
l=40) +
scale_shape(solid=F,
name="Plot",
labels=c("Plot 3",
"Plot 5",
"Plot 10")
) +
theme_bw() +
theme(legend.justification = c(1,1),
legend.position = c(1,1))
p4 <- ggplot(Coms.df$T3, aes(x = rank, y = density, colour=Com)) +
geom_point(size=2, fill = "white", aes(shape=plot)) +
xlab('Rank') +
ylab('Log of abundance') +
scale_colour_hue(name="Community",
breaks=c("AR",
"S"),
labels=c("Alive recruits",
"Survivors"),
l=40) +
scale_shape(solid=F,
name="Plot",
labels=c("Plot 4",
"Plot 8",
"Plot 12")
) +
theme_bw() +
theme(legend.justification = c(1,1),
legend.position = c(1,1))
cowplot::plot_grid(p1, p2, p3, p4, labels=c("T0", "T1", "T2", "T3"), ncol = 2, nrow = 2)
Figure 3: Rank abundance diagrams of alive recruits and survivors in 2015.
Discussion
Tropical forests are primary ecosystems in term of biodiversity and carbon storage [@lewis04]. These two aspects are key ecosystem functions, which support ecosystem services, such as timber production and climate regulation. Selective logging is increasing in tropical forests, and its long term impact is still relatively unknown. We thus need to study the impact of selective logging on tropical forest dynamic. Particularly, we need to assess resilience of tropical forest ecosystems to selective logging in order to preserve tropical forest ecosystem function and services. Paracou study design offer a perfect playground to study selective logging impact on forest dynamic. We concluded from Paracou data that 29 years after selective logging the forest didn't recover and does not seems to go toward the original equilibrium.
We found that logged forests had higher turnover rates that unlogged control, even 29 years after the log date. Moreover, turnover rates from logged plots remained stable in time with no trend of decrease. Our results are similar to findings of @osazuwa15b, at the exception that in our case turnover rates were not variable in time. It could be explain by the difference of environmental context, and we can also notice that our study contain much more time sampling.
Additionally, we found that the impact of selective logging on $\beta$-diversity mainly affect species distribution evenness instead of richness, similarly to @plumptre96. This result could explain why sometimes impact of selective logging on tropical forest biodiversity is minimized when only species richness is taken into account. Finally, our results also suggested a small trend of species turnover increase in high intensity logging forest, but this observation have to be investigated more deeply.
One limit of our study is the use of vernacular names to evaluate tree species diversity. Vernacular names induce a small discrepancy increasing evenness (Mirabel, personal communication). But the bias may affect similarly both recruits and survivors cohorts.
Majority of previous studies have investigated selective logging on tropical forests focusing a single time point, at short term, and only with stems or species richness. Few studies investigated temporal dynamics of tree species distribution. First our study show that, species distribution evenness is needed to assess selective logging impact on tropical forest. And our results suggest that models assuming mid-term recovery to selective logging are inadequate for forest management and conservation policy. Furthermore, our results question the tropical forest resilience to selective logging, implying biodiversity loss and carbon storage dysfunction. As suggested by @osazuwa15b, we need to rethink the view of benign effect of selective logging.
Acknowledgements
We would like to thanks Ariane Mirabelle, Camille Piponiot, and Bruno Hérault for their guiding in the study. The study won't have be possible without Pascal Petronelli and all people who have contributed to Paracou data collect during all the past years. Finally we would like to thank Timothy Chubb for English useful comments.
Supplementary materials
years <- c(1988:2015)
ratio <- sapply(years, function(y){
# Attributing communities
bd_y <- bd
bd_y$com <- com(tree, bd_y, y, 1987)
if(length(which(is.na(bd_y$com))) > 0){
bd_y <- bd_y[-which(is.na(bd_y$com)),]
}
# Building lists of treatment and repetition
Coms <- list(T0 = subset(bd_y, treatment == 0),
T1 = subset(bd_y, treatment == 1),
T2 = subset(bd_y, treatment == 2),
T3 = subset(bd_y, treatment == 3))
Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1),
P6 = subset(Coms$T0, plot == 6),
P11 = subset(Coms$T0, plot == 11))
Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2),
P7 = subset(Coms$T1, plot == 7),
P9 = subset(Coms$T1, plot == 9))
Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3),
P5 = subset(Coms$T2, plot == 5),
P10 = subset(Coms$T2, plot == 10))
Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4),
P8 = subset(Coms$T3, plot == 8),
P12 = subset(Coms$T3, plot == 12))
# Extracting searched beta-diversities
## T0.DS.AR
T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)})
T0.ratio <- lapply(T0.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))})
## T1.DS.AR
T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T1.df <- lapply(T1.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T1.ratio <- lapply(T1.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))})
## T2.DS.AR
T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T2.df <- lapply(T2.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T2.ratio <- lapply(T2.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))})
## T3.DS.AR
T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T3.df <- lapply(T3.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T3.ratio <- lapply(T3.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))})
return(list(
T0.ratio = T0.ratio,
T1.ratio = T1.ratio,
T2.ratio = T2.ratio,
T3.ratio = T3.ratio
))
}, simplify = F)
names(ratio) <- years
ratio <- lapply(ratio, `[`, names(ratio[[1]]))
ratio <- apply(do.call(rbind, ratio), 2, as.list)
ratio <- lapply(ratio, function(x){do.call('rbind', lapply(x, rbind))})
ratio <- lapply(ratio, function(x){row.names(x) <- years ; return(x)})
ratio <- lapply(ratio, data.frame)
library(ggplot2)
# data
ratio <- lapply(ratio, function(x){
data.frame(apply(x, 2, as.numeric), row.names = years)
})
ratio.stats <- lapply(ratio, function(x){
data.frame(year = years,
mean = apply(x, 1, mean, na.rm = T),
se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))}))
})
ratio.stats <- do.call('rbind', ratio.stats)
ratio.stats$coms <- as.factor(rep(names(ratio), each = length(years)))
# plot
ggplot(ratio.stats, aes(x = year, y = mean, colour = coms)) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se),
colour = "black", width=.1) +
geom_line() +
geom_point(size=2, shape=21, fill="white") +
xlab('year') +
ylab('Ratio') +
scale_colour_hue(name="Ratio",
breaks=c("T0.ratio",
"T1.ratio",
"T2.ratio",
"T3.ratio"),
labels=c("AR and S ratio in T0",
"AR and S ratio in T1",
"AR and S ratio in T2",
"AR and S ratio in T3"),
l=40) +
expand_limits(y=1) +
theme_bw() +
theme(legend.justification = c(1,0),
legend.position = c(1,0))
Supplementary materials 1: Ratio of number of species between alive recruits and survivors in time.
library(entropart)
years <- c(1988:2015)
q = 0 # Tsallis entropy order
diversity <- sapply(years, function(y){
# Attributing communities
bd_y <- bd
bd_y$com <- com(tree, bd_y, y, 1987)
if(length(which(is.na(bd_y$com))) > 0){
bd_y <- bd_y[-which(is.na(bd_y$com)),]
}
# Building lists of treatment and repetition
Coms <- list(T0 = subset(bd_y, treatment == 0),
T1 = subset(bd_y, treatment == 1),
T2 = subset(bd_y, treatment == 2),
T3 = subset(bd_y, treatment == 3))
Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1),
P6 = subset(Coms$T0, plot == 6),
P11 = subset(Coms$T0, plot == 11))
Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2),
P7 = subset(Coms$T1, plot == 7),
P9 = subset(Coms$T1, plot == 9))
Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3),
P5 = subset(Coms$T2, plot == 5),
P10 = subset(Coms$T2, plot == 10))
Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4),
P8 = subset(Coms$T3, plot == 8),
P12 = subset(Coms$T3, plot == 12))
# Extracting searched beta-diversities
## T0.DS.AR
T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)})
T0.mc <- lapply(T0.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T0.S.AR <- lapply(T0.mc, function(x){BetaDiversity(x, q = q)$Total})
## T1.DS.AR
T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T1.df <- lapply(T1.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T1.mc <- lapply(T1.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T1.S.AR <- lapply(T1.mc, function(x){BetaDiversity(x, q = q)$Total})
## T2.DS.AR
T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T2.df <- lapply(T2.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T2.mc <- lapply(T2.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T2.S.AR <- lapply(T2.mc, function(x){BetaDiversity(x, q = q)$Total})
## T3.DS.AR
T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T3.df <- lapply(T3.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T3.mc <- lapply(T3.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T3.S.AR <- lapply(T3.mc, function(x){BetaDiversity(x, q = q)$Total})
return(list(
T0.S.AR = T0.S.AR,
T1.S.AR = T1.S.AR,
T2.S.AR = T2.S.AR,
T3.S.AR = T3.S.AR
))
}, simplify = F)
names(diversity) <- years
diversity <- lapply(diversity, `[`, names(diversity[[1]]))
diversity <- apply(do.call(rbind, diversity), 2, as.list)
diversity <- lapply(diversity, function(x){do.call('rbind', lapply(x, rbind))})
diversity <- lapply(diversity, function(x){row.names(x) <- years ; return(x)})
diversity <- lapply(diversity, data.frame)
library(ggplot2)
# data
diversity <- lapply(diversity, function(x){
data.frame(apply(x, 2, as.numeric), row.names = years)
})
diversity.stats <- lapply(diversity, function(x){
data.frame(year = years,
mean = apply(x, 1, mean, na.rm = T),
se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))}))
})
diversity.stats <- do.call('rbind', diversity.stats)
diversity.stats$coms <- as.factor(rep(names(diversity), each = length(years)))
# plot
ggplot(diversity.stats, aes(x = year, y = mean, colour = coms)) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se),
colour = "black", width=.1) +
geom_line() +
geom_point(size=2, shape=21, fill="white") +
xlab('year') +
ylab('Beta-diversity at order 0 (Richness index)') +
scale_colour_hue(name="Beta-diversity index",
breaks=c("T0.S.AR",
"T1.S.AR",
"T2.S.AR",
"T3.S.AR"),
labels=c("AR and S in T0",
"AR and S in T1",
"AR and S in T2",
"AR and S in T3"),
l=40) +
expand_limits(y=1) +
theme_bw() +
theme(legend.justification = c(1,0),
legend.position = c(1,0))
Supplementary materials 2: Beta-diversity of order 0 (Richness index) evolution after logging. Beta-diversity was measured with Richness beta-diversity from Tsallis entropy at order 0 between alive recruits (AR) and survivors (S) in the four treatments: T0 the control, T1, T2, and T3 represented by red, green, light blue and dark blue lines respectively. Vertical bars represents the standard errors.
library(entropart)
years <- c(1988:2015)
q = 2 # Tsallis entropy order
diversity <- sapply(years, function(y){
# Attributing communities
bd_y <- bd
bd_y$com <- com(tree, bd_y, y, 1987)
if(length(which(is.na(bd_y$com))) > 0){
bd_y <- bd_y[-which(is.na(bd_y$com)),]
}
# Building lists of treatment and repetition
Coms <- list(T0 = subset(bd_y, treatment == 0),
T1 = subset(bd_y, treatment == 1),
T2 = subset(bd_y, treatment == 2),
T3 = subset(bd_y, treatment == 3))
Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1),
P6 = subset(Coms$T0, plot == 6),
P11 = subset(Coms$T0, plot == 11))
Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2),
P7 = subset(Coms$T1, plot == 7),
P9 = subset(Coms$T1, plot == 9))
Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3),
P5 = subset(Coms$T2, plot == 5),
P10 = subset(Coms$T2, plot == 10))
Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4),
P8 = subset(Coms$T3, plot == 8),
P12 = subset(Coms$T3, plot == 12))
# Extracting searched beta-diversities
## T0.DS.AR
T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)})
T0.mc <- lapply(T0.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T0.S.AR <- lapply(T0.mc, function(x){BetaDiversity(x, q = q)$Total})
## T1.DS.AR
T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T1.df <- lapply(T1.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T1.mc <- lapply(T1.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T1.S.AR <- lapply(T1.mc, function(x){BetaDiversity(x, q = q)$Total})
## T2.DS.AR
T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T2.df <- lapply(T2.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T2.mc <- lapply(T2.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T2.S.AR <- lapply(T2.mc, function(x){BetaDiversity(x, q = q)$Total})
## T3.DS.AR
T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))})
T3.df <- lapply(T3.df, function(x) {
if('AS' %in% names(x)){
x$S <- x$AS + x$DS
} else {
x$S <- x$DS
}
return(x)
})
T3.mc <- lapply(T3.df, function(x) {MetaCommunity(x[c('S', 'AR')])})
T3.S.AR <- lapply(T3.mc, function(x){BetaDiversity(x, q = q)$Total})
return(list(
T0.S.AR = T0.S.AR,
T1.S.AR = T1.S.AR,
T2.S.AR = T2.S.AR,
T3.S.AR = T3.S.AR
))
}, simplify = F)
names(diversity) <- years
diversity <- lapply(diversity, `[`, names(diversity[[1]]))
diversity <- apply(do.call(rbind, diversity), 2, as.list)
diversity <- lapply(diversity, function(x){do.call('rbind', lapply(x, rbind))})
diversity <- lapply(diversity, function(x){row.names(x) <- years ; return(x)})
diversity <- lapply(diversity, data.frame)
library(ggplot2)
# data
diversity <- lapply(diversity, function(x){
data.frame(apply(x, 2, as.numeric), row.names = years)
})
diversity.stats <- lapply(diversity, function(x){
data.frame(year = years,
mean = apply(x, 1, mean, na.rm = T),
se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))}))
})
diversity.stats <- do.call('rbind', diversity.stats)
diversity.stats$coms <- as.factor(rep(names(diversity), each = length(years)))
# plot
ggplot(diversity.stats, aes(x = year, y = mean, colour = coms)) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se),
colour = "black", width=.1) +
geom_line() +
geom_point(size=2, shape=21, fill="white") +
xlab('year') +
ylab('Beta-diversity at order 2 (Simpson index)') +
scale_colour_hue(name="Beta-diversity index",
breaks=c("T0.S.AR",
"T1.S.AR",
"T2.S.AR",
"T3.S.AR"),
labels=c("AR and S in T0",
"AR and S in T1",
"AR and S in T2",
"AR and S in T3"),
l=40) +
expand_limits(y=1) +
theme_bw() +
theme(legend.justification = c(1,0),
legend.position = c(1,0))
Supplementary materials 3: Beta-diversity of order 2 (Simpson index) evolution after logging. Beta-diversity was measured with Simpson beta-diversity from Tsallis entropy at order 1 between alive recruits (AR) and survivors (S) in the four treatments: T0 the control, T1, T2, and T3 represented by red, green, light blue and dark blue lines respectively. Vertical bars represents the standard errors.
References
sylvainschmitt/Paracou documentation built on May 30, 2019, 10:44 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Abstract
Keywords
Beta-diversity, Selective logging, Tropical evergreen forest
Introduction
Tropical logging accounts for one eighth of global timber production [@blaser11]. It represents a major threat to biodiversity [@gibson11]. Most tropical timber production originates from selective logging, the targeted harvesting of timber from commercial species in a single cutting cycle [@martin15]. Others sylvicultural practices, such as tree removal by poisoning, are sometimes also involved in selective logging.
The short-term negative impact of selective logging on biodiversity is well known [@deAvila15; @carreno12; @gibson11; @martin15]. @deAvila15 found that additional sylvicultural practices and logging damages had stronger negative effect on diversity than logging itself. @carreno12 found no direct effect of medium intensity logging on species diversity but a community drift towards 'fast' [@reich14] acquisitive species. However, even the short term impact on biodiversity is unclear, as sometimes short term increases in diversity can occur at low sylvicultural intensities [@martin15]. Much less is known about the long term impact [@osazuwa15a; @osazuwa15b].
There are a number of ways diversity can be measured. The first ways is to count the number of species, known as species richness. This measurement highlight the richness but not the evenness of the abundance distribution of species within communities. Second, the Shannon and Simpson indices takes more evenness into account, particularly the latter [@shannon49]. Finally, @tsallis88 combined these measurements with the help of entropy, i.e. the amount of uncertainty calculated from the frequency distribution of a community [@marcon14b]. @jost07 helped to clarify the concept of partitioning diversity $\gamma$, the total diversity of the metacommunity into $\beta$ and $\alpha$ diversities. But the debate is ongoing, in spite of an attempt to resolve it by @chao12.
$\alpha$-diversity, the local diversity, has been highly documented in tropical forests [@valencia94]. Whereas $\beta$-diversity, how species composition changes with distances and among communities, has been less investigated [@condit02]. Still, $\beta$-diversity is a central concept about diversity organisation inside tropical forests. $\beta$-diversity can result from species adaptations to climates or substrates [@paoli06], from limited dispersal [@qian09], or from other historical effects such as logging. Consequently $\beta$-diversity seems at least as important as $\alpha$-diversity for conservation and forest management.
In this paper we designed communities from individual cohorts with similar history. And we used $\beta$-diversity indices between those communities to study and compare turnover of logged and unlogged forests along time.
Material and methods
Study site
The study was carried out in Paracou research station ($5^\circ 18'N$ and $52^\circ 53'W$) in French Guiana. The site has a tropical wet climate, with an average annual rainfall of 2920mm and temperature of $26^\circ C$ marked with a 4 month dry season and a 8 month rainy season interrupted by a short drier period. The forest is a lowland tropical rain forest with a dominance of Fabaceae, Chrysobalanaceae, Lecythidaceae and Sapotaceae.
In 1984, twelves plots of 6.25ha were established to study the effects of different intensities of selective logging. From 1986 to 1988, four logging treatments (T0, T1, T2 and T3) were applied to the plots according to randomized block design with three replicates [@baraloto12]:
- T0 plots are the control plots were no selective logging or other sylvicultural practices were applied.
- T1 plots underwent selective felling for timber of 50-60cm diameter at breast height (dbh) trees from 58 commercial species at an intensity of about 10 trees/ha. Selective logging in T1 plots resulted in 12-33% above ground biomass (AGB) loss and $41.7 \pm 4.0 ~ ( mean \pm standard~deviation)$ species logged.
- T2 plots underwent selective felling for timber of about 30 trees/ha from same commercial species. Additionally, non-commercial species (NCS) over 50cm dbh and commercial species with major defects were eliminated by poison-girdling. Selective logging in T2 plots resulted in 33-56% AGB loss and $41.7 \pm 6.7$ species logged.
- T3 plots underwent selective felling for timber from same commercial species plus logging of NCS of 40-50cm dbh for fuelwood. Additionally, NCS over 50cm dbh and commercial species with major defects were eliminated by poison-girdling. Selective logging in T3 plots resulted in 35-56% AGB loss and $53.0 \pm 3.0$ species logged.
Floral composition
Before logging, each tree was identified by its vernacular name. Vernacular names can encompass different species and sometimes different genera. Since 1990, survivors and recruited trees have begun to be identified by species name. Consequently, trees removed in 1987 or trees that died before they could be categorized were not identified to the species level. But use of vernacular names has been shown to only induce a small discrepancy (Mirabel, personal communications). Our analysis will then use vernacular names. Censuses included 263 vernacular names.
Analysis
Time sampling
Inventories of all plots were conducted annually from 1984 to 1994, every two years from 1995 to 2012, and annually from 2012 on trees above 10 cm gbh. For our analysis we generate recruit figures for missing years by replicating that of the previous year (n-1) in order to simulate turnover measurements.
Community design
To study floral turnovers we divided the plots into communities with similar histories. We used the logging year as the start date for both logged and unlogged plots. Trees which attain 10cm dbh after this date were defined as recruits (R) whereas trees which had 10cm dbh before this date were defined as survivors (S), both in logged and unlogged plots. In each census, recruits and survivors were divided in two communities: alive trees (A) and dead trees (D). We thus obtain 4 communities: Alive Recruits (AR), Alive Survivors (AS), Dead Recruits (DR), and Dead Survivors (DS) for each census in both logged and unlogged plots (see $figure~1$).
We used this design to examine the evolution of floral turnover between the start date and alive recruits from a given census. The floral composition at start date comprised alive and dead survivors (AS and DS). Thus to study floral turnover between start date and alive recruits from a given census, we need to compare alive recruits (AR) to survivors (AS and DS).
Diversity indices
Diversity of communities can be partitioned into different level of diversity. We chose to use $\beta$-diversity, as extra information is obtained by knowing each community’s composition, compared to the mean community composition ($\alpha$-diversity). We more specifically used the Shannon diversity index. Shannon index, is an order 1 diversity index, giving a trade-off between the measurement of richness and evenness in floral composition. Shannon index of $\beta$-diversity was calculated with Tsallis generalized entropy formula [@marcon14b], and its exponential which gives the equivalent number of communities, with the following equations :
$$^1 H_{\beta} = \sum_{i}{w_{i} \sum_{s}{{p_{si} ~ ln \frac{p_{si}}{p_{s}} } }} ~~ and ~~ ^1 D_{\beta} = e^{ ^1 H_{\beta}} $$ with $^1 H_{\beta}$ the $\beta$-entropy and $^1 D_{\beta}$ the $\beta$-diversity of order 1 (Shannon) for $i$ communities of weights $w_{i}$ containing $s$ species of probability $p_{si}$ in the community $i$ and probability $p_{s}$ in the metacommunity.
Bias was reduced using the ‘best’ transformation following recommendations from @marcon14a.
All analyses were conducted using the R 3.3.0 software environment [@R] with entropart package [@entropart].
Results
Number of stems in alive recruits (AR) community from logged plots went from $2073 \pm 307$ in 1988 to $4132 \pm 301$ in 2015, representing $132 \pm 9$ and $172 \pm 6$ vernacular names respectively. Whereas, number of stems in survivors (S) from logged plots was $1780 \pm 279$ representing $113 \pm 8$ vernacular names
$\beta$-diversity of Tsallis entropy of order 1 ($\beta {treatment}$) between survivors (S) and alive recruits (AR) communities for all treatments is represented by $figure~1$. This result show that $\beta {T0}$ (unlogged) is almost stable in time with a mean of $1.08 \pm 0.018$. Consequently, a slow turnover exists on unlogged plots. There is an effect of logging represented by $\beta {logged}$. $\beta$-diversity from logged plots ($\beta {logged}$) are 4.4 times higher than $\beta {T0}$ (if we consider $\beta = 0$ as a null turnover) for the whole period. Moreover, $\beta$-diversity increases with treatment intensity ($\beta {logged} < \beta {logged} < \beta {logged}$). As unlogged plots $\beta$-diversity, the $\beta$-diversities are stable in time for the two lowest intensity treatment plots (T1 and T2). For the most intense treatment (T3), $\beta$-diversity has a tendency to increase between 1990 and 2003.
$\beta$-diversity rates remaining stable in time, log-abundance distributions should be similar from one year to the other. Consequently, we chose to look 2015 log-abundance figures. The log-abundance figures of the two community (AR and S) are relatively similar between all logged plots in 2015 (see $figure~2$). Survivors (S) distribution is always lower than the alive recruits' one. Log-abundance of logged plots are inverted between survivors (S) and alive recruits (AR) compare to unlogged plots. However, first-rank species abundance of survivors (S) in logged plots are bigger than the one of alive recruits (AR).
We can observe a pike in 1988 for unlogged plots. It is due to a smaller headcount of tree and species in the alive recruits (AR) community compare to survivors (S). It is shown by the low ratio (0.1) between the number of species of AR and S in 1988 (see $supplementary~materials~1$). Contrary to unlogged plots, the selective logging liberated space allowing many tree species to be recruited. So, alive recruits (AR) headcount is sufficient for the logged plots and they did not present this problem. The ratio for unlogged plots $supplementary~materials~2$).
We did not find clear differences in $\beta$-diversity between the four treatment for order 0 entropy (see $supplementary~materials~2$). Consequently, previous results concern mainly evenness (not present at entropy of order 0) and not richness. Still, $\beta$-diversity of order 0 is upper than 1 but mingled in time.
Additionally, if we consider more evenness by using $\beta$-diversity of order 2, results remains the same than order 1 concerning the comparison between logged and unlogged plots. But the differences between the three intensities of logging diminish (see $supplementary~materials~3$).
If we consider more evenness by using $\beta$-diversity at order 2, the results is still the same than order 1 for the comparison between logged and T0, but the differences between the three levels of log treatments disappear (see $supplementary~materials~3$).
rm(list = ls()) ; invisible(gc()) library(knitr) library(parallel) cores <- detectCores() - 1 library(Paracou) opts_chunk$set( echo = F, message = F, warning = F, include = F, fig.height = 8, fig.width = 8, cache = T, cache.lazy = F)
tree <- data.table::fread('~/Documents/BIOGET/Projet/Paracou/inst/extdata/paracou_p1_15.csv', data.table = F) # Extracting interest data tree <- tree[c('n_parcelle', 'i_arbre', 'nomPilote', 'campagne', 'Genre', 'espece', 'code_vivant')] names(tree) <- c('plot', 'id', 'vern', 'census', 'genus', 'species', 'alive') # Treatment tree$treatment <- NA tree$treatment[which(tree$plot %in% c(1, 6, 11))] <- 0 # Unlogged tree$treatment[which(tree$plot %in% c(2, 7, 9))] <- 1 # T1 tree$treatment[which(tree$plot %in% c(3, 5, 10))] <- 2 # T2 tree$treatment[which(tree$plot %in% c(4, 8, 12))] <- 3 # T3
bd <- tree[which(!duplicated(tree$id)), c('id', 'treatment', 'plot', 'vern', 'genus', 'species')] bd <- merge(bd, tree[which(tree$alive == 'FAUX'), c('id', 'census')], all.x = T) names(bd)[7] <- 'death' # Missing years bad_years <- c(1996, 1998, 2000, 2002, 2004, 2006) age <- table(tree[-which(tree$alive == 'FAUX'),'id'], tree[-which(tree$alive == 'FAUX'),'census']) age <- as.data.frame.matrix(age) age[as.character(bad_years)] <- age[as.character(bad_years - 1)] age <- rowSums(age) - 1 age <- data.frame(id = names(age), age = age) bd <- merge(bd, age) sel <- which(is.na(bd$death)) bd$birth <- NA bd$birth[-sel] <- bd$death[-sel] - bd$age[-sel] bd$birth[sel] <- 2015 - bd$age[sel] bd <- bd[which(!duplicated(bd$id)),] row.names(bd) <- bd$id
library(entropart) years <- c(1988:2015) q = 1 # Tsallis entropy order diversity <- sapply(years, function(y){ # Attributing communities bd_y <- bd bd_y$com <- com(tree, bd_y, y, 1987) if(length(which(is.na(bd_y$com))) > 0){ bd_y <- bd_y[-which(is.na(bd_y$com)),] } # Building lists of treatment and repetition Coms <- list(T0 = subset(bd_y, treatment == 0), T1 = subset(bd_y, treatment == 1), T2 = subset(bd_y, treatment == 2), T3 = subset(bd_y, treatment == 3)) Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1), P6 = subset(Coms$T0, plot == 6), P11 = subset(Coms$T0, plot == 11)) Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2), P7 = subset(Coms$T1, plot == 7), P9 = subset(Coms$T1, plot == 9)) Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3), P5 = subset(Coms$T2, plot == 5), P10 = subset(Coms$T2, plot == 10)) Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4), P8 = subset(Coms$T3, plot == 8), P12 = subset(Coms$T3, plot == 12)) # Extracting searched beta-diversities ## T0.DS.AR T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)}) T0.mc <- lapply(T0.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T0.S.AR <- lapply(T0.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T1.DS.AR T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T1.df <- lapply(T1.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T1.mc <- lapply(T1.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T1.S.AR <- lapply(T1.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T2.DS.AR T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T2.df <- lapply(T2.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T2.mc <- lapply(T2.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T2.S.AR <- lapply(T2.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T3.DS.AR T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T3.df <- lapply(T3.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T3.mc <- lapply(T3.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T3.S.AR <- lapply(T3.mc, function(x){BetaDiversity(x, q = q)$Total}) return(list( T0.S.AR = T0.S.AR, T1.S.AR = T1.S.AR, T2.S.AR = T2.S.AR, T3.S.AR = T3.S.AR )) }, simplify = F) names(diversity) <- years diversity <- lapply(diversity, `[`, names(diversity[[1]])) diversity <- apply(do.call(rbind, diversity), 2, as.list) diversity <- lapply(diversity, function(x){do.call('rbind', lapply(x, rbind))}) diversity <- lapply(diversity, function(x){row.names(x) <- years ; return(x)}) diversity <- lapply(diversity, data.frame)
library(ggplot2) # data diversity <- lapply(diversity, function(x){ data.frame(apply(x, 2, as.numeric), row.names = years) }) diversity.stats <- lapply(diversity, function(x){ data.frame(year = years, mean = apply(x, 1, mean, na.rm = T), se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))})) }) diversity.stats <- do.call('rbind', diversity.stats) diversity.stats$coms <- as.factor(rep(names(diversity), each = length(years))) # plot ggplot(diversity.stats, aes(x = year, y = mean, colour = coms)) + geom_errorbar(aes(ymin = mean - se, ymax = mean + se), colour = "black", width=.1) + geom_line() + geom_point(size=2, shape=21, fill="white") + xlab('year') + ylab('Beta-diversity of order 1 (Shannon index)') + scale_colour_hue(name="Beta-diversity index", breaks=c("T0.S.AR", "T1.S.AR", "T2.S.AR", "T3.S.AR"), labels=c("AR and S in T0", "AR and S in T1", "AR and S in T2", "AR and S in T3"), l=40) + expand_limits(y=1) + theme_bw() + theme(legend.justification = c(1,0), legend.position = c(1,0))
Figure 2: Beta-diversity of order 1 (Shannon index) evolution after logging. Beta-diversity was measured with Shannon beta-diversity from Tsallis entropy at order 2 between alive recruits (AR) and survivors (S) in the four treatments: T0 the control, T1, T2, and T3 represented by red, green, light blue and dark blue lines respectively. Vertical bars represents the standard errors.
# Attributing communities bd_y <- bd bd_y$com <- com(tree, bd_y, 2015, 1987) if(length(which(is.na(bd_y$com))) > 0){ bd_y <- bd_y[-which(is.na(bd_y$com)),] } # Building lists of treatment and repetition Coms <- list(T0 = subset(bd_y, treatment == 0), T1 = subset(bd_y, treatment == 1), T2 = subset(bd_y, treatment == 2), T3 = subset(bd_y, treatment == 3)) Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1), P6 = subset(Coms$T0, plot == 6), P11 = subset(Coms$T0, plot == 11)) Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2), P7 = subset(Coms$T1, plot == 7), P9 = subset(Coms$T1, plot == 9)) Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3), P5 = subset(Coms$T2, plot == 5), P10 = subset(Coms$T2, plot == 10)) Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4), P8 = subset(Coms$T3, plot == 8), P12 = subset(Coms$T3, plot == 12)) # Data Coms <- lapply(Coms, function(Ti){ lapply(Ti, function(x) { with(x, as.data.frame.matrix(table(vern, com))) }) }) Coms <- lapply(Coms, function(Ti){ lapply(Ti, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x[c('AR', 'S')]) }) }) Coms <- lapply(Coms, function(Ti){ lapply(Ti, function(x) { data.frame(AR = log(sort(x$AR, decreasing = T)), S = log(sort(x$S, decreasing = T)), rank = 1:length(x$AR)) }) }) Coms.df <- lapply(Coms, function(x){do.call('rbind',x)}) reps <- lapply(Coms, function(x){ unlist(lapply(x, function(y){ length(y[,1]) })) }) Coms.df <- mapply(function(x, y){ x$plot <- as.factor(rep(1:3, y)) return(x) }, x = Coms.df, y = reps, SIMPLIFY = F) Coms.df <- lapply(Coms.df, function(x){ AR <- cbind(x[-2], 'AR') names(AR) <- c('density', 'rank', 'plot', 'Com') S <- cbind(x[-1], 'S') names(S) <- names(AR) return(rbind(AR, S)) }) Coms.df <- lapply(Coms.df, function(x){ x[which(x$density > 0),] }) # plot p1 <- ggplot(Coms.df$T0, aes(x = rank, y = density, colour=Com)) + geom_point(size=2, fill = "white", aes(shape=plot)) + xlab('Rank') + ylab('Log of abundance') + scale_colour_hue(name="Community", breaks=c("AR", "S"), labels=c("Alive recruits", "Survivors"), l=40) + scale_shape(solid=F, name="Plot", labels=c("Plot 1", "Plot 6", "Plot 11") ) + theme_bw() + theme(legend.justification = c(1,1), legend.position = c(1,1)) p2 <- ggplot(Coms.df$T1, aes(x = rank, y = density, colour=Com)) + geom_point(size=2, fill = "white", aes(shape=plot)) + xlab('Rank') + ylab('Log of abundance') + scale_colour_hue(name="Community", breaks=c("AR", "S"), labels=c("Alive recruits", "Survivors"), l=40) + scale_shape(solid=F, name="Plot", labels=c("Plot 2", "Plot 7", "Plot 9") ) + theme_bw() + theme(legend.justification = c(1,1), legend.position = c(1,1)) p3 <- ggplot(Coms.df$T2, aes(x = rank, y = density, colour=Com)) + geom_point(size=2, fill = "white", aes(shape=plot)) + xlab('Rank') + ylab('Log of abundance') + scale_colour_hue(name="Community", breaks=c("AR", "S"), labels=c("Alive recruits", "Survivors"), l=40) + scale_shape(solid=F, name="Plot", labels=c("Plot 3", "Plot 5", "Plot 10") ) + theme_bw() + theme(legend.justification = c(1,1), legend.position = c(1,1)) p4 <- ggplot(Coms.df$T3, aes(x = rank, y = density, colour=Com)) + geom_point(size=2, fill = "white", aes(shape=plot)) + xlab('Rank') + ylab('Log of abundance') + scale_colour_hue(name="Community", breaks=c("AR", "S"), labels=c("Alive recruits", "Survivors"), l=40) + scale_shape(solid=F, name="Plot", labels=c("Plot 4", "Plot 8", "Plot 12") ) + theme_bw() + theme(legend.justification = c(1,1), legend.position = c(1,1)) cowplot::plot_grid(p1, p2, p3, p4, labels=c("T0", "T1", "T2", "T3"), ncol = 2, nrow = 2)
Figure 3: Rank abundance diagrams of alive recruits and survivors in 2015.
Discussion
Tropical forests are primary ecosystems in term of biodiversity and carbon storage [@lewis04]. These two aspects are key ecosystem functions, which support ecosystem services, such as timber production and climate regulation. Selective logging is increasing in tropical forests, and its long term impact is still relatively unknown. We thus need to study the impact of selective logging on tropical forest dynamic. Particularly, we need to assess resilience of tropical forest ecosystems to selective logging in order to preserve tropical forest ecosystem function and services. Paracou study design offer a perfect playground to study selective logging impact on forest dynamic. We concluded from Paracou data that 29 years after selective logging the forest didn't recover and does not seems to go toward the original equilibrium.
We found that logged forests had higher turnover rates that unlogged control, even 29 years after the log date. Moreover, turnover rates from logged plots remained stable in time with no trend of decrease. Our results are similar to findings of @osazuwa15b, at the exception that in our case turnover rates were not variable in time. It could be explain by the difference of environmental context, and we can also notice that our study contain much more time sampling.
Additionally, we found that the impact of selective logging on $\beta$-diversity mainly affect species distribution evenness instead of richness, similarly to @plumptre96. This result could explain why sometimes impact of selective logging on tropical forest biodiversity is minimized when only species richness is taken into account. Finally, our results also suggested a small trend of species turnover increase in high intensity logging forest, but this observation have to be investigated more deeply.
One limit of our study is the use of vernacular names to evaluate tree species diversity. Vernacular names induce a small discrepancy increasing evenness (Mirabel, personal communication). But the bias may affect similarly both recruits and survivors cohorts.
Majority of previous studies have investigated selective logging on tropical forests focusing a single time point, at short term, and only with stems or species richness. Few studies investigated temporal dynamics of tree species distribution. First our study show that, species distribution evenness is needed to assess selective logging impact on tropical forest. And our results suggest that models assuming mid-term recovery to selective logging are inadequate for forest management and conservation policy. Furthermore, our results question the tropical forest resilience to selective logging, implying biodiversity loss and carbon storage dysfunction. As suggested by @osazuwa15b, we need to rethink the view of benign effect of selective logging.
Acknowledgements
We would like to thanks Ariane Mirabelle, Camille Piponiot, and Bruno Hérault for their guiding in the study. The study won't have be possible without Pascal Petronelli and all people who have contributed to Paracou data collect during all the past years. Finally we would like to thank Timothy Chubb for English useful comments.
Supplementary materials
years <- c(1988:2015) ratio <- sapply(years, function(y){ # Attributing communities bd_y <- bd bd_y$com <- com(tree, bd_y, y, 1987) if(length(which(is.na(bd_y$com))) > 0){ bd_y <- bd_y[-which(is.na(bd_y$com)),] } # Building lists of treatment and repetition Coms <- list(T0 = subset(bd_y, treatment == 0), T1 = subset(bd_y, treatment == 1), T2 = subset(bd_y, treatment == 2), T3 = subset(bd_y, treatment == 3)) Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1), P6 = subset(Coms$T0, plot == 6), P11 = subset(Coms$T0, plot == 11)) Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2), P7 = subset(Coms$T1, plot == 7), P9 = subset(Coms$T1, plot == 9)) Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3), P5 = subset(Coms$T2, plot == 5), P10 = subset(Coms$T2, plot == 10)) Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4), P8 = subset(Coms$T3, plot == 8), P12 = subset(Coms$T3, plot == 12)) # Extracting searched beta-diversities ## T0.DS.AR T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)}) T0.ratio <- lapply(T0.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))}) ## T1.DS.AR T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T1.df <- lapply(T1.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T1.ratio <- lapply(T1.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))}) ## T2.DS.AR T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T2.df <- lapply(T2.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T2.ratio <- lapply(T2.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))}) ## T3.DS.AR T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T3.df <- lapply(T3.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T3.ratio <- lapply(T3.df, function(x) {length(which(x$AR > 0)) / length(which(x$S > 0))}) return(list( T0.ratio = T0.ratio, T1.ratio = T1.ratio, T2.ratio = T2.ratio, T3.ratio = T3.ratio )) }, simplify = F) names(ratio) <- years ratio <- lapply(ratio, `[`, names(ratio[[1]])) ratio <- apply(do.call(rbind, ratio), 2, as.list) ratio <- lapply(ratio, function(x){do.call('rbind', lapply(x, rbind))}) ratio <- lapply(ratio, function(x){row.names(x) <- years ; return(x)}) ratio <- lapply(ratio, data.frame)
library(ggplot2) # data ratio <- lapply(ratio, function(x){ data.frame(apply(x, 2, as.numeric), row.names = years) }) ratio.stats <- lapply(ratio, function(x){ data.frame(year = years, mean = apply(x, 1, mean, na.rm = T), se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))})) }) ratio.stats <- do.call('rbind', ratio.stats) ratio.stats$coms <- as.factor(rep(names(ratio), each = length(years))) # plot ggplot(ratio.stats, aes(x = year, y = mean, colour = coms)) + geom_errorbar(aes(ymin = mean - se, ymax = mean + se), colour = "black", width=.1) + geom_line() + geom_point(size=2, shape=21, fill="white") + xlab('year') + ylab('Ratio') + scale_colour_hue(name="Ratio", breaks=c("T0.ratio", "T1.ratio", "T2.ratio", "T3.ratio"), labels=c("AR and S ratio in T0", "AR and S ratio in T1", "AR and S ratio in T2", "AR and S ratio in T3"), l=40) + expand_limits(y=1) + theme_bw() + theme(legend.justification = c(1,0), legend.position = c(1,0))
Supplementary materials 1: Ratio of number of species between alive recruits and survivors in time.
library(entropart) years <- c(1988:2015) q = 0 # Tsallis entropy order diversity <- sapply(years, function(y){ # Attributing communities bd_y <- bd bd_y$com <- com(tree, bd_y, y, 1987) if(length(which(is.na(bd_y$com))) > 0){ bd_y <- bd_y[-which(is.na(bd_y$com)),] } # Building lists of treatment and repetition Coms <- list(T0 = subset(bd_y, treatment == 0), T1 = subset(bd_y, treatment == 1), T2 = subset(bd_y, treatment == 2), T3 = subset(bd_y, treatment == 3)) Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1), P6 = subset(Coms$T0, plot == 6), P11 = subset(Coms$T0, plot == 11)) Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2), P7 = subset(Coms$T1, plot == 7), P9 = subset(Coms$T1, plot == 9)) Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3), P5 = subset(Coms$T2, plot == 5), P10 = subset(Coms$T2, plot == 10)) Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4), P8 = subset(Coms$T3, plot == 8), P12 = subset(Coms$T3, plot == 12)) # Extracting searched beta-diversities ## T0.DS.AR T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)}) T0.mc <- lapply(T0.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T0.S.AR <- lapply(T0.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T1.DS.AR T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T1.df <- lapply(T1.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T1.mc <- lapply(T1.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T1.S.AR <- lapply(T1.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T2.DS.AR T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T2.df <- lapply(T2.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T2.mc <- lapply(T2.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T2.S.AR <- lapply(T2.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T3.DS.AR T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T3.df <- lapply(T3.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T3.mc <- lapply(T3.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T3.S.AR <- lapply(T3.mc, function(x){BetaDiversity(x, q = q)$Total}) return(list( T0.S.AR = T0.S.AR, T1.S.AR = T1.S.AR, T2.S.AR = T2.S.AR, T3.S.AR = T3.S.AR )) }, simplify = F) names(diversity) <- years diversity <- lapply(diversity, `[`, names(diversity[[1]])) diversity <- apply(do.call(rbind, diversity), 2, as.list) diversity <- lapply(diversity, function(x){do.call('rbind', lapply(x, rbind))}) diversity <- lapply(diversity, function(x){row.names(x) <- years ; return(x)}) diversity <- lapply(diversity, data.frame)
library(ggplot2) # data diversity <- lapply(diversity, function(x){ data.frame(apply(x, 2, as.numeric), row.names = years) }) diversity.stats <- lapply(diversity, function(x){ data.frame(year = years, mean = apply(x, 1, mean, na.rm = T), se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))})) }) diversity.stats <- do.call('rbind', diversity.stats) diversity.stats$coms <- as.factor(rep(names(diversity), each = length(years))) # plot ggplot(diversity.stats, aes(x = year, y = mean, colour = coms)) + geom_errorbar(aes(ymin = mean - se, ymax = mean + se), colour = "black", width=.1) + geom_line() + geom_point(size=2, shape=21, fill="white") + xlab('year') + ylab('Beta-diversity at order 0 (Richness index)') + scale_colour_hue(name="Beta-diversity index", breaks=c("T0.S.AR", "T1.S.AR", "T2.S.AR", "T3.S.AR"), labels=c("AR and S in T0", "AR and S in T1", "AR and S in T2", "AR and S in T3"), l=40) + expand_limits(y=1) + theme_bw() + theme(legend.justification = c(1,0), legend.position = c(1,0))
Supplementary materials 2: Beta-diversity of order 0 (Richness index) evolution after logging. Beta-diversity was measured with Richness beta-diversity from Tsallis entropy at order 0 between alive recruits (AR) and survivors (S) in the four treatments: T0 the control, T1, T2, and T3 represented by red, green, light blue and dark blue lines respectively. Vertical bars represents the standard errors.
library(entropart) years <- c(1988:2015) q = 2 # Tsallis entropy order diversity <- sapply(years, function(y){ # Attributing communities bd_y <- bd bd_y$com <- com(tree, bd_y, y, 1987) if(length(which(is.na(bd_y$com))) > 0){ bd_y <- bd_y[-which(is.na(bd_y$com)),] } # Building lists of treatment and repetition Coms <- list(T0 = subset(bd_y, treatment == 0), T1 = subset(bd_y, treatment == 1), T2 = subset(bd_y, treatment == 2), T3 = subset(bd_y, treatment == 3)) Coms$T0 <- list(P1 = subset(Coms$T0, plot == 1), P6 = subset(Coms$T0, plot == 6), P11 = subset(Coms$T0, plot == 11)) Coms$T1 <- list(P2 = subset(Coms$T1, plot == 2), P7 = subset(Coms$T1, plot == 7), P9 = subset(Coms$T1, plot == 9)) Coms$T2 <- list(P3 = subset(Coms$T2, plot == 3), P5 = subset(Coms$T2, plot == 5), P10 = subset(Coms$T2, plot == 10)) Coms$T3 <- list(P4 = subset(Coms$T3, plot == 4), P8 = subset(Coms$T3, plot == 8), P12 = subset(Coms$T3, plot == 12)) # Extracting searched beta-diversities ## T0.DS.AR T0.df <- lapply(Coms$T0, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T0.df <- lapply(T0.df, function(x) {x$S <- x$AS + x$DS ; return(x)}) T0.mc <- lapply(T0.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T0.S.AR <- lapply(T0.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T1.DS.AR T1.df <- lapply(Coms$T1, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T1.df <- lapply(T1.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T1.mc <- lapply(T1.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T1.S.AR <- lapply(T1.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T2.DS.AR T2.df <- lapply(Coms$T2, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T2.df <- lapply(T2.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T2.mc <- lapply(T2.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T2.S.AR <- lapply(T2.mc, function(x){BetaDiversity(x, q = q)$Total}) ## T3.DS.AR T3.df <- lapply(Coms$T3, function(x) {with(x, as.data.frame.matrix(table(vern, com)))}) T3.df <- lapply(T3.df, function(x) { if('AS' %in% names(x)){ x$S <- x$AS + x$DS } else { x$S <- x$DS } return(x) }) T3.mc <- lapply(T3.df, function(x) {MetaCommunity(x[c('S', 'AR')])}) T3.S.AR <- lapply(T3.mc, function(x){BetaDiversity(x, q = q)$Total}) return(list( T0.S.AR = T0.S.AR, T1.S.AR = T1.S.AR, T2.S.AR = T2.S.AR, T3.S.AR = T3.S.AR )) }, simplify = F) names(diversity) <- years diversity <- lapply(diversity, `[`, names(diversity[[1]])) diversity <- apply(do.call(rbind, diversity), 2, as.list) diversity <- lapply(diversity, function(x){do.call('rbind', lapply(x, rbind))}) diversity <- lapply(diversity, function(x){row.names(x) <- years ; return(x)}) diversity <- lapply(diversity, data.frame)
library(ggplot2) # data diversity <- lapply(diversity, function(x){ data.frame(apply(x, 2, as.numeric), row.names = years) }) diversity.stats <- lapply(diversity, function(x){ data.frame(year = years, mean = apply(x, 1, mean, na.rm = T), se = apply(x, 1, function(x){sd(x, na.rm = T) / sqrt(length(x))})) }) diversity.stats <- do.call('rbind', diversity.stats) diversity.stats$coms <- as.factor(rep(names(diversity), each = length(years))) # plot ggplot(diversity.stats, aes(x = year, y = mean, colour = coms)) + geom_errorbar(aes(ymin = mean - se, ymax = mean + se), colour = "black", width=.1) + geom_line() + geom_point(size=2, shape=21, fill="white") + xlab('year') + ylab('Beta-diversity at order 2 (Simpson index)') + scale_colour_hue(name="Beta-diversity index", breaks=c("T0.S.AR", "T1.S.AR", "T2.S.AR", "T3.S.AR"), labels=c("AR and S in T0", "AR and S in T1", "AR and S in T2", "AR and S in T3"), l=40) + expand_limits(y=1) + theme_bw() + theme(legend.justification = c(1,0), legend.position = c(1,0))
Supplementary materials 3: Beta-diversity of order 2 (Simpson index) evolution after logging. Beta-diversity was measured with Simpson beta-diversity from Tsallis entropy at order 1 between alive recruits (AR) and survivors (S) in the four treatments: T0 the control, T1, T2, and T3 represented by red, green, light blue and dark blue lines respectively. Vertical bars represents the standard errors.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.