Helper functions for analyzing animal telemetry data with point process models from Devin S. Johnson, Mevin B. Hooten, and Carey E. Kuhn (2013) Estimating animal resource selection from telemetry data using point process models. Journal of Animal Ecology 82:1155--1164.

stpp_rsf_helper.R

## Spatio-temporal point process RSF helper functions

###
### Create spatial polygons object from a SpatialPoints object
###
voronoiPoly <- function(layer){
  require(deldir)
  require(sp)
  crds <- layer@coords
  z <- deldir(crds[,1], crds[,2], rw=as.vector(t(bbox(layer))))
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  voronoi <- SpatialPolygonsDataFrame(SP, 
                                      data=data.frame(
                                        x=crds[,1],y=crds[,2], row.names=sapply(slot(SP, 'polygons'), function(x) slot(x, 'ID'))))
}


###
### Create a likelihood object for a Brownian motion movement model
###
makeBmLik <- function(x,y,time){
  dt <- diff(time)
  dx <- diff(x)
  dy <- diff(y)
  return(function(par){-2*(sum(dnorm(dx,sd=par*sqrt(dt),log=TRUE))+sum(dnorm(dy,sd=par*sqrt(dt), log=TRUE)))})
}


###
### Make variable kernel map for a single animal
###

makeKernHR <- function(track.data, environ.data, h, time.col, id.col){
  uid <- as.character(unique(track.data@data[,id.col]))
  kern.df <- matrix(NA, nrow(environ.data), length(uid)+1)
  colnames(kern.df) <- c(uid,"pop.kern")
  for(i in 1:length(uid)){
    track.data.i <- track.data[track.data@data[,id.col,]==uid[i],]
    delta <- diff(track.data.i@data[,time.col])
    D <- spDists(track.data.i, environ.data)
    kern.df[,i] <- apply(sapply(c(1:length(delta)), function(j, D, delta){dnorm(D[j,], 0, h*sqrt(delta[j]))}, D=D, delta=delta), 1, sum) 
    kern.df[,i] <- kern.df[,i]/sum(kern.df[,i])
  }
  kern.df[,"pop.kern"] <- apply(kern.df[,1:length(uid)], 1, sum)
  return(SpatialPixelsDataFrame(as(environ.data, "SpatialPixels"), data=as.data.frame(kern.df)))
}


###
### Simulate a weighted distribution movement/RSF model from a SpatialPixelsDataFrame
###

simWD <- function(start, times, layers, bm.sigma, coef, model){
  require(sp)
  hab <- layers
  M <- model.matrix(model, layers@data)[,-1]
  lnw <- M%*%matrix(coef,ncol=1)
  maxw <- max(lnw)
  hab@data <- data.frame(w=lnw-maxw)
  size <- 1
  track <- matrix(NA, nrow=length(times), ncol=3)
  track[,1] <- sort(times)
  delta <- diff(track[,1])
  track[1,2:3] <- c(start[1], start[2])
  while(size<length(times)){
    prop <- data.frame(x=rnorm(1,track[size,2],bm.sigma*sqrt(delta[size])), y=rnorm(1,track[size,3],bm.sigma*sqrt(delta[size])))
    coordinates(prop) <- ~x + y
    wprop <- prop %over% hab
    if((!is.na(wprop)) & runif(1,0,1)<=exp(wprop)){
      track[size+1,2:3] <- coordinates(prop)
      size=size+1
    }
  }
  colnames(track) <- c("time", "x", "y")
  track <- as.data.frame(track)
  coordinates(track) <- ~x + y
  return(track)
}


###
### Sample spatial quadrature points and find tile areas for quadrature weights
###

