A Tidyverse Approach to Alpha, Beta and Gamma Diversities

Computed and Visualized Using Formulas and Using the Package vegetarian

The Formulas

For the formulas used in this file, see the following articles:

• L. Jost, “Entropy and diversity”, Oikos, vol. 113, pp. 363–375, Jan. 2006.
• L. Jost, “Partitioning diversity into independent alpha and beta components”, Ecology, vol. 88, pp. 2427–2439, Oct. 2008.

Loading Packages

library(vegetarian)

library(tidyverse)

library(latex2exp)

# about the vegetarian package
?vegetarian

Loading the Data Set vegetarian::simesants

data(simesants)

df <- simesants

Tidying the Data Set

# adding a new column `Weight`
# `Weight` is the proportion of the total population in each `Habitat`
df <- df %>% 
  rowwise() %>% 
  mutate(Count = sum(c_across(where(is.numeric)))) %>%
  ungroup()

df <- df %>% 
  mutate(Total_Count = sum(Count))

df <- df %>% 
  mutate(Weight = Count / Total_Count)

df <- df %>% 
  select(Habitat, Weight, everything(), -Count, -Total_Count)

df

## # A tibble: 5 x 24
##   Habitat Weight Acasub Aphrud Camher Campen Crecer Crelin Dolpla Dolpus Forneo
##   <fct>    <dbl>  <int>  <int>  <int>  <int>  <int>  <int>  <int>  <int>  <int>
## 1 Hemlock  0.254      0    191      0      5      0      0      0      0      0
## 2 Hardwo…  0.238      3    132      0     16      0      0      0      0      0
## 3 Hardwo…  0.131      0     24      8      2      1      1     10      7     11
## 4 Hardwo…  0.138      0     74      0      1      0      0      0      0      0
## 5 White_…  0.239      0    102      0     16      0      0      0      0      0
## # … with 13 more variables: Forsuba <int>, Lasali <int>, Lasnea <int>,
## #   Lasspe <int>, Lasumb <int>, Myrpun <int>, Myrscu <int>, Myrsmi <int>,
## #   Stebre <int>, Steimp <int>, Stesch <int>, Tapses <int>, Temlon <int>

var_names <- df %>% 
  select(-Habitat, -Weight) %>% names() 

DF <- df %>% 
  gather(all_of(var_names), key = 'Species', value = 'Count') %>% 
  filter(Count > 0) %>%
  select(Habitat, Species, everything())

DF # gathered data set, with zero counts removed

## # A tibble: 43 x 4
##    Habitat              Species Weight Count
##    <fct>                <chr>    <dbl> <int>
##  1 Hardwood             Acasub   0.238     3
##  2 Hemlock              Aphrud   0.254   191
##  3 Hardwood             Aphrud   0.238   132
##  4 Hardwood_Pine_Swamp  Aphrud   0.131    24
##  5 Hardwood_Rocky_Slope Aphrud   0.138    74
##  6 White_Pine           Aphrud   0.239   102
##  7 Hardwood_Pine_Swamp  Camher   0.131     8
##  8 Hemlock              Campen   0.254     5
##  9 Hardwood             Campen   0.238    16
## 10 Hardwood_Pine_Swamp  Campen   0.131     2
## # … with 33 more rows

# number of habitats or communities
N <- DF %>% select(Habitat) %>% unique() %>% nrow()

Gamma Diversities

Unweighted Gamma Diversities Computed Using Fromulas and Using vegetarian::d

giml_df <- DF %>% 
  group_by(Habitat) %>% 
  mutate(Habitat_Pop = sum(Count)) %>% 
  ungroup()

giml_u_m <- giml_df %>% 
  mutate(Total_Pop = sum(Count), 
         Prop = Count / Habitat_Pop, 
         Weight = 1 / N)

giml_u_m <- giml_u_m %>% 
  select(-Habitat) %>%
  group_by(Species) %>%
  mutate(Weighted_Prop = sum(Prop * Weight)) %>% 
  ungroup() %>%
  select(Species, Weighted_Prop) %>%
  unique()

giml_u_m <- giml_u_m %>% 
  summarise(Gamma_Richness = n(),
            Gamma_Shannon = exp(-sum(Weighted_Prop * log(Weighted_Prop))),
            Gamma_Greenberg = 1 / sum(Weighted_Prop ** 2))

giml_u_m

## # A tibble: 1 x 3
##   Gamma_Richness Gamma_Shannon Gamma_Greenberg
##            <int>         <dbl>           <dbl>
## 1             22          4.91            2.44

