搜索
人大经济论坛 › 附件下载

附件下载

所在主题：
文件名:  analysisScript_KerrPindarGalpern_etal_Science_2015.txt
资料下载链接地址: https://bbs.pinggu.org/a-1832508.html
附件大小:
59.28 KB   举报本内容
  1. ## ---------------------------------------------------------------
  2. ## 标明题目（CLIMATE CHANGE IMPACTS ON BUMBLEBEES CONVERGE ACROSS CONTINENTS）
  3. ## ---------------------------------------------------------------
  4. ##
  5. ##
  6. ## 指明了R平台
  7. ##
  8. ##
  9. ## 指明了作者
  10. ## 指明了问题联系人
  11. ##
  12. ##
  13. ##
  14. ## 撰写运行说明
  15. ##
  16. ## ---------------------------------------------------------------
  17. ## 数据探索阶段及预先设置
  18. ## ---------------------------------------------------------------
  19. ## 指明文档需要的数据及存放位置

  21. rawDataPath <- "G:/自媒体/2015-07-15/DryadSubmission/"
  22. ##
  23. ##
  24. ## 设定维度和温度的最大值个数，之所以取多个最大值，是因为论文中使用的最大值是
  25. ## 5个观测最大值的均值。
  26. ##
  27. numExtremeObs <- 5
  28. ##
  29. ## 移除个别物种
  30. ## 数据探索阶段移除个别特殊类型的物种是十分必要的
  31. ## 通过正则表达式查找需要移除的物种，多个物种之间用"|"隔开形如 "affinis|occidentalis".
  32. ## 一些探索性分析:
  33. ## 1.移除具有异常值的物种
  34. ## removeSpecies <- "sporadicus|borealis|pomorum|norvegicus|frigidus|bohemicus|cryptarum|affinis|flavifrons"
  35. ## 2.移除生活在墨西哥的品种
  36. # removeSpecies <- "pensylvanicus|huntii|fervidus|nevadensis|rufocinctus"
  37. ## 3.移除种群弱化受其他因素影响的品种
  38. ## removeSpecies <- "affinis|occidentalis"
  39. removeSpecies <- FALSE #正式发表的文章中可以不用移除物种
  40. ## 设置绘图颜色
  41. ## 说明颜色将被在那些配件中使用，这里颜色用于区分欧洲和北美及置信区间带
  42. plotCol <- c("#d01c8b", "#4dac26")#点线颜色用于区分欧美
  43. ciCol <- c("#f1b6da40", "#b8e18640")#设置置信区间颜色用于区分欧美
  44. ## 设置图片保存格式，如果绘图窗口出现图片，就保存到工作目录
  45. plotType <- "JPEG"
  46. ##
  47. ##
  48. ## ---------------------------------------------------------------
  49. ## 生产文章所用的图表及分析模型
  50. ## ---------------------------------------------------------------
  51. ##
  52. ## 程序保证上一步获得R运行结果保存在工作目录下，及上一步的记过要保存在内存中，
  53. ## 中途关闭或上一步无果可能导致下一步报错终止
  54. ##
  55. ##
  56. ##
  57. ##
  58. ## 加载R包，如果缺少该包就安装
  59. ##
  60. ##
  61. if (!(suppressWarnings(require("ape")))) {
  62. install.packages("ape")
  63. require(ape)
  64. }
  65. ## suppressWarnings函数用于忽略警告（warn），如果报错suppressWarnings返回的值为FALSE

  67. if (!(suppressWarnings(require("nlme")))) {
  68. install.packages("nlme")
  69. require(nlme)
  70. }
  71. ##
  72. ##第一步，确定两个时期品种的维度、温度、海拔均值
  73. #
  74. doStep1 <- TRUE
  75. if (doStep1) {


  78. ## 载入数据
  79. if (!exists("rawData")) rawData <- read.csv(paste0(rawDataPath, "data_bombclimate.csv"), stringsAsFactors=FALSE)
  80. speciesTable <- unique(rawData[, c("species", "continent")])
  81. samplingPeriods <- unique(rawData[, c("periodFrom", "periodTo")])
  82. #exists用于检查对象是否在内存中，kmNorthEquator指从赤道向北的距离,maxPeriodAnnualMeanT年平均最高气温

  84. rawRangeShift <- list()
  85. speciesKmNorthEquator <- list()
  86. speciesElevation <- list()
  87. cnt <- 0

  89. ## 确定每个品种的均值
  90. for (iSpecies in 1:nrow(speciesTable)) {

  92. cat(speciesTable[iSpecies, "species"], "\n")#输出抽提的品种名称

  94. ## 将这一品种所有时期的数据从原始数据（rawData）中抽取出来
  95. thisSpecies <- rawData[rawData$species==speciesTable[iSpecies, "species"] &
  96. rawData$continent==speciesTable[iSpecies, "continent"], ]

  98. ## 计算这一品种的均值（不分时期）
  99. speciesKmNorthEquator[[iSpecies]] <- mean(thisSpecies$kmNorthEquator)
  100. speciesElevation[[iSpecies]] <- mean(thisSpecies$elevation)

  102. ## 分时期计算均值
  103. for (iPeriod in 1:nrow(samplingPeriods)) {

  105. ## 提取这个品种每个时期的原始数据
  106. thisSpeciesPeriod <- rawData[rawData$species==speciesTable[iSpecies, "species"] &
  107. rawData$continent==speciesTable[iSpecies, "continent"] &
  108. rawData$periodFrom==samplingPeriods[iPeriod, "periodFrom"] &
  109. rawData$periodTo==samplingPeriods[iPeriod, "periodTo"], ]

  111. thisSpeciesPeriodN <- nrow(thisSpeciesPeriod)

  113. ## 对该品种发现的地点去重，因为一个地点可能有多次记录
  114. uniqueLoc <- !duplicated(thisSpeciesPeriod$spPerLoc)

  116. ## 计算以距离定义的寒冷地区的均值和炎热地区的均值，寒冷地区指距离迟到最远的5个地区，炎热地区反之，个数可以自己定义
  117. thisLat <- thisSpeciesPeriod$kmNorthEquator[uniqueLoc][order(thisSpeciesPeriod$kmNorthEquator[uniqueLoc], decreasing=TRUE)]
  118. thisLatCold <- mean(thisLat[1:numExtremeObs], na.rm=TRUE)
  119. thisLatHot <- mean(rev(thisLat)[1:numExtremeObs], na.rm=TRUE)#rev为反转

  121. ## 计算以温度定义的炎热地区及寒冷地区
  122. thisThermalMax <- thisSpeciesPeriod$maxPeriodAnnualMeanT[uniqueLoc]
  123. thisThermalMax <- thisThermalMax[order(thisThermalMax, decreasing=TRUE)]
  124. thisThermalMin <- thisSpeciesPeriod$minPeriodAnnualMeanT[uniqueLoc]
  125. thisThermalMin <- thisThermalMin[order(thisThermalMin)]
  126. thisThermalCold <- mean(thisThermalMin[1:numExtremeObs], na.rm=TRUE)
  127. thisThermalHot <- mean(thisThermalMax[1:numExtremeObs], na.rm=TRUE)

  129. ## 计算该品种该时期的平均海拔
  130. thisMeanElevation <- mean(thisSpeciesPeriod$elevation, na.rm=TRUE)

  132. ## 计算该品种该时段的维度中间点，中间点即均值
  133. thisMeanKmNorthEquator <- mean(thisSpeciesPeriod$kmNorthEquator[uniqueLoc], na.rm=TRUE)

  135. ## 将以上计算的探索结果整理成data.frame
  136. cnt <- cnt + 1
  137. rawRangeShift[[cnt]] <-data.frame(species=speciesTable[iSpecies, "species"],
  138. continent=speciesTable[iSpecies, "continent"],
  139. periodYears=paste(samplingPeriods[iPeriod, c("periodFrom", "periodTo")], collapse="to"),
  140. period=paste0("t", iPeriod),
  141. latitudeHot=thisLatHot, latitudeCold=thisLatCold,
  142. thermalHot=thisThermalHot, thermalCold=thisThermalCold,
  143. elevation=thisMeanElevation,
  144. latitude=thisMeanKmNorthEquator,
  145. N=thisSpeciesPeriodN)

  147. }

  149. }

  151. rawRangeShift <- do.call(rbind, rawRangeShift)
  152. speciesKmNorthEquator <- do.call(rbind, speciesKmNorthEquator)
  153. speciesElevation <- do.call(rbind, speciesElevation)



  157. ## 对品种—时期均值数据框整形，以适应针对品种建模过程
  158. rangeShift <- reshape(rawRangeShift, idvar=c("species", "continent"), timevar="period", direction="wide", sep="_")
  159. ## 这个数据整形的函数用的比较古老，dcast更加高大上
  160. ## 计算每品种—时期 deltas (deltas指t2到t1; t3到t1;t4到t1变化)并将其加入数据框
  161. ##
  162. vars <- apply(expand.grid(c("latitudeHot", "latitudeCold", "thermalHot", "thermalCold", "elevation", "N"), c("t2", "t3", "t4")), 1, function(x) paste(x, collapse="_"))
  163. vars <- data.frame(tX=vars, t1=gsub("t2|t3|t4", "t1", vars), delta=paste0("delta_", vars), stringsAsFactors=FALSE)
  164. ## 这里要说一下expand.grid函数，他将所有变量下的因子进行组合形成一个data.frame,比如说第一个变量
  165. ## 第二变量有三个因子，他们组成的记过就是一个就是一个有两个变量6行数据组成的数据框

  167. delta <- lapply(1:nrow(vars), function(x) {
  168. y <- data.frame(rangeShift[, vars[x, "tX"]] - rangeShift[, vars[x, "t1"]])
  169. names(y) <- vars[x, "delta"]
  170. return(y)
  171. })
  172. rangeShift <- data.frame(rangeShift, do.call(cbind, delta))


  175. ## 将品种的维度均值和海拔均值加入数据框
  176. rangeShift <- data.frame(rangeShift, meanLatitude=speciesKmNorthEquator, meanElevation=speciesElevation)

  178. ## 将杀虫剂数据加入数据框(仅限于北美品种有杀虫剂使用数据)
  179. pesticides <- read.csv(paste0(rawDataPath, "speciesUSA_pesticides.csv"))
  180. speciesIndex <- match(pesticides$species, rangeShift$species)

  182. rangeShift$usaTotalPest_t3 <- NA
  183. rangeShift$usaTotalPest_t3[speciesIndex] <- pesticides$totalPest_t3
  184. rangeShift$usaTotalPest_t4 <- NA
  185. rangeShift$usaTotalPest_t4[speciesIndex] <- pesticides$totalPest_t4
  186. rangeShift$usaNeonics_t3 <- NA
  187. rangeShift$usaNeonics_t3[speciesIndex] <- pesticides$neonics_t3
  188. rangeShift$usaNeonics_t4 <- NA
  189. rangeShift$usaNeonics_t4[speciesIndex] <- pesticides$neonics_t4
  190. ## 作者对R的函数实在不熟悉，以上那么多句可以用一句完成
  191. ## rangeShift <- join(rangeShift, pesticides)

  193. ## 添加土地使用情况的数据（包括北美和欧洲）
  194. landuse <- read.csv(paste0(rawDataPath, "species_landuse.csv"))
  195. rangeShift <- data.frame(rangeShift, landuse[, -1])

  197. ## If removeSpecies is set, remove those species from the dataset that match the
  198. ## regular expression given above
  199. if (!is.logical(removeSpecies)) {
  200. cat("Removing:", gsub("\\|", ", ", removeSpecies), "\n")#这里gsub省略了要移除的品种，可以在“\\|”添加
  201. rangeShift <- rangeShift[!grepl(removeSpecies, rangeShift$species), ]
  202. }

  204. }

  206. ## ---------------------------------------------------------------
  207. ## 第二步
  208. ## ---------------------------------------------------------------
  209. ## 为delta维度、温度数据构建最小二乘法回归模型，
  210. ## 然后选取AIC最低的OLS模型构建phylogenetic generalized least-squares (PGLS)回归模型，用来检测
  211. ## 两个或多个相关变量在整个种群树上是否表现出显著性，例如对于一个物种，
  212. ## 如果两个变量表现出强烈相关，那与该物种比较近亲的物种也可能表现出相同的相关性
