Linking Wildbook and marked via RWildbook
Here we show how data retrieved from a Wildbook database can be formatted and analysed using the marked package by Jeff Laake (also the creator of the RMark package). In particular, we fit a simple Cormack-Jolly-Seber model to data retrieved from whaleshark.org.
First we load the RWildbook and marked package:
## Loading required package: jsonlite## Loading required package: data.table## Loading required package: lme4## Loading required package: Matrix## Loading required package: parallel## This is marked 1.2.8Next we use the searchWB() function from RWildbook to retrieve data for the first 99 whale sharks:
## Extract data for individual A-001 through A-099
vignette_2_data <- searchWB(username="username",
password="password",
baseURL ="whaleshark.org",
object="Encounter",
individualID=paste0("A-0",rep(0:9,rep(10,10)),rep(0:9,10))[-1])This results in the following data frame (not the that the time of encounter has been modified to maintain data privacy):
data("vignette_2_data",package="RWildbook")
head(vignette_2_data)[,c("individualID","location","year","month","day")]## individualID location year month day
## 1 A-098 Ningaloo Marine Park (Coral Bay) 2003 11 10
## 2 A-052 Exmouth Gulf, Back of Reef in Bundegi Sanctuary 2005 11 13
## 3 A-013 North Ningaloo 2014 11 16
## 4 A-003 north ningaloo marine park 2011 5 13
## 5 A-079 North Ningaloo 2013 6 6
## 6 A-063 coral bay ningaloo marine park 2011 7 7Now we can create capture histories by defining the start and end times of the capture occasions. For illustration, we assume that the capture occasions cover January, February, and March of each year from 1998 through 2016.
## Define start and end dates of capture occasions
start.dates1 <- paste0(1998:2016,"-01-01") #Define the start.date value
end.dates1 <- paste0(1998:2016,"-04-01") #Define the end.date value
## Format data for use in marked
markedData1 <- markedData(data = vignette_2_data,
varname_of_capturetime = "dateInMilliseconds",
varlist = c("individualID"),
start.dates = start.dates1,
end.dates = NULL,
date_format = "%Y-%m-%d",
origin = "1970-01-01")## Note: Removing 25 individuals with no observed captures.We can now fit the Cormack-Jolly-Seber model using the functions in mark to first format the data:
## Fit simple CJS model in marked
markedData1.proc=process.data(markedData1,model="cjs",begin.time=1)## 74 capture histories collapsed into 74 markedData1.ddl=make.design.data(markedData1.proc)
markedData1.cjs=crm(markedData1.proc,markedData1.ddl,model.parameters=list(Phi=list(formula=~time),p=list(formula=~time)))## Computing initial parameter estimates## Starting optimization for 36 parameters...## Number of evaluations: 100 -2lnl: 661.4810565 Number of evaluations: 200 -2lnl: 657.6593186 Number of evaluations: 300 -2lnl: 654.738069 Number of evaluations: 400 -2lnl: 654.2391281 Number of evaluations: 500 -2lnl: 653.862695 Number of evaluations: 600 -2lnl: 653.4780271 Number of evaluations: 700 -2lnl: 653.3505049 Number of evaluations: 800 -2lnl: 653.126154 Number of evaluations: 900 -2lnl: 652.916061 Number of evaluations: 1000 -2lnl: 652.85315 Number of evaluations: 1100 -2lnl: 652.7996062 Number of evaluations: 1200 -2lnl: 652.7025896 Number of evaluations: 1300 -2lnl: 652.668327 Number of evaluations: 1400 -2lnl: 652.5942031 Number of evaluations: 1500 -2lnl: 652.5536496 Number of evaluations: 1600 -2lnl: 652.5443796 Number of evaluations: 1700 -2lnl: 652.5235439 Number of evaluations: 1800 -2lnl: 652.5164447 Number of evaluations: 1900 -2lnl: 652.5126744 Number of evaluations: 2000 -2lnl: 652.5075915 Number of evaluations: 2100 -2lnl: 652.5068583 Number of evaluations: 2200 -2lnl: 652.5059795 Number of evaluations: 2300 -2lnl: 652.502895 Number of evaluations: 2400 -2lnl: 652.5016548 Number of evaluations: 2500 -2lnl: 652.4991583 Number of evaluations: 2600 -2lnl: 652.498644 Number of evaluations: 2700 -2lnl: 652.4983563 Number of evaluations: 2800 -2lnl: 652.4975948 Number of evaluations: 2900 -2lnl: 652.4968093 Number of evaluations: 3000 -2lnl: 652.4957483 Number of evaluations: 3100 -2lnl: 652.49337 Number of evaluations: 3200 -2lnl: 652.4913842 Number of evaluations: 3300 -2lnl: 652.4884967 Number of evaluations: 3400 -2lnl: 652.4840881 Number of evaluations: 3500 -2lnl: 652.4834014 Number of evaluations: 3600 -2lnl: 652.4825869 Number of evaluations: 3700 -2lnl: 652.4822268 Number of evaluations: 3800 -2lnl: 652.4820066 Number of evaluations: 3900 -2lnl: 652.4819051 Number of evaluations: 4000 -2lnl: 652.4818957 Number of evaluations: 4100 -2lnl: 652.4818741 Number of evaluations: 4200 -2lnl: 652.4818424 Number of evaluations: 4300 -2lnl: 652.5423162 Number of evaluations: 4400 -2lnl: 652.483646 Number of evaluations: 4500 -2lnl: 652.4878773 Number of evaluations: 4600 -2lnl: 652.5380763 Number of evaluations: 4700 -2lnl: 652.4848939 Number of evaluations: 4800 -2lnl: 652.4858153 Number of evaluations: 4900 -2lnl: 652.4839405 Number of evaluations: 5000 -2lnl: 652.4822173 Number of evaluations: 5100 -2lnl: 652.5428697 Number of evaluations: 5200 -2lnl: 652.4858546 Number of evaluations: 5300 -2lnl: 652.5422815 Number of evaluations: 5400 -2lnl: 652.4821965 Number of evaluations: 5500 -2lnl: 652.4819787 Number of evaluations: 5600 -2lnl: 652.482424 Number of evaluations: 5700 -2lnl: 652.5132759 Number of evaluations: 5800 -2lnl: 652.4872118 Number of evaluations: 5900 -2lnl: 652.5201534 Number of evaluations: 6000 -2lnl: 652.4848467 Number of evaluations: 6100 -2lnl: 652.4901234 Number of evaluations: 6200 -2lnl: 652.4827419 Number of evaluations: 6300 -2lnl: 652.4845256 Number of evaluations: 6400 -2lnl: 652.5361715 Number of evaluations: 6500 -2lnl: 652.4900906 Number of evaluations: 6600 -2lnl: 652.4857628 Number of evaluations: 6700 -2lnl: 652.4918907 Number of evaluations: 6800 -2lnl: 652.4857677 Number of evaluations: 6900 -2lnl: 652.5087465 Number of evaluations: 7000 -2lnl: 652.481842 Number of evaluations: 7100 -2lnl: 652.5200702 Number of evaluations: 7200 -2lnl: 652.481838 Number of evaluations: 7300 -2lnl: 652.5200059 Number of evaluations: 7400 -2lnl: 652.4856887 Number of evaluations: 7500 -2lnl: 652.5075832 Number of evaluations: 7600 -2lnl: 652.4864769 Number of evaluations: 7700 -2lnl: 652.5031411 Number of evaluations: 7800 -2lnl: 652.5360589 Number of evaluations: 7900 -2lnl: 652.4926024 Number of evaluations: 8000 -2lnl: 652.4823632 Number of evaluations: 8100 -2lnl: 652.4878746 Number of evaluations: 8200 -2lnl: 652.4823879 Number of evaluations: 8300 -2lnl: 652.4931355 Number of evaluations: 8400 -2lnl: 652.4829153 Number of evaluations: 8500 -2lnl: 652.4982733 Number of evaluations: 8600 -2lnl: 652.4828685 Number of evaluations: 8700 -2lnl: 652.489097 Number of evaluations: 8800 -2lnl: 652.4856823 Number of evaluations: 8900 -2lnl: 652.4906435 Number of evaluations: 9000 -2lnl: 652.4826364 Number of evaluations: 9100 -2lnl: 652.4847457 Number of evaluations: 9200 -2lnl: 652.4822262 Number of evaluations: 9300 -2lnl: 652.4826786 Number of evaluations: 9400 -2lnl: 652.4818565 Number of evaluations: 9500 -2lnl: 652.4821469 Number of evaluations: 9600 -2lnl: 652.4818359 Number of evaluations: 9700 -2lnl: 652.4818359 Number of evaluations: 9800 -2lnl: 652.4818359##
## Elapsed time in minutes: 0.0307and plot the estimated capture and survival probabilities:
## Plot parameter estimates
plot(1998:2015,markedData1.cjs$results$reals$Phi$estimate,
pch=16,col="green",
main="Estimated Capture and Survival Probabilities",
xlab="Occasion",ylab="Estimate",ylim=c(0,1))
points(1999:2016,markedData1.cjs$results$reals$p$estimate,pch=16,col="blue")
legend("bottomright",pch=16,col=c("green","blue"),legend=c("Survival","Capture"))