---
title: "Figures"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(root.dir = here::here()) 
```


## Helper functions
```{r}
addA = function(col, alpha = 0.25) apply(sapply(col, col2rgb)/255, 2, function(x) rgb(x[1], x[2], x[3], alpha=alpha)) 
mean_conf = function(sites, mat, col = "red", alpha = 0.1, spar = 0.10, convergence = NULL, indices = NULL, hmscObject = TRUE, pch = 15) {
  #sites2 = sites[complete.cases(mat)]
  #mat = mat[complete.cases(mat),]
  m = apply(mat, 1, function(x) mean(x, na.rm = TRUE))
  sd = apply(mat, 1, function(x) sd(x, na.rm = TRUE)/sqrt(sum(!is.na(x))))
  sites2 = sites[complete.cases(m)]
  sd = sd[complete.cases(m)]

  m = m[complete.cases(m)]

  upper = smooth.spline(y = m + sd, x = sites2, spar = spar)$y
  lower = smooth.spline(y = m - sd, x = sites2, spar = spar)$y
  polygon(c(sites2, rev(sites2)), c(upper, rev(lower)), border = NA, col =addA(col, alpha))
  sp = smooth.spline(y = m, x = sites2, spar = spar)
  lines(sp, col = col, lwd = 2.0, pch = pch, type = "b")
  
  if(!is.null(convergence)) {
    indices = indices[indices %in% which(complete.cases(convergence$result_corr_acc[,1]), TRUE)]
    for(k in 1:length(indices)){
      
      if(!hmscObject) {
         betas = as.vector(sapply(convergence$post[[indices[k]]], function(p) 1.2 <abind::abind(p$psrf.beta,along = 1L)[,1]))
        lambdas = as.vector(sapply(convergence$post[[indices[k]]], function(p) 1.2 < p$psrf.gamma[,1]))
      }else {
         betas = as.vector(sapply(convergence$post[[indices[k]]], function(p) 1.2 < p$psrf.beta[,1]))
        lambdas = unlist(sapply(convergence$post[[indices[k]]], function(p) 1.2 < gelman.diag(p$post$Lambda[[1]],multivariate = FALSE)$psrf[,1]))
      }
      
      if(any(betas,lambdas)) points(sp$x[k], sp$y[k], col = col, pch = 4,cex = 1.2)
    }
  }
}

```

```{r}
deg2rad <- function(deg) {(deg * pi) / (180)}
rad2deg <- function(rad) {(rad * 180) / (pi)}
ff = function(x){(x-min(x))/(max(x)-min(x))}

curve_text = function(pos = 0, label = "", lineSeq = 5.0, reverse = FALSE,middle = FALSE,extend = 1.1,...){
  # inspired by plotrix package
  chars = (strsplit(label,split = "")[[1]])
  char_lens = strwidth(chars)*extend
  char_angles = char_lens / lineSeq
  changrang = range(char_angles)
  char_angles[char_angles < changrang[2]/2] = changrang[2]/2
  
  if(middle & reverse) pos = pos - rad2deg(sum(char_angles)/2)
  if(middle & !reverse) pos = pos + rad2deg(sum(char_angles)/2)
  
  if(reverse) {
    angles = c(deg2rad(pos), deg2rad(pos)+cumsum(char_angles)[-length(chars)])
    angles = angles + char_angles/2
  } else {
    angles = c(deg2rad(pos), deg2rad(pos)-cumsum(char_angles)[-length(chars)])
    angles = angles - char_angles/2
  }
  
  for(i in 1:length(chars)) text(label = chars[i], 
                                 x = cos((angles[i]))*(lineSeq), 
                                 srt = rad2deg(angles[i]) - 90+ 180*reverse,
                                 y = sin((angles[i]))*(lineSeq), 
                                 xpd = NA, adj = c(0.5, 0.5),...)
  return(max(angles))
  
}


add_curve = function(p1 = coords[1,], p2 = coords[3,], n = 10, spar = 0.7, col = "black", species = TRUE, lineSeq = 5.0, lwd = 1.0) {
 xxs1 = cos(deg2rad(p1[3]))* seq(0, lineSeq, length.out = n)
 xxs2 = cos(deg2rad(p2[3]))* seq(0, lineSeq, length.out = n)
 yys1 = sin(deg2rad(p1[3]))* seq(0, lineSeq, length.out = n)
 yys2 = sin(deg2rad(p2[3]))* seq(0, lineSeq, length.out = n)
 x = c(rev(xxs1), xxs2[-1])
 y = c(rev(yys1), yys2[-1])
 m = (p1[2] - p2[2])/(p1[1] - p2[1])
 a = rad2deg(atan(m))
 a = -(a+180)
 alpha = deg2rad(a)
 alpha2 = deg2rad(-a)
 rot = matrix(c(cos((alpha)), -sin((alpha)), sin((alpha)), cos((alpha))),2,2)
 rot2 = matrix(c(cos((alpha2)), -sin((alpha2)), sin((alpha2)), cos((alpha2))),2,2)
 tt = cbind(x,y) %*% rot
 sp = smooth.spline(tt[,1], tt[,2],spar = spar,df = 6, w = c(10.0, rep(0.1,nrow(tt)-2), 10.0))
 tt2 = cbind(sp$x, sp$y)
 b = tt2 %*% rot2
 lines(b[,1], b[,2], col = col, lwd = lwd)
 
 x1 = c(cos(deg2rad(p1[3]))*(lineSeq+0.1), cos(deg2rad(p1[3]))*(lineSeq+0.3))
 x2 = c(cos(deg2rad(p2[3]))*(lineSeq+0.1), cos(deg2rad(p2[3]))*(lineSeq+0.3))
 y1 = c(sin(deg2rad(p1[3]))* (lineSeq+0.1), sin(deg2rad(p1[3]))* (lineSeq+0.3))
 y2 = c(sin(deg2rad(p2[3]))* (lineSeq+0.1), sin(deg2rad(p2[3]))* (lineSeq+0.3))
 if(species){
 segments(x0 = x1[1], x1 = x1[2], y0 = y1[1], y1 = y1[2], col = "darkgrey")
 segments(x0 = x2[1], x1 = x2[2], y0 = y2[1], y1 = y2[2],  col = "darkgrey")
 }
}


add_legend = function(cols = RColorBrewer::brewer.pal(11,"Spectral"), range = c(-1,1), lineSeq = 5.0, angles = c(110, 70)){
 angles = seq(angles[1], angles[2], length.out = length(cols)+1)
 for(i in 2:length(angles)){
  xx1 = (lineSeq+0.4)*cos( seq(deg2rad(angles[i-1]),deg2rad(angles[i]) ,length.out=50) )
  xx2 = (lineSeq+0.7)*cos( seq(deg2rad(angles[i-1]),deg2rad(angles[i]) ,length.out=50) )
  yy1 = (lineSeq+0.4)*sin( seq(deg2rad(angles[i-1]),deg2rad(angles[i]) ,length.out=50)  )
  yy2 = (lineSeq+0.7)*sin( seq(deg2rad(angles[i-1]),deg2rad(angles[i]) ,length.out=50)  )
  polygon(c(xx1, rev(xx2)), c(yy1, rev(yy2)),border = NA, col = cols[i-1], xpd = NA)
  if(i == 2 || i == length(angles)) {
    if(i ==2) label = range[1]
    else label = range[2]
    tmp_a = (angles[i-1]+angles[i])/2
    text(srt = tmp_a-90, 
         x = (lineSeq+0.99)*cos(deg2rad(tmp_a)), 
         y =  (lineSeq+0.99)*sin(deg2rad(tmp_a)), 
         xpd = NA, labels = label)
  }
 }
}

re_scale = function(x) return(cov2cor(x))

```


## Case Study 1 - Runtime Wilkinson datasets
```{r}
library(flextable)
result = readRDS("results/case_study_runtime.RDS")
gpu = sapply(result, function(rr) {sapply(rr, function(r) r$gpu$time)})
cpu = sapply(result, function(rr) {sapply(rr, function(r) r$cpu$time)})

result_mean = cbind(apply(cpu,2, mean), apply(gpu,2, mean))
# hours
result_mean/3600
# minutes
result_mean/60

result_bc = readRDS("results/case_study_runtime_BC.RDS")
bc_case = sapply(result_bc, function(rr) {sapply(rr, function(r) r[[1]])})
bc_case_mean = apply(bc_case,2, mean)
bird_speed_up = (bc_case_mean/3600)[1] / 3.8

results_min = data.frame(cbind(rownames(result_mean),round(bc_case_mean/60, 3), round(result_mean/60, 3)))


save_as_docx("minutes" = flextable(results_min), path = "results/runtimes_wilkinson_minutes.docx")
```


## Runtime
### Figure 1 - Run time
Load results:
```{r}
setwd(here::here())
cpu = readRDS(file = "./results/1_cpu_sjSDM_full.RDS")
gpu = readRDS(file = "./results/1_gpu_sjSDM_full.RDS")
gllvm = readRDS(file = "./results/1_gllvm_full.RDS")
bc = readRDS(file = "./results/1_BayesComm_full.RDS")
hmsc = readRDS(file = "./results/1_hmsc_full.RDS")
large_scale = readRDS("./results/large_scale_logit.RDS")
setup2 = large_scale$setup
load("./data_sets_full.RData")
```


```{r}
cols = RColorBrewer::brewer.pal(5, "Set2")
svg(file = paste0("./figures/Fig_1",".svg"), width = 9, height = 2.7)

