R语言：地理探测器分析以及ggplot绘制结果箱型图和热图

最新推荐文章于 2024-06-21 14:27:25 发布

曦嘻喜喜

最新推荐文章于 2024-06-21 14:27:25 发布

阅读量1.4k

点赞数 1

文章标签： r语言数学建模 java

本文链接：https://blog.csdn.net/weixin_45626690/article/details/131341843

版权

该文利用R语言进行环境数据分析，首先对数据进行预处理，包括读取、筛选、分类，然后运用地理探测器进行单因子分析和交互作用分析，得出统计量q和p值。通过箱型图展示了每个因子的中位数差异，并按中位数排序，最后绘制热图来呈现因子间的交互效果。整个过程揭示了环境因子对研究目标的影响模式和强度变化。

摘要由CSDN通过智能技术生成

单因子结果绘制为箱型图

箱型图绘制：箱体颜色按中位数值填充；同时从大到小排列；五个子图拼图，共用colorbar

library(geodetector)
library(xlsx)
library(dplyr)
library(reshape2)
library(ggplot2)
library(tidyverse)
library(Kendall)
library(classInt)
library(RColorBrewer)

setwd("E:/enviromental_data_analysis/spatiotem/paper/datasorting/data_analysis2")
current_dir <- getwd()
###5个阶段，逐一筛选各阶段文件完成每一幅子图
filenames<-list.files(current_dir, pattern = "p1\\.xlsx$")

###列名制备colnames
year<-c(1984:2013)
year<-as.character(year)
yy<-c(0)
for (y in year) {
  yy<-cbind(yy,paste0("y",y))
}

###各阶段变量使用classIntervals()函数进行自然断点分类转为类型量，为地理探测器分析做准备
lm<-c()
p1_class<-c()
for (f in filenames) {
  mydata<-read.xlsx(f,sheetName = "Sheet1")
  mydata<-mydata[,2:31]
  for (col in 1:30) {
    breaks <- classIntervals(mydata[,col], n = 6, style = "jenks")
    # 获取分类断点
    break_points <- unique(breaks$brks)
    # 使用cut()函数将数据分组
    cut_data <- cut(mydata[,col], breaks = break_points, include.lowest = TRUE,labels = FALSE)
    p1_class<-cbind(p1_class,cut_data)
    lm<-rbind(lm,paste0(strsplit(f, "_")[[1]][1],yy[col+1]))
  }
  
}
###整理表格内容
p1_class<-as.data.frame(p1_class)
names(p1_class)[1:510]<-lm

#####相应阶段因变量GPP读取
gpp<-read.xlsx("E:/enviromental_data_analysis/spatiotem/paper/datasorting/gpp/gpp_p1.xlsx",sheetName = "Sheet1")
f_sprm <- merge(gpp, fer, by.x = "pac", by.y = "pac",all = TRUE)
f_sprm<-na.omit(f_sprm)
phase1<-cbind(p1_class,f_sprm[,3:32])
write.xlsx(phase1,"E:/enviromental_data_analysis/spatiotem/paper/datasorting/geodetector2/phase1_gpp_class.xlsx")

#地理探测器-因子分析
###30年的数据，逐年循环，各变量在上述步骤中按列拼接的，所以numbers循环得到每一年的所有变量的列
###具体factor——detector参数可以help查看，案例很详细
qz<-c()
pz<-c()
for (y in 1:30) {
  numbers <- seq(y, 480+y, by = 30)
  geo<-factor_detector(510+y,numbers,phase1)
  qs<-c()
  pv<-c()
  for (l in 1:length(geo)) {
    q<-geo[[l]]["q-statistic"][1,1]
    qs<-cbind(qs,q)
    p<-geo[[l]]["p-value"][1,1]
    pv<-cbind(pv,p)
  }
  qz<-rbind(qz,qs)
  pz<-rbind(pz,pv)
}

qz<-as.data.frame(qz)
colnames(qz)<-c("ARID","SDII","Fer","Irr","Mac","Par","Pre","RM","scPDSI","Tem_10",
                "Tem_30","Tem_avg","Tem_dif","Tem_max","Tem_min","VPD","Wind")

##MK trend
library(trend)
sen<-c()
pv<-c()
for (i in 1:ncol(qz)) {
  column_data <- qz[, i]
  senn<-sens.slope(ts(column_data,start = 1984))
  sen <- cbind(sen,senn[["estimates"]])
  pv<-cbind(pv,senn[["p.value"]])
}
qtrend<-rbind(sen,pv)
write.xlsx(qtrend,"E:/enviromental_data_analysis/spatiotem/paper/datasorting/geodetector2/phase1_gpp_qtrend.xlsx")


##短数据转换为长数据
grouped_data <- qz %>% 
  gather(key = "variable", value = "value")
###分组后计算中位数
grouped_data <- grouped_data %>%
  group_by(variable) %>%
  mutate(median_value = median(value))
