Main_analyse.Rmd

title: "Main_analyse.Rmd"
author: "Isabelle Boulangeat"
date: "23/01/2020"
output: 
  html_document:
    keep_md: yes
    variant: markdown_github
editor_options: 
  chunk_output_type: console
always_allow_html: true

library(knitr)
library(kableExtra)
knitr::opts_chunk$set(echo = TRUE)
library(Jmisc)
library(tidyr)
sourceAll("R_fct")

Load data

• Load climate
• Load species cover data
• Filter out treatment plots
• filter out un-repeated plots
• filter out years > 2014 because of missing in climatologies

DATA = load.datasets()
str(DATA$cover)
str(DATA$climato)

mostRecordedSP = DATA$cover %>% group_by(Name) %>% summarise(nobs = length(Name)) %>% arrange(desc(nobs))
selection = as.vector(t(mostRecordedSP[1:30, "Name"]))
selection = selection[-which(selection %in% c("Arctagrostis latifolia", "Persicaria vivipara", "Carex aquatilis", "Cassiope tetragona", "Luzula confusa", "Luzula nivalis", "Salix arctica", "Poa arctica", "Eriophorum angustifolium", "Petasites frigidus", "Vaccinium uliginosum", "Dryas integrifolia", "Stellaria longipes", "Saxifraga cernua", "Festuca brachyphylla"))]

N.B. I have quickly selected 14 species among the 30 having the most observations, with no problem with zero-days growing season in the previous years before an observation, and no problem in the dataset... selection procedure to discuss

Prepare species specific datasets

Step 1: create_dataset:
- only subsites where the species have been seen (then ok for N=1)
- extract plot info
- build transition table
- merge with plot info
- calculate interval (years)
- tag transition types
- combine with climate (5 years earlier average)
- merge all
Step 2: adjustments
- Add anomalies
- remove intervals>=10 years
- select variables c("snowDays", "GSL", "GST", "SoilMoist", "TotalCover")
- modify SoilMoist levels (-1, 0, 1)
Step 3: prepare dataset to fit
- separate colonisation and extinction datasets
- scale data (not SoilMoist)
- prepare data for LaplacesDemon fit (optional)
- add event column for glm fit

spDATA.list = lapply(selection, function(SPECIES){
  print(SPECIES)
  spDATA = dat.species(SPECIES, DATA)
  return(spDATA)
} )
names(spDATA.list) = selection
spDATA.check = lapply(spDATA.list, check.nas)

Visualisation of transition tables

lapply(spDATA.list, function(x) table(x$spdat$trType))

Here we can see how many transitions we have for each species.

Fit models

spDATA.fit = lapply(spDATA.list, fit.species)

Summary of model evaluation

eval.models = lapply(spDATA.fit, function(glm.fit){
  eval.colo = eval.fit(glm.fit$colo, glm.fit$dat$colo, "colo.event")[-1]
  coeffs.colo = sort(abs(glm.fit$colo$coefficients[-1]), dec=T)
  vars.colo = paste(names(coeffs.colo), collapse=";")
  eval.ext = eval.fit(glm.fit$ext, glm.fit$dat$ext, "ext.event")[-1]
  coeffs.ext = sort(abs(glm.fit$ext$coefficients[-1]), dec=T)
  vars.ext = paste(names(coeffs.ext), collapse=";")
  return(list(colo.R2 = eval.colo$R2, colo.AUC = eval.colo$AUC, ext.R2 = eval.ext$R2, ext.AUC = eval.ext$AUC, vars.colo = vars.colo, vars.ext = vars.ext, coeffs.colo = coeffs.colo, coeffs.ext=coeffs.ext))
})

res.tab.eval = as.data.frame(do.call(rbind, lapply(eval.models, function(x)unlist(x[1:4]))))
res.tab.vars = as.data.frame(do.call(rbind, lapply(eval.models, function(x)unlist(x[5:6]))))
  
kable(res.tab.eval, caption = "model evaluation", digits = 2) %>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), full_width=FALSE)
# kable(res.tab.vars, caption = "selected variables")%>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), font_size = 10)

N.B. Some models do not fit well. The variables explain more or less depending on species but also on processes (extinction and colonisation).

