Commit a040a2b7 authored by Boulangeat Isabelle's avatar Boulangeat Isabelle

full analysis 14 species + README

parent 634fa733
......@@ -2,9 +2,13 @@
title: "Main_analyse.Rmd"
author: "Isabelle Boulangeat"
date: "23/01/2020"
output: html_document
output:
md_document:
variant: markdown_github
html_document:
editor_options:
chunk_output_type: console
always_allow_html: true
---
```{r setup, include=FALSE}
......@@ -29,11 +33,11 @@ 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:22, "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"))]
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 10 species among the 22 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**
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`:
......@@ -70,7 +74,7 @@ spDATA.check = lapply(spDATA.list, check.nas)
```{r 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
```{r species_specific_models, results='hide', warning=FALSE}
......@@ -95,6 +99,7 @@ res.tab.vars = as.data.frame(do.call(rbind, lapply(eval.models, function(x)unlis
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).
```{r model_plots, include = FALSE}
# x = spDATA.fit[["Betula nana"]]
......@@ -118,19 +123,20 @@ coeffs.all.list = lapply(spDATA.fit, function(x){
return(coe)
})
coeffs.all = do.call(rbind, coeffs.all.list)
str(coeffs.all)
# str(coeffs.all)
```
## Number of selection of variable per process
## Number of selected variable per process, accross all species
```{r summary_sel_vars}
data.frame(unclass(table(coeffs.all$var, coeffs.all$process))) %>%
mutate(
colonisation = cell_spec(colonisation, color = ifelse(colonisation > 5, "red", "black")),
extinction = cell_spec(extinction, color = ifelse(extinction > 5, "red", "black"))) %>%
kable(.) %>%
kable_styling("striped")
```
## Variables effect coefficients
This gives a summary of selected variables with significant impact on colonisation or extinction and their estimated coefficients.
```{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)) +
......@@ -144,8 +150,10 @@ coeffs.all[which(coeffs.all$process == "extinction" & coeffs.all$Pr...z..<0.1),]
ggtitle("Extinction process")
```
```{r lambda, include = FALSE}
## Lambda (1-ext/colo)
An overview of the performance of populations by site of by species
```{r lambda, fig-TRUE, fig.width=10}
all.sites = unique(DATA$cover[,c("SemiUniquePLOT", "YEAR", "SUBSITE", "SoilMoist", "TotalCover")])
pDATA = proj.dat.shape(all.sites, DATA$climato)
str(pDATA)
......@@ -154,21 +162,42 @@ 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))
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)
head(pred.all.together)
pred.all.together$YEAR = pDATA.all$YEAR
pred.all.together$SUBSITE = pDATA.all$SUBSITE
summary(pred.all.together)
summary(glm(lambda~YEAR, data = pred.all.together))
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 = 4)
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)
```
pred.all.together %>% group_by(SUBSITE, YEAR) %>%
summarise(lambda = mean(lambda)) %>%
ggplot(aes(YEAR, lambda, SUBSITE))+
geom_point(shape = 1)
Simplified version focusing only on the sign of lambda
```{r lambda_sign, fig=TRUE, fig.width=10}
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)
```
......
This diff is collapsed.
This diff is collapsed.
......@@ -113,10 +113,10 @@ eval.fit <- function(stepMod, data, response.var){
}
#-----------------------------------
predictions <- function(glm.fit, pred.data, selectVars){
pred.data$SoilMoist = factor(pred.data$SoilMoist)
levels(pred.data$SoilMoist) = c(-1, 0 , 1)
pred.data$SoilMoist = as.numeric(as.character(pred.data$SoilMoist))
#pred.data$SoilMoist = factor(pred.data$SoilMoist)
#levels(pred.data$SoilMoist) = c(-1, 0 , 1)
#pred.data$SoilMoist = as.numeric(as.character(pred.data$SoilMoist))
pred.data= glm.fit$dat$spdat
sc.colo = scale(glm.fit$dat$spdat[which(glm.fit$dat$spdat$trType%in%c("AA", "COLO")),selectVars])
dat.pred.colo = scale(pred.data[,selectVars], center = attr(sc.colo,"scaled:center"), scale = attr(sc.colo,"scaled:scale"))
......@@ -127,7 +127,7 @@ predictions <- function(glm.fit, pred.data, selectVars){
pred.ext = predict(glm.fit$ext,new=data.frame(dat.pred.ext),"response")
pred.lambda = 1- (pred.ext/pred.colo)
return(list(lambda = pred.lambda, sp = glm.fit$dat$sp))
return(list(lambda = pred.lambda, sp = glm.fit$dat$sp, YEAR = pred.data$YEAR, SUBSITE = pred.data$SUBSITE))
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment