Assessing the predictive ability of Demoniche

Fiona Spooner

September 6, 2017

install.packages('climates',,'http://www.rforge.net/')
install.packages("dismo")
install.packages("maps")
install.packages("mapdata")
install.packages("popbio")
install.packages("demoniche", repos="http://R-Forge.R-project.org")
library(dismo)
library(demoniche)
library(maps)
library(mapdata)
library(rgeos)
library(zoo)
library(raster)
library(mgcv)
library(randomForest)
library(rgdal)
library(doParallel)
library(ggplot2)
library(dplyr)
library(data.table)
library(reshape2)
library(taRifx)
library(plyr)
library(sp)
wd<-getwd()
genus<-"Capra"
species<-"ibex"
Europe_only<-TRUE
binomial<-paste(genus, species, sep="_")
min_lat<- 43
max_lat<- 47.9
min_lon<- 0
max_lon<- 16
bioclim_layers<-c(1,5,6,13,15,18,19)   
no_background_points<-1000
bioclim_names<-c("Bio_1_2006_2016_average", "Bio_5_2006_2016_average","Bio_6_2006_2016_average","Bio_13_2006_2016_average","Bio_15_2006_2016_average","Bio_18_2006_2016_average","Bio_19_2006_2016_average")


#Comadre

#niche values
transition_affected_niche<-c(2,11,20)  #which parts of the matrix are affected by the c(2,11,20) is  juveniles
transition_affected_env <- c(2,11,20)
transition_affected_demogr <- c(2,11,20)

transition_affected_niche<-"all"  #which parts of the matrix are affected by the c(2,11,20) is  juveniles
transition_affected_env <- "all"
transition_affected_demogr <- "all"

env_stochas_type<-"normal"   #can also be lognormal
#standard deviation of matrices
#- just doing eqaul splits for now
#proportion_initialf<- c(1/3,1/3,1/3)
#density_individuals <- 4000  #4292.32 to 16096.2 based on density being between 8 and 30 per 100 ha and the area of each cell being 53654 ha - 
carry_k<-5000

  #Things to vary
sd_seq<-0.5
LDD_seq<-c(0.9)
kern_seq<-list(c(1,1), c(19775,113000),c(13597.5, 77700), c(1000000,1000000))

var_grid<-expand.grid(sd_seq, LDD_seq, kern_seq)
colnames(var_grid)<-c("SD", "LDD", "Kern")

kern_mat<-matrix(c(rep(c(1,1), (nrow(var_grid)/ncol(var_grid))),rep(c(19775,113000), (nrow(var_grid)/ncol(var_grid))),rep(c(13597.5, 77700), (nrow(var_grid)/ncol(var_grid))),rep(c(1000000,1000000),(nrow(var_grid)/ncol(var_grid)))), ncol=2, byrow=T)

density_mine<-2671
fraction_LDD<-0.1
fraction_SDD<-0.1

reps<-12
wd<-getwd()
binomial<-paste(genus, species, sep="_")

k<-4
years<-1950:2016
spin_years<-1650:1949
lpi<-read.csv("LPI_pops_20160523_edited.csv")
species_directory<-paste(wd, binomial, sep="/")
dir.create(species_directory)
## Warning in dir.create(species_directory): 'D:\Fiona\Git_Method\Git_Method
## \Capra_ibex' already exists
sdm_folder<-paste(species_directory, "SDM_folder", sep = "/")
dir.create(sdm_folder)
## Warning in dir.create(sdm_folder): 'D:\Fiona\Git_Method\Git_Method
## \Capra_ibex\SDM_folder' already exists
demoniche_folder<-paste(species_directory, "Demoniche_Output", sep = "/")
dir.create(demoniche_folder)
## Warning in dir.create(demoniche_folder): 'D:\Fiona\Git_Method\Git_Method
## \Capra_ibex\Demoniche_Output' already exists
source("demoniche_setup_me.R")
source("demoniche_model_me.R")
source("ensemble_evaluate.R")
#source("demoniche_dispersal_function.R")

no_yrs_mine<-1 
prob_scenario<-c(0.5,0.5)    #need to check this
noise<-0.90
Bioclim Variables

Creating the 19 bioclim variables for each year between 1950 & 2016 for Europe using the E-OBS dataset - downloaded from http://www.ecad.eu/download/ensembles/download.php

rr<-brick(paste(wd, "/rr_0.25deg_reg_v15.0.nc", sep="")) #precipitation
tg<-brick(paste(wd, "/tg_0.25deg_reg_v15.0.nc", sep="")) #mean temp
tn<-brick(paste(wd, "/tn_0.25deg_reg_v15.0.nc", sep="")) #min temp
tx<-brick(paste(wd, "/tx_0.25deg_reg_v15.0.nc", sep="")) #max temp
#dates<-as.Date((gsub("X", "",names(rr))), format="%Y.%m.%d")


rr_mon<-zApply(rr, by=as.yearmon, fun = mean)
tg_mon<-zApply(tg, by=as.yearmon, fun = mean)
tn_mon<-zApply(tn, by=as.yearmon, fun = mean)
tx_mon<-zApply(tx, by=as.yearmon, fun = mean)

jans<-seq(1,804, by=12)
years<-as.character(1950:2016)

for (i in 1:length(years)){
  jan_sel<-jans[i]
  dec_sel<-jan_sel+11
  bio_vars_all<-biovars(rr_mon[[jan_sel:dec_sel]], tn_mon[[jan_sel:dec_sel]], tx_mon[[jan_sel:dec_sel]])  
  writeRaster(bio_vars_all, paste(paste(wd, "/Bioclim/", years[i],"_bioclim_variable_stack.tif", sep=""), overwrite=T))
  print(years[i])
  }
path<-"D:/Fiona/Git_Method/Git_Method/CRUTEM"
files<-list.files(path)
prec_files<-files[grepl("pre",files )]
tmx_files<-files[grepl("tmx",files )]
tmn_files<-files[grepl("tmn",files )]

prec<-stack(paste(path, "/", prec_files, sep=""))
tmin<-stack(paste(path, "/", tmn_files, sep=""))
tmax<-stack(paste(path, "/", tmx_files, sep=""))

jans<-seq(1,nlayers(prec), by=12)
years<-as.character(1954:2016)

