如何在循环中遍历所有3站点组合计算生物多样性gamma指数

实现代码 #

library(tidyverse)
set.seed(5)

# 导入示例数据集
bats_data <- structure(list(Site = c("MOZ_Codzo", "KEN_Sokoke1", "NGA_Gotel", 
                        "ZAF_Pafuri", "COD_Kitwit", "GIN_Diecke", "MWI_Mulanje2", "NER_W", 
                        "GAB_Maboumine", "COD_Kitwit", "SWZ_Lubombo", "GIN_Ziama2", "ZAF_Kruger1", 
                        "KEN_Kakamega", "NGA_Delta", "CMR_Dja", "ZAF_Pafuri", "KEN_Kakamega", 
                        "ZAF_Wakefield", "SWZ_Lubombo", "BDI_Bururi", "ZAF_Soutpansberg", 
                        "ZAF_Kruger2", "GIN_Simandou2", "GIN_Simandou1", "CMR_Maroua", 
                        "CMR_Cameroon2", "KEN_Kakamega", "SWZ_IYSIS", "CMR_Cameroon4", 
                        "ZWE_Sengwa", "TZA_Coastal2", "GIN_Ziama1", "ZAF_Kruger3", "SLE_Seli", 
                        "CIV_Tai", "GIN_Nimba1", "GIN_Simandou1", "MWI_Zomba1", "NER_W", 
                        "GIN_Ziama1", "ZAF_Soutpansberg", "LBR_Nimba3", "CIV_Tai", "TZA_Coastal2", 
                        "KEN_Sokoke1", "GNB_Cacheu", "MWI_Liwonde", "CIV_Tai", "SOM_Jubba", 
                        "ZAF_DeHoop", "ZAF_Durban", "CIV_Tai", "SWZ_IYSIS", "TZA_Usambara1", 
                        "GIN_Ziama1", "CMR_Mpem", "GIN_Nimba2", "LBR_Nimba1", "SWZ_Lubombo", 
                        "CMR_Cameroon1", "LBR_Nimba1", "COD_Epulu", "GNB_Cantanhez", 
                        "NGA_Gotel", "ZWE_Kanyati", "COD_Epulu", "CIV_Comoe", "SLE_Gola3", 
                        "GIN_Simandou1", "NGA_Edo3", "ETH_Beletta", "MWI_Liwonde", "SLE_Gola1", 
                        "BWA_Xugana", "KEN_Masalani", "CMR_Cameroon3", "ZAF_Kruger3", 
                        "ETH_Godare", "CMR_Mpem", "SLE_Gola3", "GIN_Ziama1", "ZAF_Sudwala", 
                        "CAF_Floris", "CIV_Comoe", "MLI_Bla", "TZA_Usambara1", "SLE_Seli", 
                        "GIN_Simandou2", "ZWE_Sengwa", "LBR_Nimba1", "ZWE_Kanyati", "SLE_Gola3", 
                        "ZWE_Sengwa", "SLE_Seli", "CIV_Tai", "SLE_Seli", "MWI_Mulanje1", 
                        "ZAF_Algeria", "ZAF_Kruger3"), Gen_sp = c("Mops.condylurus", 
                                                                  "Nycteris.thebaica", "Myonycteris.angolensis", "Nycteris.thebaica", 
                                                                  "Myopterus.whitleyi", "Neoromicia.tenuipinnis", "Myotis.tricolor", 
                                                                  "Scotophilus.leucogaster", "Hipposideros.ruber", "Mops.condylurus", 
                                                                  "Miniopterus.natalensis", "Mops.brachypterus", "Mops.condylurus", 
                                                                  "Glauconycteris.humeralis", "Hypsignathus.monstrosus", "Rhinolophus.alcyone", 
                                                                  "Eptesicus.hottentotus", "Taphozous.mauritianus", "Myotis.tricolor", 
                                                                  "Neoromicia.capensis", "Rhinolophus.clivosus", "Miniopterus.natalensis", 
                                                                  "Neoromicia.nana", "Epomops.buettikoferi", "Rousettus.aegyptiacus", 
                                                                  "Scotophilus.livingstonii", "Myonycteris.torquata", "Myotis.bocagii", 
                                                                  "Miniopterus.natalensis", "Megaloglossus.woermanni", "Neoromicia.nana", 
                                                                  "Nycteris.hispida", "Hipposideros.ruber", "Chaerephon.pumilus", 
                                                                  "Rhinolophus.denti", "Nycteris.grandis", "Nanonycteris.veldkampii", 
                                                                  "Kerivoula.lanosa", "Rhinolophus.blasii", "Hipposideros.ruber", 
                                                                  "Hipposideros.marisae", "Rhinolophus.clivosus", "Nycteris.arge", 
                                                                  "Chaerephon.aloysiisabaudiae", "Neoromicia.capensis", "Pipistrellus.rueppellii", 
                                                                  "Neoromicia.guineensis", "Myotis.welwitschii", "Mops.spurrelli", 
                                                                  "Epomophorus.wahlbergi", "Miniopterus.natalensis", "Myotis.tricolor", 
                                                                  "Casinycteris.ophiodon", "Scotophilus.viridis", "Rhinolophus.deckenii", 
                                                                  "Mops.condylurus", "Hipposideros.ruber", "Epomops.buettikoferi", 
                                                                  "Hypsignathus.monstrosus", "Taphozous.mauritianus", "Epomops.franqueti", 
                                                                  "Neoromicia.somalica", "Megaloglossus.woermanni", "Hypsignathus.monstrosus", 
                                                                  "Chaerephon.ansorgei", "Neoromicia.capensis", "Casinycteris.argynnis", 
                                                                  "Nycteris.intermedia", "Myonycteris.leptodon", "Kerivoula.cuprosa", 
                                                                  "Eidolon.helvum", "Myonycteris.torquata", "Neoromicia.nana", 
                                                                  "Hypsignathus.monstrosus", "Chaerephon.pumilus", "Miniopterus.arenarius", 
                                                                  "Myonycteris.angolensis", "Nycticeinops.schlieffeni", "Neoromicia.guineensis", 
                                                                  "Rousettus.aegyptiacus", "Glauconycteris.variegata", "Miniopterus.villiersi", 
                                                                  "Neoromicia.capensis", "Neoromicia.nana", "Hipposideros.ruber", 
                                                                  "Scotophilus.leucogaster", "Scotophilus.andrewreborii", "Mops.nanulus", 
                                                                  "Eidolon.helvum", "Nycteris.thebaica", "Micropteropus.pusillus", 
                                                                  "Mops.condylurus", "Hipposideros.abae", "Neoromicia.capensis", 
                                                                  "Eidolon.helvum", "Macronycteris.gigas", "Hipposideros.jonesi", 
                                                                  "Rhinolophus.fumigatus", "Cistugo.lesueuri", "Scotophilus.dinganii"
                        ), Numbers = c(3L, 14L, 1L, 2L, 2L, 1L, 10L, 3L, 10L, 10L, 6L, 
                                       2L, 9L, 1L, 4L, 3L, 1L, 2L, 1L, 2L, 1L, 2L, 4L, 1L, 5L, 227L, 
                                       4L, 33L, 14L, 1L, 14L, 1L, 1L, 6L, 1L, 4L, 29L, 1L, 187L, 2L, 
                                       1L, 6L, 1L, 2L, 3L, 5L, 1L, 1L, 12L, 39L, 2L, 1L, 80L, 5L, 1L, 
                                       1L, 8L, 6L, 1L, 5L, 64L, 1L, 34L, 1L, 1L, 36L, 1L, 1L, 14L, 5L, 
                                       7L, 1L, 1L, 5L, 12L, 29L, 31L, 21L, 1L, 8L, 1L, 1L, 2L, 2L, 19L, 
                                       1L, 1L, 6L, 3L, 1L, 1L, 17L, 1L, 53L, 4L, 1L, 6L, 4L, 5L, 32L
                        ), presence = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        1)), row.names = c(NA, -100L), class = "data.frame")

