---
title: "R Notebook"
output: html_notebook
---

1. Packages
```{r pacotes}
library(dplyr)
library(sf)
library(fitdistrplus)
library(piecewiseSEM)
library(lme4)
library(vegan)
library(MuMIn)
library(DHARMa)
library("plyr")
library(png)
library(gridExtra)
library(grid)
#library(grid)
library(ggplot2)
library(lattice)
library(ade4)
library(adespatial)
library(adegraphics)
library(spdep)
library(maptools)
library(sp)
library(plotrix)
library(DiagrammeR)
library(rbokeh)
library(htmlwidgets)
library(webshot)
library(sjPlot)
#install_phantomjs(version = "2.1.1",  baseURL = "https://github.com/wch/webshot/releases/download/v0.3.1/",  force = FALSE)
library(effects)

```

02. Database composed by datasets with standardized sampling effort among all sampling units
```{r}

data.glmm.ec<-list()

for ( i in 1:5){
 study.sel<-unique(data.glmm[[i]]$study_id)[tapply(data.glmm[[i]]$effort, data.glmm[[i]]$study_id, sd)==0]

data.glmm.ec[[i]]<-data.glmm[[i]][data.glmm[[i]]$study_id %in% study.sel,] 
}
names(data.glmm.ec)<-names(data.glmm)

```


_____________________________________________________________________________________________

                                    Richness
_____________________________________________________________________________________________


1. GLMMS update
    lapply does not work because change the name of data
    random slope of PA removed due to singularity
    
1.1. Complete 4k

```{r}
ec.rich.complete.4k<-list()

  for (i in 1:6){
    ec.rich.complete.4k[[i]]<-update(rich.complete.4k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.rich.complete.4k[[7]]<-update(rich.complete.4k[[1]]$np, .~. - (pa_4k | study_id), data=data.glmm.ec$complete)
  names(ec.rich.complete.4k)<-names(rich.complete.4k[[1]])


ec.rich.complete.4k.sf<-list()
  for (i in 1:6){
    ec.rich.complete.4k.sf[[i]]<-update(rich.complete.4k$glmm.sf[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.rich.complete.4k.sf[[7]]<-update(rich.complete.4k$glmm.sf$np,data=data.glmm.ec$complete,   .~. - (pa_4k | study_id))
  names(ec.rich.complete.4k.sf)<-names(rich.complete.4k[[1]])
  
```

1.2. Complete 2k

```{r}
ec.rich.complete.2k<-list()

  for (i in 1:6){
    ec.rich.complete.2k[[i]]<-update(rich.complete.2k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.rich.complete.2k[[7]]<-update(rich.complete.2k[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$complete)
  names(ec.rich.complete.2k)<-names(rich.complete.2k[[1]])


ec.rich.complete.2k.sf<-list()
  for (i in 1:6){
    ec.rich.complete.2k.sf[[i]]<-update(rich.complete.2k$glmm.sf[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.rich.complete.2k.sf[[7]]<-update(rich.complete.2k$glmm.sf$np,data=data.glmm.ec$complete,   .~. - (pa_2k | study_id))
  names(ec.rich.complete.2k.sf)<-names(rich.complete.2k[[1]])
  
```

1.3. preserved

```{r}
ec.rich.preserved<-list()

  for (i in 1:6){
    ec.rich.preserved[[i]]<-update(rich.preserved$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$preserved) 
  }
    ec.rich.preserved[[7]]<-update(rich.preserved[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$preserved)
  names(ec.rich.preserved)<-names(rich.preserved[[1]])
  
```


1.4. intermediate

```{r}
ec.rich.intermediated<-list()

  for (i in 1:6){
    ec.rich.intermediated[[i]]<-update(rich.intermediated$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$intermediated) 
  }
    ec.rich.intermediated[[7]]<-update(rich.intermediated[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$intermediated)
  names(ec.rich.intermediated)<-names(rich.intermediated[[1]])
  
```


1.5. relictual 2k

