1 Foraging Model
- 1.1 Dive profiles per indidivuals
- 1.2 Data Statistics before track cut
2 Data Cleaning Rules.
- 2.1 Data after track cut
- 2.2 Chains
3 Chains
4 Temporal Variation in Dive Behavior
5 Behavior Sequence
6 Spatial Prediction
7 Residual Explorations

1 Foraging Model

The proportion of time an animal is in a feeding behavioral state.

Process Model

\[Y_{i,t+1} \sim Multivariate Normal(d_{i,t},σ)\]

\[d_{i,t}= Y_{i,t} + γ_{s_{i,g,t}}*T_{i,g,t}*( Y_{i,g,t}- Y_{i,g,t-1} )\]

\[ \begin{matrix} \alpha_{i,1,1} & \beta_{i,1,1} & 1-(\alpha_{i,1,1} + \beta_{i,1,1}) \\ \alpha_{i,2,1} & \beta_{i,2,1} & 1-(\alpha_{i,2,1} + \beta_{i,2,1}) \\ \alpha_{i,3,1} & \beta_{i,3,1} & 1-(\alpha_{i,3,1} + \beta_{i,3,1}) \\ \end{matrix} \] \[logit(\phi_{Behavior}) = \alpha_{Behavior_{t-1}} \] The behavior at time t of individual i on track g is a discrete draw. \[S_{i,g,t} \sim Cat(\phi_{traveling},\phi_{foraging},\phi_{resting})\]

Dive information is a mixture model based on behavior (S)

\(\text{Average dive depth}(\psi)\) \[ DiveDepth \sim Normal(dive_{\mu_S},dive_{\tau_S})\]

if(!newModel){
  load("/Users/Ben/Dropbox/Whales/Dive/WhalePhys.RData")
  newModel=F
}

#get gps data
f<-list.files("Data/Humpback",pattern="Locations",full.names=T,recursive = T)
gdat<-lapply(f,function(x) read.csv(x,stringsAsFactors=F))
gdat<-lapply(gdat,function(x){
  x$Quality<-as.character(x$Quality)
  return(x)
}) 

gdat<-bind_rows(gdat)

#timestamp
gdat$timestamp<-as.POSIXct(gdat$Date,format="%H:%M:%S %d-%b-%Y",tz="GMT")

#
gdat<-gdat[!is.na(gdat$Latitude),]

#crop by extent
d<-SpatialPointsDataFrame(cbind(gdat$Longitude,gdat$Latitude),data=data.frame(gdat),proj4string=CRS("+proj=longlat +datum=WGS84"))

cropoly<-readShapePoly("Data/CutPolygon.shp",proj4string=CRS("+proj=longlat +datum=WGS84"))

b<-d[!is.na(d %over% cropoly)[,2],]

gdat<-b@data

#get dive data files
f<-list.files("Data/Humpback/",pattern="Behavior",full.names=T,recursive = T)
dat<-bind_rows(lapply(f,read.csv))

dat$timestamp<-as.POSIXct(dat$End,format="%H:%M:%S %d-%b-%Y",tz="GMT")

dat$Month<-months(dat$timestamp)
dat$Month<-factor(dat$Month,levels=month.name)
dat$Hour<-strftime(dat$timestamp,format="%H")
dat$Year<-years(dat$timestamp)

### for testing
gdat<-gdat %>% filter(Ptt %in% unique(dat$Ptt))  

gdat<-gdat %>% dplyr::select(Animal=Ptt,timestamp,Quality,Latitude,Longitude)
gdat$Month<-months(gdat$timestamp)
gdat$Month<-factor(gdat$Month,levels=month.name)

dive<-dat %>% filter(What=="Dive")%>% dplyr::select(Animal=Ptt,timestamp,Hour,Month,Year,DepthMax,DepthMin,DurationMax,DurationMin)
dive<-bind_rows(gdat,dive)

#order by timestamp
dive<-dive %>% arrange(timestamp)
mdat<-dive

#One very ugly timestamp that is on land.
mdat<-mdat %>% filter(!Latitude %in% '-64.5407')

1.1 Dive profiles per indidivuals

ggplot(dive[,],aes(x=timestamp,y=DepthMax)) + geom_point(size=0.1) + geom_line(size=0.1) + facet_wrap(~Animal,scales="free",ncol=2) + theme_bw() + labs(x="Date",y="Dive Depth (m)") + scale_y_reverse()

ggsave("Figures/perindividual.jpeg",height=12,width=13)
ggsave("Figures/perindividual.svg",height=12,width=13)

Dive Profiles with Argos timestamps

ggplot(dive[,],aes(x=timestamp,y=-DepthMax)) + geom_point(size=0.1) + geom_line(size=0.1) + facet_wrap(~Animal,scales="free",ncol=2) + theme_bw() + geom_point(data=dive[!is.na(dive$Latitude),],y=0,col="red",size=1,aes(x=timestamp))
ggsave("Figures/perindividual.jpeg",height=6,width=12)

Dive autocorrelation plots

acf_df<-function(x){
  bacf <- acf(x, plot = FALSE)
  bacfdf <- with(bacf, data.frame(lag, acf))
}
acplot<-dive %>% group_by(Animal) %>% arrange(timestamp) %>% filter(!is.na(DepthMax)) %>% do(acf_df(.$DepthMax))
q <- ggplot(data = acplot, mapping = aes(x = lag, y = acf)) +
       geom_hline(aes(yintercept = 0)) +
       geom_segment(mapping = aes(xend = lag, yend = 0))
