diff --git a/.DS_Store b/.DS_Store index aee14fbaf2cf0228c6ff935f308fb485d8a9fcd3..6cc5ef7076185357ee44c8454812adeecf7ebbb7 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/Main_analyse.Rmd b/Main_analyse.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..7a579fd02649c378f47b75f664a85ff6b033aeda --- /dev/null +++ b/Main_analyse.Rmd @@ -0,0 +1,93 @@ +--- +title: "Main_analyse.Rmd" +author: "Isabelle Boulangeat" +date: "23/01/2020" +output: html_document +editor_options: + chunk_output_type: console +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +library(Jmisc) +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)) +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"))] +``` + +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** + +## 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 + +```{r species_specific_datasets} +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)) + +lapply(spDATA.list, function(x) table(x$spdat$trType)) +``` + + +## Fit models +```{r species_specific_models} +spDATA.fit = lapply(spDATA.list, fit.species) +``` + +## Summary of model evaluation +```{r model_evaluation} +spDATA.fit[[1]]$colo +eval.models = lapply(spDATA.fit, function(glm.fit){ + eval.colo = eval.fit(glm.fit$colo, glm.fit$dat$colo, "colo.event")[-1] + vars.colo = paste(names(sort(abs(glm.fit$colo$coefficients[-1]), dec=T)), collapse=";") + eval.ext = eval.fit(glm.fit$ext, glm.fit$dat$ext, "ext.event")[-1] + vars.ext = paste(names(sort(abs(glm.fit$ext$coefficients[-1]), dec=T)), 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)) +}) +res.tab = do.call(rbind, eval.models) +res.tab +``` + +# ```{r model_evaluation} +# plot_model(glm.fit$colo, transform = NULL, show.values=TRUE, title = paste(SPECIES, "colonisation")) +# plot_model(glm.fit$ext, transform = NULL, show.values=TRUE, title = paste(SPECIES, "extinction")) +# plot_model(glm.fit$colo, type = "eff") +# ``` + diff --git a/Main_analyse.html b/Main_analyse.html new file mode 100644 index 0000000000000000000000000000000000000000..a40cb242d66e0989d48981844c920b2bdb86a4bc --- /dev/null +++ b/Main_analyse.html @@ -0,0 +1,771 @@ + + + + + + + + + + + + + + + +Main_analyse.Rmd + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + +
+

Load data

+ +
DATA = load.datasets()
+str(DATA$cover)
+
## 'data.frame':    32793 obs. of  43 variables:
+##  $ SUBSITE             : chr  "ALEXFIORD:CASSIOPE" "ALEXFIORD:CASSIOPE" "ALEXFIORD:CASSIOPE" "ALEXFIORD:CASSIOPE" ...
+##  $ SemiUniquePLOT      : chr  "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL" "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL" "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL" "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL" ...
+##  $ YEAR                : num  1995 1995 1995 1995 1995 ...
+##  $ Name                : chr  "Cassiope tetragona" "Papaver radicatum" "Oxyria digyna" "Luzula confusa" ...
+##  $ SITE                : chr  "ALEXFIORD" "ALEXFIORD" "ALEXFIORD" "ALEXFIORD" ...
+##  $ TRTMT               : chr  "CTL" "CTL" "CTL" "CTL" ...
+##  $ PLOT                : chr  "Cas.c.c.10new" "Cas.c.c.10new" "Cas.c.c.10new" "Cas.c.c.10new" ...
+##  $ GFNARROWwalker      : chr  "SEVER" "FORB" "FORB" "RUSH" ...
+##  $ GFNARROWarft        : chr  "ESHRUB" "FORBSV" "FORBSV" "GRAMINOID" ...
+##  $ COVER               : num  22 0 1 0 6 3 7 5 0 8 ...
+##  $ UniquePLOT          : chr  "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL_1995" "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL_1995" "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL_1995" "ALEXFIORD:CASSIOPE_Cas.c.c.10new_CTL_1995" ...
+##  $ NumSurveysPerSubsite: int  3 3 3 3 3 3 3 3 3 3 ...
+##  $ NumSurveysPerPlot   : int  3 3 3 3 3 3 3 3 3 3 ...
+##  $ MinYearSubsite      : num  1995 1995 1995 1995 1995 ...
+##  $ MaxYearSubsite      : num  2007 2007 2007 2007 2007 ...
+##  $ DurationSubsite     : num  13 13 13 13 13 13 13 13 13 13 ...
+##  $ MinYearPlot         : num  1995 1995 1995 1995 1995 ...
+##  $ MaxYearPlot         : num  2007 2007 2007 2007 2007 ...
+##  $ DurationPlot        : num  13 13 13 13 13 13 13 13 13 13 ...
+##  $ SiteOnly            : chr  "ALEXFIORD" "ALEXFIORD" "ALEXFIORD" "ALEXFIORD" ...
+##  $ SubsiteOnly         : chr  "CASSIOPE" "CASSIOPE" "CASSIOPE" "CASSIOPE" ...
+##  $ MultipleSubsites    : num  1 1 1 1 1 1 1 1 1 1 ...
+##  $ RepeatedSubsites    : num  1 1 1 1 1 1 1 1 1 1 ...
+##  $ RepeatedPlots       : num  1 1 1 1 1 1 1 1 1 1 ...
+##  $ MixedPlots          : num  0 0 0 0 0 0 0 0 0 0 ...
+##  $ Morphosp            : num  0 0 0 0 0 0 0 0 1 0 ...
+##  $ Genus               : chr  "Cassiope" "Papaver" "Oxyria" "Luzula" ...
+##  $ Family              : chr  "Ericaceae" "Papaveraceae" "Polygonaceae" "Juncaceae" ...
+##  $ SppInPlot           : num  1 1 1 1 1 1 1 1 1 1 ...
+##  $ TotalCover          : num  44 44 44 44 44 44 44 44 44 46 ...
+##  $ RelCover            : num  0.5 0 0.0227 0 0.1364 ...
+##  $ Latitude            : num  78.9 78.9 78.9 78.9 78.9 ...
+##  $ Longitude           : num  -75.8 -75.8 -75.8 -75.8 -75.8 ...
+##  $ MaxTemp             : num  70.5 70.5 70.5 70.5 70.5 ...
+##  $ MAT                 : num  -164 -164 -164 -164 -164 ...
+##  $ MinTemp             : num  -367 -367 -367 -367 -367 ...
+##  $ WarmQuarterTemp     : num  25.6 25.6 25.6 25.6 25.6 ...
+##  $ ColdQuarterTemp     : num  -321 -321 -321 -321 -321 ...
+##  $ WClimGrid           : num  57684506 57684506 57684506 57684506 57684506 ...
+##  $ CRUGrid             : num  16049 16049 16049 16049 16049 ...
+##  $ WClimGroup          : num  2 2 2 2 2 2 2 2 2 2 ...
+##  $ SoilMoist           : chr  "MOIST" "MOIST" "MOIST" "MOIST" ...
+##  $ WClimGroupSM        : num  3 3 3 3 3 3 3 3 3 3 ...
+
str(DATA$climato)
+
## 'data.frame':    257 obs. of  172 variables:
+##  $ SITE                                : Factor w/ 62 levels "ABISKO","AKUREYRI",..: 3 3 3 3 3 3 3 3 3 38 ...
+##  $ SUBSITE                             : Factor w/ 254 levels "A","ABISKODRY",..: 1 37 46 64 67 68 82 122 123 215 ...
+##  $ COMMTYPE                            : Factor w/ 5 levels "","DRY","MIXED",..: 5 5 5 5 2 2 5 2 2 4 ...
+##  $ LAT                                 : num  78.9 78.9 78.9 78.9 78.9 ...
+##  $ LONG                                : num  -75.7 -75.7 -75.7 -75.7 -75.9 ...
+##  $ Checked_Coords                      : Factor w/ 2 levels "","x": 1 1 1 1 2 2 1 2 2 2 ...
+##  $ HasCover                            : int  0 0 0 0 1 1 0 1 1 1 ...
+##  $ ELEV                                : int  30 30 30 30 540 540 30 540 540 580 ...
+##  $ AZONE                               : Factor w/ 4 levels "ALPINE","ANT",..: 3 3 3 3 3 3 3 3 3 4 ...
+##  $ PI                                  : Factor w/ 42 levels "","Blok","COOPER",..: 12 12 12 12 12 12 12 22 22 24 ...
+##  $ DomGrazer                           : Factor w/ 10 levels "","LARGE","Small",..: 8 8 8 8 8 8 8 5 5 7 ...
+##  $ GrazerIntensity                     : Factor w/ 5 levels "","LOW","high",..: 4 4 4 4 4 4 4 4 4 5 ...
+##  $ CAVM                                : Factor w/ 22 levels "","B1","B3","B4",..: 21 21 21 21 4 3 21 4 3 8 ...
+##  $ CAVMBROAD                           : Factor w/ 7 levels "","B","G","P",..: 7 7 7 7 2 2 7 2 2 3 ...
+##  $ PlotSize_m2                         : Factor w/ 16 levels "","0.04","0.0531",..: 1 1 1 1 12 12 1 16 16 6 ...
+##  $ SurveyMethod                        : Factor w/ 13 levels "","BraunBlanquet",..: 4 4 4 4 8 8 4 3 3 2 ...
+##  $ HitsPerPlot                         : Factor w/ 10 levels "","100","137",..: 1 1 1 1 2 2 1 NA NA NA ...
+##  $ Comments                            : Factor w/ 20 levels "","\"0.25 m2 plots in 1992; 0.49 m2 in 2009\"",..: 1 1 1 1 1 1 1 2 2 1 ...
+##  $ CHELSA_snow_days_1979               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1980               : int  365 365 365 365 365 365 365 365 365 298 ...
+##  $ CHELSA_snow_days_1981               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1982               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1983               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1984               : int  365 365 365 365 365 365 365 365 365 321 ...
+##  $ CHELSA_snow_days_1985               : int  365 365 365 365 365 365 365 365 365 271 ...
+##  $ CHELSA_snow_days_1986               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1987               : int  365 365 365 365 365 365 365 365 365 311 ...
+##  $ CHELSA_snow_days_1988               : int  365 365 365 365 365 365 365 365 365 312 ...
+##  $ CHELSA_snow_days_1989               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1990               : int  365 365 365 365 365 365 365 365 365 340 ...
+##  $ CHELSA_snow_days_1991               : int  365 365 365 365 365 365 365 365 365 334 ...
+##  $ CHELSA_snow_days_1992               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1993               : int  365 365 365 365 355 355 365 355 355 355 ...
+##  $ CHELSA_snow_days_1994               : int  365 365 365 365 365 365 365 365 365 351 ...
+##  $ CHELSA_snow_days_1995               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ CHELSA_snow_days_1996               : int  365 365 365 365 365 365 365 365 365 292 ...
+##  $ CHELSA_snow_days_1997               : int  365 365 365 365 365 365 365 365 365 317 ...
+##  $ CHELSA_snow_days_1998               : int  365 365 365 365 346 346 365 346 346 365 ...
+##  $ CHELSA_snow_days_1999               : int  365 365 365 365 354 354 365 354 354 308 ...
+##  $ CHELSA_snow_days_2000               : int  365 365 365 365 365 365 365 365 365 311 ...
+##  $ CHELSA_snow_days_2001               : int  365 365 365 365 365 365 365 365 365 307 ...
+##  $ CHELSA_snow_days_2002               : int  365 365 365 365 365 365 365 365 365 306 ...
+##  $ CHELSA_snow_days_2003               : int  365 365 365 365 365 365 365 365 365 303 ...
+##  $ CHELSA_snow_days_2004               : int  365 365 365 365 365 365 365 365 365 299 ...
+##  $ CHELSA_snow_days_2005               : int  365 365 365 365 316 316 365 316 316 328 ...
+##  $ CHELSA_snow_days_2006               : int  365 365 365 365 340 340 365 340 340 318 ...
+##  $ CHELSA_snow_days_2007               : int  338 338 338 338 324 324 338 324 324 354 ...
+##  $ CHELSA_snow_days_2008               : int  334 334 334 334 297 297 334 297 297 302 ...
+##  $ CHELSA_snow_days_2009               : int  303 303 303 303 298 298 303 298 298 342 ...
+##  $ CHELSA_snow_days_2010               : int  365 365 365 365 344 344 365 344 344 283 ...
+##  $ CHELSA_snow_days_2011               : int  313 313 313 313 306 306 313 306 306 365 ...
+##  $ CHELSA_snow_days_2012               : int  312 312 312 312 295 295 312 295 295 337 ...
+##  $ CHELSA_snow_days_2013               : int  365 365 365 365 365 365 365 365 365 365 ...
+##  $ GSL_1979                            : int  0 0 0 0 0 0 0 0 0 0 ...
+##  $ GSL_1980                            : int  0 0 0 0 0 0 0 0 0 61 ...
+##  $ GSL_1981                            : int  0 0 0 0 0 0 0 0 0 15 ...
+##  $ GSL_1982                            : int  0 0 0 0 0 0 0 0 0 0 ...
+##  $ GSL_1983                            : int  0 0 0 0 0 0 0 0 0 0 ...
+##  $ GSL_1984                            : int  0 0 0 0 0 0 0 0 0 50 ...
+##  $ GSL_1985                            : int  0 0 0 0 0 0 0 0 0 53 ...
+##  $ GSL_1986                            : int  0 0 0 0 0 0 0 0 0 0 ...
+##  $ GSL_1987                            : int  0 0 0 0 0 0 0 0 0 54 ...
+##  $ GSL_1988                            : int  0 0 0 0 0 0 0 0 0 50 ...
+##  $ GSL_1989                            : int  0 0 0 0 0 0 0 0 0 18 ...
+##  $ GSL_1990                            : int  0 0 0 0 0 0 0 0 0 19 ...
+##  $ GSL_1991                            : int  0 0 0 0 0 0 0 0 0 37 ...
+##  $ GSL_1992                            : int  0 0 0 0 0 0 0 0 0 0 ...
+##  $ GSL_1993                            : int  0 0 0 0 5 5 0 5 5 13 ...
+##  $ GSL_1994                            : int  0 0 0 0 0 0 0 0 0 0 ...
+##  $ GSL_1995                            : int  0 0 0 0 0 0 0 0 0 12 ...
+##  $ GSL_1996                            : int  0 0 0 0 0 0 0 0 0 70 ...
+##  $ GSL_1997                            : int  0 0 0 0 0 0 0 0 0 43 ...
+##  $ GSL_1998                            : int  0 0 0 0 0 0 0 0 0 26 ...
+##  $ GSL_1999                            : int  0 0 0 0 10 10 0 10 10 57 ...
+##  $ GSL_2000                            : int  0 0 0 0 0 0 0 0 0 57 ...
+##  $ GSL_2001                            : int  0 0 0 0 0 0 0 0 0 68 ...
+##  $ GSL_2002                            : int  0 0 0 0 0 0 0 0 0 61 ...
+##  $ GSL_2003                            : int  0 0 0 0 0 0 0 0 0 67 ...
+##  $ GSL_2004                            : int  0 0 0 0 0 0 0 0 0 55 ...
+##  $ GSL_2005                            : int  0 0 0 0 17 17 0 17 17 35 ...
+##  $ GSL_2006                            : int  0 0 0 0 17 17 0 17 17 58 ...
+##  $ GSL_2007                            : int  3 3 3 3 23 23 3 23 23 68 ...
+##  $ GSL_2008                            : int  26 26 26 26 43 43 26 43 43 60 ...
+##  $ GSL_2009                            : int  32 32 32 32 41 41 32 41 41 24 ...
+##  $ GSL_2010                            : int  0 0 0 0 10 10 0 10 10 85 ...
+##  $ GSL_2011                            : int  17 17 17 17 34 34 17 34 34 0 ...
+##  $ GSL_2012                            : int  19 19 19 19 44 44 19 44 44 0 ...
+##  $ GSL_2013                            : int  0 0 0 0 0 0 0 0 0 8 ...
+##  $ GST_1979                            : num  NA NA NA NA NA NA NA NA NA NA ...
+##  $ GST_1980                            : num  NA NA NA NA NA ...
+##  $ GST_1981                            : num  NA NA NA NA NA ...
+##  $ GST_1982                            : num  NA NA NA NA NA NA NA NA NA NA ...
+##  $ GST_1983                            : num  NA NA NA NA NA NA NA NA NA NA ...
+##  $ GST_1984                            : num  NA NA NA NA NA ...
+##  $ GST_1985                            : num  NA NA NA NA NA ...
+##  $ GST_1986                            : num  NA NA NA NA NA NA NA NA NA NA ...
+##  $ GST_1987                            : num  NA NA NA NA NA ...
+##  $ GST_1988                            : num  NA NA NA NA NA ...
+##  $ GST_1989                            : num  NA NA NA NA NA ...
+##   [list output truncated]
+
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"))]
+

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

