search

JohnsonHsieh / iNEXT

R package for interpolation and extrapolation
https://JohnsonHsieh.github.com/iNEXT
57 stars 26 forks source link

Question on observed number of species in ChaoRichness() #53

Open elegiacally opened 4 years ago

elegiacally commented 4 years ago

Good day everyone!

I've a question about ChaoRichness and how it operates:

The dataset comes from a species inventory project covering several taxa at one site. I did species estimation using ChaoRichness for each taxon as well as all taxa. When looking at the ChaoRichness values gotten from all, the number of species observed don't tally (sum of species numbers observed from individual taxon =/= observed no. of species from all taxa). Does anyone know if it's an issue with my code / or data?

Here's my code:

Any insight would be greatly appreciated! 

#Data herE: 

sa1<- structure(list(Date = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 5L, 5L, 5L, 6L, 6L, 6L), .Label = c("11/3/20", "22-Oct-19", 
"23-Oct-19", "24-Oct-19", "30/10/19", "4-Mar-20", "5/3/20"), class = "factor"), 
    Time = c(859L, 859L, 859L, 859L, 903L, 903L, 907L, 913L, 
    913L, 919L, 919L, 919L, 919L, 926L, 932L, 935L, 935L, 935L, 
    1301L, 1302L, 1302L, 1307L, 1311L, 1325L, 1324L, 1326L, 1327L, 
    1327L, 1333L, 1333L, 1333L, 900L, 900L, 900L, 900L, 902L, 
    902L, 902L, 905L, 905L, 908L, 908L, 910L, 910L, 910L, 921L, 
    921L, 921L, 921L, 923L, 926L, 927L, 929L, 929L, 929L, 931L, 
    1953L, 2005L, 2016L, 1957L, 2004L, 2008L), Scientific_name = structure(c(16L, 
    29L, 45L, 14L, 12L, 36L, 25L, 8L, 16L, 1L, 14L, 13L, 1L, 
    9L, 29L, 39L, 2L, 42L, 11L, 10L, 15L, 23L, 15L, 7L, 5L, 23L, 
    15L, 34L, 18L, 44L, 32L, 10L, 38L, 10L, 29L, 38L, 37L, 26L, 
    9L, 36L, 3L, 17L, 9L, 36L, 29L, 30L, 6L, 4L, 36L, 38L, 12L, 
    9L, 41L, 20L, 2L, 45L, 31L, 31L, 22L, 35L, 43L, 43L), .Label = c(" Pycnonotus plumosus", 
    "Acridotheres javanicus", "Aegithina tiphia", "Aethopyga siparaja", 
    "Anthreptes malacensis", "Aplonis panayensis", "Burara harisa consobrina", 
    "Cacatua goffiniana", "Callosciurus notatus", "Catopsilia pyranthe pyranthe", 
    "Catopsilia scylla cornelia", "Cinnyris jugularis", "Copsychus malabaricus", 
    "Copsychus saularis", "Delias hyparete metarete", "Dicaeum cruentatum", 
    "Dicrurus paradiseus", "Eurema sp", "Eutropis multifasciata", 
    "Gorsachius melanolophus", "Hasora sp", "Hemidactylus frenatus", 
    "Idea leuconoe clara", "Lalage nigra", "Larvivora cyane", 
    "Macronus gularis", "Merops philippinus", "Naja sumatrana", 
    "Oriolus chinensis", "Orthotomus atrogularis", "Otus lempiji", 
    "Pachliopta aristolochiae asteris", "Pavo cristatus", "Phalanta phalantha phalantha", 
    "Pitta moluccensis", "Pycnonotus goiavier", "Pycnonotus plumosus", 
    "Pycnonotus zeylanicus", "Spilopelia chinensis", "Terpsiphone affinis", 
    "Todiramphus chloris", "Troides helena cerberus", "Unidentified Fruit Bat sp", 
    "Zizula hylax pygmaea", "Zosterops simplex"), class = "factor"), 
    Common_Name = structure(c(32L, 2L, 39L, 26L, 22L, 44L, 33L, 
    41L, 32L, 23L, 26L, 43L, 23L, 30L, 2L, 36L, 18L, 8L, 25L, 
    10L, 27L, 40L, 27L, 24L, 6L, 40L, 27L, 19L, 15L, 31L, 11L, 
    10L, 37L, 10L, 2L, 37L, 23L, 29L, 30L, 44L, 9L, 16L, 30L, 
    44L, 2L, 13L, 1L, 12L, 44L, 37L, 22L, 30L, 7L, 20L, 18L, 
    39L, 38L, 38L, 35L, 4L, 42L, 42L), .Label = c("Asian Glossy Starling", 
    "Black-naped Oriole", "Blue-tailed bee-eater", "Blue-winged Pitta", 
    "Blyth's Paradise Flycatcher", "Brown-throated Sunbird", 
    "Collared Kingfisher", "Common Birdwing", "Common Iora", 
    "Common Palmfly", "Common Rose", "Crimson Sunbird", "Dark-necked Tailorbird", 
    "Equatorial Spitting Cobra", "Grass Yellow sp", "Greater Racket-tailed Drongo", 
    "Indian peafowl", "Javan Myna", "Leopard", "Malayan Night Heron", 
    "Many-lined Sun Skink", "Olive-backed Sunbird", "Olive-winged Bulbul", 
    "Orange Awlet", "Orange Emigrant", "Oriental Magpie-Robin", 
    "Painted Jezebel", "Pied triller", "Pin-striped Tit-babbler", 
    "Plantain Squirrel", "Pygmy Grass Blue", "Scarlet-backed Flowerpecker", 
    "Siberian Blue Robin", "Skipper sp", "Spiny-tailed House Gecko", 
    "Spotted Dove", "Straw-headed Bulbul", "Sunda Scops Owl", 
    "Swinhoe's White-eye", "Taiwan Tree Nymph", "Tanimbar Corella", 
    "Unidentified Fruit Bat sp", "White-rumped Shama", "Yellow-vented Bulbul"
    ), class = "factor"), Taxon = structure(c(1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 
    1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 1L, 3L, 3L), .Label = c("Bird", 
    "Butterfly", "Mammal", "Reptile"), class = "factor"), Quantity = c(1L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 7L, 1L, 
    1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 3L, 1L, 2L, 1L, 1L, 
    3L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 
    1L), Survey_method = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Incidental", 
    "Transect"), class = "factor"), Sampling_Point = structure(c(1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 6L, 6L, 
    6L, 6L, 6L, 6L, 6L, 4L, 4L, 4L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
    6L, 6L, 6L, 6L, 5L, 5L, 5L, 5L, 5L, 4L, 4L, 4L, 4L, 3L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 4L, 5L, 4L, 4L, 
    3L), .Label = c("SA_01", "SA_02", "SA_03", "SA_04", "SA_05", 
    "SA_06"), class = "factor")), row.names = c(2L, 3L, 4L, 5L, 
6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 22L, 24L, 25L, 27L, 28L, 29L, 35L, 36L, 37L, 38L, 39L, 40L, 
41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 
54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 
67L, 68L, 69L, 70L, 71L, 72L), class = "data.frame")