q + facet_wrap(~Animal)

mybbox<-make_bbox(data=mdat,lat=Latitude,lon=Longitude,f=0.1)
troy <- get_map(location = mybbox, maptype = "toner-background")

attr_troy <- attr(troy, "bb")    # save attributes from original

# change color in raster
troy[troy == "#FFFFFF"] <- "#C0C0C0"
troy[troy == "#000000"] <- "#FFFFFF"

# correct class, attributes
class(troy) <- c("ggmap", "raster")
attr(troy, "bb") <- attr_troy

#nice to have it as a function
makemap<-function(x){
  mybbox<-make_bbox(data=x,lat=Latitude,lon=Longitude,f=0.1)
troy <- get_map(location = mybbox, maptype = "toner-background")

attr_troy <- attr(troy, "bb")    # save attributes from original

# change color in raster
troy[troy == "#FFFFFF"] <- "#C0C0C0"
troy[troy == "#000000"] <- "#FFFFFF"

# correct class, attributes
class(troy) <- c("ggmap", "raster")
attr(troy, "bb") <- attr_troy
return(troy)
}

1.2 Data Statistics before track cut

mdat %>% group_by(Animal) %>% summarize(n=n(),argos=sum(!is.na(Latitude)),dive=sum(!is.na(DepthMax)))

## # A tibble: 11 x 4
##    Animal     n argos  dive
##     <int> <int> <int> <int>
##  1 131111   548   173   375
##  2 131115  1029   179   850
##  3 131116  1932   457  1475
##  4 131127  9275  2495  6780
##  5 131128   112    52    60
##  6 131130  1165   151  1014
##  7 131132  2292   589  1703
##  8 131133  8496  2077  6419
##  9 131134  2098   653  1445
## 10 131136  6844  1970  4874
## 11 154187  1866   486  1380

nrow(mdat)

## [1] 35657

#set max depth to km
#mdat$DepthMax<-mdat$DepthMax/1000

#Specifiy local time
  
#local time
mdat$localtime<-as.POSIXct(format(mdat$timestamp,tz="etc/GMT+3"))
mdat$LocalHour<-as.numeric(strftime(mdat$localtime,format="%H"))

#set duration to minutes
mdat$DurationMax<-mdat$DurationMax/60

#view data
ggmap(troy) + geom_point(data=mdat,aes(x=Longitude, y=Latitude,col=as.factor(Animal)),size=0.3)  + mytheme + labs(col="Animal")

ggsave("Figures/Map.svg")
ggsave("Figures/Map.png",height=6,width=8)

#By individual
for(x in unique(mdat$Animal)){
  mmap<-makemap(mdat %>% filter(Animal==x))
  ggmap(mmap)+geom_path(data=mdat,aes(x=Longitude, y=Latitude,group=Animal),size=.5) + geom_point(data=mdat,aes(x=Longitude, y=Latitude,col=DepthMax)) + theme_bw() + mytheme + scale_color_continuous(low="blue",high="red") + labs(col="Max Dive Depth (km)") + #facet_wrap(~Animal,ncol=2)
}

server <- function(input, output) {

# create a reactive value that will store the click position
data_of_click <- reactiveValues(clickedMarker=NULL)

dat<-mdat %>% filter(Animal=="131133",!is.na(Latitude)) %>% mutate(ID=1:nrow(.))

# Leaflet map with 2 markers
pal<-colorFactor(heat.colors(10),dat$DepthMax)

output$map <- renderLeaflet({
leaflet() %>% 
 setView(lng=mean(dat$Longitude,na.rm=T), lat =mean(dat$Latitude,na.rm=T), zoom=8) %>%
 addTiles(options = providerTileOptions(noWrap = TRUE)) %>% addPolylines(data=dat %>% arrange(Animal,timestamp),~Longitude, ~Latitude,weight=0.25) %>%
 addCircleMarkers(data=dat, ~Longitude , ~Latitude, layerId=~ID, popup=~paste(timestamp), radius=3 , fillColor="red", stroke = TRUE, fillOpacity = 0.1) 
})

# store the click
observeEvent(input$map_marker_click,{
data_of_click$clickedMarker <- input$map_marker_click
})

# Make a barplot or scatterplot depending of the selected point
output$plot=renderPlot({
my_place=dat[dat$ID %in% data_of_click$clickedMarker$id,] 
print(data_of_click$clickedMarker)
print(my_place)

p<-ggplot(mdat[mdat$Animal %in% dat$Animal,],aes(x=timestamp,y=-DepthMax)) + geom_point() + geom_line(size=0.5) + theme_bw() 
if(!is.null(my_place)){
p<-p + geom_vline(data=my_place,aes(xintercept=timestamp),size=4,col='red',alpha=0.5) 
}
return(p)
})
}


ui <- fluidPage(
br(),
leafletOutput("map", height="600px"),
plotOutput("plot", height="300px"),
br()
)

shinyApp(ui = ui, server = server)

2 Data Cleaning Rules.

All tracks start with an observed argos location.
A single track must have one argos location per time step.

##Time is the beginning of the first point.
step_length=6

sxy<-split(mdat,mdat$Animal)

