Commit 1805cbef authored by Boulangeat Isabelle's avatar Boulangeat Isabelle

Add spatial distri

MAJ Rmd
parent ebf237bd
No preview for this file type
test_model<- function(MyData, loglik_fct){
require(LaplacesDemon)
set.seed(129837)
Initial.Values <- rnorm(length(MyData$parm.names), 0, 1)
res = loglik_fct(Initial.Values, MyData)$LP
return(res)
}
first.fit <- function(MyData, loglik_fct){
set.seed(666)
Initial.Values <- rnorm(length(MyData$parm.names), 0, 1)
Fit = LaplaceApproximation(Model=loglik_fct, Initial.Values, Data=MyData, Method="HAR", Iterations = 50000)
Initial.Values <- as.initial.values(Fit)
set.seed(666)
Fit.RDMH1 <- LaplacesDemon(Model=loglik_fct, Data=MyData, Initial.Values,
Covar=NULL, Iterations=1000, Status=100, Thinning=1,
Algorithm="RDMH", Specs=NULL)
return(Fit.RDMH1)
}
second.fit <- function(MyData, ffit, loglik_fct){
Initial.Values <- as.initial.values(ffit)
set.seed(666)
sfit <- LaplacesDemon(loglik_fct, Data=MyData, Initial.Values,
Covar=ffit$Covar, Iterations=30000, Status=8333, Thinning=30,
Algorithm="RDMH", Specs=NULL)
return(sfit)
}
third.fit <- function(MyData, sfit, loglik_fct){
Initial.Values <- as.initial.values(sfit)
set.seed(666)
tfit <- LaplacesDemon(loglik_fct, Data=MyData, Initial.Values,
Covar=sfit$Covar, Iterations=830000, Status=11333, Thinning=830,
Algorithm="RDMH", Specs=NULL)
return(tfit)
}
create_dataset <- function(SPECIES, cover.dom){
merge_distribution <- function(sp.distri.path, ncell.side=3, plots, plot.diag = FALSE){
require(raster)
sp.distri = raster(sp.distri.path)
# 3x3 mean filter
sp.distri2 <- focal(sp.distri, w=matrix(1/(ncell.side*ncell.side),nrow=ncell.side,ncol=ncell.side), na.rm = TRUE)
require(sp)
plots.sp = plots[,c("SemiUniquePLOT", "Latitude", "Longitude")]
coordinates(plots.sp) = ~Longitude+Latitude
projection(plots.sp) = CRS("+proj=longlat +datum=WGS84")
plots.proj = spTransform(plots.sp, projection(sp.distri2))
distri = data.frame(SemiUniquePLOT = plots.proj$SemiUniquePLOT, distriArea = extract(sp.distri2, plots.proj, method="simple"), Longitude = plots.sp$Longitude, Latitude = plots.sp$Latitude)
# nb : some plots are not in the arctic and generate a NaN value
if(plot.diag ==TRUE){
pdf("diag_distri.pdf")
plot(sp.distri2)
points(plots.proj)
dev.off()
}
return(distri)
}
#-----------
create_dataset <- function(SPECIES, cover.dom, sp.distri.path){
# Plots
# only subsites where the species have been seen (after the PA filter)
# only subsites where the species have been seen
subsites = unique(cover.dom$SUBSITE[which(cover.dom$Name==SPECIES)])
# length(subsites)
cover.sp = cover.dom[which(cover.dom$SUBSITE %in% subsites),]
sites = unique(cover.sp[,c("SemiUniquePLOT", "YEAR")])
# nrow(sites)
# table where sp considered as present (abundance filter)
sp.table = cover.sp[which(cover.sp$Name ==SPECIES),]
# length(table(sp.table$SemiUniquePLOT))
# length(table(sp.table$SITE))
# length(table(sp.table$SUBSITE))
summary(sp.table$COVER)
sp.table.by.site = split(sp.table, sp.table$SUBSITE)
# length(sp.table.by.site)
# extract plot info
# length(unique(sites$SemiUniquePLOT))
plots = unique(cover.sp[,-which(colnames(cover.sp) %in% c("YEAR", "COVER", "Name", "Genus", "Family","GFNARROWwalker", "GFNARROWarft", "UniquePLOT", "TotalCover", "RelCover", "Morphosp"))])
# length(unique(plots$SemiUniquePLOT))
# length(unique(plots$SUBSITE))
# nrow(plots)
# build transition table # need cover.sp
sp.table.tr.list = lapply(sp.table.by.site, function(x){
......@@ -50,20 +64,12 @@ sp.table.tr.list = lapply(sp.table.by.site, function(x){
return(df)
})
sp.table.tr= do.call(rbind, sp.table.tr.list)
# head(sp.table.tr)
# nrow(sp.table.tr)
# table(sp.table.tr$st0, sp.table.tr$st1)
# merge with plot info
sp.table.tr.all = merge(sp.table.tr, plots, by = "SemiUniquePLOT")
# head(sp.table.tr.all)
nrow(sp.table.tr.all)
# intervals
sp.table.tr.all$interval = sp.table.tr.all$y1 - sp.table.tr.all$y0
# hist(sp.table.tr.all$interval, breaks = 20)
# summary(sp.table.tr.all$interval)
# explo colo/ extinctions
sp.table.tr.all$trType = "present"
......@@ -71,21 +77,10 @@ sp.table.tr.all$trType[which(sp.table.tr.all$st0==0 & sp.table.tr.all$st1==1)]="
sp.table.tr.all$trType[which(sp.table.tr.all$st0==1 & sp.table.tr.all$st1==0)]="extinction"
sp.table.tr.all$trType[which(sp.table.tr.all$st0==0 & sp.table.tr.all$st1==0)]="absent"
# par(mfrow = c(1,2))
# plot(as.factor(sp.table.tr.all$trType[sp.table.tr.all$SoilMoist=="DRY"]),sp.table.tr.all$WarmQuarterTemp[sp.table.tr.all$SoilMoist=="DRY"], notch=TRUE, col = "lightblue")
# plot(as.factor(sp.table.tr.all$trType[sp.table.tr.all$SoilMoist=="MOIST"]),sp.table.tr.all$WarmQuarterTemp[sp.table.tr.all$SoilMoist=="MOIST"], notch=TRUE, col = "lightblue")
# summary(sp.table.tr.all)
# table(sp.table.tr.all$SoilMoist)
# str(sp.table.tr.all)
#--------------------------------------------------
# combine with climate (5 years earlier average)
#---------------------------------------------------
climData = unique(sp.table.tr.all[,c("SUBSITE","y0")])
# head(climData)
# dim(climData)
# str(climData)
#snowdays
climData$snowDays = apply(climData, 1, function(x){
......@@ -111,15 +106,20 @@ climData$GST = apply(climData, 1, function(x){
return(GST)
})
# ------------------------------------------------
# distribution
#------------------------------------------------
distri = merge_distribution (sp.distri.path, ncell.side=3, plots, plot.diag = TRUE)
#--------------------------------------------------
# merge datasets
#---------------------------------------------------
# str(climData)
# str(sp.table.tr.all)
dataset = merge(climData, sp.table.tr.all[, c("SemiUniquePLOT", "st0", "st1", "y0", "y1", "N", "SUBSITE", "SITE", "PLOT", "SppInPlot", "SoilMoist","interval", "trType")], by = c("SUBSITE", "y0"))
n1 = nrow(dataset)
# str(dataset)
dataset2 = merge(dataset, distri, by = "SemiUniquePLOT")
if(!(nrow(dataset2) == n1)) stop("error merging datasets")
return(dataset)
return(dataset2)
}
File added
---
title: "test vaccinium vi"
author: "Isabelle Boulangeat"
date: "6/01/2019"
output: html_document
---
# Prepare datasets
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval=TRUE)
```
## Load climate
```{r}
temp = read.table("shrubhub/ARCfunc_timeseries_temp.txt", h=T, sep="\t", encoding = "UTF-8", comment.char="", quote = "")
precip = read.table("shrubhub/ARCfunc_timeseries_prec.txt", h=T, sep="\t", encoding = "UTF-8", comment.char="", quote = "")
climato = read.table("shrubhub/TVC_SITE_SUBSITE_UPDATED2016_snow_GSL_GST_CHELSAclimatologies.txt", h=T, sep="\t", encoding = "UTF-8", comment.char="", quote = "")
climato$SiteSubsite = paste(climato$SITE, climato$SUBSITE, sep=":")
```
## Load shrub cover
```{r, fig=TRUE}
load("coverc_sub.rdata")
hist(coverc.sub$COVER)
summary(coverc.sub$COVER)
```
## Filter out treatment plots
```{r}
cover.sub = coverc.sub[which(coverc.sub$TRTMT =="CTL"),]
```
## filter out un-repeated plots
```{r, fig=TRUE}
cover.sub = cover.sub[which(cover.sub$RepeatedPlots ==1),]
hist(asin(sqrt(cover.sub$COVER/100)), breaks =21, col = "lightblue")
summary(cover.sub$COVER)
```
## remove presences when low abundance (optional)
```{r, fig=TRUE}
#----
PA_COVER_TRESHOLD = 1
#----
cover.dom = cover.sub[which(cover.sub$COVER>PA_COVER_TRESHOLD),]
```
## create species-specific dataset
```{r}
SPECIES = "Vaccinium vitis-idaea"
# SPECIES = "Betula nana"
# SPECIES = "Salix arctica"
# SPECIES = "Salix pulchra"
# SPECIES = "Vaccinium uliginosum"
# ----
source("R_fct/shape_data.r")
spdat = create_dataset(SPECIES, cover.dom, sp.distri.path= "ArcticShrubCurrent/Vaccinium.vitisidaea_cur_pc_n.img")
```
The steps of `create_dataset`:
- only subsites where the species have been seen
- extract plot info
- build transition table
- merge with plot info
- calculate interval (years)
- tag "present" and "absent"
- combine with climate (5 years earlier average)
- extract spatial distribution
- merge all
## some checks
```{r}
table(spdat$trType)
unique(spdat[which(spdat$st0=="colonisation"&spdat$distriArea==0),c("SUBSITE", "y0", "y1", "st0", "st1", "SoilMoist")])
```
# Fit
## Cleaning-Ajustments
```{r, fig=TRUE}
hist(spdat$interval)
```
```{r}
spdat = spdat[spdat$interval<10,]
spdat = spdat[-which(is.nan(spdat$distriArea)),]
selectVars = c("snowDays", "GSL", "GST", "N", "SoilMoist")
spdat$SoilMoist = factor(spdat$SoilMoist)
levels(spdat$SoilMoist) = c(1:3)
spdat$SoilMoist = as.numeric(as.character(spdat$SoilMoist))
head(spdat)
dim(spdat)
```
## Build data object for LD fit
```{r}
require(LaplacesDemon)
source("R_fct/model_fct.r")
nvars = length(selectVars)
dat = scale(spdat[,selectVars])
parm.names = as.parm.names(list(aE=rep(0, 1+nvars), aC=rep(0, 1+nvars) ))
nbetas = length(parm.names)
MyData <- list(N = nrow(dat), dat=dat, mon.names = c("logLik") , parm.names= parm.names,
var.aE = 1:(nvars),
var.aC = 1:(nvars),
aE.pos = grep("aE", parm.names),
aC.pos = grep("aC", parm.names),
AA.pos = spdat$trType=="absent",
AP.pos = spdat$trType=="colonisation",
PA.pos = spdat$trType=="extinction",
PP.pos = spdat$trType=="present",
invlogit_ = invlogit_,
itime = spdat$interval,
density = 1,
# density = spdat$distriArea,
# lik_fct = lik_dd,
nbetas = nbetas)
MyData_N = MyData
MyData_N$density = spdat$distriArea
```
## Fit
```{r, eval=FALSE}
source("R_fct/fit_LP.r")
test_model(MyData, loglik_fct)
test_model(MyData_N, loglik_fct)
#----
ffit = first.fit(MyData, loglik_fct)
Consort(ffit)
sfit = second.fit(MyData, ffit,loglik_fct)
Consort(sfit)
tfit = second.fit(MyData, sfit,loglik_fct)
Consort(tfit)
#----
ffit_N = first.fit(MyData_N, loglik_fct)
Consort(ffit_N)
sfit_N = second.fit(MyData_N, ffit_N,loglik_fct)
Consort(sfit_N)
```
This diff is collapsed.
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