for (i in 1:length(years)){
  jan_sel<-jans[i]
  dec_sel<-jan_sel+11
  bio_vars_all<-biovars(prec[[jan_sel:dec_sel]], tmin[[jan_sel:dec_sel]], tmax[[jan_sel:dec_sel]])  
  writeRaster(bio_vars_all, paste(wd, "/Bioclim/Global/", years[i],"_bioclim_variable_stack.tif", sep=""), overwrite=T)
  
  print(years[i])
  
}

Species distribution modelling with dismo

Extracting Alpine ibex (Capra ibex) data from GBIF and with additional locations added from the LPI dataset. These points are then formatted to a SpatialPointsDataFrame.

wd<- getwd()
capra <-gbif(genus =genus, species = species)
capgeo <- subset(capra, !is.na(lon) & !is.na(lat) & (lon!="NA" & lat !="NA") | year!="NA") 
dups <- duplicated(capgeo[, c("lon", "lat")])
capg <-capgeo[!dups, ]
capg2 <- capg[capg$lon > min_lon & capg$lon< max_lon & capg$lat > min_lat & capg$lat < max_lat & capg$year>=2006, ] 
capg2<-data.frame(capg2$lon,capg2$lat)
colnames(capg2)<-c("Longitude", "Latitude")

pyr<-subset(lpi,Binomial == binomial & Specific_location==1)    #record 11470 had wrong longitude - in Russia!
pyrs<-pyr[,c("Longitude","Latitude")]
capg2<-rbind(capg2, pyrs)
capg2<-na.omit(capg2)
capg2$presence<-rep(1)
capg2$ID<-1:nrow(capg2)

capc<-as.matrix(capg2[ , c( "Longitude","Latitude", "presence")])
capc<-matrix(capc[complete.cases(capc)], ncol=3)
xy<-as.matrix(capc[,c(1,2)])
df<-data.frame(capc[,3])

write.csv(df, paste(species_directory,"/", binomial,"_gbif.csv", sep=""), row.names = FALSE)
write.csv(xy, paste(species_directory,"/",binomial, "_locs.csv", sep=""), row.names=FALSE)
library(sp)
library(raster)

df<-read.csv(paste(species_directory, "/",binomial, "_gbif.csv", sep="") )
xy<-read.csv(paste(species_directory,"/" ,binomial, "_locs.csv", sep=""))

sp<-sp:::SpatialPointsDataFrame(coords=xy, data=df)

e<-raster:::extent(sp)

e_big<-e+2

maps:::map('world',  col="light grey", fill=T, xlim=c(e_big[1],e_big[2]), ylim=c(e_big[3],e_big[4]))

points(sp, col="red", pch=20)

Creating the 2006-2016 average for each bioclim variable and then selecting out the layers we want to use in the model:

  • Bioclim 1 – Annual mean temperature

  • Bioclim 5 – Max temperature of the warmest month

  • Bioclim 6 – Min temperature of the coldest month

  • Bioclim 13 – Precipitation of the wettest month

  • Bioclim 15 – Precipitation seasonality (coefficient of variation)

  • Bioclim 18 – Precipitation of warmest quarter

  • Bioclim 19 – Precipitation of coldest quarter

Snow depth seems to be an important predictor for alpine ibex but I have struggled to find historical data on this. Existing papers are based on historical data from one weather station, which means there is no spatial variation and therefore it is not appropriate for these models.

if (Europe_only){
lf<-list.files(paste(wd, "/Bioclim/", sep=""))  
} else{
lf<-list.files(paste(wd, "/Bioclim/Global", sep=""))  
}


first<-which(grepl("2006_bioclim", lf)) #was initially 1985
last<-which(grepl("2016_bioclim", lf))

if (Europe_only){
all_years<-stack(paste(wd, "/Bioclim/",lf[first:last], sep=""))
} else{
all_years<-stack(paste(wd, "/Bioclim/Global/",lf[first:last], sep=""))
}

bios<-seq(1,nlayers(all_years), by=19)
cellStats(all_years[[bios]], stat="mean")
#creating a 1985-2016 average of each bioclim variable
for (i in 1:19){
  layers<-bios
  bio_layer<-mean(all_years[[layers]])
  if (Europe_only){
  writeRaster(bio_layer, paste(wd, "/Bioclim/Bio_",i,"_2006_2016_average.tif",sep=""), overwrite=TRUE)
  } else {
  writeRaster(bio_layer, paste(wd, "/Bioclim/Global/Bio_",i,"_2006_2016_average.tif",sep=""), overwrite=TRUE)
  }
  bios<-bios+1
  #print(layers)
}

#pred_nf<-stack(paste(wd, "/Bioclim/Bio_", bio_layer_pred,"_1985_2016_average.tif",sep="" ))

Selecting out the bioclim layers I’m interested in:

if (Europe_only){
pred_nf<-stack(paste(wd, "/Bioclim/Bio_", bioclim_layers,"_2006_2016_average.tif",sep="" ))
} else {
pred_nf<-stack(paste(wd, "/Bioclim/Global/Bio_", bioclim_layers,"_2006_2016_average.tif",sep="" ))
}

pred_crop<-crop(pred_nf, e_big)
plot(pred_crop)

Creating the testing and training datasets

Using K=4 so 75% of points are used for training and 25% for testing.

library(dismo)

set.seed(10)

group_pres<-kfold(sp, k)
#write.csv(group_pres, paste(species_directory, "k_folds_presence_2006.csv", sep="/"))
group_pres<-read.csv(paste(species_directory, "k_folds_presence_2006.csv", sep="/"))
group_pres<-group_pres[,-1]  #all presence points

sp_train<-sp[group_pres!=1,] #presence points split into training and test
sp_test<-sp[group_pres ==1,]
bg_sam<- randomPoints(pred_crop[[1]], no_background_points )


#write.csv(bg_sam, paste(species_directory, "background_random_points2.csv", sep="/"))
bg_pse<-read.csv(paste(species_directory, "background_random_points2.csv", sep="/"))
bg_pse<-bg_pse[,-1]

group_back<-kfold(bg_pse, k)
#write.csv(group_back, paste(species_directory, "k_folds_background2.csv", sep="/"))
group_back<-read.csv(paste(species_directory, "k_folds_background2.csv", sep="/"))
group_back<-group_back[,-1]
sp_backg_train <- bg_pse[group_back != 1, ] #background points split into training and test
sp_backg_test <- bg_pse[group_back == 1, ]
colnames(sp_backg_train)<-c("lon", "lat")
colnames(sp_backg_test)<-c("lon", "lat")

sp_train<-sp_train@coords #presence training points
colnames(sp_train)<-c("lon", "lat")
train <- rbind(sp_train, sp_backg_train)  #presence and background training points combined

