forked from FukamiLab/BIO202
-
Notifications
You must be signed in to change notification settings - Fork 0
/
05-C-Advanced-bayesian-example.html
841 lines (756 loc) · 45.3 KB
/
05-C-Advanced-bayesian-example.html
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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<meta name="author" content="Jes, Lizzie" />
<title>Advanced Bayesian model example</title>
<script src="site_libs/jquery-1.11.3/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/bootstrap.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="site_libs/jqueryui-1.11.4/jquery-ui.min.js"></script>
<link href="site_libs/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="site_libs/tocify-1.9.1/jquery.tocify.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
div.sourceCode { overflow-x: auto; }
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; line-height: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; } /* Keyword */
code > span.dt { color: #902000; } /* DataType */
code > span.dv { color: #40a070; } /* DecVal */
code > span.bn { color: #40a070; } /* BaseN */
code > span.fl { color: #40a070; } /* Float */
code > span.ch { color: #4070a0; } /* Char */
code > span.st { color: #4070a0; } /* String */
code > span.co { color: #60a0b0; font-style: italic; } /* Comment */
code > span.ot { color: #007020; } /* Other */
code > span.al { color: #ff0000; font-weight: bold; } /* Alert */
code > span.fu { color: #06287e; } /* Function */
code > span.er { color: #ff0000; font-weight: bold; } /* Error */
code > span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
code > span.cn { color: #880000; } /* Constant */
code > span.sc { color: #4070a0; } /* SpecialChar */
code > span.vs { color: #4070a0; } /* VerbatimString */
code > span.ss { color: #bb6688; } /* SpecialString */
code > span.im { } /* Import */
code > span.va { color: #19177c; } /* Variable */
code > span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code > span.op { color: #666666; } /* Operator */
code > span.bu { } /* BuiltIn */
code > span.ex { } /* Extension */
code > span.pp { color: #bc7a00; } /* Preprocessor */
code > span.at { color: #7d9029; } /* Attribute */
code > span.do { color: #ba2121; font-style: italic; } /* Documentation */
code > span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code > span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
</style>
<style type="text/css">
pre:not([class]) {
background-color: white;
}
</style>
<style type="text/css">
h1 {
font-size: 34px;
}
h1.title {
font-size: 38px;
}
h2 {
font-size: 30px;
}
h3 {
font-size: 24px;
}
h4 {
font-size: 18px;
}
h5 {
font-size: 16px;
}
h6 {
font-size: 12px;
}
.table th:not([align]) {
text-align: left;
}
</style>
<link rel="stylesheet" href="style.css" type="text/css" />
</head>
<body>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
code {
color: inherit;
background-color: rgba(0, 0, 0, 0.04);
}
img {
max-width:100%;
height: auto;
}
.tabbed-pane {
padding-top: 12px;
}
button.code-folding-btn:focus {
outline: none;
}
</style>
<style type="text/css">
/* padding for bootstrap navbar */
body {
padding-top: 51px;
padding-bottom: 40px;
}
/* offset scroll position for anchor links (for fixed navbar) */
.section h1 {
padding-top: 56px;
margin-top: -56px;
}
.section h2 {
padding-top: 56px;
margin-top: -56px;
}
.section h3 {
padding-top: 56px;
margin-top: -56px;
}
.section h4 {
padding-top: 56px;
margin-top: -56px;
}
.section h5 {
padding-top: 56px;
margin-top: -56px;
}
.section h6 {
padding-top: 56px;
margin-top: -56px;
}
</style>
<script>
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark it active
menuAnchor.parent().addClass('active');
// if it's got a parent navbar menu mark it active as well
menuAnchor.closest('li.dropdown').addClass('active');
});
</script>
<div class="container-fluid main-container">
<!-- tabsets -->
<script>
$(document).ready(function () {
window.buildTabsets("TOC");
});
</script>
<!-- code folding -->
<script>
$(document).ready(function () {
// move toc-ignore selectors from section div to header
$('div.section.toc-ignore')
.removeClass('toc-ignore')
.children('h1,h2,h3,h4,h5').addClass('toc-ignore');
// establish options
var options = {
selectors: "h1,h2,h3",
theme: "bootstrap3",
context: '.toc-content',
hashGenerator: function (text) {
return text.replace(/[.\\/?&!#<>]/g, '').replace(/\s/g, '_').toLowerCase();
},
ignoreSelector: ".toc-ignore",
scrollTo: 0
};
options.showAndHide = false;
options.smoothScroll = true;
// tocify
var toc = $("#TOC").tocify(options).data("toc-tocify");
});
</script>
<style type="text/css">
#TOC {
margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
position: relative;
width: 100%;
}
}
.toc-content {
padding-left: 30px;
padding-right: 40px;
}
div.main-container {
max-width: 1200px;
}
div.tocify {
width: 20%;
max-width: 260px;
max-height: 85%;
}
@media (min-width: 768px) and (max-width: 991px) {
div.tocify {
width: 25%;
}
}
@media (max-width: 767px) {
div.tocify {
width: 100%;
max-width: none;
}
}
.tocify ul, .tocify li {
line-height: 20px;
}
.tocify-subheader .tocify-item {
font-size: 0.90em;
padding-left: 25px;
text-indent: 0;
}
.tocify .list-group-item {
border-radius: 0px;
}
.tocify-subheader {
display: inline;
}
.tocify-subheader .tocify-item {
font-size: 0.95em;
}
</style>
<!-- setup 3col/9col grid for toc_float and main content -->
<div class="row-fluid">
<div class="col-xs-12 col-sm-4 col-md-3">
<div id="TOC" class="tocify">
</div>
</div>
<div class="toc-content col-xs-12 col-sm-8 col-md-9">
<div class="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html"></a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
<li>
<a href="index.html">Home</a>
</li>
<li>
<a href="00-computer-setup.html">Computer Setup</a>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W1
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="01-A-R-intro.html">Intro to R</a>
</li>
<li>
<a href="01-B-Rmarkdown-intro.html">R markdown</a>
</li>
<li>
<a href="01-C-R-workshop.html">R workshop</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W2
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="02-A-tidyr.html">ggplot2 and tidyr</a>
</li>
<li>
<a href="02-B-git-intro.html">Intro to git</a>
</li>
<li>
<a href="02-C-student-projects.html">Project introductions</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W3
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="03-A-data-exploration.html">Data exploration</a>
</li>
<li>
<a href="03-B-linear-models.html">Linear models</a>
</li>
<li>
<a href="03-C-heterogeneity.html">Heterogeneity</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W4
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="04-A-mixed-models.html">Mixed effects models</a>
</li>
<li>
<a href="04-B-binary-data.html">Binary & proportional data</a>
</li>
<li>
<a href="04-C-zero-data.html">Zero-inflated data</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W5
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="05-A-Bayesian-linear-models.html">Bayesian linear models</a>
</li>
<li>
<a href="05-B-Bayesian-priors.html">Bayesian inference with prior information</a>
</li>
<li>
<a href="05-C-Advanced-bayesian-example.html">Advanced Bayesian model example</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W6
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="06-A-unconstrained-ordination.html">Unconstrained ordination</a>
</li>
<li>
<a href="06-B-constrained-ordination.html">Constrained ordination</a>
</li>
<li>
<a href="06-C-matrix-comparison.html">Comparing multivariate data</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W8
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="08-A-mapping.html">Visualizing spatial data</a>
</li>
<li>
<a href="08-B-spatial-regression.html">Spatial regression</a>
</li>
<li>
<a href="08-C-spatial-ordination.html">Ordination approach to spatial analysis</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
W9
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="09-A-time-series.html">Time series</a>
</li>
<li>
<a href="09-B-networks.html">Network analysis</a>
</li>
<li>
<a href="09-C-occupancy-models.html">Occupancy models</a>
</li>
</ul>
</li>
</ul>
<ul class="nav navbar-nav navbar-right">
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<div class="fluid-row" id="header">
<h1 class="title toc-ignore">Advanced Bayesian model example</h1>
<h4 class="author"><em>Jes, Lizzie</em></h4>
</div>
<p><strong>Assigned Reading:</strong></p>
<blockquote>
<p>Section 14.2 in Korner-Nievergelt et al. 2015. <em>Bayesian data analysis in ecology using linear models with R, BUGS, and Stan.</em> Elsevier. <a href="http://www.sciencedirect.com/science/book/9780128013700">link</a></p>
</blockquote>
<p><strong>Optional Reading:</strong></p>
<blockquote>
<p>Chapters 16 and 17 in Korner-Nievergelt et al. 2015. <em>Bayesian data analysis in ecology using linear models with R, BUGS, and Stan.</em> Elsevier. <a href="http://www.sciencedirect.com/science/book/9780128013700">link</a></p>
<p>* These chapters contain a data analysis checklist and information on how to report analyses in scientific papers. Chapter 17 is particularly useful because it gives verbatim text for how you should describe the models that you fit.</p>
</blockquote>
<div id="key-points" class="section level3">
<h3>Key Points</h3>
<p>We want to fit a zero-inflated Poisson mixed-effects model, but the <code>hurdle()</code> and <code>zeroinfl()</code> functions from the <code>pscl</code> package don’t allow random effects. A Bayesian approach is more flexible and will allow us to fit this model.</p>
<p><strong>The Model</strong></p>
<p>Today’s example is a zero-inflated Poisson model with a random intercept in the portion of the model controlling zero-inflation. This example is taken directly from Section 14.2 of Korner-Nievergelt et al. (2015). (See above citation.)</p>
<p>In the following model, the subscript <span class="math inline">\(i\)</span> denotes a specific nestID and the subscript <span class="math inline">\(t\)</span> denotes a specific sampling year. Bold variables are vectors.</p>
<p>Create an unobserved (latent) variable <span class="math inline">\(\mathbf{z}\)</span> that equals 1 with probability <span class="math inline">\(\boldsymbol{\theta}\)</span>.</p>
<p><span class="math display">\[z_{it} \sim Bernoulli(\theta_{it})\]</span></p>
<p>If <span class="math inline">\(\mathbf{z}\)</span> is 1, then the observed variable <span class="math inline">\(\mathbf{y}\)</span> is 0. Otherwise, <span class="math inline">\(\mathbf{y}\)</span> is Poisson distributed with mean <span class="math inline">\(\boldsymbol{\lambda}\)</span> (I.e. This is a mixture model because <span class="math inline">\(\mathbf{y}\)</span> can be 0 for two different reasons.)</p>
<p><span class="math display">\[y_{it} \sim Poisson(\lambda_{it}(1 - z_{it}))\]</span></p>
<p>The log odds that <span class="math inline">\(\mathbf{z}\)</span> is 1 depends linearly on year (<span class="math inline">\(a_2\)</span>) and different nestIDs are allowed to have different effects (<span class="math inline">\(\boldsymbol{\epsilon_{nestID}}\)</span>) on the overall mean (<span class="math inline">\(a_1\)</span>)</p>
<p><span class="math display">\[logit(\theta_{it}) = a_1 + a_2 year_t + \epsilon_{nestID[i]}\]</span></p>
<p>However, the effects of different nestIDs are constrained to come from a Normal distribution with mean <span class="math inline">\(0\)</span> and variance <span class="math inline">\(\sigma_n^2\)</span>. Hence they are a “random” effect.</p>
<p><span class="math display">\[\epsilon_{nestID} \sim Norm(0, \sigma_n)\]</span></p>
<p>If <span class="math inline">\(\mathbf{z}\)</span> is <span class="math inline">\(0\)</span>, then (the log of) the expected value of <span class="math inline">\(\mathbf{y}\)</span> depends linearly on year (<span class="math inline">\(b_2\)</span>).</p>
<p><span class="math display">\[log(\lambda_{it}) = b_1 + b_2 year_t\]</span></p>
</div>
<div id="analysis-example" class="section level3">
<h3>Analysis Example</h3>
<p>In review, the zero-inflated Poisson model is a mixture of binomial and a Poisson distribution. Like all mixed models, it addresses data generated by multiple processes, as seen in life sciences, and the distributions that result from such data. Neglecting these different processes could yield biased inferences. Using multi-level models are especially useful to help reveal the different processes underlying our data.</p>
<p>Count data typically contains a high proportion of zeros and varying data, as in our example of black-stork nestlings. The number of surviving nestlings is often zero due to depradation or natural events. Furthermore, in nests that succeed, numbers of surviving nestlings differ greatly across nests. Hence, these count data (or zero-inflated data) result from two distinct processes, one producing zero counts and one producing non-zero count variance. Recall the bimodal histogram of zero-inflated data, with one peak at zero and another at a greater value. To address zero-laden bimodal distributions, zero-inflated models blend a Bernoulli model (zero vs. nonzero) with a conditional Poisson model–conditional on the Bernoulli process being nonzero. To fit zero-inflated models, we use the function <code>zeroinfl</code> from the package pscl.</p>
<p>The question posed in our example is: Did the number of black-stork nestlings surviving in Latvia decline over time? The authors use a zero-inflated Poisson model to estimate temporal trends in nest survival and, separately, the number of surviving nestlings in extant nests. The nests were repeatedly measured over 17 years.</p>
<p><br></p>
<p><strong>First steps: Read and prepare data for STAN</strong> Example data: Breeding success of Black-stork in Latvia. The data were collected and provided by Maris Stradz. They contain the number of nestlings of more than 300 black-stork nests in different years.</p>
<p><em>Read data</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(blmeco)</code></pre></div>
<pre><code>## Loading required package: MASS</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">data</span>(blackstork)
dat <-<span class="st"> </span>blackstork
<span class="co"># define nest as a factor </span>
dat$nest <-<span class="st"> </span><span class="kw">factor</span>(dat$nest)
<span class="co"># z-transform year</span>
dat$year.z <-<span class="st"> </span><span class="kw">as.numeric</span>(<span class="kw">scale</span>(dat$year))</code></pre></div>
<p><em>Prepare data for Stan</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">y <-<span class="st"> </span>dat$njuvs
N <-<span class="st"> </span><span class="kw">nrow</span>(dat)
nest <-<span class="st"> </span><span class="kw">as.numeric</span>(dat$nest)
Nnests <-<span class="st"> </span><span class="kw">nlevels</span>(dat$nest)
year <-<span class="st"> </span>dat$year.z
datax <-<span class="st"> </span><span class="kw">c</span>(<span class="st">"y"</span>, <span class="st">"N"</span>, <span class="st">"nest"</span>, <span class="st">"Nnests"</span>, <span class="st">"year"</span>)</code></pre></div>
<p><br></p>
<p><strong>Run STAN</strong></p>
<p>Next, we load the Stan code from the text file. In the book, the authors use the function <code>log_sum_exp</code> to calculate the density function of the zero count. The function <code>increment_log_prob</code> defines likelihood in Stan; however, the updated Stan language uses <code>target += u;</code> instead. Functions <code>bernoulli_log</code> and <code>poisson_log</code> define the density functions’ logarithms; they’ve also been updated to <code>bernoulli_lpdf</code> and <code>poisson_lpdf</code>.</p>
<p>There appeared to be an issue with the original code in the book, namely with “sigmanest.” I changed it to “sigma” to solve the issue.</p>
<p>Once you run the stan() function, save the model object to an RData file. Then you can load the fitted model and samples without having to run the stan() function every time you re-open the R script: load(“code/05-C-stanmod.RData”).</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(rstan)</code></pre></div>
<pre><code>## Loading required package: ggplot2</code></pre>
<pre><code>## Loading required package: StanHeaders</code></pre>
<pre><code>## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)</code></pre>
<pre><code>## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Fit the model</span>
<span class="co"># mod <- stan(file = "code/05-C.stan", data=datax, chains=3, iter=1000)</span>
<span class="co"># Save the model output</span>
<span class="co"># save(mod, file = "data/05-C-stanmod.RData")</span>
<span class="co"># Re-load the model output</span>
<span class="kw">load</span>(<span class="st">"data/05-C-stanmod.RData"</span>)
<span class="co"># View the parameter estimates</span>
<span class="kw">print</span>(mod, <span class="dt">pars =</span> <span class="kw">c</span>(<span class="st">"a"</span>, <span class="st">"b"</span>, <span class="st">"sigma"</span>))</code></pre></div>
<pre><code>## Inference for Stan model: zeroinfl.
## 3 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a[1] -1.01 0.00 0.12 -1.26 -1.09 -1.01 -0.94 -0.79 1500 1
## a[2] 0.57 0.00 0.11 0.36 0.50 0.56 0.64 0.79 1500 1
## b[1] 0.89 0.00 0.03 0.84 0.88 0.89 0.91 0.95 1500 1
## b[2] -0.05 0.00 0.02 -0.10 -0.07 -0.05 -0.04 -0.01 1500 1
## sigma 0.95 0.01 0.17 0.63 0.83 0.94 1.05 1.30 558 1
##
## Samples were drawn using NUTS(diag_e) at Thu Oct 26 16:14:33 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).</code></pre>
<p>The <code>mean</code> in the first column is the average of the simulated values for each parameter’s marginal posterior distribution. <code>se_mean</code> refers to the standard error of the simulated values (Monte Carlo uncertainty). <code>sd</code> is the simulations’ sample standard deviation, corresponding to the parameter estimates’ standard error.</p>
<p><em>Look at convergence</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">traceplot</span>(mod, <span class="st">"a"</span>)</code></pre></div>
<p><img src="images/05-C/unnamed-chunk-5-1.png" width="672" /></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">traceplot</span>(mod, <span class="st">"b"</span>)</code></pre></div>
<p><img src="images/05-C/unnamed-chunk-5-2.png" width="672" /></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">traceplot</span>(mod, <span class="st">"sigma"</span>)</code></pre></div>
<p><img src="images/05-C/unnamed-chunk-5-3.png" width="672" /></p>
<p><em>Draw plots of the parameter estimates</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(mod)</code></pre></div>
<pre><code>## 'pars' not specified. Showing first 10 parameters by default.</code></pre>
<pre><code>## ci_level: 0.8 (80% intervals)</code></pre>
<pre><code>## outer_level: 0.95 (95% intervals)</code></pre>
<p><img src="images/05-C/unnamed-chunk-6-1.png" width="672" /></p>
<p><br></p>
<p><strong>Predictive model checking</strong></p>
<p>Before interpreting results, we need to assess the model fit using predictive model checking. Accordingly, we first simulate replicated numbers of nestlings for every year and nest in the data set per the model fit.</p>
<p><em>Extract number of simulations</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">modsims <-<span class="st"> </span><span class="kw">extract</span>(mod)
nsim <-<span class="st"> </span><span class="kw">length</span>(modsims$lp_)
yrep <-<span class="st"> </span><span class="kw">matrix</span>(<span class="dt">nrow=</span><span class="kw">length</span>(y), <span class="dt">ncol=</span>nsim)
Xmat <-<span class="st"> </span><span class="kw">model.matrix</span>(~year)
for(i in <span class="dv">1</span>:nsim){
theta <-<span class="st"> </span><span class="kw">plogis</span>(Xmat%*%modsims$a[i,] +<span class="st"> </span>
<span class="st"> </span>modsims$groupefftheta[i,nest])
z <-<span class="st"> </span><span class="kw">rbinom</span>(<span class="kw">length</span>(y), <span class="dt">prob=</span>theta, <span class="dt">size=</span><span class="dv">1</span>)
lambda <-<span class="st"> </span>(<span class="dv">1</span>-z) *<span class="kw">exp</span>(Xmat%*%modsims$b[i,])
yrep[,i] <-<span class="st"> </span><span class="kw">rpois</span>(<span class="kw">length</span>(y), <span class="dt">lambda=</span>lambda)
}</code></pre></div>
<p><em>Sort the nests according to the first year for the plot</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">first <-<span class="st"> </span><span class="kw">tapply</span>(year, nest, min)
first <-<span class="st"> </span>first[<span class="kw">match</span>(nest, <span class="dv">1</span>:Nnests)]
nestnamenum <-<span class="st"> </span>first*<span class="dv">1000</span> +<span class="st"> </span>nest
nestnamenum <-<span class="st"> </span><span class="kw">as.numeric</span>(<span class="kw">factor</span>(nestnamenum))</code></pre></div>
<p><em>Visualize temporal pattern of zeros and non-zeros (Figure 14.3)</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">par</span>(<span class="dt">mfrow=</span><span class="kw">c</span>(<span class="dv">1</span>,<span class="dv">5</span>), <span class="dt">mar=</span><span class="kw">c</span>(<span class="dv">0</span>,<span class="dv">0</span>,<span class="dv">2</span>,<span class="fl">0.1</span>), <span class="dt">oma=</span><span class="kw">c</span>(<span class="dv">2</span>,<span class="dv">2</span>,<span class="dv">0</span>,<span class="dv">0</span>))
<span class="co"># First plot the real data</span>
<span class="kw">plot</span>(year, nestnamenum, <span class="dt">col=</span><span class="kw">rgb</span>(<span class="dv">1</span>,<span class="dv">0</span>,<span class="dv">0</span>,<span class="fl">0.5</span>), <span class="dt">pch=</span><span class="dv">16</span>, <span class="dt">cex=</span>dat$njuv/<span class="dv">5</span>, <span class="dt">axes=</span><span class="ot">FALSE</span>, <span class="dt">main=</span><span class="st">"observed"</span>)
<span class="kw">points</span>(year[dat$njuv==<span class="dv">0</span>], nestnamenum[dat$njuv==<span class="dv">0</span>], <span class="dt">cex=</span><span class="fl">0.3</span>) <span class="co"># add zeros</span>
<span class="co"># Add data from 4 simulated data-sets </span>
for(r in <span class="dv">1</span>:<span class="dv">4</span>){
<span class="kw">plot</span>(year, nestnamenum, <span class="dt">col=</span><span class="kw">rgb</span>(<span class="dv">1</span>,<span class="dv">0</span>,<span class="dv">0</span>,<span class="fl">0.5</span>), <span class="dt">pch=</span><span class="dv">16</span>, <span class="dt">cex=</span>yrep[,r]/<span class="dv">5</span>, <span class="dt">axes=</span><span class="ot">FALSE</span>, <span class="dt">main=</span><span class="st">"replicated"</span>)
<span class="kw">points</span>(year[yrep[,r]==<span class="dv">0</span>], nestnamenum[yrep[,r]==<span class="dv">0</span>], <span class="dt">cex=</span><span class="fl">0.3</span>) <span class="co"># add zeros </span>
}
<span class="kw">mtext</span>(<span class="st">"year"</span>, <span class="dt">outer=</span><span class="ot">TRUE</span>, <span class="dt">side=</span><span class="dv">1</span>)
<span class="kw">mtext</span>(<span class="st">"nest"</span>, <span class="dt">outer=</span><span class="ot">TRUE</span>, <span class="dt">side=</span><span class="dv">2</span>)</code></pre></div>
<p><img src="images/05-C/unnamed-chunk-9-1.png" width="672" /></p>
<p><em>Proportion of zeros for the observed data</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">mean</span>(y==<span class="dv">0</span>)</code></pre></div>
<pre><code>## [1] 0.3584071</code></pre>
<p><em>Proportion of zeros for the nsim simulated data</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">propzeros <-<span class="st"> </span><span class="kw">apply</span>(yrep, <span class="dv">2</span>, function(x) <span class="kw">mean</span>(x==<span class="dv">0</span>)) <span class="co"># </span>
<span class="kw">quantile</span>(propzeros, <span class="dt">prob=</span><span class="kw">c</span>(<span class="fl">0.025</span>, <span class="fl">0.5</span>, <span class="fl">0.975</span>))</code></pre></div>
<pre><code>## 2.5% 50% 97.5%
## 0.3238938 0.3610619 0.3978097</code></pre>
<p><em>Examine results</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">apply</span>(modsims$a, <span class="dv">2</span>, quantile, <span class="dt">prob=</span><span class="kw">c</span>(<span class="fl">0.025</span>, <span class="fl">0.5</span>, <span class="fl">0.957</span>))</code></pre></div>
<pre><code>##
## [,1] [,2]
## 2.5% -1.2574025 0.3574230
## 50% -1.0082935 0.5648238
## 95.7% -0.8142124 0.7477373</code></pre>
<p>Recall that <code>a</code> represents number of zero-survivor nests. Hence, the positive slope indicates that the number of failed nests increased over the years.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">apply</span>(modsims$b, <span class="dv">2</span>, quantile, <span class="dt">prob=</span><span class="kw">c</span>(<span class="fl">0.025</span>, <span class="fl">0.5</span>, <span class="fl">0.957</span>))</code></pre></div>
<pre><code>##
## [,1] [,2]
## 2.5% 0.8408008 -0.09805964
## 50% 0.8936011 -0.05476657
## 95.7% 0.9382964 -0.01632295</code></pre>
<p>And the number of nestlings in successful nests (indicated by <code>b</code>) decreased over time, per the estimated slope of the Poisson model’s regression coefficient (-0.05).</p>
<p><br></p>
<p><strong>Visualize results</strong></p>
<p>Three regression lines will help us understand the results: (1) the proportion of successful nests, (2) the number of nestlings in successful nests, (3) the mean number of nestlings. To draw these lines, we need to create a new data frame with predictor variable “year.” Before fitting this model, we transform this variable and create a new matrix.</p>
<p>We then extract the posterior distributions’ means for the Bernoulli coefficients (“ahat”) and Poisson model (“bhat”). To calculate the estimated proportions of surviving nests, we subtract the Bernoulli model’s fitted values from 1. (Recall that <code>zit = 1</code> denotes a dead nest.) Add the fitted value from the Poisson model, the average number of successful nests’ nestlings. We multiply the fraction of surviving nests by the mean number of nestlings for successful nests to yield the average number of nestlings per year, averaged over suriving and failed nests.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">newdat <-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">year=</span><span class="dv">1979</span>:<span class="dv">2010</span>)
newdat$year.z <-<span class="st"> </span>(newdat$year -<span class="st"> </span><span class="kw">mean</span>(dat$year))/<span class="kw">sd</span>(dat$year)
Xmat <-<span class="st"> </span><span class="kw">model.matrix</span>(~year.z, <span class="dt">data=</span>newdat)
ahat <-<span class="st"> </span><span class="kw">apply</span>(modsims$a, <span class="dv">2</span>, mean) <span class="co"># extract the estimates for a</span>
bhat <-<span class="st"> </span><span class="kw">apply</span>(modsims$b, <span class="dv">2</span>, mean) <span class="co"># extract the estimates for b</span>
newdat$propsfit <-<span class="st"> </span><span class="dv">1</span>-<span class="kw">plogis</span>(Xmat%*%ahat) <span class="co"># Proportion of nests that survived</span>
newdat$nnestfit <-<span class="st"> </span><span class="kw">exp</span>(Xmat%*%bhat)
newdat$avnnestfit<-newdat$propsfit*newdat$nnestfit</code></pre></div>
<p>Using <code>nsim</code>, we can repeat the preceding calculations with each set of model parameters–a novel MCMC iteration. And we bound every fitted value using the 2.5% and 97.5% quantiles as the lower and upper limits of a 95% CrI.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">nsim <-<span class="st"> </span><span class="kw">length</span>(modsims$lp_)
propsmat <-<span class="st"> </span><span class="kw">matrix</span>(<span class="dt">ncol=</span>nsim, <span class="dt">nrow=</span><span class="kw">nrow</span>(newdat))
nnestmat <-<span class="st"> </span><span class="kw">matrix</span>(<span class="dt">ncol=</span>nsim, <span class="dt">nrow=</span><span class="kw">nrow</span>(newdat))
avnnestmat <-<span class="st"> </span><span class="kw">matrix</span>(<span class="dt">ncol=</span>nsim, <span class="dt">nrow=</span><span class="kw">nrow</span>(newdat))
for(i in <span class="dv">1</span>:nsim){
propsmat[,i] <-<span class="st"> </span><span class="dv">1</span>-<span class="kw">plogis</span>(Xmat%*%modsims$a[i,])
nnestmat[,i] <-<span class="st"> </span><span class="kw">exp</span>(Xmat%*%modsims$b[i,])
avnnestmat[,i] <-<span class="st"> </span>propsmat[,i]*nnestmat[,i]
}
newdat$propslwr <-<span class="st"> </span><span class="kw">apply</span>(propsmat, <span class="dv">1</span>, quantile, <span class="dt">prob=</span><span class="fl">0.025</span>)
newdat$propsupr <-<span class="st"> </span><span class="kw">apply</span>(propsmat, <span class="dv">1</span>, quantile, <span class="dt">prob=</span><span class="fl">0.975</span>)
newdat$nnestlwr <-<span class="st"> </span><span class="kw">apply</span>(nnestmat, <span class="dv">1</span>, quantile, <span class="dt">prob=</span><span class="fl">0.025</span>)
newdat$nnestupr <-<span class="st"> </span><span class="kw">apply</span>(nnestmat, <span class="dv">1</span>, quantile, <span class="dt">prob=</span><span class="fl">0.975</span>)
newdat$avnnestfit <-<span class="st"> </span>newdat$propsfit*newdat$nnestfit
newdat$avnnestlwr <-<span class="st"> </span><span class="kw">apply</span>(avnnestmat, <span class="dv">1</span>, quantile, <span class="dt">prob=</span><span class="fl">0.025</span>)
newdat$avnnestupr <-<span class="st"> </span><span class="kw">apply</span>(avnnestmat, <span class="dv">1</span>, quantile, <span class="dt">prob=</span><span class="fl">0.975</span>)</code></pre></div>
<p><em>To visualize the regression lines (Figure 14.4)</em></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">par</span>(<span class="dt">mfrow=</span><span class="kw">c</span>(<span class="dv">3</span>,<span class="dv">1</span>), <span class="dt">mar=</span><span class="kw">c</span>(<span class="fl">1.5</span>,<span class="fl">5.5</span>,<span class="fl">0.1</span>,<span class="dv">1</span>), <span class="dt">oma=</span><span class="kw">c</span>(<span class="dv">3</span>,<span class="dv">0</span>,<span class="dv">0</span>,<span class="dv">0</span>))
<span class="co"># a) Number of nestlings</span>
<span class="kw">plot</span>(dat$year, dat$njuvs, <span class="dt">xlab=</span><span class="st">"Year"</span>, <span class="dt">type=</span><span class="st">"n"</span>, <span class="dt">las=</span><span class="dv">1</span>,
<span class="dt">ylab=</span><span class="st">"Number of nestlings"</span>, <span class="dt">xaxt=</span><span class="st">"n"</span>,
<span class="dt">cex.axis=</span><span class="fl">1.2</span>, <span class="dt">cex.lab=</span><span class="fl">1.4</span>, <span class="dt">ylim=</span><span class="kw">c</span>(-<span class="fl">0.4</span>, <span class="fl">5.4</span>))
<span class="kw">polygon</span>(<span class="kw">c</span>(newdat$year, <span class="kw">rev</span>(newdat$year)),
<span class="kw">c</span>(newdat$avnnestlwr, <span class="kw">rev</span>(newdat$avnnestupr)),
<span class="dt">border=</span><span class="ot">NA</span>, <span class="dt">col=</span><span class="kw">grey</span>(<span class="fl">0.5</span>))
<span class="kw">lines</span>(newdat$year, newdat$avnnestfit, <span class="dt">lwd=</span><span class="dv">2</span>)
<span class="kw">points</span>(<span class="kw">jitter</span>(dat$year), <span class="kw">jitter</span>(dat$njuvs), <span class="dt">cex=</span><span class="fl">0.5</span>)
<span class="co"># b) Number of Nestlings in successful nests</span>
<span class="kw">plot</span>(dat$year, dat$njuvs, <span class="dt">xlab=</span><span class="st">"Year"</span>, <span class="dt">type=</span><span class="st">"n"</span>,<span class="dt">las=</span><span class="dv">1</span>,
<span class="dt">ylab=</span><span class="st">"Number of nestlings</span><span class="ch">\n</span><span class="st">in successful nests"</span>, <span class="dt">xaxt=</span><span class="st">"n"</span>,
<span class="dt">cex.axis=</span><span class="fl">1.2</span>, <span class="dt">cex.lab=</span><span class="fl">1.4</span>)
<span class="kw">polygon</span>(<span class="kw">c</span>(newdat$year, <span class="kw">rev</span>(newdat$year)), <span class="kw">c</span>(newdat$nnestlwr,
<span class="kw">rev</span>(newdat$nnestupr)),
<span class="dt">border=</span><span class="ot">NA</span>, <span class="dt">col=</span><span class="kw">grey</span>(<span class="fl">0.5</span>))
<span class="kw">lines</span>(newdat$year, newdat$nnestfit, <span class="dt">lwd=</span><span class="dv">2</span>)
<span class="co"># c) Proportion of successful nests</span>
<span class="kw">plot</span>(dat$year, <span class="kw">seq</span>(<span class="dv">0</span>,<span class="dv">1</span>, <span class="dt">length=</span><span class="kw">nrow</span>(dat)), <span class="dt">xlab=</span><span class="st">"Year"</span>, <span class="dt">type=</span><span class="st">"n"</span>,
<span class="dt">ylab=</span><span class="st">"Proportion of</span><span class="ch">\n</span><span class="st">successful nests"</span>,<span class="dt">las=</span><span class="dv">1</span>,
<span class="dt">cex.axis=</span><span class="fl">1.2</span>, <span class="dt">cex.lab=</span><span class="fl">1.4</span>)
<span class="kw">polygon</span>(<span class="kw">c</span>(newdat$year, <span class="kw">rev</span>(newdat$year)), <span class="kw">c</span>(newdat$propslwr,
<span class="kw">rev</span>(newdat$propsupr)),
<span class="dt">border=</span><span class="ot">NA</span>, <span class="dt">col=</span><span class="kw">grey</span>(<span class="fl">0.5</span>))
<span class="kw">lines</span>(newdat$year, newdat$propsfit, <span class="dt">lwd=</span><span class="dv">2</span>)
<span class="kw">mtext</span>(<span class="st">"Year"</span>, <span class="dt">outer=</span><span class="ot">TRUE</span>, <span class="dt">side=</span><span class="dv">1</span>, <span class="dt">line=</span><span class="dv">2</span>, <span class="dt">cex=</span><span class="fl">1.2</span>)</code></pre></div>
<p><img src="images/05-C/unnamed-chunk-16-1.png" width="672" /></p>
<p><br></p>
</div>
<div id="discussion-questions" class="section level3">
<h3>Discussion Questions</h3>
<ul>
<li>Cold pizza from Nicole’s discussion Wednesday: What’s the golden rule when it comes to choosing optimal priors? They should not be (1) too strong, so that they are not providing us with additional information (not found in the dataset - too much overlap between prior and posterior), but also not too weak (2) to not influence the results due to lack of prior knowledge. Is there a better way than exploring different options via trail and error?</li>
<li>What are potential shortcomings of using the zero-inflated model, for instance, in the case of having a high correlation of the model parameters between the zero model and the count model? How do you address such problems?</li>
<li>What other information would be helpful to validate your zero-inflated model?</li>
<li>What’s it all about, Alfie?</li>
</ul>
<div id="updated-model-code" class="section level4">
<h4>Updated model code</h4>
<p>The 05-C.stan used old functions that will be removed from the next version of Stan (<code>increment_log_prob()</code>). The code file 05-C-updated.stan replaces these functions with their new version.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Note: Compile and run the model for a few iterations to make sure everything is working.</span>
mod <-<span class="st"> </span><span class="kw">stan</span>(<span class="dt">file =</span> <span class="st">"code/05-C-updated.stan"</span>, <span class="dt">data =</span> datax, <span class="dt">chains =</span> <span class="dv">3</span>, <span class="dt">iter =</span> <span class="dv">10</span>)</code></pre></div>
<pre><code>## In file included from C:/Program Files/R/R-3.3.1/library/BH/include/boost/config.hpp:39:0,
## from C:/Program Files/R/R-3.3.1/library/BH/include/boost/math/tools/config.hpp:13,
## from C:/Program Files/R/R-3.3.1/library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Program Files/R/R-3.3.1/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Program Files/R/R-3.3.1/library/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Program Files/R/R-3.3.1/library/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Program Files/R/R-3.3.1/library/StanHeaders/include/stan/math.hpp:4,
## from C:/Program Files/R/R-3.3.1/library/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file2c4c268f9f.cpp:8:
## C:/Program Files/R/R-3.3.1/library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL '05-C-updated' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 10 [ 10%] (Warmup)
## Iteration: 2 / 10 [ 20%] (Warmup)
## Iteration: 3 / 10 [ 30%] (Warmup)
## Iteration: 4 / 10 [ 40%] (Warmup)
## Iteration: 5 / 10 [ 50%] (Warmup)
## Iteration: 6 / 10 [ 60%] (Sampling)
## Iteration: 7 / 10 [ 70%] (Sampling)
## Iteration: 8 / 10 [ 80%] (Sampling)
## Iteration: 9 / 10 [ 90%] (Sampling)
## Iteration: 10 / 10 [100%] (Sampling)
##
## Elapsed Time: 0.084 seconds (Warm-up)
## 0.007 seconds (Sampling)
## 0.091 seconds (Total)
##
##
## SAMPLING FOR MODEL '05-C-updated' NOW (CHAIN 2).
##
## Gradient evaluation took 0.001 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 10 [ 10%] (Warmup)
## Iteration: 2 / 10 [ 20%] (Warmup)
## Iteration: 3 / 10 [ 30%] (Warmup)
## Iteration: 4 / 10 [ 40%] (Warmup)
## Iteration: 5 / 10 [ 50%] (Warmup)
## Iteration: 6 / 10 [ 60%] (Sampling)
## Iteration: 7 / 10 [ 70%] (Sampling)
## Iteration: 8 / 10 [ 80%] (Sampling)
## Iteration: 9 / 10 [ 90%] (Sampling)
## Iteration: 10 / 10 [100%] (Sampling)
##
## Elapsed Time: 0.009 seconds (Warm-up)
## 0.005 seconds (Sampling)
## 0.014 seconds (Total)
##
##
## SAMPLING FOR MODEL '05-C-updated' NOW (CHAIN 3).
##
## Gradient evaluation took 0.001 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 10 [ 10%] (Warmup)
## Iteration: 2 / 10 [ 20%] (Warmup)
## Iteration: 3 / 10 [ 30%] (Warmup)
## Iteration: 4 / 10 [ 40%] (Warmup)
## Iteration: 5 / 10 [ 50%] (Warmup)
## Iteration: 6 / 10 [ 60%] (Sampling)
## Iteration: 7 / 10 [ 70%] (Sampling)
## Iteration: 8 / 10 [ 80%] (Sampling)
## Iteration: 9 / 10 [ 90%] (Sampling)
## Iteration: 10 / 10 [100%] (Sampling)
##
## Elapsed Time: 0.017 seconds (Warm-up)
## 0.01 seconds (Sampling)
## 0.027 seconds (Total)</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Then run the same model for more iterations.</span>
<span class="co">#modfit <- stan(fit = mod, iter = 1000, chains = 3)</span></code></pre></div>
</div>
</div>
</div>
</div>
</div>
<script>
// add bootstrap table styles to pandoc tables
function bootstrapStylePandocTables() {
$('tr.header').parent('thead').parent('table').addClass('table table-condensed');
}
$(document).ready(function () {
bootstrapStylePandocTables();
});
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
script.src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>