Commit c612fdd7 authored by Kim Whoriskey's avatar Kim Whoriskey

Merge branch 'dev' into 'master'

Merging latest updates for BN

See merge request !2
parents 8e1746de 91057933
File added
Package: swim
Title: Switching Markov Movement Models
Version: 0.2.4
Date: 2017-02-10
Author: Kim Whoriskey with contributions from C. M. Albertsen
Title: Switching Movement Models
Version: 0.2.5
Date: 2020-03-20
Author: Kim Whoriskey with contributions from C. M. Albertsen
Maintainer: Kim Whoriskey <kwhoriskey@dal.ca>
Description: Fit behavioral hidden Markov models to continuous animal movement data. Currently the model provided is the Hidden Markov Movement Model, the HMMM, which estimates the behavioral state for each location of an animal track. The HMMM is appropriate for highly accurate data only.
Imports: TMB, mvtnorm, sp
Description: Fit behavioral hidden Markov models to spatially continuous animal movement data. Currently the model provided is the Hidden Markov Movement Model, the HMMM, which estimates the behavioral state for each location of an animal track. The HMMM is appropriate for highly accurate data only.
Imports: mvtnorm, TMB
Depends: ggplot2, sp
LinkingTo: TMB, RcppEigen
License: GPL-2
NeedsCompilation: yes
Packaged: 2018-06-28 13:59:07 UTC; kimwhoriskey
RoxygenNote: 7.1.0
# Export
exportPattern(".")
# Generated by roxygen2: do not edit by hand
# Import
import(
TMB,
mvtnorm,
sp
)
useDynLib(swim)
\ No newline at end of file
S3method(print,summary.swim)
S3method(print,swim)
S3method(summary,swim)
export(fitSwim)
export(genTrack)
export(plot.swim)
export(plotSwim)
useDynLib(swim)
#fit the SHMMM
fitSwim <- function(x, ...) UseMethod("swim")
fitSwim <- function(data, ts, regularize=TRUE, pars=list(logitTheta1=0, logitTheta2=0,
logitGamma1=0, logitGamma2=0,
logSdlon=0, logSdlat=0,
logA=matrix(log(1),ncol=1,nrow=2))){
#fit the HMMM
#' @export
fitSwim <- function(data, ts, regularize=TRUE, pars=list(working_theta = c(0, 0),
working_gamma = c(0, 0),
working_tau_lon=0, working_tau_lat=0,
working_A=matrix(c(log(0.75), log(0.25)),ncol=1,nrow=2))){
#data input, n X 2 vector
#ts timestep (hours)
if(any(!class(data$date) %in% c("POSIXt", "POSIXct", "POSIXlt"))){
stop("Date-time vector is not a POSIX class")
......@@ -35,34 +32,9 @@ fitSwim <- function(data, ts, regularize=TRUE, pars=list(logitTheta1=0, logitThe
}
#load TMB
requireNamespace("TMB", quietly=TRUE)
#create a list of input data. The order must match the order in the c++ file!
space = 0:1
mu = c(0.5,0.5)
dats <- list(x = t(iLoc),
# b=as.integer(space),
stateSpace=factor(space),
initDist=as.double(mu))
#create a list of input parameters. Again, order matters!
# if(pars==FALSE){
# parameters <- list(logitTheta1=0, logitTheta2=0,
# logitGamma1=0, logitGamma2=0,
# logSdlon=0, logSdlat=0,
# logA=matrix(log(1),ncol=1,nrow=2))
# } else {
parameters=pars
# }
#Then make an objective function, the negative log likelihood to minimize.
obj <- MakeADFun(dats,parameters,
obj <- TMB::MakeADFun(data=list(x=t(iLoc)),
parameters = pars,
DLL="swim")
#Now we can pass the objective function and its derivatives into any regular R optimizer.
......@@ -72,18 +44,30 @@ fitSwim <- function(data, ts, regularize=TRUE, pars=list(logitTheta1=0, logitThe
#calculate the parameter results
cat("Calculating the standard errors \n")
srep <- summary(sdreport(obj))
srep <- summary(TMB::sdreport(obj))
#calculate the latent behavioral states with the Viterbi algorithm
#get the latent behavioral states with the Viterbi algorithm
states = obj$report()$states
# get the stationary distribution
stationary <- obj$report()$delta
# get the one step ahead (forecast pseudo) residuals
pseudos <- as.data.frame(qnorm(obj$report()$pseudo[-1,]))
names(pseudos) <- c("lon", "lat")
#calculate minimized nll
nll = obj$fn()
#return the object and the parameter results
regData = data.frame(t, iLoc)
names(regData) = c("date", "lon", "lat")
rslts <- list(regData = regData, obj=obj, opt=opt, parameters=srep, states=states+1, time=time, nll=nll)
rslts <- list(regData = regData,
obj=obj, opt=opt,
parameters=srep,
states=states, stationary=stationary,
pseudos=pseudos,
time=time, nll=nll)
class(rslts) <- "swim" #set class (summary, print, and plot functions)
......
#the simulation function
genTrack = function(T=500, theta=c(0,pi), gamma=c(0.7,0.2),
Sigma=matrix(c(0.2*0.2,0,0,0.1*0.1),nrow=2),
alphas=c(0.8,0.9)){
#' Simulate Correlated Random Walk Data
#'
#' Simulate data from either a DCRWS movement process or a conditional autoregressive (step length and turning angle) process. Can simulate with or without measurement error.
#' @export
#' @param n_x Number of simulated true locations.
#' @param start_seed Starting seed
#' @param movement_process Type of process model to simulate from. "DCRWS" implies the first-difference correlated random walk, while "carHMM" is the conditional autoregressive correlated random walk.
#' @param process_pars List of relevant parameters depending on the movement process.
#' @param alpha Matrix of switching probabilities. Determines the number of states you wish to have.
#' @param start_loc Starting location in c(lon, lat). If NULL assumed to be (0,0).
#' @return Simulated data set.
genTrack <- function(n_x=500,
start_seed=42,
movement_process="DCRWS",
process_pars = list(theta = c(0, pi),
gamma = c(0.8,0.3),
sdx = rep(0.1, 2)),
alpha = matrix(c(0.95, 0.05, 0.05, 0.95), nrow=2, ncol=2),
start_loc=NULL){
#accumulators
x = matrix(NA, T, 2)
emp = matrix(NA, T-1, 2)
b = c()
alpha=c(1-alphas[1], alphas[2])
#Initialize locations randomly
x[1, ] = mvtnorm::rmvnorm(1, c(0,0), Sigma)
x[2, ] = mvtnorm::rmvnorm(1, x[1,], Sigma)
set.seed(start_seed)
#Initialize behavioral state
b[1] = rbinom(1, 1, 0.5) + 1
##########################
###### behaviours ######
##########################
#Process equation
for(i in 2:(T-1)){
#evolve states
b[i] = rbinom(1, 1, alpha[b[i-1]]) + 1
m=nrow(alpha)
B <- 1:m
b <- numeric(n_x)
b[1] <- sample(B, size=1, prob=rep(1/m, m))
for(i in 2:n_x) b[i] = sample(B, size=1, prob=alpha[b[i-1],])
##############################
###### true locations ######
##############################
if(movement_process=="DCRWS"){
X <- matrix(nrow=n_x, ncol=2)
if(is.null(start_loc)){
X[1,] <- rep(0,2)
} else {
X[1,] <- start_loc
}
X[2,] = mvtnorm::rmvnorm(1, X[1,], diag(process_pars$sdx, nrow=2))
for(i in 3:n_x){
Tr <- matrix(c(cos(process_pars$theta[b[i]]), -sin(process_pars$theta[b[i]]),
sin(process_pars$theta[b[i]]), cos(process_pars$theta[b[i]])),
nrow=2, byrow=TRUE)
X[i,] <- mvtnorm::rmvnorm(1, X[i-1,] + process_pars$gamma[b[i]] * Tr %*% (X[i-1,] - X[i-2,]), diag(process_pars$sdx))
}
} else if (movement_process=="carHMM") {
# step lengths is gamma distribution
sl <- numeric(n_x)
mu <- process_pars$mrl[b[i]]
sigma <- process_pars$sigma[b[i]]
sl[2] <- rgamma(1, (mu/sigma)^2, scale=sigma^2/mu)
for(i in 3:n_x){
mu <- process_pars$acf[b[i]]*sl[i-1] + (1-process_pars$acf[b[i]])*process_pars$mrl[b[i]]
sigma <- process_pars$sigma[b[i]]
sl[i] <- rgamma(1, shape=(mu/sigma)^2, scale=sigma^2/mu)
}
# turning angles are wrapped cauchy
# parameters should be a list of vectors of the form list(location=c(pi,0,0), rho=c(0.2, 0.6, 0.8))
ta <- numeric(n_x)
for(i in 2:(n_x-1)) ta[i] <- CircStats::rwrpcauchy(1, location=process_pars$location[b[i]], rho=process_pars$rho[b[i]])
# ta[i] <- ta[i] - (2*pi)*floor((ta[i]+pi)/(2*pi))
#evolve movement process
emp[i,1] = cos(theta[b[i]]) * (x[i,1] - x[i-1,1]) + sin(theta[b[i]]) * (x[i,2] - x[i-1,2])
emp[i,2] = -sin(theta[b[i]]) * (x[i,1] - x[i-1,1]) + cos(theta[b[i]]) * (x[i,2] - x[i-1,2])
#add random error
x[i+1,] = mvtnorm::rmvnorm(1, x[i,] + emp[i,] * gamma[b[i]], Sigma)
X.geo <- data.frame(ta=ta, sl=sl)
# turn the turning angles and step lengths into locations
# assume the first location is (0,0), and the first bearing is 0 degrees
X <- getPath(angles=ta*180/pi, type='turning', step_lengths=sl,
first_loc=start_loc, first_bearing=NULL)#*180/pi
} else {
stop("model not implemented...")
}
#Get the last behavioral state
b[T] = rbinom(1, 1, alpha[b[T-1]]) + 1
dat = data.frame(lon=x[,1], lat = x[,2], b)
if(movement_process=="DCRWS"){
rslt <- list(b=b, X=X)
} else {
rslt <- list(b=b, X=X, X.geo=X.geo)
}
return(dat)
return(rslt)
}
\ No newline at end of file
#' @useDynLib swim
dummy <- function() return(NULL)
############################################################
####################### step lengths #######################
############################################################
getStepLengths <- function(dat, units="m", convert_to_rad=TRUE, euclidean=FALSE) {
# units can take one of "m" or "km"
if(units=="m"){
R=6371e3
} else if(units=="km") {
R=6371
} else {
stop("units not supported")
}
# have to convert lon/lat into radians first
if(convert_to_rad) dat <- dat*pi/180
# initialize
n = nrow(dat)
sl <- numeric(n)
# sl[1] <- NA
for(i in 2:n){
a = sin((dat[i,2]-dat[i-1,2])/2)*sin((dat[i,2]-dat[i-1,2])/2) +
cos(dat[i-1,2])*cos(dat[i,2])*sin((dat[i,1]-dat[i-1,1])/2)*sin((dat[i,1]-dat[i-1,1])/2)
b = 2*atan2(sqrt(a), sqrt(1-a))
sl[i] = R*b
}
return(sl)
}
#############################################################
################ bearings and turning angles ################
#############################################################
getBearing <- function(loc1, loc2, euclidean=FALSE){
# takes coordinates in radians, returns bearing in radians
#lon,lat pairs
if(euclidean){
num <- loc2[2]-loc1[2]
den <- loc2[1]-loc1[1]
} else {
num <- sin(loc2[1]-loc1[1])*cos(loc2[2])
den <- cos(loc1[2])*sin(loc2[2]) - sin(loc1[2])*cos(loc2[2])*cos(loc2[1]-loc1[1])
}
bear <- atan2(num, den)
return(bear)
}
getAngles <- function(dat, type="spherical", convert_to_rad=TRUE){
# type can be one of "spherical", "spherical.adj", and "euclidean"
# have to convert lon/lat into radians first
if(convert_to_rad) dat <- dat*pi/180
# dat <- obs$X
n <- nrow(dat)
ta <- numeric(n)
# ta[1] <- ta[n] <- NA
for(i in 2:(n-1)){
if(type=="spherical.adj"){
first.bearing <- getBearing(loc1=dat[i-1,], loc2=dat[i,], euclidean=FALSE)*180/pi
second.bearing <- getBearing(loc1=dat[i,], loc2=dat[i+1,], euclidean=FALSE)*180/pi
#this one is just the intial headings of both vectors
} else {
if(type=="spherical"){
first.bearing <- 180 + getBearing(loc1=dat[i,], loc2=dat[i-1,], euclidean=FALSE)*180/pi
second.bearing <- getBearing(loc1=dat[i,], loc2=dat[i+1,], euclidean=FALSE)*180/pi
# this one is the initial heading of the second (in time) vector, and the final heading
# of the first (in time) vector
} else {
if(type=="euclidean"){
first.bearing <- getBearing(loc1=dat[i,], loc2=dat[i+1,], euclidean=TRUE)*180/pi
second.bearing <- getBearing(loc1=dat[i-1,], loc2=dat[i,], euclidean=TRUE)*180/pi
#reverse the order because the angles are counterclockwise for the euclidean bearings
} else {
stop("angle type not implemented")
}
}
}
turning.angle <- (second.bearing-first.bearing)
ta[i] <- turning.angle - 360*floor((turning.angle+180)/360)
# ensure right turns are positive, euclidean angles are calculated counterclockwise from the x-axis, while
# the angles on a sphere are calculated clockwise from North, so we have to reverse the subtraction
# relative to each other to be able to always call negative angles left turns
}
return(ta)
}
#########################################################
##################### recreate path #####################
#########################################################
getNextLoc <- function(start_loc, start_bearing, step_length){
R=6371 # radius of the Earth
# get the next latitude
next.lat = asin( sin(start_loc[2])*cos(step_length/R) +
cos(start_loc[2])*sin(step_length/R)*cos(start_bearing) )
# get the next longitude
next.lon = start_loc[1] + atan2( sin(start_bearing)*sin(step_length/R)*cos(start_loc[2]),
cos(step_length/R)-sin(start_loc[2])*sin(next.lat))
return(c(next.lon, next.lat))
}
getPath <- function(turning_angles, step_lengths,
first_loc=NULL, first_bearing=NULL){
# takes turning angles in degrees, but first.bearing and first.loc are in radians
# spits out locations in radians
# turning angles and step lengths should be as long as the path that you want
# this means that the first two turning angles and first step length won't be used
# if first location and bearing aren't specified, set to c(0,0) and 0 radians
if(is.null(first_loc)) first_loc <- c(0,0)
if(is.null(first_bearing)) first_bearing <- 0
n=length(turning_angles)
locations <- matrix(nrow=n, ncol=2)
locations[1,] <- first_loc
#call get next location
locations[2,] <- getNextLoc(start_loc=first_loc, start_bearing=first_bearing,
step_length=step_lengths[2])
for(i in 3:n){
# get the next bearing by adding the turning angle to the final bearing from x_{t-2} to x_{t-1}
next.bearing <- ((180 + getBearing(loc1=locations[i-1,], loc2=locations[i-2,], euclidean=FALSE)*180/pi) + turning_angles[i])*pi/180
# get the next location
locations[i,] <- getNextLoc(start_loc=locations[i-1,], start_bearing=next.bearing,
step_length=step_lengths[i])
}
return(locations)
}
#plot the swim fits
# this function gets the slope and intercept of a qqline for the residuals below
get_qqline <- function(vec){
y <- quantile(vec[!is.na(vec)], c(0.25, 0.75))
x <- qnorm(c(0.25, 0.75))
m <- diff(y)/diff(x)
b <- y[1] - m * x[1]
return(c(m,b))
}
#' @export
plot.swim <- function(x,...){
data = x$regData
states = x$states
# locations
bplot <- ggplot(data=x$regData, aes(x=lon, y=lat)) +
geom_path() +
geom_point(aes(x=lon, y=lat, col=as.factor(x$states)))+
guides(colour=guide_legend(title="b state")) +
ggtitle("Interpolated Path and Predicted States")
# take the residuals, get rid of the rows with Inf
if(any(is.finite(rowSums(x$pseudos)))) warning("Inf occurred at least once in the pseudoresiduals - all instances were removed.")
pseudos <- x$pseudos[is.finite(rowSums(x$pseudos)),]
plot(data$lat~data$lon, type="o", pch=20, cex=0.7,
col=ifelse(states==2, "blue", "grey"),
xlab="Longitude", ylab="Latitude")
# plot(data[states==1,]$lat~data[states==1,]$lon, type="o", pch=20, cex=0.7, col="grey",
# xlab="Longitude", ylab="Latitude")
# points(data[states==2,]$lat~data[states==2,]$lon, type="p", pch=20, cex=0.7, col="blue")
# plot the residuals
linetmp <- get_qqline(pseudos[,1]) # get qqline parameters
xresidplot <- ggplot(data=pseudos, aes(sample=lon)) +
geom_qq() +
geom_abline(slope=linetmp[1], intercept=linetmp[2]) +
ggtitle("X Coordinate Residuals")
linetmp <- get_qqline(pseudos[,2])
yresidplot <- ggplot2::ggplot(data=pseudos, aes(sample=lat)) +
geom_qq() +
geom_abline(slope=linetmp[1], intercept=linetmp[2]) +
ggtitle("Y Coordinate Residuals")
# put it all together
gridExtra::grid.arrange(bplot,
xresidplot, yresidplot,
widths=c(2,1),
layout_matrix=rbind(c(1,2),
c(1,3)))
}
#' @export
plotSwim = function(object){
requireNamespace("sp", quietly=TRUE)
data(shpfiles)
xlims = c(min(object$regData$lon)-2, max(object$regData$lon)+2) #longitude
......@@ -14,16 +12,15 @@ plotSwim = function(object){
spdf <- SpatialPointsDataFrame(coords = locs, data=locs,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
#make the spatial lines to connect the dots (locations)
line = Lines(list(Line(locs)), ID="path")
line =Lines(list(Line(locs)), ID="path")
sl = SpatialLines(list(line))
#plot everything
plot(world, col="mediumseagreen", xlim=xlims, ylim=ylims, axes=TRUE)
plot(sl, add=TRUE, col="grey")
plot(world, col="lightgrey", xlim=xlims, ylim=ylims, axes=TRUE)
plot(sl, add=TRUE, col="black")
plot(spdf[which(object$states==1),], add=TRUE, pch=20, cex=0.7, type="o",
col="grey")
col="royalblue")
plot(spdf[which(object$states==2),], add=TRUE, pch=20, cex=0.7, type="o",
col="blue")
col="tomato")
}
......
#prints the output from the summary of the swim object
print.summary.swim <- function(x, ...){
#' @export
print.summary.swim <- function(x){
cat("Switching Hidden Markov Movement Model:\n")
cat("Hidden Markov Movement Model:\n")
cat("\n Time to fit: \n")
......
#' @export
print.swim <- function(object){
cat("----------------------------------------------------------- \n the convergence was ... \n\n")
print(object$opt$message)
cat("\n\n\n")
cat("----------------------------------------------------------- \n the associated movement parameters and their standard errors are ... \n\n")
print(round(object$parameters[9:14, ], 4))
cat("\n\n\n")
cat("----------------------------------------------------------- \n the associated TPM is ... \n\n")
print(matrix(round(object$parameters[15:18, 1], 4), nrow=2, byrow=FALSE))
cat("\n\n\n")
cat("----------------------------------------------------------- \n the stationary distribution is ... \n\n")
print(object$stationary)
cat("\n\n\n")
}
#Get the Wald tests and the Confidence intervals for the parameter values
summary.swim <- function(object, ...){
#' @export
summary.swim <- function(object){
time = object$time
......@@ -30,8 +30,12 @@ summary.swim <- function(object, ...){
"sdLon", "sdLat",
"a11", "a21", "a12", "a22"))
# calculate all of the confidence intervals in the working space
lower95 = CIs[1:8,1] - 1.96*CIs[1:8,2]
upper95 = CIs[1:8,1] + 1.96*CIs[1:8,2]
# then apply a transformation to get them on the observed space
# previously just using the delta method gave me confidence intervals outside of the observed space
lower95 = append(lower95, 2*pi/(1.0+exp(-lower95[1])) - pi) #theta1
upper95 = append(upper95, 2*pi/(1.0+exp(-upper95[1])) - pi)
lower95 = append(lower95, 2*pi/(1.0+exp(-lower95[2]))) #theta2
......@@ -46,7 +50,7 @@ summary.swim <- function(object, ...){
upper95 = append(upper95, exp(upper95[6]))
lower95 = append(lower95, 1/(1+exp(-lower95[7]))) #a11, logit
upper95 = append(upper95, 1/(1+exp(-upper95[7])))
lower95 = append(lower95, 1/(1+exp(-lower95[8]))) #a21, note have to logit these ones but I'm not sure if that's correct
lower95 = append(lower95, 1/(1+exp(-lower95[8]))) #a21,
upper95 = append(upper95, 1/(1+exp(-upper95[8])))
lower95 = append(lower95, 1-upper95[15])
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/genTrack.R
\name{genTrack}
\alias{genTrack}
\title{Simulate animal track with switching behaviors}
\description{Simulates an animal path according to the code of the Hidden Markov Movement Model (HMMM). Does not incorporate measurement error.}
\usage{genTrack(T=500, theta=c(0,pi), gamma=c(0.7,0.2),
Sigma=matrix(c(0.2*0.2,0,0,0.1*0.1),nrow=2),
alphas=c(0.8,0.9))}
\arguments{
\item{T}{The length of the animal track.}
\item{theta}{Vector of two turning angles.}
\item{gamma}{Vector of two autocorrelation parameters. Must be between 0 and 1.}
\item{Sigma}{Covariance matrix for the Gaussian process error distribution.}
\item{alphas}{Vector of two switching probabilities: the first is the probability of being in state 1 and staying in state 1, and the second is the probability of being in state 2 and staying in state 2, i.e. alpha_(1,1) and alpha_(2,2) respectively. }
\title{Simulate Correlated Random Walk Data}
\usage{
genTrack(
n_x = 500,
start_seed = 42,
movement_process = "DCRWS",
process_pars = list(theta = c(0, pi), gamma = c(0.8, 0.3), sdx = rep(0.1, 2)),
alpha = matrix(c(0.95, 0.05, 0.05, 0.95), nrow = 2, ncol = 2),
start_loc = NULL
)
}
\arguments{
\item{n_x}{Number of simulated true locations.}
\details{...}
\value{A data.frame with three columns: "lon" is the simulated longitude, "lat" is the simulated latitude, and "b" is the simulated behavioral state.}
\references{...}
\item{start_seed}{Starting seed}
\author{K. Whoriskey}
\item{movement_process}{Type of process model to simulate from. "DCRWS" implies the first-difference correlated random walk, while "carHMM" is the conditional autoregressive correlated random walk.}
\examples{
\item{process_pars}{List of relevant parameters depending on the movement process.}
track = genTrack(500)
plot(track$lat~track$lon, col=ifelse(track$b==1,"red", "black"))
\item{alpha}{Matrix of switching probabilities. Determines the number of states you wish to have.}
\item{start_loc}{Starting location in c(lon, lat). If NULL assumed to be (0,0).}
}
\value{
Simulated data set.
}
\description{
Simulate data from either a DCRWS movement process or a conditional autoregressive (step length and turning angle) process. Can simulate with or without measurement error.
}
print.swim <- function(object){
cat("----------------------------------------------------------- \n the convergence was ... \n\n")
print(object$opt$message)
cat("\n\n\n")
cat("----------------------------------------------------------- \n the associated movement parameters and their standard errors are ... \n\n")
print(round(object$parameters[9:14, ], 4))
cat("\n\n\n")
cat("----------------------------------------------------------- \n the associated TPM is ... \n\n")
print(matrix(round(object$parameters[15:18, 1], 4), nrow=2, byrow=FALSE))
cat("\n\n\n")
cat("----------------------------------------------------------- \n the stationary distribution is ... \n\n")
print(object$stationary)
cat("\n\n\n")
}
This diff is collapsed.
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment