Update comments re: quadrature schemes

RavRK_0_prepare.R

@@ -43,7 +43,7 @@ createRoute <- function(route.site, dilation = 1000){
 #------
 
 # route and habitat data
-route_data <- read_csv("data/route_data_2021_02_02.csv")
+route_data <- read_csv("data/route_data_2021_02_03.csv")
 routes <- route_data %>% split(route_data$route)
 linnets <- routes %>% lapply(createRoute)
 
@@ -84,7 +84,7 @@ data <- lapply(seasons, function(seas){
 })
 
 # use spatstat::lppm to fit (full) saturated model
-# this generates a quadrat scheme and corresponding model frame for each route-season pairs (32) 
+# this generates a quadrature scheme and corresponding model frame for each route-season pairs (32) 
 # dummy spacing eps = 15 is equivalent to the default scheme, but has been set explicitly for reproducibility. 
 lppm_fits <- lapply(seasons, function(seas) with(data[[seas]], try(lppm(ravens.lpp, ~ O.dist + RK.dist + 1,eps =15))))
 model_frames <- lapply(lppm_fits, lapply, function(fit) fit %>% model.frame.lppm %>% rename(y = .mpl.Y, w = `(weights)`))

RavRK_1_fits.R

@@ -29,12 +29,12 @@ rm(list=ls())
 #-----
 # Data
 #-----
-# ravens data are the combined model matrices from (the default) linear quadrat schemes generated 
-#   for each route using spatstat's lppm (see RavRK_0_prepare.R)
-# Each row is a single quadrat in a given route in a given season 
-# y2: 0 or 1, raven not-observed/observed in quadrat
-# w: weights = 'area' of quadrat (proportional to quadrat length)
-# y1 = y2/w (intensity/rate = observed raven/metre in quadrat)
+# ravens data are the combined model matrices from the (default) linear quadrature schemes generated 
+#   for each route using spatstat's lppm (see RavRK_0_prepare.R).
+# Each row is a single quadrature point in a given route in a given season. 
+# y2: 0 or 1, raven not-observed/observed at the quadrature point.
+# w: weights = 'length' of route associated with the quadrature point
+# y1 = y2/w (intensity (or rate) = observed raven/metre in the length associated with the quadrature point)
 # O.dist (numerical covariate): distance to open farmland from quadrat (X_F)
 # RK.dist (numerical covariate): distance to roadkill from quadrat (X_R)
 # route (group factor): 8 routes where observations were made
@@ -50,7 +50,7 @@ ravensScaled <- ravens %>% mutate(O.dist = scale(O.dist), RK.dist = scale(RK.dis
 #--------------------
 
 # Point process Poisson models can be fit using standard glms, with Poisson errors and canonical log link
-# The quadrat 'sizes' (effort), w,  must incorporated into the fit. 
+# The quadrature 'weights', w,  must incorporated into the fit. 
 # Two methods:
 # 1) Use w as prior weights with response y1 = 1/w or 0 for presence/absence (non-integer counts may throw warnings but estimates are correct)
 # 2) Use log(w) as an offset with response y2 = 1 or 0 presence/absence [we use this]

RavRK_example.Rmd

@@ -47,7 +47,7 @@ m.spat <- lppm(ravens_lpp ~ 1 + X_F + X_R, eps = 15)
 m.spat %>% model.frame.lppm %>% as_tibble
 ```
 
-This generates a numerical quadrature scheme with added 'dummy points', each spaced by eps = 15 units. Each row of the model frame is a single quadrat, where the weight (`(weights)`) denotes the quadrat size (effort) and the response (`.mpl.Y`) is the observed density (0 for a dummy point, 1/weight for an observation). The covariate values are given in the remaining two columns. 
+This generates a numerical quadrature scheme which includes a point for each observation together with added 'dummy points', each spaced by eps = 15 units. The role of the quadrature scheme is to enable numerical evaluation of the model likelihood which involves (numerical) integration of the density function along the length of each route. Within the model frame, each row corresponds to a quadrature point, where the weight (`(weights)`) denotes the length of route associated with the point and the response (`.mpl.Y`) is the density (0 for a dummy point, 1/weight for an observation). The covariate values are given in the remaining two columns. 
 
 We rename the columns and add a presence/absence transformation of the response which we call `y2` (0 for a dummy point, 1 for an observation). 
 ```{r}
@@ -122,7 +122,7 @@ loo_c <- loo_compare(loo(m.stan.0, cores = 10),loo(m.stan.1, cores = 10))
 # manually stored results since loo is slow to run
 matrix(c(0.0,0.0,-6.1,5.4),2,2, byrow = T, dimnames = list(c("m.stan.1","m.stan.0"),c("elpd_diff","se_diff"))) 
 ```
-Leave-one-out cross validation, here based on predictive log densities and approximated using importance-sampling methods, is an estimate of the information-theoretic quantity Kullback-Leibler (KL) discrepancy. In this way it is similar to AIC, which is also an estimate of KL discrepancy, although cross validation is a more flexible estimator requiring less assumptions than AIC. Importantly, cross validation can be applied to compare hierarchical models with differing effect structures, providing the full set of group-level parameters are conditioned on to ensure the conditional independence of the test data points. 
+Leave-one-out cross validation, here based on predictive log densities and approximated using Pareto-smoothed importance-sampling techniques, is an estimate of the information-theoretic quantity Kullback-Leibler (KL) discrepancy. In this way, it is similar to AIC (up to a factor of $-2$), which is also an estimate of KL discrepancy, although cross validation is a more flexible estimator requiring less assumptions than AIC. Importantly, cross validation can be applied to compare hierarchical models with differing effect structures, providing the full set of group-level parameters are conditioned on to ensure the conditional independence of the test data points. 
 
 ---
 

RavRK_example_cache/html/__packages

@@ -0,0 +1,19 @@
+base
+tidyverse
+ggplot2
+tibble
+tidyr
+readr
+purrr
+dplyr
+stringr
+forcats
+spatstat.data
+spatstat.geom
+nlme
+rpart
+spatstat.core
+spatstat.linnet
+spatstat
+Rcpp
+rstanarm