sp_test<-sp_test@coords #presence test points
colnames(sp_test)<-c("lon", "lat")
test<-rbind(sp_test, sp_backg_test)    #presence and background test points combined

pb_train <- c(rep(1, nrow(sp_train)), rep(0, nrow(sp_backg_train)))  
envtrain <- extract(pred_nf, train)
envtrain <- data.frame( cbind(pa=pb_train, envtrain) ) #training points with environmental  variables

envtrain<-na.omit(envtrain)

pb_test<-c(rep(1, nrow(sp_test)), rep(0, nrow(sp_backg_test)))
envtest <- extract(pred_nf, test)
envtest <- data.frame( cbind(pa=pb_test, envtest) ) #test points with environmental  variables
envtest<-na.omit(envtest)

env_all<-rbind(envtrain, envtest) #all of the training and test data with environmental variables

env_pres<-extract(pred_nf, sp)
pa<-rep(1, nrow(env_pres))
env_pres<-data.frame(pa, env_pres)  #presence points with environmental variables
colnames(env_pres)<-c("pa", bioclim_names)
env_pres_xy<-data.frame(sp@coords, env_pres)
env_pres_xy<-na.omit(env_pres_xy)

env_back<-extract(pred_nf, bg_pse)
ap<-rep(0, nrow(env_back))
env_back<-data.frame(ap, env_back) #background points with environmental variables
colnames(env_back)<-c("pa",bioclim_names)
env_back_xy<-data.frame(bg_pse, env_back)
env_back_xy<-na.omit(env_back_xy)

group_pres<-group_pres[1:nrow(env_pres_xy)]
group_back<-group_back[1:nrow(env_back_xy)]

Bioclim

Running the Bioclim envelope model and evaluating it using three different thresholds (kappa, no omission and true skill statistic) to get AUC values.

evl_bc<- list()
tss_bc<-list()
for (i in 1:k){
  pres_train<-sp[group_pres!=i,]
  pres_test<-sp[group_pres==i,]
  backg_test<-bg_pse[group_back==i,]
  bc <- bioclim(pred_nf,pres_train)
  evl_bc[[i]] <- dismo:::evaluate(pres_test, backg_test, bc,pred_nf,type="response")#test presence, test absence, model, predictor variables 
     #print(i)
}

auc_bc <- sapply( evl_bc, function(x){slot(x, "auc")} )
print(auc_bc)
## [1] 0.8360204 0.9246301 0.8916439 0.8599503
bc_auc<-mean(auc_bc)

GAM

Running a GAM and evaluating it to get an AUC value using three different thresholds (kappa, no omission and true skill statistic) to get AUC values.

#library(mgcv)

evl_gam<- list()
tss_gam<-list()
for (i in 1:k){
  pres_train<-env_pres_xy[group_pres!=i ,-c(1,2)]
  pres_test<-env_pres_xy[(group_pres==i) ,-c(1,2)]
  back_test<-env_back_xy[(group_back==i),-c(1,2)]
  back_train<-env_back_xy[(group_back!=i),-c(1,2)]
  envtrain<-rbind(pres_train, back_train)
  gm1<-gam(pa~ s(Bio_1_2006_2016_average)+s(Bio_5_2006_2016_average)+ s(Bio_6_2006_2016_average)+ 
           s(Bio_13_2006_2016_average)+ s(Bio_15_2006_2016_average)+ s(Bio_18_2006_2016_average)+ 
           s(Bio_19_2006_2016_average), family = binomial(link = "logit"),data=envtrain)

  evl_gam[[i]] <- dismo:::evaluate(p = pres_test, a = back_test,model= gm1,type="response")
    }

auc_gam <- sapply( evl_gam, function(x){slot(x, "auc")} )
print(auc_gam)
## [1] 0.9447715 0.9569424 0.9669063 0.9315718
gam_auc<-mean(auc_gam)

Random Forest

Running a random forest model and evaluating it using three different thresholds (kappa, no omission and true skill statistic) to get AUC values.

model<-pa~  Bio_1_2006_2016_average+Bio_5_2006_2016_average+ 
  Bio_6_2006_2016_average+ Bio_13_2006_2016_average+ Bio_15_2006_2016_average+ 
  Bio_18_2006_2016_average+ Bio_19_2006_2016_average

evl_rf<- list()
for (i in 1:k){
  pres_train<-env_pres_xy[group_pres!=i,-c(1,2)]
  pres_test<-env_pres_xy[(group_pres==i) ,-c(1,2)]
  back_test<-env_back_xy[(group_back==i),-c(1,2)]
  back_train<-env_back_xy[group_back !=i,-c(1,2)]
  envtrain<-rbind(pres_train, back_train)
  rf1 <- randomForest(model, data=envtrain)
  evl_rf[[i]] <- dismo:::evaluate(pres_test, back_test, rf1,type="response")
  }

auc_rf <- sapply( evl_rf, function(x){slot(x, "auc")} )
print(auc_rf)
## [1] 0.9615924 0.9735659 0.9878706 0.9676302
rf_auc<-mean(auc_rf)

Ensemble Model

Creating a weighted (based on AUC) ensemble average suitability model for 2006-2016 from the Bioclim, GAM and Random Forest models. This will be used to predict habitat suitability for Alpine ibex for each year 1950-2016.

#total models
library(randomForest)
bc <- bioclim(pred_nf, sp)
gm1<-gam(pa~ s(Bio_1_2006_2016_average)+s(Bio_5_2006_2016_average)+ s(Bio_6_2006_2016_average)+ 
           s(Bio_13_2006_2016_average)+ s(Bio_15_2006_2016_average)+ s(Bio_18_2006_2016_average)+ 
           s(Bio_19_2006_2016_average), family = binomial(link = "logit"),data=env_all)
rf1 <- randomForest(model, data=env_all)

pb <- predict(pred_nf, bc, ext=e_big, progress='') #bioclim predict
pg <- predict(pred_nf, gm1, ext=e_big) #gam predict
pgl<-raster:::calc(pg, fun=function(x){ exp(x)/(1+exp(x))}) #backtransforming from logit space
pr <- predict(pred_nf, rf1, ext=e_big) #random forest predict

models <- stack(pb, pgl, pr)
names(models) <- c("bioclim", "gam", "random forest")
plot(models)

auc<-c(bc_auc, gam_auc, rf_auc)
w <- (auc-0.5)^2
wm <- weighted.mean( models[[c("bioclim", "gam", "random.forest")]], w)
plot(wm, main="Ensemble model")
points(sp)