giml_u_v <- df %>% 
  summarise(Gamma_Richness = d(.[, -c(1, 2)], lev = 'gamma', q = 0), 
            Gamma_Shannon = d(.[, -c(1, 2)], lev = 'gamma', q = 1), 
            Gamma_Greenberg = d(.[, -c(1, 2)], lev = 'gamma', q = 2))

giml_u_v

## # A tibble: 1 x 3
##   Gamma_Richness Gamma_Shannon Gamma_Greenberg
##            <dbl>         <dbl>           <dbl>
## 1             22          4.91            2.44

Weighted Gamma Diversities Computed Using Fromulas and Using vegetarian::d

giml_w_m <- giml_df %>% 
  mutate(Total_Pop = sum(Count), 
         Prop = Count / Habitat_Pop, 
         Weight = Habitat_Pop / Total_Pop)

giml_w_m <- giml_w_m %>%  
  select(-Habitat) %>%
  group_by(Species) %>%
  mutate(Weighted_Prop = sum(Prop * Weight)) %>% 
  ungroup() %>%
  select(Species, Weighted_Prop) %>%
  unique()

giml_w_m <- giml_w_m %>% 
  summarise(Gamma_Richness = n(),
            Gamma_Shannon = exp(-sum(Weighted_Prop * log(Weighted_Prop))),
            Gamma_Greenberg = 1 / sum(Weighted_Prop ** 2))

giml_w_m

## # A tibble: 1 x 3
##   Gamma_Richness Gamma_Shannon Gamma_Greenberg
##            <int>         <dbl>           <dbl>
## 1             22          4.32            2.17

giml_w_v <- df %>% 
  summarise(Gamma_Richness = d(.[, -c(1, 2)], lev = 'gamma', wt = .$Weight, q = 0),
            Gamma_Shannon = d(.[, -c(1, 2)], lev = 'gamma', wt = .$Weight, q = 1),
            Gamma_Greenberg = d(.[, -c(1, 2)], lev = 'gamma', wt = .$Weight, q = 2))

giml_w_v

## # A tibble: 1 x 3
##   Gamma_Richness Gamma_Shannon Gamma_Greenberg
##            <dbl>         <dbl>           <dbl>
## 1             22          4.32            2.17

Alpha Diversities

Unweighted Alpha Diversities Computed Using Fromulas and Using vegetarian::d

alep_df <- DF %>% mutate(Total_Pop = sum(Count))

alep_u_m <- alep_df %>% 
  group_by(Habitat) %>% 
  mutate(Pop = sum(Count), 
         Prop = Count / Pop,
         Total_Prop = Count / Total_Pop,
         Weight = 1 / N)

suppressMessages(alep_u_m <- alep_u_m %>% 
                   summarise(Richness = n(), 
                             Shannon = -sum(Prop * log(Prop)),
                             Greenberg = sum(Prop ** 2)) %>%
                   ungroup() %>% 
                   unique())

alep_u_m <- alep_u_m %>% 
  summarise(Alpha_Richness = mean(Richness), 
            Alpha_Shannon = exp(mean(Shannon)), 
            Alpha_Greenberg = 1 / mean(Greenberg))

alep_u_m

## # A tibble: 1 x 3
##   Alpha_Richness Alpha_Shannon Alpha_Greenberg
##            <dbl>         <dbl>           <dbl>
## 1            8.6          3.22            2.01

alep_u_v <- df %>% 
  summarise(Alpha_Richness = d(.[, -c(1, 2)], lev = 'alpha', q = 0), 
            Alpha_Shannon = d(.[, -c(1, 2)], lev = 'alpha', q = 1), 
            Alpha_Greenberg = d(.[, -c(1, 2)], lev = 'alpha', q = 2))

alep_u_v

## # A tibble: 1 x 3
##   Alpha_Richness Alpha_Shannon Alpha_Greenberg
##            <dbl>         <dbl>           <dbl>
## 1            8.6          3.22            2.01

Weighted Alpha Diversities Computed Using Fromulas and Using vegetarian::d

alep_df <- DF %>% mutate(Total_Pop = sum(Count))

alep_w_m <- alep_df %>% 
  group_by(Habitat) %>% 
  mutate(Habitat_Pop = sum(Count), 
         Prop = Count / Habitat_Pop,
         Weight = Habitat_Pop / Total_Pop,
         Richness = n(),
         Shannon = -sum(Prop * log(Prop)),
         Greenberg = sum(Prop ** 2)) %>%
  ungroup() %>%
  select(Habitat, 
         Richness, 
         Shannon, 
         Greenberg, 
         Habitat_Pop, 
         Total_Pop, 
         Weight) %>%
  unique()
        
