Main_analyse.Rmd 7.73 KB
Newer Older
1 2 3 4
---
title: "Main_analyse.Rmd"
author: "Isabelle Boulangeat"
date: "23/01/2020"
5 6
output: 
  html_document:
7 8
    keep_md: yes
    variant: markdown_github
9 10
editor_options: 
  chunk_output_type: console
11
always_allow_html: true
12 13 14
---

```{r setup, include=FALSE}
15 16
library(knitr)
library(kableExtra)
17 18
knitr::opts_chunk$set(echo = TRUE)
library(Jmisc)
19
library(tidyr)
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
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

```{r load_data}
DATA = load.datasets()
str(DATA$cover)
str(DATA$climato)

mostRecordedSP = DATA$cover %>% group_by(Name) %>% summarise(nobs = length(Name)) %>% arrange(desc(nobs))
36 37
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"))]
38 39
```

40
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**
41 42 43

## Prepare species specific datasets
Step 1: `create_dataset`:
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
<br>- only subsites where the species have been seen (then ok for N=1)
<br>- extract plot info
<br>- build transition table
<br>- merge with plot info
<br>- calculate interval (years)
<br>- tag transition types
<br>- combine with climate (5 years earlier average)
<br>- merge all
<br>Step 2: adjustments
<br>- Add anomalies
<br>- remove intervals>=10 years
<br>- select variables c("snowDays", "GSL", "GST", "SoilMoist", "TotalCover")
<br>- modify SoilMoist levels (-1, 0, 1)
<br>Step 3: prepare dataset to fit
<br>- separate colonisation and extinction datasets
<br>- scale data (not SoilMoist)
<br>- prepare data for LaplacesDemon fit (optional)
<br>- add event column for glm fit

```{r species_specific_datasets, results='hide'}
64 65 66 67 68 69
spDATA.list = lapply(selection, function(SPECIES){
  print(SPECIES)
  spDATA = dat.species(SPECIES, DATA)
  return(spDATA)
} )
names(spDATA.list) = selection
70 71
spDATA.check = lapply(spDATA.list, check.nas)
```
72

73 74
## Visualisation of transition tables
```{r transition_tables}
75 76
lapply(spDATA.list, function(x) table(x$spdat$trType))
```
77
Here we can see how many transitions we have for each species.
78 79

## Fit models
80
```{r species_specific_models, results='hide', warning=FALSE}
81 82 83 84
spDATA.fit = lapply(spDATA.list, fit.species)
```

## Summary of model evaluation
85
```{r model_evaluation, results='asis'}
86 87
eval.models = lapply(spDATA.fit, function(glm.fit){
  eval.colo = eval.fit(glm.fit$colo, glm.fit$dat$colo, "colo.event")[-1]
88 89
  coeffs.colo = sort(abs(glm.fit$colo$coefficients[-1]), dec=T)
  vars.colo = paste(names(coeffs.colo), collapse=";")
90
  eval.ext = eval.fit(glm.fit$ext, glm.fit$dat$ext, "ext.event")[-1]
91 92 93
  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))
94
})
95 96 97 98 99

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)
100
# kable(res.tab.vars, caption = "selected variables")%>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), font_size = 10)
101
```
102
N.B. Some models do not fit well. The variables explain more or less depending on species but also on processes (extinction and colonisation).
103 104 105 106 107 108

```{r model_plots, include = FALSE}
# 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"))
109 110
```

111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
## Selected variables and coefficient, sumnmary table
```{r coefficients}
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)
126
# str(coeffs.all)
127 128
```

129
## Number of selected variable per process, accross all species
130 131 132 133 134 135
```{r summary_sel_vars}
data.frame(unclass(table(coeffs.all$var, coeffs.all$process)))  %>% 
  kable(.) %>%
  kable_styling("striped")
```

136
## Variable effect coefficients
137 138 139

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

140 141 142 143 144 145 146 147 148 149 150 151 152
```{r plot_coeffs, fig = TRUE, fig.width=10}
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")
```

153 154
## Lambda (1-ext/colo)
An overview of the performance of populations by site of by species
155

156
```{r lambda, fig-TRUE, fig.width=10, warning=FALSE, fig.height=7}
157 158 159 160 161 162
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)
163
pred.tab = lapply(pred.all, function(x)data.frame(sp = x$sp, lambda = x$lambda, YEAR = x$YEAR, SUBSITE = x$SUBSITE))
164 165
pred.all.together = do.call(rbind, pred.tab)

166 167 168 169 170 171 172 173 174 175 176 177 178 179
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) +
180
  facet_wrap( ~ sp,scales = "free_x", ncol = 7)
181 182 183 184 185 186 187

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)
```
188

189
Simplified version focusing only on the sign of lambda
190
```{r lambda_sign, fig=TRUE, fig.width=10, warning=FALSE}
191
pred.all.together$lambdaSign = ifelse(pred.all.together$lambda>0, "+", "-")
192

193 194 195 196 197
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)
198 199 200
```


201