Weighted Threshold

abs_xy<-env_back_xy[,c(1,2)]

ens<-ensemble_evaluate(sp,abs_xy , wm)
thresh<-threshold(ens, stat="spec_sens")

thresh
## [1] 0.2512174

Annual Habitat Suitability Predictions

Using the ensemble model to predict suitability for each year - 1950-2016 and for each year creating a binary presence/absence map based on a three different thresholding techniques. I will go forward using true skill statistic based on Allouche 2006.

bc <- bioclim(pred_nf, sp)
gm1<-mgcv:::gam(pa~ s(Bio_1_2006_2016_average)+ s(Bio_5_2006_2016_average)+ s(Bio_6_2006_2016_average)+ s(Bio_13_2006_2016_average)+ s(Bio_15_2006_2016_average)+s(Bio_18_2006_2016_average)+ s(Bio_19_2006_2016_average), data=env_all)
rf1 <- randomForest:::randomForest(model, data=env_all)

for (i in 1:length(years)){
  
  pred_nf<-stack(paste(wd, "/Bioclim/", years[i], "_bioclim_variable_stack.tif", sep="" ))  
  pred_nf<-pred_nf[[bioclim_layers]]
  names(pred_nf)<-bioclim_names
  
  pb <- predict(pred_nf, bc, ext=e_big, progress='')
  pg <- predict(pred_nf, gm1, ext=e_big) #gam predict
  pgl<-raster:::calc(pg, fun=function(x){ exp(x)/(1+exp(x))})
  pr <- predict(pred_nf, rf1, ext=e_big)
  
  models <- stack(pb, pgl, pr)
  names(models) <- c("bioclim", "gam", "random forest")
  wm <- weighted.mean( models[[c("bioclim", "gam", "random.forest")]], w)
  
  pa_sss<-wm>thresh
 
 writeRaster(wm , paste(sdm_folder, "/weighted_ensemble_sdm_", years[i], ".tif", sep=""), overwrite=TRUE)
 writeRaster(pa_sss , paste(sdm_folder, "/pres_abs_sss_weighted_ensemble_sdm_", years[i], ".tif", sep=""), overwrite=TRUE)

  #plot(pa_sss, main=years[i])
}

Patches

Plots of the number of cells with predicted presence and number of patches (contiguous presence) over 1950-2016.

Demoniche

Formatting the population occurrence points

Formatting the LPI population data for use in Demoniche, they start as the “seed” populations. Might be better to use GBIF data for this - unrealistic that the LPI populations are the only existing populations in 1950 - alternatively a historical map of where the ibex were in 1950?

pyr<-subset(lpi, Binomial ==binomial & Specific_location==1)    #record 11470 had wrong longitude - in Russia!

#formatting the data for use in demoniche
pyrs<-pyr[,c("ID","Longitude","Latitude")]

pyrs$ID<-pyrs$ID * 100

xy<-data.frame(pyrs$Longitude, pyrs$Latitude)
coordinates(xy)<-c("pyrs.Longitude", "pyrs.Latitude")
proj4string(xy) <- CRS("+init=epsg:4326") # WGS 84
CRS.new <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
xy <- spTransform(xy, CRS.new)
xy<-data.frame(xy)

Populations<-data.frame(pyrs$ID, xy$pyrs.Longitude, xy$pyrs.Latitude)
colnames(Populations)<-c("PatchID", "X", "Y")
Populations$area<-1

pop_years<-pyr[,65:130]
pop_years[pop_years == "NULL"]<-NA

first_year<-function(x){
  pop_first<-min(which(!is.na(x)))
  pop_value<-x[pop_first]
  return(pop_value)
}
pop_values<-apply(pop_years, 1, first_year)
density_mine<-as.numeric(pop_values)

# id<-pyrs$ID*100
# lam<-rep(1,length(id))    #not sure what the value here pertains to - think it sets starting population so should use values from LPI?
# pyrxy<-SpatialPoints(pyr[,c("Longitude","Latitude")])
# 
# wd<-getwd()
# sdm<-raster(paste(sdm_folder, "/pres_abs_sss_weighted_ensemble_sdm_1950.tif", sep=""))
# e2<-extent(sdm)
# 
# r<-raster(e2, resolution=res(sdm))
# 
# rz<-rasterize(pyrxy,r,lam )
# crs(rz)<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
# rz<-projectRaster(rz, crs ="+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
# 
# rid<-rasterize(pyrxy,r,id)
# crs(rid)<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
# rid<-projectRaster(rid, crs ="+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
# 
# rz_spdf<-xyFromCell(rz, 1:ncell(rid))
# 
# rzm<-as.vector(rz)
# ridm<-as.vector(rid)
# 
# df<-data.frame(ridm,rz_spdf,rzm)
# colnames(df)<-c( "PatchID","X","Y","area")
# 
# Populations<-data.frame(na.omit(df))

Formatting the habitat suitability models

Formatting the habitat suitability model maps into vectors for use in Demoniche and plotting the average habitat suitability over time. Here we also create a fake set of habitat suitability models where the suitability is driven down at a constant rate over time.

sdm_patch_df<-data.frame(ID=1:ncell(patch))
sdm_df<-data.frame(ID=1:ncell(patch))
sdm_fake_df<-data.frame(ID=1:ncell(patch))

