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}
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
15 16
library(knitr)
library(kableExtra)
17 18
knitr::opts_chunk$set(echo = TRUE)
library(Jmisc)
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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`:
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
70 71
spDATA.check = lapply(spDATA.list, check.nas)
```
72

Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
80
```{r species_specific_models, results='hide', warning=FALSE}
81 82 83 84
spDATA.fit = lapply(spDATA.list, fit.species)
```

## Summary of model evaluation
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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]
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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]
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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
})
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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)
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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).
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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
```

Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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)
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
127 128
```

129
## Number of selected variable per process, accross all species
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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.

Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
155

156
```{r lambda, fig-TRUE, fig.width=10, warning=FALSE, fig.height=7}
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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))
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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)
```
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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, "+", "-")
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
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)
Boulangeat Isabelle's avatar
Boulangeat Isabelle committed
198 199 200
```


201