# getSpatialQuad <- function(n, type="hexagonal", track.data, environ.data, time.col, ...){
#   require(sp)
#   sp.quad <- spsample(environ.data, n=n, type=type, ...)
#   if(is.numeric(time.col)) tnms <- colnames(track.data@data)[time.col]
#   else tnms <- time.col
#   sp.quad <- SpatialPointsDataFrame(sp.quad, data=data.frame(obs.time=rep(0,nrow(sp.quad@coords))))
#   #sp.quad <- SpatialPointsDataFrame(sp.quad, data=over(sp.quad, environ.data))
#   sp.quad <- rbind(sp.quad, SpatialPointsDataFrame(SpatialPoints(track.data), data=data.frame(obs.time=track.data@data[,time.col])))
#   sp.quad@bbox <- bbox(environ.data)
#   sp.poly <- voronoiPoly(sp.quad)
#   labels <- sapply(sp.poly@polygons, FUN=function(x){return(x@ID)})
#   area <- sapply(sp.poly@polygons, FUN=function(x){return(x@area)})
#   sp.quad@data <- cbind(data.frame(ID=labels, area=area), sp.quad@data)
#   return(list(sp.quad=sp.quad, sp.quad.poly=sp.poly))
# }

# ###
# ### Sample temporal quadrature points and find interval lengths for quadrature weights
# ###
# getTempQuad <- function(n=200, track.data, time.col){
#   require(sp)
#   if(missing(time.col)) stop("Argument 'time.col' is missing please specify!")
#   times <- track.data@data[,time.col] - min(track.data@data[,time.col])
#   tm.quad <- sort(c(times, seq(min(times), max(times), length=n+2)[-c(1,n+2)]))
#   tm.quad <- tm.quad[!duplicated(tm.quad)]
#   lengths <- diff(c(tm.quad[1], sapply(c(2:length(tm.quad)), function(i,v){mean(v[(i-1):i])}, v=tm.quad), tm.quad[length(tm.quad)]))
#   lengths[2] <- sum(lengths[1:2])
#   out <- data.frame(times=tm.quad[-1], length=lengths[-1], obs=ifelse(tm.quad%in%times, 1, 0)[-1])
#   return(out)
# }