#formatting data for demoniche
for (i in 1:length(years)){
  
  #a selection of different threshold techniques for presence absence, as well as a suitability surface
  sdm<-raster(paste(sdm_folder,"/weighted_ensemble_sdm_", years[i],".tif", sep=""))   #14.6
  patch<-raster(paste(sdm_folder,"/pres_abs_sss_weighted_ensemble_sdm_", years[i],".tif", sep="")) 
  
  sdm<-projectRaster(sdm, crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")
  patch<-projectRaster(patch, crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
  #21.1%  
  #fake<-raster(paste(wd, "/Alp_SDMs/Fake_Landscape/Fake_Landscape_", years[i], ".tif", sep=""))
  if (i ==1){
    vec<-as.data.frame(sdm, xy = TRUE)
    vec_pat<-as.data.frame(patch, xy=TRUE)
    #vec_fake<-as.data.frame(fake, xy=TRUE)
  } else{
    vec<-as.data.frame(sdm)
    vec_pat<-as.data.frame(patch)
    #vec_fake<-as.data.frame(fake)
  }
  
  sdm_df<-cbind(sdm_df, vec)
  sdm_patch_df<-cbind(sdm_patch_df,vec_pat)
  #sdm_fake_df<-cbind(sdm_fake_df, vec_fake)
  #print(i)
}

sdm_df$ID[which(!is.na(df$PatchID))]<-df$PatchID[!is.na(df$PatchID)]
sdm_patch_df$ID[which(!is.na(df$PatchID))]<-df$PatchID[!is.na(df$PatchID)]
#sdm_fake_df$ID[which(!is.na(df$PatchID))]<-df$PatchID[!is.na(df$PatchID)]

niche_map_mine<-sdm_df
colnames(niche_map_mine)[1:3]<-c("gridID", "X", "Y")

patch_map_mine<-sdm_patch_df
colnames(patch_map_mine)[1:3]<-c("gridID", "X", "Y")

# fake_map_mine<-sdm_fake_df
# colnames(fake_map_mine)[1:3]<-c("gridID", "X", "Y")

niche_spin_up<-matrix(rep(niche_map_mine[,4], length(spin_years)), nrow=nrow(niche_map_mine))
niche_spin_up<-cbind(niche_map_mine[,1:3],niche_spin_up, niche_map_mine[,4:ncol(niche_map_mine)])

patch_spin_up<-matrix(rep(patch_map_mine[,4], length(spin_years)), nrow=nrow(patch_map_mine))
patch_spin_up<-cbind(patch_map_mine[,1:3],patch_spin_up, patch_map_mine[,4:ncol(patch_map_mine)]) #last patch used to be niche

# fake_spin_up<-matrix(rep(fake_map_mine[,4], 99), nrow=nrow(fake_map_mine))
# fake_spin_up<-cbind(fake_map_mine[,1:3],fake_spin_up, fake_map_mine[,4:ncol(fake_map_mine)])


col_years_short<-paste("Year_", years, sep="")
col_years<-paste("Year_", min(spin_years):max(years), sep="")

colnames(niche_map_mine)[4:length(colnames(niche_map_mine))]<-col_years_short
colnames(patch_map_mine)[4:length(colnames(patch_map_mine))]<-col_years_short
#colnames(fake_map_mine)[4:length(colnames(fake_map_mine))]<-col_years_short

colnames(niche_spin_up)[4:length(colnames(niche_spin_up))]<-col_years
colnames(patch_spin_up)[4:length(colnames(patch_spin_up))]<-col_years
#colnames(fake_spin_up)[4:length(colnames(fake_spin_up))]<-col_years

plot(years,colMeans(na.omit(niche_map_mine)[4:length(colnames(niche_map_mine))]), type="l", ylab="Mean suitability index")

niche_map_mine<-na.omit(niche_map_mine)
patch_map_mine<-na.omit(patch_map_mine)
#fake_map_mine<-na.omit(fake_map_mine)

niche_spin_up<-na.omit(niche_spin_up)
patch_spin_up<-na.omit(patch_spin_up)
#fake_spin_up<-na.omit(fake_spin_up)

#niche_map_mine[is.na(niche_map_mine)]<-0
patch_map_mine[is.na(patch_map_mine)]<-0

COMADRE

Accessing the population matrices from COMADRE - there are 5 matrices available for the alpine ibex but four are for the same area and only represent survival over two time periods in good/bad conditions.

Could I use these to calibrate the survival aspect of the matrices?

load(paste(wd, "COMADRE_v.2.0.1.RData", sep="/"))

tempMetadata<-subset(comadre$metadata, SpeciesAccepted==binomial)

keep<-as.numeric(rownames(tempMetadata))

tempMat<-comadre$mat[keep]   #MatA is matrix pop model, can be split into U, F and/or C

MatList<-list(tempMat[[1]][[1]])  #varies depending on number of matrices - need to find a way to code this better - now have five matrices available so need to sort this
AllMat<-unlist(MatList)
matrices<-matrix(AllMat, ncol=length(MatList))
colnames(matrices)<- c("Reference_matrix")

    #not yet active in demoniche

stages<-comadre$matrixClass[keep][[1]]$MatrixClassAuthor
#stages<-comadre$matrixClass[keep][[2]]$MatrixClassAuthor
#stagesf<-stages[1:3]
proportion_initial<- rep(1/length(stages), length(stages)) #Think spin up sorts this
sumweight<-rep(1, length(stages))#weight of stages  - should be equal for all mine just 
#in plants seed not included in calculating population sizes  
list_names_matrices<-colnames(matrices)
K_weight<-c(rep(1, length(stages)))  #the weight with which carrying capacity affects each stage was FALSE

carrying capacity function

lf<-list.files(sdm_folder)   

files<-lf[grepl("^weighted_ensemble_sdm_.*.tif$", lf)]

sdms<-stack(paste(sdm_folder, files, sep="/"))
sdms<-projectRaster(sdms, crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")

hsi<-extract(sdms, Populations[,c(2,3)])

link<-lin(hsi)
#sigk<-sig(hsi)
#ltk<-lt(hsi)

spin<-replicate(length(spin_years),link[,1])
link_spin<-cbind(spin, link)

colnames(link_spin)[1:length(spin_years)]<-paste("weighted_ensemble_sdm_", spin_years, sep="")

Varying standard deviation of population matrix, long distance dispersal and kernel shape

Changed the demoniche_setup function (now demoniche_setup_me.R) to change the dispersal function to a decay curve where the half-life is the median estimated dispersal distance.This was based on Bowman et al 2012 where it suggests that the maximum dispersal is 40x the square root of a species home range and median dispersal is 7x the square root of HR.

## [1] "patch_old_disp1"
## [1] "patch_old_disp1000000"
## [1] "patch_old_disp113000"
## [1] "patch_old_disp77700"
#rep_df_all<-read.csv(paste(species_directory, "kern_kap_new_params_0_5sd.csv", sep="/"))
rep_df_all<-read.csv(paste(species_directory, "my_disp_function.csv", sep="/"))
#rep_df_all<-read.csv(paste(species_directory, "my_disp_dens_function.csv", sep="/"))
#rep_df_all<-read.csv(paste(species_directory, "my_disp_only_function.csv", sep="/"))
rep_df_all$rep_id<-factor(rep_df_all$rep_id)
rep_df_all$md_id<-factor(rep_df_all$md_id)

rep_df_all$trend<-as.numeric(as.character(rep_df_all$trend))
library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:raster':
## 
##     calc
p01<-ggplot(rep_df_all, aes(x=year, y=trend, group=interaction(rep_id, md_id),colour=md_id))+
  geom_line()+ 
   theme(text = element_text(size=9))+ 
  labs(color='Maximum \ndispersal \ndistance') 
p01

rep_df_small<-rep_df_all[rep_df_all$md_id != "1000000" & rep_df_all$year>1949,]
library(ggplot2)

p01<-ggplot(rep_df_small, aes(x=year, y=trend, group=interaction(rep_id, md_id),colour=md_id))+
  geom_line()+ 
   theme(text = element_text(size=9))+ 
  labs(color='Maximum \ndispersal \ndistance') 
p01

library(reshape2)
## Warning: package 'reshape2' was built under R version 3.4.3
l<-list.files(demoniche_folder)
nf<-length(list.files(paste(demoniche_folder, l[1], sep="/")))

highfoldernames<-list.files(demoniche_folder)
highfoldernames<-highfoldernames[11:14]
lowfoldernames<-rep(1:nf, each=length(l))

foldernames<-paste(highfoldernames, lowfoldernames, sep="/")


sp_lpi<-lpi[lpi$Binomial == binomial & lpi$Specific_location ==1,]

xy<-data.frame(sp_lpi$Longitude, sp_lpi$Latitude)
coordinates(xy)<-c("sp_lpi.Longitude", "sp_lpi.Latitude")
proj4string(xy) <- CRS("+init=epsg:4326") # WGS 84
CRS.new <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
xy <- spTransform(xy, CRS.new)

convert_pop_out<-function(foldername){
  
  pop_out<-read.csv(paste(demoniche_folder ,foldername, "Reference_matrix_pop_output.csv", sep="/"), header = TRUE)
  pop_out<-pop_out[,-1]
  coordinates(pop_out) <- ~ X + Y
  gridded(pop_out) <- TRUE
  rasterDF <- stack(pop_out)
  proj4string(rasterDF)<-CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
  trends<-extract(rasterDF,xy)
  max_disp<-strsplit(foldername, "[/_]")[[1]][3]
  rep_id<-strsplit(foldername, "[/_]")[[1]][4]
  trends_df<-data.frame(sp_lpi$ID,max_disp, rep_id,trends)
}

demoniche_pop_out<-lapply(foldernames, convert_pop_out)
df <- do.call("rbind", demoniche_pop_out)
dfm<-as.matrix(df)

lambda<-function(x){

  l10<-10^diff(log10(as.numeric(x[4:length(x)])))
  
}

dft<-t(apply(dfm,1,lambda))

df_lambda<-data.frame(dfm[,1:3],dft)

colnames(df_lambda)[4:ncol(df_lambda)]<-colnames(dfm)[5:ncol(dfm)]

melt_df<-melt(df, id=1:3)
melt_df$year<-as.numeric(gsub("Year_", "", melt_df$variable))

melt_lambda<-melt(df_lambda, id=1:3)
melt_lambda$year<-as.numeric(gsub("Year_", "", melt_lambda$variable))

melt_short<-melt_df[melt_df$year>spin_years[length(spin_years)] ,]
melt_short$sp_lpi.ID<-as.factor(melt_short$sp_lpi.ID)

melt_lambda_short<-melt_lambda[melt_lambda$year>years[1] ,]
melt_lambda_short$sp_lpi.ID<-as.factor(melt_lambda_short$sp_lpi.ID)

ggplot(melt_short, aes(x= year, y=value, group=interaction(rep_id, max_disp), colour= sp_lpi.ID))+
  geom_line()+
  facet_grid(max_disp~ sp_lpi.ID)

ggplot(melt_lambda_short, aes(x= year, y=value, group=interaction(rep_id, max_disp), colour= sp_lpi.ID))+
  geom_line()+
  facet_grid(max_disp~ sp_lpi.ID)
## Warning: Removed 55440 rows containing missing values (geom_path).

# 
# melt_short<-melt_short[melt_short$max_disp == "113000" ,]
# 
# melt_lambda_short<-melt_lambda_short[melt_lambda_short$max_disp == "113000",]

melt_short<-melt_short[melt_short$max_disp == "disp113000" ,]

melt_lambda_short<-melt_lambda_short[melt_lambda_short$max_disp == "disp113000",]
library(taRifx)
library(plyr)
library(mgcv)

pops<-sp_lpi[,c(1,65:130)]
colnames(pops)[2:ncol(pops)]<-paste("Year", 1950:2015, sep="_")
pops[pops=="NULL"]<-NA
pops$rep_id<-"Observed"
pops$md_id<-"Observed"

popsm<-as.matrix(pops)

gam_lpi<-function(x){
   #subsetting the population data by each population 
  spid = x[2:(length(x)-2)]                     #subsetting only the dates
  names(spid)<-1950:2015              #renaming the date column names as R doesn't like numbered column names
  spid<-as.numeric(spid)
  pop_datab <- (!is.na(spid) )
  points = sum(pop_datab)
  id<-x[1]
  Date<-1950:2015
  spidt<-destring(t(spid))
  time<-length(min(which(!is.na(spidt))):max(which(!is.na(spidt))))
  missing<-time-points
  
  Year<-Date[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
  Population<-spidt[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
  Population[Population == 0] <- mean(Population, na.rm=TRUE)*0.01 #if a population is zero one year thhis is replaced with 1% of the average population estimate - because you can log zeros
  
  df<-data.frame(Year,Population)
  
  #not sure what this does - adding a constant of 1 so that logging doesn't go weird?
  if (sum(na.omit(df$Population<1))>0) {
    df$Population<-df$Population+1
  } 
    
  
  if (points >=6) {           
    PopN = df$Population
    if (length(na.omit(PopN)) >=6) {
      SmoothParm = round(length(na.omit(PopN))/2)    
    } else {
      SmoothParm=3
    }
    
    mg2<-mgcv:::gam(PopN ~ s(Year, k=SmoothParm), fx=TRUE)
    pv2 <- predict(mg2,df,type="response",se=TRUE) 
    R_sq2<-summary(mg2)$r.sq
    model<-1
    pv2$fit[pv2$fit <= 0] <- NA
    

    lambda2<-pv2$fit

    ial<-data.frame(id, Year,lambda2)
   
    colnames(ial)<-c("ID", "Year", "Abundance")
  }

  return(ial)
}

gam_lpi_r<-apply(popsm,  1, gam_lpi)
gam_r<-do.call( "rbind", gam_lpi_r)

fill<-data.frame(rep(pops$ID, each=length(1950:2015)), 1950:2015)
colnames(fill)<-c("ID", "Year")

all_year_ab<-join(fill, gam_r, type="right")

all_year_ab$max_disp<-"Observed"
all_year_ab$rep_id<-"Observed"

colnames(all_year_ab)[1:3]<-c("sp_lpi.ID", "year", "value")
mldab<-melt_short[,-4]
all_year_ab$sp_lpi.ID<-as.factor(all_year_ab$sp_lpi.ID)
all_year_ab$max_disp<-as.factor(all_year_ab$max_disp)
all_year_ab$rep_id<-as.factor(all_year_ab$rep_id)
all_year_ab$value<-as.numeric(all_year_ab$value)
all_year_ab$year<-as.numeric(all_year_ab$year)


both_df_ab<-rbind(mldab, all_year_ab)

both_df_ab[both_df_ab$sp_lpi.ID == "  539",]$sp_lpi.ID<-539
library(mgcv)

gam_demon<-function(id){
#for (i in 1:length(unique(melt_short$sp_lpi.ID))){
  
  #id<-unique(melt_short$sp_lpi.ID)[i]
  
  df<-melt_short[melt_short$sp_lpi.ID ==id,]
  PopN = df$value
  Year = df$year
  max_disp = df$max_disp
  rep_id = df$rep_id
  SmoothParm = round(length(na.omit(PopN))/2)    

  mg2<-mgcv:::gam(PopN ~ s(Year, k=15), fx=TRUE)
  pv2<-unique(fitted.values(mg2))
  
  ial<-data.frame(id, unique(Year),pv2)
   
  colnames(ial)<-c("sp_lpi.ID", "Year", "Abundance")
  #print(id)
  return(ial)
    }



gam_demon_r<-lapply(unique(melt_short$sp_lpi.ID), gam_demon)
gam_r<-do.call( "rbind", gam_demon_r)
library(mgcv)

gam_demon<-function(id){
  df<-melt_lambda_short[melt_lambda_short$sp_lpi.ID ==id,]
  PopN = df$value
  Year = df$year
  max_disp = df$max_disp
  rep_id = df$rep_id
  SmoothParm = round(length(na.omit(PopN))/2)    

  mg2<-mgcv:::gam(PopN ~ s(Year, k=15), fx=TRUE)
  pv2<-unique(fitted.values(mg2))
  
  ial<-data.frame(id, unique(Year),pv2)
   
  colnames(ial)<-c("sp_lpi.ID", "Year", "Abundance")
  return(ial)
    }

gam_demon_r_lambda<-lapply(unique(melt_lambda_short$sp_lpi.ID), gam_demon)
gam_r_lambda<-do.call( "rbind", gam_demon_r_lambda)
library(ggplot2)


ggplot(mldab, aes(x= year, y=value, group=interaction(rep_id, max_disp, sp_lpi.ID), colour= sp_lpi.ID))+
  geom_line(colour="grey")+
  geom_line(data=both_df_ab[both_df_ab$max_disp=="Observed",], aes(x=year, y=value), colour="red")+
  geom_line(data =  gam_r, aes(x =Year, y = Abundance, group=sp_lpi.ID), colour = "blue" )+
  facet_grid(.~ sp_lpi.ID)

pops<-sp_lpi[,c(1,65:130)]
colnames(pops)[2:ncol(pops)]<-paste("Year", 1950:2015, sep="_")
pops[pops=="NULL"]<-NA
pops$rep_id<-"Observed"
pops$md_id<-"Observed"

library(taRifx)
popsm<-as.matrix(pops)

gam_lpi<-function(x){
   #subsetting the population data by each population 
  spid = x[2:(length(x)-2)]                     #subsetting only the dates
  names(spid)<-1950:2015              #renaming the date column names as R doesn't like numbered column names
  spid<-as.numeric(spid)
  pop_datab <- (!is.na(spid) )
  points = sum(pop_datab)
  id<-x[1]
  Date<-1950:2015
  spidt<-destring(t(spid))
  time<-length(min(which(!is.na(spidt))):max(which(!is.na(spidt))))
  missing<-time-points
  
  Year<-Date[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
  Population<-spidt[min(which(!is.na(spidt))):max(which(!is.na(spidt)))]
  Population[Population == 0] <- mean(Population, na.rm=TRUE)*0.01 #if a population is zero one year thhis is replaced with 1% of the average population estimate - because you can log zeros
  
  df<-data.frame(Year,Population)
  
  #not sure what this does - adding a constant of 1 so that logging doesn't go weird?
  if (sum(na.omit(df$Population<1))>0) {
    df$Population<-df$Population+1
  } 
    
  
  if (points >=6) {           
    PopN = log10(df$Population)
    if (length(na.omit(PopN)) >=6) {
      SmoothParm = round(length(na.omit(PopN))/2)    
    } else {
      SmoothParm=3
    }
    
    mg2<-mgcv:::gam(PopN ~ s(Year, k=SmoothParm), fx=TRUE)
    pv2 <- predict(mg2,df,type="response",se=TRUE) 
    R_sq2<-summary(mg2)$r.sq
    model<-1
    pv2$fit[pv2$fit <= 0] <- NA
    

    lambda2<-diff(pv2$fit)

    ial<-data.frame(id, Year[-length(Year)], 10^lambda2)
   
    colnames(ial)<-c("ID", "Year", "R")
  }

  return(ial)
}

gam_lpi_r<-apply(popsm,  1, gam_lpi)
gam_r<-do.call( "rbind", gam_lpi_r)

fill<-data.frame(rep(pops$ID, each=length(1950:2015)), 1950:2015)
colnames(fill)<-c("ID", "Year")

all_year_r<-join(fill, gam_r, type="right")

all_year_r$max_disp<-"Observed"
all_year_r$rep_id<-"Observed"

colnames(all_year_r)[1:3]<-c("sp_lpi.ID", "year", "value")
mld<-melt_lambda_short[,-4]
all_year_r$sp_lpi.ID<-as.factor(all_year_r$sp_lpi.ID)
all_year_r$max_disp<-as.factor(all_year_r$max_disp)
all_year_r$rep_id<-as.factor(all_year_r$rep_id)
all_year_r$value<-as.numeric(all_year_r$value)
all_year_r$year<-as.numeric(all_year_r$year)


both_df<-rbind(mld, all_year_r)

need to add in x axis. do plots separately for each of the dispersal options?

library(ggplot2)

both_df<-both_df[both_df$max_disp != "1000",]

ggplot(both_df, aes(x= year, y=value, group=interaction(rep_id, max_disp,sp_lpi.ID), colour= sp_lpi.ID))+
  geom_line(colour="grey")+
  geom_line(data=both_df[both_df$max_disp=="Observed",], aes(x=year, y=value), colour="red")+
  geom_line(data =  gam_r_lambda, aes(x =Year, y = Abundance, group=sp_lpi.ID), colour = "blue" )+
  facet_grid(.~ sp_lpi.ID)

patch<-stack(paste(sdm_folder, "/pres_abs_sss_weighted_ensemble_sdm_", years,".tif", sep=""))
plot(patch[[s]])
patch<-projectRaster(patch, crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")


raster_read<-function(x){
b<-read.csv(paste(demoniche_folder ,"patch_old_disp113000/",x,"/Reference_matrix_pop_output.csv", sep="/"))
b<-b[,-1]
bef_ab_map<-rasterize(b[,c(1,2)],patch,b[,-c(1,2)])
return(bef_ab_map)
}

all_rep_stack<-lapply(1:12, raster_read)

big_stack<-stack(unlist(all_rep_stack))

year_seq<-seq(1, nlayers(big_stack),(nlayers(big_stack)/length(all_rep_stack)))


all_years_stack<-stack()
year_seq<-seq(1, nlayers(big_stack),(nlayers(big_stack)/length(all_rep_stack)))
for(i in 1:((nlayers(big_stack)/length(all_rep_stack)))-1){
  year_avg<-mean(big_stack[[year_seq]])
  all_years_stack<-stack(all_years_stack, year_avg)
  year_seq<-year_seq+1
  print(year_seq[1])
}

# 
# plot(all_years_stack[[1]])
# plot(all_years_stack[[100:167]])
# names(all_years_stack)<-1850:2016

writeRaster(all_years_stack, "Alpine_Ibex_Projection_disp_old.tif", overwrite=T)

ignore vvv

Using simulated habitat data

Basing the models on environmental data where the trend is known and constant - to better understand the impacts of the other variables.

plot(years_short,colMeans(na.omit(fake_map_mine)[4:length(colnames(fake_map_mine))]), type="l", ylab="Mean fake suitability index")

#setwd("C:Users/Fiona/Documents/Demoniche")

sd_seq<-c(0.5, 0.75,1, 1.25, 1.5)
#matrices_var<-matrix(0.45, ncol = 1, nrow = nrow(matrices), dimnames = list(NULL, "sd")) 
#k_seq<-c(1000,2000,4000,8000,16000,32000)
LDD_seq<-c(0.1, 0.5, 0.9)
kern_seq<-list(c(0.5,1,1,1),c(0.5,1,1,5), c(1,1,1,1), c(1,1,1,5))

var_grid<-expand.grid(sd_seq, LDD_seq, kern_seq)
colnames(var_grid)<-c("SD", "LDD", "Kern")

kern_mat<-matrix(c(rep(c(0.5,1,1,1), (nrow(var_grid)/length(unique(kern_seq)))),rep(c(0.5,1,1,5), (nrow(var_grid)/length(unique(kern_seq)))),rep(c(1,1,1,1),(nrow(var_grid)/length(unique(kern_seq)))),rep(c(1,1,1,5), (nrow(var_grid)/length(unique(kern_seq))))), ncol=4, byrow=T)


for (s in 1:nrow(var_grid)){

  print(paste (s, " out of ", nrow(var_grid) ), sep="")
  
  SD<-var_grid[s,1]
  LDD<-var_grid[s,2]
  
  
  matrices_var<-matrix(SD, ncol = 1, nrow = nrow(matrices), dimnames = list(NULL, "sd")) 

start.time <- Sys.time()

dir.create(paste("D:/Fiona/Git_Method/Git_Method/Demoniche_Repetitions/Fake_Three_Vars_kern/Kern_", s,"LDD_",LDD,"SD_", SD,sep=""))

reps<-5

rep_demoniche<-function(i){
  library(demoniche)
  source("demoniche_model_me.R")
 demoniche_setup(modelname = "Capra_k_16000",Populations = Populations, Nichemap = fake_spin_up,
                matrices = matrices,matrices_var = matrices_var, prob_scenario = prob_scenario,
                  stages = stages, proportion_initial = proportion_initial,
                  density_individuals = 400,  #number of individuals at the start of each population - can be a vector of length(Populations)
                  fraction_LDD = LDD, fraction_SDD = 0.5,
                  dispersal_constants = kern_mat[i,],
                  transition_affected_niche = "all",
                  transition_affected_demogr = transition_affected_demogr,
                  transition_affected_env=transition_affected_env,
                  env_stochas_type = env_stochas_type,
                  no_yrs = no_yrs_mine, K=2000, Kweight = K_weight, 
                  sumweight =sumweight)
c_ibex_k_16000 <- demoniche_model_me(modelname = "Capra_k_16000", Niche = TRUE, 
                                     Dispersal = TRUE, repetitions = 1,
                                     foldername = paste("Demoniche_Repetitions/Fake_Three_Vars_kern/Kern_", s,"LDD_",LDD,"SD_", SD,"/rep",i,sep=""))
}


#lapply(1:reps, rep_demoniche)

library(doParallel)

if (Sys.info()["nodename"] == "FIONA-PC"){
 cl <- makeCluster(4) 
} else {
 cl <- makeCluster(2)
}


registerDoParallel(cl)
foreach(i=(1:reps)) %dopar% rep_demoniche(i)
stopCluster(cl)

end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken 
}
  • Scaling population matrices to habitat suitability – only have one stage matrix from a particular point in time

  • Other matrices are just transition matrices but available in good and bad conditions, use these to great upper and lower bounds of survival? – need to check location, date and existence of these matrices

  • Noise is currently not supported but a value of 1 imputes white noise, and values higher or lower than this give positive or negative temporal autocorrelation – more important when you have multiple matrices?

  • Transition_affected_niche has more impact when specified values are inputted – if(is.numeric(transition_affected_niche){}

  • Transition affected niche pulls out parts of the matrix that are multiplied by the niche value – so if just presence/absence then it is either no change or 0

  • Potential problem with the way I have seeded the model – only ibex in the original populations whereas it might be more accurate to seed ibex more widely and then just sample from the original population sites.

  • Calculate dispersal kernel from gps data?

  • Density_individuals = is the starting populations for each locations that is “seeded” in Populations

  • Proportion_initial is how density individual is divided between stages, e.g all juveniles or evenly spread etc

  • Need to check if you can run multiple matrices_var in one run