```{r}
ec.rich.relictual.2k<-list()

  for (i in 1:6){
    ec.rich.relictual.2k[[i]]<-update(rich.relictual.2k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$relictual2k) 
  }
    ec.rich.relictual.2k[[7]]<-update(rich.relictual.2k[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$relictual2k)
  names(ec.rich.relictual.2k)<-names(rich.relictual.2k[[1]])
  
```

1.6. relictual 1k

```{r}
ec.rich.relictual.1k<-list()

  for (i in 1:6){
    ec.rich.relictual.1k[[i]]<-update(rich.relictual.1k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$relictual1k) 
  }
    ec.rich.relictual.1k[[7]]<-update(rich.relictual.1k[[1]]$np, .~. - (pa_1k | study_id), data=data.glmm.ec$relictual1k)
  names(ec.rich.relictual.1k)<-names(rich.relictual.1k[[1]])
  
```


2. Creating directory and files

```{r}
grad.name<-c("Relictual_1k","Relictual_2k","Intermediated","Preserved","Complete_2k_base", "Complete_2k_sf", "Complete_4k_base", "Complete_4k_sf")

dir.rich<- paste( "rich_", grad.name, "/", sep="")

lapply(dir.rich, dir_create)

```


3. pSEM

3.1. relictual 1k

```{r}
ec.rich.base.relic.1k<-psem.3H(g1.lands= ec.rich.relictual.1k$for.land,
                             g2.lands= ec.rich.relictual.1k$gen.land,
                             g1.sps=   ec.rich.relictual.1k$for.sps,
                             g2.sps=   ec.rich.relictual.1k$gen.sps,
                             frag=     ec.rich.relictual.1k$np,
                             data=data.glmm.ec$relictual1k,
                             dir=dir.rich[1], #change directory
                             direction=c("rich_gen<-np_1k", "rich_for<-np_1k"))
```


3.2. relictual 2k

```{r}
ec.rich.base.relic.2k<-psem.3H(g1.lands= ec.rich.relictual.2k$for.land,
                             g2.lands= ec.rich.relictual.2k$gen.land,
                             g1.sps=   ec.rich.relictual.2k$for.sps,
                             g2.sps=   ec.rich.relictual.2k$gen.sps,
                             frag=     ec.rich.relictual.2k$np,
                             data=data.glmm.ec$relictual2k,
                             dir=dir.rich[2], #change directory
                             direction=c("rich_gen<-np_2k", "rich_for<-np_2k"))
```


3.3. Intermediated

```{r}
ec.rich.base.intermediated<-psem.3H(g1.lands= ec.rich.intermediated$for.land,
                             g2.lands= ec.rich.intermediated$gen.land,
                             g1.sps=   ec.rich.intermediated$for.sps,
                             g2.sps=   ec.rich.intermediated$gen.sps,
                             frag=     ec.rich.intermediated$np,
                             data=data.glmm.ec$intermediated,
                             dir=dir.rich[3], #change directory
                             direction=c("rich_gen<-np_2k", "rich_for<-np_2k"))
```


3.4. Preserved

```{r}
ec.rich.base.preserved<-psem.3H(g1.lands= ec.rich.preserved$for.land,
                             g2.lands= ec.rich.preserved$gen.land,
                             g1.sps=   ec.rich.preserved$for.sps,
                             g2.sps=   ec.rich.preserved$gen.sps,
                             frag=     ec.rich.preserved$np,
                             data=data.glmm.ec$preserved,
                             dir=dir.rich[4], #change directory
                             direction=c("rich_gen<-np_2k", "rich_for<-np_2k"))
```


3.5. Complete 2k

a. without autocorrelation control

```{r}
ec.rich.base.complete.2k<-psem.3H(g1.lands=  ec.rich.complete.2k$for.land,
                             g2.lands= ec.rich.complete.2k$gen.land,
                             g1.sps=   ec.rich.complete.2k$for.sps,
                             g2.sps=   ec.rich.complete.2k$gen.sps,
                             frag=     ec.rich.complete.2k$np,
                             data=data.glmm.ec$complete,
                             dir=dir.rich[5], #change directory
                             direction=c("rich_gen<-np_2k", "rich_for<-np_2k"))
```