###
### Cross spatial and temporal quadratures to obtain 3d locations and volumes
###
SpatTempQuadrature <- function(track.data, environ.data, time.col, time.int, tile.dim){
  require(sp)
  require(rgeos)
  
  environ.data$cellID <- 1:nrow(environ.data)
  # make quad times
  times <- track.data@data[,time.col] - min(track.data@data[,time.col])
  qt <- seq(0, ceiling(max(times)/time.int)*time.int, time.int)
  qt <- sort(c(times, qt))
  qt <- qt[!duplicated(qt)]
  lengths <- diff(c(qt[1], sapply(c(2:length(qt)), function(i,v){mean(v[(i-1):i])}, v=qt), qt[length(qt)]))
  lengths[2] <- sum(lengths[1:2])
  temp.quad <- data.frame(times=qt[-1], length=lengths[-1], obs=ifelse(qt%in%times, 1, 0)[-1])
  temp.quad$x <- temp.quad$y <- NA
  temp.quad$x[temp.quad$obs==1] <- coordinates(track.data)[-1,1]
  temp.quad$y[temp.quad$obs==1] <- coordinates(track.data)[-1,2]
  
  # Get disperion boundries
  delta <- diff(times)
  d1 <- spDists(track.data)
  pwd <- d1[col(d1)==(row(d1)+1)]
  disp.per.time <- pwd/delta
  cut.vls <- quantile(delta, probs=seq(0,1,0.2))
  cut.vls[6] <- cut.vls[6]*1.1
  cut.vls[1] <- 0
  delta.cat <- cut(delta, cut.vls)
  disp.rad.lookup <- aggregate(disp.per.time, list(delta.cat), FUN=max)
  disp.rad.lookup$x <- disp.rad.lookup$x/tile.dim[1]
  disp.rad.lookup$y <- (disp.rad.lookup$x*tile.dim[1])/tile.dim[2]
  
  # Make spatial quad points at temporal quad points
  cat("\nConstructing spatial x temporal quadrature data\n")
  cat("Please be patient...\n")
  mu <- track.data@coords[1,]
  time.last <- 0
  df <- NULL
  quad.time <- temp.quad$time
  time.last <- 0
  pb <- txtProgressBar(min = 0, max = nrow(temp.quad), style = 3)
  for(i in 1:nrow(temp.quad)){
    delta.last <- quad.time[i]-time.last
    rad.x <- ceiling(disp.rad.lookup$x[findInterval(delta.last, cut.vls)]*delta.last)
    rad.y <- ceiling(disp.rad.lookup$y[findInterval(delta.last, cut.vls)]*delta.last)
    x.grid <- c(seq(mu[1]-rad.x*tile.dim[1], mu[1]-tile.dim[1], tile.dim[1]), mu[1], seq(mu[1]+tile.dim[1], mu[1]+rad.x*tile.dim[1], tile.dim[1]))
    y.grid <-  c(seq(mu[2]-rad.y*tile.dim[2], mu[2]-tile.dim[2], tile.dim[2]), mu[2], seq(mu[2]+tile.dim[2], mu[2]+rad.y*tile.dim[2], tile.dim[2]))
    #grd.bbox <- matrix(c(mu[1]-(rad.x+0.5)*tile.dim[1], mu[2]-(rad.y+0.5)*tile.dim[2], mu[1]+(rad.x+0.5)*tile.dim[1], mu[2]+(rad.y+0.5)*tile.dim[2]), 2, 2)
    #colnames(grd.bbox) <- c("min", "max")
    #rownames(grd.bbox) <- c("x","y")
    spqt <- as(SpatialPoints(expand.grid(x=x.grid, y=y.grid), proj4string=CRS(proj4string(environ.data))), "SpatialPixels")
    spqtPoly <- as(spqt, "SpatialPolygons")
    if(temp.quad$obs[i]==1) spqt <- rbind(SpatialPoints(temp.quad[i,c("x","y")], proj4string=CRS(proj4string(environ.data))), as(spqt,"SpatialPoints"))
    else spqt <- as(spqt, "SpatialPoints")
    
    #spqt@bbox <- grd.bbox    
    out <- as(spqt, "data.frame")
    out$t <- temp.quad$time[i]
    out$delta.last <- delta.last
    out$area <- prod(tile.dim)
    out$length<- temp.quad$length[i]
    out$volume <- out$area*temp.quad$length[i]
    out$response <- rep(0,length(out$area))
    if(temp.quad$obs[i]==1){
      out$area <- out$area/c(2, table(spqt%over%spqtPoly))
      out$response[1] <- 1
    } 
    out$response <- out$response/out$volume[1]  
    cov <- spqt %over% environ.data
    out  <- cbind(out, cbind(dist=spDistsN1(spqt, mu), cov))
    out$bmKern <- -0.5*out$dist^2/out$delta.last
    out <- out[!is.na(cov[,1]),]
    df <- rbind(df, out)
    if(temp.quad$obs[i]==1){
      mu <- as.numeric(temp.quad[i,c("x","y")])
      time.last <- temp.quad$time[i]
    }
    setTxtProgressBar(pb, i)
  }
  close(pb)
  class(df) <- c("stQuad", class(df))
  return(df)
}

###
### Create spatial quadrature data
###
SpatQuadrature <- function(track.data, environ.data){
  require(sp)
  require(rgeos)
    
  # Make spatial quad points at temporal quad points
  cat("\nConstructing spatial quadrature data\n")
  cat("Please be patient...\n")
  saPoly <- as(environ.data, "SpatialPolygons")
  spGrid <- SpatialPointsDataFrame(as(environ.data, "SpatialPoints"), data=data.frame(obs=rep(0,dim(environ.data)[1])))
  spTrack <- SpatialPointsDataFrame(as(track.data, "SpatialPoints"), data=data.frame(obs=rep(1,dim(track.data)[1])))
  spq <- rbind(spTrack, spGrid)
  quad.data <- cbind(coordinates(spq), cellID=spq %over% saPoly, obs=spq@data$obs, spq %over% environ.data) 
  wts <- as.vector(prod(environ.data@grid@cellsize)/(spq %over% aggregate(spq, by=saPoly, FUN=length))[,1])
  quad.data$wts <- wts
  quad.data$response <- quad.data$obs/quad.data$wts
  class(quad.data) <- c("spatQuad", class(quad.data))
  names(quad.data)[1:2] <- c("x","y")
  return(quad.data)
}
  


