Day 2 ready

psolymos · Mar 18, 2021 · 5132aff · 5132aff
1 parent e7d9403
commit 5132aff
Show file tree

Hide file tree

Showing 9 changed files with 173 additions and 84 deletions.
diff --git a/README.md b/README.md
@@ -33,8 +33,8 @@ Each day will consist of 3 sessions, roughly one hour each, with short breaks in
 | **Day 2** |  **Behavioral complexities**  |    |
 |  |  1. Statistical assumptions and nuisance variables  |  [Slides](https://peter.solymos.org/qpad-workshop/day2-1-intro.pdf)  |
 |  |  2. Behavioral complexities  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-2-behavior.html)  |
-|  |  3. Removal modeling techniques 1.  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-3-removal.html)  |
-|  |  3. Removal modeling techniques 2.  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-3-mixtures.html)  |
+|  |  3. Removal modeling techniques  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-3-removal.html)  |
+|  |  3. Finite mixture models  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-4-mixtures.html)  |
 | **Day 3** |  **The detection process** |    |
 |  |  1. Distance sampling |    |
 |  |  2. Calibrating population density  |    |

diff --git a/day2-1-intro.Rmd b/day2-1-intro.Rmd
@@ -28,6 +28,7 @@ Day 1
 Day 2
 
 - Behavioral complexities
+- Removal models and assumptions
 
 Day 3
 
@@ -64,6 +65,16 @@ LOCAL copies will not be tracked and overwritten by git. You can copy new files
 
 ***
 
+# Update bSims
+
+Patched the Shiny apps, please update!
+
+```{r eval=FALSE}
+remotes::install_github("psolymos/bSims")
+```
+
+***
+
 # What is detectability?
 
 In the most colloquial terms, $\delta$ is the probability

diff --git a/day2-3-removal.Rmd b/day2-3-removal.Rmd
@@ -91,6 +91,8 @@ We use the `fitdistr` function to fit an exponential distribution to these event
 ```{r}
 phi <- 0.5
 (phi_hat <- fitdistr(v1$t, "exponential")$estimate)
+# which is the same as 1/mean(v1$t)
+1/mean(v1$t)
 ```
 
 The estimate is close to the true value of 0.5. Let's plot this
@@ -122,6 +124,23 @@ legend("bottomright", bty="n", lty=c(1,1,1,1), col=c(1,2,4,3),
   legend=c("Empirical", "Expected", "Estimated", "Binned"))
 ```
 
+Fitting survival model to the timt-to-event data
+
+```{r}
+library(survival)
+t21 <- v1$t
+y01 <- ifelse(is.na(t21), 0, 1)
+## censoring at max when we have nondetections (not here but in general)
+t21[is.na(t21)] <- attr(v1, "tlim")[2]
+#time cannot be 0, so we use a small number instead
+t21[t21 == 0] <- 0.001
+## survival object
+sv <- Surv(t21, y01)
+m <- survreg(sv ~ 1, dist="exponential")
+1/exp(coef(m))
+```
+
+
 ## Removal model
 
 The time-removal model, originally developed for estimating wildlife and fish abundances from mark-recapture studies, was later reformulated for avian surveys with the goal of improving estimates of bird abundance by accounting for the availability bias inherent in point-count data. The removal model applied to point-count surveys estimates the probability that a bird is available for detection as a function of the average number of detectable cues that an individual bird gives per minute (singing rate, $\phi$), and the known count duration ($t$).

diff --git a/day2-4-mixtures.Rmd b/day2-4-mixtures.Rmd
@@ -353,47 +353,77 @@ This result tells us mean abundance after correcting for availability bias, but
 
 We'll address these problems next week. Let's just circle back to the assumptions.
 
+## Exercise 1
 
+What other mechanisms can lead to heterogeneity in behavior?
 
+Use the `run_app("bsimsHER")` Shiny app to explore:
 
+- find "edge cases"
+- copy `bsims_all()` calls from Shiny
 
-## Further issues
 
-Stratify the landscape, habitat related behavior (mixture)
+## Exercise 2
 
-## Time to 1st detection info
+How does over/under counting influence estimated vocalization rates?
 
-ARU, lots of hits, not always clear how many inds in recording
-Use tt1 (use the ABMI SM data set)
+(Hint: use the `perception` argument.)
 
-## Methodological diffs
+```{r eval=FALSE}
+library(bSims)
 