b. with autocorrelation control (Moran Eigenvector Maps - MEM)

```{r}
ec.rich.sf.complete.2k<-psem.3H(g1.lands=  ec.rich.complete.2k$for.land,
                             g2.lands= ec.rich.complete.2k$gen.land,
                             g1.sps=   ec.rich.complete.2k$for.sps,
                             g2.sps=   ec.rich.complete.2k$gen.sps,
                             frag=     ec.rich.complete.2k.sf$np,
                             data=data.glmm.ec$complete,
                             dir=dir.rich[6], #change directory
                             direction=c("rich_gen<-np_2k", "rich_for<-np_2k"))
```


3.6. Complete 4k

a. without autocorrelation control

```{r}
ec.rich.base.complete.4k<-psem.3H(g1.lands=  ec.rich.complete.4k$for.land,
                             g2.lands= ec.rich.complete.4k$gen.land,
                             g1.sps=   ec.rich.complete.4k$for.sps,
                             g2.sps=   ec.rich.complete.4k$gen.sps,
                             frag=     ec.rich.complete.4k$np,
                             data=data.glmm.ec$complete,  dir=dir.rich[7], #change directory
                             direction=c("rich_gen<-np_4k", "rich_for<-np_4k"))
```

b. with autocorrelation control (Moran Eigenvector Maps - MEM)

```{r}
ec.rich.sf.complete.4k<-psem.3H(g1.lands=  ec.rich.complete.4k$for.land,
                             g2.lands= ec.rich.complete.4k$gen.land,
                             g1.sps=   ec.rich.complete.4k$for.sps,
                             g2.sps=   ec.rich.complete.4k$gen.sps,
                             frag=     ec.rich.complete.4k.sf$np,
                             data=data.glmm.ec$complete,
                             dir=dir.rich[8], #change directory
                             direction=c("rich_gen<-np_4k", "rich_for<-np_4k"))
```


4. Table

```{r}
ec.rich.psem.sel<-rbind(ec.rich.base.complete.4k$aic.tab, ec.rich.sf.complete.4k$aic.tab,
                     ec.rich.base.complete.2k$aic.tab, ec.rich.sf.complete.2k$aic.tab,
                     ec.rich.base.preserved$aic.tab, ec.rich.base.intermediated$aic.tab,
                     ec.rich.base.relic.2k$aic.tab, ec.rich.base.relic.1k$aic.tab)
write.table(ec.rich.psem.sel, "ec_rich_psem_sel.txt")
View(ec.rich.psem.sel)

ec.rich.psem.coef<-rbind(ec.rich.base.complete.4k$psem$env.fil$psem.summary$coefficients,
                         ec.rich.sf.complete.4k$psem$env.fil$psem.summary$coefficients,
                         ec.rich.base.complete.2k$psem$env.fil$psem.summary$coefficients,
                         ec.rich.sf.complete.2k$psem$env.fil$psem.summary$coefficients,
                         ec.rich.base.preserved$psem$comp.Out$psem.summary$coefficients,
                         ec.rich.base.intermediated$psem$comp.lib$psem.summary$coefficients,
                         ec.rich.base.relic.2k$psem$comp.lib$psem.summary$coefficients,
                         ec.rich.base.relic.1k$psem$comp.lib$psem.summary$coefficients)
write.table(ec.rich.psem.coef, "ec_rich_psem_coef.txt")
View(ec.rich.psem.coef)

ec.rich.psem.R2<-rbind(ec.rich.base.complete.4k$psem$env.fil$psem.summary$R2, 
                         ec.rich.sf.complete.4k$psem$env.fil$psem.summary$R2,
                         ec.rich.base.complete.2k$psem$env.fil$psem.summary$R2,
                         ec.rich.sf.complete.2k$psem$env.fil$psem.summary$R2,
                         ec.rich.base.preserved$psem$comp.Out$psem.summary$R2,
                         ec.rich.base.intermediated$psem$comp.lib$psem.summary$R2,
                         ec.rich.base.relic.2k$psem$comp.lib$psem.summary$R2,
                         ec.rich.base.relic.1k$psem$comp.lib$psem.summary$R2)
write.table(ec.rich.psem.R2, "ec_rich_psem_R2.txt")
View(ec.rich.psem.R2)

```


