sims_workflow.Rmd

#---
title: "sims_workflow"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(boot)
library(ggplot2)
library(plyr)

```

Get parameter file and time stamp

```{r}

rm(list=ls())

parameterDF <- read.csv("parameters_moreRangeChangeSims.csv")
#parameteDF is a file with the list of the all the parameters and their values for each simulation run

#decie on which run
type = 22

systemTime <- gsub(" ","_", Sys.time(),fixed=TRUE)
systemTime <- gsub(":","-",systemTime)

```

Load in functions

```{r}

source('R/getParams.R')
source('R/getOccuHistory.R')
source('R/getRepeatVisits.R')
source('R/getSpartaFormat.R')
source('R/getBugsData.R')
source('R/runModel.R')
source('R/helperFunctions.R')

```


Create some replicated worlds

```{r, echo=FALSE, results='hide'}

getParams(parameterDF,mysimNu=type)
Obs <- replicate(300,getSims(),simplify=FALSE)
save(Obs,file=paste0("output/Obs_",type,"_",systemTime,".RData"))
getObsSummaries(Obs)
trueSummaries <- gettrueSummaries(combineObs(Obs))

```

Run the basic model

```{r, message=FALSE, warning=FALSE}

out<- lapply(Obs,function(x)getModels(x,type=type,myModel))
#save(out,file=paste0("output/out_",type,"_",systemTime,".RData"))

```

Compare real and predicted intercept (psi1) and slope (psi change from last to first)

```{r}

checkResults(out,Obs)

```

Get summaries over parameters, model types

```{r}

print(modelSummaries <- getModelSummaries(out))

```

Compare to the analysis of raw data

```{r, message=FALSE, warning=FALSE}

outRaw <- ldply(Obs,function(x)rawdataAnalysis(x))
#save(outRaw,file=paste0("output/outRaw_",type,"_",systemTime,".RData"))
rawSummaries <- getrawSummaries(outRaw)

```

plot comparison
```{r}

plotComparison(modelSummaries,rawSummaries,trueSummaries)
#500/250

```

plot comparison of differences

```{r}
plotDifference(Obs,out,outRaw)

```

Power analysis

```{r}

powerAnalysis(out,outRaw)

```

Later comparison

```{r}

type <- 19

#load outputs
files <- list.files("output/")
setwd("output")
load(files[grepl(paste("Obs",type,sep="_"),files)])
load(files[grepl(paste("out",type,sep="_"),files)])
load(files[grepl(paste("outRaw",type,sep="_"),files)])

#summarise
modelSummaries <- getModelSummaries(out)
rawSummaries <- getrawSummaries(outRaw)
trueSummaries <- gettrueSummaries(combineObs(Obs))

#plot
plotComparison(modelSummaries,rawSummaries,trueSummaries)
plotDifference(Obs,out,outRaw)

```


(1) Autocorrelation

(a) probability of dependence between sequential visits
- at a proportion of sites
- autocorrelation of order p with certain probability

#fix = indicate previous observation as an indicator in the obs model

```{r}

source('R/scenario_autocorrelation.R')

runSims <- function(Obs,propSites=0.5,autoProb=0.5){
  #reduce visits in non-atlas years
  Obs_auto <-applyAutocorrelation(Obs,propSites,autoProb)
  return(Obs_auto)
}
Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="Autocorrelation",myModel))

#fixed OD model
outFixed <- lapply(Obs_auto,function(x)fixAutocorrelation(x))

```

(b) Not reporting after first observation

- at prop p of sites

# fix = remove records after first detection??

```{r}

source('R/scenario_autocorrelation.R')

runSims <- function(Obs,propSites=0.5){
  #reduce visits in non-atlas years
  Obs_auto <-applyFullAutocorrelation(Obs,propSites)
  return(Obs_auto)
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="FullAutocorrelation",myModel))

#fixed OD model
outFixed <- lapply(Obs_auto,function(x)fixFullAutocorrelation(x))

```

(2) Atlas scenarios

(a) Pulse in recording - increased visits in one year (atlas year) 

#fix might not be needed

```{r}

source('R/scenario_pulsedvisits.R')

runSims <- function(Obs,nStandardVisits=3){
  #reduce visits in non-atlas years
  Obs_pulsed <- pulsedVisits(Obs, nStandardVisits)
  #return summary data
  return(Obs_pulsed)
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw <- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit OD models
out <- lapply(Obs_auto,function(x)getModels(x,bias="Pulsed",myModel))

```

(b) Pulse in recording - fewer sites are visited outside of atlas sites

#as long as unbiased should be ok

```{r}

source('R/scenario_atlassites.R')

runSims <- function(Obs,nStandardSites=10){
  #reduce visits in non-atlas years
  Obs_spread <- reduceSites(Obs,nStandardSites)
  #return summary data
  return(Obs_spread)
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))
#fit OD models
out <- lapply(Obs_auto,function(x)getModels(x,bias="Spread",myModel))

```

(c) Pulse in recording - Lower quality grid cells less likely to be visited outside of atlas years

```{r}

source('R/scenario_pulsedvisits.R')

runSims <- function(Obs,lowProb=0.25){
  #reduce visits in non-atlas years
  Obs_spread <- pulsedSpreadVisits(Obs,SCovariate,lowProb)
  #return summary data
  return(Obs_spread)
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit OD models
out <- lapply(Obs_auto,function(x)getModels(x,bias="SpreadLower",myModel))

```

(2) Effects of habitat degradation on effort ("car park")

(a) fewer visits to low quality sites (site quality static)

```{r}

source('R/scenario_frequencyvisits.R')

runSims <- function(Obs,type="static",siteEffect=0.15){
  Obs_reduced <-applyReducedVisits(Obs,type=type,siteEffect=siteEffect)
  return(Obs_reduced)
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="ReducedQ",myModel))

