Renyi Diversity Profiles with vegan, BiodiversityR and ggplot2
Roeland Kindt
21/11/2020
1 Packages needed
library(BiodiversityR) # also loads vegan
library(ggplot2)
library(ggsci)
library(readxl)2 Introduction
The main objective of this document is to give some examples of how data on Renyi diversity profiles, obtained via vegan and BiodiversityR, can be plotted via ggplot2. I expect that researchers and students will modify the scripts given here to create their own graphs. This obviously requires some working knowledge of the ggplot2 package, for which I will not give details here on specific functions, arguments and choices that I selected (as an introduction to ggplot2, I would especially recommend a two-part online workshop given by Thomas Lin Pedersen).
I will also not dive deeply in the theory of diversity profiles. Most important details and key references are available from Chapter 5 - Analysis of diversity of the Tree Diversity Analysis book that I wrote with Ric Coe. This manual on community ecology and biodiversity analysis can be downloaded for free and is also the recommended citation for BiodiversityR.
Another method by which somebody may become familiar with plotting diversity profiles via ggplot2 is from the Graphical User Interface of BiodiversityR, available from function BiodiversityRGUI
In the two examples below, the same work flow is used of:
Step 1. Obtain the species accumulation result via function renyicomp
Step 2. Extract data for plotting via function renyicomp.long
Step 3. Use the extracted data to generate a plot with ggplot2
The theme for ggplot2 plotting will be modified as follows, resulting in a more empty plotting canvas than the default for ggplot2.
BioR.theme <- theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line("gray25"),
text = element_text(size = 12),
axis.text = element_text(size = 10, colour = "gray25"),
axis.title = element_text(size = 14, colour = "gray25"),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
legend.key = element_blank())3 A very brief explanation about diversity profiles
The reason that diversity profile methods are needed to rank communities by diversity can shown by the following example that I adapted from a keynote address that I gave many years ago at a biodiversity congress.
Community.A <- data.frame(t(c(10, 6, 4, 1)))
Community.B <- data.frame(t(c(17, rep(1, 7))))
Communities <- data.frame(Community=c(rep("Community A", 4), rep("Community B", 8)),
Abundance=c(10, 6, 4, 1, 17, rep(1, 7)))
Communities$Species <- factor(c(c(1:4), c(1:8)))Community A and B are significantly different in species richness (the number of species) and evenness:
ggplot(data=Communities, aes(x=Species, y=Abundance, colour=Species, label=Abundance)) +
geom_col(aes(fill=Community),
position = position_dodge2(width=1.0, preserve="single"), colour="black") +
geom_text(aes(y=Abundance+1),
position = position_dodge2(width=0.5, preserve="single"), colour="black", size=5) +
BioR.theme +
scale_fill_npg() +
facet_grid( ~ Community, scales = "free_x", space = "free_x")
Yet the Shannon diversity of these communities is very similar:
diversity(Community.A, index="shannon")## [1] 1.172066diversity(Community.B, index="shannon")## [1] 1.171194Whereas the richness of Community B is clearly higher (8 species > 4 species), the Simpson diversity index is higher for Community A.
diversity(Community.A, index="invsimpson")## [1] 2.882353diversity(Community.B, index="invsimpson")## [1] 1.945946Therefore, it is not possible to rank these two communities by diversity.
Community.AB <- data.frame(array(0, dim=c(2, 8)))
rownames(Community.AB) <- c("Community A", "Community B")
Community.AB[1, 1:4] <- Community.A
Community.AB[2, ] <- Community.B
renyi.example <- renyiresult(Community.AB, scales=c(0, 0.25, 0.5, 1, 2, 4, 8, Inf), method="s")
renyi.example## 0 0.25 0.5 1 2 4 8
## Community A 1.386294 1.319182 1.261737 1.172066 1.0586070 0.9411298 0.8454568
## Community B 2.079442 1.874799 1.639995 1.171194 0.6657482 0.4597594 0.3941034
## Inf
## Community A 0.7419373
## Community B 0.3448405The Renyi diversity profile (calculated above) allows to compare the diversity of communities. If the profile for community X is above the profile of Community Y over the entire range of their profiles, then community X is more diverse that community Y. If profiles cross over, as is the case with communities A and B, then it is not possible to order the communities by diversity.
plotgg0 <- ggplot(data=renyi.long(renyi.example), aes(x=Scales, y=Diversity, group=Grouping)) +
scale_x_discrete() +
scale_y_continuous(sec.axis = dup_axis(name=NULL)) +
geom_line(aes(colour=Grouping), size=2) +
geom_point(aes(colour=Grouping, shape=Grouping), alpha=0.8, size=5) +
BioR.theme +
scale_colour_npg() +
labs(x=expression(alpha), y="Diversity", colour="Community", shape="Community")
plotgg0
The Renyi diversity profile at scale 0 reflects species richness, at scale 1 the Shannon index, at scale 2 the Simpson index and at scale Inf the Berger-Parker index (the dominance of the most abundant species):
log(specnumber(Community.A)) - renyi.example[1, "0"]## [1] 0diversity(Community.A, index="shannon") - renyi.example[1, "1"]## [1] 0log(diversity(Community.A, index="invsimpson")) - renyi.example[1, "2"]## [1] 0log(sum(Community.A)/max(Community.A)) - renyi.example[1, "Inf"]## [1] 0A community with a perfectly even species distribution will have a horizontal diversity profile:
Community.C <- data.frame(t(rep(17, 8)))
renyiresult(Community.C)## 0 0.25 0.5 1 2 4 8 Inf
## all 2.079442 2.079442 2.079442 2.079442 2.079442 2.079442 2.079442 2.079442Two communities that have the same evenness …
Community.D <- data.frame(t(c(10, 10, 6, 6, 4, 4, 1, 1)))
Communities <- data.frame(Community=c(rep("Community A", 4), rep("Community D", 8)),
Abundance=c(10, 6, 4, 1, 10, 10, 6, 6, 4, 4, 1, 1))
Communities$Species <- factor(c(c(1:4), c(1:8)))
Community.AD <- data.frame(array(0, dim=c(2, 8)))
rownames(Community.AD) <- c("Community A", "Community D")
Community.AD[1, 1:4] <- Community.A
Community.AD[2, ] <- Community.D
ggplot(data=Communities, aes(x=Species, y=Abundance, colour=Species, label=Abundance)) +
geom_col(aes(fill=Community),
position = position_dodge2(width=1.0, preserve="single"), colour="black") +
geom_text(aes(y=Abundance+1),
position = position_dodge2(width=0.5, preserve="single"), colour="black", size=5) +
BioR.theme +
scale_fill_npg() +
facet_grid( ~ Community, scales = "free_x", space = "free_x")
… will have parallel diversity profiles.
renyi.example <- renyiresult(Community.AD, scales=c(0, 0.25, 0.5, 1, 2, 4, 8, Inf), method="s")
plotgg0 <- ggplot(data=renyi.long(renyi.example), aes(x=Scales, y=Diversity, group=Grouping)) +
scale_x_discrete() +
scale_y_continuous(sec.axis = dup_axis(name=NULL)) +
geom_line(aes(colour=Grouping), size=2) +
geom_point(aes(colour=Grouping, shape=Grouping), alpha=0.8, size=5) +
BioR.theme +
scale_colour_npg() +
labs(x=expression(alpha), y="Diversity", colour="Community", shape="Community")
plotgg0
renyi.example[2, ] - renyi.example[1, ]## 0 0.25 0.5 1 2 4
## Community D 0.6931472 0.6931472 0.6931472 0.6931472 0.6931472 0.6931472
## 8 Inf
## Community D 0.6931472 0.6931472log(2) # the difference reflects the difference in species richness only## [1] 0.69314724 Example 1: Dune meadow data
Typical for community ecology analysis in vegan and BiodiversityR, example data are included in two data sets, the community data set and the environmental data set. It is expected that the rows of these data sets correspond to the same sample sites.
The data in this example are data that were used as case study in Data Analysis in Community and Landscape Ecology and also in the Tree Diversity Analysis manual.
Data set dune is a community data set, where variables (columns) typically correspond to different species and data represents abundance of each species. Species names were abbreviated to eight characters, with for example Agrostol representing Agrostis stolonifera.
data(dune)
str(dune)## 'data.frame': 20 obs. of 30 variables:
## $ Achimill: num 1 3 0 0 2 2 2 0 0 4 ...
## $ Agrostol: num 0 0 4 8 0 0 0 4 3 0 ...
## $ Airaprae: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Alopgeni: num 0 2 7 2 0 0 0 5 3 0 ...
## $ Anthodor: num 0 0 0 0 4 3 2 0 0 4 ...
## $ Bellpere: num 0 3 2 2 2 0 0 0 0 2 ...
## $ Bromhord: num 0 4 0 3 2 0 2 0 0 4 ...
## $ Chenalbu: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Cirsarve: num 0 0 0 2 0 0 0 0 0 0 ...
## $ Comapalu: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Eleopalu: num 0 0 0 0 0 0 0 4 0 0 ...
## $ Elymrepe: num 4 4 4 4 4 0 0 0 6 0 ...
## $ Empenigr: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Hyporadi: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Juncarti: num 0 0 0 0 0 0 0 4 4 0 ...
## $ Juncbufo: num 0 0 0 0 0 0 2 0 4 0 ...
## $ Lolipere: num 7 5 6 5 2 6 6 4 2 6 ...
## $ Planlanc: num 0 0 0 0 5 5 5 0 0 3 ...
## $ Poaprat : num 4 4 5 4 2 3 4 4 4 4 ...
## $ Poatriv : num 2 7 6 5 6 4 5 4 5 4 ...
## $ Ranuflam: num 0 0 0 0 0 0 0 2 0 0 ...
## $ Rumeacet: num 0 0 0 0 5 6 3 0 2 0 ...
## $ Sagiproc: num 0 0 0 5 0 0 0 2 2 0 ...
## $ Salirepe: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Scorautu: num 0 5 2 2 3 3 3 3 2 3 ...
## $ Trifprat: num 0 0 0 0 2 5 2 0 0 0 ...
## $ Trifrepe: num 0 5 2 1 2 5 2 2 3 6 ...
## $ Vicilath: num 0 0 0 0 0 0 0 0 0 1 ...
## $ Bracruta: num 0 0 2 2 2 6 2 2 2 2 ...
## $ Callcusp: num 0 0 0 0 0 0 0 0 0 0 ...Data set dune.env is an environmental data set, where variables (columns) correspond to different descriptors (typically continuous and categorical variables) of the sample sites. One of the variables is Management, a categorical variable that describes different management categories, coded as BF (an abbreviation for biological farming), HF (hobby farming), NM (nature conservation management) and SF (standard farming).
data(dune.env)
summary(dune.env)## A1 Moisture Management Use Manure
## Min. : 2.800 1:7 BF:3 Hayfield:7 0:6
## 1st Qu.: 3.500 2:4 HF:5 Haypastu:8 1:3
## Median : 4.200 4:2 NM:6 Pasture :5 2:4
## Mean : 4.850 5:7 SF:6 3:4
## 3rd Qu.: 5.725 4:3
## Max. :11.5004.1 Step 1
First we obtain the BiodiversityR::renyicomp result. Similar to the calculation of species accumulation curves, the result reflects the diversity for the same number of sites for each level of a categorical variable. This is achieved by calculating a diversity accumulation surface for each of these levels first, then selecting the profiles for the number of sites = the number of sites of the factor levels with fewest sites (see details in the Tree Diversity Manual or from Kindt et al. 2006.
Renyi.1 <- renyicomp(dune, y=dune.env, factor='Management',
scales=c(0, 0.25, 0.5, 1, 2, 4, 8, Inf), permutations=1000, plotit=F)
Renyi.1## , , = mean
##
## scale
## Management 0 0.25 0.5 1 2 4 8
## BF 2.772589 2.716195 2.662749 2.566756 2.420204 2.252630 2.105191
## HF 2.947666 2.905714 2.866355 2.795894 2.685510 2.545415 2.402173
## NM 2.775906 2.731391 2.688070 2.606328 2.468467 2.295470 2.153456
## SF 2.784369 2.716481 2.652115 2.536473 2.360986 2.169637 2.022988
## scale
## Management Inf
## BF 1.871802
## HF 2.164092
## NM 1.950664
## SF 1.822491
##
## , , = stdev
##
## scale
## Management 0 0.25 0.5 1 2 4
## BF NA NA NA NA NA NA
## HF 0.08982381 0.08922659 0.08804295 0.08442440 0.07550271 0.06368178
## NM 0.14075376 0.13467084 0.12842097 0.11566838 0.09546145 0.09483038
## SF 0.08845501 0.08975961 0.09036065 0.08954927 0.08610534 0.09013808
## scale
## Management 8 Inf
## BF NA NA
## HF 0.06429136 0.07361085
## NM 0.11395439 0.11491759
## SF 0.09737947 0.10290850
##
## , , = min
##
## scale
## Management 0 0.25 0.5 1 2 4 8
## BF NA NA NA NA NA NA NA
## HF 2.708050 2.667574 2.631770 2.572529 2.489872 2.395696 2.307250
## NM 2.302585 2.286839 2.270992 2.239186 2.176726 2.068560 1.944138
## SF 2.639057 2.568654 2.504130 2.378479 2.189765 2.005640 1.851827
## scale
## Management Inf
## BF NA
## HF 2.055725
## NM 1.725510
## SF 1.677646
##
## , , = max
##
## scale
## Management 0 0.25 0.5 1 2 4 8
## BF NA NA NA NA NA NA NA
## HF 3.044522 3.002965 2.963682 2.892405 2.776243 2.632510 2.523664
## NM 2.944439 2.892382 2.839436 2.736108 2.598566 2.439474 2.355136
## SF 2.944439 2.881677 2.817994 2.692785 2.495236 2.337517 2.247298
## scale
## Management Inf
## BF NA
## HF 2.302585
## NM 2.151762
## SF 2.100061
##
## , , = Qnt 0.025
##
## scale
## Management 0 0.25 0.5 1 2 4 8
## BF NA NA NA NA NA NA NA
## HF 2.708050 2.667574 2.631770 2.572529 2.489872 2.395696 2.307250
## NM 2.302585 2.286839 2.270992 2.239186 2.176726 2.068560 1.944138
## SF 2.639057 2.568654 2.504130 2.378479 2.189765 2.005640 1.851827
## scale
## Management Inf
## BF NA
## HF 2.055725
## NM 1.725510
## SF 1.677646
##
## , , = Qnt 0.975
##
## scale
## Management 0 0.25 0.5 1 2 4 8
## BF NA NA NA NA NA NA NA
## HF 3.044522 3.002965 2.963682 2.892405 2.776243 2.632510 2.523664
## NM 2.944439 2.892382 2.839436 2.736108 2.598566 2.439474 2.355136
## SF 2.944439 2.881677 2.817994 2.692785 2.495236 2.337517 2.247298
## scale
## Management Inf
## BF NA
## HF 2.302585
## NM 2.151762
## SF 2.1000614.2 Step 2
To plot data via ggplot2, diversity profile data is rendered in the long format via function BiodiversityR::renyicomp.long. The label.freq option modifies the variable labelit that can be used to select how often a symbol will be plotted along the species accumulation curve (but here we keep the default).
renyi.long1 <- renyicomp.long(Renyi.1, label.freq=1)
head(renyi.long1)## Grouping Obs Scales Diversity Stdev Min Min.1 Min.2 LWR UPR labelit
## 1 BF 1 0 2.772589 NA NA NA NA NA NA TRUE
## 2 BF 2 0.25 2.716195 NA NA NA NA NA NA TRUE
## 3 BF 3 0.5 2.662749 NA NA NA NA NA NA TRUE
## 4 BF 4 1 2.566756 NA NA NA NA NA NA TRUE
## 5 BF 5 2 2.420204 NA NA NA NA NA NA TRUE
## 6 BF 6 4 2.252630 NA NA NA NA NA NA TRUE4.3 Step 3
Now we can generate the plot…
plotgg1 <- ggplot(data=renyi.long1, aes(x=Scales, y=Diversity, ymax=UPR, ymin=LWR)) +
scale_x_discrete() +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_line(data=renyi.long1,
aes(x=Obs, colour=Grouping),
size=2) +
geom_point(data=subset(renyi.long1, labelit==TRUE),
aes(colour=Grouping, shape=Grouping),
size=5) +
geom_ribbon(data=renyi.long1,
aes(x=Obs, colour=Grouping, fill=after_scale(alpha(colour, 0.1))),
show.legend=FALSE) +
BioR.theme +
scale_colour_npg() +
labs(x=expression(alpha), y = "Diversity", colour = "Management", shape = "Management")This is the resulting graph, showing a lot of overlap in the confidence intervals, but otherwise a clear order of average diversity.
plotgg1## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
To only show the average diversity profile, use the following script.
plotgg1 <- ggplot(data=renyi.long1, aes(x=Scales, y=Diversity, ymax=UPR, ymin=LWR)) +
scale_x_discrete() +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_line(data=renyi.long1,
aes(x=Obs, colour=Grouping),
size=2) +
geom_point(data=subset(renyi.long1, labelit==TRUE),
aes(colour=Grouping, shape=Grouping),
alpha=0.8, size=5) +
BioR.theme +
scale_colour_npg() +
labs(x=expression(alpha), y = "Diversity", colour = "Management", shape = "Management")
plotgg1
5 Example 2: Tree inventories from Panama
The second data sets document tree species composition and site characteristics for 1-ha plots sampled in Panama. The data set can be be downloaded from the supplementary data provided by Condit et al. 2002. Beta-diversity in tropical forest trees. Science 295: 666-669. (The same data set is used in a Rpub where I showed how to generate ordination diagrams with vegan, BiodiversityR and ggplot2 ).
As I had some problems with accessing the downloaded Excel file directly with readxl, I saved the file in the .xlsx format on a local directory in my computer.
Condit.down <- "C:\\BiodiversityR_test\\Condit.xlsx"The following script reads and formats the community data set.
Condit.1 <- as.data.frame(readxl::read_excel(Condit.down, sheet=1))
spec.columns <- as.character(Condit.1[, 1])
Condit.spec <- as.data.frame(t(Condit.1[, -1]))
names(Condit.spec) <- spec.columnsThe community data set has 100 sites (rows) and 802 species (columns):
dim(Condit.spec)## [1] 100 802Next the environmental data set is read:
Condit.2 <- as.data.frame(readxl::read_excel(Condit.down, sheet=2))## New names:
## * `` -> ...8site.rows <- as.character(Condit.2[, 1])
Condit.env <- Condit.2[, 2:7]
rownames(Condit.env) <- site.rows
Condit.env$Age.cat <- factor(Condit.env$age, levels=c("1", "2", "3"))In the script above, I created the new variable: Age.cat, a categorical variable corresponding to age category.
summary(Condit.env)## EW coord NS coord precip elev
## Min. :600714 Min. : 962862 Min. :1887 Min. : 10.0
## 1st Qu.:625929 1st Qu.:1011569 1st Qu.:2516 1st Qu.:120.0
## Median :626404 Median :1011819 Median :2530 Median :120.0
## Mean :631536 Mean :1013061 Mean :2581 Mean :150.4
## 3rd Qu.:637892 3rd Qu.:1012908 3rd Qu.:2530 3rd Qu.:120.0
## Max. :688165 Max. :1045987 Max. :4002 Max. :830.0
## NA's :2 NA's :2
## age geology Age.cat
## Min. :1.00 Length:100 1 :18
## 1st Qu.:2.00 Class :character 2 :15
## Median :3.00 Mode :character 3 :65
## Mean :2.48 NA's: 2
## 3rd Qu.:3.00
## Max. :3.00
## NA's :2Some further changes were made to the data sets, first to ensure that they had the same rownames (I checked first that the order of sites was the same).
check.datasets(Condit.spec, Condit.env)## Warning: rownames for community and environmental datasets are different# same rownames
rownames(Condit.spec) <- rownames(Condit.env)A next change was to remove cases with missing values.
Condit.env <- Condit.env[complete.cases(Condit.env), ]
Condit.spec <- same.sites(Condit.spec, Condit.env)
check.datasets(Condit.spec, Condit.env)## OKAnd a final change was to only use one (plot B27) of the 1-ha subplots of a 50-ha plot, in order to have a more even spatial distribution among the sites.
rows.include <- c(28, 51:98)
rownames(Condit.spec)[rows.include]## [1] "B27" "p1" "p2" "p3" "p4" "p5" "p6" "p7" "p8" "p9" "p10" "p11"
## [13] "p12" "p13" "p14" "p15" "p16" "p17" "p18" "p19" "p20" "p21" "p22" "p23"
## [25] "p24" "p25" "p26" "p27" "p28" "p29" "p30" "p31" "p32" "p33" "p34" "p35"
## [37] "p36" "p37" "p38" "p39" "C1" "C2" "C3" "C4" "S0" "S1" "S2" "S3"
## [49] "S4"Condit.spec2 <- Condit.spec[rows.include, ]
Condit.env2 <- Condit.env[rows.include, ]
check.datasets(Condit.spec2, Condit.env2)## OK5.1 Step 1
Now that we have a community and environmental data set in the proper format, we can proceed with step 1.
Renyi.2 <- renyicomp(Condit.spec2, y=Condit.env2, factor='Age.cat',
scales=c(0, 0.25, 0.5, 1, 2, 4, 8, Inf), permutations=1000, plotit=F)
Renyi.2## , , = mean
##
## scale
## Age.cat 0 0.25 0.5 1 2 4 8 Inf
## 1 5.931706 5.659103 5.383055 4.880992 4.227764 3.689357 3.340273 2.960063
## 2 5.680173 5.409075 5.151660 4.690496 3.934532 3.102367 2.675142 2.340815
## 3 6.243888 5.985644 5.708610 5.157543 4.392376 3.844386 3.537326 3.156063
##
## , , = stdev
##
## scale
## Age.cat 0 0.25 0.5 1 2 4
## 1 0.08435777 0.08481410 0.08326040 0.08151900 0.09481565 0.12182411
## 2 NA NA NA NA NA NA
## 3 0.05165209 0.06230859 0.07295866 0.08742688 0.08238485 0.06258282
## scale
## Age.cat 8 Inf
## 1 0.14835522 0.15634426
## 2 NA NA
## 3 0.05915276 0.06850305
##
## , , = min
##
## scale
## Age.cat 0 0.25 0.5 1 2 4 8 Inf
## 1 5.616771 5.339160 5.064835 4.576979 3.879414 3.287813 2.906058 2.544122
## 2 NA NA NA NA NA NA NA NA
## 3 6.035481 5.751046 5.453703 4.894488 4.189870 3.718710 3.391563 2.983457
##
## , , = max
##
## scale
## Age.cat 0 0.25 0.5 1 2 4 8 Inf
## 1 6.025866 5.760286 5.494676 5.023316 4.414129 3.956729 3.669467 3.334138
## 2 NA NA NA NA NA NA NA NA
## 3 6.326149 6.092742 5.854880 5.379010 4.661448 4.083298 3.795388 3.441388
##
## , , = Qnt 0.025
##
## scale
## Age.cat 0 0.25 0.5 1 2 4 8 Inf
## 1 5.676754 5.414149 5.161365 4.694434 4.031643 3.425239 3.042117 2.673368
## 2 NA NA NA NA NA NA NA NA
## 3 6.124683 5.842034 5.544218 4.969651 4.241360 3.743419 3.442612 3.036144
##
## , , = Qnt 0.975
##
## scale
## Age.cat 0 0.25 0.5 1 2 4 8 Inf
## 1 6.016157 5.748156 5.480288 5.002559 4.386654 3.893356 3.594226 3.251207
## 2 NA NA NA NA NA NA NA NA
## 3 6.317165 6.080198 5.825519 5.316432 4.563523 3.985795 3.659206 3.3160525.2 Step 2
Step 2 is the same as before:
renyi.long2 <- renyicomp.long(Renyi.2, label.freq=1)
head(renyi.long2)## Grouping Obs Scales Diversity Stdev Min Min.1 Min.2 LWR
## 1 1 1 0 5.931706 0.08435777 5.616771 5.616771 5.616771 5.676754
## 2 1 2 0.25 5.659103 0.08481410 5.339160 5.339160 5.339160 5.414149
## 3 1 3 0.5 5.383055 0.08326040 5.064835 5.064835 5.064835 5.161365
## 4 1 4 1 4.880992 0.08151900 4.576979 4.576979 4.576979 4.694434
## 5 1 5 2 4.227764 0.09481565 3.879414 3.879414 3.879414 4.031643
## 6 1 6 4 3.689357 0.12182411 3.287813 3.287813 3.287813 3.425239
## UPR labelit
## 1 6.016157 TRUE
## 2 5.748156 TRUE
## 3 5.480288 TRUE
## 4 5.002559 TRUE
## 5 4.386654 TRUE
## 6 3.893356 TRUE5.3 Step 3
And step 3 is also similar.
plotgg2 <- ggplot(data=renyi.long2, aes(x=Scales, y=Diversity, ymax=UPR, ymin=LWR)) +
scale_x_discrete() +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_line(data=renyi.long2,
aes(x=Obs, colour=Grouping),
size=2) +
geom_point(aes(colour=Grouping, shape=Grouping),
size=5) +
geom_ribbon(data=renyi.long2,
aes(x=Obs, colour=Grouping, fill=after_scale(alpha(colour, 0.1))),
show.legend=FALSE) +
BioR.theme +
scale_colour_npg() +
labs(x=expression(alpha), y = "Diversity", colour = "Age", shape = "Age")This is now the resulting plot.
plotgg2## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
6 Session Information
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252
## [2] LC_CTYPE=English_United Kingdom.1252
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.1252
##
## attached base packages:
## [1] tcltk stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] readxl_1.3.1 ggsci_2.9 ggplot2_3.3.2
## [4] BiodiversityR_2.12-2 rgl_0.100.54 vegan3d_1.1-2
## [7] vegan_2.5-6 lattice_0.20-41 permute_0.9-5
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_1.4-1 ellipsis_0.3.1
## [4] class_7.3-17 rio_0.5.16 htmlTable_2.0.0
## [7] base64enc_0.1-3 rstudioapi_0.11 farver_2.0.3
## [10] splines_4.0.2 knitr_1.28 effects_4.1-4
## [13] Formula_1.2-3 jsonlite_1.6.1 nloptr_1.2.2.1
## [16] cluster_2.1.0 Rcmdr_2.6-2 png_0.1-7
## [19] shiny_1.4.0.2 compiler_4.0.2 backports_1.1.7
## [22] Matrix_1.2-18 fastmap_1.0.1 survey_4.0
## [25] later_1.1.0.1 tcltk2_1.2-11 acepack_1.4.1
## [28] htmltools_0.5.0 tools_4.0.2 gtable_0.3.0
## [31] glue_1.4.1 dplyr_1.0.2 Rcpp_1.0.4.6
## [34] carData_3.0-4 cellranger_1.1.0 vctrs_0.3.4
## [37] nlme_3.1-148 crosstalk_1.1.0.1 xfun_0.15
## [40] stringr_1.4.0 openxlsx_4.1.5 lme4_1.1-23
## [43] mime_0.9 miniUI_0.1.1.1 lifecycle_0.2.0
## [46] RcmdrMisc_2.7-0 statmod_1.4.34 MASS_7.3-51.6
## [49] zoo_1.8-8 scales_1.1.1 hms_0.5.3
## [52] promises_1.1.1 parallel_4.0.2 sandwich_2.5-1
## [55] RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3
## [58] gridExtra_2.3 rpart_4.1-15 latticeExtra_0.6-29
## [61] stringi_1.4.6 nortest_1.0-4 e1071_1.7-3
## [64] checkmate_2.0.0 boot_1.3-25 zip_2.0.4
## [67] manipulateWidget_0.10.1 rlang_0.4.8 pkgconfig_2.0.3
## [70] evaluate_0.14 purrr_0.3.4 htmlwidgets_1.5.1
## [73] labeling_0.3 tidyselect_1.1.0 magrittr_1.5
## [76] R6_2.4.1 generics_0.1.0 Hmisc_4.4-0
## [79] DBI_1.1.0 pillar_1.4.4 haven_2.3.1
## [82] foreign_0.8-80 withr_2.2.0 mgcv_1.8-31
## [85] survival_3.1-12 scatterplot3d_0.3-41 abind_1.4-5
## [88] nnet_7.3-14 tibble_3.0.1 crayon_1.3.4
## [91] car_3.0-8 relimp_1.0-5 rmarkdown_2.3
## [94] jpeg_0.1-8.1 grid_4.0.2 data.table_1.12.8
## [97] forcats_0.5.0 digest_0.6.25 webshot_0.5.2
## [100] xtable_1.8-4 httpuv_1.5.4 munsell_0.5.0
## [103] mitools_2.4