_____________________________________________________________________________________________

                                    Abundance
_____________________________________________________________________________________________


1. GLMMS update

1.1. Complete 4k

```{r}
ec.abund.complete.4k<-list()

  for (i in 1:6){
    ec.abund.complete.4k[[i]]<-update(abund.complete.4k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.abund.complete.4k[[7]]<-update(abund.complete.4k[[1]]$np, .~. - (pa_4k | study_id), data=data.glmm.ec$complete)
  names(ec.abund.complete.4k)<-names(abund.complete.4k[[1]])


ec.abund.complete.4k.sf<-list()
  for (i in 1:6){
    ec.abund.complete.4k.sf[[i]]<-update(abund.complete.4k$glmm.sf[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.abund.complete.4k.sf[[7]]<-update(abund.complete.4k$glmm.sf$np,data=data.glmm.ec$complete,   .~. - (pa_4k | study_id))
  names(ec.abund.complete.4k.sf)<-names(abund.complete.4k[[1]])
  
```


1.2. Complete 2k

```{r}
ec.abund.complete.2k<-list()

  for (i in 1:6){
    ec.abund.complete.2k[[i]]<-update(abund.complete.2k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.abund.complete.2k[[7]]<-update(abund.complete.2k[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$complete)
  names(ec.abund.complete.2k)<-names(abund.complete.2k[[1]])


ec.abund.complete.2k.sf<-list()
  for (i in 1:6){
    ec.abund.complete.2k.sf[[i]]<-update(abund.complete.2k$glmm.sf[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$complete) 
  }
    ec.abund.complete.2k.sf[[7]]<-update(abund.complete.2k$glmm.sf$np,data=data.glmm.ec$complete,   .~. - (pa_2k | study_id))
  names(ec.abund.complete.2k.sf)<-names(abund.complete.2k[[1]])
  
```

1.3. preserved

```{r}
ec.abund.preserved<-list()

  for (i in 1:6){
    ec.abund.preserved[[i]]<-update(abund.preserved$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$preserved) 
  }
    ec.abund.preserved[[7]]<-update(abund.preserved[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$preserved)
  names(ec.abund.preserved)<-names(abund.preserved[[1]])
  
```

1.4. intermediated

```{r}
ec.abund.intermediated<-list()

  for (i in 1:6){
    ec.abund.intermediated[[i]]<-update(abund.intermediated$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$intermediated) 
  }
    ec.abund.intermediated[[7]]<-update(abund.intermediated[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$intermediated)
  names(ec.abund.intermediated)<-names(abund.intermediated[[1]])
  
```

1.5. relictual 2k

```{r}
ec.abund.relictual.2k<-list()

  for (i in 1:6){
    ec.abund.relictual.2k[[i]]<-update(abund.relictual.2k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$relictual2k) 
  }
    ec.abund.relictual.2k[[7]]<-update(abund.relictual.2k[[1]]$np, .~. - (pa_2k | study_id), data=data.glmm.ec$relictual2k)
  names(ec.abund.relictual.2k)<-names(abund.relictual.2k[[1]])
  
```

1.6. relictual 1k