total.matrix <- xtabs(Quantity~Scientific_name+Sampling_Point, data=sa1)
total.df <- as.data.frame.matrix(total.matrix)
(total.df1 <- ifelse(total.df>0, 1, 0)) #change to incidence matrix
total.df2<- list(total.df1)
(se.total<-ChaoRichness(total.df2, datatype="incidence_raw"))

#PER TAXON

#Split by taxon
(sa.pertaxon <- split(sa1, sa1$Taxon))
names(sa.pertaxon) #all the taxa recorded; Mammal, Bird, Butterfly, Reptile

#Generate vectors 
(sa2 <-lapply(sa.pertaxon, function(x) xtabs(Quantity ~ Scientific_name + Sampling_Point, droplevels(x))))
(sa3 <- lapply(sa2, function(x) ifelse(x>0, 1, 0))) #change to incidence matrix
sa4 <- lapply(sa3, function (x) list(x))

#chao richness
se.results <- lapply(sa4, function(x) {
  tryCatch(ChaoRichness(x, datatype="incidence_raw"),
           error = function(e) e) 
})

(se.ok <- !sapply(se.results, inherits, "error"))  #determine which ran alright 

#   Bird Butterfly    Mammal   Reptile 
#   TRUE      TRUE      TRUE     TRUE 
#note that reptile got a warning message, but it still worked 

se.all<-se.results[se.ok]
(se.all2 <-do.call("rbind",se.all))

#combine in one table
row.names(se.total) <- "All Taxa"
(se.table.all<- (rbind(se.all2,se.total)) )