-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathWeek-12-lab.Rmd
92 lines (58 loc) · 5.55 KB
/
Week-12-lab.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
Week 12 Lab
========================================================
Before we even get to mark recapture models, we need to introduce a couple of tricks in JAGS that can be used to handle non-standard distributions.
The 'zeros' trick
-------------------
The ‘zeroes trick’ is JAGS goes as follows. Let’s say we want to model
$$
Y \sim G(\theta)
$$
but $G()$ is not a distribution that is built into JAGS. We can achieve the same effect (i.e. have the same likelihood enter the model) if we instead say that we observed a value of $Y=0$ from the following model
$$
Y \sim Pois(\lambda)
$$
where
$$
\lambda = -\mbox{log}(g(y,\theta))
$$
and $g(y,\theta)$ is the likelihood associated with distribution $G$. In other words, the probability of getting a 0 from $\mbox{Pois}(-\mbox{log}(g(y,\theta))$ is just $g(y,\theta)$, which is the likelihood we wanted in the first place.
(Reminder, the pmf of a Poisson is $\frac{\lambda^{y}}{y!}e^{-\lambda}$ so the probability of getting $y=0$ is $e^{-\lambda}$.)
In practice, we add a constant $C$ to get $\lambda = -\mbox{log}(g(y,\theta))+C$ to ensure that $\lambda > 0$.)
Let’s take for example the Cauchy distribution
$$
f(y|\theta) = \frac{1}{\pi}\frac{1}{1+(y-\theta)^{2}}
$$
The Cauchy is an evil distribution, as you all know. Nevertheless, we can use the zeros trick to sample from it as follows:
```{r eval=FALSE}
Zeros[i]<-0
lambda[i]<- -log(1/(1+pow(y[i]-theta,2)))
Zeros[i]~dpois(lambda[i])
```
Notice that I just need the portion of the pdf that involves the data and the parameters because all that matters is that I get something proportional to the correct likelihood, and so I have dropped the constants involved in $f(y|\theta)$.
The 'ones' trick
-------------------
The ‘ones trick’ is similar, but it uses a dummy variable of 1 drawn from a Binomial distribution:
$$
Y \sim \mbox{Bern}(p) \\
p = g(y,\theta)
$$
The code for the Cauchy looks like
```{r eval=FALSE}
Ones[i]<-1
p[i]<- 1/(1+pow(y[i]-theta,2))
Ones[i]~dbern(p[i])
```
##Initial values
Sometimes, the hardest part about getting JAGS to run successfully is supplying it with good initial values. For mark-recapture models I suggest trying NAs up to the first capture, 1 through the period of known survival (from marking to last resighting) and 0 for the period where its fate is unknown (after the last resighting).
##First, a warm up
Let’s look through some [code](https://github.com/hlynch/Bayesian2020/tree/master/_data/Week12LabWarmup.R) for occupancy modeling that will serve as a warm up to the Lab and the Problem Set. This code simulates the detection of an amphibian (usually detected through calling) at 100 sites assuming 30 surveys each breeding season, and detection probabilities that vary as a function of humidity. (This is a bit of a contrived example. We assume 30 replicate surveys at each site, assuming a closed population and assuming that each site has fixed humidity. Nevertheless, these assumptions simplify the data simulation and are easily relaxed.)
(Note that this warm up exercise involves occupancy of an animal species at different sites, whereas the lab itself involves mark-recapture of individual animals over several years. Occupancy and mark-recapture are related problems, but the actual model code will be different.)
##Fitting mark-recapture models
Mark-recapture models are a huge field of modeling, and so we will only just barely scratch the surface on what can be done. Nevertheless, they are a common application of Bayesian models, and they will give us some practice writing Bayesian models in a wider variety of applications.
We are going to start with a fairly straightforward dataset on the survival of the European Dipper, which was described in the paper by McCarthy and Masters (2005) that you read for class this week. I have posted the survival data [here](https://github.com/hlynch/Bayesian2020/tree/master/_data/dipper_data.R). (While the code used by McCarthy and Masters is online as part of the supplement, no peaking! You’ll learn more working through it from the beginning.)
Write a mark-recapture model to estimate survival (constant over time) and resighting rate (also constant over time) for the Dipper. We’ll start with an uninformative prior, and we have time can go back and discuss how McCarthy and Masters generated a more informative prior if we have time.
**Exercise #1**: Use the Method #1 (Brute force) approach described in lecture on Monday. This is the most straightforward and the easiest to code.
**Exercise #2**: Use the Method #2 (Modelling the entire capture history) approach described on Monday. You’ll have to think carefully how to represent the likelihood for the entire recapture history. (Hint: You may need the ‘ones’ or ‘zeroes’ trick.)
##FAQ
**Question**: Why not turn the 0s at the end of each row to NA, since you did not detect anything?
**Answer**: The answer is that until you fit the model, you don't know how informative (or not) a 0 at the end of the time series is. Lets say you fit the model and the model estimates detection probability is 0.999. Then (conditional on the animal being alive), you are almost certain to detect it, which means a non-detection in the last year is almost certainly a death, and thats highly informative for the model. On the other hand, if detection probability is 0.001, then the 0 at the end has very little value to the model, and it could be set to NA with little change to the posteriors. So, _a priori_, the model actually doesn't know how "informative" the non-detections are in the period at the end where the animal could be dead or alive. (edited)