```{r}
ec.abund.relictual.1k<-list()

  for (i in 1:6){
    ec.abund.relictual.1k[[i]]<-update(abund.relictual.1k$glmm.base[[i]], .~. -effort - (effort | study_id), data=data.glmm.ec$relictual1k) 
  }
    ec.abund.relictual.1k[[7]]<-update(abund.relictual.1k[[1]]$np, .~. - (pa_1k | study_id), data=data.glmm.ec$relictual1k)
  names(ec.abund.relictual.1k)<-names(abund.relictual.1k[[1]])
  
```


02. creating directory and files

```{r}
grad.name<-c("Relictual_1k","Relictual_2k","Intermediated","Preserved","Complete_2k_base", "Complete_2k_sf", "Complete_4k_base", "Complete_4k_sf")

dir.abund<- paste("abund_", grad.name, "/", sep="")

lapply(dir.abund, dir_create)

```


03. pSEM

3.1. relictual 1k

```{r}
ec.abund.base.relic.1k<-psem.3H(g1.lands= ec.abund.relictual.1k$for.land,
                             g2.lands= ec.abund.relictual.1k$gen.land,
                             g1.sps=   ec.abund.relictual.1k$for.sps,
                             g2.sps=   ec.abund.relictual.1k$gen.sps,
                             frag=     ec.abund.relictual.1k$np,
                             data=data.glmm.ec$relictual1k,
                             dir=dir.abund[1], #change directory
                             direction=c("abund_gen<-np_1k", "abund_for<-np_1k"))
```

3.2. relictual 2k

```{r}
ec.abund.base.relic.2k<-psem.3H(g1.lands= ec.abund.relictual.2k$for.land,
                             g2.lands= ec.abund.relictual.2k$gen.land,
                             g1.sps=   ec.abund.relictual.2k$for.sps,
                             g2.sps=   ec.abund.relictual.2k$gen.sps,
                             frag=     ec.abund.relictual.2k$np,
                             data=data.glmm.ec$relictual2k,
                             dir=dir.abund[2], #change directory
                             direction=c("abund_gen<-np_2k", "abund_for<-np_2k"))
```

3.3. Intermediated

```{r}
ec.abund.base.intermediated<-psem.3H(g1.lands= ec.abund.intermediated$for.land,
                             g2.lands= ec.abund.intermediated$gen.land,
                             g1.sps=   ec.abund.intermediated$for.sps,
                             g2.sps=   ec.abund.intermediated$gen.sps,
                             frag=     ec.abund.intermediated$np,
                             data=data.glmm.ec$intermediated,
                             dir=dir.abund[3], #change directory
                             direction=c("abund_gen<-np_2k", "abund_for<-np_2k"))
```

3.4. Preserved

```{r}
ec.abund.base.preserved<-psem.3H(g1.lands= ec.abund.preserved$for.land,
                             g2.lands= ec.abund.preserved$gen.land,
                             g1.sps=   ec.abund.preserved$for.sps,
                             g2.sps=   ec.abund.preserved$gen.sps,
                             frag=     ec.abund.preserved$np,
                             data=data.glmm.ec$preserved,
                             dir=dir.abund[4], #change directory
                             direction=c("abund_gen<-np_2k", "abund_for<-np_2k"))
```

3.5. Complete 2k

a. without autocorrelation control

```{r}
ec.abund.base.complete.2k<-psem.3H(g1.lands=  ec.abund.complete.2k$for.land,
                             g2.lands= ec.abund.complete.2k$gen.land,
                             g1.sps=   ec.abund.complete.2k$for.sps,
                             g2.sps=   ec.abund.complete.2k$gen.sps,
                             frag=     ec.abund.complete.2k$np,
                             data=data.glmm.ec$complete,
                             dir=dir.abund[5], #change directory
                             direction=c("abund_gen<-np_2k", "abund_for<-np_2k"))
```

b. with autocorrelation control (Moran Eigenvector Maps - MEM)

