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.1Next 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
## 1 A-098 Ningaloo Marine Park (Coral Bay) 2003 11
## 2 A-052 Exmouth Gulf, Back of Reef in Bundegi Sanctuary 2005 11
## 3 A-013 North Ningaloo 2014 11
## 4 A-003 north ningaloo marine park 2011 5
## 5 A-079 North Ningaloo 2013 6
## 6 A-063 coral bay ningaloo marine park 2011 7
## day
## 1 10
## 2 13
## 3 16
## 4 13
## 5 6
## 6 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)))## Fitting model## 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.513276
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.5087466
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.0414and 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"))