mean_values <- grouped_data %>% group_by(variable) %>% summarise(mean_value = mean(value))
median_values <- grouped_data %>% group_by(variable) %>% summarise(median_value = median(value))
###为分面绘制准备
grouped_data$phase<-c("Phase_1")
###boxplot
# 使用ggplot绘制箱型图，并根据均值大小指定渐变色填充
theme_set(theme_bw(base_size = 16, base_family = "serif"))
phase1<-ggplot(grouped_data, aes(x = reorder(variable, -median_value), y = value,fill = median_value)) +
  geom_boxplot()+
  scale_fill_gradientn(colors = rev(brewer.pal(11, "RdBu")),
                       aesthetics = "fill",
                       name = "q(median)",
                       guide = guide_colorbar(title.position = "top",
                                              title.hjust = 0.5,
                                              label.position = "right")) +
  guides(fill = guide_colorbar(title.position = "top",
                               title.hjust = 0.5,
                               label.position = "right")) +
  facet_wrap(~phase)+
  labs(title ="", 
       x=NULL, 
       y="q")+
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_text(angle=30, hjust = 1,size = 12, color = "black",family = 'serif',face = "bold"),
        axis.text.y= element_text(size = 12, color = "black",family = 'serif',face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(family = 'serif',size = 14, face = "bold"),
        legend.text = element_text(family = 'serif',size = 14, color = "black"),
        legend.title = element_text(family = 'serif',size = 14, color = "black"),
        strip.text = element_text(family = 'serif',size=14, lineheight=1,face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5))

###其它阶段此处省略
#拼图
library(ggpubr)
plot_grid <- ggarrange(phase1, phase2, phase3, phase4, phase5, nrow = 3, ncol = 2, common.legend = TRUE, legend = "right")
# Display the combined plot
print(plot_grid)

ggsave(filename = "E:/enviromental_data_analysis/spatiotem/paper/figure2/geo_individual.tiff", 
       width = 16,
       height = 12, 
       units = "in",
       dpi = 600)

双因子交互结果绘制为热图

热图绘制：结果标在图上

###上述数据制备过程一致，此处不再重复，最终分析的数据为制备好的phase1（包括自变量和因变量）
###30年的数据，逐年循环求取，以中位数值绘图分析
qz<-c()
for (y in 1:30) {
  numbers <- seq(y, 480+y, by = 30)
  geo<-interaction_detector(510+y,numbers,phase1)
  geo<-as.data.frame(geo)
  v<-geo[,3]
  v<-as.numeric(v)
  qz<-cbind(qz,v)
}

write.xlsx(qz,file=paste0("E:/enviromental_data_analysis/spatiotem/paper/datasorting/geodetector2/","phase1_inter_q30_gpp.xlsx"))
qzm<-apply(qz, 1, median)
# qzme<-apply(qz, 1, mean)
geo<-geo[,-3]
geo$value<-qzm
geo_inter<-geo
names(geo_inter)[3]<-"values"
geo_inter$values<-as.numeric(geo_inter$values)
# geo_inter_2<-dcast(geo_inter, V1 ~ V2, value.var = "values",fun.aggregate = mean)
write.xlsx(geo_inter,file=paste0("E:/enviromental_data_analysis/spatiotem/paper/datasorting/geodetector2/","phase1_inter_q30median_gpp.xlsx"))
##MK trend
library(trend)
sen<-c()
pv<-c()
for (i in 1:ncol(qz)) {
  column_data <- qz[, i]
  senn<-sens.slope(ts(column_data,start = 1984))
  sen <- cbind(sen,senn[["estimates"]])
  pv<-cbind(pv,senn[["p.value"]])
}
qtrend<-rbind(sen,pv)
write.xlsx(qtrend,"E:/enviromental_data_analysis/spatiotem/paper/datasorting/geodetector2/phase1_inter_qtrend_gpp.xlsx")

###绘热图
theme_set(theme_bw(base_size = 14, base_family = "serif"))

x_labels <- c("ARID","SDII","Fer","Irr","Mac","Par","Pre","RM","scPDSI","Tem_10",
              "Tem_30","Tem_avg","Tem_dif","Tem_max","Tem_min","VPD","Wind"
)
y_labels <- c("ARID","SDII","Fer","Irr","Mac","Par","Pre","RM","scPDSI","Tem_10",
              "Tem_30","Tem_avg","Tem_dif","Tem_max","Tem_min","VPD","Wind")

geo_inter$label <- sprintf("%.2f", geo_inter$values)

phase1<-ggplot(geo_inter, aes(x = V1, y = V2, fill = values)) + 
  geom_tile() + 
  scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(11, "RdBu"))) +
  # 修改刻度大小和颜色
  theme(axis.text.x = element_text(angle=30, hjust = 1,size = 12, color = "black"),
        axis.text.y = element_text(size = 12, color = "black"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.text = element_text(size = 12, color = "black"),
        legend.title = element_text(size = 12, color = "black")) +
  # 调整大小
  theme(plot.title = element_text(size = 20, hjust = 0.5, face = "bold"))+
  scale_x_discrete(labels = x_labels) +
  scale_y_discrete(labels = y_labels)+
  geom_text(aes(label = label), color = "black", size = 3, hjust = 0.5, vjust = 0.5)##标签结果数据

###拼图方式同上