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: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 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:
- 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)
} )
names(spDATA.list) = selection
spDATA.check = lapply(spDATA.list, check.nas)

Visualisation of transition tables

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 
## 
## $`Salix reticulata`
## 
##   AA COLO  EXT   PP 
##  491   14   24  161 
## 
## $`Geum rossii`
## 
##   AA COLO  EXT   PP 
##   16    4    8  153 
## 
## $`Empetrum nigrum`
## 
##   AA COLO  EXT   PP 
##  448   18   23  133 
## 
## $`Bistorta bistortoides`
## 
##   AA COLO  EXT   PP 
##   43   20   22   96 
## 
## $`Hierochloe alpina`
## 
##   AA COLO  EXT   PP 
##  603   39   36   61

Here we can see how many transitions we have for each species.

Fit models

spDATA.fit = lapply(spDATA.list, fit.species)

Summary of model evaluation

eval.models = lapply(spDATA.fit, function(glm.fit){
  eval.colo = eval.fit(glm.fit$colo, glm.fit$dat$colo, "colo.event")[-1]
  coeffs.colo = sort(abs(glm.fit$colo$coefficients[-1]), dec=T)
  vars.colo = paste(names(coeffs.colo), collapse=";")
  eval.ext = eval.fit(glm.fit$ext, glm.fit$dat$ext, "ext.event")[-1]
  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))
})

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)
model evaluation
colo.R2 colo.AUC ext.R2 ext.AUC
Vaccinium vitis-idaea 0.42 0.89 0.03 0.61
Carex bigelowii 0.31 0.79 0.08 0.64
Ledum palustre 0.14 0.67 0.04 0.60
Betula nana 0.29 0.85 0.18 0.88
Eriophorum vaginatum 0.02 0.53 0.07 0.70
Salix pulchra 0.18 0.75 0.03 0.59
Eriophorum chamissonis 0.22 0.74 0.06 0.65
Persicaria bistorta 0.11 0.67 0.08 0.63
Dupontia fisheri 0.24 0.80 0.25 0.80
Rubus chamaemorus 0.00 0.50 0.00 0.50
Salix reticulata 0.18 0.75 0.21 0.73
Geum rossii 0.25 0.70 0.28 0.86
Empetrum nigrum 0.16 0.77 0.11 0.67
Bistorta bistortoides 0.20 0.71 0.04 0.59
Hierochloe alpina 0.18 0.77 0.28 0.77
# 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).

Selected variables and coefficient, sumnmary table

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)
# str(coeffs.all)

Number of selected variable per process, accross all species

data.frame(unclass(table(coeffs.all$var, coeffs.all$process)))  %>% 
  kable(.) %>%
  kable_styling("striped")
colonisation extinction
GSL 11 6
GSL:SoilMoist 5 3
GST 9 6
GST:SoilMoist 5 2
SoilMoist 13 10
TotalCover 11 7
TotalCover:SoilMoist 9 4
snowDays 10 10
snowDays:SoilMoist 7 4

Variable effect coefficients

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

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")

Lambda (1-ext/colo)

An overview of the performance of populations by site of by species

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

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)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'