suppressMessages(alep_w_m <- alep_w_m %>% 
                   summarise(Alpha_Richness = mean(Richness),
                             Alpha_Shannon = exp(sum(Shannon * Weight)),
                             Alpha_Greenberg = sum(Weight ** 2) / sum(Weight ** 2 * Greenberg )))

alep_w_m

## # A tibble: 1 x 3
##   Alpha_Richness Alpha_Shannon Alpha_Greenberg
##            <dbl>         <dbl>           <dbl>
## 1            8.6          2.98            1.77

alep_w_v <- df %>% 
  summarise(Alpha_Richness = d(.[, -c(1, 2)], lev = 'alpha', wt = .$Weight, q = 0), 
            Alpha_Shannon = d(.[, -c(1, 2)], lev = 'alpha', wt = .$Weight, q = 1),
            Alpha_Greenberg = d(.[, -c(1, 2)], lev = 'alpha', wt = .$Weight, q = 2))

alep_w_v

## # A tibble: 1 x 3
##   Alpha_Richness Alpha_Shannon Alpha_Greenberg
##            <dbl>         <dbl>           <dbl>
## 1            8.6          2.98            1.77

Beta Diversities

Unweighted Beta Diversities Computed Using Fromulas and Using vegetarian::d

bet_u_m <- giml_u_m / alep_u_m 

names(bet_u_m) <- c('Beta_Richness', 'Beta_Shannon', 'Beta_Greenberg')

bet_u_v <- giml_u_v / alep_u_v

names(bet_u_v) <- c('Beta_Richness', 'Beta_Shannon', 'Beta_Greenberg')

bet_u_m

##   Beta_Richness Beta_Shannon Beta_Greenberg
## 1       2.55814     1.523506       1.212183

bet_u_v

##   Beta_Richness Beta_Shannon Beta_Greenberg
## 1       2.55814     1.523506       1.212183

# alternatively
beta_u_v <- df %>% 
  summarise(Beta_Richness = d(.[, -c(1, 2)], lev = 'beta', q = 0), 
            Beta_Shannon = d(.[, -c(1, 2)], lev = 'beta', q = 1), 
            Beta_Greenberg = d(.[, -c(1, 2)], lev = 'beta', q = 2))

beta_u_v

## # A tibble: 1 x 3
##   Beta_Richness Beta_Shannon Beta_Greenberg
##           <dbl>        <dbl>          <dbl>
## 1          2.56         1.52           1.21

Weighted Beta Diversities Computed Using Fromulas and Using vegetarian::d

bet_w_m <- giml_w_m / alep_w_m

names(bet_w_m) <- c('Beta_Richness', 'Beta_Shannon', 'Beta_Greenberg')

bet_w_v <- giml_w_v / alep_w_v

names(bet_w_v) <- c('Beta_Richness', 'Beta_Shannon', 'Beta_Greenberg')

bet_w_m

##   Beta_Richness Beta_Shannon Beta_Greenberg
## 1       2.55814     1.449162       1.227822

bet_w_v

##   Beta_Richness Beta_Shannon Beta_Greenberg
## 1       2.55814     1.449162       1.227822

# alternatively
beta_u_v <- df %>% 
  summarise(Beta_Richness = d(.[, -c(1, 2)], lev = 'beta', wt = .$Weight, q = 0), 
            Beta_Shannon = d(.[, -c(1, 2)], lev = 'beta', wt = .$Weight, q = 1), 
            Beta_Greenberg = d(.[, -c(1, 2)], lev = 'beta', wt = .$Weight, q = 2))

beta_u_v

## # A tibble: 1 x 3
##   Beta_Richness Beta_Shannon Beta_Greenberg
##           <dbl>        <dbl>          <dbl>
## 1          2.56         1.45           1.23

Visualizing Diversities

# defining two ranges: (0, 1) and (1, 10)
range_1 <- seq(0.001, 1, .01)

range_2 <- seq(1.001, 5, .01)

Visualizing Unweighted Gamma Diversity Computed Using Fromulas

qsum1_g_u_m <- NULL

