ctmcmove-package: ctmcmove

Description Details Author(s) References Examples

Description

Software to facilitates taking movement data in xyt format and pairing it with raster covariates within a continuous time Markov chain (CTMC) framework. As described in Hanks et al. (2015) <DOI:10.1214/14-AOAS803> , this allows flexible modeling of movement in response to covariates (or covariate gradients) with model fitting possible within a Poisson GLM framework.

Details

Typical work flow for analysis of telemetry / GPS movement data:

1. Fit a quasi-continuous path model to telemetry xyt data. The ctmcmove package facilitates this through the "mcmc.fmove" function.

2. Create or import raster layers (from package "raster") for each covariate.

3. Impute a quasi-continuous path (done jointly with model fitting in the "mcmc.fmove" function.

4. Turn the quasi-continuous path into a CTMC discrete-space path using the "path2ctmc" command.

5. Turn discrete-space path into Poisson GLM format using the "ctmc2glm" command.

6. Repeat #3 - #5 multiple times (M times). Stack together the response "z", model matrix "X", and offset "tau" elements from each imputed path.

7. Fit a Poisson GLM model to the stacked data with response "z", model matrix "X", offset "log(tau)", and weights for each row equal to "1/M".

7 (alternate). Alternately, multiple imputation could be used, as described in Hanks et al., (2015).

Author(s)

Ephraim M. Hanks

Maintainer: Ephraim M. Hanks

References

Hanks, E. M.; Hooten, M. B. & Alldredge, M. W. Continuous-time Discrete-space Models for Animal Movement The Annals of Applied Statistics, 2015, 9, 145-165

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

Hooten, M. B.; Johnson, D. S.; Hanks, E. M. & Lowry, J. H. Agent-Based Inference for Animal Movement and Selection Journal of Agricultural, Biological, and Environmental Statistics, 2010, 15, 523-538

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
## Not run: 

##
## 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/4)
## 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="CAR",n.mcmc=400,a=1,r=1,num.paths.save=30)
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
##
##########################################################################

P=20

plot(sst,col=grey.colors(100))
for(i in 1:P){
    points(out$pathlist[[i]]$xy,col=i,type="l",lwd=2)
}
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,method="LinearInterp")
## alternate method, useful if you have impassible barriers, but slower
## ctmc=path2ctmc(path$xy,path$t,int,method="ShortestPath")

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.list=list()
glm.list[[1]]=ctmc2glm(ctmc,loc.stack,grad.stack)

str(glm.list)

for(i in 2:P){
    cat(i," ")
    path=out$pathlist[[i]]
    ctmc=path2ctmc(path$xy,path$t,int,method="LinearInterp")
    glm.list[[i]]=ctmc2glm(ctmc,loc.stack,grad.stack)
}

## remove transitions that are nearly instantaneous
##  (These are essentially outliers in the following regression analyses)
for(i in 1:P){
    idx.0=which(glm.list[[i]]$tau<10^-5)
    if(length(idx.0)>0){
        glm.list[[i]]=glm.list[[i]][-idx.0,]
    }
    glm.list[[i]]$t=glm.list[[i]]$t-min(glm.list[[i]]$t)
}


##
## Stack the P imputations together
##

glm.data=glm.list[[1]]
for(i in 2:P){
    glm.data=rbind(glm.data,glm.list[[i]])
}

str(glm.data)

##########################################################################
##
## 6. Fit Poisson GLM
##    (here we are fitting all "M" paths simultaneously,
##     giving each one a weight of "1/M")
##
##########################################################################

fit.SWL=glm(z~cha+npp+sst+crw+d2r+sst.loc,
        weights=rep(1/P,nrow(glm.data)),family="poisson",offset=log(tau),data=glm.data)
summary(fit.SWL)

beta.hat.SWL=coef(fit.SWL)
beta.se.SWL=summary(fit.SWL)$coef[,2]

##########################################################################
##
## 6. Fit Poisson GLM
##    (here we are fitting using Multiple Imputation
##
##########################################################################

## Fit each path individually
glm.fits=list()
for(i in 1:P){
    glm.fits[[i]]=glm(z~cha+npp+sst+crw+d2r+sst.loc,
        family="poisson",offset=log(tau),data=glm.list[[i]])
}

## get point estimates and sd estimates using Rubin's MI combining rules
beta.hat.mat=integer()
beta.se.mat=integer()
for(i in 1:P){
    beta.hat.mat=rbind(beta.hat.mat,coef(glm.fits[[i]]))
    beta.se.mat=rbind(beta.se.mat,summary(glm.fits[[i]])$coef[,2])
}

beta.hat.mat
beta.se.mat

## E(beta) = E_paths(E(beta|path))
beta.hat.MI=apply(beta.hat.mat,2,mean)
beta.hat.MI

## Var(beta) = E_paths(Var(beta|path))+Var_paths(E(beta|path))
beta.var.MI=apply(beta.se.mat^2,2,mean)+apply(beta.hat.mat,2,var)
beta.se.MI=sqrt(beta.var.MI)

cbind(beta.hat.MI,beta.se.MI)

