Tyler D. Rudolph
Estimation of the weighted and un-weighted bivariate kernel probability density function for spatially explicit1 population-level processes given repeated measures of individual sample units. May be used to estimate the utilization distribution of a wildlife or fish population monitored via satellite (e.g. GPS collar transmitters) or acoustic (e.g. PTT tags & receptor stations) data telemetry while accounting for individual variation in point sample size, dispersion and quality, thereby generating unbiased estimates of population distribution.
Two weighting strata are possible:
-
Individual sample observations (e.g. GPS relocations) may be weighted to reflect variation in positional accuracy (e.g. PDOP, HDOP) and/or datum quality;
-
Individual kernel distributions may be weighted to reflect variation in input data sufficiency, representativity and/or reliability.
if(!require('devtools')) install.packages('devtools')
devtools::install_github('https://github.com/tylerdrudolph/wmKDE')
library(wmKDE)
if(!require('dplyr')) install.packages('dplyr')
library(dplyr)
if(!require('sf')) install.packages('sf')
library(sf)
if(!require('terra')) install.packages('terra')
library(terra)
options(warn = -1)Case 1) Two individual point patterns with different centres of activity but identical dispersion (sd = 1) and sample size (n=500)
locs <- bind_rows(data.frame(id = 1, col=8, x = rnorm(n=500), y = rnorm(n=500)),
data.frame(id = 2, col=4, x = rnorm(n=500, mean=-3.5), y = rnorm(n=500, mean=-3.5))) %>%
st_as_sf(., coords=c('x','y'), crs=st_crs(6624))
## basic kernel (relocations pooled)
bkern1 <- wmKDE(locs, avg=F, spatres=0.1, )
## mean kernel utilization distribution (identical result)
mkern1 <- wmKDE(locs, id='id', ncores=1, avg=T, spatres=0.1)## plot side by side
par(mfrow=c(1,2))
plot(bkern1$wmKDE$iso, main='Pooled Kernel \n n(ID1) = 500; n(ID2) = 500')
plot(st_geometry(locs), pch=20, cex = 0.6, col=locs$col, add=T)
plot(mkern1$isocontours %>% st_geometry, border=grDevices::heat.colors(nrow(bkern1$isocontours)), add=T)
plot(mkern1$wmKDE, main='Mean Kernel \n n(ID1) = 500; n(ID2) = 500')
plot(st_geometry(locs), pch=20, cex = 0.6, col=locs$col, add=T)
plot(mkern1$isocontours %>% st_geometry, border=grDevices::heat.colors(nrow(mkern1$isocontours)), add=T) 
## 3-D volumes (identical)
persp(bkern1$wmKDE, main = 'Pooled Kernel \n n(ID1) = 500; n(ID2) = 500')
persp(mkern1$wmKDE, main = 'Mean Kernel \n n(ID1) = 500; n(ID2) = 500')
Case 2) Two individual point patterns with different centres of activity and unequal sample size
locs2 <- bind_rows(data.frame(id = 1, col=8, x = rnorm(n=500), y = rnorm(n=500)),
data.frame(id = 2, col=4, x = rnorm(n=100, mean=-3.5), y = rnorm(n=100, mean=-3.5))) %>%
st_as_sf(., coords=c('x','y'), crs=st_crs(32198))
## pooled kernel
bkern2 <- wmKDE(locs2, avg=F, spatres=0.1, verbose=F)
## mean kernel (different results)
mkern2 <- wmKDE(locs2, id='id', ncores=1, avg=T, spatres=0.1, verbose=F)
## plot estimates side by side
par(mfrow=c(1,2))
plot(bkern2$wmKDE$iso, main='Pooled Kernel \n n(ID1) = 500; n(ID2) = 100', legend=F)
plot(st_geometry(locs2), pch=20, cex = 0.6, col=locs2$col, add=T)
plot(bkern2$isocontours %>% st_geometry, border=heat.colors(nrow(bkern2$isocontours)), add=T)
plot(mkern2$wmKDE, main='Mean Kernel \n n(ID1) = 500; n(ID2) = 100', legend=F)
plot(st_geometry(locs2), pch=20, cex = 0.6, col=locs2$col, add=T)
plot(mkern2$isocontours %>% st_geometry, border=heat.colors(nrow(mkern2$isocontours)), add=T) 
## 3-D volumes (the mean UD correctly assigns equal weight to each UD)
persp(bkern2$wmKDE, main = 'Pooled Kernel \n n(ID1) = 500; n(ID2) = 100')
persp(mkern2$wmKDE, main = 'Mean Kernel \n n(ID1) = 500; n(ID2) = 100')
Case 3) Two individual point patterns with different centres of activity and identical sample sizes, but different degrees of spread (dispersion)
locs3 <- bind_rows(data.frame(id = 1, col=8, x = rnorm(n=500, sd=0.25), y = rnorm(n=500, sd=0.25)),
data.frame(id = 2, col=4, x = rnorm(n=500, mean=-3.5), y = rnorm(n=500, mean=-3.5))) %>%
st_as_sf(., coords=c('x','y'), crs=st_crs(6624))
## pooled kernel
bkern3 <- wmKDE(locs3, avg=F, spatres=0.1, verbose=F)
## mean kernels with and without individual rescaling
mkern3a <- wmKDE(locs3, id='id', ncores=1, avg=T, zscale=F, spatres=0.1, verbose=F)
mkern3b <- wmKDE(locs3, id='id', ncores=1, avg=T, zscale=T, spatres=0.1, verbose=F)
## plot resulting estimates side by side
par(mfrow=c(1,3))
plot(bkern3$wmKDE, legend=F, main = 'Pooled Kernel \n sd(ID1) = 0.25; sd(ID2) = 2')
plot(st_geometry(locs3), cex = 0.6, pch=20, col=locs$col, add=T)
plot(bkern3$isocontours %>% st_geometry, border=heat.colors(nrow(bkern3$isocontours)), add=T)
plot(mkern3a$wmKDE, legend=F, main = 'Mean Kernel w/o rescaling \n sd(ID1) = 0.25; sd(ID2) = 2')
plot(st_geometry(locs3), cex = 0.6, pch=20, col=locs$col, add=T)
plot(mkern3a$isocontours %>% st_geometry, border=heat.colors(nrow(mkern3a$isocontours)), add=T)
plot(mkern3b$wmKDE, legend=F, main = 'Mean Kernel w/rescaling \n sd(ID1) = 0.25; sd(ID2) = 2')
plot(st_geometry(locs3), cex = 0.6, pch=20, col=locs$col, add=T)
plot(mkern3b$isocontours %>% st_geometry, border=heat.colors(nrow(mkern3b$isocontours)), add=T) 
## 3-D volumes (rescaling correctly assigns equal weight to each UD)
par(mfrow=c(1,3))
persp(bkern3$wmKDE, main = 'Pooled Kernel \n sd(ID1) = 0.25; sd(ID2) = 2')
persp(mkern3a$wmKDE, main = 'Mean Kernel w/o rescaling \n sd(ID1) = 0.25; sd(ID2) = 2')
persp(mkern3b$wmKDE, main = 'Mean Kernel w/rescaling \n sd(ID1) = 0.25; sd(ID2) = 2')
Case 4) Two individual point patterns with different centres of activity and degree of spread (dispersion) and identical sample size, with and without spatial weights
locs4 <- dplyr::bind_rows(data.frame(id = 1, w = 0.025, col=8, x = rnorm(n=400, mean=1.75, sd=0.25), y = rnorm(n=400, mean=1.75, sd=0.25)),
data.frame(id = 1, w = 1, col=8, x=rnorm(n=100, sd=0.8), y = rnorm(n=100, sd=0.8)),
data.frame(id = 2, w = 1, col=4, x = rnorm(n=500, mean=-3.5), y = rnorm(n=500, mean=-3.5))) %>%
st_as_sf(., coords=c('x','y'), crs=st_crs(6624))
## 4 a) basic kernel
bkern4 <- wmKDE(locs4, avg=F, spatres=0.1, verbose=F)
## 4 b) mean kernel, no rescaling
mkern4a <- wmKDE(locs4, avg=T, id='id', spatres=0.1, ncores=1, verbose=F)
## 4 c) mean kernel w/spatial weights
mkern4b <- wmKDE(locs4, avg=T, id='id', spw = 'w', spatres=0.1, ncores=1, verbose=F)
## 4 d) plot estimates side by side
par(mfrow=c(1,3))
plot(bkern4$wmKDE, legend=F, main = 'Pooled Kernel')
plot(st_geometry(locs4), cex = 0.6, pch=20, col=locs4$col, add=T)
plot(bkern4$isocontours %>% st_geometry, border=heat.colors(nrow(bkern4$isocontours)), add=T)
plot(mkern4a$wmKDE, legend=F, main = 'Mean Kernel w/o spatial weights')
plot(st_geometry(locs4), cex = 0.6, pch=20, col=locs$col, add=T)
plot(mkern4a$isocontours %>% st_geometry, border=heat.colors(nrow(mkern4a$isocontours)), add=T)
plot(mkern4b$wmKDE, legend=F, main = 'Mean Kernel w/spatial weights (w=0.025)')
plot(st_geometry(locs4), cex = 0.6, pch=20, col=locs$col, add=T)
plot(mkern4b$isocontours %>% st_geometry, border=heat.colors(nrow(mkern4b$isocontours)), add=T) 
## 3-D
par(mfrow=c(1,3))
persp(bkern4$wmKDE, main = 'Pooled Kernel')
persp(mkern4a$wmKDE, main = 'Mean Kernel w/o spatial weights')
persp(mkern4b$wmKDE, main = 'Mean Kernel w/spatial weights')
Footnotes
-
Although current implementation requires a geo-spatially referenced vector feature input (sf) object, this method can (and may eventually) be generalized to any comparable bivariate distribution problem. ↩