#Cut into tracks
#time diff function
timed<-function(d,step_length){
  
  #Order and startmwith a valid observation
  d<-d %>% arrange(timestamp)
  startpoint<-min(which(!is.na(d$Latitude)))
  d<-d[startpoint:nrow(d),]
  d$interval <- cut(d$timestamp,breaks = seq(min(d$timestamp), max(d$timestamp)+step_length*3600, as.difftime(step_length, units="hours")))
  
  #If there are no argos observations, remove interval.
  remove_interval<-d %>% group_by(interval) %>% summarize(n=sum(!is.na(Latitude))) %>% filter(n==0)
  d<-d %>% filter(!interval %in% remove_interval$interval )
  
  #No empty dive intervals
  remove_dive<-d %>% group_by(interval) %>% summarize(n=sum(!is.na(DepthMax))) %>% filter(n==0)
  d<-d %>% filter(!interval %in% remove_dive$interval )
  
  #Can we have empty intervals?
  d<-d %>% filter(!is.na(interval)) %>% droplevels()
  
  #which intervals are more than 12 hours apart
  Track=c()
  Track[1]<-1
  counter=1
  
  for(x in 2:length(unique(d$interval))){
    difft<-as.numeric(difftime(as.POSIXct(unique(d$interval)[x]),as.POSIXct(unique(d$interval)[x-1]),units="hours"))
    if(difft>step_length){
      counter=counter+1
    } 
    Track[x]<-counter
  }
  
  d<-d %>% inner_join(data.frame(interval=unique(d$interval),Track))
  
  #First position in each track must be a valid position, recode 
  d<-d %>% group_by(Track) %>% filter(timestamp >= min(timestamp[!is.na(Latitude)]))
  
  #get jStep and interval time
  d<-d %>% group_by(interval) %>% mutate(j=difftime(timestamp,as.POSIXct(interval,tz="GMT"),units="hours")/step_length,jStep=1:n())
  
  #set step and remove tracks less than 2 steps, a bit ugly to maintain subgroup order.
  d<-d %>% group_by(Track) %>% mutate(step=as.numeric(as.factor(as.character(interval)))) %>% filter(max(step) > 2)
  
  #reset track numbering
  #need more than 3 jSteps
    #d<-d %>% group_by(Track,step) %>% filter(max(jStep) > 3)

  #remove tracks that are shorter than 12 hours
  track_time<-d %>% group_by(Track) %>% summarize(mt=difftime(max(as.POSIXct(timestamp)),min(as.POSIXct(timestamp)),units="hours")) %>% filter(mt>=12) %>% .$Track
  
  return(d)
  }

sxy<-lapply(sxy,timed,step_length=step_length)
mdat<-bind_rows(sxy)

#recode tracks and steps to make ordinal
mdat<-mdat %>% group_by(Animal) %>% mutate(Track=as.numeric(as.factor(Track)))

2.1 Data after track cut

mdat %>% group_by(Animal) %>% summarize(n=n())

## # A tibble: 11 x 2
##    Animal     n
##     <int> <int>
##  1 131111   454
##  2 131115  1005
##  3 131116  1892
##  4 131127  8224
##  5 131128    99
##  6 131130   483
##  7 131132  1545
##  8 131133  6873
##  9 131134  1870
## 10 131136  6092
## 11 154187  1836

mdat %>% group_by(Animal) %>% summarize(Tracks=length(unique(Track)))

## # A tibble: 11 x 2
##    Animal Tracks
##     <int>  <int>
##  1 131111      1
##  2 131115      2
##  3 131116      1
##  4 131127     29
##  5 131128      1
##  6 131130      2
##  7 131132      5
##  8 131133     34
##  9 131134      6
## 10 131136     15
## 11 154187      2

ggplot(mdat[,],aes(x=timestamp,y=-DepthMax)) + geom_point(size=0.1) + geom_line(size=0.1) + facet_wrap(~Animal,scales="free",ncol=2) + theme_bw() + geom_point(data=dive[!is.na(dive$Latitude),],y=0,col="red",size=1,aes(x=timestamp))

###################################################
#filter by two whales for the moment
#mdat<-mdat %>% filter(Animal %in% c("131132","131127"))
#mdat<-mdat %>% filter(Animal %in% c("131115","131116"))

#remake map
troy<-makemap(mdat)
####################################################

#refactor animal
mdat$jAnimal<-as.numeric(as.factor(mdat$Animal))

#unique track code
mdat<-mdat %>% ungroup() %>% mutate(TrackID=as.numeric(as.factor(paste(Animal,Track,sep="_"))))

##Split into format for JAGS
#Cast time array
mdat<-mdat %>% mutate(j=as.numeric(j))
j<-reshape2::acast(mdat,TrackID~step~jStep,value.var="j")

#how many observations per individual in each step
idx<-mdat %>% group_by(TrackID,step) %>% summarize(n=n())
colnames(idx)<-c("Track","step","jStep")
idx<-reshape2::acast(data=idx,Track~step)

#tracks
tracks<-length(unique(mdat$TrackID))

#steps per track
steps<-mdat %>% group_by(TrackID) %>% summarize(n=length(unique(step))) %>% .$n

#obs array
obs<-reshape2::melt(mdat,measure.vars=c("Longitude","Latitude"))
obs<-reshape2::acast(obs,TrackID~step~jStep~variable)
obs[!is.finite(obs)]<-NA

#argos class array
mdat$argos.lc<-factor(mdat$Quality,levels=c(3,2,1,0,"A","B"))
mdat$numargos<-as.numeric(mdat$argos.lc)
obs_class<-reshape2::acast(mdat,TrackID~step~jStep,value.var="numargos")
#set interpolated observations to having lowest class of argos error
obs_class[!is.finite(obs_class)]<-6

#average dive depth array
maxdive<-reshape2::acast(mdat,TrackID~step~jStep,value.var="DepthMax")

#fill the empty values
maxdive[!is.finite(maxdive)]<-NA

