##
## Example of using a CTMC model for movement
##
## Steps:
## 1. Fit Quasi-Continuous Path Model to telemetry data (done using Buderman et al 2015)
## 2. Create covariate raster objects (the CTMC will be on the raster
## grid cells)
## 3. Impute a quasi-continuous path
## 4. Turn quasi-continuous path into a CTMC discrete-space path
## 5. Turn discrete-space path into latent Poisson GLM format
## 6. Fit a Poisson GLM model to the data
##
library(ctmcmove)
data(seal)
xyt=seal$locs[,3:1]
head(xyt)
plot(xyt[,1:2],type="b")
xy=xyt[,-3]
x=xyt[,1]
y=xyt[,2]
t=xyt[,3]
########################
##########################################################################
##
## 1. Fit functional movement model to telemetry data
##
##########################################################################
library(fda)
## Define the knots of the spline expansion.
##
## Problems with fitting the functional movement model can often be fixed by
## varying the spacing of the knots.
knots = seq(min(t),max(t),by=1/6)
## create B-spline basis vectors used to approximate the path
b=create.bspline.basis(c(min(t),max(t)),breaks=knots,norder=3)
## define the sequence of times on which to sample the imputed path
tpred=seq(min(t),max(t),by=1/24/60)
## Fit latent Gaussian model using MCMC
out=mcmc.fmove(xy,t,b,tpred,QQ="CAR2",n.mcmc=400,a=1,r=1,num.paths.save=20)
str(out)
## plot 3 imputed paths
plot(xy,type="b")
points(out$pathlist[[1]]$xy,col="red",type="l")
points(out$pathlist[[2]]$xy,col="blue",type="l")
points(out$pathlist[[3]]$xy,col="green",type="l")
##########################################################################
##
## 2. Creating rasters
##
##########################################################################
cov.df=seal$cov.df
str(cov.df)
NN=sqrt(nrow(cov.df$X))
sst=matrix(seal$cov.df$X$sst,NN,byrow=TRUE)
sst=sst[NN:1,]
sst=raster(sst,xmn=min(seal$cov.df$X$x),xmx=max(seal$cov.df$X$x),
ymn=min(seal$cov.df$X$y),ymx=max(seal$cov.df$X$y))
crs(sst)="+proj=longlat +datum=WGS84"
plot(sst)
chA=matrix(seal$cov.df$X$chA,NN,byrow=TRUE)
chA=chA[NN:1,]
chA=raster(chA,xmn=min(seal$cov.df$X$x),xmx=max(seal$cov.df$X$x),
ymn=min(seal$cov.df$X$y),ymx=max(seal$cov.df$X$y))
crs(chA)="+proj=longlat +datum=WGS84"
pro=matrix(seal$cov.df$X$pro,NN,byrow=TRUE)
pro=pro[NN:1,]
npp=raster(pro,xmn=min(seal$cov.df$X$x),xmx=max(seal$cov.df$X$x),
ymn=min(seal$cov.df$X$y),ymx=max(seal$cov.df$X$y))
crs(npp)="+proj=longlat +datum=WGS84"
int=sst
values(int) <- 1
d2r=int
rookery.cell=cellFromXY(int,xyt[1,1:2])
values(d2r)=NA
values(d2r)[rookery.cell]=0
d2r=distance(d2r)
grad.stack=stack(sst,chA,npp,d2r)
names(grad.stack) <- c("sst","cha","npp","d2r")
plot(sst)
points(xyt[,1:2],type="b")
plot(grad.stack)
##########################################################################
##
## 3 Impute Quasi-Continuous Paths
##
##########################################################################
plot(sst)
for(i in 1:10){
points(out$pathlist[[i]]$xy,col=i,type="l")
}
points(xyt[,1:2],type="b",pch=20,cex=2,lwd=2)
##########################################################################
##
## 4. Turn continuous space path into a CTMC discrete space path
##
##########################################################################
path=out$pathlist[[1]]
ctmc=path2ctmc(path$xy,path$t,int)
str(ctmc)
##########################################################################
##
## 5. Turn CTMC discrete path into latent Poisson GLM data
##
##########################################################################
loc.stack=stack(int,sst)
names(loc.stack) <- c("Intercept","sst.loc")
glm.data=ctmc2glm(ctmc,loc.stack,grad.stack)
str(glm.data)
summary(glm.data)
## now repeat 3-5 for 10 imputations
for(i in 2:10){
path=out$pathlist[[i]]
ctmc=path2ctmc(path$xy,path$t,int)
glm.data=rbind(glm.data,ctmc2glm(ctmc,loc.stack,grad.stack))
}
str(glm.data)
## remove transitions that are nearly instantaneous
summary(glm.data)
idx.0=which(glm.data$tau<10^-5)
idx.0
if(length(idx.0)>0){
glm.data=glm.data[-idx.0,]
}
summary(glm.data)
##########################################################################
##
## 6. Fit Poisson GLM and Poisson GAM with time-varying coefficients
## (here we are fitting all "M" paths simultaneously,
## giving each one a weight of "1/M")
##
##########################################################################
fit=glm(z~cha+npp+sst+crw+d2r+sst.loc,
weights=rep(1/10,nrow(glm.data)),family="poisson",offset=log(tau),data=glm.data)
summary(fit)
##
## 6. (Alternate) We can use any software which fits Poisson glm data.
## The following uses "gam" in package "mgcv" to fit a time-varying
## effect of "d2r" using penalized regression splines. The result
## is similar to that found in:
##
## Hanks, E.; Hooten, M.; Johnson, D. & Sterling, J. Velocity-Based
## Movement Modeling for Individual and Population Level Inference
## PLoS ONE, Public Library of Science, 2011, 6, e22795
##
library(mgcv)
fit=gam(z~cha+npp+crw+s(t,by=-d2r),
weights=rep(1/10,nrow(glm.data)),family="poisson",offset=log(tau),data=glm.data)
summary(fit)
plot(fit)
abline(h=0,col="red")Run the code above in your browser using DataLab