Simulated_Data_Ex.Rmd

---
title: "ClusterExcercise"
output: pdf_document
date: "2023-06-20"
---

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

## Group Size, Detection Probability Bias

Package Prep and Library

```{r packages, warning=FALSE}
library(tidyverse)
library(Distance)
library(knitr)
library(mrds)
#install.packages("devtools")
 library(devtools)
  install_github("DistanceDevelopment/dssd", build_vignettes = TRUE)
  install_github("DistanceDevelopment/dsims", build_vignettes = TRUE)
```

Data import and formatting

```{r import}
library(Distance)
data("sikadeer")

```

Working with Dsims

```{r}
library(dsims)

# Find the file path to the example shapefile in dssd
shapefile.name <- system.file("extdata", "StAndrew.shp", package = "dssd")
# Create the survey region object
region <- make.region(region.name = "St Andrews bay",
                      shape = shapefile.name,
                      units = "m")
plot(region)
```

```{r}
# We first create a flat density grid
density <- make.density(region = region,
                        x.space = 500,
                        constant = 1)

# Now we can add some high and low points to give some spatial variability
density <- add.hotspot(object = density,
                       centre = c(-170000, 6255000),
                       sigma = 8000,
                       amplitude = 4)

density <- add.hotspot(object = density,
                       centre = c(-160000, 6275000),
                       sigma = 6000,
                       amplitude = 4)

density <- add.hotspot(object = density,
                       centre = c(-155000, 6260000),
                       sigma = 3000,
                       amplitude = 2)

density <- add.hotspot(object = density,
                       centre = c(-150000, 6240000),
                       sigma = 10000,
                       amplitude = -0.9)

density <- add.hotspot(object = density,
                       centre = c(-155000, 6285000),
                       sigma = 10000,
                       amplitude = -1)

# I will choose to plot in km rather than m (scale = 0.001)
plot(density, region, scale = 0.001)
```

```{r}
density.data <- density@density.surface[[1]]
head(density.data)
```

```{r}
library(mgcv)
fit.gam <- gam(density ~ s(x,y), data = density.data, family = gaussian(link="log"))

# Use the gam object to create a density object
gam.density <- make.density(region = region,
                            x.space = 500,
                            fitted.model = fit.gam)

plot(gam.density, region, scale = 0.001)
```

```{r}
# Create a covariate list describing the distribution of cluster sizes
covariates <- list(size = list(distribution = "ztruncpois", mean = 3))

# Define the population description
pop.desc <- make.population.description(region = region,
                                        density = gam.density,
                                        covariates = covariates,
                                        N = 300,
                                        fixed.N = TRUE)
```

```{r}
parallel.design <- make.design(region = region, 
                               design = "systematic",
                               spacing = 2500,
                               edge.protocol = "minus",
                               design.angle = 90,
                               truncation = 750)

zigzag.design <- make.design(region = region, 
                             design = "eszigzag",
                             spacing = 2233,
                             edge.protocol = "minus",
                             design.angle = 0,
                             bounding.shape = "convex.hull",
                             truncation = 750)

p.survey <- generate.transects(parallel.design)
plot(region, p.survey)

z.survey <- generate.transects(zigzag.design)
plot(region, z.survey)
```

```{r}
# Define the covariate parameters on the log scale
cov.param <- list(size = log(1.08))

# Create the detectability description
detect <- make.detectability(key.function = "hn",
                             scale.param = 300,
                             cov.param = cov.param,
                             truncation = 750)

# Plot the simulation detection functions
plot(detect, pop.desc)
```

```{r}
analyses <- make.ds.analysis(dfmodel = list(~1, ~1, ~size),
                             key = c("hn", "hr", "hn"),
                             truncation = 750,
                             er.var = "R2",
                             criteria = "AIC")
```

```{r}
sim.parallel <- make.simulation(reps = 5,
                                design = parallel.design,
                                population.description = pop.desc,
                                detectability = detect,
                                ds.analysis = analyses)

sim.zigzag <- make.simulation(reps = 5,
                              design = zigzag.design,
                              population.description = pop.desc,
                              detectability = detect,
                              ds.analysis = analyses)
```

```{r}
# Generate a single instance of a survey: a population, set of transects 
# and the resulting distance data
eg.parallel.survey <- run.survey(sim.parallel)

# Plot it to view a summary
plot(eg.parallel.survey, region)
```

```{r}
# Generate a single instance of a survey: a population, set of transects 
# and the resulting distance data
eg.zigzag.survey <- run.survey(sim.zigzag)

# Plot it to view a summary
plot(eg.zigzag.survey, region)
```

```{r}
# Running the simulations
sim.parallel <- run.simulation(sim.parallel)
sim.zigzag <- run.simulation(sim.zigzag)
```

```{r}
summary(sim.zigzag)

oldparams <- par(mfrow = c(1,2))
histogram.N.ests(sim.parallel)
histogram.N.ests(sim.zigzag)

par(oldparams)
```