sink(“Bayesian/NestedDive.jags”) cat(" model{

pi <- 3.141592653589

#for each if 6 argos class observation error

for(x in 1:6){

##argos observation error##
argos_prec[x,1:2,1:2] <- argos_cov[x,,]

#Constructing the covariance matrix
argos_cov[x,1,1] <- argos_sigma[x]
argos_cov[x,1,2] <- 0
argos_cov[x,2,1] <- 0
argos_cov[x,2,2] <- argos_alpha[x]
}

for(i in 1:ind){
for(g in 1:tracks[i]){

## Priors for first true location
#for lat long
y[i,g,1,1:2] ~ dmnorm(argos[i,g,1,1,1:2],argos_prec[1,1:2,1:2])

#First movement - random walk.
y[i,g,2,1:2] ~ dmnorm(y[i,g,1,1:2],iSigma)

###First Behavioral State###
state[i,g,1] ~ dcat(lambda[]) ## assign state for first obs

#Process Model for movement
for(t in 2:(steps[i,g]-1)){

#Behavioral State at time T
phi[i,g,t,1] <- alpha[state[i,g,t-1]] 
phi[i,g,t,2] <- 1-phi[i,g,t,1]
state[i,g,t] ~ dcat(phi[i,g,t,])

#Turning covariate
#Transition Matrix for turning angles
T[i,g,t,1,1] <- cos(theta[state[i,g,t]])
T[i,g,t,1,2] <- (-sin(theta[state[i,g,t]]))
T[i,g,t,2,1] <- sin(theta[state[i,g,t]])
T[i,g,t,2,2] <- cos(theta[state[i,g,t]])

#Correlation in movement change
d[i,g,t,1:2] <- y[i,g,t,] + gamma[state[i,g,t]] * T[i,g,t,,] %*% (y[i,g,t,1:2] - y[i,g,t-1,1:2])

#Gaussian Displacement in location
y[i,g,t+1,1:2] ~ dmnorm(d[i,g,t,1:2],iSigma)

}

#Final behavior state
phi[i,g,steps[i,g],1] <- alpha[state[i,g,steps[i,g]-1]] 
phi[i,g,steps[i,g],2] <- 1-phi[i,g,steps[i,g],1]
state[i,g,steps[i,g]] ~ dcat(phi[i,g,steps[i,g],])

##  Measurement equation - irregular observations
# loops over regular time intervals (t)    

for(t in 2:steps[i,g]){

#first substate
sub_state[i,g,t,1] ~ dcat(sub_lambda[state[i,g,t]])

# loops over observed dive within interval t

for(u in 2:idx[i,g,t]){ 

#Substate, resting or foraging dives?
sub_phi[i,g,t,u,1] <- sub_alpha[state[i,g,t],sub_state[i,g,t,u-1]] 
sub_phi[i,g,t,u,2] <- 1-sub_phi[i,g,t,u,1]
sub_state[i,g,t,u] ~ dcat(sub_phi[i,g,t,u,])
}

# loops over observed locations within interval t
for(u in 1:idx[i,g,t]){ 
zhat[i,g,t,u,1:2] <- (1-j[i,g,t,u]) * y[i,g,t-1,1:2] + j[i,g,t,u] * y[i,g,t,1:2]

#for each lat and long
#argos error
argos[i,g,t,u,1:2] ~ dmnorm(zhat[i,g,t,u,1:2],argos_prec[argos_class[i,g,t,u],1:2,1:2])

#for each dive depth
#dive depth at time t
divedepth[i,g,t,u] ~ dnorm(depth_mu[state[i,g,t],sub_state[i,g,t,u]],depth_tau[state[i,g,t],sub_state[i,g,t,u]])T(0,)

#Assess Model Fit

#Fit dive discrepancy statistics - comment out for memory savings
eval[i,g,t,u] ~ dnorm(depth_mu[state[i,g,t],sub_state[i,g,t,u]],depth_tau[state[i,g,t],sub_state[i,g,t,u]])T(0,)
E[i,g,t,u]<-pow((divedepth[i,g,t,u]-eval[i,g,t,u]),2)/(eval[i,g,t,u])

dive_new[i,g,t,u] ~ dnorm(depth_mu[state[i,g,t],sub_state[i,g,t,u]],depth_tau[state[i,g,t],sub_state[i,g,t,u]])T(0,)
Enew[i,g,t,u]<-pow((dive_new[i,g,t,u]-eval[i,g,t,u]),2)/(eval[i,g,t,u])

}
}
}
}

###Priors###

#Process Variance
iSigma ~ dwish(R,2)
Sigma <- inverse(iSigma)

##Mean Angle
tmp[1] ~ dbeta(10, 10)
tmp[2] ~ dbeta(10, 10)

# prior for theta in 'traveling state'
theta[1] <- (2 * tmp[1] - 1) * pi

# prior for theta in 'foraging state'    
theta[2] <- (tmp[2] * pi * 2)

##Move persistance
# prior for gamma (autocorrelation parameter)
#from jonsen 2016

##Behavioral States
#Movement autocorrelation
gamma[1] ~ dbeta(10,2)  
gamma[2] ~ dbeta(2,10)  

#Transition Intercepts
alpha[1] ~ dbeta(1,1)
alpha[2] ~ dbeta(1,1)

#Probability of init behavior switching 
lambda[1] ~ dbeta(1,1)
lambda[2] <- 1 - lambda[1]

#Probability of init subbehavior switching 
sub_lambda[1] ~ dbeta(1,1)
sub_lambda[2] <- 1 - sub_lambda[1]

#Dive Priors
#Foraging dives
depth_mu[2,1]  ~ dunif(50,250)
depth_sigma[1] ~  dunif(0,90)
depth_tau[2,1] <- 1/pow(depth_sigma[1],2) 

#Resting Dives
depth_mu[2,2] ~ dunif(0,30)
depth_sigma[2] ~ dunif(0,20)
depth_tau[2,2] <- 1/pow(depth_sigma[2],2) 

#Traveling Dives
depth_mu[1,1] ~ dunif(0,100)
depth_sigma[3] ~ dunif(0,20)
depth_tau[1,1] <- 1/pow(depth_sigma[3],2) 

#Dummy traveling substate
depth_mu[1,2]<-0
depth_tau[1,2]<-0.01

#Sub states
#Traveling has no substate
sub_alpha[1,1]<-1 
sub_alpha[1,2]<-0

#ARS has two substates, foraging and resting
#Foraging probability
sub_alpha[2,1] ~ dbeta(1,1)
sub_alpha[2,2] ~ dbeta(1,1)

##Argos priors##
#longitudinal argos precision, from Jonsen 2005, 2016, represented as precision not sd

#by argos class
argos_sigma[1] <- 11.9016
argos_sigma[2] <- 10.2775
argos_sigma[3] <- 1.228984
argos_sigma[4] <- 2.162593
argos_sigma[5] <- 3.885832
argos_sigma[6] <- 0.0565539

#latitidunal argos precision, from Jonsen 2005, 2016
argos_alpha[1] <- 67.12537
argos_alpha[2] <- 14.73474
argos_alpha[3] <- 4.718973
argos_alpha[4] <- 0.3872023
argos_alpha[5] <- 3.836444
argos_alpha[6] <- 0.1081156

}"
,fill=TRUE)

