发帖
雷达卡
- ## ---------------------------------------------------------------
- ## 标明题目（CLIMATE CHANGE IMPACTS ON BUMBLEBEES CONVERGE ACROSS CONTINENTS）
- ## ---------------------------------------------------------------
- ##
- ##
- ## 指明了R平台
- ##
- ##
- ## 指明了作者
- ## 指明了问题联系人
- ##
- ##
- ##
- ## 撰写运行说明
- ##
- ## ---------------------------------------------------------------
- ## 数据探索阶段及预先设置
- ## ---------------------------------------------------------------
- ## 指明文档需要的数据及存放位置
- rawDataPath <- "G:/自媒体/2015-07-15/DryadSubmission/"
- ##
- ##
- ## 设定维度和温度的最大值个数，之所以取多个最大值，是因为论文中使用的最大值是
- ## 5个观测最大值的均值。
- ##
- numExtremeObs <- 5
- ##
- ## 移除个别物种
- ## 数据探索阶段移除个别特殊类型的物种是十分必要的
- ## 通过正则表达式查找需要移除的物种，多个物种之间用"|"隔开形如 "affinis|occidentalis".
- ## 一些探索性分析:
- ## 1.移除具有异常值的物种
- ## removeSpecies <- "sporadicus|borealis|pomorum|norvegicus|frigidus|bohemicus|cryptarum|affinis|flavifrons"
- ## 2.移除生活在墨西哥的品种
- # removeSpecies <- "pensylvanicus|huntii|fervidus|nevadensis|rufocinctus"
- ## 3.移除种群弱化受其他因素影响的品种
- ## removeSpecies <- "affinis|occidentalis"
- removeSpecies <- FALSE #正式发表的文章中可以不用移除物种
- ## 设置绘图颜色
- ## 说明颜色将被在那些配件中使用，这里颜色用于区分欧洲和北美及置信区间带
- plotCol <- c("#d01c8b", "#4dac26")#点线颜色用于区分欧美
- ciCol <- c("#f1b6da40", "#b8e18640")#设置置信区间颜色用于区分欧美
- ## 设置图片保存格式，如果绘图窗口出现图片，就保存到工作目录
- plotType <- "JPEG"
- ##
- ##
- ## ---------------------------------------------------------------
- ## 生产文章所用的图表及分析模型
- ## ---------------------------------------------------------------
- ##
- ## 程序保证上一步获得R运行结果保存在工作目录下，及上一步的记过要保存在内存中，
- ## 中途关闭或上一步无果可能导致下一步报错终止
- ##
- ##
- ##
- ##
- ## 加载R包，如果缺少该包就安装
- ##
- ##
- if (!(suppressWarnings(require("ape")))) {
- install.packages("ape")
- require(ape)
- }
- ## suppressWarnings函数用于忽略警告（warn），如果报错suppressWarnings返回的值为FALSE
- if (!(suppressWarnings(require("nlme")))) {
- install.packages("nlme")
- require(nlme)
- }
- ##
- ##第一步，确定两个时期品种的维度、温度、海拔均值
- #
- doStep1 <- TRUE
- if (doStep1) {
-
-
- ## 载入数据
- if (!exists("rawData")) rawData <- read.csv(paste0(rawDataPath, "data_bombclimate.csv"), stringsAsFactors=FALSE)
- speciesTable <- unique(rawData[, c("species", "continent")])
- samplingPeriods <- unique(rawData[, c("periodFrom", "periodTo")])
- #exists用于检查对象是否在内存中，kmNorthEquator指从赤道向北的距离,maxPeriodAnnualMeanT年平均最高气温
-
- rawRangeShift <- list()
- speciesKmNorthEquator <- list()
- speciesElevation <- list()
- cnt <- 0
-
- ## 确定每个品种的均值
- for (iSpecies in 1:nrow(speciesTable)) {
-
- cat(speciesTable[iSpecies, "species"], "\n")#输出抽提的品种名称
-
- ## 将这一品种所有时期的数据从原始数据（rawData）中抽取出来
- thisSpecies <- rawData[rawData$species==speciesTable[iSpecies, "species"] &
- rawData$continent==speciesTable[iSpecies, "continent"], ]
-
- ## 计算这一品种的均值（不分时期）
- speciesKmNorthEquator[[iSpecies]] <- mean(thisSpecies$kmNorthEquator)
- speciesElevation[[iSpecies]] <- mean(thisSpecies$elevation)
-
- ## 分时期计算均值
- for (iPeriod in 1:nrow(samplingPeriods)) {
-
- ## 提取这个品种每个时期的原始数据
- thisSpeciesPeriod <- rawData[rawData$species==speciesTable[iSpecies, "species"] &
- rawData$continent==speciesTable[iSpecies, "continent"] &
- rawData$periodFrom==samplingPeriods[iPeriod, "periodFrom"] &
- rawData$periodTo==samplingPeriods[iPeriod, "periodTo"], ]
-
- thisSpeciesPeriodN <- nrow(thisSpeciesPeriod)
-
- ## 对该品种发现的地点去重，因为一个地点可能有多次记录
- uniqueLoc <- !duplicated(thisSpeciesPeriod$spPerLoc)
-
- ## 计算以距离定义的寒冷地区的均值和炎热地区的均值，寒冷地区指距离迟到最远的5个地区，炎热地区反之，个数可以自己定义
- thisLat <- thisSpeciesPeriod$kmNorthEquator[uniqueLoc][order(thisSpeciesPeriod$kmNorthEquator[uniqueLoc], decreasing=TRUE)]
- thisLatCold <- mean(thisLat[1:numExtremeObs], na.rm=TRUE)
- thisLatHot <- mean(rev(thisLat)[1:numExtremeObs], na.rm=TRUE)#rev为反转
-
- ## 计算以温度定义的炎热地区及寒冷地区
- thisThermalMax <- thisSpeciesPeriod$maxPeriodAnnualMeanT[uniqueLoc]
- thisThermalMax <- thisThermalMax[order(thisThermalMax, decreasing=TRUE)]
- thisThermalMin <- thisSpeciesPeriod$minPeriodAnnualMeanT[uniqueLoc]
- thisThermalMin <- thisThermalMin[order(thisThermalMin)]
- thisThermalCold <- mean(thisThermalMin[1:numExtremeObs], na.rm=TRUE)
- thisThermalHot <- mean(thisThermalMax[1:numExtremeObs], na.rm=TRUE)
-
- ## 计算该品种该时期的平均海拔
- thisMeanElevation <- mean(thisSpeciesPeriod$elevation, na.rm=TRUE)
-
- ## 计算该品种该时段的维度中间点，中间点即均值
- thisMeanKmNorthEquator <- mean(thisSpeciesPeriod$kmNorthEquator[uniqueLoc], na.rm=TRUE)
-
- ## 将以上计算的探索结果整理成data.frame
- cnt <- cnt + 1
- rawRangeShift[[cnt]] <- data.frame(species=speciesTable[iSpecies, "species"],
- continent=speciesTable[iSpecies, "continent"],
- periodYears=paste(samplingPeriods[iPeriod, c("periodFrom", "periodTo")], collapse="to"),
- period=paste0("t", iPeriod),
- latitudeHot=thisLatHot, latitudeCold=thisLatCold,
- thermalHot=thisThermalHot, thermalCold=thisThermalCold,
- elevation=thisMeanElevation,
- latitude=thisMeanKmNorthEquator,
- N=thisSpeciesPeriodN)
-
- }
-
- }
-
- rawRangeShift <- do.call(rbind, rawRangeShift)
- speciesKmNorthEquator <- do.call(rbind, speciesKmNorthEquator)
- speciesElevation <- do.call(rbind, speciesElevation)
-
-
-
- ## 对品种—时期均值数据框整形，以适应针对品种建模过程
- rangeShift <- reshape(rawRangeShift, idvar=c("species", "continent"), timevar="period", direction="wide", sep="_")
- ## 这个数据整形的函数用的比较古老，dcast更加高大上
- ## 计算每品种—时期 deltas (deltas指t2到t1; t3到t1;t4到t1变化)并将其加入数据框
- ##
- vars <- apply(expand.grid(c("latitudeHot", "latitudeCold", "thermalHot", "thermalCold", "elevation", "N"), c("t2", "t3", "t4")), 1, function(x) paste(x, collapse="_"))
- vars <- data.frame(tX=vars, t1=gsub("t2|t3|t4", "t1", vars), delta=paste0("delta_", vars), stringsAsFactors=FALSE)
- ## 这里要说一下expand.grid函数，他将所有变量下的因子进行组合形成一个data.frame,比如说第一个变量
- ## 第二变量有三个因子，他们组成的记过就是一个就是一个有两个变量6行数据组成的数据框
-
- delta <- lapply(1:nrow(vars), function(x) {
- y <- data.frame(rangeShift[, vars[x, "tX"]] - rangeShift[, vars[x, "t1"]])
- names(y) <- vars[x, "delta"]
- return(y)
- })
- rangeShift <- data.frame(rangeShift, do.call(cbind, delta))
-
-
- ## 将品种的维度均值和海拔均值加入数据框
- rangeShift <- data.frame(rangeShift, meanLatitude=speciesKmNorthEquator, meanElevation=speciesElevation)
-
- ## 将杀虫剂数据加入数据框(仅限于北美品种有杀虫剂使用数据)
- pesticides <- read.csv(paste0(rawDataPath, "speciesUSA_pesticides.csv"))
- speciesIndex <- match(pesticides$species, rangeShift$species)
-
- rangeShift$usaTotalPest_t3 <- NA
- rangeShift$usaTotalPest_t3[speciesIndex] <- pesticides$totalPest_t3
- rangeShift$usaTotalPest_t4 <- NA
- rangeShift$usaTotalPest_t4[speciesIndex] <- pesticides$totalPest_t4
- rangeShift$usaNeonics_t3 <- NA
- rangeShift$usaNeonics_t3[speciesIndex] <- pesticides$neonics_t3
- rangeShift$usaNeonics_t4 <- NA
- rangeShift$usaNeonics_t4[speciesIndex] <- pesticides$neonics_t4
- ## 作者对R的函数实在不熟悉，以上那么多句可以用一句完成
- ## rangeShift <- join(rangeShift, pesticides)
-
- ## 添加土地使用情况的数据（包括北美和欧洲）
- landuse <- read.csv(paste0(rawDataPath, "species_landuse.csv"))
- rangeShift <- data.frame(rangeShift, landuse[, -1])
-
- ## If removeSpecies is set, remove those species from the dataset that match the
- ## regular expression given above
- if (!is.logical(removeSpecies)) {
- cat("Removing:", gsub("\\|", ", ", removeSpecies), "\n")#这里gsub省略了要移除的品种，可以在“\\|”添加
- rangeShift <- rangeShift[!grepl(removeSpecies, rangeShift$species), ]
- }
-
- }
- ## ---------------------------------------------------------------
- ## 第二步
- ## ---------------------------------------------------------------
- ## 为delta维度、温度数据构建最小二乘法回归模型，
- ## 然后选取AIC最低的OLS模型构建phylogenetic generalized least-squares (PGLS)回归模型，用来检测
- ## 两个或多个相关变量在整个种群树上是否表现出显著性，例如对于一个物种，
- ## 如果两个变量表现出强烈相关，那与该物种比较近亲的物种也可能表现出相同的相关性
上述是Science近日发表的文章的R脚本特点如下：
1.规范，包括书写、命名、格式规范
2.全程一个脚本，完成所有分析任务，完整性和整合性明显
3.作者用的R包很少，但这也限制了作者代码的间接性。
4.细节和代码串联的方法值得学习
5.个别小错误，比如将cbind写成了c
更多数据和关于phylogenetic generalized least-squares (PGLS)回归模型的脚本请关注我们
相关数据及代码：http://pan.baidu.com/s/1i3rNTul 密码：微信索取
关于我们，关注理性与文艺，用数据创作内容性的精致阅读，这里是数据分析挖掘人员与文艺青年的集结地，不做培训，不做鼓吹，只踏踏实实的做一个又一个数据驱动的文章，并设计机器人减轻数据分析的负担，无论你感兴趣还是想参与都可以关注，请加微信公众号大音如霜
analysisScript_KerrPindarGalpern_etal_Science_2015.txt
(59.28 KB)
提升卡
置顶卡
沉默卡
变色卡
抢沙发
千斤顶
显身卡
京公网安备 11010802022788号