+
+
+

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)
+} )
+
## [1] "Vaccinium vitis-idaea"
+## [1] "Carex bigelowii"
+## [1] "Ledum palustre"
+## [1] "Betula nana"
+## [1] "Eriophorum vaginatum"
+## [1] "Salix pulchra"
+## [1] "Eriophorum chamissonis"
+## [1] "Persicaria bistorta"
+## [1] "Dupontia fisheri"
+## [1] "Rubus chamaemorus"
+
names(spDATA.list) = selection
+(spDATA.check = lapply(spDATA.list, check.nas))
+
## [1] "Vaccinium vitis-idaea"
+## [1] "Carex bigelowii"
+## [1] "Ledum palustre"
+## [1] "Betula nana"
+## [1] "Eriophorum vaginatum"
+## [1] "Salix pulchra"
+## [1] "Eriophorum chamissonis"
+## [1] "Persicaria bistorta"
+## [1] "Dupontia fisheri"
+## [1] "Rubus chamaemorus"
+
## $`Vaccinium vitis-idaea`
+## [1] 0
+## 
+## $`Carex bigelowii`
+## [1] 0
+## 
+## $`Ledum palustre`
+## [1] 0
+## 
+## $`Betula nana`
+## [1] 0
+## 
+## $`Eriophorum vaginatum`
+## [1] 0
+## 
+## $`Salix pulchra`
+## [1] 0
+## 
+## $`Eriophorum chamissonis`
+## [1] 0
+## 
+## $`Persicaria bistorta`
+## [1] 0
+## 
+## $`Dupontia fisheri`
+## [1] 0
+## 
+## $`Rubus chamaemorus`
+## [1] 0
+
lapply(spDATA.list, function(x) table(x$spdat$trType))
+
## $`Vaccinium vitis-idaea`
+## 
+##   AA COLO  EXT   PP 
+##  168   18   13  647 
+## 
+## $`Carex bigelowii`
+## 
+##   AA COLO  EXT   PP 
+##  194   53   38  490 
+## 
+## $`Ledum palustre`
+## 
+##   AA COLO  EXT   PP 
+##  147   22   17  513 
+## 
+## $`Betula nana`
+## 
+##   AA COLO  EXT   PP 
+##  191   19    5  510 
+## 
+## $`Eriophorum vaginatum`
+## 
+##   AA COLO  EXT   PP 
+##  199   34   27  419 
+## 
+## $`Salix pulchra`
+## 
+##   AA COLO  EXT   PP 
+##  287   50   31  376 
+## 
+## $`Eriophorum chamissonis`
+## 
+##   AA COLO  EXT   PP 
+##   81   40   46  251 
+## 
+## $`Persicaria bistorta`
+## 
+##   AA COLO  EXT   PP 
+##  296   77   60  183 
+## 
+## $`Dupontia fisheri`
+## 
+##   AA COLO  EXT   PP 
+##  214   29   38  218 
+## 
+## $`Rubus chamaemorus`
+## 
+##   AA COLO  EXT   PP 
+##  315   34   22  191
+
+
+

Fit models

+
spDATA.fit = lapply(spDATA.list, fit.species)
+
## [1] "Vaccinium vitis-idaea"
+## [1] "Carex bigelowii"
+
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
## [1] "Ledum palustre"
+## [1] "Betula nana"
+
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
## [1] "Eriophorum vaginatum"
+## [1] "Salix pulchra"
+
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
## [1] "Eriophorum chamissonis"
+## [1] "Persicaria bistorta"
+
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
+## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+
## [1] "Dupontia fisheri"
+## [1] "Rubus chamaemorus"
+
+
+

Summary of model evaluation

+
spDATA.fit[[1]]$colo
+
## 
+## Call:  glm(formula = colo.event ~ GSL + GST + TotalCover + SoilMoist + 
+##     GSL:SoilMoist + GST:SoilMoist + TotalCover:SoilMoist, family = "binomial", 
+##     data = spDATA$colo)
+## 
+## Coefficients:
+##          (Intercept)                   GSL                   GST  
+##              -3.9995               -0.9393                3.2678  
+##           TotalCover             SoilMoist         GSL:SoilMoist  
+##               0.9502               -1.4958               -2.5655  
+##        GST:SoilMoist  TotalCover:SoilMoist  
+##               6.5223                3.2555  
+## 
+## Degrees of Freedom: 185 Total (i.e. Null);  178 Residual
+## Null Deviance:       118.3 
+## Residual Deviance: 77.27     AIC: 93.27
+
eval.models = lapply(spDATA.fit, function(glm.fit){
+  eval.colo = eval.fit(glm.fit$colo, glm.fit$dat$colo, "colo.event")[-1]
+  vars.colo = paste(names(sort(abs(glm.fit$colo$coefficients[-1]), dec=T)), collapse=";")
+  eval.ext = eval.fit(glm.fit$ext, glm.fit$dat$ext, "ext.event")[-1]
+  vars.ext = paste(names(sort(abs(glm.fit$ext$coefficients[-1]), dec=T)), 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))
+})
+res.tab = do.call(rbind, eval.models)
+res.tab
+
##                        colo.R2      colo.AUC  ext.R2       ext.AUC  
+## Vaccinium vitis-idaea  0.420457     0.8941799 0.03035965   0.6103912
+## Carex bigelowii        0.3131377    0.7919665 0.08320889   0.6446294
+## Ledum palustre         0.1403581    0.6742424 0.03557624   0.6017659
+## Betula nana            0.2892933    0.8483053 0.1751081    0.8796078
+## Eriophorum vaginatum   0.02042669   0.5316287 0.0713155    0.6987978
+## Salix pulchra          0.1812621    0.7474913 0.02601676   0.5888384
+## Eriophorum chamissonis 0.2176915    0.7401235 0.064756     0.6487095
+## Persicaria bistorta    0.1146254    0.6679756 0.08382495   0.6345628
+## Dupontia fisheri       0.2435445    0.7993877 0.2461762    0.8039594
+## Rubus chamaemorus      2.351904e-16 0.5       4.573867e-16 0.5      
+##                        vars.colo                                                                                
+## Vaccinium vitis-idaea  "GST:SoilMoist;GST;TotalCover:SoilMoist;GSL:SoilMoist;SoilMoist;TotalCover;GSL"          
+## Carex bigelowii        "GSL;snowDays:SoilMoist;snowDays;TotalCover:SoilMoist;SoilMoist;TotalCover"              
+## Ledum palustre         "SoilMoist;TotalCover:SoilMoist;GST;TotalCover"                                          
+## Betula nana            "TotalCover:SoilMoist;snowDays;GSL;SoilMoist;TotalCover"                                 
+## Eriophorum vaginatum   "SoilMoist"                                                                              
+## Salix pulchra          "GSL;snowDays;SoilMoist;TotalCover:SoilMoist;TotalCover"                                 
+## Eriophorum chamissonis "snowDays:SoilMoist;GSL:SoilMoist;TotalCover:SoilMoist;SoilMoist;TotalCover;snowDays;GSL"
+## Persicaria bistorta    "SoilMoist;snowDays:SoilMoist;GST:SoilMoist;snowDays;GSL;GST"                            
+## Dupontia fisheri       "GSL:SoilMoist;snowDays;snowDays:SoilMoist;GST;GSL;SoilMoist;TotalCover"                 
+## Rubus chamaemorus      ""                                                                                       
+##                        vars.ext                                                                                 
+## Vaccinium vitis-idaea  "SoilMoist"                                                                              
+## Carex bigelowii        "SoilMoist;TotalCover:SoilMoist;snowDays;TotalCover"                                     
+## Ledum palustre         "snowDays"                                                                               
+## Betula nana            "TotalCover;snowDays"                                                                    
+## Eriophorum vaginatum   "GSL;snowDays;TotalCover;GST"                                                            
+## Salix pulchra          "TotalCover:SoilMoist;SoilMoist;TotalCover"                                              
+## Eriophorum chamissonis "snowDays:SoilMoist;GSL:SoilMoist;GSL;GST;snowDays;SoilMoist"                            
+## Persicaria bistorta    "SoilMoist;snowDays;GSL"                                                                 
+## Dupontia fisheri       "snowDays;snowDays:SoilMoist;GSL;GSL:SoilMoist;TotalCover:SoilMoist;TotalCover;SoilMoist"
+## Rubus chamaemorus      ""
+
+
+

{r model_evaluation} # plot_model(glm.fit$colo, transform = NULL, show.values=TRUE, title = paste(SPECIES, "colonisation")) # plot_model(glm.fit$ext, transform = NULL, show.values=TRUE, title = paste(SPECIES, "extinction")) # plot_model(glm.fit$colo, type = "eff") #