复制代码

上述是Science近日发表的文章的R脚本特点如下：
1.规范，包括书写、命名、格式规范
2.全程一个脚本，完成所有分析任务，完整性和整合性明显
3.作者用的R包很少，但这也限制了作者代码的间接性。
4.细节和代码串联的方法值得学习
5.个别小错误，比如将cbind写成了c
更多数据和关于phylogenetic generalized least-squares (PGLS)回归模型的脚本请关注我们
相关数据及代码：http://pan.baidu.com/s/1i3rNTul密码：微信索取
关于我们，关注理性与文艺，用数据创作内容性的精致阅读，这里是数据分析挖掘人员与文艺青年的集结地，不做培训，不做鼓吹，只踏踏实实的做一个又一个数据驱动的文章，并设计机器人减轻数据分析的负担，无论你感兴趣还是想参与都可以关注，请加微信公众号大音如霜


    熟悉论坛请点击新手指南
下载说明
1、论坛支持迅雷和网际快车等p2p多线程软件下载，请在上面选择下载通道单击右健下载即可。
2、论坛会定期自动批量更新下载地址,所以请不要浪费时间盗链论坛资源,盗链地址会很快失效。
3、本站为非盈利性质的学术交流网站,鼓励和保护原创作品，拒绝未经版权人许可的上传行为。本站如接到版权人发出的合格侵权通知，将积极的采取必要措施；同时，本站也将在技术手段和能力范围内，履行版权保护的注意义务。
(如有侵权，欢迎举报)
二维码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

京ICP备16021002号-2 京B2-20170662号 京公网安备 11010802022788号 论坛法律顾问：王进律师 知识产权保护声明   免责及隐私声明

GMT+8, 2025-12-31 20:56