Load package

pacman::p_load(tidyverse,magrittr,DT)

Create function

b.plot <- function(df){
    ggplot(df, aes(x=simu))+facet_wrap( ~ carte)+
        geom_density(fill="#159858", colour="#159858",alpha=.5)+
        geom_vline(data = d2, aes(xintercept = lower),colour="blue") +
        geom_text(data = d2, aes(x= lower,y=1000), label="2.5%", colour="blue") + 
        geom_vline(data = d2, aes(xintercept = upper),colour="blue")+
        geom_text(data = d2, aes(x= upper,y=1000), label="97.5%", colour="blue") + 
        geom_text(data = d2, aes(x= obs,y=900), label="Obs", colour="red", size=5) + 
        geom_segment(data =d2, mapping=aes(x=obs, y=800, xend=obs, yend=0), arrow = arrow(), size=1, colour="red")+
        scale_x_continuous(limits = c(-0.005, 0.035))+
        theme_bw()
}

Analysis

Load results from MATLAB

setwd("D:/GIS/doc_Elo/R_figure")
dat.vd <- read.csv("results.csv") #vd =  Vincent Deblauwe
datatable(dat.vd,class="compact")

1.

Create dat2 from data.vd and add 0.001 to P-value to make P-value > 0

a <- 0.001
dat2 <- dat.vd
dat2$p_elv <- dat2$p_ele + a
dat2$p_slope <- dat2$p_slope + a
dat2$p_rel.el <- dat2$p_rel.el + a
datatable(dat2,class="compact")

2.

Added the sign of the correlation
Create 2 new columns p_Slope_Sign and p_RelEl_Sign
if value r < 0 ==> P-value * -1

dat2$p_Slope_Sign <- vector(nrow(dat2),mode="numeric")
dat2$p_RelEl_Sign <- vector(nrow(dat2),mode="numeric")

for (i in 1:nrow(dat2)) {
    ifelse(dat2[i,"r_rel.el"] < 0, x <- -1, x <- 1)
    dat2[i,"p_RelEl_Sign"] <- dat2[i,"p_rel.el"] * x
    ifelse(dat2[i,"r_slope"] < 0, x <- -1, x <- 1)
    dat2[i,"p_Slope_Sign"] <- dat2[i,"p_slope"] * x
}
datatable(dat2,class="compact")

3.

Convert P-value to visualisation by log(p_value)
Create 2 new columns p_Slope_Graph and p_RelEl_Graph

dat2$p_Slope_Graph <- vector(nrow(dat2),mode="numeric")
dat2$p_RelEl_Graph <- vector(nrow(dat2),mode="numeric")
for (i in 1:nrow(dat2)) {
    ifelse(dat2[i,"r_rel.el"] < 0, x <- -1, x <- 1)
    dat2[i,"p_RelEl_Graph"] <- abs(log(dat2[i,"p_rel.el"])) * x
    ifelse(dat2[i,"r_slope"] < 0, x <- -1, x <- 1)
    dat2[i,"p_Slope_Graph"] <- abs(log(dat2[i,"p_slope"])) * x
}
datatable(dat2,class="compact")

Plot with 2 columns log scale p_Slope_Graph and p_RelEl_Graph

my.color <-c("blue","red")
plot(dat2$p_RelEl_Graph,dat2$p_Slope_Graph,pch=c(17,19),col=my.color,
     xlim= c(-4,2), ylim=c(-8,4),xlab="Relative elevation", ylab="Slope")
text(dat2$p_RelEl_Graph,dat2$p_Slope_Graph, labels = dat2$Species,pos=1,col=my.color)
abline(v=0,h=0,col="gray")

We can see C. obtusa presence in lowland, valley more than C. sciadophylla

Analysis with data simulation 999 times and observer value

Get observer value from data.vd

obs.data<-dat.vd %>% dplyr::select(starts_with("r_")) %>% t() %>% as.data.frame()
colnames(obs.data) <- dat.vd$Species
obs.data[order(rownames(obs.data)),]
##          C_obtusa C_sciadophylla
## r_ele      0.0805         0.0892
## r_rel.el  -0.0788        -0.0119
## r_slope   -0.0728        -0.1615
datatable(obs.data,class="compact")

For C. obtusa

Load 999 simulations values of C. obtusa

obtusa_simu <- read.csv("obtusa_simulation.csv") %>% gather(key=carte,value=simu)

Test two-tail at 2.5 % and 97.5 %

d2 <- obtusa_simu %>%
    group_by(carte) %>%
    summarize(lower = quantile(simu, probs = .025),
              upper = quantile(simu, probs = .975)) %>% 
    mutate(obs=obs.data$C_obtusa^2,
           test_2_tail = case_when(obs>=upper~"+",
                                   obs<=lower~"-",
                                   obs>lower&obs<upper~"n"))

datatable(d2,class="compact")

Plot

(g1 <- b.plot(obtusa_simu))

For C. sciadophylla

Load 999 simulations values of C. sciadophylla

sciado_simu <- read.csv("sciado_simulation.csv") %>% gather(key=carte,value=simu)

Test two-tail at 2.5 % and 97.5 %

d2 <- sciado_simu %>%
    group_by(carte) %>%
    summarize(lower = quantile(simu, probs = .025),
              upper = quantile(simu, probs = .975)) %>% 
    mutate(obs=obs.data$C_sciadophylla^2,
           test_2_tail = case_when(obs>=upper~"+",
                                   obs<=lower~"-",
                                   obs>lower&obs<upper~"n"))
datatable(d2,class="compact")

Plot

(g2 <- b.plot(sciado_simu))

library(ggpubr)
ggarrange(g1,g2,
          labels = c("C. obtusa", "C. sciado"),
          nrow=2)

There is a positive relationship between C. obtusa and the slope. And positive relationship between C. sciadophylla and relative elevation