Vingette update

vignettes/Marked_Point_Process.Rmd

@@ -60,7 +60,10 @@ euc <- euc[unlist(st_intersects(boundary, euc)),]
 
 class(trees)
 
-mesh = inla.mesh.2d(boundary = boundary, max.edge = 20)
+mesh = inla.mesh.2d(boundary = boundary, 
+                    max.edge = c(10, 20), 
+                    offset = c(15,20),
+                    cutoff = 2)
 mesh$crs <- st_crs(proj)
 
 ggplot() + 
@@ -74,77 +77,97 @@ ggplot() +
   ggtitle('Plot showing the food value index at each site')
 
 
-```
-
-```{r, analysis_of_data, eval = FALSE, echo= FALSE}
-
-data(euc) ##will add this in the future when data is on archive
-
 ```
 
 Lets first create a model for the tree locations only: in this example, we will assume that these locations are treated as *present-only* data.
 
 ```{r, points_only,fig.width=7, fig.height=5}
 
-points <- intModel(euc, Coordinates = c('x', 'y'), 
-                  Projection = proj, Mesh = mesh)
+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)
 
-plot(pointsPredictions)
+pointsPlot <- plot(pointsPredictions, variable = 'mean', 
+                   plot = FALSE)
+
+pointsPlot$predictions$predictions$mean + 
+  gg(euc)
+
 ```
 
 Now lets add the marks to this framework by specifying the name of the response variable of the marks with the argument `markNames`, and the family of the marks with the argument `markFamily`. Doing so will model each marks as a separate observation process, based on the family specified.
 
 ```{r, include_marks,fig.width=7, fig.height=5}
 
-marks <- intModel(euc, Coordinates = c('x', 'y'), Projection = proj,
-                  markNames = c('food', 'koala'), markFamily = c('gamma', 'poisson'),
-                  Mesh = mesh)
+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)
+                           mesh = mesh, marks = 'food', spatial = TRUE,
+                           fun = 'exp')
 
 koalaPredictions <- predict(marksModel, mask = boundary, 
                             mesh = mesh, marks = 'koala', spatial = TRUE)
 
-plot(foodPredictions)
-plot(koalaPredictions)
+plot(foodPredictions, variable = c('mean', 'sd'))
+plot(koalaPredictions, variable = c('mean', 'sd'))
 ```
 
 For the second mark model we only include the *food* mark, but this time we use a *log-linear* model with additive Gaussian noise. This is specified using the `.$updateFormula` slot function, and then by adding the scaling component to the *inlabru* components with the `.$addComponents` slot function to ensure that it is actually estimated. Moreover we assume *penalizing complexity* priors for the two spatial effects, as well as specify *bru_max_iter* in the *options* argument to keep the time to estimate down.
 
 ```{r, marks_add_scaling,fig.width=7, fig.height=5}
 
-marks2 <- intModel(euc, Coordinates = c('x', 'y'), Projection = proj,
-                  markNames = 'food', markFamily = 'gaussian',
-                  Mesh = mesh, pointsSpatial = 'individual')
+marks2 <- startMarks(euc, Boundary = boundary,
+                     Projection = proj,
+                     markNames = 'food', 
+                     markFamily = 'gaussian',
+                     Mesh = mesh)
 
-marks2$updateFormula(markName = 'food', 
-      newFormula = ~ exp(food_intercept + (euc_spatial + 1e-6)*scaling + food_spatial))
+marks2$updateFormula(Mark = 'food', 
+      newFormula = ~ exp(food_intercept + (shared_spatial + 1e-6)*scaling + food_spatial))
 
-marks2$changeComponents(addComponent = 'scaling')
+marks2$changeComponents(addComponent = 'scaling(1)')
 
-marks2$specifySpatial(datasetName = 'euc', 
-                      prior.sigma = c(0.1, 0.01), 
-                      prior.range = c(10, 0.01))
+marks2$specifySpatial(sharedSpatial = TRUE, 
+                      prior.sigma = c(1, 0.01), 
+                      prior.range = c(120, 0.01))
 
 marks2$specifySpatial(Mark = 'food',  
-                      prior.sigma = c(0.1, 0.01), 
-                      prior.range = c(10, 0.01))
+                      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 + (euc_spatial + 1e-6)*scaling + food_spatial))
+    formula =  ~ (food_intercept + (shared_spatial + 1e-6)*scaling + food_spatial))
 
 plot(predsMarks2)
 ```

vignettes/Spatiotemporal_example.Rmd

@@ -39,20 +39,18 @@ as well as define some objects required by the model to run.
 
 ```{r, Alabama_map}
 
-proj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
+proj <- "+proj=utm +zone=17 +datum=WGS84 +units=km"#"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
 
 AL <- USAboundaries::us_states(states = "Alabama", resolution = 'high')
 AL <- as(AL, "sf")