-ARU vs human on p
+phi <- 0.5
+B <- 10
+perc <- seq(0.5, 1.5, 0.1)
+
+l <- expand_list(
+  abund_fun = list(identity),
+  duration = 10,
+  vocal_rate = phi,
+  tau = Inf,
+  tint = list(c(3, 5, 10)),
+  perception = perc)
+str(l[1:2])
+
+## a list of bsims_all objects
+## $settings() $new(), $replicate(B, cl)
+b <- lapply(l, bsims_all)
+
+## repeat the runs B times for each setting
+s <- lapply(b, function(z) {
+  z$replicate(B, cl=4)
+})
+
+## removal model
+phi_hat <- t(sapply(s, function(r) sapply(r, estimate)["phi",]))
+
+matplot(perc, phi_hat, lty=1, type="l", col="grey", ylim=c(0, max(phi_hat)))
+lines(perc, apply(phi_hat, 1, median), lwd=2)
+abline(h=phi)
+
+matplot(perc, 1-exp(-1*phi_hat), lty=1, type="l", col="grey", ylim=c(0,1))
+lines(perc, 1-exp(-1*apply(phi_hat, 1, median)), lwd=2)
+abline(h=1-exp(-1*phi), lty=2)
+```
 
-## Std data ove time
+This is how perceived individual ID is deduced using locations:
 
-TSSR/ JDAY effects, but migration?
+```{r eval=FALSE}
+set.seed(1)
+x <- bsims_all(density=0.1)$new()
+perception <- 0.75
 
-## bSims/Shiny
+z <- get_events(x)
+z <- z[!duplicated(z$i),]
+dim(z)
+hc <- hclust(dist(cbind(z$x, z$y)), method="ward.D2")
 
-We can collect all our settings into a `bsims_all` call
+h <- length(unique(z$i)) * perception
+z$j <- cutree(hc, k=min(nrow(z), max(1, round(h))))
 
-```{r}
-xall <- bsims_all(
-  extent=10,
-  road=0.25, edge=0.5,
-  density=c(1, 1, 0),
-  vocal_rate=phi,
-  move_rate=1, movement=0.2,
-  tau=0.8,
-  tint=c(3,5,10), 
-  rint=c(0.5, 1, 1.5))
-xall
+plot(hc)
+table(true=z$i, perceived=z$j)
+plot(z$x, z$y, pch=z$j, col=z$j)
 ```
 
-This call does not evaluate the expression, but it creates a 'closure' with all the info inside to create independent realizations (i.e. none of the layers will match across the runs)
-
-```{r}
-xall$new()
-```
 
diff --git a/docs/README.md b/docs/README.md
@@ -33,8 +33,8 @@ Each day will consist of 3 sessions, roughly one hour each, with short breaks in
 | **Day 2** |  **Behavioral complexities**  |    |
 |  |  1. Statistical assumptions and nuisance variables  |  [Slides](https://peter.solymos.org/qpad-workshop/day2-1-intro.pdf)  |
 |  |  2. Behavioral complexities  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-2-behavior.html)  |
