diff --git a/.Rproj.user/8641BD5D/sdb/s-C3ACBEE/449C785F b/.Rproj.user/8641BD5D/sdb/s-C3ACBEE/449C785F index 4dc3ec1..a99d7e9 100644 --- a/.Rproj.user/8641BD5D/sdb/s-C3ACBEE/449C785F +++ b/.Rproj.user/8641BD5D/sdb/s-C3ACBEE/449C785F @@ -1,14 +1,14 @@ { "collab_server" : "", - "contents" : "---\ntitle: \"Foraging Model\"\nauthor: \"Ben Weinstein\"\ndate: \"May 30, 2017\"\noutput: \n html_document:\n toc: true\n number_sections: true\n theme: united\n---\n\n#Abstract\n\nForaging behavior shapes animal health through the balance of the resources expended to obtain nutrients and the energetic rewards of captured resources. The health of individuals will influence reproductive success, fitness and ultimately population viability. Despite the key role of individual health in evolutionary ecology, it can be challenging to connect foraging behavior and energetics in free living marine organisms. Combining data on foraging strategies, feeding behavior and physiology, we estimate body mass for humpback whales as they forage near the Antarctic Peninsula. Humpback whales require extreme energetic storage to facilitate the yearly migration between austral foraging areas to equatorial breeding grounds. Based on tagged individuals, we model the foraging effort, duration and energetics of animals throughout the austral season. We then connect these measurements to studies of whale feeding physiology and biomechanics. While predictive modeling of individual behavior will be naturally complex, we utilize hierarchical Bayesian models to propagate uncertainty throughout the model. By testing a variety of model assumptions surrounding energetics and foraging rate, we can understand the sensitivity of our predictions to parameter bounds. Our model forms the basis for field tests of animal body condition and further exploration of individual behavior and animal health. Future work combining the environmental impacts of changing ocean conditions on foraging success will inform conservation strategies in a changing marine ecosystem.\n\n\n##Time Feeding\n\nThe proportion of time an animal is in a feeding behavioral state.\n\n(@eq1)\n\n\n*Process Model*\n\n$$Y_{i,t+1} \\sim Multivariate Normal(d_{i,t},σ)$$\n\n$$d_{i,t}= Y_{i,t} + γ_(s_(i,g,t) )*T_(i,g,t)*( Y_(i,g,t)- Y_(i,g,t-1) )$$\n\n$$\n\\begin{matrix}\n \\alpha_{i,1,1} & 1-\\alpha_{i,1,1} \\\\\n \\alpha_{i,2,1} & 1-\\alpha_{i,2,1} \\\\\n\\end{matrix}\n$$\n\n\n$$logit(\\phi_{traveling}) = \\alpha_{Behavior_{t-1}}$$\nThe behavior at time t of individual i on track g is a discrete draw.\n$$S_{i,g,t} \\sim Cat(\\phi_{traveling},\\phi_{foraging})$$\n\nDive information is a mixture model based on behavior (S)\n\n$Average dive depth (\\psi)$\n$$ \\psi \\sim Normal(dive_{\\mu_S},dive_{\\tau_S})$$\n\n$\\text{Average number of dives }(\\omega)$\n$$ \\omega \\sim Poisson(\\lambda_S)$$\n\n```{r,echo=F,warning=F,message=F}\nlibrary(tidyr)\nlibrary(ggplot2)\nlibrary(maptools)\nlibrary(raster)\nlibrary(data.table)\n\nlibrary(dplyr)\nlibrary(stringr)\nlibrary(chron)\nlibrary(R2jags)\nlibrary(knitr)\n\nnewModel=T\n\nopts_chunk$set(echo=F,warning=F,message=F,fig.width=11)\n\nmytheme<-theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank(),axis.title.x=element_blank(),axis.title.y=element_blank(),panel.grid=element_blank())\n```\n\n```{r}\n#get dive data\n\nf<-list.files(\"C:/Users/Ben/Dropbox/Whales/Data/Humpback\",pattern=\"Locations\",full.names=T,recursive = T)\n\n#just grab five individuals for now\ngdat<-lapply(f,function(x) read.csv(x,stringsAsFactors=F))\ngdat<-lapply(gdat,function(x){\n x$Quality<-as.character(x$Quality)\n return(x)\n}) \n\ngdat<-bind_rows(gdat)\n\n#timestamp\ngdat$Date<-as.POSIXct(gdat$Date,format=\"%H:%M:%S %d-%b-%Y\",tz=\"Etc/GMT+3\")\n\n#get data files\nf<-list.files(\"C:/Users/Ben/Dropbox/Whales/Data/Humpback/\",pattern=\"Behavior\",full.names=T,recursive = T)\n dat<-bind_rows(lapply(f,read.csv))\n dat$GMTtime<-as.POSIXct(dat$End,format=\"%H:%M:%S %d-%b-%Y\",tz=\"Etc/GMT+3\")\n \n#local time\ndat$timestamp<-as.POSIXct(format(dat$GMTtime,tz=\"Etc/GMT+3\"))\n\ndat$Month<-months(dat$timestamp)\ndat$Month<-factor(dat$Month,levels=month.name)\ndat$Hour<-strftime(dat$timestamp,format=\"%H\")\ndat$Year<-years(dat$timestamp)\n \n#create unique messages\nindices<-which(dat$What==\"Message\")\ncounter=1\ndat$ID<-NA\nfor (x in 1:(length(indices)-1)){\ndat[indices[x]:indices[x+1],\"ID\"]<-counter\ncounter=counter+1\n}\n\ndive<-dat %>% filter(What==\"Dive\")%>% dplyr::select(Animal=Ptt,timestamp,Hour,Month,Year, ID,Start,DepthMax,DepthMin,DurationMax)\ndive<-dive[!is.na(dive$Month),]\n \n#remove duplicate data\ndive<-dive %>% arrange(Animal,timestamp) \ndive<-dive[!duplicated(dive),]\n```\n\n```{r}\n##Merge with geographic data and format for JAGS\nmessages<-dat %>% filter(dat$What==\"Message\")\nmessages$timestamp<-as.POSIXct(messages$End,format=\"%H:%M:%S %d-%b-%Y\",tz=\"Etc/GMT+3\")\n\n#look up the time interval that best fits\nsetDT(gdat) ## convert to data.table by reference\nsetDT(messages) ## same\n\nsetkey(messages, Ptt,timestamp) ## set the column to perform the join on\nsetkey(gdat,Ptt, Date) ## same as above\n\nans = gdat[messages, roll=Inf] ## perform rolling join\nans<-as.data.frame(ans)\n\nmessage_join<-ans %>% select(Date,Animal=Ptt,Date,Quality,Latitude,Longitude,ID)\n\nmdat<-merge(dat,message_join,by=\"ID\")\nmdat<-mdat %>% filter(What==\"Dive\")\n\n#The maxiumum dive depth based on the geographic message\nmdat<-mdat %>% select(ID,Animal,Latitude,Longitude,timestamp,Date,Month,Hour,Year,DepthMax,Quality) \n\nmdat<-mdat %>% group_by(ID,Animal) %>% slice(which.max(DepthMax)) %>% arrange(Animal,timestamp) %>% filter(!is.na(Latitude))\n\n#crop by extent\nd<-SpatialPointsDataFrame(cbind(mdat$Longitude,mdat$Latitude),data=data.frame(mdat),proj4string=CRS(\"+proj=longlat +datum=WGS84\"))\n\ncropoly<-readShapePoly(\"Data/CutPolygon.shp\",proj4string=CRS(\"+proj=longlat +datum=WGS84\"))\n\nb<-d[!is.na(d %over% cropoly)[,2],]\n\nmdat<-b@data\n\n#view data\nggplot(data=mdat)+geom_path(aes(x=Longitude, y=Latitude,group=Animal),size=.5) + geom_point(aes(x=Longitude, y=Latitude,col=DepthMax))+ borders(fill=\"grey90\") + coord_cartesian(ylim = range(mdat$Latitude),xlim=range(mdat$Longitude)) + theme_bw() + mytheme + scale_color_continuous(low=\"blue\",high=\"red\") + labs(col=\"Max Dive Depth(m)\")\n\n##Time is the beginning of the first point.\nstep_length=6\n\nsxy<-split(mdat,mdat$Animal)\n\n#time diff function\ntimed<-function(d,step_length){\n d$j[1]<-0\n for (x in 2:nrow(d)){\n d$j[x]<-as.numeric(difftime(as.POSIXct(d$timestamp[x]),as.POSIXct(d$timestamp[x-1]),units=\"mins\"))/(step_length*60)\n }\n \n #Split out track endings\n ends<-c(1,which(d$j>1),nrow(d))\n\n for(w in 2:length(ends)){\n d[ends[w-1]:ends[w],\"Track\"]<-w-1\n }\n \n #remove tracks that are shorter than three days\n track_time<-d %>% group_by(Track) %>% summarize(mt=difftime(max(as.POSIXct(timestamp)),min(as.POSIXct(timestamp)),units=\"hours\")) %>% filter(mt>=24) %>% .$Track\n \n d<-d[d$Track %in% track_time,]\n \n #renumber the tracks\n d$Track<-as.numeric(as.factor(d$Track))\n return(d)\n }\n\nsxy<-lapply(sxy,timed,step_length=step_length)\n\n#Format matrices for jags\nmdat<-bind_rows(sxy)\n\nsxy<-split(mdat,list(mdat$Animal,mdat$Track),drop=TRUE)\n\nsxy<-lapply(sxy,function(x){\n#How many observations in each step length segment\nx$step<-as.numeric(cut(as.POSIXct(x$timestamp),paste(step_length,\"hours\")))\nreturn(x)\n})\n\nmdat<-bind_rows(sxy)\n\n###################################################\n#filter by two whales for the moment\nmdat<-mdat %>% filter(Animal %in% c(\"131127\",\"131136\"))\n####################################################\n\n#refactor animal\nmdat$Animal<-as.numeric(as.factor(mdat$Animal))\n\n##Split into format for JAGS\n\n#total number of steps per track/animal\nsteps_all<-mdat %>% group_by(Animal,Track) %>% summarize(n=length(unique(step)))\n\n# give each step a label\nmdat<-mdat %>% group_by(Animal,Track,step) %>% mutate(jStep=1:n())\n\n#Cast time array\nj<-reshape2::acast(mdat,Animal~Track~step~jStep,value.var=\"j\")\n\n#how many observations per individual in each step\nmdat$step<-factor(mdat$step,levels=1:max(steps_all$n))\nidx<-reshape2::melt(table(mdat$Animal,mdat$Track,mdat$step))\ncolnames(idx)<-c(\"Animal\",\"Track\",\"step\",\"jStep\")\nidx<-reshape2::acast(data=idx,Animal~Track~step)\n\n#Individuals\nind=length(unique(mdat$Animal))\n\n#tracks per indivudal\ntracks<-mdat %>% group_by(Animal) %>% summarize(tracks=length(unique(Track))) %>% .$tracks\n\n#steps per track\nsteps<-reshape2::acast(steps_all,Animal~Track,value.var=\"n\")\n\n#obs array\nobs<-reshape2::melt(mdat,measure.vars=c(\"Longitude\",\"Latitude\"))\nobs<-reshape2::acast(obs,Animal~Track~step~jStep~variable,fun.aggregate = mean)\n\n#argos class array\nmdat$argos.lc<-factor(mdat$Quality,levels=c(3,2,1,0,\"A\",\"B\"))\nmdat$numargos<-as.numeric(mdat$argos.lc)\nobs_class<-reshape2::acast(mdat,Animal~Track~step~jStep,value.var=\"numargos\",fun.aggregate = min)\n\nobs_class[!is.finite(obs_class)]<-NA\n\n#average dive depth array\nmaxdive<-reshape2::acast(mdat,Animal~Track~step~jStep,value.var=\"DepthMax\",fun.aggregate = mean)\n\n#fill the empty values\nmaxdive[!is.finite(maxdive)]<-NA\n\n#average number of dives\ndivecount<-reshape2::acast(mdat,Animal~Track~step,value.var=\"DepthMax\",fun.aggregate = length)\n\n#no dives, fill empty\ndivecount[divecount==0]<-NA\n\n```\n\n```{r,child=\"Bayesian/Diving.R\",eval=T}\n```\n\n```{r,eval=T}\n#source jags file\nsource(\"Bayesian/Diving.R\")\n\n#prior cov shape\nR <- diag(c(1,1))\ndata=list(divecount=divecount,dive=maxdive,argos=obs,steps=steps,R=R,ind=ind,j=j,idx=idx,tracks=tracks,argos_class=obs_class)\n\n#paramters to track\npt<-c(\"alpha\",\"depth_mu\",\"depth_tau\",\"dive_new\",\"state\",\"lambda\")\n\nif(newModel){\n system.time(diving<-jags.parallel(model.file = \"Bayesian/Diving.jags\",data=data,n.chains=2,parameters.to.save=pt,n.iter=10000,n.burnin=9500,n.thin=2,DIC=FALSE))\n}\n\n```\n\n##Chains\n```{r,eval=T}\n#bind chains\npc_dive<-reshape2::melt(diving$BUGSoutput$sims.array)\ncolnames(pc_dive)<-c(\"Draw\",\"chain\",\"par\",\"value\")\n\n#extract parameter name\npc_dive$parameter<-data.frame(str_match(pc_dive$par,\"(\\\\w+)\"))[,-1]\n\n#Extract index\nsplitpc<-split(pc_dive,pc_dive$parameter)\n\n#single index\nsplitpc[c(\"alpha\",\"depth_mu\",\"depth_tau\",\"lambda\")]<-lapply(splitpc[c(\"alpha\",\"depth_mu\",\"depth_tau\",\"lambda\")],function(x){\n sv<-data.frame(str_match(x$par,\"(\\\\w+)\\\\[(\\\\d+)]\"))[,3]\n pc<-data.frame(x,Behavior=sv)\n return(pc)\n}) \n\n#3 index\nsplitpc[c(\"state\")]<-lapply(splitpc[c(\"state\")],function(x){\n sv<-data.frame(str_match(x$par,\"(\\\\w+)\\\\[(\\\\d+),(\\\\d+),(\\\\d+)\"))[,3:5]\n colnames(sv)<-c(\"Animal\",\"Track\",\"step\")\n pc<-data.frame(x,sv)\n return(pc)\n}) \n\n#4 index\nsplitpc[c(\"dive_new\")]<-lapply(splitpc[c(\"dive_new\")],function(x){\n sv<-data.frame(str_match(x$par,\"(\\\\w+)\\\\[(\\\\d+),(\\\\d+),(\\\\d+),(\\\\d+)]\"))[,3:6]\n colnames(sv)<-c(\"Animal\",\"Track\",\"step\",\"jStep\")\n pc<-data.frame(x,sv)\n return(pc)\n}) \n#bind all matrices back together\npc_dive<-bind_rows(splitpc)\nrm(splitpc)\n\nggplot(pc_dive[!pc_dive$parameter %in% c(\"dive_new\",\"state\"),],aes(x=Draw,y=value,col=as.factor(chain))) + geom_line() + facet_wrap(~par,scales=\"free\")\n\n```\n\n## Empirical distribution of foraging \n\n```{r}\nalpha<-pc_dive %>% filter(parameter==\"alpha\")\nalpha$Behavior<-as.character(alpha$Behavior)\nalpha$Behavior[alpha$Behavior==\"1\"]<-\"Traveling->Traveling\"\nalpha$Behavior[alpha$Behavior==\"2\"]<-\"Foraging->Traveling\"\n\nggplot(alpha,aes(x=value)) + geom_histogram() + facet_wrap(~Behavior) + labs(x=\"Probability of Transition\")\n\n#label previous and prediction state, treat traveling as a success and foraging as a failure in a binomial trial\nalpha[alpha$Behavior==\"Traveling->Traveling\",\"Previous_State\"]<-1\nalpha[alpha$Behavior==\"Foraging->Traveling\",\"Previous_State\"]<-0\n```\n\n### Time foraging as a function of time\n\n* 12 time step\n\n```{r}\nforaging_time<-function(alpha,draws,step){\n states<-c()\n \n #initial state is traveling\n states[1]<-(1)\n \n for(x in 2:draws){\n \n #probability of staying, or switching, to traveling state\n travel_prob=alpha %>% filter(Previous_State==states[x-1]) %>% sample_n(1) %>% .$value\n \n #pick next state\n states[x]<-rbinom(1,1,travel_prob)\n }\n \n #total hours of foraging\n total_foraging<-sum(states==0) * step\nreturn(total_foraging)\n}\n\nstates<-foraging_time(alpha=alpha,draw=20,step=12)\n```\n\nThis allows us to specify the hours foraging while capturing our uncertainty in behavioral states\n\nFor example, if we wanted to know how many hours an animal was predicted to forage on average over a 10 day span. We can draw 100 simulations given the posterior distributions of our empirical model.\n\n```{r}\n\ntest_time<-sapply(1:1000,function(x){\n foraging_time(alpha=alpha,draw=20,step=12)\n})\npaste(\"Mean foraging hours:\", mean(test_time),sep=\" \")\n\nqplot(test_time) +ggtitle(\"Foraging Time (10 Days)\") +xlab(\"Hours\")\n```\n\n\n##Distribution of average max dive depths\n\n```{r}\ndivep<-pc_dive %>% filter(parameter %in% c(\"depth_mu\",\"depth_tau\"))\n\n## Simulate dive depth\ndepth_mu1<-pc_dive %>% filter(parameter %in% c(\"depth_mu\"),Behavior==1)\ndepth_mu2<-pc_dive %>% filter(parameter %in% c(\"depth_mu\"),Behavior==2)\ndepth_tau1<-pc_dive %>% filter(parameter %in% c(\"depth_tau\"),Behavior==1)\ndepth_tau2<-pc_dive %>% filter(parameter %in% c(\"depth_tau\"),Behavior==2)\n\ntravel_dives<-data.frame(MaxDepth=rnorm(1e5,depth_mu1$value,sqrt(1/depth_tau1$value)),Behavior=\"Traveling\")\nforage_dives<-data.frame(MaxDepth=rnorm(1e5,depth_mu2$value,sqrt(1/depth_tau2$value)),Behavior=\"Foraging\")\npred_dives<-bind_rows(travel_dives,forage_dives)\nggplot(pred_dives) + geom_density(alpha=0.5,aes(x=MaxDepth,fill=Behavior)) + geom_density(data=mdat,aes(x=DepthMax),col=\"black\",size=1) + theme_bw() \n\nggplot(pc_dive[pc_dive$parameter==\"dive_new\",])+ geom_density(aes(x=value,fill=par)) + geom_density(data=mdat,aes(x=DepthMax),col=\"black\",size=1.5) + theme_bw() \n```\n\n##Distribution of average dive counts\n\n```{r}\n\n## Simulate dive counts\ncount_lambda1<-pc_dive %>% filter(parameter %in% c(\"lambda\"),Behavior==1)\ncount_lambda2<-pc_dive %>% filter(parameter %in% c(\"lambda\"),Behavior==2)\n\ntravel_counts<-data.frame(Counts=rpois(1e5,count_lambda1$value),Behavior=\"Traveling\")\nforage_counts<-data.frame(Counts=rpois(1e5,count_lambda2$value),Behavior=\"Foraging\")\n\n#calculate observed\nobs_count<-mdat %>% group_by(Animal,Track,step) %>% summarize(Counts=n())\n\npred_counts<-bind_rows(travel_counts,forage_counts)\nggplot(pred_counts) + geom_density(alpha=0.5,aes(x=Counts,fill=Behavior)) + geom_density(data=obs_count,aes(x=Counts),col=\"black\",size=1) + theme_bw() \n```\n\n## Posterior checks\n\n* Compare observed dive depths by inferred behavior against generated data\n\n```{r,eval=F}\npred_dives<-pc_dive %>% filter(parameter==\"dive_new\") %>% group_by(Animal,Track,step,jStep) %>% summarize(Pred_DepthMax=mean(value)) %>% ungroup() %>% mutate(Animal=as.numeric(as.character(Animal)),Track=as.numeric(as.character(Track)),jStep=as.numeric(as.character(jStep)))\n\npredframe<-pred_dives %>% inner_join(mdat,by=c(\"Animal\",\"Track\",\"step\",\"jStep\")) %>% select(Animal,Latitude,Longitude,Track,step,jStep,timestamp,Predicted=Pred_DepthMax,Observed=DepthMax)\n\n#merge with state\npredframe<-pc_dive %>% filter(parameter==\"state\") %>% group_by(Animal,Track,step) %>% summarize(state=names(sort(table(value)))[1]) %>% ungroup() %>% mutate(Animal=as.numeric(as.character(Animal)),Track=as.numeric(as.character(Track))) %>% inner_join(predframe,by=c(\"Animal\",\"Track\",\"step\"))\n\npredframe<-reshape2::melt(predframe, measure.vars=c(\"Predicted\",\"Observed\"),value.name=\"DepthMax\")\n\nggplot(data=predframe) + geom_histogram(aes(x=DepthMax,fill=variable),alpha=0.8) + theme_bw()\n\nggplot(data=predframe) + geom_density(aes(x=DepthMax,fill=variable),alpha=0.8) \n\nggplot(data=predframe) + geom_density(aes(x=DepthMax,fill=variable),alpha=0.8) + facet_wrap(~state) + ggtitle(\"By Behavior\")\n\nggplot(data=predframe) + geom_density(aes(x=DepthMax,fill=variable),alpha=0.8) + facet_wrap(~Animal) + ggtitle(\"By Animal\")\n\n```\n\n```{r}\nsave.image(\"WhalePhys.RData\")\n```\n", + "contents" : "---\ntitle: \"Foraging Model\"\nauthor: \"Ben Weinstein\"\ndate: \"May 30, 2017\"\noutput: \n html_document:\n toc: true\n number_sections: true\n theme: united\n---\n\n#Abstract\n\nForaging behavior shapes animal health through the balance of the resources expended to obtain nutrients and the energetic rewards of captured resources. The health of individuals will influence reproductive success, fitness and ultimately population viability. Despite the key role of individual health in evolutionary ecology, it can be challenging to connect foraging behavior and energetics in free living marine organisms. Combining data on foraging strategies, feeding behavior and physiology, we estimate body mass for humpback whales as they forage near the Antarctic Peninsula. Humpback whales require extreme energetic storage to facilitate the yearly migration between austral foraging areas to equatorial breeding grounds. Based on tagged individuals, we model the foraging effort, duration and energetics of animals throughout the austral season. We then connect these measurements to studies of whale feeding physiology and biomechanics. While predictive modeling of individual behavior will be naturally complex, we utilize hierarchical Bayesian models to propagate uncertainty throughout the model. By testing a variety of model assumptions surrounding energetics and foraging rate, we can understand the sensitivity of our predictions to parameter bounds. Our model forms the basis for field tests of animal body condition and further exploration of individual behavior and animal health. Future work combining the environmental impacts of changing ocean conditions on foraging success will inform conservation strategies in a changing marine ecosystem.\n\n\n##Time Feeding\n\nThe proportion of time an animal is in a feeding behavioral state.\n\n(@eq1)\n\n\n*Process Model*\n\n$$Y_{i,t+1} \\sim Multivariate Normal(d_{i,t},σ)$$\n\n$$d_{i,t}= Y_{i,t} + γ_(s_(i,g,t) )*T_(i,g,t)*( Y_(i,g,t)- Y_(i,g,t-1) )$$\n\n$$\n\\begin{matrix}\n \\alpha_{i,1,1} & 1-\\alpha_{i,1,1} \\\\\n \\alpha_{i,2,1} & 1-\\alpha_{i,2,1} \\\\\n\\end{matrix}\n$$\n\n\n$$logit(\\phi_{traveling}) = \\alpha_{Behavior_{t-1}}$$\nThe behavior at time t of individual i on track g is a discrete draw.\n$$S_{i,g,t} \\sim Cat(\\phi_{traveling},\\phi_{foraging})$$\n\nDive information is a mixture model based on behavior (S)\n\n$Average dive depth (\\psi)$\n$$ \\psi \\sim Normal(dive_{\\mu_S},dive_{\\tau_S})$$\n\n$\\text{Average number of dives }(\\omega)$\n$$ \\omega \\sim Poisson(\\lambda_S)$$\n\n```{r,echo=F,warning=F,message=F}\nlibrary(tidyr)\nlibrary(ggplot2)\nlibrary(maptools)\nlibrary(raster)\nlibrary(data.table)\n\nlibrary(dplyr)\nlibrary(stringr)\nlibrary(chron)\nlibrary(R2jags)\nlibrary(knitr)\n\n\nnewModel=F\n\nopts_chunk$set(echo=F,warning=F,message=F,fig.width=11)\n\nmytheme<-theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank(),axis.title.x=element_blank(),axis.title.y=element_blank(),panel.grid=element_blank())\n```\n\n```{r}\n\nif(!newModel){\n load(\"WhalePhys.RData\")\n}\n\n#get dive data\n\nf<-list.files(\"C:/Users/Ben/Dropbox/Whales/Data/Humpback\",pattern=\"Locations\",full.names=T,recursive = T)\n\n#just grab five individuals for now\ngdat<-lapply(f,function(x) read.csv(x,stringsAsFactors=F))\ngdat<-lapply(gdat,function(x){\n x$Quality<-as.character(x$Quality)\n return(x)\n}) \n\ngdat<-bind_rows(gdat)\n\n#timestamp\ngdat$Date<-as.POSIXct(gdat$Date,format=\"%H:%M:%S %d-%b-%Y\",tz=\"Etc/GMT+3\")\n\n#get data files\nf<-list.files(\"C:/Users/Ben/Dropbox/Whales/Data/Humpback/\",pattern=\"Behavior\",full.names=T,recursive = T)\n dat<-bind_rows(lapply(f,read.csv))\n dat$GMTtime<-as.POSIXct(dat$End,format=\"%H:%M:%S %d-%b-%Y\",tz=\"Etc/GMT+3\")\n \n#local time\ndat$timestamp<-as.POSIXct(format(dat$GMTtime,tz=\"Etc/GMT+3\"))\n\ndat$Month<-months(dat$timestamp)\ndat$Month<-factor(dat$Month,levels=month.name)\ndat$Hour<-strftime(dat$timestamp,format=\"%H\")\ndat$Year<-years(dat$timestamp)\n \n#create unique messages\nindices<-which(dat$What==\"Message\")\ncounter=1\ndat$ID<-NA\nfor (x in 1:(length(indices)-1)){\ndat[indices[x]:indices[x+1],\"ID\"]<-counter\ncounter=counter+1\n}\n\ndive<-dat %>% filter(What==\"Dive\")%>% dplyr::select(Animal=Ptt,timestamp,Hour,Month,Year, ID,Start,DepthMax,DepthMin,DurationMax)\ndive<-dive[!is.na(dive$Month),]\n \n#remove duplicate data\ndive<-dive %>% arrange(Animal,timestamp) \ndive<-dive[!duplicated(dive),]\n```\n\n```{r}\n##Merge with geographic data and format for JAGS\nmessages<-dat %>% filter(dat$What==\"Message\")\nmessages$timestamp<-as.POSIXct(messages$End,format=\"%H:%M:%S %d-%b-%Y\",tz=\"Etc/GMT+3\")\n\n#look up the time interval that best fits\nsetDT(gdat) ## convert to data.table by reference\nsetDT(messages) ## same\n\nsetkey(messages, Ptt,timestamp) ## set the column to perform the join on\nsetkey(gdat,Ptt, Date) ## same as above\n\nans = gdat[messages, roll=Inf] ## perform rolling join\nans<-as.data.frame(ans)\n\nmessage_join<-ans %>% select(Date,Animal=Ptt,Date,Quality,Latitude,Longitude,ID)\n\nmdat<-merge(dat,message_join,by=\"ID\")\nmdat<-mdat %>% filter(What==\"Dive\")\n\n#The maxiumum dive depth based on the geographic message\nmdat<-mdat %>% select(ID,Animal,Latitude,Longitude,timestamp,Date,Month,Hour,Year,DepthMax,Quality) \n\nmdat<-mdat %>% group_by(ID,Animal) %>% slice(which.max(DepthMax)) %>% arrange(Animal,timestamp) %>% filter(!is.na(Latitude))\n\n#crop by extent\nd<-SpatialPointsDataFrame(cbind(mdat$Longitude,mdat$Latitude),data=data.frame(mdat),proj4string=CRS(\"+proj=longlat +datum=WGS84\"))\n\ncropoly<-readShapePoly(\"Data/CutPolygon.shp\",proj4string=CRS(\"+proj=longlat +datum=WGS84\"))\n\nb<-d[!is.na(d %over% cropoly)[,2],]\n\nmdat<-b@data\n\n#view data\nggplot(data=mdat)+geom_path(aes(x=Longitude, y=Latitude,group=Animal),size=.5) + geom_point(aes(x=Longitude, y=Latitude,col=DepthMax))+ borders(fill=\"grey90\") + coord_cartesian(ylim = range(mdat$Latitude),xlim=range(mdat$Longitude)) + theme_bw() + mytheme + scale_color_continuous(low=\"blue\",high=\"red\") + labs(col=\"Max Dive Depth(m)\")\n\n##Time is the beginning of the first point.\nstep_length=6\n\nsxy<-split(mdat,mdat$Animal)\n\n#time diff function\ntimed<-function(d,step_length){\n d$j[1]<-0\n for (x in 2:nrow(d)){\n d$j[x]<-as.numeric(difftime(as.POSIXct(d$timestamp[x]),as.POSIXct(d$timestamp[x-1]),units=\"mins\"))/(step_length*60)\n }\n \n #Split out track endings\n ends<-c(1,which(d$j>1),nrow(d))\n\n for(w in 2:length(ends)){\n d[ends[w-1]:ends[w],\"Track\"]<-w-1\n }\n \n #remove tracks that are shorter than three days\n track_time<-d %>% group_by(Track) %>% summarize(mt=difftime(max(as.POSIXct(timestamp)),min(as.POSIXct(timestamp)),units=\"hours\")) %>% filter(mt>=24) %>% .$Track\n \n d<-d[d$Track %in% track_time,]\n \n #renumber the tracks\n d$Track<-as.numeric(as.factor(d$Track))\n return(d)\n }\n\nsxy<-lapply(sxy,timed,step_length=step_length)\n\n#Format matrices for jags\nmdat<-bind_rows(sxy)\n\nsxy<-split(mdat,list(mdat$Animal,mdat$Track),drop=TRUE)\n\nsxy<-lapply(sxy,function(x){\n#How many observations in each step length segment\nx$step<-as.numeric(cut(as.POSIXct(x$timestamp),paste(step_length,\"hours\")))\nreturn(x)\n})\n\nmdat<-bind_rows(sxy)\n\n###################################################\n#filter by two whales for the moment\nmdat<-mdat %>% filter(Animal %in% c(\"131127\",\"131136\"))\n####################################################\n\n#refactor animal\nmdat$Animal<-as.numeric(as.factor(mdat$Animal))\n\n##Split into format for JAGS\n\n#total number of steps per track/animal\nsteps_all<-mdat %>% group_by(Animal,Track) %>% summarize(n=length(unique(step)))\n\n# give each step a label\nmdat<-mdat %>% group_by(Animal,Track,step) %>% mutate(jStep=1:n())\n\n#Cast time array\nj<-reshape2::acast(mdat,Animal~Track~step~jStep,value.var=\"j\")\n\n#how many observations per individual in each step\nmdat$step<-factor(mdat$step,levels=1:max(steps_all$n))\nidx<-reshape2::melt(table(mdat$Animal,mdat$Track,mdat$step))\ncolnames(idx)<-c(\"Animal\",\"Track\",\"step\",\"jStep\")\nidx<-reshape2::acast(data=idx,Animal~Track~step)\n\n#Individuals\nind=length(unique(mdat$Animal))\n\n#tracks per indivudal\ntracks<-mdat %>% group_by(Animal) %>% summarize(tracks=length(unique(Track))) %>% .$tracks\n\n#steps per track\nsteps<-reshape2::acast(steps_all,Animal~Track,value.var=\"n\")\n\n#obs array\nobs<-reshape2::melt(mdat,measure.vars=c(\"Longitude\",\"Latitude\"))\nobs<-reshape2::acast(obs,Animal~Track~step~jStep~variable,fun.aggregate = mean)\n\n#argos class array\nmdat$argos.lc<-factor(mdat$Quality,levels=c(3,2,1,0,\"A\",\"B\"))\nmdat$numargos<-as.numeric(mdat$argos.lc)\nobs_class<-reshape2::acast(mdat,Animal~Track~step~jStep,value.var=\"numargos\",fun.aggregate = min)\n\nobs_class[!is.finite(obs_class)]<-NA\n\n#average dive depth array\nmaxdive<-reshape2::acast(mdat,Animal~Track~step~jStep,value.var=\"DepthMax\",fun.aggregate = mean)\n\n#fill the empty values\nmaxdive[!is.finite(maxdive)]<-NA\n\n#average number of dives\ndivecount<-reshape2::acast(mdat,Animal~Track~step,value.var=\"DepthMax\",fun.aggregate = length)\n\n#no dives, fill empty\ndivecount[divecount==0]<-NA\n\n```\n\n```{r,child=\"Bayesian/Diving.R\",eval=T}\n```\n\n```{r,eval=T}\n#source jags file\nsource(\"Bayesian/Diving.R\")\n\n#prior cov shape\nR <- diag(c(1,1))\ndata=list(divecount=divecount,dive=maxdive,argos=obs,steps=steps,R=R,ind=ind,j=j,idx=idx,tracks=tracks,argos_class=obs_class)\n\n#paramters to track\npt<-c(\"alpha\",\"depth_mu\",\"depth_tau\",\"dive_new\",\"state\",\"lambda\")\n\nif(newModel){\n system.time(diving<-jags.parallel(model.file = \"Bayesian/Diving.jags\",data=data,n.chains=2,parameters.to.save=pt,n.iter=10000,n.burnin=9500,n.thin=2,DIC=FALSE))\n}\n\n```\n\n##Chains\n```{r,eval=T}\n#bind chains\npc_dive<-reshape2::melt(diving$BUGSoutput$sims.array)\ncolnames(pc_dive)<-c(\"Draw\",\"chain\",\"par\",\"value\")\n\n#extract parameter name\npc_dive$parameter<-data.frame(str_match(pc_dive$par,\"(\\\\w+)\"))[,-1]\n\n#Extract index\nsplitpc<-split(pc_dive,pc_dive$parameter)\n\n#single index\nsplitpc[c(\"alpha\",\"depth_mu\",\"depth_tau\",\"lambda\")]<-lapply(splitpc[c(\"alpha\",\"depth_mu\",\"depth_tau\",\"lambda\")],function(x){\n sv<-data.frame(str_match(x$par,\"(\\\\w+)\\\\[(\\\\d+)]\"))[,3]\n pc<-data.frame(x,Behavior=sv)\n return(pc)\n}) \n\n#3 index\nsplitpc[c(\"state\")]<-lapply(splitpc[c(\"state\")],function(x){\n sv<-data.frame(str_match(x$par,\"(\\\\w+)\\\\[(\\\\d+),(\\\\d+),(\\\\d+)\"))[,3:5]\n colnames(sv)<-c(\"Animal\",\"Track\",\"step\")\n pc<-data.frame(x,sv)\n return(pc)\n}) \n\n#4 index\nsplitpc[c(\"dive_new\")]<-lapply(splitpc[c(\"dive_new\")],function(x){\n sv<-data.frame(str_match(x$par,\"(\\\\w+)\\\\[(\\\\d+),(\\\\d+),(\\\\d+),(\\\\d+)]\"))[,3:6]\n colnames(sv)<-c(\"Animal\",\"Track\",\"step\",\"jStep\")\n pc<-data.frame(x,sv)\n return(pc)\n}) \n#bind all matrices back together\npc_dive<-bind_rows(splitpc)\nrm(splitpc)\n\nggplot(pc_dive[!pc_dive$parameter %in% c(\"dive_new\",\"state\"),],aes(x=Draw,y=value,col=as.factor(chain))) + geom_line() + facet_wrap(~par,scales=\"free\")\n\n```\n\n## Empirical distribution of foraging \n\n```{r}\nalpha<-pc_dive %>% filter(parameter==\"alpha\")\nalpha$Behavior<-as.character(alpha$Behavior)\nalpha$Behavior[alpha$Behavior==\"1\"]<-\"Traveling->Traveling\"\nalpha$Behavior[alpha$Behavior==\"2\"]<-\"Foraging->Traveling\"\n\nggplot(alpha,aes(x=value)) + geom_histogram() + facet_wrap(~Behavior) + labs(x=\"Probability of Transition\")\n\n#label previous and prediction state, treat traveling as a success and foraging as a failure in a binomial trial\nalpha[alpha$Behavior==\"Traveling->Traveling\",\"Previous_State\"]<-1\nalpha[alpha$Behavior==\"Foraging->Traveling\",\"Previous_State\"]<-0\n```\n\n### Time foraging as a function of time\n\n* 12 time step\n\n```{r}\nforaging_time<-function(alpha,draws,step){\n states<-c()\n \n #initial state is traveling\n states[1]<-(1)\n \n for(x in 2:draws){\n \n #probability of staying, or switching, to traveling state\n travel_prob=alpha %>% filter(Previous_State==states[x-1]) %>% sample_n(1) %>% .$value\n \n #pick next state\n states[x]<-rbinom(1,1,travel_prob)\n }\n \n #total hours of foraging\n total_foraging<-sum(states==0) * step\nreturn(total_foraging)\n}\n\nstates<-foraging_time(alpha=alpha,draw=20,step=12)\n```\n\nThis allows us to specify the hours foraging while capturing our uncertainty in behavioral states\n\nFor example, if we wanted to know how many hours an animal was predicted to forage on average over a 10 day span. We can draw 100 simulations given the posterior distributions of our empirical model.\n\n```{r}\n\ntest_time<-sapply(1:1000,function(x){\n foraging_time(alpha=alpha,draw=20,step=12)\n})\npaste(\"Mean foraging hours:\", mean(test_time),sep=\" \")\n\nqplot(test_time) +ggtitle(\"Foraging Time (10 Days)\") +xlab(\"Hours\")\n```\n\n\n##Distribution of average max dive depths\n\n```{r}\ndivep<-pc_dive %>% filter(parameter %in% c(\"depth_mu\",\"depth_tau\"))\n\n## Simulate dive depth\ndepth_mu1<-pc_dive %>% filter(parameter %in% c(\"depth_mu\"),Behavior==1)\ndepth_mu2<-pc_dive %>% filter(parameter %in% c(\"depth_mu\"),Behavior==2)\ndepth_tau1<-pc_dive %>% filter(parameter %in% c(\"depth_tau\"),Behavior==1)\ndepth_tau2<-pc_dive %>% filter(parameter %in% c(\"depth_tau\"),Behavior==2)\n\n#how many of each state to draw?\nstatecount<-pc_dive %>% filter(parameter==\"state\") %>% group_by(Animal,Track,step) %>% summarize(state=names(sort(table(value)))[1]) %>% group_by(state) %>% summarize(n=n()) %>% data.frame()\n\ntravel_dives<-data.frame(MaxDepth=rnorm(statecount[1,\"n\"],depth_mu1$value,sqrt(1/depth_tau1$value)),Behavior=\"Traveling\")\nforage_dives<-data.frame(MaxDepth=rnorm(statecount[2,\"n\"],depth_mu2$value,sqrt(1/depth_tau2$value)),Behavior=\"Foraging\")\npred_dives<-bind_rows(travel_dives,forage_dives)\nggplot(pred_dives) + geom_density(alpha=0.5,aes(x=MaxDepth,fill=Behavior)) + geom_density(data=mdat,aes(x=DepthMax),col=\"black\",size=1) + theme_bw() \n\nggplot(pc_dive[pc_dive$parameter==\"dive_new\",])+ geom_density(aes(x=value),col=\"red\") + geom_density(data=mdat,aes(x=DepthMax),col=\"black\",size=1.5) + theme_bw() \n```\n\n##Distribution of average dive counts\n\n```{r}\n\n## Simulate dive counts\ncount_lambda1<-pc_dive %>% filter(parameter %in% c(\"lambda\"),Behavior==1)\ncount_lambda2<-pc_dive %>% filter(parameter %in% c(\"lambda\"),Behavior==2)\n\nstatecount<-pc_dive %>% filter(parameter==\"state\") %>% group_by(Animal,Track,step) %>% summarize(state=names(sort(table(value)))[1]) %>% group_by(state) %>% summarize(n=n()) %>% data.frame()\n\ntravel_counts<-data.frame(Counts=rpois(statecount[1,\"n\"],count_lambda1$value),Behavior=\"Traveling\")\nforage_counts<-data.frame(Counts=rpois(statecount[2,\"n\"],count_lambda2$value),Behavior=\"Foraging\")\n\n#calculate observed\nobs_count<-mdat %>% group_by(Animal,Track,step) %>% summarize(Counts=n())\n\npred_counts<-bind_rows(travel_counts,forage_counts)\nggplot(pred_counts) + geom_density(alpha=0.5,aes(x=Counts,fill=Behavior)) + geom_density(data=obs_count,aes(x=Counts),col=\"black\",size=1) + theme_bw() \n```\n\n## Posterior checks\n\n* Compare observed dive depths by inferred behavior against generated data\n\n```{r,eval=F}\npred_dives<-pc_dive %>% filter(parameter==\"dive_new\") %>% group_by(Animal,Track,step,jStep) %>% summarize(Pred_DepthMax=mean(value)) %>% ungroup() %>% mutate(Animal=as.numeric(as.character(Animal)),Track=as.numeric(as.character(Track)),jStep=as.numeric(as.character(jStep)))\n\npredframe<-pred_dives %>% inner_join(mdat,by=c(\"Animal\",\"Track\",\"step\",\"jStep\")) %>% select(Animal,Latitude,Longitude,Track,step,jStep,timestamp,Predicted=Pred_DepthMax,Observed=DepthMax)\n\n#merge with state\npredframe<-pc_dive %>% filter(parameter==\"state\") %>% group_by(Animal,Track,step) %>% summarize(state=names(sort(table(value)))[1]) %>% ungroup() %>% mutate(Animal=as.numeric(as.character(Animal)),Track=as.numeric(as.character(Track))) %>% inner_join(predframe,by=c(\"Animal\",\"Track\",\"step\"))\n\npredframe<-reshape2::melt(predframe, measure.vars=c(\"Predicted\",\"Observed\"),value.name=\"DepthMax\")\n\nggplot(data=predframe) + geom_histogram(aes(x=DepthMax,fill=variable),alpha=0.8) + theme_bw()\n\nggplot(data=predframe) + geom_density(aes(x=DepthMax,fill=variable),alpha=0.8) \n\nggplot(data=predframe) + geom_density(aes(x=DepthMax,fill=variable),alpha=0.8) + facet_wrap(~state) + ggtitle(\"By Behavior\")\n\nggplot(data=predframe) + geom_density(aes(x=DepthMax,fill=variable),alpha=0.8) + facet_wrap(~Animal) + ggtitle(\"By Animal\")\n\n```\n\n```{r}\nsave.image(\"WhalePhys.RData\")\n```\n", "created" : 1496200015945.000, "dirty" : false, "encoding" : "UTF-8", "folds" : "", - "hash" : "2874219066", + "hash" : "1129738221", "id" : "449C785F", - "lastKnownWriteTime" : 1496785911, - "last_content_update" : 1496785911480, + "lastKnownWriteTime" : 1496803268, + "last_content_update" : 1496803268262, "path" : "~/WhalePhys/RunModel.Rmd", "project_path" : "RunModel.Rmd", "properties" : { diff --git a/.gitignore b/.gitignore index 45225f3..0769a39 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ *.RData +.Rproj.user diff --git a/RunModel.Rmd b/RunModel.Rmd index d2e7df5..153e872 100644 --- a/RunModel.Rmd +++ b/RunModel.Rmd @@ -60,7 +60,8 @@ library(chron) library(R2jags) library(knitr) -newModel=T + +newModel=F opts_chunk$set(echo=F,warning=F,message=F,fig.width=11) @@ -68,6 +69,11 @@ mytheme<-theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.tex ``` ```{r} + +if(!newModel){ + load("WhalePhys.RData") +} + #get dive data f<-list.files("C:/Users/Ben/Dropbox/Whales/Data/Humpback",pattern="Locations",full.names=T,recursive = T) @@ -385,12 +391,15 @@ depth_mu2<-pc_dive %>% filter(parameter %in% c("depth_mu"),Behavior==2) depth_tau1<-pc_dive %>% filter(parameter %in% c("depth_tau"),Behavior==1) depth_tau2<-pc_dive %>% filter(parameter %in% c("depth_tau"),Behavior==2) -travel_dives<-data.frame(MaxDepth=rnorm(1e5,depth_mu1$value,sqrt(1/depth_tau1$value)),Behavior="Traveling") -forage_dives<-data.frame(MaxDepth=rnorm(1e5,depth_mu2$value,sqrt(1/depth_tau2$value)),Behavior="Foraging") +#how many of each state to draw? +statecount<-pc_dive %>% filter(parameter=="state") %>% group_by(Animal,Track,step) %>% summarize(state=names(sort(table(value)))[1]) %>% group_by(state) %>% summarize(n=n()) %>% data.frame() + +travel_dives<-data.frame(MaxDepth=rnorm(statecount[1,"n"],depth_mu1$value,sqrt(1/depth_tau1$value)),Behavior="Traveling") +forage_dives<-data.frame(MaxDepth=rnorm(statecount[2,"n"],depth_mu2$value,sqrt(1/depth_tau2$value)),Behavior="Foraging") pred_dives<-bind_rows(travel_dives,forage_dives) ggplot(pred_dives) + geom_density(alpha=0.5,aes(x=MaxDepth,fill=Behavior)) + geom_density(data=mdat,aes(x=DepthMax),col="black",size=1) + theme_bw() -ggplot(pc_dive[pc_dive$parameter=="dive_new",])+ geom_density(aes(x=value,fill=par)) + geom_density(data=mdat,aes(x=DepthMax),col="black",size=1.5) + theme_bw() +ggplot(pc_dive[pc_dive$parameter=="dive_new",])+ geom_density(aes(x=value),col="red") + geom_density(data=mdat,aes(x=DepthMax),col="black",size=1.5) + theme_bw() ``` ##Distribution of average dive counts @@ -401,8 +410,10 @@ ggplot(pc_dive[pc_dive$parameter=="dive_new",])+ geom_density(aes(x=value,fill= count_lambda1<-pc_dive %>% filter(parameter %in% c("lambda"),Behavior==1) count_lambda2<-pc_dive %>% filter(parameter %in% c("lambda"),Behavior==2) -travel_counts<-data.frame(Counts=rpois(1e5,count_lambda1$value),Behavior="Traveling") -forage_counts<-data.frame(Counts=rpois(1e5,count_lambda2$value),Behavior="Foraging") +statecount<-pc_dive %>% filter(parameter=="state") %>% group_by(Animal,Track,step) %>% summarize(state=names(sort(table(value)))[1]) %>% group_by(state) %>% summarize(n=n()) %>% data.frame() + +travel_counts<-data.frame(Counts=rpois(statecount[1,"n"],count_lambda1$value),Behavior="Traveling") +forage_counts<-data.frame(Counts=rpois(statecount[2,"n"],count_lambda2$value),Behavior="Foraging") #calculate observed obs_count<-mdat %>% group_by(Animal,Track,step) %>% summarize(Counts=n()) diff --git a/RunModel.html b/RunModel.html index b3da0f6..3a59950 100644 --- a/RunModel.html +++ b/RunModel.html @@ -127,16 +127,13 @@

May 30, 2017