## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE,
  warning = FALSE,
  message = FALSE
)


## ----setup--------------------------------------------------------------------
#  library(spatstat)
#  library(PointedSDMs)
#  library(sf)
#  library(sp)
#  library(ggplot2)
#  library(inlabru)
#  library(INLA)

## ----load_data----------------------------------------------------------------
#  
#  data(Koala)
#  eucTrees <- Koala$eucTrees
#  boundary <- Koala$boundary
#  

## ----clean_data, echo = FALSE,fig.width=7, fig.height=5-----------------------
#  
#  proj <- "+init=epsg:27700"
#  
#  boundary <- as(boundary, 'sf')
#  st_crs(boundary) <- proj
#  
#  euc <- st_as_sf(x = eucTrees,
#                  coords = c('E', 'N'),
#                  crs = proj)
#  
#  euc$food <- euc$FOOD/1000
#  euc <- euc[unlist(st_intersects(boundary, euc)),]
#  
#  class(trees)
#  
#  mesh = inla.mesh.2d(boundary = boundary,
#                      max.edge = c(10, 20),
#                      offset = c(15,20),
#                      cutoff = 2)
#  mesh$crs <- st_crs(proj)
#  
#  ggplot() +
#    geom_sf(data = st_boundary(boundary)) +
#    geom_sf(data = euc, aes(color = koala)) +
#    ggtitle('Plot showing number of koalas at each site')
#  
#  ggplot() +
#    geom_sf(data = st_boundary(boundary)) +
#    geom_sf(data = euc, aes(color = food)) +
#    ggtitle('Plot showing the food value index at each site')
#  
#  

## ----points_only,fig.width=7, fig.height=5------------------------------------
#  
#  points <- startISDM(euc, Boundary = boundary,
#                      Projection = proj,
#                      Mesh = mesh)
#  
#  points$specifySpatial(sharedSpatial = TRUE,
#                        prior.range = c(120, 0.1),
#                        prior.sigma = c(1, 0.1))
#  
#  pointsModel <- fitISDM(points, options = list(control.inla = list(int.strategy = 'eb')))
#  
#  pointsPredictions <- predict(pointsModel, mask = boundary,
#                               mesh = mesh, predictor = TRUE)
#  
#  pointsPlot <- plot(pointsPredictions, variable = 'mean',
#                     plot = FALSE)
#  
#  pointsPlot$predictions$predictions$mean +
#    gg(euc)
#  

## ----include_marks,fig.width=7, fig.height=5----------------------------------
#  
#  marks <- startMarks(euc, Boundary = boundary,
#                      Projection = proj,
#                      markNames = c('food', 'koala'),
#                      markFamily = c('gamma', 'poisson'),
#                      Mesh = mesh)
#  
#  marks$specifySpatial(sharedSpatial = TRUE,
#                        prior.range = c(120, 0.1),
#                        prior.sigma = c(1, 0.1))
#  
#  marks$specifySpatial(Mark = 'koala',
#                       prior.range = c(120, 0.1),
#                       prior.sigma = c(1, 0.1))
#  
#  marks$specifySpatial(Mark = 'food',
#                       prior.range = c(60, 0.1),
#                       prior.sigma = c(1, 0.1))
#  
#  marksModel <- fitISDM(marks, options = list(control.inla = list(int.strategy = 'eb'),
#                                               safe = TRUE))
#  
#  foodPredictions <- predict(marksModel, mask = boundary,
#                             mesh = mesh, marks = 'food', spatial = TRUE,
#                             fun = 'exp')
#  
#  koalaPredictions <- predict(marksModel, mask = boundary,
#                              mesh = mesh, marks = 'koala', spatial = TRUE)
#  
#  plot(foodPredictions, variable = c('mean', 'sd'))
#  plot(koalaPredictions, variable = c('mean', 'sd'))

## ----marks_add_scaling,fig.width=7, fig.height=5------------------------------
#  
#  marks2 <- startMarks(euc, Boundary = boundary,
#                       Projection = proj,
#                       markNames = 'food',
#                       markFamily = 'gaussian',
#                       Mesh = mesh)
#  
#  marks2$updateFormula(Mark = 'food',
#        newFormula = ~ exp(food_intercept + (shared_spatial + 1e-6)*scaling + food_spatial))
#  
#  marks2$changeComponents(addComponent = 'scaling(1)')
#  
#  marks2$specifySpatial(sharedSpatial = TRUE,
#                        prior.sigma = c(1, 0.01),
#                        prior.range = c(120, 0.01))
#  
#  marks2$specifySpatial(Mark = 'food',
#                        prior.sigma = c(1, 0.01),
#                        prior.range = c(120, 0.01))
#  
#  marksModel2 <- fitISDM(marks2, options = list(control.inla = list(int.strategy = 'eb'),
#                                                bru_max_iter = 2, safe = TRUE))
#  
#  predsMarks2 <- predict(marksModel2, mask = boundary, mesh = mesh,
#      formula =  ~ (food_intercept + (shared_spatial + 1e-6)*scaling + food_spatial))
#  
#  plot(predsMarks2)