[R] 最优参数地理探测器

## install and library the pacakge
#install.packages("GD")
library("GD")
library(tidyverse)
library(openxlsx)

# set road
setwd("C:/Users/mingxu/Desktop")
getwd()

#read data
zhizidata <- read.xlsx("zhizi.xlsx")

## 设置最佳离散化参数
## optional methods: equal, natural, quantile, geometric, sd and manual
discmethod <- c("equal","natural","quantile","geometric", "sd")
discitv <- c(3:6)

## "gdm" function
## veg and soil are categorical variables,
## bio_14 ,bio_10 ,Elev and Human are continuous variables.
zhizigdm <- gdm(Crocetin ~ bio_14 + bio_10 + Elev + Human + veg + soil,
               continuous_variable = c("bio_14", "bio_10", "Elev", "Human"),
               data = zhizidata,
               discmethod = discmethod, discitv = discitv) 
zhizigdm
plot(zhizigdm, sig = TRUE)



#因子探测
## multiple variables including continuous variables
discmethod <- c("equal","natural","quantile","geometric","sd")
discitv <- c(3:6)

#deal variables
data.continuous <- zhizidata[, c(1:5)]
odc1 <- optidisc(Crocetin ~ ., data = data.continuous, discmethod, discitv) 
data.continuous <- do.call(cbind, lapply(1:4, function(x)
  data.frame(cut(data.continuous[, 2:5][, x], unique(odc1[[x]]$itv), include.lowest = TRUE))))
# add stratified data to explanatory variables
zhizidata[, c(2:5)] <- data.continuous

g3 <- gd(Crocetin ~ ., data = zhizidata)
g3
print(g3)
plot(g3, sig = FALSE) 


#风险探测
#按变量划分的风险平均值：
## multiple variables inclusing continuous variables
rm2 <- riskmean(Crocetin ~ ., data = zhizidata)
rm2
plot(rm2)

#风险矩阵：
## multiple variables inclusing continuous variables
gr2 <- gdrisk(Crocetin ~ ., data = zhizidata)
gr2
plot(gr2)

#交互检测器
## multiple variables inclusing continuous variables
gi2 <- gdinteract(Crocetin ~ ., data = zhizidata)
gi2
plot(gi2)
#生态探测器：gdeco
## multiple variables inclusing continuous variables
gd3 <- gdeco(Crocetin ~ ., data = zhizidata)
gd3
plot(gd3)

#空间单位的大小效应比较
ndvilist <- list(bio_14, bio_10 , Elev , Human , veg , soil)
su <- c(20,30,40,50) ## sizes of spatial units
## "gdm" function
gdlist <- lapply(ndvilist, function(x){
  gdm(NDVIchange ~ Climatezone + Mining + Tempchange + GDP,
      continuous_variable = c("Tempchange", "GDP"),
      data = x, discmethod = "quantile", discitv = 6)
})
sesu(gdlist, su) ## size effects of spatial units