-|  |  3. Removal modeling techniques 1.  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-3-removal.html)  |
-|  |  3. Removal modeling techniques 2.  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-3-mixtures.html)  |
+|  |  3. Removal modeling techniques  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-3-removal.html)  |
+|  |  3. Finite mixture models  |  [Notes](https://peter.solymos.org/qpad-workshop/day2-4-mixtures.html)  |
 | **Day 3** |  **The detection process** |    |
 |  |  1. Distance sampling |    |
 |  |  2. Calibrating population density  |    |

diff --git a/docs/day2-1-intro.pdf b/docs/day2-1-intro.pdf
diff --git a/docs/day2-2-behavior.html b/docs/day2-2-behavior.html
@@ -405,7 +405,7 @@ <h2>Prerequisites</h2>
 ##      of deldir() was changed to (simply) &quot;plot&quot;.
 ##  
 ##      See the help for deldir() and plot.deldir().</code></pre>
-<pre><code>## bSims 0.2-2   2019-12-17      ow-ow-ow-a-la</code></pre>
+<pre><code>## bSims 0.2-3   2021-03-15      kicky-chew</code></pre>
 <pre class="r"><code>library(detect)             # multinomial models</code></pre>
 <pre><code>## Loading required package: Formula</code></pre>
 <pre><code>## Loading required package: stats4</code></pre>

diff --git a/docs/day2-3-removal.html b/docs/day2-3-removal.html
diff --git a/docs/day2-4-mixtures.html b/docs/day2-4-mixtures.html
@@ -384,7 +384,7 @@ <h4 class="author">Peter Solymos <a href="mailto:[email protected]" class="ema
 ##      of deldir() was changed to (simply) &quot;plot&quot;.
 ##  
 ##      See the help for deldir() and plot.deldir().</code></pre>
-<pre><code>## bSims 0.2-2   2019-12-17      ka-ka-ka-kowp-kowp-kowp</code></pre>
+<pre><code>## bSims 0.2-3   2021-03-15      chew-chew-chew</code></pre>
 <pre class="r"><code>library(detect)             # multinomial models</code></pre>
 <pre><code>## Loading required package: Formula</code></pre>
 <pre><code>## Loading required package: stats4</code></pre>
@@ -785,57 +785,69 @@ <h2>Estimating abundance</h2>
 <p>This result tells us mean abundance after correcting for availability bias, but we don’t know what area was effectively sampled, and detection of individuals given availability is probably less than 1 because this happens to be a real data set and it is guaranteed that humans in the forest cannot detect birds that are very far (say &gt; 500 m away).</p>
 <p>We’ll address these problems next week. Let’s just circle back to the assumptions.</p>
 </div>
-<div id="further-issues" class="section level2">
-<h2>Further issues</h2>
-<p>Stratify the landscape, habitat related behavior (mixture)</p>
+<div id="exercise-1" class="section level2">
+<h2>Exercise 1</h2>
+<p>What other mechanisms can lead to heterogeneity in behavior?</p>
+<p>Use the <code>run_app(&quot;bsimsHER&quot;)</code> Shiny app to explore:</p>
+<ul>
+<li>find “edge cases”</li>
+<li>copy <code>bsims_all()</code> calls from Shiny</li>
+</ul>
 </div>
-<div id="time-to-1st-detection-info" class="section level2">
-<h2>Time to 1st detection info</h2>
-<p>ARU, lots of hits, not always clear how many inds in recording Use tt1 (use the ABMI SM data set)</p>
-</div>
-<div id="methodological-diffs" class="section level2">
-<h2>Methodological diffs</h2>
-<p>ARU vs human on p</p>
-</div>
-<div id="std-data-ove-time" class="section level2">
-<h2>Std data ove time</h2>
-<p>TSSR/ JDAY effects, but migration?</p>
-</div>
-<div id="bsimsshiny" class="section level2">
-<h2>bSims/Shiny</h2>
-<p>We can collect all our settings into a <code>bsims_all</code> call</p>
-<pre class="r"><code>xall &lt;- bsims_all(
-  extent=10,
-  road=0.25, edge=0.5,
-  density=c(1, 1, 0),
-  vocal_rate=phi,
-  move_rate=1, movement=0.2,
-  tau=0.8,
-  tint=c(3,5,10), 
-  rint=c(0.5, 1, 1.5))
-xall</code></pre>
-<pre><code>## bSims wrapper object with settings:
-##   extent    : 10
-##   road      : 0.25
-##   edge      : 0.5
-##   density   : 1, 1, 0
-##   vocal_rate: 0.5
-##   move_rate : 1
-##   movement  : 0.2
-##   tau       : 0.8
-##   tint      : 3, 5, 10
-##   rint      : 0.5, 1, 1.5</code></pre>
-<p>This call does not evaluate the expression, but it creates a ‘closure’ with all the info inside to create independent realizations (i.e. none of the layers will match across the runs)</p>
-<pre class="r"><code>xall$new()</code></pre>
-<pre><code>## bSims transcript
-##   1 km x 1 km
-##   stratification: HER
-##   total abundance: 106
-##   duration: 10 min
-##   detected: 4 heard
-##   1st event detected by breaks:
-##     [0, 3, 5, 10 min]
-##     [0, 50, 100, 150 m]</code></pre>
+<div id="exercise-2" class="section level2">
+<h2>Exercise 2</h2>
+<p>How does over/under counting influence estimated vocalization rates?</p>
+<p>(Hint: use the <code>perception</code> argument.)</p>
+<pre class="r"><code>library(bSims)
+
+phi &lt;- 0.5
+B &lt;- 10
+perc &lt;- seq(0.5, 1.5, 0.1)
+
+l &lt;- expand_list(
+  abund_fun = list(identity),
+  duration = 10,
+  vocal_rate = phi,
+  tau = Inf,
+  tint = list(c(3, 5, 10)),
+  perception = perc)
+str(l[1:2])
+
+## a list of bsims_all objects
+## $settings() $new(), $replicate(B, cl)
+b &lt;- lapply(l, bsims_all)
+
+## repeat the runs B times for each setting
+s &lt;- lapply(b, function(z) {
+  z$replicate(B, cl=4)
+})
+
+## removal model
+phi_hat &lt;- t(sapply(s, function(r) sapply(r, estimate)[&quot;phi&quot;,]))
+
+matplot(perc, phi_hat, lty=1, type=&quot;l&quot;, col=&quot;grey&quot;, ylim=c(0, max(phi_hat)))
+lines(perc, apply(phi_hat, 1, median), lwd=2)
+abline(h=phi)
+
+matplot(perc, 1-exp(-1*phi_hat), lty=1, type=&quot;l&quot;, col=&quot;grey&quot;, ylim=c(0,1))
+lines(perc, 1-exp(-1*apply(phi_hat, 1, median)), lwd=2)
+abline(h=1-exp(-1*phi), lty=2)</code></pre>
+<p>This is how perceived individual ID is deduced using locations:</p>
+<pre class="r"><code>set.seed(1)
+x &lt;- bsims_all(density=0.1)$new()
+perception &lt;- 0.75
+
+z &lt;- get_events(x)
+z &lt;- z[!duplicated(z$i),]
+dim(z)
+hc &lt;- hclust(dist(cbind(z$x, z$y)), method=&quot;ward.D2&quot;)
+
+h &lt;- length(unique(z$i)) * perception
+z$j &lt;- cutree(hc, k=min(nrow(z), max(1, round(h))))
+
+plot(hc)
+table(true=z$i, perceived=z$j)
+plot(z$x, z$y, pch=z$j, col=z$j)</code></pre>
 </div>