library(tidyverse)
library(readxl)
library(janitor)5 Join Functions
In this chapter, we will learn how to merge information from another dataset (using left_join(), full_join()) and how to filter our data based on another dataset (semi_join(), anti_join()). We will also use our skills from the previous chapter to calculate various statistics about our data, such as community means and community-weighted means of selected traits or indicator values (using left_join(), group_by(), and summarise()).
In the picture below, you can see simple visualisation of two datasets joined by different approaches. The resulting dataset is depicted in grey.
Credits: https://r4ds.hadley.nz/joins.html
5.1 Matching selected information with left_join()
In most cases, the data you need for analyses is not organised in one file, but rather in multiple files. For example, we have measured characteristics, known as traits, for each species in the region, and we have species lists for all the sites. To calculate the mean values of selected traits for each site, we must first append the information to individual records. For this, we will use the function called left_join(), which belongs to the group of mutating joins available in the tidyverse package dplyr.
This function will automatically look for columns with exactly matching names, or we can specify the corresponding variables directly using the by argument. It will append all the columns available in the second “lifeforms” dataset, but only for rows where there is an exact match. Although we have information for other species, such as Abies alba or Carex sylvatica, these do not have matching rows in the first dataset named “data”, and the information will therefore not be used. This means that left_join() is only selecting rows that are present in the first dataset. If the information is missing in the matching file, for example, if there is no information for Trientalis europaea in the lifeforms dataset, the resulting file will contain NA.
Let’s train a bit. Start a new script for this chapter again.
We will import the data from the Forest understory again.
spe <- read_csv("data/forest_understory/Axmanova-Forest-spe.csv")Rows: 2277 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Species
dbl (3): RELEVE_NR, Layer, CoverPerc
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.We will also import data with some traits, such as plant height and Ellenberg-type indicator values, which indicate demands on environmental factors, including soil reaction, light, and nutrients. Lower values indicate lower demands, meaning that species can tolerate low (acidic) pH levels or low nutrient levels. High values suggest that the species requires or grows in habitats with high (basic) soil reaction and high nutrient availability.
traits <- read_csv("data/forest_understory/traits.csv")Rows: 316 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Taxon
dbl (6): Height, SeedMass, EIV_Light, EIV_Moisture, EIV_Reaction, EIV_Nutrients
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.If we try to merge these two files directly, it does not work. We will get a warning explaining that there are no common variables.
spe %>% left_join(traits)Error in `left_join()`:
! `by` must be supplied when `x` and `y` have no common variables.
ℹ Use `cross_join()` to perform a cross-join.So there are more options for what we can do. We check the data and find out that the corresponding/matching information is stored once as ‘Species’ (in the data) and once as ‘Taxon’ (in traits).
tibble(spe)# A tibble: 2,277 × 4
RELEVE_NR Species Layer CoverPerc
<dbl> <chr> <dbl> <dbl>
1 1 Quercus petraea agg. 1 63
2 1 Acer campestre 1 18
3 1 Crataegus laevigata 4 2
4 1 Cornus mas 4 8
5 1 Lonicera xylosteum 4 3
6 1 Galium sylvaticum 6 3
7 1 Carex digitata 6 3
8 1 Melica nutans 6 2
9 1 Polygonatum odoratum 6 2
10 1 Geum urbanum 6 2
# ℹ 2,267 more rowstibble(traits)# A tibble: 316 × 7
Taxon Height SeedMass EIV_Light EIV_Moisture EIV_Reaction EIV_Nutrients
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Quercus pe… 25 975. 6 5 5 4
2 Acer campe… 7.5 54.3 5 5 7 6
3 Crataegus … 5.5 276. 5 5 7 5
4 Cornus mas 4 488. 6 4 8 5
5 Lonicera x… 2 164. 4 5 7 6
6 Galium syl… 0.85 1.10 5 5 6 5
7 Carex digi… 0.2 1.44 4 5 6 6
8 Melica nut… 0.45 2.49 4 5 6 4
9 Polygonatu… 0.45 264. 6 3 7 3
10 Geum urban… 0.55 2.05 5 5 6 7
# ℹ 306 more rowsWe can decide to rename the ‘Species’ column to ‘Taxon’ and use this column for matching.
spe %>% rename(Taxon = Species) %>%
left_join(traits, by = "Taxon")# A tibble: 2,277 × 10
RELEVE_NR Taxon Layer CoverPerc Height SeedMass EIV_Light EIV_Moisture
<dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 Quercus pet… 1 63 25 975. 6 5
2 1 Acer campes… 1 18 7.5 54.3 5 5
3 1 Crataegus l… 4 2 5.5 276. 5 5
4 1 Cornus mas 4 8 4 488. 6 4
5 1 Lonicera xy… 4 3 2 164. 4 5
6 1 Galium sylv… 6 3 0.85 1.10 5 5
7 1 Carex digit… 6 3 0.2 1.44 4 5
8 1 Melica nuta… 6 2 0.45 2.49 4 5
9 1 Polygonatum… 6 2 0.45 264. 6 3
10 1 Geum urbanum 6 2 0.55 2.05 5 5
# ℹ 2,267 more rows
# ℹ 2 more variables: EIV_Reaction <dbl>, EIV_Nutrients <dbl>Or we can even skip the by argument, as there is only one common matching column.
spe %>% rename(Taxon = Species) %>%
left_join(traits)Joining with `by = join_by(Taxon)`# A tibble: 2,277 × 10
RELEVE_NR Taxon Layer CoverPerc Height SeedMass EIV_Light EIV_Moisture
<dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 Quercus pet… 1 63 25 975. 6 5
2 1 Acer campes… 1 18 7.5 54.3 5 5
3 1 Crataegus l… 4 2 5.5 276. 5 5
4 1 Cornus mas 4 8 4 488. 6 4
5 1 Lonicera xy… 4 3 2 164. 4 5
6 1 Galium sylv… 6 3 0.85 1.10 5 5
7 1 Carex digit… 6 3 0.2 1.44 4 5
8 1 Melica nuta… 6 2 0.45 2.49 4 5
9 1 Polygonatum… 6 2 0.45 264. 6 3
10 1 Geum urbanum 6 2 0.55 2.05 5 5
# ℹ 2,267 more rows
# ℹ 2 more variables: EIV_Reaction <dbl>, EIV_Nutrients <dbl>We can also skip the rename() step and directly say in the left_join() function what is matching with what. It will retain ‘Species’ as the outcome, because that is how we named it in the primary data file.
spe %>% left_join(traits, by = c("Species" = "Taxon"))# A tibble: 2,277 × 10
RELEVE_NR Species Layer CoverPerc Height SeedMass EIV_Light EIV_Moisture
<dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 Quercus pet… 1 63 25 975. 6 5
2 1 Acer campes… 1 18 7.5 54.3 5 5
3 1 Crataegus l… 4 2 5.5 276. 5 5
4 1 Cornus mas 4 8 4 488. 6 4
5 1 Lonicera xy… 4 3 2 164. 4 5
6 1 Galium sylv… 6 3 0.85 1.10 5 5
7 1 Carex digit… 6 3 0.2 1.44 4 5
8 1 Melica nuta… 6 2 0.45 2.49 4 5
9 1 Polygonatum… 6 2 0.45 264. 6 3
10 1 Geum urbanum 6 2 0.55 2.05 5 5
# ℹ 2,267 more rows
# ℹ 2 more variables: EIV_Reaction <dbl>, EIV_Nutrients <dbl>Note that left_join() is a great help. However, you should always check before, or know the data well. If you have multiple matching rows, you can create a lot of unwanted mess with left_join(). For example, if you do not use a unique ID, but a descriptor such as forest type. Please read the warning messages carefully. :-)
*Right_join works almost the same as the left_join. You just have to turn the logic. It keeps all the observations in the matching file.
5.2 full_join() vs. bind_rows()
Sometimes it is needed to match both datasets exactly as they are. First, we will look at the full_join. In contrast to left_join(), full_join() will also add additional rows that appear in the other dataset.
Credits: https://r4ds.hadley.nz/joins.html
We will upload the forest env data and make two subsets to train a bit.
env <- read_csv("data/forest_understory/Axmanova-Forest-env.csv")Rows: 65 Columns: 13
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): ForestTypeName
dbl (12): RELEVE_NR, ForestType, Biomass, pH_KCl, Canopy_E3, Radiation, Heat...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.In the first subset, we will use only the first two forest types, namely oak and oak-hornbeam forests and the following selection of variables.
forest12 <- env %>% filter(ForestType %in% c(1, 2)) %>%
select(PlotID = RELEVE_NR, ForestType, ForestTypeName, pH_KCl, Biomass)In the second subset, we will use one forest type shared with the previous subset and two others: oak-hornbeam, ravine, and alluvial forests, with a slightly different selection of variables.
forest234 <- env %>%
filter(ForestType %in% c(2, 3, 4))%>%
select(PlotID = RELEVE_NR, ForestType, ForestTypeName, pH_KCl, Canopy_E3)Now, left_join() the data and check the preview (use view()) to see more of the dataset. We have added the extra information available in the second dataset to our selection of plots and variables in the first dataset, but the number of rows remained the same as it was in the forest1 selection i.e. 44 rows of the oak and oak-hornbeam forest types.
forest12 %>%
left_join(forest234)Joining with `by = join_by(PlotID, ForestType, ForestTypeName, pH_KCl)`# A tibble: 44 × 6
PlotID ForestType ForestTypeName pH_KCl Biomass Canopy_E3
<dbl> <dbl> <chr> <dbl> <dbl> <dbl>
1 1 2 oak hornbeam forest 5.28 12.8 80
2 2 1 oak forest 3.24 9.9 NA
3 3 1 oak forest 4.01 15.2 NA
4 4 1 oak forest 3.76 16 NA
5 5 1 oak forest 3.50 20.7 NA
6 6 1 oak forest 3.8 46.4 NA
7 7 1 oak forest 3.48 49.2 NA
8 8 2 oak hornbeam forest 3.68 48.7 85
9 9 2 oak hornbeam forest 4.24 13.8 80
10 10 1 oak forest 4.00 79.1 NA
# ℹ 34 more rowsIn the next step, we will combine both subsets using a full_join(). What is the output? We will obtain a dataset containing 65 rows and all available information, encompassing all four forest types. However, since Forest type 1 is only present in the first subset, the corresponding rows will include only the information available there (forest type and pH). Type 2 is present in both datasets. Therefore, it will contain all information from both datasets (forest type, pH, biomass, and canopy). Types 3 and 4 will have information only from the second dataset (forest type, pH, canopy). Use view() to see the complete merged dataset.
forest12 %>%
full_join(forest234)Joining with `by = join_by(PlotID, ForestType, ForestTypeName, pH_KCl)`# A tibble: 65 × 6
PlotID ForestType ForestTypeName pH_KCl Biomass Canopy_E3
<dbl> <dbl> <chr> <dbl> <dbl> <dbl>
1 1 2 oak hornbeam forest 5.28 12.8 80
2 2 1 oak forest 3.24 9.9 NA
3 3 1 oak forest 4.01 15.2 NA
4 4 1 oak forest 3.76 16 NA
5 5 1 oak forest 3.50 20.7 NA
6 6 1 oak forest 3.8 46.4 NA
7 7 1 oak forest 3.48 49.2 NA
8 8 2 oak hornbeam forest 3.68 48.7 85
9 9 2 oak hornbeam forest 4.24 13.8 80
10 10 1 oak forest 4.00 79.1 NA
# ℹ 55 more rowsNote that the joining worked well, since we had the same names of variables in both datasets. For example, if the soil reaction were called soil_pH in one and pH_KCl in the other dataset, it will treat them as different variables, as you can see below:
forest12 %>%
rename(soil_pH = pH_KCl) %>%
full_join(forest234)Joining with `by = join_by(PlotID, ForestType, ForestTypeName)`# A tibble: 65 × 7
PlotID ForestType ForestTypeName soil_pH Biomass pH_KCl Canopy_E3
<dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
1 1 2 oak hornbeam forest 5.28 12.8 5.28 80
2 2 1 oak forest 3.24 9.9 NA NA
3 3 1 oak forest 4.01 15.2 NA NA
4 4 1 oak forest 3.76 16 NA NA
5 5 1 oak forest 3.50 20.7 NA NA
6 6 1 oak forest 3.8 46.4 NA NA
7 7 1 oak forest 3.48 49.2 NA NA
8 8 2 oak hornbeam forest 3.68 48.7 3.68 85
9 9 2 oak hornbeam forest 4.24 13.8 4.24 80
10 10 1 oak forest 4.00 79.1 NA NA
# ℹ 55 more rows*Tip: if you are merging data from different sources and you get a similar mess, check the function coalesce() https://dplyr.tidyverse.org/reference/coalesce.html.
forest12 %>%
rename(soil_pH = pH_KCl) %>%
full_join(forest234) %>%
mutate(pH=coalesce(soil_pH,pH_KCl)) %>%
select(PlotID, pH, soil_pH, pH_KCl) %>%
filter(PlotID <5 | PlotID>120)Joining with `by = join_by(PlotID, ForestType, ForestTypeName)`# A tibble: 16 × 4
PlotID pH soil_pH pH_KCl
<dbl> <dbl> <dbl> <dbl>
1 1 5.28 5.28 5.28
2 2 3.24 3.24 NA
3 3 4.01 4.01 NA
4 4 3.76 3.76 NA
5 121 5.76 5.76 5.76
6 124 4.89 4.89 4.89
7 126 4.88 4.88 4.88
8 122 6.68 NA 6.68
9 123 6.91 NA 6.91
10 125 4.99 NA 4.99
11 128 7.03 NA 7.03
12 129 5.7 NA 5.7
13 130 5.74 NA 5.74
14 131 5.22 NA 5.22
15 132 5.56 NA 5.56
16 127 5.56 NA 5.56Imagine a situation where data were sampled in different years/months or by different colleagues, or in different regions separately. If these subsets have exactly the same structure, we can use the bind_rows() function to merge them, one below the other.
forest1 <- env %>%
filter(ForestType == 1) %>%
select(PlotID = RELEVE_NR, ForestType, ForestTypeName, pH_KCl, Biomass)
forest2 <- env %>%
filter(ForestType == 2)%>%
select(PlotID = RELEVE_NR, ForestType, ForestTypeName, pH_KCl, Biomass)Bind and check the whole dataset with view().
forest1 %>%
bind_rows(forest2)# A tibble: 44 × 5
PlotID ForestType ForestTypeName pH_KCl Biomass
<dbl> <dbl> <chr> <dbl> <dbl>
1 2 1 oak forest 3.24 9.9
2 3 1 oak forest 4.01 15.2
3 4 1 oak forest 3.76 16
4 5 1 oak forest 3.50 20.7
5 6 1 oak forest 3.8 46.4
6 7 1 oak forest 3.48 49.2
7 10 1 oak forest 4.00 79.1
8 11 1 oak forest 3.83 7.4
9 28 1 oak forest 3.56 29.7
10 31 1 oak forest 3.20 14.9
# ℹ 34 more rowsNote that the rbind() function in base R works similarly, but less effectively and it is designed for data frames.
rm(forest1, forest2, forest12, forest234) #remove unwanted data from the environment5.3 *Inner_join
Inner_join compares the two files and keeps only rows that are common in both files. Therefore, it is easy to lose information, and we should be careful about this approach. In the example below, I have a dataset with forest types 1, 2 and 3 and information about biomass and another dataset with forest types 2, 3 and 4 and information about pH. In the merged file, inner_join keeps only the common rows with the same key= PlotID. Therefore, the result will include only vegetation types 2 and 3, but with all information available in both files (i.e., biomass and pH).
forest123 <- env %>%
filter(ForestType %in% c (1,2,3)) %>%
select(PlotID = RELEVE_NR, ForestType, ForestTypeName, Biomass)
forest234 <- env %>%
filter(ForestType %in% c (2,3,4))%>%
select(PlotID = RELEVE_NR, ForestType, ForestTypeName, pH_KCl)
forest<- forest123 %>% inner_join(forest234)Joining with `by = join_by(PlotID, ForestType, ForestTypeName)`tibble(forest)# A tibble: 38 × 5
PlotID ForestType ForestTypeName Biomass pH_KCl
<dbl> <dbl> <chr> <dbl> <dbl>
1 1 2 oak hornbeam forest 12.8 5.28
2 8 2 oak hornbeam forest 48.7 3.68
3 9 2 oak hornbeam forest 13.8 4.24
4 14 2 oak hornbeam forest 14.2 4.07
5 16 2 oak hornbeam forest 26.2 3.60
6 18 2 oak hornbeam forest 72.2 4.62
7 29 2 oak hornbeam forest 37.9 3.61
8 32 2 oak hornbeam forest 19.2 3.35
9 36 2 oak hornbeam forest 3.3 3.38
10 41 2 oak hornbeam forest 26.9 3.34
# ℹ 28 more rowsWith the lines below, I remove specified objects that I no longer need from the environment. rm is a shortcut for remove.
rm(forest, forest123, forest234)5.4 Filtering joins
Sometimes we need to filter the data according to another dataset. For example, we have reviewed the environmental file and determined that for a follow-up analysis, we only need the alluvial forest data. Therefore, we filtered rows in the env file.
forest4 <- env %>%
filter(ForestType == 4)%>%
select(PlotID = RELEVE_NR,ForestType, ForestTypeName, pH_KCl, Biomass) %>%
print()# A tibble: 11 × 5
PlotID ForestType ForestTypeName pH_KCl Biomass
<dbl> <dbl> <chr> <dbl> <dbl>
1 101 4 alluvial forest 6.67 91.1
2 103 4 alluvial forest 7.12 114.
3 104 4 alluvial forest 7.14 188.
4 110 4 alluvial forest 5.82 126.
5 111 4 alluvial forest 7.01 84.8
6 113 4 alluvial forest 6.61 74.5
7 125 4 alluvial forest 4.99 123.
8 128 4 alluvial forest 7.03 121.
9 130 4 alluvial forest 5.74 238
10 132 4 alluvial forest 5.56 287.
11 127 4 alluvial forest 5.56 176 However, we also need to filter the rows in the species data where there is no information about forest types. The only common identifier is the PlotID (i.e. sample/site ID).
spe %>% rename(PlotID = RELEVE_NR) %>%
semi_join(forest4)Joining with `by = join_by(PlotID)`# A tibble: 530 × 4
PlotID Species Layer CoverPerc
<dbl> <chr> <dbl> <dbl>
1 101 Fraxinus excelsior 1 38
2 103 Alnus glutinosa 1 63
3 104 Alnus glutinosa 1 63
4 110 Alnus glutinosa 1 38
5 111 Alnus glutinosa 1 38
6 113 Alnus glutinosa 1 63
7 125 Alnus glutinosa 1 38
8 128 Corylus avellana 1 38
9 130 Acer pseudoplatanus 1 38
10 132 Alnus glutinosa 1 38
# ℹ 520 more rowsWith the distinct() function, you can check if the IDs are the same.
spe %>% rename(PlotID = RELEVE_NR) %>%
semi_join(forest4) %>%
distinct(PlotID)Joining with `by = join_by(PlotID)`# A tibble: 11 × 1
PlotID
<dbl>
1 101
2 103
3 104
4 110
5 111
6 113
7 125
8 128
9 130
10 132
11 127Sometimes you may need the opposite - to keep only rows/IDs that are not in the other file. For this, you should use the anti_join() function. For example, you may want to remove a particular group of species from your list, such as annuals or non-vascular plants (mosses and lichens).
To show you how anti_join works, I will prepare an example of the list of herbs and filter the dataset of traits.
herbs <- spe %>%
filter(Layer==6) %>%
distinct(Species)First, I will filter the dataset of traits to keep only the data for the herbs with semi_join i.e. I will keep the same species.
traits %>% semi_join(herbs, by=c ("Taxon"= "Species"))# A tibble: 269 × 7
Taxon Height SeedMass EIV_Light EIV_Moisture EIV_Reaction EIV_Nutrients
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Galium syl… 0.85 1.10 5 5 6 5
2 Carex digi… 0.2 1.44 4 5 6 6
3 Melica nut… 0.45 2.49 4 5 6 4
4 Polygonatu… 0.45 264. 6 3 7 3
5 Geum urban… 0.55 2.05 5 5 6 7
6 Anemone sp… NA NA NA NA NA NA
7 Viola mira… 0.2 1.7 4 5 8 5
8 Hieracium … 0.4 0.4 5 5 5 4
9 Platanther… 0.35 NA 6 5 6 4
10 Convallari… 0.2 134. 4 5 5 4
# ℹ 259 more rowsNow I will filter the dataset of traits to keep everything except the herbs with anti_join. I will keep only the species not in the list of herbs.
traits %>% anti_join(herbs, by=c ("Taxon"= "Species"))# A tibble: 47 × 7
Taxon Height SeedMass EIV_Light EIV_Moisture EIV_Reaction EIV_Nutrients
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Quercus pe… 25 975. 6 5 5 4
2 Acer campe… 7.5 54.3 5 5 7 6
3 Crataegus … 5.5 276. 5 5 7 5
4 Cornus mas 4 488. 6 4 8 5
5 Lonicera x… 2 164. 4 5 7 6
6 Euonymus v… 2.5 NA 5 4 7 6
7 Fraxinus e… 27.5 58.4 4 6 7 7
8 Crataegus … NA NA NA NA NA NA
9 Ligustrum … 3.5 35.4 6 4 7 5
10 Rosa speci… NA NA NA NA NA NA
# ℹ 37 more rows5.5 Community means and community weighted means with summarise()
summarise() reduces the information to a single value per group, as defined by the group_by() function. We also used summarise() with across() to calculate the minimum, mean, median, maximum, or sum of values in the env data. A simple example is below.
env %>%
group_by(ForestTypeName) %>%
summarise(meanBiomass = mean(Biomass))# A tibble: 4 × 2
ForestTypeName meanBiomass
<chr> <dbl>
1 alluvial forest 148.
2 oak forest 26.7
3 oak hornbeam forest 44.4
4 ravine forest 79.1Here, we will use joins to add new possibilities for how to use summarise.
summarise() is often used for obtaining so-called community means or community-weighted means. For example, we have traits for individual species, and we want to calculate a mean for each site and compare it.
spe %>%
left_join(traits, by = c("Species" = "Taxon")) %>%
summarise(meanHeight = mean(Height), .by = RELEVE_NR)# A tibble: 65 × 2
RELEVE_NR meanHeight
<dbl> <dbl>
1 1 NA
2 2 6.00
3 3 3.50
4 4 8.14
5 5 3.82
6 6 NA
7 7 NA
8 8 3.52
9 9 8.25
10 10 NA
# ℹ 55 more rowsSo far, we have had data with no missing values, but you will most likely encounter situations where there are some NAs. For this, you will need to use the argument na.rm = TRUE, which indicates that you want to skip the NAs from counting. Otherwise, a single NA would mean that the result will also be NA. The same issue occurred here, so we need to address it.
spe %>%
left_join(traits, by = c("Species" = "Taxon")) %>%
summarise(meanHeight = mean(Height, na.rm = T), .by = RELEVE_NR)# A tibble: 65 × 2
RELEVE_NR meanHeight
<dbl> <dbl>
1 1 4.60
2 2 6.00
3 3 3.50
4 4 8.14
5 5 3.82
6 6 3.3
7 7 4.25
8 8 3.52
9 9 8.25
10 10 3.53
# ℹ 55 more rowsWe can also consider using abundance to weight the result. This is called the community-weighted mean, and it assigns higher importance to species that are more abundant and lower importance to those that are rare. Here, the abundance is approximated as a percentage cover in the plot.
spe %>%
left_join(traits, by = c("Species" = "Taxon")) %>%
summarise(meanHeight = mean(Height, na.rm = T),
meanHeight_weighted = weighted.mean(Height, CoverPerc, na.rm = TRUE),
.by = RELEVE_NR)# A tibble: 65 × 3
RELEVE_NR meanHeight meanHeight_weighted
<dbl> <dbl> <dbl>
1 1 4.60 12.9
2 2 6.00 16.3
3 3 3.50 16.7
4 4 8.14 17.2
5 5 3.82 15.3
6 6 3.3 9.36
7 7 4.25 12.9
8 8 3.52 11.2
9 9 8.25 14.6
10 10 3.53 14.9
# ℹ 55 more rowsYou can see that the results are pretty different! Note that we are in a forest, so we now compare the heights of the dominant trees and the smaller, less abundant herbs. Especially for the CWM, it plays a huge role. Therefore, we will now calculate the CWM only for the herb layer (6).
spe %>%
left_join(traits, by = c("Species" = "Taxon")) %>%
filter(Layer == 6) %>%
summarise(meanHeight = mean(Height, na.rm = T),
meanHeight_weighted = weighted.mean(Height, CoverPerc, na.rm = TRUE),
.by = RELEVE_NR)# A tibble: 65 × 3
RELEVE_NR meanHeight meanHeight_weighted
<dbl> <dbl> <dbl>
1 1 0.379 0.367
2 2 0.496 0.483
3 3 0.430 0.404
4 4 0.598 0.593
5 5 0.563 0.527
6 6 0.562 0.561
7 7 0.491 0.522
8 8 0.452 0.657
9 9 0.492 0.481
10 10 0.539 0.555
# ℹ 55 more rowsSimilarly, we can calculate community means across multiple variables simultaneously. Here, we will try the Ellenberg-type indicator values (abbreviated as EIV), which measure species demands for particular factors. The higher the value, the higher the demands or affinity to habitats with higher values of these environmental factors. We will apply them solely based on the presence of species; no weights will be used. Therefore, we do not need to care that much about the differences among layers. However, we will keep each species only once (a list of unique species for each plot), so we will remove information about layers first and use distinct(). An alternative is to add filter(Layer == 6) %>% if we want to focus on the herb layer only.
spe %>%
left_join(traits, by = c("Species" = "Taxon")) %>%
select(-c(Layer, CoverPerc)) %>%
distinct() %>%
group_by(PlotID=RELEVE_NR) %>%
summarise(across(starts_with("EIV"), ~mean(.x, na.rm = TRUE)))# A tibble: 65 × 5
PlotID EIV_Light EIV_Moisture EIV_Reaction EIV_Nutrients
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1 5.08 4.64 6.54 4.95
2 2 4.88 4.81 4.94 4.31
3 3 4.87 4.8 4.93 4.4
4 4 5.3 4.5 5.7 4.25
5 5 6.07 4.14 5.57 3.86
6 6 5.72 4.5 6.39 5.22
7 7 5.8 4.4 6.25 4.9
8 8 5.17 4.61 5.43 4.65
9 9 5.32 4.5 5.77 4.41
10 10 6.38 3.69 6.19 3.88
# ℹ 55 more rowsWhen you become confident with joins and summarise(), you can even do something like this: We are working with env data, but at the same time, we want to use the summarised information from species data. We can calculate it directly in the pipeline without stepping outside, as there is no need to calculate and save intermediate steps separately. We will train this in chapters 10 and 11.
env %>%
left_join(spe %>%
group_by(RELEVE_NR) %>%
count(name = "Richness")) %>%
select(PlotID = RELEVE_NR, ForestType, ForestTypeName, Richness)Joining with `by = join_by(RELEVE_NR)`# A tibble: 65 × 4
PlotID ForestType ForestTypeName Richness
<dbl> <dbl> <chr> <int>
1 1 2 oak hornbeam forest 44
2 2 1 oak forest 17
3 3 1 oak forest 16
4 4 1 oak forest 22
5 5 1 oak forest 15
6 6 1 oak forest 22
7 7 1 oak forest 23
8 8 2 oak hornbeam forest 24
9 9 2 oak hornbeam forest 28
10 10 1 oak forest 19
# ℹ 55 more rows5.6 Exercises
1. Prepare a script from this chapter with your notes.
2. Do a bit more training with Forest understory env data. Upload the env file called Axmanova-Forest-env.csv and create two subsets saved in your environment: forest1 with variables PlotID = RELEVE_NR, ForestType, ForestTypeName and forest2 with variables ForestType, pH_KCl, Biomass. Now use left_join() to match forest2 to forest1. What happened?
3. Get back to the Forest understory species file (Axmanova-Forest-spe.csv). In the data, you can see which species occur in the tree layer (=1) >> prepare a table where for each plotID you will have in another column a single name of the tree dominant, i.e. the tree with the highest cover (for simplicity, specify in the slice_max(n= 1, with_ties = F)>> What is the most common dominant? In how many plots?
4. Append this list of dominants to env data >> calculate counts of plots with each of the dominant for each forest type >> arrange the data. *Create a table in a wide format of forest types x dominants.
5. Use Forest understory species and traits >> calculate community mean and community weighted mean for SeedMass.
6. Take the file with traits and filter it to contain only values for dominants from the exercise 2.
7. Upload Penguins data by loading the library(palmerpenguins) >> check names >> calculate mean weight and standard deviation of weight for different penguin species (group_by, summarise). >> Do you remember how to visualise the weight data for each species with ggplot?
8. Stay with Penguins. In Chapter 3.4 / Bar charts, we used a ggplot to visualise how many individuals are on each island and to which species they belong. Now, try to summarise the information into a table format with islands as rows and species as names, where values will be counts of individuals.
9. *Upload the data from the Lepidoptera folder “env_MSVejnoha_transposed.csv”. First, load the env data and check its structure and assign it to the lepidoptera_env. >> Check the histogram of the cover of trees across sites >> prepare a new variable called tree_shading corresponding to low/medium/high cover of trees (select the thresholds based on the histogram) and save it in the data >> filter out the group with only low cover and assign it to lepidoptera_env_selection >>
Import the species data from the same folder “spe_matrix_MSVejnoha.csv” which contains for each site the list of taxa and corresponding counts of individuals >> change it to long format (pivot_longer) >> and, using filtering join, filter this data to lepidoptera_spe_selection, where you will have only the same sites as in the env selection.
5.7 Further reading
Joins - chapter in R for Data Science: https://r4ds.hadley.nz/joins.html
Summarise: https://dplyr.tidyverse.org/reference/summarise.html
Summarise multiple columns: https://dplyr.tidyverse.org/reference/summarise_all.html