# x = spDATA.fit[["Betula nana"]]
# layout(matrix(c(1,1), ncol = 2))
# plot_model(list(x$colo, x$ext), transform = NULL, show.values=TRUE, title = paste(x$dat$sp, "colonisation"))
# plot_model(x$ext, transform = NULL, show.values=TRUE, title = paste(x$dat$sp, "extinction"))

Selected variables and coefficient, sumnmary table

coeffs.all.list = lapply(spDATA.fit, function(x){
  coe.colo = data.frame(summary(x$colo)$coefficients)
  coe.colo[,"sp"] = x$dat$sp
  coe.colo[,"var"] = rownames(coe.colo)
  coe.colo[,"process"] = "colonisation"
  coe.ext = data.frame(summary(x$ext)$coefficients)
  coe.ext[,"sp"] = x$dat$sp
  coe.ext[,"var"] = rownames(coe.ext)
  coe.ext[,"process"] = "extinction"
  coe = rbind(coe.colo[-1,], coe.ext[-1,])
  return(coe)
  })
coeffs.all = do.call(rbind, coeffs.all.list)
# str(coeffs.all)

Number of selected variable per process, accross all species

data.frame(unclass(table(coeffs.all$var, coeffs.all$process)))  %>% 
  kable(.) %>%
  kable_styling("striped")

Variable effect coefficients

This gives a summary of selected variables with significant impact on colonisation or extinction and their estimated coefficients.

coeffs.all[which(coeffs.all$process == "colonisation" & coeffs.all$Pr...z..<0.1),] %>%
  ggplot(aes(sp, Estimate, fill = var)) +
  geom_col(position = "dodge") +
  theme_bw() + facet_wrap(~sp,scales = "free_x", ncol=4) +
  ggtitle("Colonisation process")
coeffs.all[which(coeffs.all$process == "extinction" & coeffs.all$Pr...z..<0.1),] %>%
  ggplot(aes(sp, Estimate, fill = var)) +
  geom_col(position = "dodge") +
  theme_bw() + facet_wrap(~sp,scales = "free_x", ncol=4) +
  ggtitle("Extinction process")

Lambda (1-ext/colo)

An overview of the performance of populations by site of by species

all.sites = unique(DATA$cover[,c("SemiUniquePLOT", "YEAR", "SUBSITE", "SoilMoist", "TotalCover")])
pDATA = proj.dat.shape(all.sites, DATA$climato)
selectVars = c("snowDays", "GSL", "GST", "SoilMoist", "TotalCover")
pDATA.all = na.omit(pDATA[,c(selectVars,"SemiUniquePLOT", "YEAR", "SUBSITE") ])

pred.all = lapply(spDATA.fit, predictions, pred.data = pDATA.all, selectVars=selectVars)
pred.tab = lapply(pred.all, function(x)data.frame(sp = x$sp, lambda = x$lambda, YEAR = x$YEAR, SUBSITE = x$SUBSITE))
pred.all.together = do.call(rbind, pred.tab)

pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
  ggplot(aes(YEAR, lambda)) +
  geom_point(shape = 1) +
  geom_smooth() +
  facet_wrap( ~ sp, ncol = 4)

pred.all.together$SITE =  unlist( lapply(as.character(pred.all.together$SUBSITE), function(x){
  strsplit(x, ":")[[1]][[1]]
}))

pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
  ggplot(aes(SITE, lambda, sp)) + 
  theme(axis.text.x = element_text(angle=90, hjust=1)) +
  geom_boxplot() + geom_hline(yintercept = 0) +
  facet_wrap( ~ sp,scales = "free_x", ncol = 7)

pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
  ggplot(aes(SITE, lambda, sp)) +
  theme(axis.text.x = element_text(angle=90, hjust=1)) +
  geom_boxplot() + ylim (-1, 1) + geom_hline(yintercept = 0) +
  facet_wrap( ~ sp,scales = "free_x", ncol = 7)

Simplified version focusing only on the sign of lambda

pred.all.together$lambdaSign = ifelse(pred.all.together$lambda>0, "+", "-")

pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
  ggplot(aes(x=SITE, y=lambdaSign, fill= lambdaSign )) +
  theme(axis.text.x = element_text(angle=90, hjust=1)) +
  geom_bar(stat="identity")  + 
  facet_wrap( ~ sp,scales = "free_x", ncol = 7)