class: center, middle, title-slide # Abundance - live demos ### Aurélien Besnard for the team ### last updated: 2022-03-23 --- class: center, middle background-size: cover # Live demo #1: Capture-recapture for closed populations <style type="text/css">pre {font-size: 18px}</style> --- # Example with a single group ## Load package ```r library(RMark) ``` --- ## Read and prepare the data ```r mouse <- convert.inp("dat/mouse.inp", group.df = NULL, covariates = NULL) ``` Inspect first lines of the dataset ```r head(mouse) ``` ``` ch freq 1 1000 1 2 0011 1 3 0010 1 4 0001 1 5 1001 1 6 0011 1 ``` ```r dim(mouse)[1] # number of individuals captured at least once ``` ``` [1] 82 ``` --- ## Fit a first batch of models To use RMark, we need three steps : - Prepare the data - Define the models to fit - Fit the models --- ## Prepare the data ```r mouse.proc <- process.data(mouse, begin.time = 1, model = "FullHet") mouse.ddl <- make.design.data(mouse.proc) ``` `model = FullHet` refers to the most complex model form, allowing simpler models to be fitted as specific cases --- ## Define models (1) Use a function that do three things: - Specify the effects in the model to be fitted - Create a list of all models to be fitted - Prepare them for Mark Note that RMark considers capture probability `\(p\)` and recapture probability `\(c\)` to be different To fit a model vith `\(p\)` and `\(c\)` equal, use `share = TRUE` --- ## Define models (2) ```r run.mouse <- function() { # List of effects p.dot <- list(formula = ~ 1, share = TRUE) # M0 : p constant over time p.dot.behav <- list(formula = ~ 1) # Mb : p and c differ p.time <- list(formula = ~ time, share = TRUE) # Mt : p varies over time p.h <- list(formula = ~ mixture, share = TRUE) # Mh : p is heterogeneous between individuals (Pledger's mixture model) p.time.behav <- list(formula = ~ time) # Mtb p.h.behav <- list(formula = ~ mixture) # Mbh p.h.time <- list(formula = ~ time + mixture, share = TRUE) # Mth p.h.time.behav <- list(formula = ~ mixture + time) # Mtbh ## Build model list mouse.model.list <- create.model.list("FullHet") ## Prepare list to send to Mark mouse.results <- mark.wrapper(mouse.model.list,data = mouse.proc,ddl = mouse.ddl) ## Grab results return(mouse.results) } ``` --- # Fit models Fit at once all models prepared with the function above ```r mouse.results <- run.mouse() ``` ``` Output summary for FullHet model Name : pi(~1)p(~1)c()f0(~1) Npar : 3 (unadjusted=2) -2lnL: -181.4414 AICc : -175.3674 (unadjusted=-177.40452) Beta estimate se lcl ucl pi:(Intercept) 3.110286e-05 0.0000000 3.110286e-05 3.110286e-05 p:(Intercept) -1.236656e+00 0.1811077 -1.591627e+00 -8.816851e-01 f0:(Intercept) 3.823569e+00 0.3178170 3.200647e+00 4.446490e+00 Real Parameter pi mixture:1 0.5000078 Real Parameter p 1 2 3 4 mixture:1 0.2250185 0.2250185 0.2250185 0.2250185 mixture:2 0.2250185 0.2250185 0.2250185 0.2250185 Real Parameter c 2 3 4 mixture:1 0.2250185 0.2250185 0.2250185 mixture:2 0.2250185 0.2250185 0.2250185 Real Parameter f0 1 45.76725 Output summary for FullHet model Name : pi(~1)p(~1)c(~1)f0(~1) Npar : 4 (unadjusted=3) -2lnL: -189.5136 AICc : -181.3897 (unadjusted=-183.43951) Beta estimate se lcl ucl pi:(Intercept) 1.135781e-05 575.4448600 -1127.8719000 1127.8720000 p:(Intercept) -3.404630e-01 0.2668737 -0.8635354 0.1826095 c:(Intercept) -1.437966e+00 0.1936421 -1.8175049 -1.0584278 f0:(Intercept) 2.333576e+00 0.6268405 1.1049689 3.5621836 Real Parameter pi mixture:1 0.5000028 Real Parameter p 1 2 3 4 mixture:1 0.415697 0.415697 0.415697 0.415697 mixture:2 0.415697 0.415697 0.415697 0.415697 Real Parameter c 2 3 4 mixture:1 0.1918605 0.1918605 0.1918605 mixture:2 0.1918605 0.1918605 0.1918605 Real Parameter f0 1 10.31476 Output summary for FullHet model Name : pi(~1)p(~mixture)c()f0(~1) Npar : 4 (unadjusted=3) -2lnL: -181.4414 AICc : -173.3176 (unadjusted=-175.36737) Beta estimate se lcl ucl pi:(Intercept) -4.204911e+00 0.0000000 -4.204911 -4.204911 p:(Intercept) -1.236690e+00 4.2986234 -9.661992 7.188612 p:mixture2 3.492275e-05 4.3589156 -8.543440 8.543510 f0:(Intercept) 3.823568e+00 0.3178171 3.200646 4.446489 Real Parameter pi mixture:1 0.0147027 Real Parameter p 1 2 3 4 mixture:1 0.2250126 0.2250126 0.2250126 0.2250126 mixture:2 0.2250187 0.2250187 0.2250187 0.2250187 Real Parameter c 2 3 4 mixture:1 0.2250126 0.2250126 0.2250126 mixture:2 0.2250187 0.2250187 0.2250187 Real Parameter f0 1 45.76721 Output summary for FullHet model Name : pi(~1)p(~mixture)c(~1)f0(~1) Npar : 5 (unadjusted=4) -2lnL: -189.5136 AICc : -179.3272 (unadjusted=-181.38975) Beta estimate se lcl ucl pi:(Intercept) 2.204661e-01 748.1448200 -1466.143400 1466.584300 p:(Intercept) -3.405065e-01 4.3084070 -8.784984 8.103971 p:mixture2 9.660131e-05 9.6601463 -18.933790 18.933984 c:(Intercept) -1.437966e+00 0.1936421 -1.817505 -1.058428 f0:(Intercept) 2.333577e+00 0.6268410 1.104969 3.562186 Real Parameter pi mixture:1 0.5548943 Real Parameter p 1 2 3 4 mixture:1 0.4156864 0.4156864 0.4156864 0.4156864 mixture:2 0.4157099 0.4157099 0.4157099 0.4157099 Real Parameter c 2 3 4 mixture:1 0.1918605 0.1918605 0.1918605 mixture:2 0.1918605 0.1918605 0.1918605 Real Parameter f0 1 10.31477 Output summary for FullHet model Name : pi(~1)p(~time + mixture)c()f0(~1) Npar : 7 (unadjusted=6) -2lnL: -190.6236 AICc : -176.2736 (unadjusted=-178.36193) Beta estimate se lcl ucl pi:(Intercept) -4.288145e+00 9379.1667000 -18387.455000 18378.8790000 p:(Intercept) -7.217527e-01 4.6257544 -9.788232 8.3447261 p:time2 -8.256117e-01 0.3038091 -1.421077 -0.2301459 p:time3 -6.683295e-01 0.2947945 -1.246127 -0.0905322 p:time4 -5.712429e-01 0.2898297 -1.139309 -0.0031767 p:mixture2 5.124827e-07 4.6822722 -9.177253 9.1772542 f0:(Intercept) 3.769986e+00 0.3229069 3.137089 4.4028841 Real Parameter pi mixture:1 0.0135444 Real Parameter p 1 2 3 4 mixture:1 0.3270071 0.1754673 0.1993946 0.2153462 mixture:2 0.3270072 0.1754673 0.1993947 0.2153463 Real Parameter c 2 3 4 mixture:1 0.1754673 0.1993946 0.2153462 mixture:2 0.1754673 0.1993947 0.2153463 Real Parameter f0 1 43.37948 Output summary for FullHet model Name : pi(~1)p(~mixture + time)c(~1)f0(~1) Npar : 8 (unadjusted=4) -2lnL: -195.6955 AICc : -179.2441 (unadjusted=-187.57165) Beta estimate se lcl ucl pi:(Intercept) -0.4231560 0.0000000 -0.4231560 -0.4231560 p:(Intercept) 3.0166405 575.8003600 -1125.5521000 1131.5854000 p:mixture2 -4.3837581 885.6926900 -1740.3415000 1731.5740000 p:time2 0.9319872 0.0000000 0.9319872 0.9319872 p:time3 1.8760688 0.0000000 1.8760688 1.8760688 p:time4 20.0474640 3540.1957000 -6918.7363000 6958.8312000 c:(Intercept) -1.4379663 0.1936421 -1.8175049 -1.0584278 f0:(Intercept) -18.7515600 3524.2770000 -6926.3346000 6888.8315000 Real Parameter pi mixture:1 0.3957618 Real Parameter p 1 2 3 4 mixture:1 0.9533203 0.9810836 0.9925548 1 mixture:2 0.2030859 0.3929019 0.6245606 1 Real Parameter c 2 3 4 mixture:1 0.1918605 0.1918605 0.1918605 mixture:2 0.1918605 0.1918605 0.1918605 Real Parameter f0 1 7.182917e-09 Output summary for FullHet model Name : pi(~1)p(~time)c()f0(~1) Npar : 6 -2lnL: -190.6236 AICc : -178.3619 Beta estimate se lcl ucl pi:(Intercept) -3.328945e-05 1773.6801000 -3476.413000 3476.4129000 p:(Intercept) -7.217525e-01 0.2525860 -1.216821 -0.2266840 p:time2 -8.256115e-01 0.3038092 -1.421078 -0.2301454 p:time3 -6.683295e-01 0.2947947 -1.246127 -0.0905318 p:time4 -5.712428e-01 0.2898300 -1.139310 -0.0031760 f0:(Intercept) 3.769987e+00 0.3229070 3.137089 4.4028844 Real Parameter pi mixture:1 0.4999917 Real Parameter p 1 2 3 4 mixture:1 0.3270072 0.1754673 0.1993947 0.2153463 mixture:2 0.3270072 0.1754673 0.1993947 0.2153463 Real Parameter c 2 3 4 mixture:1 0.1754673 0.1993947 0.2153463 mixture:2 0.1754673 0.1993947 0.2153463 Real Parameter f0 1 43.37949 Output summary for FullHet model Name : pi(~1)p(~time)c(~1)f0(~1) Npar : 7 (unadjusted=4) -2lnL: -195.6955 AICc : -181.3455 (unadjusted=-187.57165) Beta estimate se lcl ucl pi:(Intercept) -1.126492e-04 1.448318e+03 -2.838703e+03 2838.7029000 p:(Intercept) -5.398854e-07 2.208630e-01 -4.328921e-01 0.4328910 p:time2 -3.448392e-01 3.863555e-01 -1.102096e+00 0.4124177 p:time3 5.108269e-01 4.759814e-01 -4.220965e-01 1.4437504 p:time4 2.130196e+01 1.100090e+04 -2.154047e+04 21583.0690000 c:(Intercept) -1.437966e+00 1.936421e-01 -1.817505e+00 -1.0584276 f0:(Intercept) -2.666935e+01 2.872209e+04 -5.632197e+04 56268.6320000 Real Parameter pi mixture:1 0.4999718 Real Parameter p 1 2 3 4 mixture:1 0.4999999 0.4146343 0.6250002 1 mixture:2 0.4999999 0.4146343 0.6250002 1 Real Parameter c 2 3 4 mixture:1 0.1918605 0.1918605 0.1918605 mixture:2 0.1918605 0.1918605 0.1918605 Real Parameter f0 1 2.61607e-12 ``` --- # Inspect outputs ```r mouse.results ``` Look at column `model` to get details of effects included --- # Inspect outputs ```r mouse.results ``` ``` model npar AICc DeltaAICc weight 2 pi(~1)p(~1)c(~1)f0(~1) 4 -181.3897 0.00000000 0.328826111 8 pi(~1)p(~time)c(~1)f0(~1) 7 -181.3455 0.04425099 0.321630567 4 pi(~1)p(~mixture)c(~1)f0(~1) 5 -179.3272 2.06249639 0.117246772 6 pi(~1)p(~mixture + time)c(~1)f0(~1) 8 -179.2441 2.14566165 0.112471319 7 pi(~1)p(~time)c()f0(~1) 6 -178.3619 3.02781323 0.072357742 5 pi(~1)p(~time + mixture)c()f0(~1) 7 -176.2736 5.11613099 0.025469040 1 pi(~1)p(~1)c()f0(~1) 3 -175.3674 6.02237506 0.016189154 3 pi(~1)p(~mixture)c()f0(~1) 4 -173.3176 8.07214000 0.005809294 Deviance 2 11.72881 8 5.54690 4 11.72881 6 5.54690 7 10.61878 5 10.61878 1 19.80095 3 19.80095 ``` --- # How do we identify models? The models appear in the same order as we specify them in the building model function: ```r names(mouse.results) ``` ``` [1] "p.dot" "p.dot.behav" "p.h" "p.h.behav" [5] "p.h.time" "p.h.time.behav" "p.time" "p.time.behav" [9] "model.table" ``` --- # Get estimates of a specific model We called `\(M_0\)` "p.dot" in the building model function This is the name used in the object storing the results ```r mouse.results$p.dot$results$real ``` ``` estimate se lcl ucl fixed note pi g1 m1 0.5000078 0.0000000 0.5000078 0.5000078 p g1 t1 m1 0.2250185 0.0315825 0.1691551 0.2928287 f0 g1 a0 t1 45.7672480 14.5456110 24.9175180 84.0629860 ``` + `pi` is not relevant here as the model does not include heterogeneity + `p` is the capture and recapture probability + `f0` is the estimated number of missed individuals (never captured) --- # How to get abundance estimates? Abundance `\(N\)` is a derived parameter estimated using 'f0' and the number of caught individuals To get its estimate: ```r mouse.results$p.dot$results$derived ``` ``` $`N Population Size` estimate lcl ucl 1 128 107 166 ``` --- # Model with two groups ## Read and convert the data Read the file, do not forget to grab the information on sex of individuals ```r mouse <- convert.inp("dat/mouse_groups.inp", group.df = data.frame(sex = c("M","F")), covariates = NULL) ``` --- # Look at the data ```r head(mouse) ``` ``` ch freq sex 1:1 0110 1 M 1:2 0101 1 M 1:3 1100 1 M 1:4 1100 1 M 1:5 1100 1 M 1:6 1010 1 M ``` --- # Prepare the data Same steps as before but specify that we have some groups, here 'sex' ```r mouse.proc <- process.data(mouse, begin.time = 1, model = "FullHet", groups = "sex") mouse.ddl <- make.design.data(mouse.proc) ``` --- ## Define models to be fitted ```r run.mouse <- function() { p.dot <- list(formula = ~ 1, share = TRUE) p.dot.behav <- list(formula = ~ 1) p.time <- list(formula = ~ time, share = TRUE) p.h <- list(formula = ~ mixture, share = TRUE) p.time.behav <- list(formula = ~ time) p.h.behav <- list(formula = ~ mixture) p.h.time <- list(formula = ~ time + mixture, share = TRUE) p.h.time.behav <- list(formula = ~ mixture + time) # Mbsex - model with sex on capture probability AND recapture probability p.sex.behav <- list(p = list(formula = ~ sex),c = list(formula = ~ sex)) ## Build model list mouse.model.list <- create.model.list("FullHet") ## Prepare list to send to Mark mouse.results <- mark.wrapper(mouse.model.list, data = mouse.proc, ddl = mouse.ddl) ## Grab results return(mouse.results) } ``` --- # Fit all models ```r mouse.results <- run.mouse() ``` ``` Output summary for FullHet model Name : pi(~1)p(~1)c()f0(~1) Npar : 3 (unadjusted=2) -2lnL: -55.76693 AICc : -49.70241 (unadjusted=-51.734755) Beta estimate se lcl ucl pi:(Intercept) 5.316308e-05 0.0000000 5.316308e-05 5.316308e-05 p:(Intercept) -6.141761e-01 0.1352648 -8.792952e-01 -3.490571e-01 f0:(Intercept) 2.262456e+00 0.3546695 1.567303e+00 2.957608e+00 Real Parameter pi Group:sexF mixture:1 0.5000133 Group:sexM mixture:1 0.5000133 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.3511072 0.3511072 0.3511072 0.3511072 mixture:2 0.3511072 0.3511072 0.3511072 0.3511072 Group:sexM 1 2 3 4 mixture:1 0.3511072 0.3511072 0.3511072 0.3511072 mixture:2 0.3511072 0.3511072 0.3511072 0.3511072 Real Parameter c Group:sexF 2 3 4 mixture:1 0.3511072 0.3511072 0.3511072 mixture:2 0.3511072 0.3511072 0.3511072 Group:sexM 2 3 4 mixture:1 0.3511072 0.3511072 0.3511072 mixture:2 0.3511072 0.3511072 0.3511072 Real Parameter f0 Group:sexF 1 9.606651 Group:sexM 1 9.606651 Output summary for FullHet model Name : pi(~1)p(~1)c(~1)f0(~1) Npar : 4 -2lnL: -64.6307 AICc : -56.52289 Beta estimate se lcl ucl pi:(Intercept) 4.472987e-05 465.7962700 -912.9606700 912.9607600 p:(Intercept) 5.183910e-02 0.2198099 -0.3789883 0.4826665 c:(Intercept) -8.092194e-01 0.1491105 -1.1014759 -0.5169629 f0:(Intercept) 8.247932e-01 0.7956631 -0.7347066 2.3842929 Real Parameter pi Group:sexF mixture:1 0.5000112 Group:sexM mixture:1 0.5000112 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.5129569 0.5129569 0.5129569 0.5129569 mixture:2 0.5129569 0.5129569 0.5129569 0.5129569 Group:sexM 1 2 3 4 mixture:1 0.5129569 0.5129569 0.5129569 0.5129569 mixture:2 0.5129569 0.5129569 0.5129569 0.5129569 Real Parameter c Group:sexF 2 3 4 mixture:1 0.3080569 0.3080569 0.3080569 mixture:2 0.3080569 0.3080569 0.3080569 Group:sexM 2 3 4 mixture:1 0.3080569 0.3080569 0.3080569 mixture:2 0.3080569 0.3080569 0.3080569 Real Parameter f0 Group:sexF 1 2.281409 Group:sexM 1 2.281409 Output summary for FullHet model Name : pi(~1)p(~mixture)c()f0(~1) Npar : 4 (unadjusted=3) -2lnL: -55.76693 AICc : -47.65911 (unadjusted=-49.70241) Beta estimate se lcl ucl pi:(Intercept) -1.679831e+00 991.0703400 -1944.177700 1940.818100 p:(Intercept) -6.141753e-01 1.2492707 -3.062746 1.834395 p:mixture2 -9.787118e-07 1.4734283 -2.887921 2.887919 f0:(Intercept) 2.262456e+00 0.3546694 1.567304 2.957608 Real Parameter pi Group:sexF mixture:1 0.1571178 Group:sexM mixture:1 0.1571178 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.3511073 0.3511073 0.3511073 0.3511073 mixture:2 0.3511071 0.3511071 0.3511071 0.3511071 Group:sexM 1 2 3 4 mixture:1 0.3511073 0.3511073 0.3511073 0.3511073 mixture:2 0.3511071 0.3511071 0.3511071 0.3511071 Real Parameter c Group:sexF 2 3 4 mixture:1 0.3511073 0.3511073 0.3511073 mixture:2 0.3511071 0.3511071 0.3511071 Group:sexM 2 3 4 mixture:1 0.3511073 0.3511073 0.3511073 mixture:2 0.3511071 0.3511071 0.3511071 Real Parameter f0 Group:sexF 1 9.606651 Group:sexM 1 9.606651 Output summary for FullHet model Name : pi(~1)p(~mixture)c(~1)f0(~1) Npar : 5 (unadjusted=4) -2lnL: -64.6307 AICc : -54.46854 (unadjusted=-56.522887) Beta estimate se lcl ucl pi:(Intercept) -4.280296e-01 1079.9895000 -2117.2076000 2116.3515000 p:(Intercept) 5.184060e-02 0.7289000 -1.3768035 1.4804847 p:mixture2 -2.098418e-06 1.1479404 -2.2499654 2.2499612 c:(Intercept) -8.092194e-01 0.1491105 -1.1014759 -0.5169629 f0:(Intercept) 8.247922e-01 0.7956630 -0.7347072 2.3842917 Real Parameter pi Group:sexF mixture:1 0.3945969 Group:sexM mixture:1 0.3945969 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.5129572 0.5129572 0.5129572 0.5129572 mixture:2 0.5129567 0.5129567 0.5129567 0.5129567 Group:sexM 1 2 3 4 mixture:1 0.5129572 0.5129572 0.5129572 0.5129572 mixture:2 0.5129567 0.5129567 0.5129567 0.5129567 Real Parameter c Group:sexF 2 3 4 mixture:1 0.3080569 0.3080569 0.3080569 mixture:2 0.3080569 0.3080569 0.3080569 Group:sexM 2 3 4 mixture:1 0.3080569 0.3080569 0.3080569 mixture:2 0.3080569 0.3080569 0.3080569 Real Parameter f0 Group:sexF 1 2.281407 Group:sexM 1 2.281407 Output summary for FullHet model Name : pi(~1)p(~time + mixture)c()f0(~1) Npar : 7 (unadjusted=6) -2lnL: -62.91561 AICc : -48.61126 (unadjusted=-50.687965) Beta estimate se lcl ucl pi:(Intercept) 1.902482e-01 0.0000000 0.1902482 0.1902482 p:(Intercept) -1.851382e-01 0.5629945 -1.2886073 0.9183310 p:time2 -4.468512e-01 0.2747448 -0.9853509 0.0916486 p:time3 -6.925080e-01 0.2817889 -1.2448142 -0.1402019 p:time4 -5.669672e-01 0.2778927 -1.1116369 -0.0222976 p:mixture2 3.101711e-06 1.1467735 -2.2476730 2.2476792 f0:(Intercept) 2.217709e+00 0.3611861 1.5097840 2.9256336 Real Parameter pi Group:sexF mixture:1 0.5474191 Group:sexM mixture:1 0.5474191 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.4538472 0.3470596 0.2936658 0.3203627 mixture:2 0.4538480 0.3470603 0.2936664 0.3203634 Group:sexM 1 2 3 4 mixture:1 0.4538472 0.3470596 0.2936658 0.3203627 mixture:2 0.4538480 0.3470603 0.2936664 0.3203634 Real Parameter c Group:sexF 2 3 4 mixture:1 0.3470596 0.2936658 0.3203627 mixture:2 0.3470603 0.2936664 0.3203634 Group:sexM 2 3 4 mixture:1 0.3470596 0.2936658 0.3203627 mixture:2 0.3470603 0.2936664 0.3203634 Real Parameter f0 Group:sexF 1 9.186259 Group:sexM 1 9.186259 Output summary for FullHet model Name : pi(~1)p(~mixture + time)c(~1)f0(~1) Npar : 8 (unadjusted=4) -2lnL: -71.99011 AICc : -55.59774 (unadjusted=-63.882296) Beta estimate se lcl ucl pi:(Intercept) 1.1550120 3.408509e+02 -6.669128e+02 669.2227900 p:(Intercept) 0.1342399 1.850896e+01 -3.614333e+01 36.4118070 p:mixture2 0.1527493 0.000000e+00 1.527493e-01 0.1527493 p:time2 -0.2150763 0.000000e+00 -2.150763e-01 -0.2150763 p:time3 0.8151702 0.000000e+00 8.151702e-01 0.8151702 p:time4 20.5093940 1.612136e+04 -3.157735e+04 31618.3680000 c:(Intercept) -0.8092195 1.491105e-01 -1.101476e+00 -0.5169630 f0:(Intercept) -19.8673230 3.132458e+03 -6.159485e+03 6119.7505000 Real Parameter pi Group:sexF mixture:1 0.7604252 Group:sexM mixture:1 0.7604252 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.5335097 0.4798019 0.7209965 1 mixture:2 0.5712589 0.5179705 0.7506645 1 Group:sexM 1 2 3 4 mixture:1 0.5335097 0.4798019 0.7209965 1 mixture:2 0.5712589 0.5179705 0.7506645 1 Real Parameter c Group:sexF 2 3 4 mixture:1 0.3080568 0.3080568 0.3080568 mixture:2 0.3080568 0.3080568 0.3080568 Group:sexM 2 3 4 mixture:1 0.3080568 0.3080568 0.3080568 mixture:2 0.3080568 0.3080568 0.3080568 Real Parameter f0 Group:sexF 1 2.353593e-09 Group:sexM 1 2.353593e-09 Output summary for FullHet model Name : pi(~1)p(~sex)c(~sex)f0(~1) Npar : 6 (unadjusted=5) -2lnL: -73.42423 AICc : -61.19659 (unadjusted=-63.262067) Beta estimate se lcl ucl pi:(Intercept) 8.545351e-06 2048.0178000 -4014.1150000 4014.1150000 p:(Intercept) 3.072807e-01 0.3077730 -0.2959545 0.9105158 p:sexM -3.723105e-01 0.3142225 -0.9881865 0.2435655 c:(Intercept) -4.367175e-01 0.1979751 -0.8247487 -0.0486863 c:sexM -8.222374e-01 0.3082500 -1.4264074 -0.2180673 f0:(Intercept) 6.097045e-01 0.8913361 -1.1373142 2.3567232 Real Parameter pi Group:sexF mixture:1 0.5000021 Group:sexM mixture:1 0.5000021 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.5762214 0.5762214 0.5762214 0.5762214 mixture:2 0.5762214 0.5762214 0.5762214 0.5762214 Group:sexM 1 2 3 4 mixture:1 0.4837483 0.4837483 0.4837483 0.4837483 mixture:2 0.4837483 0.4837483 0.4837483 0.4837483 Real Parameter c Group:sexF 2 3 4 mixture:1 0.3925234 0.3925234 0.3925234 mixture:2 0.3925234 0.3925234 0.3925234 Group:sexM 2 3 4 mixture:1 0.2211538 0.2211538 0.2211538 mixture:2 0.2211538 0.2211538 0.2211538 Real Parameter f0 Group:sexF 1 1.839888 Group:sexM 1 1.839888 Output summary for FullHet model Name : pi(~1)p(~time)c()f0(~1) Npar : 6 (unadjusted=5) -2lnL: -62.91561 AICc : -50.68797 (unadjusted=-52.753445) Beta estimate se lcl ucl pi:(Intercept) -0.0000120239 450.5053600 -882.9905300 882.9905100 p:(Intercept) -0.1851365000 0.2181576 -0.6127254 0.2424524 p:time2 -0.4468513000 0.2747445 -0.9853506 0.0916479 p:time3 -0.6925082000 0.2817886 -1.2448139 -0.1402025 p:time4 -0.5669675000 0.2778924 -1.1116367 -0.0222983 f0:(Intercept) 2.2177088000 0.3611859 1.5097844 2.9256333 Real Parameter pi Group:sexF mixture:1 0.499997 Group:sexM mixture:1 0.499997 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.4538476 0.34706 0.2936661 0.320363 mixture:2 0.4538476 0.34706 0.2936661 0.320363 Group:sexM 1 2 3 4 mixture:1 0.4538476 0.34706 0.2936661 0.320363 mixture:2 0.4538476 0.34706 0.2936661 0.320363 Real Parameter c Group:sexF 2 3 4 mixture:1 0.34706 0.2936661 0.320363 mixture:2 0.34706 0.2936661 0.320363 Group:sexM 2 3 4 mixture:1 0.34706 0.2936661 0.320363 mixture:2 0.34706 0.2936661 0.320363 Real Parameter f0 Group:sexF 1 9.186259 Group:sexM 1 9.186259 Output summary for FullHet model Name : pi(~1)p(~time)c(~1)f0(~1) Npar : 7 (unadjusted=4) -2lnL: -71.99011 AICc : -57.68577 (unadjusted=-63.882296) Beta estimate se lcl ucl pi:(Intercept) -3.009316e-04 3.548308e+03 -6.954684e+03 6954.6836000 p:(Intercept) 1.706246e-01 2.070354e-01 -2.351648e-01 0.5764140 p:time2 -2.171396e-01 3.686967e-01 -9.397853e-01 0.5055060 p:time3 8.102057e-01 5.215654e-01 -2.120624e-01 1.8324739 p:time4 2.167078e+01 2.451764e+04 -4.803291e+04 48076.2530000 c:(Intercept) -8.092236e-01 1.491106e-01 -1.101480e+00 -0.5169669 f0:(Intercept) -2.511835e+01 1.816978e+04 -3.563788e+04 35587.6450000 Real Parameter pi Group:sexF mixture:1 0.4999248 Group:sexM mixture:1 0.4999248 Real Parameter p Group:sexF 1 2 3 4 mixture:1 0.542553 0.4883733 0.727273 1 mixture:2 0.542553 0.4883733 0.727273 1 Group:sexM 1 2 3 4 mixture:1 0.542553 0.4883733 0.727273 1 mixture:2 0.542553 0.4883733 0.727273 1 Real Parameter c Group:sexF 2 3 4 mixture:1 0.308056 0.308056 0.308056 mixture:2 0.308056 0.308056 0.308056 Group:sexM 2 3 4 mixture:1 0.308056 0.308056 0.308056 mixture:2 0.308056 0.308056 0.308056 Real Parameter f0 Group:sexF 1 1.233788e-11 Group:sexM 1 1.233788e-11 ``` --- # Inspect results ```r mouse.results ``` ``` model npar AICc DeltaAICc weight 7 pi(~1)p(~sex)c(~sex)f0(~1) 6 -61.19659 0.000000 0.726572672 9 pi(~1)p(~time)c(~1)f0(~1) 7 -57.68577 3.510822 0.125578082 2 pi(~1)p(~1)c(~1)f0(~1) 4 -56.52289 4.673699 0.070209862 6 pi(~1)p(~mixture + time)c(~1)f0(~1) 8 -55.59774 5.598844 0.044208468 4 pi(~1)p(~mixture)c(~1)f0(~1) 5 -54.46854 6.728045 0.025136377 8 pi(~1)p(~time)c()f0(~1) 6 -50.68797 10.508621 0.003796304 1 pi(~1)p(~1)c()f0(~1) 3 -49.70241 11.494177 0.002319264 5 pi(~1)p(~time + mixture)c()f0(~1) 7 -48.61126 12.585327 0.001344034 3 pi(~1)p(~mixture)c()f0(~1) 4 -47.65911 13.537477 0.000834936 Deviance 7 45.45361 9 46.88773 2 54.24714 6 46.88773 4 54.24714 8 55.96223 1 63.11091 5 55.96223 3 63.11091 ``` --- # Get estimates ```r mouse.results$p.sex.behav$results$real ``` ``` estimate se lcl ucl fixed note pi gF m1 0.5000021 512.0044500 5.562732e-309 1.0000000 p gF t1 m1 0.5762214 0.0751552 4.265467e-01 0.7131057 p gM t1 m1 0.4837483 0.0579550 3.728809e-01 0.5962388 c gF t2 m1 0.3925234 0.0472069 3.047566e-01 0.4878308 c gM t2 m1 0.2211538 0.0406964 1.516078e-01 0.3109110 f0 gF a0 t1 1.8398876 1.6399582 4.110492e-01 8.2354779 ``` + `pi` not relevant as no heterogeneity was fitted + `p gF` capture probability of females / `p gM` capture probability of males + `c gF` recapture probability of females / `c gM` recapture probability of males + `f0 gF` number of females never caugth / `f0 gM` number of males never caugth --- # Get abundance estimates ```r mouse.results$p.sex.behav$results$derived ``` ``` $`N Population Size` estimate lcl ucl 1 45.83989 44.41105 52.23548 2 51.83989 50.41105 58.23548 ``` + Line 1 - Population size of Females + Line 2 - Population size of Males --- # Clean temporary files created by RMark ```r rm(list = ls(all = TRUE)) cleanup(ask = FALSE) ``` --- class: center, middle background-size: cover # Live demo #2: Distance sampling <style type="text/css">pre {font-size: 18px}</style> --- # Load Distance package ```r library(Distance) ``` --- # Read the data (1) + An example of data provided by Eric Rexstad in Distance package + Line transect data of winter wrens in Scotland ```r data(wren_lt) ``` --- # Explore the data (1) ```r wren_lt ``` ``` Region.Label Area Sample.Label Effort object distance Study.Area 1 Montrave 33.2 1 0.416 5 15 Montrave 4 2 Montrave 33.2 1 0.416 6 80 Montrave 4 3 Montrave 33.2 1 0.416 7 35 Montrave 4 4 Montrave 33.2 1 0.416 8 55 Montrave 4 5 Montrave 33.2 1 0.416 12 12 Montrave 4 6 Montrave 33.2 1 0.416 13 75 Montrave 4 7 Montrave 33.2 1 0.416 14 23 Montrave 4 8 Montrave 33.2 1 0.416 15 70 Montrave 4 9 Montrave 33.2 1 0.416 16 75 Montrave 4 10 Montrave 33.2 1 0.416 17 56 Montrave 4 11 Montrave 33.2 2 0.802 25 20 Montrave 4 12 Montrave 33.2 2 0.802 26 24 Montrave 4 13 Montrave 33.2 2 0.802 27 10 Montrave 4 14 Montrave 33.2 2 0.802 28 60 Montrave 4 15 Montrave 33.2 2 0.802 29 65 Montrave 4 16 Montrave 33.2 2 0.802 30 0 Montrave 4 17 Montrave 33.2 2 0.802 41 40 Montrave 4 18 Montrave 33.2 2 0.802 42 32 Montrave 4 19 Montrave 33.2 2 0.802 43 20 Montrave 4 20 Montrave 33.2 2 0.802 44 60 Montrave 4 21 Montrave 33.2 2 0.802 45 40 Montrave 4 22 Montrave 33.2 3 0.802 56 33 Montrave 4 23 Montrave 33.2 3 0.802 57 5 Montrave 4 24 Montrave 33.2 3 0.802 58 35 Montrave 4 25 Montrave 33.2 3 0.802 59 15 Montrave 4 26 Montrave 33.2 3 0.802 60 70 Montrave 4 27 Montrave 33.2 3 0.802 61 5 Montrave 4 28 Montrave 33.2 3 0.802 73 48 Montrave 4 29 Montrave 33.2 3 0.802 74 32 Montrave 4 30 Montrave 33.2 3 0.802 75 0 Montrave 4 31 Montrave 33.2 3 0.802 76 33 Montrave 4 32 Montrave 33.2 3 0.802 77 20 Montrave 4 33 Montrave 33.2 3 0.802 78 35 Montrave 4 34 Montrave 33.2 4 0.598 87 5 Montrave 4 35 Montrave 33.2 4 0.598 88 80 Montrave 4 36 Montrave 33.2 4 0.598 89 80 Montrave 4 37 Montrave 33.2 4 0.598 90 55 Montrave 4 38 Montrave 33.2 4 0.598 98 70 Montrave 4 39 Montrave 33.2 4 0.598 99 75 Montrave 4 40 Montrave 33.2 4 0.598 100 22 Montrave 4 41 Montrave 33.2 4 0.598 101 85 Montrave 4 42 Montrave 33.2 4 0.598 102 75 Montrave 4 43 Montrave 33.2 4 0.598 103 35 Montrave 4 44 Montrave 33.2 4 0.598 104 40 Montrave 4 45 Montrave 33.2 5 0.700 115 35 Montrave 4 46 Montrave 33.2 5 0.700 116 10 Montrave 4 47 Montrave 33.2 5 0.700 117 70 Montrave 4 48 Montrave 33.2 5 0.700 118 70 Montrave 4 49 Montrave 33.2 5 0.700 119 75 Montrave 4 50 Montrave 33.2 5 0.700 132 60 Montrave 4 51 Montrave 33.2 5 0.700 133 30 Montrave 4 52 Montrave 33.2 5 0.700 134 25 Montrave 4 53 Montrave 33.2 5 0.700 135 80 Montrave 4 54 Montrave 33.2 5 0.700 136 80 Montrave 4 55 Montrave 33.2 5 0.700 137 40 Montrave 4 56 Montrave 33.2 5 0.700 138 85 Montrave 4 57 Montrave 33.2 5 0.700 139 55 Montrave 4 58 Montrave 33.2 5 0.700 140 65 Montrave 4 59 Montrave 33.2 6 0.802 149 20 Montrave 4 60 Montrave 33.2 6 0.802 150 30 Montrave 4 61 Montrave 33.2 6 0.802 151 40 Montrave 4 62 Montrave 33.2 6 0.802 152 25 Montrave 4 63 Montrave 33.2 6 0.802 153 45 Montrave 4 64 Montrave 33.2 6 0.802 154 60 Montrave 4 65 Montrave 33.2 6 0.802 162 25 Montrave 4 66 Montrave 33.2 6 0.802 163 35 Montrave 4 67 Montrave 33.2 6 0.802 164 25 Montrave 4 68 Montrave 33.2 6 0.802 165 10 Montrave 4 69 Montrave 33.2 6 0.802 166 2 Montrave 4 70 Montrave 33.2 6 0.802 167 35 Montrave 4 71 Montrave 33.2 7 0.786 174 27 Montrave 4 72 Montrave 33.2 7 0.786 175 35 Montrave 4 73 Montrave 33.2 7 0.786 176 65 Montrave 4 74 Montrave 33.2 7 0.786 177 45 Montrave 4 75 Montrave 33.2 7 0.786 178 80 Montrave 4 76 Montrave 33.2 7 0.786 179 30 Montrave 4 77 Montrave 33.2 7 0.786 187 25 Montrave 4 78 Montrave 33.2 7 0.786 188 70 Montrave 4 79 Montrave 33.2 7 0.786 189 30 Montrave 4 80 Montrave 33.2 7 0.786 190 55 Montrave 4 81 Montrave 33.2 7 0.786 191 60 Montrave 4 82 Montrave 33.2 7 0.786 192 40 Montrave 4 83 Montrave 33.2 7 0.786 193 75 Montrave 4 84 Montrave 33.2 8 0.810 199 70 Montrave 4 85 Montrave 33.2 8 0.810 200 23 Montrave 4 86 Montrave 33.2 8 0.810 201 55 Montrave 4 87 Montrave 33.2 8 0.810 202 60 Montrave 4 88 Montrave 33.2 8 0.810 203 75 Montrave 4 89 Montrave 33.2 8 0.810 204 30 Montrave 4 90 Montrave 33.2 8 0.810 205 65 Montrave 4 91 Montrave 33.2 8 0.810 211 45 Montrave 4 92 Montrave 33.2 8 0.810 212 20 Montrave 4 93 Montrave 33.2 8 0.810 213 55 Montrave 4 94 Montrave 33.2 8 0.810 214 30 Montrave 4 95 Montrave 33.2 9 0.770 217 45 Montrave 4 96 Montrave 33.2 9 0.770 218 10 Montrave 4 97 Montrave 33.2 9 0.770 219 75 Montrave 4 98 Montrave 33.2 9 0.770 223 30 Montrave 4 99 Montrave 33.2 9 0.770 224 70 Montrave 4 100 Montrave 33.2 9 0.770 225 35 Montrave 4 101 Montrave 33.2 10 0.408 229 55 Montrave 4 102 Montrave 33.2 10 0.408 230 23 Montrave 4 103 Montrave 33.2 10 0.408 233 85 Montrave 4 104 Montrave 33.2 11 0.078 234 25 Montrave 4 105 Montrave 33.2 12 0.094 237 90 Montrave 4 106 Montrave 33.2 12 0.094 238 40 Montrave 4 107 Montrave 33.2 12 0.094 240 50 Montrave 4 108 Montrave 33.2 12 0.094 241 90 Montrave 4 109 Montrave 33.2 13 0.408 246 40 Montrave 4 110 Montrave 33.2 13 0.408 247 0 Montrave 4 111 Montrave 33.2 13 0.408 248 100 Montrave 4 112 Montrave 33.2 13 0.408 249 45 Montrave 4 113 Montrave 33.2 13 0.408 250 20 Montrave 4 114 Montrave 33.2 13 0.408 256 40 Montrave 4 115 Montrave 33.2 13 0.408 257 20 Montrave 4 116 Montrave 33.2 13 0.408 258 45 Montrave 4 117 Montrave 33.2 13 0.408 259 45 Montrave 4 118 Montrave 33.2 14 0.542 265 80 Montrave 4 119 Montrave 33.2 14 0.542 266 5 Montrave 4 120 Montrave 33.2 14 0.542 267 35 Montrave 4 121 Montrave 33.2 14 0.542 268 35 Montrave 4 122 Montrave 33.2 14 0.542 269 10 Montrave 4 123 Montrave 33.2 14 0.542 273 50 Montrave 4 124 Montrave 33.2 14 0.542 274 30 Montrave 4 125 Montrave 33.2 14 0.542 275 20 Montrave 4 126 Montrave 33.2 14 0.542 276 35 Montrave 4 127 Montrave 33.2 14 0.542 277 60 Montrave 4 128 Montrave 33.2 14 0.542 278 15 Montrave 4 129 Montrave 33.2 14 0.542 279 10 Montrave 4 130 Montrave 33.2 15 0.472 286 22 Montrave 4 131 Montrave 33.2 15 0.472 287 15 Montrave 4 132 Montrave 33.2 15 0.472 288 30 Montrave 4 133 Montrave 33.2 15 0.472 289 40 Montrave 4 134 Montrave 33.2 15 0.472 294 15 Montrave 4 135 Montrave 33.2 15 0.472 295 35 Montrave 4 136 Montrave 33.2 15 0.472 296 16 Montrave 4 137 Montrave 33.2 15 0.472 297 15 Montrave 4 138 Montrave 33.2 15 0.472 298 25 Montrave 4 139 Montrave 33.2 16 0.378 301 60 Montrave 4 140 Montrave 33.2 16 0.378 303 35 Montrave 4 141 Montrave 33.2 17 0.354 312 10 Montrave 4 142 Montrave 33.2 17 0.354 313 55 Montrave 4 143 Montrave 33.2 17 0.354 314 12 Montrave 4 144 Montrave 33.2 17 0.354 318 35 Montrave 4 145 Montrave 33.2 17 0.354 319 30 Montrave 4 146 Montrave 33.2 17 0.354 320 70 Montrave 4 147 Montrave 33.2 17 0.354 321 25 Montrave 4 148 Montrave 33.2 18 0.400 326 10 Montrave 4 149 Montrave 33.2 18 0.400 327 60 Montrave 4 150 Montrave 33.2 18 0.400 335 10 Montrave 4 151 Montrave 33.2 18 0.400 336 60 Montrave 4 152 Montrave 33.2 18 0.400 337 40 Montrave 4 153 Montrave 33.2 18 0.400 338 50 Montrave 4 154 Montrave 33.2 19 0.040 340 40 Montrave 4 155 Montrave 33.2 19 0.040 341 80 Montrave 4 156 Montrave 33.2 19 0.040 343 25 Montrave 4 ``` --- # Explore the data (2) ```r dim(wren_lt) ``` ``` [1] 156 7 ``` ```r sum(is.na(wren_lt$distance)) # transect with no observations at all ``` ``` [1] 0 ``` --- # Explore the data (3) ```r hist(wren_lt$distance, xlab = "Distance (m)", main = "Winter wren line transects") ``` --- # Explore the data (3) <img src="abundance-tutorial_files/figure-html/unnamed-chunk-26-1.svg" style="display: block; margin: auto;" /> --- # Specify units Choose the units in which density is to be reported - Distance_units - units of measure for perpendicular/radial distances - Effort_units - units of measure for effort (NULL for point transects) - Area_units -units of measure for the study area. ```r conversion.factor <- convert_units("meter", "kilometer","hectare") ``` --- # Fit a simple detection function model ```r wren.hn <- ds(data = wren_lt, key = "hn", adjustment = NULL, convert_units = conversion.factor) ``` --- # Fit a simple detection function model ```r summary(wren.hn) ``` ``` Summary for distance analysis Number of observations : 156 Distance range : 0 - 100 Model : Half-normal key function AIC : 1418.188 Detection function parameters Scale coefficient(s): estimate se (Intercept) 4.105816 0.1327744 Estimate SE CV Average p 0.685037 0.05678866 0.08289868 N in covered region 227.724931 21.47288594 0.09429308 Summary statistics: Region Area CoveredArea Effort n k ER se.ER cv.ER 1 Montrave 33.2 193.2 9.66 156 19 16.14907 1.226096 0.07592366 Abundance: Label Estimate se cv lcl ucl df 1 Total 39.13286 4.399026 0.1124126 31.30227 48.92235 74.24692 Density: Label Estimate se cv lcl ucl df 1 Total 1.1787 0.1325008 0.1124126 0.9428394 1.473565 74.24692 ``` --- # Explore results (1) ```r cutpoints <- c(0,5,10,15,20,30,40,50,65,80,100) plot(wren.hn, breaks = cutpoints, main = "Half normal model, wren line transects") ``` --- # Explore results (1) <img src="abundance-tutorial_files/figure-html/unnamed-chunk-31-1.svg" style="display: block; margin: auto;" /> --- # Explore results (2) ```r str(wren.hn) ``` ``` List of 3 $ ddf :List of 15 ..$ call : language ddf(dsmodel = ~cds(key = "hn", formula = ~1), data = data, method = "ds", meta.data = meta.data, control = control) ..$ data :'data.frame': 156 obs. of 7 variables: .. ..$ Region.Label: Factor w/ 1 level "Montrave": 1 1 1 1 1 1 1 1 1 1 ... .. ..$ Sample.Label: chr [1:156] "1" "1" "1" "1" ... .. ..$ object : int [1:156] 5 6 7 8 12 13 14 15 16 17 ... .. ..$ distance : int [1:156] 15 80 35 55 12 75 23 70 75 56 ... .. ..$ Study.Area : chr [1:156] "Montrave 4" "Montrave 4" "Montrave 4" "Montrave 4" ... .. ..$ observer : num [1:156] 1 1 1 1 1 1 1 1 1 1 ... .. ..$ detected : num [1:156] 1 1 1 1 1 1 1 1 1 1 ... ..$ model : symbol dsmodel ..$ meta.data :List of 7 .. ..$ width : int 100 .. ..$ point : logi FALSE .. ..$ binned : logi FALSE .. ..$ mono : logi TRUE .. ..$ mono.strict: logi TRUE .. ..$ left : num 0 .. ..$ int.range : num [1:2] 0 100 ..$ control :List of 19 .. ..$ optimx.method: chr "nlminb" .. ..$ showit : num 0 .. ..$ estimate : logi TRUE .. ..$ refit : logi TRUE .. ..$ nrefits : num 25 .. ..$ initial : logi NA .. ..$ lowerbounds : logi NA .. ..$ upperbounds : logi NA .. ..$ limit : logi TRUE .. ..$ parscale : logi TRUE .. ..$ maxiter : num 12 .. ..$ standardize : logi TRUE .. ..$ mono.points : num 20 .. ..$ mono.tol : num 1e-08 .. ..$ mono.delta : num 1e-07 .. ..$ debug : logi FALSE .. ..$ nofit : logi FALSE .. ..$ optimx.maxit : num 300 .. ..$ silent : logi FALSE ..$ method : chr "ds" ..$ ds :List of 17 .. ..$ p1 : num 4.11 .. ..$ value : num 708 .. ..$ fevals : num 1 .. ..$ gevals : num 1 .. ..$ niter : num 1 .. ..$ kkt1 : logi TRUE .. ..$ kkt2 : logi TRUE .. ..$ xtime : num 0.001 .. ..$ par : num 4.11 .. ..$ message : chr "" .. ..$ hessian : num [1, 1] 72.3 .. ..$ aux :List of 28 .. .. ..$ maxit : num 300 .. .. ..$ lower : num 1.72 .. .. ..$ upper : num 5.16 .. .. ..$ setlower : logi FALSE .. .. ..$ setupper : logi FALSE .. .. ..$ point : logi FALSE .. .. ..$ int.range : num [1:2] 0 100 .. .. ..$ showit : num 0 .. .. ..$ integral.numeric: NULL .. .. ..$ breaks : NULL .. .. ..$ maxiter : num 12 .. .. ..$ refit : logi TRUE .. .. ..$ nrefits : num 25 .. .. ..$ parscale : logi TRUE .. .. ..$ mono : logi FALSE .. .. ..$ mono.strict : logi TRUE .. .. ..$ binned : logi FALSE .. .. ..$ width : int 100 .. .. ..$ standardize : logi TRUE .. .. ..$ mono.points : num 20 .. .. ..$ mono.tol : num 1e-08 .. .. ..$ mono.delta : num 1e-07 .. .. ..$ debug : logi FALSE .. .. ..$ nofit : logi FALSE .. .. ..$ left : num 0 .. .. ..$ silent : logi FALSE .. .. ..$ optim.history : logi [1:3] NA NA NA .. .. ..$ ddfobj :List of 4 .. .. .. ..$ type : chr "hn" .. .. .. ..$ scale :List of 3 .. .. .. .. ..$ formula : chr "~1" .. .. .. .. ..$ dm : num [1:156, 1] 1 1 1 1 1 1 1 1 1 1 ... .. .. .. .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. .. .. .. ..$ : NULL .. .. .. .. .. .. ..$ : chr "(Intercept)" .. .. .. .. ..$ parameters: num 4.11 .. .. .. ..$ xmat :'data.frame': 156 obs. of 5 variables: .. .. .. .. ..$ object : int [1:156] 5 6 7 8 12 13 14 15 16 17 ... .. .. .. .. ..$ distance: int [1:156] 15 80 35 55 12 75 23 70 75 56 ... .. .. .. .. ..$ binned : logi [1:156] FALSE FALSE FALSE FALSE FALSE FALSE ... .. .. .. .. ..$ observer: num [1:156] 1 1 1 1 1 1 1 1 1 1 ... .. .. .. .. ..$ detected: num [1:156] 1 1 1 1 1 1 1 1 1 1 ... .. .. .. ..$ intercept.only: logi TRUE .. ..$ model :List of 1 .. .. ..$ scalemodel: NULL .. ..$ converge : num 0 .. ..$ bounded : logi FALSE .. ..$ bounds :List of 4 .. .. ..$ lower : num 1.72 .. .. ..$ upper : num 5.16 .. .. ..$ setlower: logi FALSE .. .. ..$ setupper: logi FALSE .. ..$ optim.history: num [1:3] 0 -708.09 4.11 .. ..- attr(*, "details")=List of 1 .. .. ..$ nhatend: num [1, 1] 72.3 .. ..- attr(*, "maximize")= logi FALSE .. ..- attr(*, "npar")= int 1 .. ..- attr(*, "class")= chr [1:2] "mcds" "ds" ..$ par : num 4.11 ..$ lnl : num -708 ..$ hessian : num [1, 1] 56.7 ..$ dsmodel : chr [1:2] "~" "cds(key = \"hn\", formula = ~1)" ..$ criterion : num 1418 ..$ fitted : Named num [1:156] 0.685 0.685 0.685 0.685 0.685 ... .. ..- attr(*, "names")= chr [1:156] "5" "6" "7" "8" ... ..$ Nhat : num 228 ..$ name.message: chr "half-normal key function" ..- attr(*, "class")= chr [1:2] "ds" "ddf" $ dht :List of 1 ..$ individuals:List of 8 .. ..$ bysample :'data.frame': 19 obs. of 9 variables: .. .. ..$ Region : Factor w/ 1 level "Montrave": 1 1 1 1 1 1 1 1 1 1 ... .. .. ..$ Sample : Factor w/ 19 levels "1","10","11",..: 1 2 3 4 5 6 7 8 9 10 ... .. .. ..$ Effort : num [1:19] 0.416 0.408 0.078 0.094 0.408 0.542 0.472 0.378 0.354 0.4 ... .. .. ..$ Sample.Area: num [1:19] 8.32 8.16 1.56 1.88 8.16 ... .. .. ..$ Area : num [1:19] 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 ... .. .. ..$ n : num [1:19] 10 3 1 4 9 12 9 2 7 6 ... .. .. ..$ Nhat : num [1:19] 2.509 0.753 0.251 1.003 2.258 ... .. .. ..$ Nchat : num [1:19] 14.6 4.38 1.46 5.84 13.14 ... .. .. ..$ Dhat : num [1:19] 1.755 0.537 0.936 3.106 1.61 ... .. ..$ summary :'data.frame': 1 obs. of 9 variables: .. .. ..$ Region : Factor w/ 1 level "Montrave": 1 .. .. ..$ Area : num 33.2 .. .. ..$ CoveredArea: num 193 .. .. ..$ Effort : num 9.66 .. .. ..$ n : num 156 .. .. ..$ k : num [1(1d)] 19 .. .. ..$ ER : num 16.1 .. .. ..$ se.ER : num 1.23 .. .. ..$ cv.ER : num 0.0759 .. ..$ N :'data.frame': 1 obs. of 7 variables: .. .. ..$ Label : chr "Total" .. .. ..$ Estimate: num 39.1 .. .. ..$ se : num 4.4 .. .. ..$ cv : num 0.112 .. .. ..$ lcl : num 31.3 .. .. ..$ ucl : num 48.9 .. .. ..$ df : num 74.2 .. ..$ D :'data.frame': 1 obs. of 7 variables: .. .. ..$ Label : chr "Total" .. .. ..$ Estimate: num 1.18 .. .. ..$ se : num 0.133 .. .. ..$ cv : num 0.112 .. .. ..$ lcl : num 0.943 .. .. ..$ ucl : num 1.47 .. .. ..$ df : num 74.2 .. ..$ average.p : num 0.685 .. ..$ cormat : num [1, 1] 1 .. ..$ vc :List of 3 .. .. ..$ total : num [1, 1] 19.4 .. .. ..$ detection:List of 2 .. .. .. ..$ variance: num [1, 1] 10.5 .. .. .. ..$ partial : num [1, 1] -24.4 .. .. ..$ er : num [1, 1] 8.83 .. ..$ Nhat.by.sample:'data.frame': 19 obs. of 9 variables: .. .. ..$ Region.Label: Factor w/ 1 level "Montrave": 1 1 1 1 1 1 1 1 1 1 ... .. .. ..$ Label : chr [1:19] "Montrave1" "Montrave10" "Montrave11" "Montrave12" ... .. .. ..$ Area : num [1:19] 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 ... .. .. ..$ Sample.Label: Factor w/ 19 levels "1","10","11",..: 1 2 3 4 5 6 7 8 9 10 ... .. .. ..$ Effort.x : num [1:19] 0.416 0.408 0.078 0.094 0.408 0.542 0.472 0.378 0.354 0.4 ... .. .. ..$ Nhat : num [1:19] 2.509 0.753 0.251 1.003 2.258 ... .. .. ..$ n : num [1:19] 10 3 1 4 9 12 9 2 7 6 ... .. .. ..$ CoveredArea : num [1:19] 193 193 193 193 193 ... .. .. ..$ Effort.y : num [1:19] 9.66 9.66 9.66 9.66 9.66 9.66 9.66 9.66 9.66 9.66 ... ..- attr(*, "ER_var")= chr [1:3] "R2" "TRUE" "FALSE" ..- attr(*, "class")= chr "dht" $ call: language ds(data = wren_lt, key = "hn", adjustment = NULL, convert.units = conversion.factor) - attr(*, "class")= chr "dsmodel" ``` --- # Get estimates - detection function ```r wren.hn$ddf ``` ``` Distance sampling analysis object Summary for ds object Number of observations : 156 Distance range : 0 - 100 AIC : 1418.188 Detection function: Half-normal key function Detection function parameters Scale coefficient(s): estimate se (Intercept) 4.105816 0.1327744 Estimate SE CV Average p 0.685037 0.05678866 0.08289868 N in covered region 227.724931 21.47288594 0.09429308 ``` --- # Get estimates - density ```r wren.hn$dht ``` ``` Abundance and density estimates from distance sampling Variance : R2, N/L Summary statistics: Region Area CoveredArea Effort n k ER se.ER cv.ER 1 Montrave 33.2 193.2 9.66 156 19 16.14907 1.226096 0.07592366 Abundance: Label Estimate se cv lcl ucl df 1 Total 39.13286 4.399026 0.1124126 31.30227 48.92235 74.24692 Density: Label Estimate se cv lcl ucl df 1 Total 1.1787 0.1325008 0.1124126 0.9428394 1.473565 74.24692 ``` --- # Fit some more models ```r wren.unif.cos <- ds(wren_lt, key = "unif", adjustment = "cos", convert_units = conversion.factor) wren.hr.poly <- ds(wren_lt, key = "hr", adjustment = "poly", convert_units = conversion.factor) ``` --- # Compare these models ```r AIC(wren.hn, wren.hr.poly, wren.unif.cos) ``` ``` df AIC wren.hn 1 1418.188 wren.hr.poly 2 1412.133 wren.unif.cos 3 1416.433 ``` --- # Goodness-Of-Fit tests .pull-left[ ```r z <- gof_ds(wren.hn) ``` <img src="abundance-tutorial_files/figure-html/unnamed-chunk-38-1.svg" style="display: block; margin: auto;" /> ] .pull-right[ ``` Goodness of fit results for ddf object Distance sampling Cramer-von Mises test (unweighted) Test statistic = 0.389692 p-value = 0.0769337 ``` ] --- class: center, middle background-size: cover # Live demo #3: N-mixture models <style type="text/css">pre {font-size: 18px}</style> --- # Ocellated lizard (*Timon lepidus*) <img src="img/ocelle.jpg" width="70%" style="display: block; margin: auto;" /> --- # Quadrat sampling in a Oleron's Island <img src="img/Diapositive35.JPG" width="80%" style="display: block; margin: auto;" /> --- # Load unmarked package ```r library(unmarked) library(AICcmodavg) ``` --- # Read the data (1) ```r dat <- read.csv2('dat/ocellatedlezard.csv') ``` --- # Explore the data (1) ```r head(dat) ``` ``` quadrat Y1 Y2 Y3 burrow habitat temp1 temp2 temp3 p1 p2 p3 1 1 0 2 0 0 dune 21.2 24.8 21.6 April May June 2 2 1 2 5 5 dune 21.7 23.8 23.8 April May June 3 3 0 3 1 1 dune 22.7 22.5 24.2 April May June 4 4 0 1 2 2 dune 23.5 22.4 18.3 April May June 5 5 1 2 2 5 dune 24.2 24.8 20.1 April May June 6 6 1 1 3 10 dune 23.9 26.3 23.9 April May June ``` ```r dim(dat)[1] # number of quadrats ``` ``` [1] 70 ``` --- # Explore the data (2) Look at the detection on the quadrats ```r y <- dat[,2:4] head(y) ``` ``` Y1 Y2 Y3 1 0 2 0 2 1 2 5 3 0 3 1 4 0 1 2 5 1 2 2 6 1 1 3 ``` --- # Explore the data (3) Look at the detection on the quadrats ```r colSums(y) ``` ``` Y1 Y2 Y3 21 47 44 ``` ```r rowSums(y) ``` ``` [1] 2 8 4 3 5 5 4 2 0 3 0 3 0 1 1 2 0 0 2 2 4 4 4 3 2 1 0 1 0 0 0 0 0 3 2 0 2 0 [39] 1 1 0 1 3 0 0 2 3 3 1 4 2 3 1 2 1 3 2 0 0 0 1 0 0 0 0 0 2 0 0 3 ``` --- # Prepare the data for unmarked (1) Collect the covariates at the site (that do not change over field sessions) ```r covsite <- dat[,5:6] # collect site covariates head(covsite) # look at covariates ``` ``` burrow habitat 1 0 dune 2 5 dune 3 1 dune 4 2 dune 5 5 dune 6 10 dune ``` --- # Prepare the data for unmarked (2) Collect the covariates of the field sessions ```r temp <- dat[,7:9] # collect temperature at each field session month <- dat[,10:12] # collect month of each field session obscov <- list(temp = temp, month = month) # make a list with session covariates str(obscov) # look at the covariates ``` ``` List of 2 $ temp :'data.frame': 70 obs. of 3 variables: ..$ temp1: num [1:70] 21.2 21.7 22.7 23.5 24.2 23.9 15.6 23.3 23.2 18.9 ... ..$ temp2: num [1:70] 24.8 23.8 22.5 22.4 24.8 26.3 22.4 24.9 26.3 26 ... ..$ temp3: num [1:70] 21.6 23.8 24.2 18.3 20.1 23.9 23.4 21.6 23.8 27.4 ... $ month:'data.frame': 70 obs. of 3 variables: ..$ p1: chr [1:70] "April" "April" "April" "April" ... ..$ p2: chr [1:70] "May" "May" "May" "May" ... ..$ p3: chr [1:70] "June" "June" "June" "June" ... ``` --- # Prepare the data for unmarked (3) Build the object for unmarked function ```r count <- unmarkedFramePCount(y, siteCovs = covsite, obsCovs = obscov) ``` --- # Fit models Fit a constant model just to check ```r res0 <- pcount(~1 ~1, count, K = 100, mixture = c("P")) ``` --- # Explore the returned object ```r str(res0) ``` ``` Formal class 'unmarkedFitPCount' [package "unmarked"] with 15 slots ..@ K : num 100 ..@ mixture : chr "P" ..@ fitType : chr "pcount" ..@ call : language pcount(formula = ~1 ~ 1, data = count, K = 100, mixture = c("P")) ..@ formula :Class 'formula' language ~1 ~ 1 .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> ..@ data :Formal class 'unmarkedFramePCount' [package "unmarked"] with 5 slots .. .. ..@ y : int [1:70, 1:3] 0 1 0 0 1 1 1 0 0 1 ... .. .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. .. ..$ : NULL .. .. .. .. ..$ : chr [1:3] "Y1" "Y2" "Y3" .. .. ..@ obsCovs :'data.frame': 210 obs. of 2 variables: .. .. .. ..$ temp : num [1:210] 21.2 24.8 21.6 21.7 23.8 23.8 22.7 22.5 24.2 23.5 ... .. .. .. ..$ month: Factor w/ 3 levels "April","June",..: 1 3 2 1 3 2 1 3 2 1 ... .. .. ..@ siteCovs:'data.frame': 70 obs. of 2 variables: .. .. .. ..$ burrow : int [1:70] 0 5 1 2 5 10 0 0 0 5 ... .. .. .. ..$ habitat: Factor w/ 2 levels "clairiere","dune": 2 2 2 2 2 2 2 2 2 2 ... .. .. ..@ mapInfo : NULL .. .. ..@ obsToY : num [1:3, 1:3] 1 0 0 0 1 0 0 0 1 ..@ sitesRemoved : int(0) ..@ estimates :Formal class 'unmarkedEstimateList' [package "unmarked"] with 1 slot .. .. ..@ estimates:List of 2 .. .. .. ..$ state:Formal class 'unmarkedEstimate' [package "unmarked"] with 9 slots .. .. .. .. .. ..@ name : chr "Abundance" .. .. .. .. .. ..@ short.name : chr "lam" .. .. .. .. .. ..@ estimates : Named num 0.392 .. .. .. .. .. .. ..- attr(*, "names")= chr "(Intercept)" .. .. .. .. .. ..@ covMat : num [1, 1] 0.0429 .. .. .. .. .. ..@ fixed : int 1 .. .. .. .. .. ..@ covMatBS : NULL .. .. .. .. .. ..@ invlink : chr "exp" .. .. .. .. .. ..@ invlinkGrad : chr "exp" .. .. .. .. .. ..@ randomVarInfo: list() .. .. .. ..$ det :Formal class 'unmarkedEstimate' [package "unmarked"] with 9 slots .. .. .. .. .. ..@ name : chr "Detection" .. .. .. .. .. ..@ short.name : chr "p" .. .. .. .. .. ..@ estimates : Named num -0.573 .. .. .. .. .. .. ..- attr(*, "names")= chr "(Intercept)" .. .. .. .. .. ..@ covMat : num [1, 1] 0.0952 .. .. .. .. .. ..@ fixed : int 1 .. .. .. .. .. ..@ covMatBS : NULL .. .. .. .. .. ..@ invlink : chr "logistic" .. .. .. .. .. ..@ invlinkGrad : chr "logistic.grad" .. .. .. .. .. ..@ randomVarInfo: list() ..@ AIC : num 388 ..@ opt :List of 6 .. ..$ par : num [1:2] 0.392 -0.573 .. ..$ value : num 192 .. ..$ counts : Named int [1:2] 14 6 .. .. ..- attr(*, "names")= chr [1:2] "function" "gradient" .. ..$ convergence: int 0 .. ..$ message : NULL .. ..$ hessian : num [1:2, 1:2] 68.8 37.6 37.6 31 ..@ negLogLike : num 192 ..@ nllFun :function (parms) ..@ bootstrapSamples: NULL ..@ covMatBS : NULL ..@ TMB : NULL ``` --- # Examine model results ```r summary(res0) ``` ``` Call: pcount(formula = ~1 ~ 1, data = count, K = 100, mixture = c("P")) Abundance (log-scale): Estimate SE z P(>|z|) 0.392 0.207 1.89 0.0585 Detection (logit-scale): Estimate SE z P(>|z|) -0.573 0.308 -1.86 0.063 AIC: 388.2797 Number of sites: 70 optim convergence code: 0 optim iterations: 14 Bootstrap iterations: 0 ``` --- # Get estimates on "natural" scale ```r lambda <- backTransform(res0, type = "state") lambda ``` ``` Backtransformed linear combination(s) of Abundance estimate(s) Estimate SE LinComb (Intercept) 1.48 0.306 0.392 1 Transformation: exp ``` ```r confint(lambda) ``` ``` 0.025 0.975 0.99 2.2 ``` --- # Get estimates on "natural" scale ```r detection <- backTransform(res0, type = "det") detection ``` ``` Backtransformed linear combination(s) of Detection estimate(s) Estimate SE LinComb (Intercept) 0.36 0.0711 -0.573 1 Transformation: logistic ``` ```r confint(detection) ``` ``` 0.025 0.975 0.24 0.51 ``` --- # Check identifiability (1) ```r res0a <- pcount(~1 ~1, count, K = 100, mixture = c("P")) res0b <- pcount(~1 ~1, count, K = 150, mixture = c("P")) res0c <- pcount(~1 ~1, count, K = 200, mixture = c("P")) ``` --- # Check identifiability (2) ```r print(res0a@AIC, digits = 5) ``` ``` [1] 388.28 ``` ```r print(res0b@AIC, digits = 5) ``` ``` [1] 388.28 ``` ```r print(res0c@AIC, digits = 5) ``` ``` [1] 388.28 ``` + If AIC is stable, no problem + If AIC decrease then identifiability issues, try: - Other distributions (but may pose the same problem) - Add some covariates and then check on the best model --- # Fit some models with covariates Fit a model with month as covariate on detection ```r res_month <- pcount(~month ~1, count, K = 100, mixture = c("P")) ``` --- # Examine model results ```r summary(res_month) ``` ``` Call: pcount(formula = ~month ~ 1, data = count, K = 100, mixture = c("P")) Abundance (log-scale): Estimate SE z P(>|z|) 0.256 0.177 1.45 0.147 Detection (logit-scale): Estimate SE z P(>|z|) (Intercept) -1.20 0.310 -3.85 0.000116 monthJune 1.14 0.339 3.37 0.000745 monthMay 1.27 0.344 3.70 0.000214 AIC: 373.7345 Number of sites: 70 optim convergence code: 0 optim iterations: 24 Bootstrap iterations: 0 ``` --- # Compare models ```r res0@AIC ``` ``` [1] 388.2797 ``` ```r res_month@AIC ``` ``` [1] 373.7345 ``` --- # Examine detection probability over months ```r newData <- data.frame(month = c('April','May','June')) pred <- predict(res_month, type = "det", newdata = newData, appendData = TRUE) pred ``` ``` Predicted SE lower upper month 1 0.23 0.055 0.14 0.36 April 2 0.52 0.091 0.35 0.69 May 3 0.49 0.087 0.32 0.65 June ``` --- ## A model with rabbit burrows effect on abundance ```r res_burrow <- pcount(~month ~burrow, count, K = 100, mixture = c("P")) ``` --- # Compare models ```r res_month@AIC ``` ``` [1] 373.7345 ``` ```r res_burrow@AIC ``` ``` [1] 370.7005 ``` --- # Plot burrow effect on abundance ```r newData2 <- data.frame(burrow = seq(min(dat$burrow), max(dat$burrow), by = 1)) E.p <- predict(res_burrow, type = "state", newdata = newData2, appendData = TRUE) plot(Predicted ~ burrow, E.p, type = "l", ylim = c(0,10), xlab = "Number of burrow", ylab = "Expected abundance") lines(lower ~ burrow, E.p, type = "l", col = gray(0.5)) lines(upper ~ burrow, E.p, type = "l", col = gray(0.5)) ``` --- # Plot burrow effect on abundance <img src="abundance-tutorial_files/figure-html/unnamed-chunk-64-1.svg" style="display: block; margin: auto;" /> --- # Goodness-of-fit tests <img src="abundance-tutorial_files/figure-html/unnamed-chunk-65-1.svg" style="display: block; margin: auto;" /> --- # Goodness-of-fit tests ```r print(obs.boot, digits.vals = 4, digits.chisq = 4) ``` ``` Chi-square goodness-of-fit for N-mixture model of 'unmarkedFitPCount' class Observed chi-square statistic = 194.9798 Number of bootstrap samples = 100 P-value = 0.78 Quantiles of bootstrapped statistics: 0% 25% 50% 75% 100% 153.2 195.6 207.0 226.6 257.1 Estimate of c-hat = 0.9302 ``` --- # Model selection ```r res1 <- pcount(~month ~burrow, count, K = 100, mixture = c("P")) res2 <- pcount(~month + burrow ~1, count, K = 100, mixture = c("P")) res3 <- pcount(~month + burrow ~burrow, count, K = 100, mixture = c("P")) res4 <- pcount(~month ~habitat, count, K = 100, mixture = c("P")) res5 <- pcount(~month + habitat ~1, count, K = 100, mixture = c("P")) res6 <- pcount(~month + habitat ~habitat, count, K = 100, mixture = c("P")) ``` --- # Model selection ```r fms <- fitList( "lambda(.) p(.)" = res0, "lambda(burrow) p(month)" = res1, "lambda(.) p(month+burrow)" = res2, "lambda(burrow) p(month+burrow)" = res3, "lambda(habitat) p(month)" = res4, "lambda(.) p(month+habitat)" = res5, "lambda(habitat) p(month+habitat)" = res6 ) ``` --- # Model selection ```r ms <- modSel(fms) ms ``` ``` nPars AIC delta AICwt cumltvWt lambda(burrow) p(month+burrow) 6 370.32 0.00 4.4e-01 0.44 lambda(burrow) p(month) 5 370.70 0.38 3.6e-01 0.80 lambda(.) p(month+burrow) 5 373.61 3.29 8.4e-02 0.88 lambda(.) p(month+habitat) 5 374.50 4.18 5.4e-02 0.94 lambda(habitat) p(month) 5 374.88 4.57 4.5e-02 0.98 lambda(habitat) p(month+habitat) 6 376.48 6.16 2.0e-02 1.00 lambda(.) p(.) 2 388.28 17.96 5.5e-05 1.00 ```