-
Notifications
You must be signed in to change notification settings - Fork 0
/
Fig_5_accu_curve.R
165 lines (142 loc) · 5.08 KB
/
Fig_5_accu_curve.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
# Install packages----
# install.packages("remotes")
# remotes::install_github("KaiHsiangHu/iNEXT.3D")
# vignette: https://cran.r-project.org/web/packages/iNEXT/vignettes/Introduction.pdf
# library ----
library(tidyverse)
library(patchwork)
library(iNEXT.3D)
library(viridis)
library(gt)
# Data----
Site_dat <- read.csv(
"biomass_data.csv",
header = T,
fill = TRUE,
sep = ",",
na.strings = c("", " ", "NA", "NA ", "na", "NULL")
)
# Data wrangling----
Site_dat <- Site_dat %>%
mutate(
Treatment = as.factor(Treatment),
Life_form = as.factor(Life_form),
Functional_groups = as.factor(Functional_groups)
)
# Analysis----
Site_prep_iNext <- Site_dat %>%
mutate(species = Sci_name) %>%
arrange(Site, Treatment) %>%
select(-c(Palatability, Sci_name)) %>%
mutate(
Treatment = case_when(Treatment == "ab" ~ "Control", # Cymbopogon present fire present
Treatment == "bgpnf" ~ "CPFA", # Cymbopogon present fire absent
Treatment == "bgrnf" ~ "CAFA" # Cymbopogon absent fire absent
),
) %>%
mutate(Cymb_status = case_when(
!species == "Cymbopogon sp." ~ "Other Sp.",
TRUE ~ as.character(species)
),) %>%
mutate(pres = as.numeric(1)) %>%
mutate(Treatment = factor(Treatment)) %>% # to order treatments in the plot
mutate(Treatment = fct_relevel(Treatment, c("Control","CPFA","CAFA")))
Site.list <- Site_prep_iNext %>%
split(.$Treatment)
Site.matrix.list <- purrr::map(
Site.list,
~ .x %>%
select(species, Site, pres) %>%
distinct() %>%
spread(key = Site, value = pres) %>%
replace(is.na(.), 0) %>%
column_to_rownames(var = "species")
)
# Taxonomic diversity
TD_treat_out <-
iNEXT3D(
data = Site.matrix.list,
diversity = 'TD',
q = c(0, 1, 2),
datatype = 'incidence_raw',
size = c(1:60),
nboot = 0
)
TD_treat_out$DataInfo
# Assemblage = the treatment or groups
# T = Reference sample size
# U = Total number of incidents
# S.obs = Observed species richness
# SC = Sample coverage
TD_treat_out$AsyEst # to see the asymptotic diversity estimates
# Make df for ploting----
Site.TD.df <- TD_treat_out %>%
purrr::pluck("iNextEst", "size_based")
Site_info <- Site_prep_iNext %>%
distinct() %>% mutate(Assemblage = as.character(Treatment))
Site.hill.TD <- Site.TD.df %>% left_join(Site_info) %>%
mutate(Order.q = case_when(Order.q == "0" ~ "q = 0", # q=0 species richness
Order.q == "1" ~ "q = 1", # q=1 Shannon diversity
Order.q == "2" ~ "q = 2")) %>% # q=2 Simpson diversity or evenness
filter(!Order.q == "q = 1")
df.point <-
Site.hill.TD[which(Site.hill.TD$Method == "Observed"), ]
df.line <-
Site.hill.TD[which(Site.hill.TD$Method != "Observed"), ]
df.line$Method <- factor(df.line$Method,
c("Rarefaction", "Extrapolation"))
# Plot----
treatment_colors <- c("Control" = "#3b5d4d", "CPFA" = "#c5af99","CAFA" = "#ffd365")
fig_6 <- ggplot(Site.hill.TD ,
aes(x = nt, y = qD, color = Treatment)) +
facet_wrap( ~ Order.q) +
geom_point(aes(),
shape = 1,
size = 3,
data = df.point) +
geom_line(aes(linetype = Method), lwd = 0.75, data = df.line) +
labs(
x = "Number of sites",
y = "Taxonomic diversity",) +
scale_color_manual(values = treatment_colors) +
theme_bw(base_size = 12) +
theme(legend.text = element_text(size = 8)) +
guides(col = guide_legend(ncol = 15)) +
theme(plot.title = element_text(size = 14, hjust = 0.5)) +
theme(
panel.grid.major = element_line(colour = "gray86", size = 0.1),
panel.background = element_rect(fill = "white")
)
fig_6
(
fig_6 +
theme(plot.caption = element_text(
size = 8, face = "italic",
hjust = 0.0
)) +theme(
legend.position = c(1.3, .95),
legend.justification = c('right', 'top')
) +
theme(legend.background = element_rect(fill = NA))
)
ggsave('fig_6_acu_curve.jpg',
width = 10,
height = 6,
dpi = 300)
TD_treat_out$DataInfo # S.obs= Observed species richness will be the same, this is the q=0 in the plot
TD_treat_out$AsyEst
iNextEststimates <- as.data.frame(TD_treat_out$AsyEst) %>% rename(Treatment = Assemblage) %>% mutate(Treatment= as.factor(Treatment))%>% mutate(Treatment= fct_relevel(Treatment, c('Control', 'CPFA', 'CAFA'))) %>% arrange(Treatment)
iNextEststimates %>% mutate(Treatment= fct_relevel(Treatment, c('Control', 'CPFA', 'CAFA'))) %>% arrange(Treatment)
# iNext table
TableS7 <- iNextEststimates %>% select(Treatment, Diversity, Observed, Estimator) %>%
filter(Diversity!="Shannon diversity") %>%
mutate_if(is.numeric, round, 2) %>%
gt()%>%
tab_options(column_labels.font.size = 11,
table.font.size = 10,
column_labels.font.weight = "bold") %>%
tab_header(subtitle = '', 'Accumulation of species richness and evenness') %>%
opt_table_font(default_fonts()) %>% # Fonts: Roboto Mono,IBM Plex Mono, Red Hat Mono
opt_table_outline(style = "solid", width = px(2))
TableS7
TableS7 %>% gtsave('TableS7 (iNext_AsyEst).png', expand = 5) # expand to set white space