sink()

#source jags file
source("Bayesian/NestedDiveRagged.R")

#prior cov shape
R <- diag(c(1,1))
data=list(divedepth=maxdive,argos=obs,steps=steps,R=R,j=j,idx=idx,tracks=tracks,argos_class=obs_class)

#paramters to track
pt<-c("alpha","sub_alpha","gamma","depth_mu","depth_tau","state", "sub_state")

if(newModel){
 system.time(diving<-jags(model.file = "Bayesian/NestedDive_NoInd.jags",data=data,n.chains=2,parameters.to.save=pt,n.iter=30000,n.burnin=28000,n.thin=4,DIC=FALSE,parallel=T,codaOnly = c("state", "sub_state")))
  }

2.2 Chains

print("model complete")

## [1] "model complete"

#bind chains
names(diving$samples)<-1:length(diving$samples)
pc_dive<-melt(lapply(diving$samples, as.data.frame))
pc_dive<-pc_dive %>% dplyr::select(chain=L1,variable,value) %>% group_by(chain,variable) %>% mutate(Draw=1:n()) %>% dplyr::select(Draw,chain,par=variable,value)

#extract parameter name
pc_dive$parameter<-data.frame(str_match(pc_dive$par,"(\\w+)"))[,-1]

#Extract index
splitpc<-split(pc_dive,pc_dive$parameter)

#single index
splitpc[c("alpha","gamma")]<-lapply(splitpc[c("alpha","gamma")],function(x){
    sv<-data.frame(str_match(x$par,"(\\w+)\\[(\\d+)]"))[,3]
    pc<-data.frame(x,Behavior=sv)
    return(pc)
}) 

#double index
splitpc[c("depth_mu","depth_tau","sub_alpha")]<-lapply(splitpc[c("depth_mu","depth_tau","sub_alpha")],function(x){
    sv<-data.frame(str_match(x$par,"(\\w+)\\[(\\d+),(\\d+)"))[,3:4]
    colnames(sv)<-c("Behavior","sub_state")
    pc<-data.frame(x,sv)
    return(pc)
}) 

#3 index
splitpc[c("state")]<-lapply(splitpc[c("state")],function(x){
    sv<-data.frame(str_match(x$par,"(\\w+)\\[(\\d+),(\\d+)]"))[,3:4]
    colnames(sv)<-c("TrackID","step")
    pc<-data.frame(x,sv)
    return(pc)
}) 

#4 index
splitpc[c("sub_state")]<-lapply(splitpc[c("sub_state")],function(x){
    sv<-data.frame(str_match(x$par,"(\\w+)\\[(\\d+),(\\d+),(\\d+)]"))[,3:5]
    colnames(sv)<-c("TrackID","step","jStep")
    pc<-data.frame(x,sv)
    return(pc)
}) 

#bind all matrices back together
pc_dive<-bind_rows(splitpc)
rm(splitpc)

3 Chains

ggplot(pc_dive[!pc_dive$parameter %in% c("dive_new","state","sub_state"),],aes(x=Draw,y=value,col=as.factor(chain))) + geom_line() + facet_wrap(~par,scales="free")

#posteriors
ggplot(pc_dive[!pc_dive$parameter %in% c("dive_new","state","sub_state"),],aes(x=value,fill=as.factor(chain))) + geom_histogram() + facet_wrap(~par,scales="free")

#sum table
sumtable<-pc_dive %>% filter(!parameter %in% c("dive_new","state","eval","sub_state","E","Enew")) %>% group_by(parameter,Behavior,sub_state) %>% summarize(mean=round(mean(value),3),upper=round(quantile(value,0.95),3),lower=round(quantile(value,0.05),3)) 

write.csv(sumtable,"Figures/sumtable.csv",row.names = F)