##
## compare estimates from MI and Stacked Weighted Likelihood approach
##

## standardize regression coefficients by multiplying by the SE of the X matrix
sds=apply(model.matrix(fit.SWL),2,sd)
sds[1]=1

## plot MI and SWL regression coefficients
par(mfrow=c(1,2))
plot(beta.hat.MI*sds,beta.hat.SWL*sds,main="(a) Coefficient Estimates",
xlab="Weighted Likelihood Coefficient",
ylab="Multiple Imputation Coefficient",pch=20,cex=2)
abline(0,1,col="red")
plot(log(beta.se.MI),log(beta.se.SWL),
main="(b) Estimated log(Standard Errors)",xlab="Weighted Likelihood log(SE)",
 ylab="Multiple Imputation log(SE)",pch=20,cex=2)
abline(0,1,col="red")


###########################################################################
##
## 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+sst.loc+s(t,by=-d2r),
        weights=rep(1/P,nrow(glm.data)),family="poisson",offset=log(tau),data=glm.data)
summary(fit)

plot(fit)
abline(h=0,col="red")





############################################################
##
## Overview Plot
##
############################################################


## pdf("sealfig.pdf",width=8.5,height=8.85)
par(mfrow=c(3,3))
##
plot(sst,col=(terrain.colors(30)),main="(a) Sea Surface Temperature")
points(xyt[1,1:2]-c(0,.05),type="p",pch=17,cex=2,col="red")
points(xyt[,1:2],type="b",pch=20,cex=.75,lwd=1)
##
plot(d2r/1000,col=(terrain.colors(30)),main="(b) Distance to Rookery")
points(xyt[1,1:2]-c(0,.05),type="p",pch=17,cex=2,col="red")
points(xyt[,1:2],type="b",pch=20,cex=.75,lwd=1)
##
image(sst,col=rev(terrain.colors(30)),main="(c) Imputed Functional Paths",xlab="",ylab="")
for(i in 1:5){
    ## points(out$pathlist[[i]]$xy,col=i+1,type="l",lwd=3)
    points(out$pathlist[[i]]$xy,col=i+1,type="l",lwd=2)
}
points(xyt[,1:2],type="p",pch=20,cex=.75,lwd=1)
##
ee=extent(c(188.5,190.5,58.4,59.1))
sst.crop=crop(sst,ee)
bg=sst.crop
values(bg)=NA
for(i in c(2)){
    values(bg)[cellFromXY(bg,out$pathlist[[i]]$xy)] <- 1
}
image(sst.crop,col=(terrain.colors(30)),xlim=c(188.85,190.2),
ylim=c(58.5,59),main="(d) CTMC Path",xlab="",ylab="")
image(bg,col="blue",xlim=c(188.85,190.2),ylim=c(58.5,59),add=TRUE)
for(i in c(2)){
    points(out$pathlist[[i]]$xy,col=i,type="l",lwd=3)
}
points(xyt[,1:2],type="b",pch=20,cex=2,lwd=2)
##
image(sst.crop,col=(terrain.colors(30)),xlim=c(189.62,189.849),
ylim=c(58.785,58.895),main="(e) CTMC Model Detail",xlab="",ylab="")
abline(v=189.698+res(sst)[1]*c(-1,0,1,2))
abline(h=58.823+res(sst)[2]*c(-1,0,1,2))
##
plot(fit,main="(f) Time-Varying Response to Rookery",shade=TRUE,
shade.col="orange",lwd=3,rug=F,xlab="Day of Trip",
ylab="Coefficient of Distance To Rookery")
abline(h=0,col="red")
##



###############################################
##
## Get UD (following Kenady et al 2017+)
##
###############################################

RR=get.rate.matrix(fit.SWL,loc.stack,grad.stack)
UD=get.UD(RR,method="lu")
ud.rast=sst
values(ud.rast) <- as.numeric(UD)
plot(ud.rast)


###############################################
##
## Get shortest path and current maps (following Brennan et al 2017+)
##
###############################################

library(gdistance)

## create a dummy transition layer from a raster.
## make sure the "directions" argument matches that used in path2ctmc
## also make sure to add the "symm=FALSE" argument
trans=transition(sst,mean,directions=4,symm=FALSE)
## now replace the transition object with the "rate" matrix
## so "conductance" values are "transition rates"
transitionMatrix(trans) <- RR
str(trans)

##
## now calculate least cost paths using "shortestPath" from gdistance
##

## pick start and end locations
plot(sst)
st=c(185,59.5)
en=c(190,57.3)

st.cell=cellFromXY(sst,st)
en.cell=cellFromXY(sst,en)

## shortest path
sp=shortestPath(trans,st,en,output="SpatialLines")
plot(sst,main="Shortest Path (SST in background)")
lines(sp,col="brown",lwd=7)



##
## Now calculate "current maps" that show space use of random walkers
## moving between two given locations.
##
## gdistance's "passage" function allows for asymmetric transition rates
##

passage.gdist=passage(trans,st,en,theta=.001,totalNet="net")
plot((passage.gdist))



## End(Not run)

ctmcmove documentation built on April 20, 2018, 5:03 p.m.