```

(b) lower detection probability in lower quality sites (site quality static)

```{r}

source('R/scenario_lowquality.R')

runSims <- function(Obs,meanSiteDetEffect=0.15){
  
  #specify site effect on detection probability
  SiteDetEffects <- rep(meanSiteDetEffect,NSpecies)
  
  #get new repeat visits
  Obs_lq <- getLQSims(Obs)
  
  return(Obs_lq)
  
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="ReducedD",myModel))

```

(c) as(a) except with site quality of some sites declining (positive interaction term)
fewer visits over time in low quality sites

```{r}

source('R/scenario_frequencyvisits.R')

runSims <- function(Obs,type="dynamic",timeEffect=0.1){
  
  #reduce visits in non-atlas years
  Obs_reduced <-applyReducedVisits(Obs,type=type,timeEffect=timeEffect)
  
  return(Obs_reduced)
  
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="ReducedQD",myModel))

```

(d) as(b) except with site quality of some sites declining (positive interaction term)
detection probability declines over time in low quality sites

```{r}

source('R/scenario_lowquality.R')


runSims <- function(Obs,timeEffect=0.1){
  
  #low quality sites are declining
  IntDetEffects <<- rep(timeEffect,NSpecies)

  #get new repeat visits
  Obs_lq <- getLQSims(Obs)
  
  return(Obs_lq)
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="ReducedDD",myModel))

```

(3) CS observer bias

(a) Bias towards specialists

- specialist trait covary with abundance (OccProb or Lambda)
- assume this affects species' detection probability 

```{r}

source('R/scenario_specialistbias.R')

runSims <- function(Obs,propSites=1){

  Obs <- getbiasedSims(propSites)
  
  return(Obs)
  
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="SpecialistSp",myModel))


```

(b) Bias towards declining species

- assume some species are declining
- assume this affects species' detection probability 
- at prop p of sites

```{r}

source('R/scenario_specialistbias.R')

runSims <- function(Obs,biasEffect=0.15,propSites=1){
  
  Obs <- getbiaseddeclineSims(biasEffect,propSites)
  
  #get new repeat visits
  Obs_lq <- getLQSims(Obs)
  
  return(Obs_lq)
  
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out <- lapply(Obs_auto,function(x)getModels(x,bias="DeclineSp",myModel))

```

(4) Probability of site visitation is higher if species previously seen there

#not work

```{r}

source('R/scenario_frequencyvisits.R')

runSims <- function(Obs){
  
  #reduce visits in non-atlas years
  Obs_reduced <-applyBiasedVisits(Obs)
  
  return(Obs_reduced)
  
}

Obs_auto <- lapply(Obs,runSims)

#raw data analysis
outRaw<- ldply(Obs_auto,function(x)rawdataAnalysis(x))

#fit standard OD model
out_reduced <- lapply(Obs_auto,function(x)getModels(x,bias="ReducedV",myModel))

```

<!-- (5) Mismatch between scale of dynamics and observation -->

<!-- ```{r} -->
<!-- source('C:/Users/db40fysa/Nextcloud/sMon-Analyses/MoreRangeChangeSims/R/scenario_spatialmismatch.R') -->

<!-- runSims <- function(x){ -->
<!--    Occ <- getOccuHistory(OccProb, spatialbeta=Seffects, interactionbeta=Ieffects, spatialcovariates=Scovariate,temporalbeta=Teffects, temporalcovariates=Tcovariate) -->
<!--   Obs <- getRepeatVisits(Occ, NVisits, DetProb, SiteDetEffects,YearDetEffects,IntDetEffects, Scovariate,Tcovariate) -->

<!--   #apply scenario -->
<!--  Obs_grouped <- groupSite(NSites=NSites, Obs=Obs) -->
<!--  Obs_reduced <- reduceVisits(Obs=Obs_grouped) -->

<!--   #with original dataset -->
<!--   modelData <- getSpartaFormat(Obs,focalSpecies=focalSpecies) -->
<!--   bugs.data <- getBugsData(modelData) -->
<!--   modelSummary <- runModel(bugs.data,modelData,focalSpecies) -->
<!--   modelSummary$Type <- "Standard" -->

<!--   #with fixed dataset -->
<!--   #modelData_grouped <- getSpartaFormat(Obs=Obs_grouped,focalSpecies=focalSpecies) -->
<!--   #bugs.data_grouped <- getBugsData(modelData=modelData_grouped) -->
<!--   #modelSummary_grouped <- runModel(bugs.data=bugs.data_grouped,modelData=modelData_grouped,focalSpecies) -->
<!--   #modelSummary_grouped$Type <- "grouped" -->

<!--   #with fixed dataset -->
<!--   modelData_reduced <- getSpartaFormat(Obs=Obs_reduced,focalSpecies=focalSpecies) -->
<!--   bugs.data_reduced <- getBugsData(modelData=modelData_reduced) -->
<!--   modelSummary_reduced <- runModel(bugs.data=bugs.data_reduced, -->
<!--                                    modelData=modelData_reduced, -->
<!--                                    focalSpecies,myModel) -->
<!--   modelSummary_reduced$Type <- "Reduced" -->

<!--   return(rbind(modelSummary,modelSummary_reduced)) -->

<!-- } -->
<!-- ``` -->