+
+ + + + +
+ + + + + + + + + + + + + + + diff --git a/R_fct/analyse_fct.r b/R_fct/analyse_fct.r new file mode 100644 index 0000000000000000000000000000000000000000..11d13b7028cd1b5b92607a5213821e603c1ef890 --- /dev/null +++ b/R_fct/analyse_fct.r @@ -0,0 +1,118 @@ +require(dplyr) +require(MASS) +require(fmsb) +require(ROCR) +require(LaplacesDemon) +#----------------------------------- + +load.datasets <- function(){ + 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("coverc_sub.rdata") + + cover.sub = coverc.sub[which(coverc.sub$TRTMT =="CTL"),] + cover.sub = cover.sub[which(cover.sub$RepeatedPlots ==1),] + + cover.dom = cover.sub[which(cover.sub$YEAR<2014),] + + return(list(cover=cover.dom, climato=climato)) +} + +#----------------------------------- + +dat.species <- function(SPECIES, DATA, Bayesian=TRUE){ + spdat = create_dataset(SPECIES, DATA$cover, DATA$climato) + spdat$GSL_anom = spdat$GSL-spdat$GSL_av + spdat$GST_anom = spdat$GST-spdat$GST_av + spdat$snowDays_anom = spdat$snowDays-spdat$snow_days_av + spdat$trType = factor(spdat$trType) + levels(spdat$trType) <- c("AA", "COLO", "EXT", "PP") + spdat = spdat[spdat$interval<10,] + + selectVars = c("snowDays", "GSL", "GST", "SoilMoist", "TotalCover") + spdat$SoilMoist = factor(spdat$SoilMoist) + levels(spdat$SoilMoist) = c(-1,0,1) + spdat$SoilMoist = as.numeric(as.character(spdat$SoilMoist)) + + dat.ext = scale(spdat[which(spdat$trType%in%c("PP", "EXT")),selectVars]) + dat.ext[,"SoilMoist"] = spdat[which(spdat$trType%in%c("PP", "EXT")),"SoilMoist"] + + dat.colo = scale(spdat[which(spdat$trType%in%c("AA", "COLO")),selectVars]) + dat.colo[,"SoilMoist"] = spdat[which(spdat$trType%in%c("AA", "COLO")),"SoilMoist"] + + if(Bayesian){ + nvars = length(selectVars) + parm.names.ext = as.parm.names(list(aE=rep(0, 1+nvars))) + nbetas.ext = length(parm.names.ext) + MyData.ext <- list(N = nrow(dat.ext), dat=dat.ext, mon.names = c("logLik") , parm.names= parm.names.ext, + var.aE = 1:(nvars), + aE.pos = grep("aE", parm.names.ext), + PA.pos = spdat$trType[which(spdat$trType%in%c("PP", "EXT"))]=="EXT", + PP.pos = spdat$trType[which(spdat$trType%in%c("PP", "EXT"))]=="PP", + invlogit_ = invlogit_, + itime = spdat$interval[which(spdat$trType%in%c("PP", "EXT"))], + nbetas = nbetas.ext) + + parm.names.colo = as.parm.names(list(aC=rep(0, 1+nvars))) + nbetas.colo = length(parm.names.colo) + MyData.colo<- list(N = nrow(dat.colo), dat=dat.colo, mon.names = c("logLik") , parm.names= parm.names.colo, + var.aC = 1:(nvars), + aC.pos = grep("aC", parm.names.colo), + AP.pos = spdat$trType[which(spdat$trType%in%c("AA", "COLO"))]=="COLO", + AA.pos = spdat$trType[which(spdat$trType%in%c("AA", "COLO"))]=="AA", + invlogit_ = invlogit_, + density = rep(1, nrow(dat.colo)), + itime = spdat$interval[which(spdat$trType%in%c("AA", "COLO"))], + nbetas = nbetas.colo) + + LD = list(colo = MyData.colo, ext = MyData.ext) + } else (LD = NULL) + + dat.colo = data.frame(dat.colo) + dat.colo$colo.event = as.numeric(MyData.colo$AP.pos) + + dat.ext = data.frame(dat.ext) + dat.ext$ext.event = as.numeric(MyData.ext$PA.pos) + + return(list(spdat = spdat,LD = LD, colo = dat.colo, ext = dat.ext, sp = SPECIES)) +} + +#----------------------------------- +check.nas <- function(spDATA){ + print(spDATA$sp) + return(sum(is.na(spDATA$spdat$GST))) +} +#----------------------------------- + +fit.species <- function(spDATA){ + print(spDATA$sp) + + colo.mod = glm(colo.event~ (snowDays + GSL + GST + TotalCover) * SoilMoist, family = "binomial", data = spDATA$colo) + stepMod.colo = stepAIC(colo.mod, direction = "both", trace=0) + + ext.mod = glm(ext.event~ (snowDays + GSL + GST + TotalCover) * SoilMoist, family = "binomial", data = spDATA$ext) + stepMod.ext = stepAIC(ext.mod, direction = "both", trace=0) + + return(list(colo = stepMod.colo, ext = stepMod.ext, dat = spDATA)) +} + +#----------------------------------- + +eval.fit <- function(stepMod, data, response.var){ + pred = predict(stepMod,new=data,"response") + R2 = NagelkerkeR2(stepMod)$R2 + perf = performance(prediction(pred, data[,response.var]), "auc") + AUC = perf@y.values[[1]] + return(list(pred = pred, R2=R2, AUC=AUC)) +} + + +#----------------------------------- + +fit.species.LD <- function(spDATA){ + +} diff --git a/R_fct/fit_LP.r b/R_fct/fit_LP.r index 5f082c4bb015b6f3ffceb90b0f84923ce80733fb..99014324f89893bd8903a21e0a1199e0b5181bea 100644 --- a/R_fct/fit_LP.r +++ b/R_fct/fit_LP.r @@ -1,3 +1,4 @@ + test_model<- function(MyData, loglik_fct){ require(LaplacesDemon) set.seed(129837) diff --git a/R_fct/graphs.r b/R_fct/graphs.r new file mode 100644 index 0000000000000000000000000000000000000000..b0d93212030e6942622a955b7854505a0f99ff99 --- /dev/null +++ b/R_fct/graphs.r @@ -0,0 +1,41 @@ +require(ggplot2) +require(reshape2) +require(sjPlot) +plot.effects.ld <- function(Fit.mcmc, vars, nmetapars=1, colos = c("lightblue", "steelblue", "grey50", "tomato", "orange")){ + all.vars = rep(vars,nmetapars) + df = as.data.frame(Fit.mcmc$Posterior1) + df.long= melt(df) + df.long$par = unlist(lapply(df.long$variable, substr, start=1, stop=2)) + df.long$index = unlist(lapply(df.long$variable, gsub, pattern="[^(0-9)]", replacement="")) + df.long2 = df.long[df.long$index>1,] + df.long2$var = factor(df.long2$index) + levels(df.long2$var) = all.vars + levels(df.long2$var) = vars + df.long2$par = factor(df.long2$par) + levels(df.long2$par) + # = c("aC", "aD", "bC", "bD", "d") + p1 = ggplot(df.long2, aes(x = var, y = value, fill = var)) + + geom_violin(trim= FALSE) + + geom_boxplot(alpha=0.3, width=0.2, outlier.size = 0) + + facet_wrap(~ par, ncol=3, dir="v", scales="free") + + theme(axis.title.x=element_blank(), + axis.text.x=element_blank(), + axis.ticks.x=element_blank()) + p2 = p1 + geom_hline(yintercept=0, linetype="dashed", color = "black") + + labs(title="Posterior distributions", x="variable", y="effect") + + scale_fill_manual(name="variables", + labels=levels(df.long2$var) , + values = colos) + + return(p2) +} + +#------------------------- +plot.effects.glm <- function(SPECIES, glm.fit, colos = c("lightblue", "steelblue", "grey50", "tomato", "orange")){ + coeff.colo = summary(glm.fit$colo)$coefficients[-1,1] + vars.colo = rownames(glm.fit$colo$coeff)[-1] + plot_model(glm.fit$colo, transform = NULL, show.values=TRUE, title = SPECIES) + plot_model(glm.fit$colo, transform = NULL, show.values=TRUE, title = SPECIES) + + return(barplot(coeff.colo, col = colos)) +} diff --git a/R_fct/model_fct.r b/R_fct/model_fct.r index fa3c3cecb007173d1ed91080375e91d8f9bb0b73..a657ce4c9e307f6f42cae5a9132595b10ebb45fd 100644 --- a/R_fct/model_fct.r +++ b/R_fct/model_fct.r @@ -51,3 +51,71 @@ invlogit_ <- function(x, itime) proba_itime = 1 - (1 - proba)^(itime) return(proba_itime) } + + +ext_model <- function(parm, Data) +{ + require(data.table) + aE = parm[Data$aE.pos] + # fit transitions + logit_alphaE = tcrossprod(aE, as.matrix(cbind(rep(1, Data$N), Data$dat[,Data$var.aE]))) + + # compute transitions accounting for interval time and logit transformation + # ! might be NaN because exp(bigNumber) + aE.fct=as.vector(Data$invlogit_(logit_alphaE, Data$itime)) + + # LIKELIHOOD + lik = rep(NA, nrow(Data$dat)) + lik[which(Data$PA.pos)]= (aE.fct) [Data$PA.pos] + lik[which(Data$PP.pos)]= (1-aE.fct) [Data$PP.pos] + + lik[lik==0] = .Machine$double.xmin + loglik = sum(log(lik)) + + + # PRIOR + # prior cauchy following Gelman 2008 + # intercepts different of slope + beta.aE = dcauchy(aE[1], 1, 10, log = TRUE) + sum(dcauchy(aE[-1], 0, 2.5, log = TRUE)) + + # POSTERIOR + LP = loglik + beta.aE + + return(list(LP = LP, Dev = -2*loglik, Monitor = LP, yhat = lik,parm = parm)) +} + +colo_model <- function(parm, Data) +{ + require(data.table) + + aC = parm[Data$aC.pos] + + # fit transitions + logit_alphaC = tcrossprod(aC, as.matrix(cbind(rep(1, Data$N), Data$dat[,Data$var.aC]))) + + # compute transitions accounting for interval time and logit transformation + # ! might be NaN because exp(bigNumber) + aC.fct=as.vector(Data$invlogit_(logit_alphaC, Data$itime)) + + N = Data$density + + # LIKELIHOOD + lik = rep(NA, nrow(Data$dat)) + + lik[which(Data$AA.pos)]= (1-(aC.fct*N)) [Data$AA.pos] + lik[which(Data$AP.pos)]= (N*aC.fct) [Data$AP.pos] + + lik[lik==0] = .Machine$double.xmin + loglik = sum(log(lik)) + + + # PRIOR + # prior cauchy following Gelman 2008 + # intercepts different of slope + beta.aC = dcauchy(aC[1], 1, 10, log = TRUE) + sum(dcauchy(aC[-1], 0, 2.5, log = TRUE)) + + # POSTERIOR + LP = loglik + beta.aC + + return(list(LP = LP, Dev = -2*loglik, Monitor = LP, yhat = lik,parm = parm)) +} diff --git a/R_fct/shape_data.r b/R_fct/shape_data.r index 4bab8d9978ab0f75d6e9f614e478cdcd944bbaa6..78fb51af13b5c9a754bd5da6ab115d35335c1543 100644 --- a/R_fct/shape_data.r +++ b/R_fct/shape_data.r @@ -1,72 +1,71 @@ - -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) -} +# +# 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 = unique(plots[,c("SUBSITE", "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(SUBSITE = plots.proj$SUBSITE, 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){ +create_dataset <- function(SPECIES, cover.dom, climato){ # Plots # only subsites where the species have been seen subsites = unique(cover.dom$SUBSITE[which(cover.dom$Name==SPECIES)]) cover.sp = cover.dom[which(cover.dom$SUBSITE %in% subsites),] -sites = unique(cover.sp[,c("SemiUniquePLOT", "YEAR")]) +sites = unique(cover.sp[,c("SemiUniquePLOT", "YEAR", "SUBSITE")]) # table where sp considered as present (abundance filter) -sp.table = cover.sp[which(cover.sp$Name ==SPECIES),] +sp.table = cover.sp[which(cover.sp$Name ==SPECIES),c("SemiUniquePLOT", "YEAR", "COVER")] +sp.table$PA = ifelse(sp.table$COVER>0, 1, 0) -sp.table.by.site = split(sp.table, sp.table$SUBSITE) +# long table with presences and absences +long.tab = merge(sites, sp.table, by = c("SemiUniquePLOT", "YEAR"), all.x=TRUE) +long.tab$PA[is.na(long.tab$PA)] = 0 -# 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"))]) +#sp.table.by.site = split(cover.sp, cover.sp$SemiUniquePLOT) # build transition table # need cover.sp -sp.table.tr.list = lapply(sp.table.by.site, function(x){ -# x = sp.table.by.site[[1]] - # print(x$SemiUniquePLOT) - p = cover.sp[which(cover.sp$SUBSITE%in%x$SUBSITE),] # all plots - cover = merge(p[which(p$Name==SPECIES),c("SemiUniquePLOT", "YEAR","COVER")], - unique(p[,c("SemiUniquePLOT", "YEAR")]), all=TRUE) - cover[which(is.na(cover$COVER)), "COVER"]=0 +sp.table.tr.list = lapply(unique(long.tab$SemiUniquePLOT), function(x){ +# x = unique(long.tab$SemiUniquePLOT)[10] + p = long.tab[which(long.tab$SemiUniquePLOT==x),] + # p = cover.sp[which(cover.sp$SemiUniquePLOT%in%x$SemiUniquePLOT),] # all years + # cover = merge(p[which(p$Name==SPECIES),c("SemiUniquePLOT", "YEAR","COVER")], + # unique(p[,c("SemiUniquePLOT", "YEAR")]), all=TRUE) + # #cover[which(is.na(cover$COVER)), "COVER"]=0 ntransitions = length(unique(p$YEAR))-1 - nplots = length(unique(p$SemiUniquePLOT)) df = data.frame(SemiUniquePLOT = rep(unique(p$SemiUniquePLOT), each =ntransitions), - st0 = rep(0,ntransitions*nplots), - st1 = rep(0,ntransitions*nplots), - y0 =rep(sort(unique(p$YEAR))[-length(unique(p$YEAR))],length(unique(p$PLOT))), - y1 =rep(sort(unique(p$YEAR))[-1],length(unique(p$PLOT))), - N = rep(by(cover[,"COVER"], cover[,"YEAR"], mean)[1:ntransitions], nplots) + y0 =sort(unique(p$YEAR))[-length(unique(p$YEAR))], + y1 =sort(unique(p$YEAR))[-1] ) - df$st0[which(paste(df$SemiUniquePLOT, df$y0, sep="_")%in%x$UniquePLOT)] = 1 - df$st1[which(paste(df$SemiUniquePLOT, df$y1, sep="_")%in%x$UniquePLOT)] = 1 + df$st0 =merge(df, p[, c("YEAR", "PA")], by.x = "y0", by.y = "YEAR")$PA + df$st1 =merge(df, p[, c("YEAR", "PA")], by.x = "y1", by.y = "YEAR")$PA + df return(df) }) sp.table.tr= do.call(rbind, sp.table.tr.list) + # merge with plot info -sp.table.tr.all = merge(sp.table.tr, plots, by = "SemiUniquePLOT") +sp.table.tr.all = unique(merge(sp.table.tr, cover.sp[, c("SemiUniquePLOT", "YEAR", "SoilMoist", "TotalCover")], by.x = c("SemiUniquePLOT", "y0"),by.y= c("SemiUniquePLOT", "YEAR"))) # intervals sp.table.tr.all$interval = sp.table.tr.all$y1 - sp.table.tr.all$y0 @@ -77,49 +76,100 @@ 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" + +# extract plot info + #-------------------------------------------------- # combine with climate (5 years earlier average) #--------------------------------------------------- -climData = unique(sp.table.tr.all[,c("SUBSITE","y0")]) -#snowdays -climData$snowDays = apply(climData, 1, function(x){ - # print(as.numeric(x[2])) - subClim = climato[which(climato$SiteSubsite==x[1]),] - snowDays = mean(as.numeric(subClim[sapply(1:5, function(k){grep(paste("snow_days_",as.numeric(x[2])-k+1,sep=""), names(subClim), value=FALSE)})])) +climData = unique(sites[,c("SUBSITE","YEAR")]) + +grepClim = function(climData, climato, variable, nyears){ + clim = apply(climData, 1, function(x){ + subClim = climato[which(climato$SiteSubsite==x[[1]]),] + meanClim = mean(as.numeric(subClim[sapply(1:nyears, function(k){grep(paste(variable,as.numeric(x[[2]])-k+1,sep=""), names(subClim), value=FALSE)})])) + return(meanClim) + }) + return(clim) +} + +climData$snowDays = grepClim(climData, climato, "snow_days_", 5) +climData$GSL = grepClim(climData, climato, "GSL_", 5) +climData$GST = grepClim(climData, climato, "GST_", 5) + +# ------------------------------------------------ +# average +#------------------------------------------------ + +mean.snowDays = sapply(climato$SiteSubsite, function(x){ + subClim = climato[which(climato$SiteSubsite==x),] + snowDays = mean(as.numeric(subClim[sapply(1979:2013, function(k){grep(paste("snow_days_",k,sep=""), names(subClim), value=FALSE)})])) return(snowDays) }) - -#growing season length -climData$GSL = apply(climData, 1, function(x){ - # print(as.numeric(x[2])) - subClim = climato[which(climato$SiteSubsite==x[1]),] - GSL = mean(as.numeric(subClim[sapply(1:5, function(k){grep(paste("GSL_",as.numeric(x[2])-k+1,sep=""), names(subClim), value=FALSE)})])) +mean.GSL = sapply(climato$SiteSubsite, function(x){ + subClim = climato[which(climato$SiteSubsite==x),] + GSL = mean(as.numeric(subClim[sapply(1979:2013, function(k){grep(paste("GSL_",k,sep=""), names(subClim), value=FALSE)})])) return(GSL) }) - -#growing season temperature -climData$GST = apply(climData, 1, function(x){ - # print(as.numeric(x[2])) - subClim = climato[which(climato$SiteSubsite==x[1]),] - GST = mean(as.numeric(subClim[sapply(1:5, function(k){grep(paste("GST_",as.numeric(x[2])-k+1,sep=""), names(subClim), value=FALSE)})])) +mean.GST = sapply(climato$SiteSubsite, function(x){ + subClim = climato[which(climato$SiteSubsite==x),] + GST = mean(as.numeric(subClim[sapply(1979:2013, function(k){grep(paste("GST_",k,sep=""), names(subClim), value=FALSE)})])) return(GST) }) +mean.snowDays = unique(data.frame(snow_days_av=mean.snowDays, SUBSITE = names(mean.snowDays))) +mean.GSL = unique(data.frame(GSL_av=mean.GSL, SUBSITE = names(mean.GSL))) +mean.GST = unique(data.frame(GST_av=mean.GST, SUBSITE = names(mean.GST))) +require(dplyr) +all_means <- mean.snowDays %>% inner_join(mean.GSL, by = "SUBSITE") %>% inner_join(mean.GST, by ="SUBSITE") + # ------------------------------------------------ # distribution #------------------------------------------------ -distri = merge_distribution (sp.distri.path, ncell.side=3, plots, plot.diag = TRUE) +#distri = merge_distribution (sp.distri.path, ncell.side=3, plots, plot.diag = TRUE) #-------------------------------------------------- # merge datasets #--------------------------------------------------- +dataset = merge(climData, sites, by=c("SUBSITE", "YEAR")) +dataset2 = merge(dataset, sp.table.tr.all,by.x = c("SemiUniquePLOT", "YEAR"), by.y = c("SemiUniquePLOT", "y0"), all.y=TRUE) +#nrow(dataset2) +#nrow(dataset) -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) +#dataset2 = merge(dataset, distri, by = "SUBSITE", all.x=TRUE, all.y=FALSE) +#if(!(nrow(dataset2) == n1)) stop("error merging datasets") -dataset2 = merge(dataset, distri, by = "SemiUniquePLOT") -if(!(nrow(dataset2) == n1)) stop("error merging datasets") +dataset3 = merge(dataset2, all_means, by = "SUBSITE", all.x=TRUE, all.y=FALSE) +#if(!(nrow(dataset3) == n1)) stop("error merging datasets") + +return(dataset3) +} -return(dataset2) +##------------- +proj.dat.shape <- function(sites, year=2012, climato){ + proj.dat = sites + #snowdays + proj.dat$snowDays = apply(proj.dat , 1, function(x){ + # print(as.numeric(x[2])) + subClim = climato[which(climato$SiteSubsite==x[1]),] + snowDays = mean(as.numeric(subClim[sapply(1:5, function(k){grep(paste("snow_days_",as.numeric(year)-k+1,sep=""), names(subClim), value=FALSE)})])) + return(snowDays) + }) + proj.dat$GST = apply(proj.dat , 1, function(x){ + subClim = climato[which(climato$SiteSubsite==x[1]),] + GST = mean(as.numeric(subClim[sapply(1:5, function(k){grep(paste("GST_",as.numeric(year)-k+1,sep=""), names(subClim), value=FALSE)})])) + return(GST) + }) + proj.dat$GSL = apply(proj.dat , 1, function(x){ + subClim = climato[which(climato$SiteSubsite==x[1]),] + GSL = mean(as.numeric(subClim[sapply(1:5, function(k){grep(paste("GSL_",as.numeric(year)-k+1,sep=""), names(subClim), value=FALSE)})])) + return(GSL) + }) + proj.dat$GST = apply(proj.dat , 1, function(x){ + subClim = climato[which(climato$SiteSubsite==x[1]),] + GST = mean(as.numeric(subClim[sapply(1:5, function(k){grep(paste("GST_",as.numeric(year)-k+1,sep=""), names(subClim), value=FALSE)})])) + return(GST) + }) + return(proj.dat) } diff --git a/diag_distri.pdf b/diag_distri.pdf deleted file mode 100644 index 554a3a44b842fefc7ca45d3aa0b2490e3561b591..0000000000000000000000000000000000000000 Binary files a/diag_distri.pdf and /dev/null differ diff --git a/outputs/sfit.rds b/outputs/sfit.rds new file mode 100644 index 0000000000000000000000000000000000000000..2d6fa19c26c27861bd6d2f572334b9ed7c61ca5e Binary files /dev/null and b/outputs/sfit.rds differ diff --git a/outputs/sfit_BN.rds b/outputs/sfit_BN.rds new file mode 100644 index 0000000000000000000000000000000000000000..0e61300f34d943464c137a9ff1255323398fb4c0 Binary files /dev/null and b/outputs/sfit_BN.rds differ diff --git a/outputs/sfit_VVI.rds b/outputs/sfit_VVI.rds new file mode 100644 index 0000000000000000000000000000000000000000..97f157187803b758d16267f81955f4bca8fd12d9 Binary files /dev/null and b/outputs/sfit_VVI.rds differ diff --git a/outputs/tfit-.rds b/outputs/tfit-.rds new file mode 100644 index 0000000000000000000000000000000000000000..751123866f204db4f8e4800483e7232d46a30e7d Binary files /dev/null and b/outputs/tfit-.rds differ diff --git a/outputs/tfit.rds b/outputs/tfit.rds new file mode 100644 index 0000000000000000000000000000000000000000..234e2a73c1dc8d3f57790aff9af4fe06de7a6598 Binary files /dev/null and b/outputs/tfit.rds differ diff --git a/outputs/tfit_BN.RData b/outputs/tfit_BN.RData new file mode 100644 index 0000000000000000000000000000000000000000..ea1ec4b932f980e9e386eff3166885effea1f0fd Binary files /dev/null and b/outputs/tfit_BN.RData differ diff --git a/outputs/tfit_BN.rds b/outputs/tfit_BN.rds new file mode 100644 index 0000000000000000000000000000000000000000..5b78ff92a5afca7c4c82bd32cbdbe7ffbaa34697 Binary files /dev/null and b/outputs/tfit_BN.rds differ diff --git a/outputs/tfit_colo_BN.RData b/outputs/tfit_colo_BN.RData new file mode 100644 index 0000000000000000000000000000000000000000..d0b21d7a9892e35a8ceba66772e5109cd681db5c Binary files /dev/null and b/outputs/tfit_colo_BN.RData differ diff --git a/outputs/tfit_ext_BN.RData b/outputs/tfit_ext_BN.RData new file mode 100644 index 0000000000000000000000000000000000000000..114831e3691e6a77aa4f556e22eb8232477ace4a Binary files /dev/null and b/outputs/tfit_ext_BN.RData differ diff --git a/test-analyse-dyn.Rmd b/test-analyse-dyn.Rmd index 60e787ccfb83b504241e69bdff89a84dd8d26ca1..4a469e4468c598c55d107a6c79f7ee5589355f43 100644 --- a/test-analyse-dyn.Rmd +++ b/test-analyse-dyn.Rmd @@ -1,14 +1,18 @@ --- -title: "test vaccinium vi" +title: "test arctic analyse" author: "Isabelle Boulangeat" date: "6/01/2019" output: html_document +editor_options: + chunk_output_type: console --- # Prepare datasets ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, eval=TRUE) +library(dplyr) +library(ggplot2) ``` ## Load climate @@ -18,9 +22,12 @@ precip = read.table("shrubhub/ARCfunc_timeseries_prec.txt", h=T, sep="\t", encod 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=":") +# climato_zeroNA = as.matrix(climato[,grep("GST", names(climato))]) +# climato_zeroNA[which(is.na(climato_zeroNA[]))] = 0 +# climato[,grep("GST", names(climato))] = climato_zeroNA ``` -## Load shrub cover +## Load species cover data ```{r, fig=TRUE} load("coverc_sub.rdata") hist(coverc.sub$COVER) @@ -41,20 +48,17 @@ hist(asin(sqrt(cover.sub$COVER/100)), breaks =21, col = "lightblue") summary(cover.sub$COVER) ``` -## remove presences when low abundance (optional) + +## filter out years > 2014 because of missing in climatologies ```{r, fig=TRUE} -#---- -PA_COVER_TRESHOLD = 1 -#---- -cover.dom = cover.sub[which(cover.sub$COVER>PA_COVER_TRESHOLD),] +cover.dom = cover.sub[which(cover.sub$YEAR<2014),] ``` - ## create species-specific dataset ```{r} SPECIES = "Vaccinium vitis-idaea" -# SPECIES = "Betula nana" +SPECIES = "Betula nana" # SPECIES = "Salix arctica" # SPECIES = "Salix pulchra" # SPECIES = "Vaccinium uliginosum" @@ -62,42 +66,127 @@ SPECIES = "Vaccinium vitis-idaea" source("R_fct/shape_data.r") -spdat = create_dataset(SPECIES, cover.dom, sp.distri.path= "ArcticShrubCurrent/Vaccinium.vitisidaea_cur_pc_n.img") +spdat = create_dataset(SPECIES, cover.dom, climato=climato) +#spdat = create_dataset(SPECIES, cover.dom, sp.distri.path= "ArcticShrubCurrent/Betula.nana_cur_pc_n.img", climato=climato) ``` The steps of `create_dataset`: -- only subsites where the species have been seen +- 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 "present" and "absent" +- tag transition types - combine with climate (5 years earlier average) -- extract spatial distribution - merge all +## anomalies +```{r} +spdat$GSL_anom = spdat$GSL-spdat$GSL_av +spdat$GST_anom = spdat$GST-spdat$GST_av +spdat$snowDays_anom = spdat$snowDays-spdat$snow_days_av +``` + ## some checks ```{r} table(spdat$trType) -unique(spdat[which(spdat$st0=="colonisation"&spdat$distriArea==0),c("SUBSITE", "y0", "y1", "st0", "st1", "SoilMoist")]) +nrow(spdat) +``` +Here it means that 617 successive observations of presences are outside of the species distribution range!! +```{r} +#unique(spdat[which(spdat$trType=="present"&spdat$distriArea==0),c("SUBSITE", "y0", "y1", "st0", "st1", "SoilMoist", "distriArea")]) +``` + +# Data exploration + +```{r, fig=TRUE, fig.width=10} +spdat$trType = factor(spdat$trType) +levels(spdat$trType) <- c("AA", "COLO", "EXT", "PP") + +bp <- ggplot(na.omit(spdat), aes(fill=trType, y=snowDays, x=trType)) + + geom_point() + + ggtitle("Snow") + + facet_wrap(~SUBSITE) + + theme(legend.position="none") + + xlab("") +bp +``` + +```{r, fig=TRUE, fig.width=10} +bp <- ggplot(na.omit(spdat), aes(fill=trType, y=GST, x=trType)) + + geom_point() + + ggtitle("GST") + + facet_wrap(~SUBSITE) + + theme(legend.position="none") + + xlab("") +bp +``` + +```{r, fig=TRUE, fig.width=10} +bp <- ggplot(na.omit(spdat), aes(fill=trType, y=GSL, x=trType)) + + geom_point() + + ggtitle("GSL") + + facet_wrap(~SUBSITE) + + theme(legend.position="none") + + xlab("") +bp +``` + +```{r, fig=TRUE, fig.width=10} +bp <- ggplot(na.omit(spdat), aes(fill=trType, y=SoilMoist, x=trType)) + + geom_point() + + ggtitle("GSL") + + facet_wrap(~SUBSITE) + + theme(legend.position="none") + + xlab("") +bp +``` + +Events observed total +```{r, fig=TRUE, fig.width=10} +tab = spdat %>% + group_by(SUBSITE) %>% count(trType) + +bp <- ggplot(tab, aes(fill=trType, y=n, x=trType)) + + geom_bar(position="dodge", stat="identity") + + ggtitle("Events") + + facet_wrap(~SUBSITE) + + theme(legend.position="none") + + xlab("") +bp ``` +Events observed change +```{r, fig=TRUE, fig.width=10} +tab2 = tab[which(tab$trType%in%c("COLO", "EXT")),] +bp <- ggplot(tab2, aes(fill=trType, y=n, x=trType)) + + geom_bar(position="dodge", stat="identity") + + ggtitle("Events") + + facet_wrap(~SUBSITE) + + theme(legend.position="none") + + xlab("") +bp +``` + + # 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 = unique(spdat[, -which(colnames(spdat)=="TotalCover")]) +#spdat = na.omit(spdat) +# spdat = spdat[-which(is.nan(spdat$distriArea)),] +# spdat = spdat[-which(spdat$distriArea==0),] +selectVars = c("snowDays", "GSL", "GST", "SoilMoist", "TotalCover") spdat$SoilMoist = factor(spdat$SoilMoist) -levels(spdat$SoilMoist) = c(1:3) +levels(spdat$SoilMoist) = c(-1,0,1) spdat$SoilMoist = as.numeric(as.character(spdat$SoilMoist)) head(spdat) dim(spdat) @@ -110,48 +199,170 @@ require(LaplacesDemon) source("R_fct/model_fct.r") nvars = length(selectVars) -dat = scale(spdat[,selectVars]) +dat.ext = scale(spdat[which(spdat$trType%in%c("PP", "EXT")),selectVars]) +dat.ext[,"SoilMoist"] = spdat[which(spdat$trType%in%c("PP", "EXT")),"SoilMoist"] -parm.names = as.parm.names(list(aE=rep(0, 1+nvars), aC=rep(0, 1+nvars) )) -nbetas = length(parm.names) +parm.names.ext = as.parm.names(list(aE=rep(0, 1+nvars))) +nbetas.ext = length(parm.names.ext) -MyData <- list(N = nrow(dat), dat=dat, mon.names = c("logLik") , parm.names= parm.names, +MyData.ext <- list(N = nrow(dat.ext), dat=dat.ext, mon.names = c("logLik") , parm.names= parm.names.ext, var.aE = 1:(nvars), +aE.pos = grep("aE", parm.names.ext), +PA.pos = spdat$trType[which(spdat$trType%in%c("PP", "EXT"))]=="EXT", +PP.pos = spdat$trType[which(spdat$trType%in%c("PP", "EXT"))]=="PP", +invlogit_ = invlogit_, +itime = spdat$interval[which(spdat$trType%in%c("PP", "EXT"))], +nbetas = nbetas.ext) + +#-- +dat.colo = scale(spdat[which(spdat$trType%in%c("AA", "COLO")),selectVars]) +dat.colo[,"SoilMoist"] = spdat[which(spdat$trType%in%c("AA", "COLO")),"SoilMoist"] + +parm.names.colo = as.parm.names(list(aC=rep(0, 1+nvars))) +nbetas.colo = length(parm.names.colo) + +MyData.colo<- list(N = nrow(dat.colo), dat=dat.colo, mon.names = c("logLik") , parm.names= parm.names.colo, 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", +aC.pos = grep("aC", parm.names.colo), +AP.pos = spdat$trType[which(spdat$trType%in%c("AA", "COLO"))]=="COLO", +AA.pos = spdat$trType[which(spdat$trType%in%c("AA", "COLO"))]=="AA", invlogit_ = invlogit_, -itime = spdat$interval, -density = 1, -# density = spdat$distriArea, -# lik_fct = lik_dd, -nbetas = nbetas) +density = rep(1, nrow(dat.colo)), +itime = spdat$interval[which(spdat$trType%in%c("AA", "COLO"))], +nbetas = nbetas.colo) -MyData_N = MyData -MyData_N$density = spdat$distriArea ``` ## Fit +```{r, eval=TRUE, fig=TRUE} +summary(dat.colo) +dat.colo = data.frame(dat.colo) +dat.colo$colo.event = as.numeric(MyData.colo$AP.pos) +colo.mod = glm(colo.event~ snowDays + TotalCover, family = "binomial", data = dat.colo) + +colo.mod = glm(colo.event~ (snowDays + GSL + GST + TotalCover) * SoilMoist, data = dat.colo) +require(MASS) +stepMod.colo = stepAIC(colo.mod, direction = "both") +summary(stepMod.colo) + +pred = predict(stepMod.colo,new=dat.colo,"response") +require(fmsb) +(R2 = NagelkerkeR2(stepMod.colo)$R2) +require(ROCR) +perf = performance(prediction(pred, dat.colo$colo.event), "auc") +(AUC = perf@y.values[[1]]) +(coeff = summary(stepMod.colo)$coefficients) +(vars = rownames(coeff)[-1]) +(effect = coeff[-1,1]) +barplot(effect) +# pval = coeff[-1,4] +``` + +```{r, fig=TRUE} + +summary(dat.ext) +dat.ext = data.frame(dat.ext) +dat.ext$ext.event = as.numeric(MyData.ext$PA.pos) +ext.mod = glm(ext.event~ TotalCover, family = "binomial", data = dat.ext) +summary(ext.mod) + +require(MASS) +stepMod.ext = stepAIC(ext.mod, direction = "both") +summary(stepMod.ext) + +pred = predict(stepMod.ext,new=dat.ext,"response") +require(fmsb) +(R2 = NagelkerkeR2(stepMod.ext)$R2) +require(ROCR) +perf = performance(prediction(pred, dat.ext$ext.event), "auc") +(AUC = perf@y.values[[1]]) +(coeff = summary(stepMod.ext)$coefficients) +(vars = rownames(coeff)[-1]) +(effect = coeff[-1,1]) +barplot(effect) +``` ```{r, eval=FALSE} source("R_fct/fit_LP.r") test_model(MyData, loglik_fct) -test_model(MyData_N, loglik_fct) +test_model(MyData.ext, ext_model) +test_model(MyData.colo, colo_model) +#save.image("test-analyse.RData") +#0-- +load("test-analyse.RData") +library(LaplacesDemon) #---- ffit = first.fit(MyData, loglik_fct) Consort(ffit) sfit = second.fit(MyData, ffit,loglik_fct) Consort(sfit) -tfit = second.fit(MyData, sfit,loglik_fct) +tfit = third.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) +save.image(file= "outputs/tfit_BN.RData") +#--- +ffit = first.fit(MyData.ext, ext_model) +Consort(ffit) +sfit = second.fit(MyData.ext, ffit,ext_model) +Consort(sfit) +tfit = third.fit(MyData.ext, sfit,ext_model) +Consort(tfit) +save.image(file= "outputs/tfit_ext_BN.RData") +#--- +ffit = first.fit(MyData.colo, colo_model) +Consort(ffit) +sfit = second.fit(MyData.colo, ffit,colo_model) +Consort(sfit) +tfit = third.fit(MyData.colo, sfit,colo_model) +Consort(tfit) +save.image(file= "outputs/tfit_colo_BN.RData") + +``` + +## plot variables +```{r, fig=TRUE, fig.width=8} +load("outputs/tfit_BN.RData") +source("R_fct/graphs.r") +efCE = plot.effects(tfit, selectVars) +efCE +``` + +```{r, fig=TRUE, fig.width=8} +load("outputs/tfit_ext_BN.RData") +eff_ext = plot.effects(tfit, selectVars) +eff_ext +``` +```{r, fig=TRUE, fig.width=8} +load("outputs/tfit_colo_BN.RData") +eff_colo = plot.effects(tfit, selectVars) +eff_colo +``` +# ## compute lambda + ```{r, fig=TRUE, eval=FALSE} +# mean.posterior.aC = Fit.mcmc$Summary2[1:6, "Mean"] +# mean.posterior.aE = Fit.mcmc$Summary2[7:12, "Mean"] +# library(dplyr) +# proj.dat = unique(spdat[, c("SUBSITE", "N", "SoilMoist", "y1")]) +# Nlast = unlist(lapply(split(proj.dat, proj.dat$SUBSITE), function(x){x[which.max(x$y1),"N"]})) +# proj.dat = merge(unique(proj.dat[,-c(2,4)]), data.frame(SUBSITE = names(Nlast), N = Nlast)) +# proj.dat_all = proj.dat.shape(proj.dat, year = 2012, climato) +# proj.dat_all$SoilMoist = factor(proj.dat_all$SoilMoist) +# levels(proj.dat_all$SoilMoist) = c(1:3) +# proj.dat_all$SoilMoist = as.numeric(as.character(proj.dat_all$SoilMoist)) +# proj.dat_all = merge(proj.dat_all, unique(spdat[,c("distriArea","SUBSITE")]), by = "SUBSITE") +# +# proj.dat_dat = proj.dat_all[,vars] +# proj.dat_dat = scale(proj.dat_dat, center = attributes(dat)$`scaled:center`, scale = attributes(dat)$`scaled:scale`) +# logit_alphaC = tcrossprod(mean.posterior.aC, as.matrix(cbind(rep(1, nrow(proj.dat)), proj.dat_dat))) +# logit_alphaE = tcrossprod(mean.posterior.aE, as.matrix(cbind(rep(1, nrow(proj.dat)), proj.dat_dat))) +# +# aC.fct=as.vector(invlogit_(logit_alphaC, 1)) +# aE.fct=as.vector(invlogit_(logit_alphaE, 1)) +# +# proj.dat_all$lambda = 1-aE.fct/aC.fct +# +# bp <- ggplot(na.omit(proj.dat_all), aes(x=distriArea, y=lambda, group=distriArea)) + +# geom_boxplot(alpha=0.3, width=0.2, outlier.size = 0) + +# geom_point() +# bp ``` diff --git a/test-analyse-dyn.html b/test-analyse-dyn.html index 792797e24e34695940ff87ac874b10212a7e8f93..23afd596efd30a1f0d83986b74625b3b007e2559 100644 --- a/test-analyse-dyn.html +++ b/test-analyse-dyn.html @@ -13,7 +13,7 @@ -test vaccinium vi +test arctic analyse @@ -177,7 +177,7 @@ summary { -