-st_crs(AL) <- proj
+AL <- st_transform(AL, proj)
 
 mesh <- inla.mesh.2d(boundary = inla.sp2segment(AL[1]), 
-                     cutoff = 0.1,
-                     max.edge = c(0.2, 0.8), 
-                     offset = c(0.1, 0.2),
+                     cutoff = 0.1 * 50,
+                     max.edge = c(0.2, 0.8) * 70, 
+                     offset = c(0.1, 0.2) * 150,
                      crs = st_crs(proj))
 
-
-
 ```
 
 The first dataset we consider is obtained from the North American Breeding Bird Survey. This dataset may be loaded directly from the package, and contains observations of the species between 2015 and 2017. This dataset is treated as replicate present-absent, where every point is assumed to be a visited site (or *route*).
@@ -71,21 +69,24 @@ eBird2015 <- spocc::occ(
   query = 'Colinus virginianus',
   from = 'gbif',
   date = c("2015-01-01", "2015-12-31"),
-  geometry = AL
+   geometry = st_bbox(st_transform(AL,
+                '+proj=longlat +datum=WGS84 +no_defs'))
 )$gbif
 
 eBird2016 <- spocc::occ(
   query = 'Colinus virginianus',
   from = 'gbif',
   date = c("2016-01-01", "2016-12-31"),
-  geometry = AL
+ geometry = st_bbox(st_transform(AL,
+                '+proj=longlat +datum=WGS84 +no_defs'))
 )$gbif
 
 eBird2017 <- spocc::occ(
   query = 'Colinus virginianus',
   from = 'gbif',
   date = c("2017-01-01", "2017-12-31"),
-  geometry = AL
+ geometry = st_bbox(st_transform(AL,
+                '+proj=longlat +datum=WGS84 +no_defs'))
 )$gbif
 
 eBird <- data.frame(eBird2015$data[[1]]) %>% 
@@ -95,28 +96,35 @@ eBird <- data.frame(eBird2015$data[[1]]) %>%
 
 eBird <- st_as_sf(x = eBird,
                   coords = c('longitude', 'latitude'),
-                  crs = proj)
+                  crs =  '+proj=longlat +datum=WGS84 +no_defs')
 
 eBird$Year <- eBird$year
 
+eBird <- st_transform(eBird, proj)
+
 eBird <- eBird[unlist(st_intersects(AL, eBird)),]
 
 ```
 
-We then get onto the model description, which in this case includes a shared spatial field between the two datasets. This shared spatial field is characterized by an *ar1* process. To add this structure into the model, we specify the parameter *temporalModel* in the function `intModel` appropriately, Furthermore we specified the hyper parameters for both the random field and the temporal effect.
+We then get onto the model description, which in this case includes a shared spatial field between the two datasets. This shared spatial field is characterized by an *ar1* process. To add this structure into the model, we specify the parameter *temporalModel* in the function `startISDM` appropriately, Furthermore we specified the hyper parameters for both the random field and the temporal effect and the priors for the intercepts.
 
 ```{r, setup_model}
 
-hyperParams <- list(rho = list(prior = "pc.prec", param = c(0.01, 0.01)))
+hyperParams <- list(model = 'ar1', 
+                    hyper = list(rho = list(prior = "pc.prec", param = c(0.9, 0.1))))
 
-modelSetup <- intModel(eBird, BBSColinusVirginianus,
-                       Coordinates = c('Longitude', 'Latitude'),                                                    temporalName = 'Year',
+modelSetup <- startISDM(eBird, BBSColinusVirginianus,
+                       temporalName = 'Year',
+                       Boundary = AL,
                        Projection =  proj, Mesh = mesh, 
-                       responsePA =  'NPres', trialsPA = 'Ntrials',
-                       temporalModel = list(model = 'ar1', hyper = hyperParams))
+                       responsePA =  'NPres', trialsPA = 'Ntrials')
+
+modelSetup$specifySpatial(sharedSpatial = TRUE, prior.sigma = c(1, 0.5), 
+                          prior.range = c(100, 0.5))
+
+modelSetup$specifyRandom(temporalModel = hyperParams)
 
-modelSetup$specifySpatial(sharedSpatial = TRUE, prior.sigma = c(0.2, 0.01), 
-                          prior.range = c(0.4, 0.01))
+modelSetup$priorsFixed(Effect = 'intercept', mean.linear = 0, prec.linear = 0.001)
 
 ```
 
@@ -141,7 +149,8 @@ Next we run the model, using the function `fitISDM`. Due to time considerations,
 ```{r, run_model}
 
 mod <- fitISDM(modelSetup,
-                options = list(control.inla = list(int.strategy = 'eb')))
+                options = list(control.inla = list(int.strategy = 'eb',
+                                                   diagonal = 1)))
  
 ```