# 筛选出现在超过3个站点的物种
bats3 <- bats_data %>% count(Gen_sp) %>% filter(n > 3)
bat_names <- unique(bats3$Gen_sp)

# 初始化结果存储数据框
theta.output.bats.3 <- data.frame(Species = character(),
                                  combo_id = integer(),
                                  gamma = integer(),
                                  stringsAsFactors = FALSE)

# 遍历每个目标物种
for (sp in bat_names) {
  # 提取该物种出现的所有站点
  sp_sites <- unique(bats_data$Site[bats_data$Gen_sp == sp])
  # 生成所有3站组合
  all_combo <- combn(sp_sites, 3, simplify = FALSE)
  
  # 遍历每个组合计算gamma
  for (c_id in 1:length(all_combo)) {
    combo_sites <- all_combo[[c_id]]
    # 筛选对应站点的所有物种记录
    sites1 <- subset(bats_data, Site %in% combo_sites)
    gamma <- nrow(distinct(sites1, Gen_sp))
    # 写入结果
    theta.output.bats.3 <- rbind(theta.output.bats.3,
                                 data.frame(Species = sp,
                                            combo_id = c_id,
                                            gamma = gamma))
  }
}

# 查看结果
head(theta.output.bats.3)

核心改动说明 #

• 替换了原来的随机抽样逻辑，改用combn()生成目标物种所有出现站点的3站不重复组合
• 新增了组合层级的循环，每个组合单独计算gamma多样性，实现每个物种对应所有组合都有独立结果
• 输出改为结构化数据框，列依次为物种名、组合编号、对应gamma值，方便后续统计分析

内容的提问来源于stack exchange，提问作者adkane