```{r}
ec.abund.sf.complete.2k<-psem.3H(g1.lands=  ec.abund.complete.2k$for.land,
                             g2.lands= ec.abund.complete.2k$gen.land,
                             g1.sps=   ec.abund.complete.2k$for.sps,
                             g2.sps=   ec.abund.complete.2k$gen.sps,
                             frag=     ec.abund.complete.2k.sf$np,
                             data=data.glmm.ec$complete,
                             dir=dir.abund[6], #change directory
                             direction=c("abund_gen<-np_2k", "abund_for<-np_2k"))
```

3.6. Complete 4k

a. without autocorrelation control

```{r}
ec.abund.base.complete.4k<-psem.3H(g1.lands=  ec.abund.complete.4k$for.land,
                             g2.lands= ec.abund.complete.4k$gen.land,
                             g1.sps=   ec.abund.complete.4k$for.sps,
                             g2.sps=   ec.abund.complete.4k$gen.sps,
                             frag=     ec.abund.complete.4k$np,
                             data=data.glmm.ec$complete,
                             dir=dir.abund[7], #change directory
                             direction=c("abund_gen<-np_4k", "abund_for<-np_4k"))
```

b. with autocorrelation control (Moran Eigenvector Maps - MEM)

```{r}
ec.abund.sf.complete.4k<-psem.3H(g1.lands=  ec.abund.complete.4k$for.land,
                             g2.lands= ec.abund.complete.4k$gen.land,
                             g1.sps=   ec.abund.complete.4k$for.sps,
                             g2.sps=   ec.abund.complete.4k$gen.sps,
                             frag=     ec.abund.complete.4k.sf$np,
                             data=data.glmm.ec$complete,
                             dir=dir.abund[8], #change directory
                             direction=c("abund_gen<-np_4k", "abund_for<-np_4k"))
```


4. Table

```{r}
ec.abund.psem.sel<-rbind(ec.abund.base.complete.4k$aic.tab, ec.abund.sf.complete.4k$aic.tab,
                     ec.abund.base.complete.2k$aic.tab, ec.abund.sf.complete.2k$aic.tab,
                     ec.abund.base.preserved$aic.tab, ec.abund.base.intermediated$aic.tab,
                     #ec.abund.base.relic.2k$aic.tab, 
                     ec.abund.base.relic.1k$aic.tab)
write.table(ec.abund.psem.sel, "ec_abund_psem_sel.txt")
View(ec.abund.psem.sel)

ec.abund.psem.coef<-rbind(ec.abund.base.complete.4k$psem$env.fil$psem.summary$coefficients,
                      ec.abund.sf.complete.4k$psem$env.fil$psem.summary$coefficients,
                      ec.abund.base.complete.2k$psem$env.fil$psem.summary$coefficients,
                      ec.abund.sf.complete.2k$psem$env.fil$psem.summary$coefficients,
                      ec.abund.base.preserved$psem$comp.Out$psem.summary$coefficients,
                      ec.abund.base.intermediated$psem$comp.lib$psem.summary$coefficients,
                      #ec.abund.base.relic.2k$psem$comp.lib$psem.summary$coefficients,
                      ec.abund.base.relic.1k$psem$comp.Out$psem.summary$coefficients)
write.table(ec.abund.psem.coef, "ec_abund_psem_coef.txt")
View(ec.abund.psem.coef)

ec.abund.psem.R2<-rbind(ec.abund.base.complete.4k$psem$env.fil$psem.summary$R2,
                      ec.abund.sf.complete.4k$psem$env.fil$psem.summary$R2,
                      ec.abund.base.complete.2k$psem$env.fil$psem.summary$R2,
                      ec.abund.sf.complete.2k$psem$env.fil$psem.summary$R2,
                      ec.abund.base.preserved$psem$comp.Out$psem.summary$R2,
                      ec.abund.base.intermediated$psem$comp.lib$psem.summary$R2,
                      #ec.abund.base.relic.2k$psem$comp.lib$psem.summary$R2,
                      ec.abund.base.relic.1k$psem$comp.Out$psem.summary$R2)
write.table(ec.abund.psem.R2, "ec_abund_psem_R2.txt")
View(ec.abund.psem.R2)

```