giml <- giml_df %>% 
  mutate(Total_Pop = sum(Count), 
         Prop = Count / Habitat_Pop, 
         Weight = 1 / N)

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- giml %>%  
      select(-Habitat) %>%
      group_by(Species) %>%
      mutate(Weighted_Prop = sum(Prop * Weight)) %>% 
      ungroup() %>%
      select(Species, Weighted_Prop) %>%
      unique()
    
    df2 <- df2 %>% 
      summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1/ (1 - q)), 
                q = q)
    
    qsum1_g_u_m <- rbind(qsum1_g_u_m, df2)
    
}

print(qsum1_g_u_m %>% 
        ggplot(aes(x = q, y = Giml_Manual)) + 
        geom_line( color = 'blue') +
        geom_hline(yintercept = giml_u_m$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_u_m$Gamma_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color= 'yellow') + 
        labs(title = 'Unweighted Gamma Diversity Computed Using Fromulas ') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

qsum2_g_u_m <- NULL

giml <- giml_df %>% 
  mutate(Total_Pop = sum(Count), 
         Prop = Count / Habitat_Pop,
         Weight = 1 / N)

for (q in range_2) {
  
    df2 <- NULL
    
    df2 <- giml %>%  
      select(-Habitat) %>%
      group_by(Species) %>%
      mutate(Weighted_Prop = sum(Prop * Weight)) %>% 
      ungroup() %>%
      select(Species, Weighted_Prop) %>%
      unique()
    
    df2 <- df2 %>% 
      summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1 / (1 - q)),
                q = q)
    
    qsum2_g_u_m <- rbind(qsum2_g_u_m, df2)
    
}

print(qsum2_g_u_m %>% 
        ggplot(aes(x = q, y = Giml_Manual)) + 
        geom_line( color = 'blue') +
        geom_hline(yintercept = giml_u_m$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_u_m$Gamma_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color= 'yellow') + 
        labs(title = 'Unweighted Gamma Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

Visualizing Weighted Gamma Diversity Computed Using Fromulas

qsum1_g_w_m <- NULL

giml <- giml_df %>% 
  mutate(Total_Pop = sum(Count), 
         Prop = Count / Habitat_Pop, 
         Weight = Habitat_Pop / Total_Pop)

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- giml %>%  
      select(-Habitat) %>%
      group_by(Species) %>%
      mutate(Weighted_Prop = sum(Prop * Weight)) %>% 
      ungroup() %>%
      select(Species, Weighted_Prop) %>%
      unique()
    
    df2 <- df2 %>% 
      summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1/ (1 - q)), 
                q = q)
    
    qsum1_g_w_m <- rbind(qsum1_g_w_m, df2)
    
}

print(qsum1_g_w_m %>% 
        ggplot(aes(x = q, y = Giml_Manual)) + 
        geom_line( color = 'blue') +
        geom_hline(yintercept = giml_w_m$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_w_m$Gamma_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color= 'yellow') + 
        labs(title = 'Weighted Gamma Diversity Computed Using Fromulas ') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

qsum2_g_w_m <- NULL

giml <- giml_df %>% 
  mutate(Total_Pop = sum(Count), 
         Prop = Count / Habitat_Pop, 
         Weight = Habitat_Pop / Total_Pop)

for (q in range_2) {
  
    df2 <- NULL
    
    df2 <- giml %>%  
      select(-Habitat) %>%
      group_by(Species) %>%
      mutate(Weighted_Prop = sum(Prop * Weight)) %>% 
      ungroup() %>%
      select(Species, Weighted_Prop) %>%
      unique()
    
    df2 <- df2 %>% 
      summarise(Giml_Manual = sum(Weighted_Prop ** q) ** (1 / (1 - q)),
                q = q)
    
    qsum2_g_w_m <- rbind(qsum2_g_w_m, df2)
    
}

print(qsum2_g_w_m %>% 
        ggplot(aes(x = q, y = Giml_Manual)) + 
        geom_line(color = 'blue') +
        geom_hline(yintercept = giml_w_m$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_w_m$Gamma_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color= 'yellow') + 
        labs(title = 'Weighted Gamma Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

Visualizing Unweighted Gamma Diversity Computed Using vegetarian::d

qsum1_g_u_v <- NULL

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'gamma', wt = 1 / N, q = q), q = q)
    
    names(df2) <- c('Giml_Vegetarian', 'q')
    
    qsum1_g_u_v <- rbind(qsum1_g_u_v, df2)
    
}