###
### Fit RSF model with STPP
###

stppRSF <- function(model, Q.data, use.gam=FALSE, ...){
  if(!inherits(Q.data, "stQuad")) stop("\nQ.data is not an object of class 'stQuad'! See function 'SpatTempQuadrature'\n")
  if(!use.gam){
    fit <- glm(model, data=Q.data, family=quasi(link=log, variance=mu), weights=volume, ...)
    fit$aic <- 2*(deviance(fit)/2 + sum(log(fit$prior.weights[fit$y!=0])) + sum(fit$y!=0)) + 2*length(fit$coef)
  }
  else{
    fit <- gam(model, data=Q.data, family=quasi(link=log, variance=mu), weights=volume, ...)
    fit$aic <- 2*(deviance(fit)/2 + sum(log(fit$prior.weights[fit$y!=0])) + sum(fit$y!=0)) + 2*sum(fit$edf)
  }
  fit$maxLogLik <- -(deviance(fit)/2 + sum(log(fit$prior.weights[fit$y!=0])) + sum(fit$y!=0))
  class(fit) <- c("stppRSF.fit", class(fit))
  return(fit)
}

###
### Fit RSF model with spatial process
###

spatRSF <- function(model, Q.data, use.gam=FALSE, ...){
  if(!inherits(Q.data, "spatQuad")) stop("\nQ.data is not an object of class 'spatQuad'! See function 'SpatQuadrature'\n")
  if(!use.gam){
    fit <- glm(model, data=Q.data, family=quasi(link=log, variance=mu), weights=wts, ...)
    fit$aic <- 2*(deviance(fit)/2 + sum(log(fit$prior.weights[fit$y!=0])) + sum(fit$y!=0)) + 2*length(fit$coef)
  }
  else{
    fit <- gam(model, data=Q.data, family=quasi(link=log, variance=mu), weights=wts, ...)
    fit$aic <- 2*(deviance(fit)/2 + sum(log(fit$prior.weights[fit$y!=0])) + sum(fit$y!=0)) + 2*sum(fit$edf)
  }
  fit$maxLogLik <- -(deviance(fit)/2 + sum(log(fit$prior.weights[fit$y!=0])) + sum(fit$y!=0))
  class(fit) <- c("spatRSF.fit", class(fit))
  return(fit)
}

summaryRSF <- function(object){
  if(!(inherits(object, "spatRSF.fit") | inherits(object, "stppRSF.fit"))) stop("\nThis is a summery function for 'stppRSF.fit' and 'spatRSF.fit' objects only/n")
  return(summary(object, disp=1))
}


###
### STPP Likelihood calculation
###
n2logLikStpp <- function(model, Q.data, lower=-Inf, upper=Inf, makeFUN=TRUE){
  mm <- model.matrix(model, Q.data)
  y <- as.vector(Q.data$response)
  if(inherits(Q.data,"spatQuad")) w <- as.vector(Q.data$wts)
  else if(inherits(Q.data,"stQuad")) w <- as.vector(Q.data$volume)
  else stop("Q.data argumaent is not an 'stQuad' or 'spatQuad' data object!\n")
  return(
    function(par){
      if(any(par<lower)) val <- Inf
      else if(any(par>upper)) val <- Inf
      else{
        lnL <- as.vector(mm%*%par)
        val <- -2*crossprod(y*lnL - exp(lnL), w)
      }
      return(val)
      }
    )
}