#Animal lookup table
pc_dive<-mdat %>% ungroup() %>%  dplyr::select(Animal,TrackID,Track)  %>% distinct() %>% mutate(TrackID=as.factor(TrackID)) %>% full_join(pc_dive)

4 Temporal Variation in Dive Behavior

#Take the most common estimate of behavior
Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

#Combine posterior summary
mdat$Animal<-as.factor(mdat$Animal)
mdat$Track<-as.factor(mdat$Track)
mdat$step<-as.factor(mdat$step)

state_est<-pc_dive %>% group_by(Animal,Track,step,parameter) %>%  filter(parameter %in% c("state","sub_state")) %>% summarize(Behavior=Mode(value)) %>% ungroup()  %>% spread(parameter,Behavior) 

state_est %>% group_by(state,sub_state) %>% summarize(n=n())

## # A tibble: 3 x 3
## # Groups:   state [?]
##   state sub_state     n
##   <dbl>     <dbl> <int>
## 1     1         1   232
## 2     2         1   632
## 3     2         2   115

state_est$Animal<-as.factor(state_est$Animal)
state_est$Track<-as.factor(state_est$Track)
state_est$step<-as.factor(state_est$step)

state_est<-state_est %>% inner_join(mdat)

blookup<-data.frame(state=c(1,1,2,2),sub_state=c(1,2,1,2),Behavior=c("Traveling","Dummy","Foraging","Resting"))
state_est<-state_est %>% inner_join(blookup)

ggplot(state_est %>% filter(!is.na(DepthMax)),aes(x=timestamp,y=DepthMax,col=Behavior,group=Track)) + geom_point() + geom_line(size=0.1,aes(group=Track)) + facet_wrap(~Animal,scales="free",ncol=1) + theme_bw() + scale_y_reverse() + labs(x="Date",y="Dive Depth (m)") + scale_color_manual(values=c("violetred3","pink","steelblue"))

ggsave("Figures/TemporalBehavior.svg",height=12,width=10)
ggsave("Figures/TemporalBehavior.jpeg",height=12,width=10)

Grouped by stage

ggplot(state_est %>% filter(!is.na(DepthMax)),aes(x=timestamp,y=-DepthMax,col=as.factor(state),group=step)) + geom_point() + geom_line(size=0.25,aes(group=Track)) + facet_wrap(~Animal,scales="free",ncol=1) + theme_bw()

5 Behavior Sequence

ggplot(state_est,aes(x=timestamp,y=Behavior)) + geom_line(aes(group=1))  + geom_point() + facet_wrap(~Animal,scales="free",ncol=1) + theme_bw()

5.1 Proportion of behaviors

state_est %>% group_by(Behavior) %>% filter(is.na(Latitude)) %>% summarize(n=length(jStep)) %>% mutate(p=n/sum(n)) %>% dplyr::select(-n) %>% spread(Behavior,p)

## # A tibble: 1 x 3
##   Foraging   Resting Traveling
## *    <dbl>     <dbl>     <dbl>
## 1 0.613093 0.1559115 0.2309955

prop<-state_est %>% group_by(Animal,Behavior) %>% filter(is.na(Latitude)) %>% summarize(n=length(jStep)) %>% mutate(p=n/sum(n)) %>% dplyr::select(-n) %>% spread(Behavior,p)

state_est %>% group_by(Behavior,Month) %>% summarize(LQuartile=quantile(DepthMax,0.25,na.rm=T),HQuartile=quantile(DepthMax,0.75,na.rm=T),max=max(DepthMax,na.rm=T),s=sum(DepthMax>400,na.rm = T),mean=mean(DepthMax,na.rm=T),n=n())

## # A tibble: 18 x 8
## # Groups:   Behavior [?]
##     Behavior    Month LQuartile HQuartile   max     s      mean     n
##       <fctr>    <chr>     <dbl>     <dbl> <dbl> <int>     <dbl> <int>
##  1  Foraging    April      54.0    275.50 499.5    43 168.30276  2212
##  2  Foraging February      29.0    245.50 347.5     0 143.94693   512
##  3  Foraging     July      70.5    181.50 363.5     0 128.04462   718
##  4  Foraging     June      98.5    213.50 543.5    24 158.44698  6069
##  5  Foraging    March     125.5    315.50 403.5     1 223.10590  1635
##  6  Foraging      May      90.5    307.50 499.5   284 195.20576  8144
##  7   Resting    April      19.5     57.00 371.5     0  64.43393   391
##  8   Resting February      17.5     47.00 331.5     0  52.49114   776
##  9   Resting     July      19.0     47.50 299.5     0  56.92308   174
## 10   Resting     June      19.0     73.50 395.5     0  62.92764   650
## 11   Resting    March      17.0     38.00 387.5     0  54.02968  1098
## 12   Resting      May      19.5     58.50 395.5     0  58.90488  1114
## 13 Traveling    April      20.5     51.00 110.5     0  37.46882   640
## 14 Traveling February      21.5     38.00 157.5     0  33.13429   386
## 15 Traveling     July      19.0     30.00 157.5     0  27.76562   161
## 16 Traveling     June      21.5     60.50 209.5     0  42.71062  1732
## 17 Traveling    March      18.0     43.25 173.5     0  33.89600   621
## 18 Traveling      May      25.0     64.50 371.5     0  46.77831  3340

5.2 Simulate dive depths

#how many of each state to draw?
statecount<-state_est %>% group_by(Animal,Behavior) %>% filter(is.na(Latitude)) %>% summarize(n=length(jStep))