print(qsum1_g_u_v %>% 
        ggplot(aes(x = q, y = Giml_Vegetarian)) + 
        geom_line(color = 'blue') +
        geom_hline(yintercept = giml_u_v$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_u_v$Gamma_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color= 'yellow') + 
        labs(title = 'Unweighted Gamma Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

qsum2_g_u_v <- NULL

for (q in range_2) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'gamma', wt = 1 / N, q = q), q = q)
    
    names(df2) <- c('Giml_Vegetarian', 'q')
    
    qsum2_g_u_v <- rbind(qsum2_g_u_v, df2)
    
}

print(qsum2_g_u_v %>% 
        ggplot(aes(x = q, y = Giml_Vegetarian)) + 
        geom_line(color = 'blue') +
        geom_hline(yintercept = giml_u_v$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_u_v$Gamma_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color= 'yellow') + 
        labs(title = 'Unweighted Gamma Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

Visualizing Weighted Gamma Diversity Computed Using vegetarian::d

qsum1_g_w_v <- NULL

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'gamma', wt = df$Weight, q = q), q = q)
    
    names(df2) <- c('Giml_Vegetarian', 'q')
    
    qsum1_g_w_v <- rbind(qsum1_g_w_v, df2)
    
}

print(qsum1_g_w_v %>% 
        ggplot(aes(x = q, y = Giml_Vegetarian)) + 
        geom_line(color = 'blue') +
        geom_hline(yintercept = giml_w_v$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_w_v$Gamma_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color= 'yellow') + 
        labs(title = 'Weighted Gamma Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

qsum2_g_w_v <- NULL

for (q in range_2) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'gamma', wt = df$Weight, q = q), q = q)
    
    names(df2) <- c('Giml_Vegetarian', 'q')
    
    qsum2_g_w_v <- rbind(qsum2_g_w_v, df2)
    
}

print(qsum2_g_w_v %>% 
        ggplot(aes(x = q, y = Giml_Vegetarian)) + 
        geom_line(color = 'blue') +
        geom_hline(yintercept = giml_w_v$Gamma_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = giml_w_v$Gamma_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color= 'yellow') + 
        labs(title = 'Weighted Gamma Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\gamma$ Diversity'))) 

alep <- DF %>% mutate(Total_Pop = sum(Count))

Visualizing Unweighted Alpha Diversity Computed Using Fromulas

qsum1_a_u_m <- NULL

for (q in range_1) {
  
    df2 <- NULL

    df2 <- alep %>%
      group_by(Habitat) %>% 
      mutate(Habitat_Pop = sum(Count), 
             Prop = Count / Habitat_Pop,
             Weight = 1 / N,
             Smallqsum = sum((Prop * Weight) ** q),
             q = q) %>%
      ungroup() %>%
      select(Habitat, Smallqsum, Weight, q) %>%
      unique()
    
    suppressMessages(df2 <- df2 %>% 
                       group_by(q) %>% 
                       summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), 
                                 q = q))
    
    qsum1_a_u_m <- rbind(qsum1_a_u_m, df2)
    
}

print(qsum1_a_u_m %>% 
        ggplot(aes(x = q, y = Alep_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = alep_u_m$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = alep_u_m$Alpha_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Unweighted Alpha Diversity Computed Using Fromulas ') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity'))) 

qsum2_a_u_m <- NULL

for (q in range_2) {
    
    df2 <- NULL
    
    df2 <- alep %>%
      group_by(Habitat) %>% 
      mutate(Habitat_Pop = sum(Count), 
             Prop = Count / Habitat_Pop,
             Weight = 1 / N,
             Smallqsum = sum((Prop * Weight) ** q),
             q = q) %>%
      ungroup() %>%
      select(Habitat, Smallqsum, Weight, q) %>%
      unique()
    
    suppressMessages(df2 <- df2 %>% 
                       group_by(q) %>% 
                       summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), 
                                 q = q))
    
    qsum2_a_u_m <- rbind(qsum2_a_u_m, df2)
}

print(qsum2_a_u_m %>% 
        ggplot(aes(x = q, y = Alep_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = alep_u_m$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = alep_u_m$Alpha_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Unweighted Alpha Diversity Computed Using Fromulas ') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity')))  

Visualizing Weighted Alpha Diversity Computed Using Fromulas

qsum1_a_w_m <- NULL

for (q in range_1) {
  
    df2 <- NULL

    df2 <- alep %>%
      group_by(Habitat) %>% 
      mutate(Habitat_Pop = sum(Count), 
             Prop = Count / Habitat_Pop,
             Weight = Habitat_Pop / Total_Pop,
             Smallqsum = sum((Prop * Weight) ** q),
             q = q) %>%
      ungroup() %>%
      select(Habitat, Smallqsum, Weight, q) %>%
      unique()
    
    suppressMessages(df2 <- df2 %>% 
                       group_by(q) %>% 
                       summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), 
                                 q = q))
    
    qsum1_a_w_m <- rbind(qsum1_a_w_m, df2)
    
}