test vaccinium vi

+

test arctic analyse

Isabelle Boulangeat

6/01/2019

@@ -192,10 +192,13 @@ summary { 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=":") +climato$SiteSubsite = paste(climato$SITE, climato$SUBSITE, sep=":") +# climato_zeroNA = as.matrix(climato[,grep("GST", names(climato))]) +# climato_zeroNA[which(is.na(climato_zeroNA[]))] = 0 +# climato[,grep("GST", names(climato))] = climato_zeroNA -
-

Load shrub cover

+
+

Load species cover data

load("coverc_sub.rdata")
 hist(coverc.sub$COVER)

@@ -216,17 +219,14 @@ hist(asin(sqrt(cover.sub$COVER/100)), breaks =21, col = "lightblue")## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.000 0.000 2.000 5.979 7.000 100.000
-
-

remove presences when low abundance (optional)

-
#----
-PA_COVER_TRESHOLD = 1
-#----
-cover.dom = cover.sub[which(cover.sub$COVER>PA_COVER_TRESHOLD),]
+
+

filter out years > 2014 because of missing in climatologies

+
cover.dom = cover.sub[which(cover.sub$YEAR<2014),]

create species-specific dataset

SPECIES = "Vaccinium vitis-idaea"
-# SPECIES = "Betula nana"
+SPECIES = "Betula nana"
 # SPECIES = "Salix arctica"
 # SPECIES = "Salix pulchra"
 # SPECIES = "Vaccinium uliginosum"