pred_dives<-list()
for(x in 1:nrow(statecount)){
  state_index<-blookup[ blookup$Behavior %in% statecount$Behavior[x],"state"]
  sub_state_index<-blookup[ blookup$Behavior %in%  statecount$Behavior[x],"sub_state"]

  depth_tau<-pc_dive %>% filter(parameter %in% c("depth_tau"),Behavior==state_index,sub_state==sub_state_index)
  depth_mu<-pc_dive %>% filter(parameter %in% c("depth_mu"),Behavior==state_index,sub_state==sub_state_index)

  travel_dives<-data.frame(Animal=statecount[x,"Animal"],DepthMax=rtruncnorm(n=as.numeric(statecount[x,"n"]),mean=depth_mu$value,sd=1/sqrt(depth_tau$value),a=0),Behavior=statecount[x,"Behavior"])
  pred_dives[[x]]<-travel_dives
}

pred_dives<-bind_rows(pred_dives)

ggplot(pred_dives) + geom_histogram(alpha=0.8,aes(x=DepthMax,fill=Behavior))  + theme_bw() + labs(x="Dive Depth (m)") +facet_wrap(~Behavior,scales="free")

ggplot(pred_dives) + geom_density(alpha=0.8,aes(x=DepthMax,fill=Behavior))  + theme_bw() + labs(x="Dive Depth (m)") + scale_fill_manual(values=c("violetred3","pink","steelblue"))

ggsave("Figures/DiveDensity.jpg",height=4,width=7)

ggplot(pred_dives)  + geom_histogram(data=mdat,aes(x=DepthMax),col="black")+ geom_histogram(alpha=0.8,aes(x=DepthMax,fill=Behavior))  + theme_bw() + labs(x="Dive Depth (m)")

ggsave("Figures/DiveHist.jpg",height=4,width=7)

ggplot(pred_dives)  +  geom_density(alpha=0.8,aes(x=DepthMax),col='red')  + theme_bw() + labs(x="Dive Depth (m)") +geom_density(data=mdat,aes(x=DepthMax),col="black",size=1.5) + theme_bw() + labs(x="Dive Depth (m)")

ggsave("Figures/DivePredict.jpg",height=4,width=7) 

ggplot(pred_dives)  +  geom_density(alpha=0.8,aes(x=DepthMax),col='red')  + theme_bw() + labs(x="Dive Depth (m)") +geom_density(data=mdat,aes(x=DepthMax),col="black",size=1.5) + facet_wrap(~Animal) + theme_bw() + labs(x="Dive Depth (m)")

5.3 Diel

ggplot(data=state_est) + geom_boxplot(aes(x=as.factor(LocalHour),y=DepthMax,fill=Behavior)) + ylab("Dive Depth (m)")+ labs(x="Hour (Local GMT+3)")

ggsave("Figures/Diel.svg")
ggsave("Figures/Diel.png",height=5,width=8,unit="in")

ggplot(data=state_est) + geom_boxplot(aes(x=as.factor(LocalHour),y=DepthMax,fill=Behavior)) + ylab("Dive Depth (m)") + facet_wrap(~Month) + labs(x="Hour (Local GMT+3)")

ggsave("Figures/DielbyMonth.svg")
ggsave("Figures/DielbyMonth.png",height=5,width=10,unit="in")

5.4 Month

state_est$Month<-factor(state_est$Month,levels=month.name)
ggplot(data=state_est) + geom_boxplot(aes(x=Month,y=DepthMax,fill=Behavior)) + ylab("Dive Depth (m)")

ggsave("Figures/Month.svg")
ggsave("Figures/Month.png",height=5,width=8,unit="in")

ggplot(data=state_est) + geom_density(aes(x=DepthMax,fill=Month)) + ylab("Dive Depth (m)") + facet_wrap(~Behavior,scale="free",ncol=1)

6 Spatial Prediction

state_est<-state_est %>% arrange(Animal,Track,timestamp)
ggmap(troy) +geom_path(data=state_est %>% filter(!is.na(Latitude)),aes(x=Longitude, y=Latitude,group=paste(Animal,Track)),size=0.5) + geom_point(data=state_est,aes(x=Longitude, y=Latitude,col=Behavior),size=0.5) + mytheme

ggmap(troy) +geom_path(data=state_est %>% filter(!is.na(Latitude)),aes(x=Longitude, y=Latitude,group=paste(Animal,Track)),size=0.5) + geom_point(data=state_est,aes(x=Longitude, y=Latitude,col=Behavior),size=0.5) + mytheme +facet_wrap(~Animal)

#Show nested behaviors
state_est$TopLayer<-NA
state_est[state_est$state==1,"TopLayer"]<-"Traveling"
state_est[state_est$state==2,"TopLayer"]<-"Area-restricted Search"
state_est$SecondLayer<-state_est$Behavior
state_est<-gather(state_est,key="Layer",value="estimate",TopLayer:SecondLayer)

state_est$Layer<-factor(state_est$Layer,levels = c("TopLayer","SecondLayer"),labels=c("Horizontal Movement","Vertical Movement"))
ggmap(troy) +geom_path(data=state_est %>% filter(!is.na(Latitude)),aes(x=Longitude, y=Latitude,group=paste(Animal,Track)),size=0.1) + geom_point(data=state_est,aes(x=Longitude, y=Latitude,col=estimate),size=0.7) + mytheme +facet_wrap(~Layer) + scale_color_manual(values=c("tomato","violetred3","pink","steelblue")) + labs(col="Behavior") + theme(legend.position="bottom")

