2025-12-27
#Importing the library vegan- library(vegan)
## Loading required package: permute
library(iNEXT) library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ── ## ✔ dplyr 1.1.4 ✔ readr 2.1.6 ## ✔ forcats 1.0.1 ✔ stringr 1.6.0 ## ✔ ggplot2 4.0.1 ✔ tibble 3.3.0 ## ✔ lubridate 1.9.4 ✔ tidyr 1.3.2 ## ✔ purrr 1.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ## ✖ dplyr::filter() masks stats::filter() ## ✖ dplyr::lag() masks stats::lag() ## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2) library(reshape2)
## ## Attaching package: 'reshape2' ## ## The following object is masked from 'package:tidyr': ## ## smiths
set.seed(2025)
#stimulating abundance data for 3 habitats and 10 species
species<- paste0("Sp", 1:10)
habitats<- c( "Forest", "Grasslands", "Wetland")
data<- data.frame(
Habitat = rep(habitats, each = 10),
Species = rep(species, time = 3),
Abundance = c(
rpois(10, lambda = 30), #Forest
rpois(10, lambda = 15), #Grassland
rpois(10, lambda = 20) #Wetland
)
)
#Coverting to community matrix (species as columns)
comm<- reshape2::dcast(data, Habitat ~ Species, value.var = "Abundance")
row.names(comm)<- comm$Habitat
comm<- comm[, -1]
comm
## Sp1 Sp10 Sp2 Sp3 Sp4 Sp5 Sp6 Sp7 Sp8 Sp9 ## Forest 33 33 30 34 36 32 29 32 29 28 ## Grasslands 13 20 8 20 11 14 14 16 9 20 ## Wetland 24 22 23 14 14 14 18 23 15 16
# Species Richness specnumber(comm)
## Forest Grasslands Wetland ## 10 10 10
#calculating Shannon and Simpson Indices diversity(comm, index = "shannon") # Shannon diversity (H')
## Forest Grasslands Wetland ## 2.299665 2.258591 2.278567
diversity(comm, index = "simpson") # Simpson diversity (D)
## Forest Grasslands Wetland ## 0.8994152 0.8914150 0.8951596
#calculating Hill Number hill_numbers<- data.frame( Habitat = rownames(comm), D0 = specnumber(comm), D1 = exp(diversity(comm,"shannon")), D2 = 1/ diversity(comm, "simpson") ) hill_numbers
## Habitat D0 D1 D2 ## Forest Forest 10 9.970843 1.111834 ## Grasslands Grasslands 10 9.569596 1.121812 ## Wetland Wetland 10 9.762679 1.117119
#probability of interspecific ENCOUNTER (PIE) MEASURES THE PROBABILITY
calc_pie<- function(x) {
N <- sum(x)
p <- x/N
(N / (N-1)*(1-sum(p**2)))
}
PIE <- apply(comm,1, calc_pie)
PIE
## Forest Grasslands Wetland ## 0.9022704 0.8976054 0.9000781
#Spicies Accumulation and Rarefaction curves #Spicies Accumulation curves help visualize effort vs richness #Spicies Accumulation sac <- specaccum(comm, method = "exact")
## Warning in cor(x > 0): the standard deviation is zero
#Plotting plot(sac, ci.type = "polygon", col= "blue")
# iNEXT-based visualization for interpolation / Extrapolation inext_res <- iNEXT(comm, q= c(0,1, 2), datatype = "abundance") ggiNEXT(inext_res, type = 1, facet.var = "Order.q") + theme_minimal() + labs(title = "Hill Numbers (q=0, 1, 2): Interpolation and Extrapolation")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the iNEXT package.
## Please report the issue at <https://github.com/AnneChao/iNEXT/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The shape palette can deal with a maximum of 6 discrete values because more ## than 6 becomes difficult to discriminate ## ℹ you have requested 10 values. Consider specifying shapes manually if you need ## that many of them.
## Warning: Removed 12 rows containing missing values or values outside the scale range ## (`geom_point()`).
#we stimulate capture histories for 20 individuals over 5 sampling occasion set.seed(2025) n_individuals <- 20 n_occasions <- 5 capture_prob <- 0.4 capture_hist <- matrix(rbinom(n_individuals * n_occasions, 1, capture_prob),nrow = n_individuals, n_occasions) capture_hist
## [,1] [,2] [,3] [,4] [,5] ## [1,] 1 0 0 0 1 ## [2,] 0 1 1 1 0 ## [3,] 0 0 1 1 1 ## [4,] 0 1 1 0 0 ## [5,] 1 0 1 0 0 ## [6,] 0 1 1 1 0 ## [7,] 1 1 1 1 0 ## [8,] 0 0 0 0 0 ## [9,] 1 1 0 0 0 ## [10,] 0 0 0 0 0 ## [11,] 0 0 0 0 1 ## [12,] 1 0 1 1 1 ## [13,] 1 0 0 1 0 ## [14,] 0 0 0 0 1 ## [15,] 0 1 0 0 0 ## [16,] 1 0 0 0 1 ## [17,] 0 0 0 1 1 ## [18,] 0 0 0 1 1 ## [19,] 1 1 1 1 1 ## [20,] 0 1 1 0 0
n1 <- sum(capture_hist[,1]==1) n2 <- sum(capture_hist[,2]==1) m2 <- sum(capture_hist[,1]==1 & capture_hist[,2]==1) n2
## [1] 8
x<- c(1,2,3,4) sort(x, decreasing =TRUE)
## [1] 4 3 2 1