print(qsum1_a_w_m %>% 
        ggplot(aes(x = q, y = Alep_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = alep_w_m$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = alep_w_m$Alpha_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Weighted Alpha Diversity Computed Using Fromulas ') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity'))) 

qsum2_a_w_m <- NULL

for (q in range_2) {
    
    df2 <- NULL
    
    df2 <- alep %>%
      group_by(Habitat) %>% 
      mutate(Habitat_Pop = sum(Count), 
             Prop = Count / Habitat_Pop,
             Weight = Habitat_Pop / Total_Pop,
             Smallqsum = sum((Prop * Weight) ** q),
             q = q) %>%
      ungroup() %>%
      select(Habitat, Smallqsum, Weight, q) %>%
      unique()
    
    suppressMessages(df2 <- df2 %>% 
                       group_by(q) %>% 
                       summarise(Alep_Manual = (sum(Smallqsum) / sum(Weight ** q)) ** (1 / (1 - q)), 
                                 q = q))
    
    qsum2_a_w_m <- rbind(qsum2_a_w_m, df2)
    
}

print(qsum2_a_w_m %>% 
        ggplot(aes(x = q, y = Alep_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = alep_w_m$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = alep_w_m$Alpha_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Weighted Alpha Diversity Computed Using Fromulas ') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity')))  

Visualizing Unweighted Alpha Diversity Computed Using vegetarian::d

qsum1_a_u_v <- NULL

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'alpha', q = q), q = q)
    
    names(df2) <- c('Alep_Vegetarian', 'q')
    
    qsum1_a_u_v <- rbind(qsum1_a_u_v, df2)
    
}

print(qsum1_a_u_v %>% 
        ggplot(aes(x = q, y = Alep_Vegetarian)) + 
        geom_line(color = 'blue') +
        geom_hline(yintercept = alep_u_v$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = alep_u_v$Alpha_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color= 'yellow') + 
        labs(title = 'Unweighted Alpha Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity'))) 

qsum2_a_u_v <- NULL

for (q in range_2) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'alpha', q = q), q = q)
    
    names(df2) <- c('Alep_Vegetarian', 'q')
    
    qsum2_a_u_v <- rbind(qsum2_a_u_v, df2)
    
}

print(qsum2_a_u_v %>% 
        ggplot(aes(x = q, y = Alep_Vegetarian)) + 
        geom_line(color = 'blue') +
        geom_hline(yintercept = alep_u_v$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color= 'green') + 
        geom_hline(yintercept = alep_u_v$Alpha_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color= 'yellow') + 
        labs(title = 'Unweighted Alpha Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity'))) 