@@ -234,59 +234,128 @@ cover.dom = cover.sub[which(cover.sub$COVER>PA_COVER_TRESHOLD),]
source("R_fct/shape_data.r") -spdat = create_dataset(SPECIES, cover.dom, sp.distri.path= "ArcticShrubCurrent/Vaccinium.vitisidaea_cur_pc_n.img") -
## Loading required package: raster
-
## Loading required package: sp
-

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

+spdat = create_dataset(SPECIES, cover.dom, climato=climato) +#spdat = create_dataset(SPECIES, cover.dom, sp.distri.path= "ArcticShrubCurrent/Betula.nana_cur_pc_n.img", climato=climato) +

The steps of 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

+
+
+

anomalies

+
spdat$GSL_anom = spdat$GSL-spdat$GSL_av
+spdat$GST_anom = spdat$GST-spdat$GST_av
+spdat$snowDays_anom = spdat$snowDays-spdat$snow_days_av

some checks

table(spdat$trType)
## 
 ##       absent colonisation   extinction      present 
-##          205           49           28          639
-
unique(spdat[which(spdat$st0=="colonisation"&spdat$distriArea==0),c("SUBSITE", "y0", "y1", "st0", "st1", "SoilMoist")])
-
## [1] SUBSITE   y0        y1        st0       st1       SoilMoist
-## <0 rows> (or 0-length row.names)
+## 199 20 6 515 +
nrow(spdat)
+
## [1] 740
+

Here it means that 617 successive observations of presences are outside of the species distribution range!!

+
#unique(spdat[which(spdat$trType=="present"&spdat$distriArea==0),c("SUBSITE", "y0", "y1", "st0", "st1", "SoilMoist", "distriArea")])
+
+
+

Data exploration

+
spdat$trType = factor(spdat$trType)
+levels(spdat$trType) <- c("AA", "COLO", "EXT", "PP")
+
+bp <- ggplot(na.omit(spdat), aes(fill=trType, y=snowDays, x=trType)) +
+    geom_point() +
+    ggtitle("Snow") +
+    facet_wrap(~SUBSITE) +
+    theme(legend.position="none") +
+    xlab("")
+bp
+