par(mfrow = c(1,3), mgp = c(2.1,0.6,0), mar = c(3.3, 2.3, 2,0.7), oma = c(0.0,1.3,0.0,0))
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2
site_plot = c("0.1" = "A", "0.3" = "B", "0.5" = "C")
## A
for(i in (as.character(unique(number)))){
  lineT = rep(1,5) #1:5
  names(lineT) =  (as.character(unique(number)))
  if( i == "0.1") {
    ylab = "runtime in minutes"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(-1,3600), xaxt = "n", main = "",yaxt = "n",xlab = "", ylab = ylab, xpd = NA, font = 1)
  text(0.05*(nrow(setup)), 3600*1.09, pos = 2, labels = site_plot[i], font = 2, xpd = NA, cex = 1.5)
  tt = seq(2, (3600), length.out = 5)
  tt = (tt)
  axis(1, at = seq(1, nrow(setup), by = 3)+1.5, labels = rep("", 9))
  axis(2, at = round((tt),1),labels = c(round(tt/60,0)), las = 2)
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -320,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -670,pos = 1, xpd = NA, labels =unique(setup$sites))
  if(i == "0.1"){
  text(x = 2, y = -370-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
  text(x = 2, y = -720-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
  }
    cat(i, "\n")
    X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
    X = seq(1, nrow(setup), by = 3)+1.5
    ind = which(as.character(number) == i & setup$env == e, TRUE)
    mean_conf(X, (cpu$result_time[as.character(number) == i & setup$env == e, ]), col = cols[1], pch = 2)
    mean_conf(X, (gpu$result_time[as.character(number) == i & setup$env == e, ]), col = cols[2], pch = 17)
    mean_conf(X, (gllvm$result_time[as.character(number) == i & setup$env == e, ]), col = cols[3], pch = 19)
    mean_conf(X, (bc$result_time[as.character(number) == i & setup$env == e, ]), col = cols[4], pch = 15)
    mean_conf(X, (hmsc$result_time[as.character(number) == i & setup$env == e, ]/2), col = cols[5], pch = 4)
    if(i == "0.1") legend("topleft", legend = c("CPU-sjSDM - 100 MC samples", "GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))
}
dev.off()

## B
svg(file = paste0("./figures/Fig_2",".svg"), width = 6, height = 4.5)
par(mfrow = c(1,1), mgp = c(2.1,0.6,0), mar = c(4, 2, 1,0.7), oma = c(0.0,1.3,0.0,0))
plot(NULL, NULL, xlim = c(0,1), ylim = c(0.0,5), xaxt = "n",xlab = "",main = "",yaxt = "n",  ylab = "runtime in minutes", xpd = NA,yaxs="i")
coords = matrix(NA, 9, 2)

coords[1:3,1] = c(0, 0.07+0.02, 2*0.07+ 2*0.02)
coords[1:3,2] = c(0, 0.07+0.02, 2*0.07+ 2*0.02)+0.07

coords[4:6, 1] = coords[1:3,1]+0.125+0.25
coords[4:6, 2] = coords[1:3,2]+0.125+0.25

coords[7:9, 1] = coords[1:3,1]+2*0.125+2*0.25
coords[7:9, 2] = coords[1:3,2]+2*0.125+2*0.25

counter = 1
site_labels = c("5000" = "5,000", "15000" = "15,000", "30000" = "30,000")
species_labels = c("300" = "300", "500" = "500", "1000" = "1,000")
xx = unique(setup2$sites)
tt = 0:5
for(j in xx){
  for(i in unique(setup2$species)) {
      
      mu = mean(large_scale$result_time[setup2$species == i & setup2$sites == j]/60)
      rect(xleft = coords[counter, 1], xright = coords[counter,2], ybottom = 0, ytop = mu,
           col = cols[2]) 
      
      se = sd(large_scale$result_time[setup2$species == i & setup2$sites == j]/60)/sqrt(10)
      segments(x0 = mean(coords[counter,]), x1 = mean(coords[counter,]),
              y0 = mu-se, y1=mu+se)
      segments(x0 = mean(coords[counter,])+0.02, x1=mean(coords[counter,])-0.02,
               y0 = mu+se, y1=mu+se)
      segments(x0 = mean(coords[counter,])+0.02, x1=mean(coords[counter,])-0.02,
               y0 = mu-se, y1=mu-se)
          
   
      if(counter==1)axis(2, at = tt,labels = c(round(tt,1)), las = 2, lwd = 0, lwd.ticks = 1, font = 1)
      if(i == "300") {
        if(j == 15000) d = 0.5
        else if(j == 30000) d = 0.2
        else d = 0.5
        text(x = coords[counter,2] + 0.04, y = max(large_scale$result_time[setup2$species == i & setup2$sites == j]/60) + d ,pos = 3, labels = paste0(site_labels[as.character(j)], " sites"))
      }
      
      text(x= coords[counter,1], y = -0.1, pos = 1, srt = 45, labels = species_labels[as.character(i)], xpd = NA)
      counter = counter + 1
  }
}
text(0.5, y = -0.70, pos = 1, label = "Species", xpd = NA)
    legend("topleft",horiz = FALSE, legend = c("GPU-sjSDM"), col = cols[2], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.2)
dev.off()

```



### Figure S1 - Log Runtime
```{r}
svg(file = paste0("./figures/Fig_S1",".svg"), width = 9, height = 2.7)
par(mfrow = c(1,3), mgp = c(3.1,0.6,0), mar = c(3.3, 2.8, 2,0.7), oma = c(0.0,1.3,0.0,0))
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2
site_plot = c("0.1" = "A", "0.3" = "B", "0.5" = "C")
## A
for(i in (as.character(unique(number)))){
  lineT = rep(1,5) #1:5
  names(lineT) =  (as.character(unique(number)))
  if( i == "0.1") {
    ylab = "log runtime in minutes"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(-1.5,11), xaxt = "n", main = "",yaxt = "n",xlab = "", ylab = ylab, xpd = NA, font = 1)
  text(0.05*(nrow(setup)), 11*1.09, pos = 2, labels = site_plot[i], font = 2, xpd = NA, cex = 1.5)
  tt = seq((-2), ((9)), length.out = 5)
  tt = (tt)
  axis(1, at = seq(1, nrow(setup), by = 3)+1.5, labels = rep("", 9))
  axis(2, at = round((tt),3),labels = round(exp(tt)/60,2), las = 2)
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -2.2,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -3.3,pos = 1, xpd = NA, labels =unique(setup$sites))
  if(i == "0.1"){
  text(x = 2, y = -2.3-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
  text(x = 2, y = -3.4-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
  }
    cat(i, "\n")
    X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
    X = seq(1, nrow(setup), by = 3)+1.5
    ind = which(as.character(number) == i & setup$env == e, TRUE)
    mean_conf(X, log(cpu$result_time[as.character(number) == i & setup$env == e, ]), col = cols[1], pch = 2)
    mean_conf(X, log(gpu$result_time[as.character(number) == i & setup$env == e, ]), col = cols[2], pch = 17)
    mean_conf(X, log(gllvm$result_time[as.character(number) == i & setup$env == e, ]), col = cols[3], pch = 19)
    mean_conf(X, log(bc$result_time[as.character(number) == i & setup$env == e, ]), col = cols[4], pch = 15)
    mean_conf(X, log(hmsc$result_time[as.character(number) == i & setup$env == e, ]/2), col = cols[5], pch = 4)

    if(i == "0.1") legend("topleft", legend = c("CPU-sjSDM - 100 MC samples", "GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))
}
dev.off()

```


## Inference
### Figure 3 - Covariance accuracy
Load results:
```{r}
setwd(here::here())
if(!exists("gllvmSparse")) gllvmSparse = readRDS(file = "./results/6_gllvm_sparse.RDS")
if(!exists("gpuSparse")) gpuSparse = readRDS(file = "./results/6_gpu_sjSDM_sparse2.RDS")
if(!exists("bcSparse")) bcSparse = readRDS(file = "./results/6_BayesComm_sparse.RDS")
if(!exists("hmscSparse")) hmscSparse = readRDS(file = "./results/6_hmsc_sparse.RDS")
load("./data_sets_sparse_95.RData")
```

Post-hoc calculation of TSS:
```{r}
transform_pred = function(cov, t=0.001) {
  cov = cov[lower.tri(cov, diag = FALSE)]
  cov[abs(cov) < t] = 0.0
  cov[abs(cov) > t] = 1.0
  return(cov)
}

tss = function(true, pred, t = 0.01) {
  pp = as.factor(transform_pred(pred,t))
  tt = as.factor(transform_pred(true,t))
  levels(pp) = levels(tt) = c(0, 1)
  cf = caret::confusionMatrix(pp, tt)
  return(sum(cf$byClass[1:2])-1)
}

results_sparse = 
  lapply(list(gllvmSparse, gpuSparse, bcSparse, hmscSparse), function(rr) {
    res = array(NA, c(4, 27, 5))
    counter = 1
    for(i in 1:27) {
      for(j in 1:5){
        for(tre in c(1:4)) {
          true = data_sets[[counter]]$sim$correlation
          pred = rr$post[[i]][[j]]$correlation
          try({
  
            if(is.list(pred)) {
              supportLevel = 0.95
              pred = ((pred[[1]]$support>supportLevel) + (pred[[1]]$support<(1-supportLevel))>0)*pred[[1]]$mean
            }
            
            if(any(diag(pred) < 0.001)) {
              pred = pred + t(pred)
              diag(pred) = 1.0
            }
            pred = cov2cor(pred)
            
            res[tre, i,j] =tss(true, pred, t = c(0.005, 0.01, 0.05, 0.1)[tre])
          },silent = TRUE)
        }
        counter=counter+1
      }
    }
    return(res)
  })

gllvmSparse$result_corr_acc = results_sparse[[1]]
gpuSparse$result_corr_acc = results_sparse[[2]]
bcSparse$result_corr_acc = results_sparse[[3]]
hmscSparse$result_corr_acc = results_sparse[[4]]
```


```{r}
cols = RColorBrewer::brewer.pal(5, "Set2")
svg(file = paste0("./figures/Fig_3",".svg"), width = 9, height = 4.9)
par(mfrow = c(2,3), mgp = c(2.1,0.6,0), mar = c(0.6, 2.2, 1.4,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2

## A
counter = 1
for(i in (as.character(unique(number)))){
e = 5L
  lineT = rep(1,5) #1:5
  names(lineT) =  (as.character(unique(number)))
  if( i == "0.1") {
    ylab = "Accuracy"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,27-0.5), ylim = c(0.5,1), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA)
  axis(2, las = 2)
  text(0.05*(nrow(setup)-4), 1*1.05, pos = 2, labels = LETTERS[1:3][counter], font = 2, xpd = NA, cex = 1.5)
    cat(i, "\n")
     X = seq(1, by = 3L,to = nrow(setup))+1.5
    ind = which(as.character(number) == i & setup$env == e, TRUE)
    mean_conf(X, cpu$result_corr_acc[as.character(number) == i & setup$env == e, ], col = cols[1], pch = 2)
    mean_conf(X, gpu$result_corr_acc[as.character(number) == i & setup$env == e, ], col = cols[2], pch = 17)
    mean_conf(X, gllvm$result_corr_acc[as.character(number) == i & setup$env == e, ], col = cols[3], pch = 19)
    mean_conf(X, bc$result_corr_acc[as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
    mean_conf(X, hmsc$result_corr_acc[as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
    if(i == "0.1") legend("topleft", legend = c("CPU-sjSDM - 100 MC samples", "GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))
    counter = counter + 1
  }

# B Sparse
e = 5L
lineT = rep(1,5) #1:5
names(lineT) =  (as.character(unique(number)))
counter = 1
for(i in (as.character(unique(number)))){
  if( i == "0.1") {
    ylab = "TSS"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.0,1), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA, xaxs = "r")
  text(0.05*(nrow(setup)-4), 1*1.10, pos = 2, labels = LETTERS[4:6][counter], font = 2, xpd = NA, cex = 1.5)
  axis(1, at = seq(1, 27, by = 3)+1.5, labels = rep("", 9))
  axis(2, las = 2)
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.066,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.166,pos = 1, xpd = NA, labels =unique(setup$sites))
  if(i == "0.1"){
  text(x = 2, y = -0.079-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
  text(x = 2, y = -0.179-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
  }
  cat(i, "\n")
  X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
  X = seq(1, nrow(setup), by = 3)+1.5
  ind = which(as.character(number) == i & setup$env == e, TRUE)
  mean_conf(X, gpuSparse$result_corr_acc[2,as.character(number) == i & setup$env == e, ], col = cols[2], pch = 17)
  mean_conf(X, gllvmSparse$result_corr_acc[2,as.character(number) == i & setup$env == e, ], col = cols[3], pc = 19)
  mean_conf(X, bcSparse$result_corr_acc[2,as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
  mean_conf(X, hmscSparse$result_corr_acc[2,as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
  
  if(i == "0.1") legend("topleft", legend = c("GPU-sjSDM - 2,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols[-1], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(17, 19, 15, 4))
  counter = counter + 1
}
dev.off()

```

### Figure S3 - Sparse matrices, different TSS

```{r}
svg(file = paste0("./figures/Fig_S3",".svg"), width = 9.0, height = 8.9)
par(mfrow = c(4,3), mgp = c(2.1,0.6,0), mar = c(0.6, 1.3, 0.3,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2
# B Sparse
e = 5L
lineT = rep(1,5) #1:5
for(k in 1:4) {
  names(lineT) =  (as.character(unique(number)))
  for(i in (as.character(unique(number)))){
    if( i == "0.1") {
      ylab = "TSS"
    } else {
      ylab = ""
    }
    plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.0,1), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA, xaxs = "r")
   # title(paste0(as.numeric(i)*100, "% species"), line = 2, xpd = NA)
    if(i == "0.1") text(-0.05*(nrow(setup)-4), 1*1.05, pos = 2, labels = LETTERS[1:4][k], font = 2, xpd = NA, cex = 1.5)
    
    if(k == 4) {
      axis(1, at = seq(1, 27, by = 3)+1.5, labels = rep("", 9))
      axis(2, las = 2)
    
      text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.066,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
      text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.166,pos = 1, xpd = NA, labels =unique(setup$sites))
      if(i == "0.1"){
      text(x = 2, y = -0.079-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
      text(x = 2, y = -0.179-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
      }
    }
    cat(i, "\n")
    X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
    X = seq(1, nrow(setup), by = 3)+1.5
    ind = which(as.character(number) == i & setup$env == e, TRUE)
    mean_conf(X, gpuSparse$result_corr_acc[k,as.character(number) == i & setup$env == e, ], col = cols[2], pch = 17)
    mean_conf(X, gllvmSparse$result_corr_acc[k,as.character(number) == i & setup$env == e, ], col = cols[3], pc = 19)
    mean_conf(X, bcSparse$result_corr_acc[k,as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
    mean_conf(X, hmscSparse$result_corr_acc[k,as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
    
    if(i == "0.1") legend("topleft", legend = c("GPU-sjSDM - 2,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols[-1], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(17, 19, 15, 4))
  
  }
}
dev.off()
```



### Figure S2 - Covariance RMSE
Calculate RMSE:
```{r}
cpu = readRDS(file = "./results/1_cpu_sjSDM_full.RDS")
gpu = readRDS(file = "./results/1_gpu_sjSDM_full.RDS")
gllvm = readRDS(file = "./results/1_gllvm_full.RDS")
bc = readRDS(file = "./results/1_BayesComm_full.RDS")
hmsc = readRDS(file = "./results/1_hmsc_full.RDS")

load("data_sets_full.RData")
counter = 1
result_corr_rmse = array(NA, dim = c(5L, 27L, 5L))
rmse = function(obs, est) return(sqrt(mean((obs[lower.tri(obs)] - est[lower.tri(obs)])^2)))
for(i in 1:27) {
  for(j in 1:5) {
    obs = data_sets[[counter]]$sim$correlation
    result_corr_rmse[1, i, j] = rmse(obs, cov2cor(cpu$post[[i]][[j]]$correlation))
    result_corr_rmse[2, i, j] = rmse(obs, cov2cor(gpu$post[[i]][[j]]$correlation))
    gllvm_rmse = try({rmse(obs, cov2cor(gllvm$post[[i]][[j]]$correlation))}, silent=TRUE)
    if("try-error" %in% class(gllvm_rmse)) gllvm_rmse = NA
    result_corr_rmse[3, i, j] = gllvm_rmse
    bc_rmse = (bc$post[[i]][[j]]$correlation) + t((bc$post[[i]][[j]]$correlation))
    diag(bc_rmse) = 1.0
    result_corr_rmse[4, i, j] = rmse(obs, bc_rmse)
    result_corr_rmse[5, i, j] = rmse(obs, (hmsc$post[[i]][[j]]$correlation))
    counter = counter+1
  }
}

load("data_sets_sparse_95.RData")
gllvmSparse = readRDS(file = "./results/6_gllvm_sparse.RDS")
gpuSparse = readRDS(file = "./results/6_gpu_sjSDM_sparse2.RDS")
bcSparse = readRDS(file = "./results/6_BayesComm_sparse.RDS")
hmscSparse = readRDS(file = "./results/6_hmsc_sparse.RDS")
supportLevel = 0.95
counter = 1
result_corr_rmse_sparse = array(NA, dim = c(4L, 27L, 5L))
rmse = function(obs, est) return(sqrt(mean((obs[lower.tri(obs)] - est[lower.tri(obs)])^2)))
for(i in 1:27) {
  for(j in 1:5) {
    obs = data_sets[[counter]]$sim$correlation
    result_corr_rmse_sparse[1, i, j] = rmse(obs, cov2cor(gpuSparse$post[[i]][[j]]$correlation))
    gllvm_rmse = try({rmse(obs, cov2cor(gllvmSparse$post[[i]][[j]]$correlation))}, silent=TRUE)
    if("try-error" %in% class(gllvm_rmse)) gllvm_rmse = NA
    result_corr_rmse_sparse[2, i, j] = gllvm_rmse
    result_corr_rmse_sparse[3, i, j] = rmse(obs, (bcSparse$post[[i]][[j]]$correlation))
    
    pred = hmscSparse$post[[i]][[j]]$correlation
    pred = ((pred[[1]]$support>supportLevel) + (pred[[1]]$support<(1-supportLevel))>0)*pred[[1]]$mean
    result_corr_rmse_sparse[4, i, j] = rmse(obs, (pred))
    counter = counter+1
  }
}
```

Plot RMSE:
```{r}
cols = RColorBrewer::brewer.pal(5, "Set2")

svg(file = paste0("./figures/Fig_S2",".svg"), width = 9, height = 4.9)
par(mfrow = c(2,3), mgp = c(2.1,0.6,0), mar = c(0.6, 1.6, 0.3,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2
pchs = c(2, 17, 19, 15, 4)

## A
for(i in (as.character(unique(number)))){
e = 5L
  lineT = rep(1,5) #1:5
  names(lineT) =  (as.character(unique(number)))
  if( i == "0.1") {
    ylab = "RMSE"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,27-0.5), ylim = c(0.0,0.8), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA)
  axis(2, las = 2)
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 0.8*1.05, pos = 2, labels = "A", font = 2, xpd = NA, cex = 1.5)
    cat(i, "\n")
     X = seq(1, by = 3L,to = nrow(setup))+1.5
    ind = which(as.character(number) == i & setup$env == e, TRUE)
    for(l in 1:5){
      mean_conf(X, result_corr_rmse[l,,][as.character(number) == i & setup$env == e, ], 
                col = cols[l], pch = pchs[l], spar = 0.0)
    }

  if(i == "0.1") legend("topleft", legend = c("CPU-sjSDM - 100 MC samples", "GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))
  }

pchs2 = pchs[-1]
# B Sparse
e = 5L
lineT = rep(1,5) #1:5
names(lineT) =  (as.character(unique(number)))
for(i in (as.character(unique(number)))){
  if( i == "0.1") {
    ylab = "RMSE"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.0,0.2), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA, xaxs = "r")
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 1*1.05, pos = 2, labels = "B", font = 2, xpd = NA, cex = 1.5)
  axis(1, at = seq(1, 27, by = 3)+1.5, labels = rep("", 9))
  axis(2, las = 2)
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.066*0.2,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.166*0.17,pos = 1, xpd = NA, labels =unique(setup$sites))
  if(i == "0.1"){
  text(x = 2, y = -0.079*0.2-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
  text(x = 2, y = -0.179*0.17-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
  }
  
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 0.2*1.05, pos = 2, labels = "B", font = 2, xpd = NA, cex = 1.5)
  cat(i, "\n")
  X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
  X = seq(1, nrow(setup), by = 3)+1.5
  ind = which(as.character(number) == i & setup$env == e, TRUE)
    for(l in 1:4){
      mean_conf(X, result_corr_rmse_sparse[l,,][as.character(number) == i & setup$env == e, ], 
                col = cols[l+1], pch = pchs2[l], spar = 0.0)
    }
  legend_pos = "topleft"
  if(i == "0.1") legend_pos = "bottomleft"
  if(i=="0.1") legend(legend_pos, legend = c("GPU-sjSDM - 2,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols[-1], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(17, 19, 15, 4))
}

dev.off()

```




### Figure S4 - Env Acc

```{r}
cols = RColorBrewer::brewer.pal(5, "Set2")

svg(file = paste0("./figures/Fig_S4",".svg"), width = 9, height = 9)
par(mfrow = c(4,3), mgp = c(2.1,0.6,0), mar = c(0.6, 1.3, 0.3,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2
## A
for(i in (as.character(unique(number)))){
  lineT = rep(1,5) #1:5
  names(lineT) =  (as.character(unique(number)))
  if( i == "0.1") {
    ylab = "Env accuracy"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.0,1.0), xaxt = "n", main = "",yaxt = "n",xlab = "", ylab = ylab, xpd = NA, font = 1)
  if(i == "0.1") text(-0.05*(nrow(setup)), 1.0*1.09, pos = 2, labels = "A", font = 2, xpd = NA, cex = 1.5)
  tt = seq(0.0, (1.0), length.out = 5)
  axis(2,las = 2)
  cat(i, "\n")
  X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
  X = seq(1, nrow(setup), by = 3)+1.5
  ind = which(as.character(number) == i & setup$env == e, TRUE)
  mean_conf(X, (cpu$result_env[as.character(number) == i & setup$env == e, ]), col = cols[1], pch = 2)
  mean_conf(X, (gpu$result_env[as.character(number) == i & setup$env == e, ]), col = cols[2], pch = 17)
  mean_conf(X, (gllvm$result_env[as.character(number) == i & setup$env == e, ]), col = cols[3], pch = 19)
  mean_conf(X, (bc$result_env[as.character(number) == i & setup$env == e, ]), col = cols[4], pch = 15)
  mean_conf(X, (hmsc$result_env[as.character(number) == i & setup$env == e, ]), col = cols[5], pch = 4)
  if(i == "0.1") legend("bottomleft", legend = c("CPU-sjSDM - 100 MC samples", "GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))
}

## B
e = 5L
lineT = rep(1,5) #1:5
names(lineT) =  (as.character(unique(number)))
for(i in (as.character(unique(number)))){
  if( i == "0.1") {
    ylab = "Env RMSE"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.0,1), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA, xaxs = "r")
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 1*1.10, pos = 2, labels = "B", font = 2, xpd = NA, cex = 1.5)
  axis(2, las = 2)
  cat(i, "\n")
  X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
  X = seq(1, nrow(setup), by = 3)+1.5
  ind = which(as.character(number) == i & setup$env == e, TRUE)
  mean_conf(X, cpu$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[1], pch = 2)
  mean_conf(X, gpu$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[2], pch = 17)
  mean_conf(X, gllvm$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[3], pch = 19)
  mean_conf(X, bc$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
  mean_conf(X, hmsc$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
  if(i == "0.1") legend("topleft", legend = c("CPU-sjSDM - 100 MC samples", "GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))
}



## C
e = 5L
lineT = rep(1,5) #1:5
names(lineT) =  (as.character(unique(number)))
for(i in (as.character(unique(number)))){
  if( i == "0.1") {
    ylab = "Env accuracy"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.5,1.0), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA, xaxs = "r")
  # title(paste0(as.numeric(i)*100, "% species"), line = 2, xpd = NA)
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 1*1.06, pos = 2, labels = "C", font = 2, xpd = NA, cex = 1.5)
  axis(2, las = 2)
  cat(i, "\n")
  X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
  X = seq(1, nrow(setup), by = 3)+1.5
  mean_conf(X, gpuSparse$result_env[as.character(number) == i & setup$env == e, ], col = cols[2], pch = 17)
  mean_conf(X, gllvmSparse$result_env[as.character(number) == i & setup$env == e, ], col = cols[3], pch = 19)
  mean_conf(X, bcSparse$result_env[as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
  mean_conf(X, hmscSparse$result_env[as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
  if(i == "0.1") legend("bottomleft", legend = c( "GPU-sjSDM - 2,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols[-1], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c( 17, 19, 15, 4))
}

# D

e = 5L
lineT = rep(1,5) #1:5
names(lineT) =  (as.character(unique(number)))
for(i in (as.character(unique(number)))){
  if( i == "0.1") {
    ylab = "Env RMSE"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.0,1.3), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA, xaxs = "r")
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 1.3*1.06, pos = 2, labels = "D", font = 2, xpd = NA, cex = 1.5)
  axis(1, at = seq(1, 27, by = 3)+1.5, labels = rep("", 9))
  axis(2, las = 2)
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.1,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.22,pos = 1, xpd = NA, labels =unique(setup$sites))
  if(i == "0.1"){
  text(x = 2, y = -0.11-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
  text(x = 2, y = -0.23-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
  }
  cat(i, "\n")
  X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
  X = seq(1, nrow(setup), by = 3)+1.5
  ind = which(as.character(number) == i & setup$env == e, TRUE)
  mean_conf(X, gpuSparse$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[2], pch = 17)
  mean_conf(X, gllvmSparse$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[3], pch = 19)
  mean_conf(X, bcSparse$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
  mean_conf(X, hmscSparse$result_rmse_env[as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
  if(i == "0.1") legend("topleft", legend = c("GPU-sjSDM - 2,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols[-1], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c( 17, 19, 15, 4))
  
}

dev.off()

```

### Table large scale accuracies
```{r}
cov_acc = apply(large_scale$result_corr_acc, 1, mean)
env_acc = apply(large_scale$result_env, 1, mean)
env_rmse = apply(large_scale$result_rmse_env, 1, mean)
write.csv2(round(rbind(cov_acc, env_acc, env_rmse),3), "results/large_scale_results.csv")
```


### Figure S6 - Covariance behaviour
Load results:
```{r}
setwd(here::here())
gpu_beh = readRDS(file = "./results/gpu_sjSDM_behaviour_sites.RDS")
gllvm_beh = readRDS(file = "./results/gllvm_behaviour_sites.RDS")
bc_beh = readRDS(file = "./results/bc_behaviour_sites.RDS")
hmsc_beh = readRDS(file = "./results/hmsc_behaviour_sites.RDS")

```


```{r}
svg(file = paste0("./figures/Fig_S6",".svg"), width = 9, height = 4.3)
cols = RColorBrewer::brewer.pal(5, "Set2")
xx = 1:25
sites = seq(50,to = 350, length.out = 7)
par(mfrow = c(1,2), mar = c(2.1, 3, 1.6, 1), mgp = c(2.4, 1, 0), oma = c(1.4,1.3,0.3,0))

plot(NULL, NULL, xlim = c(min(sites), max(sites)), ylim = c(0.5, 1.0), ylab = "accuracy", xlab = "Sites", yaxt = "n", xpd = NA, main = "", xaxt = "n")
axis(1, at = sites[seq(1, length(sites), by = 2)], labels = sites[seq(1, length(sites), by = 2)])
text(x = 27, y = 1.05, pos = 2, labels = "A", xpd = NA, cex = 1.5, font = 2)
title("covariance accuracy",line = 1, xpd = NA)
axis(2, las = 2)
ind = 1:nrow(gpu_beh$result_corr_acc)
mean_conf(sites,mat = gpu_beh$result_corr_acc, col = cols[2], pch = 17)
mean_conf(sites,gllvm_beh$result_corr_acc, col = cols[3], pch = 19)
mean_conf(sites,bc_beh$result_corr_acc,col=cols[4], pch = 15)
mean_conf(sites,hmsc_beh$result_corr_acc,col=cols[5],pch = 4)
legend("topleft", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols[-1], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c( 17, 19, 15, 4))
plot(NULL, NULL, xlim = c(min(sites), max(sites)), ylim = c(0.0, 1.0), ylab = "RMSE", xlab = "Sites", yaxt = "n", xpd = NA, main = "", xaxt = "n")
axis(1, at = sites[seq(1, length(sites), by = 2)], labels = sites[seq(1, length(sites), by = 2)])
text(x = 27, y = 1.1, pos = 2, labels = "B", xpd = NA, cex = 1.5, font = 2)
title("env RMSE",line = 1, xpd = NA)
axis(2, las = 2)
mean_conf(sites, gpu_beh$result_rmse_env, col = cols[2], pch = 17)
mean_conf(sites, gllvm_beh$result_rmse_env, col = cols[3], pch = 19)
mean_conf(sites, bc_beh$result_rmse_env, col = cols[4], pch = 15)
mean_conf(sites, hmsc_beh$result_rmse_env, col = cols[5], pch = 4)
dev.off()

```


## Predictive performance
### Figure 4 - Predictive Performance
Load results:
```{r}
setwd(here::here())
load("./data_sets_full.RData")
library(sjSDM)
p_gpu = p_cpu = p_gllvm = p_bc = p_hmsc= matrix(NA, 27, 5)
acc = function(t, p) Metrics::auc(t, p)
counter = 1
for(i in 1:27) {
  for(j in 1:5) {
    p_gpu[i,j] = acc( gpu$auc[[i]][[j]]$true, gpu$auc[[i]][[j]]$pred )
    p_cpu[i,j] = acc( cpu$auc[[i]][[j]]$true, cpu$auc[[i]][[j]]$pred )
    if(!is.null(gllvm$auc[[i]][[j]]$pred )) p_gllvm[i,j] = acc( gllvm$auc[[i]][[j]]$true, gllvm$auc[[i]][[j]]$pred )
    p_bc[i,j] = acc( bc$auc[[i]][[j]]$true, bc$auc[[i]][[j]]$pred )
    p_hmsc[i,j] = acc( hmsc$auc[[i]][[j]]$true, hmsc$auc[[i]][[j]]$pred )
    counter = counter+1
  }
}
```



```{r}
cols = RColorBrewer::brewer.pal(5, "Set2")
svg(file = paste0("./figures/Fig_4",".svg"), width = 9, height = 2.7)
e = 5L
lineT = rep(1,5) #1:5
names(lineT) =  (as.character(unique(number)))
par(mfrow = c(1,3), mgp = c(2.1,0.6,0), mar = c(3.3, 2.3, 2,0.7), oma = c(0.0,1.3,0.0,0))
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2
site_plot = c("0.1" = "A", "0.3" = "B", "0.5" = "C")
ylim = c(0.7,0.8)
for(i in (as.character(unique(number)))){
  if( i == "0.1") {
    ylab = "AUC"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1,nrow(setup)), ylim = ylim, xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, main = "", xaxs = "i", xpd = NA)
  text(0.05*(nrow(setup))-1, 0.81, pos = 2, labels = site_plot[i], font = 2, xpd = NA, cex = 1.5)
  axis(1, at = seq(1, nrow(setup), by = 3)+1.5, labels = rep("", 9))
  axis(2, las = 2)
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = 0.693,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = 0.685,pos = 1, xpd = NA, labels = unique(setup$sites))
  if(i == "0.1"){
  text(x = 2, y = 0.692-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
  text(x = 2, y = 0.684-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
  }
  cat(i, "\n")
  X = seq(1, by = 3L,to = nrow(setup))+1.5
  mean_conf(X, p_cpu[as.character(number) == i & setup$env == e, ], col = cols[1], pch = 2)
  mean_conf(X, p_gpu[as.character(number) == i & setup$env == e, ], col = cols[2], pch = 17)
  mean_conf(X, rbind(p_gllvm, matrix(NA, length(number) - nrow(p_gllvm) ,5))[as.character(number) == i & setup$env == e, ], col = cols[3], pch = 19)
  mean_conf(X, rbind(p_bc, matrix(NA, length(number) - nrow(p_bc) ,5))[as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
  mean_conf(X, rbind(p_hmsc, matrix(NA, length(number) - nrow(p_hmsc) ,5))[as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
 if(i == "0.5")legend("bottomright", legend = c("CPU-sjSDM - 100 MC samples", "GPU-sjSDM - 1,000 MC samples", "gllvm", "BayesComm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))

}
dev.off()
```

## Convergence checks
### Figure S5
```{r}
if(!exists("bc")) bc = readRDS(file = "./results/BayesCommDiag.RDS")
if(!exists("hmsc")) hmsc = readRDS(file = "./results/hmsc.RDS")
if(!exists("bcSparse")) bcSparse = readRDS(file = "./results/sparse_bc.RDS")
if(!exists("hmscSparse")) hmscSparse = readRDS(file = "./results/sparse_hmsc2.RDS")
remove_last_factor = function(lambda) {
  fact_numbers = as.numeric(sapply(strsplit(names(lambda), ","), function(X) gsub("[^0-9.]", "", X[2])))
  return(lambda[fact_numbers != max(fact_numbers)])
}
bc_conv = hmsc_conv = bc_sparse_conv = hmsc_sparse_conv = matrix(NA, 27L, 5L)
for(i in 1:27) {
  for(j in 1:5) {
    bc_conv[i,j] = mean(c(bc$post[[i]][[j]]$psrf.gamma, bc$post[[i]][[j]]$beta.conv))
    bc_sparse_conv[i,j] = mean(c(bcSparse$post[[i]][[j]]$psrf.gamma, bcSparse$post[[i]][[j]]$beta.conv))
    
    hmsc_conv[i,j] = mean(c(hmsc$post[[i]][[j]]$beta.conv, remove_last_factor( hmsc$post[[i]][[j]]$lambda.conv )))
    hmsc_sparse_conv[i,j] = mean(c(hmscSparse$post[[i]][[j]]$beta.conv, remove_last_factor( hmscSparse$post[[i]][[j]]$lambda.conv )))
  }
}
cols = RColorBrewer::brewer.pal(5, "Set2")

svg(file = paste0("./figures/Fig_S5",".svg"), width = 9, height = 4.9)

par(mfrow = c(2,3), mgp = c(2.1,0.6,0), mar = c(0.6, 1.3, 0.3,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
number = setup$species
spar = 0.5
e = 5L
lwd = 2.2

## A
for(i in (as.character(unique(number)))){
e = 5L
  lineT = rep(1,5) #1:5
  names(lineT) =  (as.character(unique(number)))
  if( i == "0.1") {
    ylab = "PSRF > 1.2"
  } else {
    ylab = "" 
  }
  plot(NULL, NULL, xlim = c(1+1.5,27-0.5), ylim = c(0.0,1), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA)
  axis(2, las = 2)
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 1*1.05, pos = 2, labels = "A", font = 2, xpd = NA, cex = 1.5)
    cat(i, "\n")
     X = seq(1, by = 3L,to = nrow(setup))+1.5
    ind = which(as.character(number) == i & setup$env == e, TRUE)
    mean_conf(X, bc_conv[as.character(number) == i & setup$env == e, ], col = cols[4], pch = 15)
    mean_conf(X, hmsc_conv[as.character(number) == i & setup$env == e, ], col = cols[5], pch = 4)
  if(i == "0.1")legend("topleft", legend = c("BayesComm", "Hmsc"), col = cols[4:5], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))
  }


# B Sparse
e = 5L
lineT = rep(1,5) #1:5
names(lineT) =  (as.character(unique(number)))
for(i in (as.character(unique(number)))){
  if( i == "0.1") {
    ylab = "PSRF > 1.2"
  } else {
    ylab = ""
  }
  plot(NULL, NULL, xlim = c(1+1.5,nrow(setup)-0.5), ylim = c(0.0,1), xaxt = "n", yaxt = "n", xlab = "", ylab = ylab, xpd = NA, xaxs = "r")
  if(i == "0.1") text(-0.05*(nrow(setup)-4), 1*1.05, pos = 2, labels = "B", font = 2, xpd = NA, cex = 1.5)
  axis(1, at = seq(1, 27, by = 3)+1.5, labels = rep("", 9))
  axis(2, las = 2)
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.066,pos = 1, xpd = NA, labels = unique(setup$sites)* as.numeric(i))
  text(x = seq(1, nrow(setup), by = 3)+1.5, y = -0.166,pos = 1, xpd = NA, labels =unique(setup$sites))
  if(i == "0.1"){
  text(x = 2, y = -0.079-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Species:", pos = 2, xpd = NA)
  text(x = 2, y = -0.179-strheight(unique(setup$sites)* as.numeric(i))[1], labels = "Sites:", pos = 2, xpd = NA)
  }
  cat(i, "\n")
  X = (1:nrow(setup))[as.character(number) == i & setup$env == e]+1.5
  X = seq(1, nrow(setup), by = 3)+1.5
  ind = which(as.character(number) == i & setup$env == e, TRUE)
  mean_conf(X, bc_sparse_conv[as.character(number) == i & setup$env == e, ], col = cols[4], pch = 19)
  mean_conf(X, hmsc_sparse_conv[as.character(number) == i & setup$env == e, ], col = cols[5], pch = 15)
  
  if(i == "0.1")legend("topleft", legend = c( "BayesComm", "Hmsc"), col = cols[-1][3:4], bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 17, 19, 15, 4))

}
dev.off()
```


## Simulating from LVM-JSDM
### Figure S7 - LVM Env RMSE
```{r}
load("results/LVMsimulation_scenarios.RData")

# 1data_10_E 
# 2data_10_EL 
# 3data_10_EL_51
# 4data_10_EL_15
# 
# 5data_50_E   
# 6data_50_EL 
# 7data_50_EL_51  
# 8data_50_EL_15  
# 
# 9data_100_E    
# 10data_100_EL 
# 11data_100_EL_51  
# 12data_100_EL_15  
svg(file = paste0("./figures/Fig_S7",".svg"), width = 9, height = 7)
cols = RColorBrewer::brewer.pal(5, "Set2")
# 1 3 5
cols = cols[c(2, 3, 5)]
ylim = c(0, 0.8)
species_label = c("")
par(mfrow = c(3,3), mgp = c(2.1,0.6,0), mar = c(0.6, 1.3, 0.3,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
#### EL ####
counter = 1
ylab = "Env RMSE"
for( i in c(1, 5, 9)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab,xlab="", xpd = NA)
  ylab = ""
  text(3, y = 1.04*ylim[2], pos=3, label = species_label[counter], xpd = NA, font =1)
  if(i == 1) text(-1, y = 1.04*ylim[2], pos=3, label = "A", xpd = NA, font =2, cex = 1.2)
  counter = counter+1
  mean_conf(0:5, t(cbind(sjSDM_res[[i]]$rmse[,1], sjSDM_res[[i+1]]$rmse)) , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(0:5, t(cbind(gllvm_res[[i]]$rmse[,1], gllvm_res[[i+1]]$rmse)) , spar = 0.0, col = cols[2], pch = 19)
  mean_conf(0:5,  t(cbind(Hmsc_res[[i]]$rmse[,1],  Hmsc_res[[i+1]]$rmse)) , spar = 0.4, col = cols[3], pch = 4)
  if(i == 1)legend("topleft", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
}

#### 51 ####
ylab = "Env RMSE"
for( i in c(3, 7, 11)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
  ylab = ""
  if(i == 3) text(-1, y = 1.04*ylim[2], pos=3, label =  "B", xpd = NA, font =2, cex = 1.2)
  if(i == 3) legend("topleft", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  mean_conf(1:5, sjSDM_res[[i]]$rmse , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(1:5, gllvm_res[[i]]$rmse , spar = 0.0, col = cols[2], pch = 19)
  mean_conf(1:5,  Hmsc_res[[i]]$rmse , spar = 0.4, col = cols[3], pch = 4)
}


### 15 #####
ylab = "Env RMSE"
for( i in c(4, 8, 12)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
  if(i == 4) text(-1, y = 1.04*ylim[2], pos=3, label =  "C", xpd = NA, font =2, cex = 1.2)
  text(2.5, y = -0.15*ylim[2], pos=1, label =  "Number of latent variables", xpd = NA, font =1, cex = 1.2)
  ylab = ""
    if(i == 4)legend("topleft", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  axis(1, at = 0:5, labels = 0:5)
  mean_conf(1:5, sjSDM_res[[i]]$rmse , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(1:5, gllvm_res[[4]]$rmse , spar = 0.0, col = cols[2], pch = 19)
  mean_conf(1:5,  Hmsc_res[[i]]$rmse , spar = 0.4, col = cols[3], pch = 4)
}

dev.off()
```


### Figure S8 - LVM Cov accuracy
```{r}
svg(file = paste0("./figures/Fig_S8",".svg"), width = 9, height = 7)
### Cov acc ###
ylim = c(0, 1)
par(mfrow = c(3,3), mgp = c(2.1,0.6,0), mar = c(0.6, 1.3, 0.3,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
#### EL ####
counter = 1
ylab = "Cov accuracy"
for( i in c(1, 5, 9)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
  text(3, y = 1.04*ylim[2], pos=3, label = species_label[counter], xpd = NA, font =1)
  ylab = ""
  if(i == 1) text(-1, y = 1.04*ylim[2], pos=3, label = "A", xpd = NA, font =2, cex = 1.2)
  counter = counter+1
  if(i == 1) legend("bottomright", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  mean_conf(0:5, t(cbind(sjSDM_res[[i]]$cov[,1], sjSDM_res[[i+1]]$cov)) , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(0:5, t(cbind(gllvm_res[[i]]$cov[,1], gllvm_res[[i+1]]$cov)) , spar = 0.4, col = cols[2], pch = 19)
  mean_conf(0:5,  t(cbind(Hmsc_res[[i]]$cov[,1],  Hmsc_res[[i+1]]$cov)) , spar = 0.4, col = cols[3], pch = 4)
}

#### 51 ####
ylab = "Cov accuracy"
for( i in c(3, 7, 11)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
  if(i == 3) text(-1, y = 1.04*ylim[2], pos=3, label =  "B", xpd = NA, font =2, cex = 1.2)
  ylab = ""
  if(i == 3) legend("bottomright", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  mean_conf(1:5, sjSDM_res[[i]]$cov , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(1:5, gllvm_res[[i]]$cov , spar = 0.4, col = cols[2], pch = 19)
  mean_conf(1:5,  Hmsc_res[[i]]$cov , spar = 0.4, col = cols[3], pch = 4)
}


### 15 #####
ylab = "Cov accuracy"
for( i in c(4, 8, 12)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
  if(i == 4) text(-1, y = 1.04*ylim[2], pos=3, label =  "C", xpd = NA, font =2, cex = 1.2)
  ylab = ""
  text(2.5, y = -0.15*ylim[2], pos=1, label =  "Number of latent variables", xpd = NA, font =1, cex = 1.2)
  if(i == 4) legend("bottomright", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  axis(1, at = 0:5, labels = 0:5)
  mean_conf(1:5, sjSDM_res[[i]]$cov , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(1:5, gllvm_res[[4]]$cov , spar = 0.4, col = cols[2], pch = 19)
  mean_conf(1:5,  Hmsc_res[[i]]$cov , spar = 0.4, col = cols[3], pch = 4)
}
dev.off()
```



### Figure S9 - LVM Cov RMSE
```{r}

svg(file = paste0("./figures/Fig_S9",".svg"), width = 9, height = 7)
### Cov rmse ###
ylim = c(0, 2)
par(mfrow = c(3,3), mgp = c(2.1,0.6,0), mar = c(0.6, 1.3, 0.3,0.7), oma = c(2.6,3.0,2.2,0)+0.3)
#### EL ####
counter = 1
ylab = "Cov RMSE"
for( i in c(1, 5, 9)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
    text(3, y = 1.04*ylim[2], pos=3, label = species_label[counter], xpd = NA, font =1)
    ylab = ""
  if(i == 1) text(-1, y = 1.04*ylim[2], pos=3, label = "A", xpd = NA, font =2, cex = 1.2)
  counter = counter+1
  if(i == 1) legend("topleft", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  mean_conf(0:5, t(cbind(sjSDM_res[[i]]$cov_rmse[,1], sjSDM_res[[i+1]]$cov_rmse)) , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(0:5, t(cbind(gllvm_res[[i]]$cov_rmse[,1], gllvm_res[[i+1]]$cov_rmse)) , spar = 0.4, col = cols[2], pch = 19)
  mean_conf(0:5,  t(cbind(Hmsc_res[[i]]$cov_rmse[,1],  Hmsc_res[[i+1]]$cov_rmse)) , spar = 0.4, col = cols[3], pch = 4)
}

#### 51 ####
ylab = "Cov RMSE"
for( i in c(3, 7, 11)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
      if(i == 3) text(-1, y = 1.04*ylim[2], pos=3, label =  "B", xpd = NA, font =2, cex = 1.2)
  ylab = ""
  if(i == 3) legend("topleft", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  mean_conf(1:5, sjSDM_res[[i]]$cov_rmse , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(1:5, gllvm_res[[i]]$cov_rmse , spar = 0.4, col = cols[2], pch = 19)
  mean_conf(1:5,  Hmsc_res[[i]]$cov_rmse , spar = 0.4, col = cols[3], pch = 4)
}


### 15 #####
ylab = "Cov RMSE"
ylim = c(0, 30)
for( i in c(4, 8, 12)){
  plot(NULL, NULL, xlim = c(0, 5), ylim = ylim, las = 1, xaxt = "n", ylab = ylab, xlab="", xpd = NA)
    if(i == 4) text(-1, y = 1.04*ylim[2], pos=3, label =  "C", xpd = NA, font =2, cex = 1.2)
  ylab = ""
  text(2.5, y = -0.15*ylim[2], pos=1, label =  "Number of latent variables", xpd = NA, font =1, cex = 1.2)
  if(i == 4) legend("topleft", legend = c("GPU-sjSDM - 1,000 MC samples", "gllvm", "Hmsc"), col = cols, bty="n", lty = 1,  seg.len = 0.5, lwd = 2, xpd = NA, cex = 1.0, x.intersp = 1.0, pch = c(2, 19, 4))
  axis(1, at = 0:5, labels = 0:5)
  mean_conf(1:5, sjSDM_res[[i]]$cov_rmse , spar = 0.4, col = cols[1], pch = 2)
  mean_conf(1:5, gllvm_res[[4]]$cov_rmse , spar = 0.4, col = cols[2], pch = 19)
  mean_conf(1:5,  Hmsc_res[[i]]$cov_rmse , spar = 0.4, col = cols[3], pch = 4)
}
dev.off()
```


## Case study - environmental DNA
### Figure 5
```{r}
load("data/eDNA/SDM_data.RData")
library(sjSDM)

occ = ttab_ds[[7]]
dim(occ)
rates = apply(occ, 2, mean)
occ = occ[,rates > 2/125]
dim(occ)

env_scaled = scale(env2)
results = readRDS("results/5_fungi_models.RDS")
base_sigma = re_scale(getCov(results$base))
tuned_sigma = re_scale(getCov(results$best2))

svg(file = paste0("./figures/Fig_5",".svg"), width = 7.0, height = 4.6) 
layout(matrix(c(1,2,rep(3,4)), 2,3 ,byrow = F), c(1.0,1.0,1.0), c(1,1))
number = 20
#par(mfrow = c(2,2))
par( mar = c(1,2,2.1,2)+0.3)
#A
sigma = base_sigma[order(apply(occ, 2, sum)), order(apply(occ, 2, sum))]
sigmas = sigma[upper.tri(sigma)]
upper = order(sigmas, decreasing = TRUE)[1:number]
lower = order(sigmas, decreasing = FALSE)[1:number]
cuts = cut(sigmas, breaks = seq(-1,1,length.out = 21))
to_plot = (1:length(sigmas) %in% upper) | (1:length(sigmas) %in% lower)
levels(cuts) = viridis::viridis(20)
cuts = as.character(cuts)
n = ncol(sigma)
lineSeq = 4.7
nseg = 100
plot(NULL, NULL, xlim = c(-5,5), ylim =c(-5,5),pty="s", axes = F, xlab = "", ylab = "")
text(x = -6, y = 5.7, pos = 3, xpd = NA, labels = "A", font = 2, cex = 1.5)

xx = lineSeq*cos( seq(0,2*pi, length.out=nseg) )
yy = lineSeq*sin( seq(0,2*pi, length.out=nseg) )
polygon(xx,yy, col= "white", border = "black", lty = 1, lwd = 1)
angles = seq(0,355,length.out = n+1)[1:(n)]
xx = cos(deg2rad(angles))*lineSeq
yy = sin(deg2rad(angles))*lineSeq

counter = 1
coords = cbind(xx, yy, angles)
for(i in 1:n) {
  for(j in i:n){
    if(i!=j) {
      if(to_plot[counter]) add_curve(coords[i,], coords[j,], col = cuts[counter], n = 5, lineSeq = lineSeq)
      counter = counter + 1
    }
  }
}

lineSeq = 5.0
occ_logs = log(sort(apply(occ, 2, sum)))
cuts = cut(occ_logs, breaks = 10)
cols = viridis::magma(10) #colfunc(5)
levels(cuts) = cols
for(i in 1:length(occ_logs)){
  p1 = coords[i,]
 x1 = c(cos(deg2rad(p1[3]))*(lineSeq+0.1), cos(deg2rad(p1[3]))*(lineSeq+0.3))
 y1 = c(sin(deg2rad(p1[3]))* (lineSeq+0.1), sin(deg2rad(p1[3]))* (lineSeq+0.3))
 segments(x0 = x1[1], x1 = x1[2], y0 = y1[1], y1 = y1[2], col = as.character(cuts[i]), lend = 1)
}
add_legend(viridis::viridis(20), angles = c(140,110))
text(cos(deg2rad(123))*(lineSeq+0.7), sin(deg2rad(123))*(lineSeq+0.7)+0.5, labels = "correlation", pos = 2, xpd = NA)
add_legend(cols = cols, range = c(3, 112), angles = c(70,40))
text(cos(deg2rad(53))*(lineSeq+0.7), sin(deg2rad(55))*(lineSeq+0.7)+0.5, labels = "Sp. abundance", pos = 4, xpd = NA)
### arrows
segments(x0 = cos(deg2rad(-1))*(lineSeq-0.2), x1 = cos(deg2rad(-1))*(lineSeq+0.9), 
         y0 = sin(deg2rad(-1))*(lineSeq-0.2), y1 = sin(deg2rad(-1))*(lineSeq+0.9), xpd = NA)
segments(x0 = cos(deg2rad(356))*(lineSeq-0.2), x1 = cos(deg2rad(356))*(lineSeq+0.9), 
         y0 = sin(deg2rad(356))*(lineSeq-0.2), y1 = sin(deg2rad(356))*(lineSeq+0.9), xpd = NA)
# first
angles = seq(150,195,length.out = n+1)[1:(n)]
xx = cos(deg2rad(angles))*(lineSeq+0.6)
yy = sin(deg2rad(angles))*(lineSeq+0.6)
lines(xx, yy, xpd = NA)
end = curve_text(195+3, "Species",lineSeq = lineSeq+0.6,reverse = TRUE)

# second
angles = seq(rad2deg(end)+3,rad2deg(end)+45+8,length.out = n+1)[1:(n)]
xx = cos(deg2rad(angles))*(lineSeq+0.6)
yy = sin(deg2rad(angles))*(lineSeq+0.6)
lines(xx, yy, xpd = NA)
arrow_angle = max(angles)-2.8
polygon(x = c(cos(deg2rad(arrow_angle))*(lineSeq+0.5), cos(deg2rad(arrow_angle))*(lineSeq+0.7), cos(deg2rad(max(angles)))*(lineSeq+0.6), cos(deg2rad(arrow_angle))*(lineSeq+0.5)),
        y = c(sin(deg2rad(arrow_angle))*(lineSeq+0.5), sin(deg2rad(arrow_angle))*(lineSeq+0.7), sin(deg2rad(max(angles)))*(lineSeq+0.6), sin(deg2rad(arrow_angle))*(lineSeq+0.5)),col = "black", xpd = NA)


##############

sigma = tuned_sigma[order(apply(occ, 2, sum)), order(apply(occ, 2, sum))]
sigmas = sigma[upper.tri(sigma)]
upper = order(sigmas, decreasing = TRUE)[1:number]
lower = order(sigmas, decreasing = FALSE)[1:number]
cuts = cut(sigmas, breaks = seq(-1,1,length.out = 21))
to_plot = (1:length(sigmas) %in% upper) | (1:length(sigmas) %in% lower)
levels(cuts) = viridis::viridis(20)
cuts = as.character(cuts)
n = ncol(sigma)
lineSeq = 4.7
nseg = 100
plot(NULL, NULL, xlim = c(-5,5), ylim =c(-5,5),pty="s", axes = F, xlab = "", ylab = "")
text(x = -6, y = 5.7, pos = 3, xpd = NA, labels = "B", font = 2, cex = 1.5)

xx = lineSeq*cos( seq(0,2*pi, length.out=nseg) )
yy = lineSeq*sin( seq(0,2*pi, length.out=nseg) )
polygon(xx,yy, col= "white", border = "black", lty = 1, lwd = 1)
angles = seq(0,355,length.out = n+1)[1:(n)]
xx = cos(deg2rad(angles))*lineSeq
yy = sin(deg2rad(angles))*lineSeq

counter = 1
coords = cbind(xx, yy, angles)
for(i in 1:n) {
  for(j in i:n){
    if(i!=j) {
      if(to_plot[counter]) add_curve(coords[i,], coords[j,], col = cuts[counter], n = 5, lineSeq = lineSeq)
      counter = counter + 1
      #cat(counter, "\n")
    }
  }
}

lineSeq = 5.0
occ_logs = log(sort(apply(occ, 2, sum)))
cuts = cut(occ_logs, breaks = 10)
cols = viridis::magma(10) #colfunc(5)
levels(cuts) = cols
for(i in 1:length(occ_logs)){
  p1 = coords[i,]
 x1 = c(cos(deg2rad(p1[3]))*(lineSeq+0.1), cos(deg2rad(p1[3]))*(lineSeq+0.3))
 y1 = c(sin(deg2rad(p1[3]))* (lineSeq+0.1), sin(deg2rad(p1[3]))* (lineSeq+0.3))
 segments(x0 = x1[1], x1 = x1[2], y0 = y1[1], y1 = y1[2], col = as.character(cuts[i]), lend = 1)
}
add_legend(viridis::viridis(20), angles = c(140,110))
text(cos(deg2rad(123))*(lineSeq+0.7), sin(deg2rad(123))*(lineSeq+0.7)+0.5, labels = "correlation", pos = 2, xpd = NA)
add_legend(cols = cols, range = c(3, 112), angles = c(70,40))
text(cos(deg2rad(53))*(lineSeq+0.7), sin(deg2rad(55))*(lineSeq+0.7)+0.5, labels = "Sp. abundance", pos = 4, xpd = NA)
### arrows
segments(x0 = cos(deg2rad(-1))*(lineSeq-0.2), x1 = cos(deg2rad(-1))*(lineSeq+0.9), 
         y0 = sin(deg2rad(-1))*(lineSeq-0.2), y1 = sin(deg2rad(-1))*(lineSeq+0.9), xpd = NA)
segments(x0 = cos(deg2rad(356))*(lineSeq-0.2), x1 = cos(deg2rad(356))*(lineSeq+0.9), 
         y0 = sin(deg2rad(356))*(lineSeq-0.2), y1 = sin(deg2rad(356))*(lineSeq+0.9), xpd = NA)
# first
angles = seq(150,195,length.out = n+1)[1:(n)]
xx = cos(deg2rad(angles))*(lineSeq+0.6)
yy = sin(deg2rad(angles))*(lineSeq+0.6)
lines(xx, yy, xpd = NA)
end = curve_text(195+3, "Species",lineSeq = lineSeq+0.6,reverse = TRUE)

# second
angles = seq(rad2deg(end)+3,rad2deg(end)+45+8,length.out = n+1)[1:(n)]
xx = cos(deg2rad(angles))*(lineSeq+0.6)
yy = sin(deg2rad(angles))*(lineSeq+0.6)
lines(xx, yy, xpd = NA)
arrow_angle = max(angles)-2.8
polygon(x = c(cos(deg2rad(arrow_angle))*(lineSeq+0.5), cos(deg2rad(arrow_angle))*(lineSeq+0.7), cos(deg2rad(max(angles)))*(lineSeq+0.6), cos(deg2rad(arrow_angle))*(lineSeq+0.5)),
        y = c(sin(deg2rad(arrow_angle))*(lineSeq+0.5), sin(deg2rad(arrow_angle))*(lineSeq+0.7), sin(deg2rad(max(angles)))*(lineSeq+0.6), sin(deg2rad(arrow_angle))*(lineSeq+0.5)),col = "black", xpd = NA)

# C


effects = coef(results$best1)[[1]][,-1]
within = apply(effects, 1, function(e) which.max(abs(e)))
max_eff = sapply(1:length(within), function(i) effects[i,within[i]])
effect_comb = cbind(max_eff, within)
effect_comb_ind = order(within, max_eff)

sigma = tuned_sigma[effect_comb_ind, effect_comb_ind]
sigmas = sigma[upper.tri(sigma)]
upper = order(sigmas, decreasing = TRUE)[1:number]
lower = order(sigmas, decreasing = FALSE)[1:number]
cuts = cut(sigmas, breaks = seq(-1,1,length.out = 21))
to_plot = 1:length(sigmas) %in% upper | 1:length(sigmas) %in% lower
levels(cuts) = viridis::viridis(20)
cuts = as.character(cuts)
n = ncol(sigma)
lineSeq = 3.5
nseg = 100
plot(NULL, NULL, xlim = c(-5,5), ylim =c(-5,5),pty="s", axes = F, xlab = "", ylab = "")
text(x = -5, y = 5.53, pos = 3, xpd = NA, labels = "C", font = 2, cex = 1.5)

xx = lineSeq*cos( seq(0,2*pi, length.out=nseg) )
yy = lineSeq*sin( seq(0,2*pi, length.out=nseg) )
polygon(xx,yy, col= "white", border = "black", lty = 1, lwd = 1)
angles = seq(0,360,length.out = n+1)[1:(n)]
xx = cos(deg2rad(angles))*lineSeq
yy = sin(deg2rad(angles))*lineSeq

counter = 1
coords = cbind(xx, yy, angles)
for(i in 1:n) {
  for(j in i:n){
    if(i!=j) {
      if(to_plot[counter]) add_curve(coords[i,], coords[j,], col = cuts[counter], n = 5, species = TRUE, lineSeq = lineSeq, lwd = 1.3)
      counter = counter + 1
    }
  }
}
cols = RColorBrewer::brewer.pal(8,"Dark2")
lineSeq = 4.2
d = 1.0

effect_comb2 = effect_comb
effect_comb2[,1] = ff(effect_comb[,1])
effect_comb2 = cbind(effect_comb2, effect_comb[,1])
for(i in sort(unique(within))) {
  sub = coords[within[effect_comb_ind] == i,]
  sub_eff = effect_comb2[effect_comb_ind,][effect_comb2[effect_comb_ind,2] == i,]
  from = sub[1,3]
  to = sub[nrow(sub),3]

  
  x = c((lineSeq+d*(sub_eff[,1]))*cos(deg2rad(sub[,3]) ), 
        rev((lineSeq+d/2)*cos(deg2rad(sub[,3]))))
  
  y = c((lineSeq+d*(sub_eff[,1]))*sin(deg2rad(sub[,3])),
        rev((lineSeq+d/2)*sin(deg2rad(sub[,3]))))
  
  names(env2)
  angleName = (from+to)/2
  if(angleName > 180) reverse = TRUE
  else reverse = FALSE
  curve_text(angleName, label = names(env2)[i],reverse = reverse,lineSeq = 4.1, middle = TRUE, extend = 1.17, col = cols[i])
  #y
  polygon(x, y, xpd = NA,col = cols[i])
  text(srt = 0, 
         x = (lineSeq+0.1+d)*cos(deg2rad(sub[1,3]+4)), 
         y =  (lineSeq+0.1+d)*sin(deg2rad(sub[1,3]+4)), 
         xpd = NA, labels = round(min(sub_eff[,3]), 2), col = cols[i], cex = 0.9)
  
   text(srt = 0, 
         x = (lineSeq+0.1+d)*cos(deg2rad(sub[nrow(sub),3]-4)), 
         y =  (lineSeq+0.1+d)*sin(deg2rad(sub[nrow(sub),3]-4)), 
         xpd = NA, labels = round(max(sub_eff[,3]), 2), col = cols[i], cex = 0.9)
}
rec_cols = viridis::viridis(11)

x = seq(3,5, length.out = 12)
for(i in 1:length(rec_cols)){
  rect(xleft = x[i], xright = x[i+1], ybottom = -5, ytop = -5+diff(x)[1], col = rec_cols[i], xpd = NA, border = NA)
}
text(x[1],-5.2, labels = -1)
text(x[11],-5.2, labels = +1)
text(mean(x),-5.4, labels = "correlation", xpd = NA)

dev.off()
```




## Process based simulation
### Figure S10
```{r}
plot_truth = function(scen, lwd = 2.0, col = "darkred") {
  ranges = par("usr")
  x = ranges[1:2]
  y = ranges[3:4]
  d1 = (-x[1]+x[2])/11
  if(scen == "A") return(diag(1.0, 12))
  if(scen == "B") {
    for(i in 0:10) {
      rect(x[1]+i*d1,x[1]+(i)*d1 , x[1]+(i+1)*d1, x[1]+(i+1)*d1 , lwd = lwd, border = col)
    }
    true = diag(1.0, 12)
    for(i in 2:12) {
      true[i, i-1] = 1.0
    }
    return(true)
  }
  if(scen == "C") {
    for(i in 6:11) {
      rect(x[1]+i*d1,x[1]+(i)*d1 , x[1]+(i+1)*d1, x[1]+(i+1)*d1  , lwd = lwd, border = col)
    }
      true = diag(1.0, 12)
      for(i in 8:12) {
        true[i, i-1] = 1.0
      }
      return(true)
  }
  if(scen == "D") {
    
    for(i in 0:10) {
      if(!i %in% c(3, 7)) rect(x[1]+i*d1,x[1]+(i)*d1 , x[1]+(i+1)*d1, x[1]+(i+1)*d1 , lwd = lwd, border = col, lty = )
    }
    true = diag(1.0, 12)
    for(i in 2:12) {
      if(!i %in% c(5, 9)) true[i, i-1] = 1.0
    }
    return(true)
  }
}
results_1 = readRDS("results/7_process_based_sjSDM_BC.RDS")
results_2 = readRDS("results/7_process_based_gllvm_Hmsc.RDS")
sims = readRDS("data_process_based.RDS")
cols = colorRampPalette(c( "white", "#24526e"))(39)
results = lapply(1:7, function(i) { lapply(1:5, function(j) c(results_1[[i]][[j]], results_2[[i]][[j]]))})

r2 = function(a, b) mean(sqrt((a-b)^2))


svg("figures/Fig_S10.svg", width = 10, height = 12)
par(mfrow = c(7, 4), mar = c(1, 1, 1,1), oma = c(3, 3, 3, 2))
scen = c("A", "A", "B", "B", "C", "A", "D")
for(i in 1:7) {
  for(j in 1:4) {
      tmp = matrix(0.0, 12, 12)
      for(k in 1:5) {
        c1 = results[[i]][[k]][[j]]$correlation
        if(sum(diag(c1)) < 0.1) diag(c1) = 1.0
        c1= abs(cov2cor(c1))
        c1[upper.tri(c1)] = 0.0
        tmp = tmp + c1
        diag(tmp) = 0.0
      }
      tmp = tmp/5
      image(tmp[2:12,1:11 ], col = cols, breaks = seq(0, 1, length.out=40), xaxt="n", yaxt="n", xaxs="i", yaxs = "i", useRaster = TRUE)
      true = plot_truth(scen[i])
      rmse = r2( true[lower.tri(true, diag = FALSE)], tmp[lower.tri(tmp, diag = FALSE)])
      text(-0.04, pos = 4, y = 0.5, labels = paste0("RMSE: ", round(rmse, 3)))
      
      if(j == 1) text(x = -0.1, y = 1.1, pos = 2, labels = LETTERS[i], xpd = NA, cex = 1.2, font = 2)
      if(i == 1) text(x = 0.5, y = 1.1, pos = 3, labels = c("sjSDM", "BayesComm", "Hmsc", "gllvm")[j], xpd = NA, font = 2)
      
      if(j == 4) {
          text(y = seq(from = 0, to = 1, length.out = 11), pos = 4,x = 1.07, xpd = NA, labels = paste0("sp", 2:12))
      }
      
      if(i == 7) {
          text(x = seq(from = 0, to = 1, length.out = 11)-0.01, y = -0.2, xpd = NA, 
               srt = 45, labels = paste0("sp", 1:11))
      }
  }
}
dev.off()

```