Visualizing Weighted Alpha Diversity Computed Using `vegetarian::d``

qsum1_a_w_v <- NULL

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'alpha', wt = df$Weight, q = q), q = q)
    
    names(df2) <- c('Alep_Vegetarian', 'q')
    
    qsum1_a_w_v <- rbind(qsum1_a_w_v, df2)
    
}

print(qsum1_a_w_v %>% 
        ggplot(aes(x = q, y = Alep_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = alep_w_v$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = alep_w_v$Alpha_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Weighted Alpha Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity'))) 

qsum2_a_w_v <- NULL

for (q in range_2) {
    
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'alpha', wt = df$Weight, q = q), q = q)
    
    names(df2) <- c('Alep_Vegetarian', 'q')
    
    qsum2_a_w_v <- rbind(qsum2_a_w_v, df2)
    
}

print(qsum2_a_w_v %>% 
        ggplot(aes(x = q, y = Alep_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = alep_w_v$Alpha_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = alep_w_v$Alpha_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Weighted Alpha Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity')))  

Visualizing Unweighted Beta Diversity Computed Using Fromulas

qsum1_b_u_m <- inner_join(qsum1_g_u_m, qsum1_a_u_m, by = 'q')

qsum1_b_u_m <- qsum1_b_u_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)

print(qsum1_b_u_m %>%
        ggplot(aes(x = q, y = Bet_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_u_m$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_u_m$Beta_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Unweighted Beta Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

qsum2_b_u_m <- inner_join(qsum2_g_u_m, qsum2_a_u_m, by = 'q')

qsum2_b_u_m <- qsum2_b_u_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)

print(qsum2_b_u_m %>% 
        ggplot(aes(x = q, y = Bet_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_u_m$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_u_m$Beta_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Unweighted Beta Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

# alternatively
qsum1_b_u_v <- NULL

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'beta', q = q), q = q)
    
    names(df2) <- c('Bet_Vegetarian', 'q')
    
    qsum1_b_u_v <- rbind(qsum1_b_u_v, df2)
    
}

print(qsum1_b_u_v %>% 
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_u_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_u_v$Beta_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Weighted Alpha Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity'))) 

# alternatively
qsum2_b_u_v <- NULL

for (q in range_2) {
    
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'beta', q = q), q = q)
    
    names(df2) <- c('Bet_Vegetarian', 'q')
    
    qsum2_b_u_v <- rbind(qsum2_b_u_v, df2)
    
}

print(qsum2_b_u_v %>% 
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_u_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_u_v$Beta_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Weighted Beta Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

Visualizing Weighted Beta Diversity Computed Using Fromulas

qsum1_b_w_m <- inner_join(qsum1_g_w_m, qsum1_a_w_m, by = 'q')

qsum1_b_w_m <- qsum1_b_w_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)

print(qsum1_b_w_m %>% 
        ggplot(aes(x = q, y = Bet_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_w_m$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_w_m$Beta_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Weighted Beta Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

qsum2_b_w_m <- inner_join(qsum2_g_w_m, qsum2_a_w_m, by = 'q')

qsum2_b_w_m <- qsum2_b_w_m %>% mutate(Bet_Manual = Giml_Manual / Alep_Manual)

print(qsum2_b_w_m %>% 
        ggplot(aes(x = q, y = Bet_Manual)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_w_m$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_w_m$Beta_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Weighted Beta Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

Visualizing Unweighted Beta Diversity Computed Using vegetarian::d

qsum1_b_u_v <- inner_join(qsum1_g_u_v, qsum1_a_u_v, by = 'q')

qsum1_b_u_v <- qsum1_b_u_v %>% mutate(Bet_Vegetarian = Giml_Vegetarian / Alep_Vegetarian)

print(qsum1_b_u_v %>% 
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_u_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_u_v$Beta_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Unweighted Beta Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

qsum2_b_u_v <- inner_join(qsum2_g_u_v, qsum2_a_u_v, by = 'q')

qsum2_b_u_v <- qsum2_b_u_v %>% mutate(Bet_Vegetarian = Giml_Vegetarian / Alep_Vegetarian)

print(qsum2_b_u_v %>%
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_u_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_u_v$Beta_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Unweighted Beta Diversity Computed Using Fromulas') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

Visualizing Weighted Beta Diversity Computed Using vegetarian::d

qsum1_b_w_v <- inner_join(qsum1_g_w_v, qsum1_a_w_v, by = 'q')

qsum1_b_w_v <- qsum1_b_w_v %>% mutate(Bet_Vegetarian = Giml_Vegetarian / Alep_Vegetarian)

print(qsum1_b_w_v %>% 
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_w_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_w_v$Beta_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Weighted Beta Diversity Computed Using `vegetarian::d`') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

qsum2_b_w_v <- inner_join(qsum2_g_w_v, qsum2_a_w_v, by = 'q')

qsum2_b_w_v <- qsum2_b_w_v %>% mutate(Bet_Vegetarian = Giml_Vegetarian / Alep_Vegetarian)

print(qsum2_b_w_v %>% 
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_w_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_w_v$Beta_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Weighted Beta Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

# alternatively
qsum1_b_w_v <- NULL

for (q in range_1) {
  
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'beta', wt = df$Weight, q = q), q = q)
    
    names(df2) <- c('Bet_Vegetarian', 'q')
    
    qsum1_b_w_v <- rbind(qsum1_b_w_v, df2)
    
}

print(qsum1_b_w_v %>% 
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_w_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_w_v$Beta_Richness, color = 'yellow') +
        geom_vline(xintercept = 0, color = 'yellow') + 
        labs(title = 'Weighted Beta Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\alpha$ Diversity'))) 

# alternatively
qsum2_b_w_v <- NULL

for (q in range_2) {
    
    df2 <- NULL
    
    df2 <- data.frame(d(df[, -c(1, 2)], lev = 'beta', wt = df$Weight, q = q), q = q)
    
    names(df2) <- c('Bet_Vegetarian', 'q')
    
    qsum2_b_w_v <- rbind(qsum2_b_w_v, df2)
    
}

print(qsum2_b_w_v %>% 
        ggplot(aes(x = q, y = Bet_Vegetarian)) + 
        geom_line( color = 'blue') + 
        geom_hline(yintercept = bet_w_v$Beta_Shannon, color = 'green') +
        geom_vline(xintercept = 1, color = 'green') + 
        geom_hline(yintercept = bet_w_v$Beta_Greenberg, color = 'yellow') +
        geom_vline(xintercept = 2, color = 'yellow') + 
        labs(title = 'Weighted Beta Diversity Computed Using vegetarian::d') +
        xlab(TeX('$q$')) +
        ylab(TeX('$\\beta$ Diversity')))  

MacArthur’s Homogeneity Measure

Order 0 (Richness)

# using formulas
alep_u_m$Alpha_Richness / giml_u_m$Gamma_Richness 

## [1] 0.3909091

# using vegetarian::M.homog
M.homog(df[, -c(1, 2)], q = 0)

## [1] 0.3909091

# alternatively
alep_u_v$Alpha_Richness / giml_u_v$Gamma_Richness 

## [1] 0.3909091

Order 1 (Shannon)

# using formulas
alep_u_m$Alpha_Shannon / giml_u_m$Gamma_Shannon

## [1] 0.6563805

# using vegetarian::M.homog
M.homog(df[, -c(1, 2)])

## [1] 0.6563805

# using vegetarian
alep_u_v$Alpha_Shannon / giml_u_v$Gamma_Shannon

## [1] 0.6563805

Order 2 (Greenberg)

# using formulas
alep_u_m$Alpha_Greenberg / giml_u_m$Gamma_Greenberg

## [1] 0.8249579

# using vegetarian::M.homog
M.homog(df[, -c(1, 2)], q = 2)

## [1] 0.8249579

# alternatively
alep_u_v$Alpha_Greenberg / giml_u_v$Gamma_Greenberg

## [1] 0.8249579

Relative Homogeneity

Unweighted

# using formulas
homog_u_m <- data.frame(list((1 / bet_u_m$Beta_Richness - 1 / N) / (1 - 1 / N),
                             (1 / bet_u_m$Beta_Shannon - 1 / N) / (1 - 1 / N), 
                             (1 / bet_u_m$Beta_Greenberg - 1 / N) / (1 - 1 / N)))

names(homog_u_m) <- c('Order 0 Homogeneity', 
                      'Order 1 Homogeneity', 
                      'Order 2 Homogeneity')

homog_u_m

##   Order 0 Homogeneity Order 1 Homogeneity Order 2 Homogeneity
## 1           0.2386364           0.5704757           0.7811974

# using vegetarian::Rel.homog
Rel.homog(df[, -c(1, 2)])

## [1] 0.5704757

Weighted

# using formulas
d_1_w_m <- exp(-sum(df$Weight * log(df$Weight))) 

(1 / bet_w_m$Beta_Shannon - 1 / d_1_w_m) / (1 - 1 / d_1_w_m)

## [1] 0.6087685

(1 / bet_u_m$Beta_Shannon - 1 / d_1_w_m) / (1 - 1 / d_1_w_m)

## [1] 0.5662637

# using vegetarian::Rel.homog
Rel.homog(df[, -c(1, 2)], wt = df$Weight)

## [1] 0.5662637

# is there a bug in vegetarian::Rel.homog?

Turnover

Order 0 (Richness)

# using formulas
(bet_u_m$Beta_Richness - 1) / (N - 1)

## [1] 0.3895349

# using vegetarian::turnover
turnover(df[, -c(1, 2)], q = 0)

## [1] 0.3895349

# alternatively
(bet_u_v$Beta_Richness - 1) / (N - 1)

## [1] 0.3895349

Order 1 (Shannon)

# using formulas
(bet_u_m$Beta_Shannon - 1) / (N - 1)

## [1] 0.1308766

# using vegetarian::turnover
turnover(df[, -c(1, 2)])

## [1] 0.1308766

# alternatively
(bet_u_v$Beta_Shannon - 1) / (N - 1)

## [1] 0.1308766

Order 2 (Greenberg)

# using formulas
(bet_u_m$Beta_Greenberg - 1) / (N - 1)

## [1] 0.05304575

# using vegetarian::turnover
turnover(df[, -c(1, 2)], q = 2)

## [1] 0.05304575

# alternatively
(bet_u_v$Beta_Greenberg - 1) / (N - 1)

## [1] 0.05304575