ggsave("Figures/NestedMap.jpeg",height=8,width=11)
ggsave("Figures/NestedMap.svg",height=8,width=11)

7 Residual Explorations

7.1 Diel

ggplot(state_est,aes(x=LocalHour,y=DepthMax,col=Behavior)) + geom_smooth(method = "gam", formula = y ~ s(x, k = 5)) + facet_wrap(~Animal)

byhour<-state_est %>% group_by(Animal,LocalHour,Behavior) %>% summarize(n=n()) %>% mutate(prop=n/sum(n))
ggplot(byhour,aes(x=LocalHour,y=prop,col=Behavior)) + geom_line()  + geom_point() + facet_wrap(~Animal,scales="free",ncol=1) + theme_bw() + scale_y_continuous(labels=scales::percent)

Lines connect individuals

byhour<-state_est %>% group_by(LocalHour,Behavior) %>% summarize(n=n()) %>% mutate(prop=n/sum(n))
ggplot(byhour,aes(x=LocalHour,y=prop,col=Behavior)) + geom_line()  + geom_point() + theme_bw() + scale_y_continuous("Frequency of Behavior",labels=scales::percent) + labs(x="Hour") + scale_color_manual(values=c("violetred3","pink","steelblue"))

ggsave("Figures/DielFreq.svg",height=6,width=9)
ggsave("Figures/DielFreq.png",height=6,width=9)

7.2 Monthly

Lines connect individuals

bymonth<-state_est %>% group_by(Animal,Month,Behavior) %>% summarize(n=n()) %>% mutate(prop=n/sum(n)) %>% filter(!is.na(Month)) %>% filter(Behavior %in% c("Foraging","Resting"))
bymonth$Month<-factor(bymonth$Month,levels=month.name)
ggplot(bymonth,aes(x=Month,y=prop,col=Behavior)) + geom_line(aes(group=paste(Animal,Behavior)))  + geom_point() + theme_bw() + scale_y_continuous("Frequency of Behavior",labels=scales::percent)

ggplot(bymonth,aes(x=Month,y=prop,fill=Behavior)) + geom_boxplot() + theme_bw() + scale_y_continuous("Frequency of Behavior",labels=scales::percent) + scale_fill_manual(values=c("violetred3","pink"))

ggsave("Figures/MonthFreq.svg",height=6,width=9)
ggsave("Figures/MonthFreq.png",height=6,width=9)

7.3 Posterior Checks

The goodness of fit is a measured as chi-squared. The expected value is compared to the observed value of the actual data. In addition, a replicate dataset is generated from the posterior predicted intensity. Better fitting models will have lower discrepancy values and be Better fitting models are smaller values and closer to the 1:1 line. A perfect model would be 0 discrepancy. This is unrealsitic given the stochasticity in the sampling processes. Rather, its better to focus on relative discrepancy. In addition, a model with 0 discrepancy would likely be seriously overfit and have little to no predictive power.

fitstat<-pc_dive %>% filter(parameter %in% c("E","Enew")) %>% group_by(parameter,Draw,chain) %>% summarize(fit=mean(value))

fitstat<-dcast(fitstat,Draw+chain~parameter,value.var="fit")

ymin<-min(c(fitstat$E,fitstat$Enew)) - min(c(fitstat$E,fitstat$Enew)) * .1
ymax<-max(c(fitstat$E,fitstat$Enew)) + max(c(fitstat$E,fitstat$Enew)) * .1
ggplot(fitstat,aes(x=E,y=Enew)) + geom_point() + theme_bw() + labs(x="Discrepancy of observed data",y="Discrepancy of replicated data") + geom_abline() + coord_fixed() + ylim(ymin=ymin,ymax=ymax) + xlim(xmin=ymin,xmax=ymax)

ggsave("Figures/PosteriorCheck.jpg",height = 4,width = 4)
fitstat %>% group_by() %>% summarize(mean(E),var(Enew))

#By Animal

fitstat<-pc_dive %>% filter(parameter %in% c("E","Enew")) %>% group_by(parameter,Draw,chain,Animal) %>% summarize(fit=mean(value))

fitstat<-dcast(fitstat,Animal+Draw+chain~parameter,value.var="fit")

ymin<-min(c(fitstat$E,fitstat$Enew)) - min(c(fitstat$E,fitstat$Enew)) * .1
ymax<-max(c(fitstat$E,fitstat$Enew)) + max(c(fitstat$E,fitstat$Enew)) * .1
ggplot(fitstat,aes(x=E,y=Enew)) + geom_point() + theme_bw() + labs(x="Discrepancy of observed data",y="Discrepancy of replicated data") + geom_abline() + coord_fixed() + facet_wrap(~Animal,scale="free")

ggsave("Figures/PosteriorCheckAnimal.jpg",height = 4,width = 4)
fitstat %>% group_by() %>% summarize(mean(E),var(Enew))

#save.image("/Users/Ben/Dropbox/Whales/Dive/WhalePhys.RData")

Foraging Model

Ben Weinstein

2018-01-25 10:05:00

1 Foraging Model

1.1 Dive profiles per indidivuals

1.2 Data Statistics before track cut

2 Data Cleaning Rules.

2.1 Data after track cut

2.2 Chains

3 Chains

4 Temporal Variation in Dive Behavior

5 Behavior Sequence

5.1 Proportion of behaviors

5.2 Simulate dive depths

5.3 Diel

5.4 Month

6 Spatial Prediction

7 Residual Explorations

7.1 Diel

7.2 Monthly

7.3 Posterior Checks