-
Notifications
You must be signed in to change notification settings - Fork 0
/
index.html
6361 lines (6331 loc) · 736 KB
/
index.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
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.2.269">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<meta name="author" content="MH Manuel Haqiqatkhah">
<title>Skewness and staging: Does the floor effect induce bias in multilevel AR(1) models?</title>
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1.6em;
vertical-align: middle;
}
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { display: inline-block; line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
pre.numberSource code > span
{ position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
{ content: counter(source-line);
position: relative; left: -1em; text-align: right; vertical-align: baseline;
border: none; display: inline-block;
-webkit-touch-callout: none; -webkit-user-select: none;
-khtml-user-select: none; -moz-user-select: none;
-ms-user-select: none; user-select: none;
padding: 0 4px; width: 4em;
color: #aaaaaa;
}
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa; padding-left: 4px; }
div.sourceCode
{ }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { color: #008000; } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { color: #008000; font-weight: bold; } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
div.csl-bib-body { }
div.csl-entry {
clear: both;
}
.hanging div.csl-entry {
margin-left:2em;
text-indent:-2em;
}
div.csl-left-margin {
min-width:2em;
float:left;
}
div.csl-right-inline {
margin-left:2em;
padding-left:1em;
}
div.csl-indent {
margin-left: 2em;
}
</style>
<script src="Code documentations_files/libs/clipboard/clipboard.min.js"></script>
<script src="Code documentations_files/libs/quarto-html/quarto.js"></script>
<script src="Code documentations_files/libs/quarto-html/popper.min.js"></script>
<script src="Code documentations_files/libs/quarto-html/tippy.umd.min.js"></script>
<script src="Code documentations_files/libs/quarto-html/anchor.min.js"></script>
<link href="Code documentations_files/libs/quarto-html/tippy.css" rel="stylesheet">
<link href="Code documentations_files/libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" class="quarto-color-scheme" id="quarto-text-highlighting-styles">
<link href="Code documentations_files/libs/quarto-html/quarto-syntax-highlighting-dark.css" rel="prefetch" class="quarto-color-scheme quarto-color-alternate" id="quarto-text-highlighting-styles">
<script src="Code documentations_files/libs/bootstrap/bootstrap.min.js"></script>
<link href="Code documentations_files/libs/bootstrap/bootstrap-icons.css" rel="stylesheet">
<link href="Code documentations_files/libs/bootstrap/bootstrap.min.css" rel="stylesheet" class="quarto-color-scheme" id="quarto-bootstrap" data-mode="light">
<link href="Code documentations_files/libs/bootstrap/bootstrap-dark.min.css" rel="prefetch" class="quarto-color-scheme quarto-color-alternate" id="quarto-bootstrap" data-mode="dark">
<link href="Code documentations_files/libs/quarto-contrib/fontawesome6-0.1.0/all.css" rel="stylesheet">
<link href="Code documentations_files/libs/quarto-contrib/fontawesome6-0.1.0/latex-fontsize.css" rel="stylesheet">
<script src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-chtml-full.js" type="text/javascript"></script>
</head>
<body>
<div id="quarto-content" class="page-columns page-rows-contents page-layout-article toc-left">
<div id="quarto-sidebar-toc-left" class="sidebar toc-left">
<nav id="TOC" role="doc-toc" class="toc-active">
<h2 id="toc-title">Table of contents</h2>
<ul>
<li><a href="#sec-introduction" id="toc-sec-introduction" class="nav-link active" data-scroll-target="#sec-introduction"><span class="toc-section-number">0</span> Introduction</a></li>
<li><a href="#sec-components" id="toc-sec-components" class="nav-link" data-scroll-target="#sec-components"><span class="toc-section-number">1</span> Core components</a>
<ul>
<li><a href="#sec-simulation-component" id="toc-sec-simulation-component" class="nav-link" data-scroll-target="#sec-simulation-component"><span class="toc-section-number">1.1</span> The <em>Simulation</em> component</a>
<ul>
<li><a href="#sec-dgm-specifications" id="toc-sec-dgm-specifications" class="nav-link" data-scroll-target="#sec-dgm-specifications"><span class="toc-section-number">1.1.1</span> Data-generating models specifications</a></li>
<li><a href="#sec-general-dgm-wrappers" id="toc-sec-general-dgm-wrappers" class="nav-link" data-scroll-target="#sec-general-dgm-wrappers"><span class="toc-section-number">1.1.2</span> General DGM wrappers</a></li>
<li><a href="#sec-dataset-generation" id="toc-sec-dataset-generation" class="nav-link" data-scroll-target="#sec-dataset-generation"><span class="toc-section-number">1.1.3</span> Dataset generation</a></li>
</ul></li>
<li><a href="#sec-analysis-component" id="toc-sec-analysis-component" class="nav-link" data-scroll-target="#sec-analysis-component"><span class="toc-section-number">1.2</span> The <em>Analysis</em> component</a></li>
<li><a href="#sec-harvesting-component" id="toc-sec-harvesting-component" class="nav-link" data-scroll-target="#sec-harvesting-component"><span class="toc-section-number">1.3</span> The <em>Harvesting</em> component</a></li>
<li><a href="#sec-reporting-component" id="toc-sec-reporting-component" class="nav-link" data-scroll-target="#sec-reporting-component"><span class="toc-section-number">1.4</span> The <em>Reporting</em> component</a>
<ul>
<li><a href="#sec-report-clean-harvest" id="toc-sec-report-clean-harvest" class="nav-link" data-scroll-target="#sec-report-clean-harvest"><span class="toc-section-number">1.4.1</span> Cleaning the harvested dataset</a></li>
<li><a href="#sec-make-distribution-plots" id="toc-sec-make-distribution-plots" class="nav-link" data-scroll-target="#sec-make-distribution-plots"><span class="toc-section-number">1.4.2</span> Making distribution plots of simulated datasets</a></li>
<li><a href="#sec-report-estimation-plots" id="toc-sec-report-estimation-plots" class="nav-link" data-scroll-target="#sec-report-estimation-plots"><span class="toc-section-number">1.4.3</span> Making plots of estimation results</a></li>
<li><a href="#sec-report-caterpillar-plots" id="toc-sec-report-caterpillar-plots" class="nav-link" data-scroll-target="#sec-report-caterpillar-plots"><span class="toc-section-number">1.4.4</span> Making plots of parameter estimates coverage</a></li>
<li><a href="#making-plots-for-dgm-time-series" id="toc-making-plots-for-dgm-time-series" class="nav-link" data-scroll-target="#making-plots-for-dgm-time-series"><span class="toc-section-number">1.4.5</span> Making plots for DGM time series</a></li>
</ul></li>
</ul></li>
<li><a href="#sec-pipeline" id="toc-sec-pipeline" class="nav-link" data-scroll-target="#sec-pipeline"><span class="toc-section-number">2</span> Pipeline</a>
<ul>
<li><a href="#sec-bookkeeping" id="toc-sec-bookkeeping" class="nav-link" data-scroll-target="#sec-bookkeeping"><span class="toc-section-number">2.1</span> Book-keeping</a></li>
<li><a href="#sec-pipeline-functions" id="toc-sec-pipeline-functions" class="nav-link" data-scroll-target="#sec-pipeline-functions"><span class="toc-section-number">2.2</span> Pipeline functions</a></li>
<li><a href="#sec-pipeline-code" id="toc-sec-pipeline-code" class="nav-link" data-scroll-target="#sec-pipeline-code"><span class="toc-section-number">2.3</span> Pipeline code</a></li>
</ul></li>
<li><a href="#sec-simulation-results" id="toc-sec-simulation-results" class="nav-link" data-scroll-target="#sec-simulation-results"><span class="toc-section-number">3</span> Simulation results</a>
<ul>
<li><a href="#plots-main-results" id="toc-plots-main-results" class="nav-link" data-scroll-target="#plots-main-results"><span class="toc-section-number">3.1</span> Plots main results</a></li>
<li><a href="#plots-of-mcmc-parameter-estimates-coverage" id="toc-plots-of-mcmc-parameter-estimates-coverage" class="nav-link" data-scroll-target="#plots-of-mcmc-parameter-estimates-coverage"><span class="toc-section-number">3.2</span> Plots of MCMC parameter estimates coverage</a></li>
</ul></li>
<li><a href="#sec-additional-analyses-and-figures" id="toc-sec-additional-analyses-and-figures" class="nav-link" data-scroll-target="#sec-additional-analyses-and-figures"><span class="toc-section-number">4</span> Additional figures and analyses</a>
<ul>
<li><a href="#sec-additional-DGM-time-series-plots" id="toc-sec-additional-DGM-time-series-plots" class="nav-link" data-scroll-target="#sec-additional-DGM-time-series-plots"><span class="toc-section-number">4.1</span> DGM time series plots</a></li>
<li><a href="#sec-profiles-of-simulated-datasets" id="toc-sec-profiles-of-simulated-datasets" class="nav-link" data-scroll-target="#sec-profiles-of-simulated-datasets"><span class="toc-section-number">4.2</span> Profiles of the simulated datasets</a></li>
<li><a href="#sec-cogito-data-analysis" id="toc-sec-cogito-data-analysis" class="nav-link" data-scroll-target="#sec-cogito-data-analysis"><span class="toc-section-number">4.3</span> COGITO data analysis</a></li>
</ul></li>
</ul>
</nav>
</div>
<div id="quarto-margin-sidebar" class="sidebar margin-sidebar">
</div>
<main class="content" id="quarto-document-content">
<header id="title-block-header" class="quarto-title-block default">
<div class="quarto-title">
<div class="quarto-title-block"><div><h1 class="title">Skewness and staging: Does the floor effect induce bias in multilevel AR(1) models?</h1><button type="button" class="btn code-tools-button dropdown-toggle" id="quarto-code-tools-menu" data-bs-toggle="dropdown" aria-expanded="false"><i class="bi"></i> Code</button><ul class="dropdown-menu dropdown-menu-end" aria-labelelledby="quarto-code-tools-menu"><li><a id="quarto-show-all-code" class="dropdown-item" href="javascript:void(0)" role="button">Show All Code</a></li><li><a id="quarto-hide-all-code" class="dropdown-item" href="javascript:void(0)" role="button">Hide All Code</a></li><li><hr class="dropdown-divider"></li><li><a id="quarto-view-source" class="dropdown-item" href="javascript:void(0)" role="button">View Source</a></li></ul></div></div>
<p class="subtitle lead">Reproducible codes</p>
</div>
<div class="quarto-title-meta">
<div>
<div class="quarto-title-meta-heading">Author</div>
<div class="quarto-title-meta-contents">
<p>MH Manuel Haqiqatkhah </p>
</div>
</div>
</div>
</header>
<section id="sec-introduction" class="level1" data-number="0">
<h1 data-number="0"><span class="header-section-number">0</span> Introduction</h1>
<p>This document contains the reproducible code for the manuscript <em>Skewness and staging: Does the floor effect induce bias in multilevel AR(1) models?</em> by <span class="citation" data-cites="haqiqatkhah_2022_SkewnessStagingDoes">Haqiqatkhah, Ryan, and Hamaker (<a href="#ref-haqiqatkhah_2022_SkewnessStagingDoes" role="doc-biblioref">2022</a>)</span>. Please cite this document with the following info:</p>
<blockquote class="blockquote">
<p>Haqiqatkhah, M. M., Ryan, O., & Hamaker, E. L. (2022). <em>Skewness and staging: Does the floor effect induce bias in multilevel AR(1) models?</em> PsyArXiv. <a href="https://doi.org/10.31234/osf.io/myuvr" class="uri">https://doi.org/10.31234/osf.io/myuvr</a></p>
</blockquote>
<p>In this study, we simulated multilevel data from three data generating mechanisms (DGMs), namely, the AR(1), <span class="math inline">\(\chi^2\)</span>AR(1), BinAR(1), and PoDAR(1) models with different parameter sets. The simulation was conducted using the following modular pipeline design, inspired by Bien’s R package <code>simulator</code> <span class="citation" data-cites="bien_2016_SimulatorEngineStreamline">(<a href="#ref-bien_2016_SimulatorEngineStreamline" role="doc-biblioref">2016</a>)</span>, consisting of the following <strong>components</strong>:</p>
<ol type="A">
<li><em>Simulation</em>: generating the datasets</li>
<li><em>Analysis</em>: modeling the data</li>
<li><em>Harvesting</em>: collecting the relevant parameter estimates</li>
<li><em>Reporting</em>: making tables and plots</li>
</ol>
<p>And the components were placed in a <strong>pipeline</strong>, that managed:</p>
<ol type="1">
<li>Making the simulation design matrix that include all relevant conditions</li>
<li>Book-keeping data files belonging to each replication of each condition</li>
<li>Performing simulations in batch</li>
<li>Performing Analyses in batch</li>
<li>Collecting the data in batch</li>
</ol>
<p>This document is structured as follows. In <a href="#sec-components">Section 1</a>, we explain the four components and the functions used therein. Then, in <a href="#sec-pipeline">Section 2</a> we explain the wrapper functions used in the pipeline, and show how the pipeline was—and can be—executed. Then, in <a href="#sec-simulation-results">Section 3</a>, we discuss how the harvested data was used to make the figures used in the paper (and others that were not included). Finally, in <a href="#sec-additional-analyses-and-figures">Section 4</a> we discuss how supplementary plots were made and how the empirical data analysis [on COGITO data; <span class="citation" data-cites="schmiedek_2010_HundredDaysCognitive">Schmiedek, Lövdén, and Lindenberger (<a href="#ref-schmiedek_2010_HundredDaysCognitive" role="doc-biblioref">2010</a>)</span>] was done. Note that although the codes provided here are cleaned as much as possible, they are not necessarily succinctly written; some functions were written to accommodate the most general functionalities which turned out to be not necessary for the simulation study.</p>
<p>The <a href="https://github.com/psyguy/skewness-staging">GitHub <i class="fa-brands fa-github" aria-label="github"></i> repository of this study</a> contains all of the code necessary to run the study (in the <code>scripts</code> folder <a href="https://github.com/psyguy/skewness-staging/tree/main/scripts"><i class="fa-solid fa-globe" aria-label="globe"></i></a>). Furthermore, because running the whole study would take a huge amount of time (it took us more than 88 days on a 24-core server and produced 2.69 TB of data), the the raw and summarized datasets of the estimated parameters are provided in <code>simulation-files/harvests</code> <a href="https://github.com/psyguy/skewness-staging/tree/main/simulation-files/harvest-files"><i class="fa-solid fa-globe" aria-label="globe"></i></a>. To give an impression of the simulated datasets and the analysis output files, the <code>.rds</code> files of the first replication of the study are also provided in <code>simulation-files</code> folder <a href="https://github.com/psyguy/skewness-staging/tree/main/simulation-files"><i class="fa-solid fa-globe" aria-label="globe"></i></a>, and the <em>M</em>plus files of the first replication of the <span class="math inline">\(N=100, T=100\)</span> condition is provided in the <code>Mplus-files</code> folder <a href="https://github.com/psyguy/skewness-staging/tree/main/Mplus-files"><i class="fa-solid fa-globe" aria-label="globe"></i></a>. Finally, all the figures that can be generated based on the study results are provided in the <code>figures</code> folder <a href="https://github.com/psyguy/skewness-staging/tree/main/figures"><i class="fa-solid fa-globe" aria-label="globe"></i></a>. The figures are also referenced individually in <a href="#sec-simulation-results">Section 3</a> and <a href="#sec-additional-analyses-and-figures">Section 4</a> of this document.</p>
<p>To replicate the study from the scratch, you should first either clone the repository (using <code>git clone https://github.com/psyguy/skewness-staging.git</code>) or <a href="https://github.com/psyguy/skewness-staging/archive/refs/heads/main.zip">download the repository as a zip file <i class="fa-solid fa-download" aria-label="download"></i></a> and extract it on your machine. Then you can sequentially run the <code>.R</code> files you find in the <code>scripts</code> folder. Note that you would need to have a licensed version of <em>M</em>plus version 8.6 <span class="citation" data-cites="muthen_2017_MplusUserGuide">(<a href="#ref-muthen_2017_MplusUserGuide" role="doc-biblioref">Muthén and Muthén 2017</a>)</span> on your machine to run the whole study. Instead of running the scripts separately, if you have <a href="https://quarto.org/docs/get-started/">Quarto installed <i class="fa-solid fa-globe" aria-label="globe"></i></a>, you can also compile <code>Code documentations.qmd</code> located in the root directory using Quarto after setting the following variables to <code>TRUE</code>:</p>
<div class="cell" data-hash="Code-documentations_cache/html/unnamed-chunk-1_9371825e8b9a4f2a95545c5bbf41783c">
<div class="sourceCode cell-code" id="cb1"><pre class="sourceCode numberSource r number-lines code-with-copy"><code class="sourceCode r"><span id="cb1-1"><a href="#cb1-1"></a>study_RunPipeline <span class="ot"><-</span> <span class="cn">FALSE</span></span>
<span id="cb1-2"><a href="#cb1-2"></a>study_MakeMainPlots <span class="ot"><-</span> <span class="cn">FALSE</span></span>
<span id="cb1-3"><a href="#cb1-3"></a>study_MakeAdditionalPlots <span class="ot"><-</span> <span class="cn">FALSE</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Finally, in case you want to change the scripts (e.g., to run a smaller portion of the simulation, or try other parameters, etc.), you should look up the <code>kintr</code> parameters called in each chunk (with <code><<some_param>></code>) and find the corresponding code (under <code>## @knitr some_param</code>) in the <code>scripts</code> folder.</p>
</section>
<section id="sec-components" class="level1" data-number="1">
<h1 data-number="1"><span class="header-section-number">1</span> Core components</h1>
<section id="sec-simulation-component" class="level2" data-number="1.1">
<h2 data-number="1.1" class="anchored" data-anchor-id="sec-simulation-component"><span class="header-section-number">1.1</span> The <em>Simulation</em> component</h2>
<p>The Simulation component consists of three sets of functions:</p>
<ol type="i">
<li>Functions that implement the DGMs and generate univariate (<span class="math inline">\(N=1\)</span>) time series of length <span class="math inline">\(T\)</span> from the parameters given to them;</li>
<li>Wrappers that interface the DGM functions;</li>
<li>A wrapper to generate datasets (consisting on <span class="math inline">\(N\)</span> time series of length <span class="math inline">\(T\)</span>) with a given DGM</li>
</ol>
<section id="sec-dgm-specifications" class="level3" data-number="1.1.1">
<h3 data-number="1.1.1" class="anchored" data-anchor-id="sec-dgm-specifications"><span class="header-section-number">1.1.1</span> Data-generating models specifications</h3>
<p>First we define functions for each data-generating models (DGMs) that can produce univariate, single-subject (<span class="math inline">\(N=1\)</span>) time series of desired length <span class="math inline">\(T\)</span> (default: <code>T = 100</code>) with the two canonical parameters and a given random seed (default: <code>seed = 0</code>). All model(-implied) parameters are saved in a list (called <code>pa</code>).</p>
<p>For each model, the first observation (<span class="math inline">\(X_1\)</span>) is randomly drawn from the model-implied marginal distribution, to eliminate the need for removing the burn-in window in the beginning of the data. After the data is generated, in case the argument <code>only.ts</code> is set to be <code>TRUE</code>, the raw data (as a vector of length <code>T</code>) is returned. Otherwise, the function calculates empirical dynamic (<span class="math inline">\(\phi\)</span>) and marginal (<span class="math inline">\(\mu\)</span>, <span class="math inline">\(\sigma^2\)</span>, and <span class="math inline">\(\gamma\)</span>) parameters based on the simulated data, and save it in a list (<code>Empirical.Parameters</code>). Furthermore, two <span class="math inline">\(\mathrm{\LaTeX{}}\)</span>-ready strings (<code>Model.Description</code> and <code>Model.Description.Short</code>) are made which include a summary of the model parameters (that can be used, e.g., in plots). Finally, in case <code>only.ts != TRUE</code>, the function returns a list consisting of the time series (stored in <code>x</code>), verbal description of the dataset (<code>Model.Description</code> and <code>Model.Description.Short</code>), theoretical (i.e., model-implied) parameters (<code>Model.Parameters</code>), and empirical (i.e., sample) estimated parameters (<code>Empirical.Parameters</code>).</p>
<div class="panel-tabset">
<ul class="nav nav-tabs" role="tablist"><li class="nav-item" role="presentation"><a class="nav-link active" id="tabset-1-1-tab" data-bs-toggle="tab" data-bs-target="#tabset-1-1" role="tab" aria-controls="tabset-1-1" aria-selected="true" href="">The AR(1) model</a></li><li class="nav-item" role="presentation"><a class="nav-link" id="tabset-1-2-tab" data-bs-toggle="tab" data-bs-target="#tabset-1-2" role="tab" aria-controls="tabset-1-2" aria-selected="false" href="">The <span class="math inline">\(\chi^2\)</span>AR(1) model</a></li><li class="nav-item" role="presentation"><a class="nav-link" id="tabset-1-3-tab" data-bs-toggle="tab" data-bs-target="#tabset-1-3" role="tab" aria-controls="tabset-1-3" aria-selected="false" href="">The BinAR(1) model</a></li><li class="nav-item" role="presentation"><a class="nav-link" id="tabset-1-4-tab" data-bs-toggle="tab" data-bs-target="#tabset-1-4" role="tab" aria-controls="tabset-1-4" aria-selected="false" href="">The PoDAR(1) model</a></li></ul>
<div class="tab-content">
<div id="tabset-1-1" class="tab-pane active" role="tabpanel" aria-labelledby="tabset-1-1-tab">
<p>The canonical parameters of the AR(1) model with normally distributed residuals (which we referred to as <code>NAR(1)</code> in the simulation) are the autoregressive parameter <span class="math inline">\(\phi\)</span> (default: <code>phi = 0.4</code>), mean <span class="math inline">\(\mu\)</span> (default: <code>Mean = 50</code>), and the marginal variance <span class="math inline">\(\sigma^2\)</span> (default: <code>var.marginal = 4</code>). Based on the marginal variance, the residual variance (<code>var.resid</code>) is calculated via <span class="math inline">\(\sigma^2_\epsilon = \sigma^2 (1 - \phi^2)\)</span>.</p>
<section id="construction" class="level5" data-number="1.1.1.0.1">
<h5 data-number="1.1.1.0.1" class="anchored" data-anchor-id="construction"><span class="header-section-number">1.1.1.0.1</span> Construction</h5>
<p>The time series is constructed by first generating a zero-centered time series <span class="math inline">\(\tilde X_t\)</span> (<code>x_cent</code>). To do so, first the initial observation in the time series (<code>x_cent[1]</code>) is sampled from normal distribution with mean zero and a variance equal to the marginal variance of the model:</p>
<p><span class="math display">\[
\tilde X_1 \sim \mathcal{N}(0, \sigma^2)
\]</span></p>
<p>Then, the remainder of the time series is generated using the definition of the AR(1) model (note that here the residual variance is used in the normal distribution):</p>
<p><span class="math display">\[
\begin{aligned}
\tilde X_{t} &= \phi \tilde X_{t-1} + \epsilon_{t} \\
\epsilon_{t} &\sim \mathcal{N}(0, {\sigma^2_{\epsilon}})
\end{aligned}
\]</span></p>
<p>Finally, the mean is added to the centered zero-centered time series to reach the final time series with mean <span class="math inline">\(\mu\)</span>:</p>
<p><span class="math display">\[
X_t = \tilde X_t + \mu
\]</span></p>
</section>
<section id="code" class="level5" data-number="1.1.1.0.2">
<h5 data-number="1.1.1.0.2" class="anchored" data-anchor-id="code"><span class="header-section-number">1.1.1.0.2</span> Code</h5>
<div class="cell" data-hash="Code-documentations_cache/html/unnamed-chunk-4_9f104b3390a1abb4d1d7294dbc3e3b75">
<details>
<summary>Click to expand the code</summary>
<div class="sourceCode cell-code" id="cb2"><pre class="sourceCode numberSource r number-lines code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1"></a>dgm_nar <span class="ot"><-</span> <span class="cf">function</span>(...){</span>
<span id="cb2-2"><a href="#cb2-2"></a></span>
<span id="cb2-3"><a href="#cb2-3"></a> pa <span class="ot"><-</span> <span class="fu">list</span>(...)</span>
<span id="cb2-4"><a href="#cb2-4"></a></span>
<span id="cb2-5"><a href="#cb2-5"></a> <span class="cf">if</span>(<span class="fu">is.list</span>(pa<span class="sc">$</span>pa)) pa <span class="ot"><-</span> pa<span class="sc">$</span>pa</span>
<span id="cb2-6"><a href="#cb2-6"></a></span>
<span id="cb2-7"><a href="#cb2-7"></a> <span class="do">## setting default seed if not given</span></span>
<span id="cb2-8"><a href="#cb2-8"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>phi)) pa<span class="sc">$</span>phi <span class="ot"><-</span> <span class="fl">0.4</span></span>
<span id="cb2-9"><a href="#cb2-9"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>Mean)) pa<span class="sc">$</span>Mean <span class="ot"><-</span> <span class="dv">50</span></span>
<span id="cb2-10"><a href="#cb2-10"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>var.marginal)) pa<span class="sc">$</span>var.marginal <span class="ot"><-</span> <span class="dv">4</span></span>
<span id="cb2-11"><a href="#cb2-11"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>var.resid)) pa<span class="sc">$</span>var.resid <span class="ot"><-</span> pa<span class="sc">$</span>var.marginal <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>)</span>
<span id="cb2-12"><a href="#cb2-12"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>k)) pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb2-13"><a href="#cb2-13"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>T)) pa<span class="sc">$</span>T <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb2-14"><a href="#cb2-14"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>seed)) pa<span class="sc">$</span>seed <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb2-15"><a href="#cb2-15"></a></span>
<span id="cb2-16"><a href="#cb2-16"></a> <span class="do">## Making sure var.marginal and var.resid are correctly related</span></span>
<span id="cb2-17"><a href="#cb2-17"></a> pa<span class="sc">$</span>var.marginal <span class="ot"><-</span> pa<span class="sc">$</span>var.resid <span class="sc">/</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>)</span>
<span id="cb2-18"><a href="#cb2-18"></a></span>
<span id="cb2-19"><a href="#cb2-19"></a> <span class="do">### first make time series centered around zero</span></span>
<span id="cb2-20"><a href="#cb2-20"></a> <span class="do">## drawing the first sample x_cent_1</span></span>
<span id="cb2-21"><a href="#cb2-21"></a> <span class="fu">set.seed</span>(pa<span class="sc">$</span>seed)</span>
<span id="cb2-22"><a href="#cb2-22"></a> x_cent <span class="ot"><-</span> <span class="fu">rep</span>(<span class="cn">NA</span>, pa<span class="sc">$</span>T)</span>
<span id="cb2-23"><a href="#cb2-23"></a> x_cent[<span class="dv">1</span>] <span class="ot"><-</span> <span class="fu">rnorm</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb2-24"><a href="#cb2-24"></a> <span class="at">mean =</span> <span class="dv">0</span>,</span>
<span id="cb2-25"><a href="#cb2-25"></a> <span class="at">sd =</span> <span class="fu">sqrt</span>(pa<span class="sc">$</span>var.marginal)</span>
<span id="cb2-26"><a href="#cb2-26"></a> )</span>
<span id="cb2-27"><a href="#cb2-27"></a></span>
<span id="cb2-28"><a href="#cb2-28"></a> <span class="do">## making the rest of the centered time series</span></span>
<span id="cb2-29"><a href="#cb2-29"></a> <span class="cf">for</span> (t <span class="cf">in</span> <span class="dv">2</span><span class="sc">:</span>pa<span class="sc">$</span>T){</span>
<span id="cb2-30"><a href="#cb2-30"></a></span>
<span id="cb2-31"><a href="#cb2-31"></a> x_cent[t] <span class="ot"><-</span> pa<span class="sc">$</span>phi<span class="sc">*</span>x_cent[t<span class="dv">-1</span>] <span class="sc">+</span> <span class="fu">rnorm</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb2-32"><a href="#cb2-32"></a> <span class="at">mean =</span> <span class="dv">0</span>,</span>
<span id="cb2-33"><a href="#cb2-33"></a> <span class="at">sd =</span> <span class="fu">sqrt</span>(pa<span class="sc">$</span>var.resid)</span>
<span id="cb2-34"><a href="#cb2-34"></a> )</span>
<span id="cb2-35"><a href="#cb2-35"></a> }</span>
<span id="cb2-36"><a href="#cb2-36"></a> <span class="do">## adding the mean to the centered time series</span></span>
<span id="cb2-37"><a href="#cb2-37"></a> x <span class="ot"><-</span> x_cent <span class="sc">+</span> pa<span class="sc">$</span>Mean</span>
<span id="cb2-38"><a href="#cb2-38"></a></span>
<span id="cb2-39"><a href="#cb2-39"></a> <span class="do">## quick output of raw time series without book-keeping variables/parameters</span></span>
<span id="cb2-40"><a href="#cb2-40"></a> <span class="cf">if</span>(<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>only.ts))</span>
<span id="cb2-41"><a href="#cb2-41"></a> <span class="cf">if</span>(pa<span class="sc">$</span>only.ts<span class="sc">==</span><span class="cn">TRUE</span>) <span class="fu">return</span>(x)</span>
<span id="cb2-42"><a href="#cb2-42"></a></span>
<span id="cb2-43"><a href="#cb2-43"></a> Empirical.Parameters <span class="ot">=</span> <span class="fu">list</span>(<span class="at">Mean =</span> <span class="fu">mean</span>(x),</span>
<span id="cb2-44"><a href="#cb2-44"></a> <span class="at">Variance =</span> <span class="fu">var</span>(x),</span>
<span id="cb2-45"><a href="#cb2-45"></a> <span class="at">Skewness =</span> moments<span class="sc">::</span><span class="fu">skewness</span>(x),</span>
<span id="cb2-46"><a href="#cb2-46"></a> <span class="at">AR =</span> <span class="fu">acf</span>(x, <span class="at">lag.max =</span> <span class="dv">1</span>, <span class="at">plot =</span> <span class="cn">FALSE</span>)<span class="sc">$</span>acf[<span class="dv">2</span>]</span>
<span id="cb2-47"><a href="#cb2-47"></a> )</span>
<span id="cb2-48"><a href="#cb2-48"></a></span>
<span id="cb2-49"><a href="#cb2-49"></a> <span class="do">## making a LaTeX-ready list description of the model</span></span>
<span id="cb2-50"><a href="#cb2-50"></a> Model.Description <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"</span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb2-51"><a href="#cb2-51"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb2-52"><a href="#cb2-52"></a> <span class="st">"("</span>,</span>
<span id="cb2-53"><a href="#cb2-53"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb2-54"><a href="#cb2-54"></a> <span class="st">")"</span>,</span>
<span id="cb2-55"><a href="#cb2-55"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb2-56"><a href="#cb2-56"></a> <span class="dv">0</span>,</span>
<span id="cb2-57"><a href="#cb2-57"></a> <span class="st">"("</span>,</span>
<span id="cb2-58"><a href="#cb2-58"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb2-59"><a href="#cb2-59"></a> <span class="st">")"</span>,</span>
<span id="cb2-60"><a href="#cb2-60"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">phi = "</span>,</span>
<span id="cb2-61"><a href="#cb2-61"></a> <span class="fu">round</span>(pa<span class="sc">$</span>phi,<span class="dv">2</span>),</span>
<span id="cb2-62"><a href="#cb2-62"></a> <span class="st">"("</span>,</span>
<span id="cb2-63"><a href="#cb2-63"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>AR,<span class="dv">2</span>),</span>
<span id="cb2-64"><a href="#cb2-64"></a> <span class="st">")"</span>,</span>
<span id="cb2-65"><a href="#cb2-65"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">sigma^2_{marginal} = "</span>,</span>
<span id="cb2-66"><a href="#cb2-66"></a> <span class="fu">round</span>(pa<span class="sc">$</span>var.marginal,<span class="dv">2</span>),</span>
<span id="cb2-67"><a href="#cb2-67"></a> <span class="st">"("</span>,</span>
<span id="cb2-68"><a href="#cb2-68"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Variance,<span class="dv">2</span>),</span>
<span id="cb2-69"><a href="#cb2-69"></a> <span class="st">")"</span>,</span>
<span id="cb2-70"><a href="#cb2-70"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; T = "</span>,</span>
<span id="cb2-71"><a href="#cb2-71"></a> pa<span class="sc">$</span>T,</span>
<span id="cb2-72"><a href="#cb2-72"></a> <span class="st">"$"</span>)</span>
<span id="cb2-73"><a href="#cb2-73"></a></span>
<span id="cb2-74"><a href="#cb2-74"></a> Model.Description.Short <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"$</span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb2-75"><a href="#cb2-75"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Mean, <span class="dv">2</span>),</span>
<span id="cb2-76"><a href="#cb2-76"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">sigma^2_{</span><span class="sc">\\</span><span class="st">epsilon} = "</span>,</span>
<span id="cb2-77"><a href="#cb2-77"></a> <span class="fu">round</span>(pa<span class="sc">$</span>var.resid, <span class="dv">2</span>),</span>
<span id="cb2-78"><a href="#cb2-78"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">phi = "</span>,</span>
<span id="cb2-79"><a href="#cb2-79"></a> <span class="fu">round</span>(pa<span class="sc">$</span>phi, <span class="dv">2</span>),</span>
<span id="cb2-80"><a href="#cb2-80"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">rightarrow"</span>,</span>
<span id="cb2-81"><a href="#cb2-81"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">sigma^2 = "</span>,</span>
<span id="cb2-82"><a href="#cb2-82"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Variance,<span class="dv">2</span>),</span>
<span id="cb2-83"><a href="#cb2-83"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb2-84"><a href="#cb2-84"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb2-85"><a href="#cb2-85"></a> <span class="st">"$"</span>)</span>
<span id="cb2-86"><a href="#cb2-86"></a></span>
<span id="cb2-87"><a href="#cb2-87"></a></span>
<span id="cb2-88"><a href="#cb2-88"></a> <span class="do">## making the output object</span></span>
<span id="cb2-89"><a href="#cb2-89"></a> output <span class="ot"><-</span> <span class="fu">list</span>(<span class="at">x =</span> x,</span>
<span id="cb2-90"><a href="#cb2-90"></a> <span class="at">Model.Description =</span> Model.Description,</span>
<span id="cb2-91"><a href="#cb2-91"></a> <span class="at">Model.Description.Short =</span> Model.Description.Short,</span>
<span id="cb2-92"><a href="#cb2-92"></a> <span class="at">Model.Parameters =</span> pa,</span>
<span id="cb2-93"><a href="#cb2-93"></a> <span class="at">Empirical.Parameters =</span> Empirical.Parameters)</span>
<span id="cb2-94"><a href="#cb2-94"></a></span>
<span id="cb2-95"><a href="#cb2-95"></a> <span class="fu">return</span>(output)</span>
<span id="cb2-96"><a href="#cb2-96"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
</section>
</div>
<div id="tabset-1-2" class="tab-pane" role="tabpanel" aria-labelledby="tabset-1-2-tab">
<p>The canonical parameters of the <span class="math inline">\(\chi^2\)</span>AR(1) model (which we referred to as <code>ChiAR(1)</code> in the simulation) are the autoregressive parameter <span class="math inline">\(\phi\)</span> (default: <code>phi = 0.4</code>), and degrees of freedom <span class="math inline">\(\nu\)</span> (default: <code>nu = 3</code>). We set the intercept to zero (<code>c = 0</code>).<a href="#fn1" class="footnote-ref" id="fnref1" role="doc-noteref"><sup>1</sup></a></p>
<section id="construction-1" class="level5" data-number="1.1.1.0.3">
<h5 data-number="1.1.1.0.3" class="anchored" data-anchor-id="construction-1"><span class="header-section-number">1.1.1.0.3</span> Construction</h5>
<p>Similar to the AR(1) model, we need to sample the first observation of the <span class="math inline">\(\chi^2\)</span>AR(1) model from its marginal distribution. However, since this model does not have a closed-form marginal distribution, as an approximation, we instead sample <code>x[1]</code> from a <span class="math inline">\(\chi^2\)</span> distribution with <span class="math inline">\(\nu\)</span> degrees of freedom:</p>
<p><span class="math display">\[
X_1 \sim \chi^2(\nu)
\]</span></p>
<p>Then, we generate the remainder of the time series using the definition of the <span class="math inline">\(\chi^2\)</span>AR(1) model:</p>
<p><span class="math display">\[
\begin{aligned}
X_{t} &= c + \phi X_{t-1} + a_{t} \\
a_{t} &\sim \chi^2(\nu).
\end{aligned}
\]</span></p>
</section>
<section id="code-1" class="level5" data-number="1.1.1.0.4">
<h5 data-number="1.1.1.0.4" class="anchored" data-anchor-id="code-1"><span class="header-section-number">1.1.1.0.4</span> Code</h5>
<div class="cell" data-hash="Code-documentations_cache/html/unnamed-chunk-5_6f5049e009b4ac4e5e973cf814c18e34">
<details>
<summary>Click to expand the code</summary>
<div class="sourceCode cell-code" id="cb3"><pre class="sourceCode numberSource r number-lines code-with-copy"><code class="sourceCode r"><span id="cb3-1"><a href="#cb3-1"></a>dgm_chiar <span class="ot"><-</span> <span class="cf">function</span>(...){</span>
<span id="cb3-2"><a href="#cb3-2"></a></span>
<span id="cb3-3"><a href="#cb3-3"></a> pa <span class="ot"><-</span> <span class="fu">list</span>(...)</span>
<span id="cb3-4"><a href="#cb3-4"></a></span>
<span id="cb3-5"><a href="#cb3-5"></a> <span class="cf">if</span>(<span class="fu">is.list</span>(pa<span class="sc">$</span>pa)) pa <span class="ot"><-</span> pa<span class="sc">$</span>pa</span>
<span id="cb3-6"><a href="#cb3-6"></a></span>
<span id="cb3-7"><a href="#cb3-7"></a> <span class="do">## setting default seed if not given</span></span>
<span id="cb3-8"><a href="#cb3-8"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>phi)) pa<span class="sc">$</span>phi <span class="ot"><-</span> <span class="fl">0.4</span></span>
<span id="cb3-9"><a href="#cb3-9"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>nu)) pa<span class="sc">$</span>nu <span class="ot"><-</span> <span class="dv">3</span></span>
<span id="cb3-10"><a href="#cb3-10"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>c)) pa<span class="sc">$</span>c <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb3-11"><a href="#cb3-11"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>k)) pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb3-12"><a href="#cb3-12"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>T)) pa<span class="sc">$</span>T <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb3-13"><a href="#cb3-13"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>seed)) pa<span class="sc">$</span>seed <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb3-14"><a href="#cb3-14"></a></span>
<span id="cb3-15"><a href="#cb3-15"></a> <span class="do">## drawing the first sample x_1</span></span>
<span id="cb3-16"><a href="#cb3-16"></a> <span class="fu">set.seed</span>(pa<span class="sc">$</span>seed)</span>
<span id="cb3-17"><a href="#cb3-17"></a> x <span class="ot"><-</span> <span class="fu">rep</span>(<span class="cn">NA</span>, pa<span class="sc">$</span>T)</span>
<span id="cb3-18"><a href="#cb3-18"></a> x[<span class="dv">1</span>] <span class="ot"><-</span> <span class="fu">rchisq</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb3-19"><a href="#cb3-19"></a> <span class="at">df =</span> pa<span class="sc">$</span>Mean)</span>
<span id="cb3-20"><a href="#cb3-20"></a></span>
<span id="cb3-21"><a href="#cb3-21"></a> <span class="do">## making the rest of the time series</span></span>
<span id="cb3-22"><a href="#cb3-22"></a> <span class="cf">for</span> (t <span class="cf">in</span> <span class="dv">2</span><span class="sc">:</span>pa<span class="sc">$</span>T){</span>
<span id="cb3-23"><a href="#cb3-23"></a></span>
<span id="cb3-24"><a href="#cb3-24"></a> x[t] <span class="ot"><-</span> pa<span class="sc">$</span>c <span class="sc">+</span> pa<span class="sc">$</span>phi<span class="sc">*</span>x[t<span class="dv">-1</span>] <span class="sc">+</span> <span class="fu">rchisq</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb3-25"><a href="#cb3-25"></a> <span class="at">df =</span> pa<span class="sc">$</span>nu)</span>
<span id="cb3-26"><a href="#cb3-26"></a> }</span>
<span id="cb3-27"><a href="#cb3-27"></a></span>
<span id="cb3-28"><a href="#cb3-28"></a> <span class="do">## quick output of raw time series without book-keeping variables/parameters</span></span>
<span id="cb3-29"><a href="#cb3-29"></a> <span class="cf">if</span>(<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>only.ts))</span>
<span id="cb3-30"><a href="#cb3-30"></a> <span class="cf">if</span>(pa<span class="sc">$</span>only.ts<span class="sc">==</span><span class="cn">TRUE</span>) <span class="fu">return</span>(x)</span>
<span id="cb3-31"><a href="#cb3-31"></a></span>
<span id="cb3-32"><a href="#cb3-32"></a> Empirical.Parameters <span class="ot">=</span> <span class="fu">list</span>(<span class="at">Mean =</span> <span class="fu">mean</span>(x),</span>
<span id="cb3-33"><a href="#cb3-33"></a> <span class="at">Variance =</span> <span class="fu">var</span>(x),</span>
<span id="cb3-34"><a href="#cb3-34"></a> <span class="at">Skewness =</span> moments<span class="sc">::</span><span class="fu">skewness</span>(x),</span>
<span id="cb3-35"><a href="#cb3-35"></a> <span class="at">AR =</span> <span class="fu">acf</span>(x, <span class="at">lag.max =</span> <span class="dv">1</span>, <span class="at">plot =</span> <span class="cn">FALSE</span>)<span class="sc">$</span>acf[<span class="dv">2</span>]</span>
<span id="cb3-36"><a href="#cb3-36"></a> )</span>
<span id="cb3-37"><a href="#cb3-37"></a></span>
<span id="cb3-38"><a href="#cb3-38"></a> <span class="do">## making a LaTeX-ready list description of the model</span></span>
<span id="cb3-39"><a href="#cb3-39"></a> Model.Description <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"$</span><span class="sc">\\</span><span class="st">chi^2AR(1):"</span>, <span class="co"># \\; with",</span></span>
<span id="cb3-40"><a href="#cb3-40"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb3-41"><a href="#cb3-41"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb3-42"><a href="#cb3-42"></a> <span class="st">"("</span>,</span>
<span id="cb3-43"><a href="#cb3-43"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb3-44"><a href="#cb3-44"></a> <span class="st">")"</span>,</span>
<span id="cb3-45"><a href="#cb3-45"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb3-46"><a href="#cb3-46"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb3-47"><a href="#cb3-47"></a> <span class="st">"("</span>,</span>
<span id="cb3-48"><a href="#cb3-48"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb3-49"><a href="#cb3-49"></a> <span class="st">")"</span>,</span>
<span id="cb3-50"><a href="#cb3-50"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">phi = "</span>,</span>
<span id="cb3-51"><a href="#cb3-51"></a> <span class="fu">round</span>(pa<span class="sc">$</span>phi,<span class="dv">2</span>),</span>
<span id="cb3-52"><a href="#cb3-52"></a> <span class="st">"("</span>,</span>
<span id="cb3-53"><a href="#cb3-53"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>AR,<span class="dv">2</span>),</span>
<span id="cb3-54"><a href="#cb3-54"></a> <span class="st">")"</span>,</span>
<span id="cb3-55"><a href="#cb3-55"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">nu = "</span>,</span>
<span id="cb3-56"><a href="#cb3-56"></a> <span class="fu">round</span>(pa<span class="sc">$</span>nu,<span class="dv">2</span>),</span>
<span id="cb3-57"><a href="#cb3-57"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; c = "</span>,</span>
<span id="cb3-58"><a href="#cb3-58"></a> <span class="fu">round</span>(pa<span class="sc">$</span>c,<span class="dv">3</span>),</span>
<span id="cb3-59"><a href="#cb3-59"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; T = "</span>,</span>
<span id="cb3-60"><a href="#cb3-60"></a> pa<span class="sc">$</span>T,</span>
<span id="cb3-61"><a href="#cb3-61"></a> <span class="st">"$"</span>)</span>
<span id="cb3-62"><a href="#cb3-62"></a></span>
<span id="cb3-63"><a href="#cb3-63"></a> Model.Description.Short <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"$c = "</span>,</span>
<span id="cb3-64"><a href="#cb3-64"></a> <span class="fu">round</span>(pa<span class="sc">$</span>c,<span class="dv">2</span>),</span>
<span id="cb3-65"><a href="#cb3-65"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">nu = "</span>,</span>
<span id="cb3-66"><a href="#cb3-66"></a> <span class="fu">round</span>(pa<span class="sc">$</span>nu, <span class="dv">2</span>),</span>
<span id="cb3-67"><a href="#cb3-67"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">phi = "</span>,</span>
<span id="cb3-68"><a href="#cb3-68"></a> <span class="fu">round</span>(pa<span class="sc">$</span>phi, <span class="dv">2</span>),</span>
<span id="cb3-69"><a href="#cb3-69"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">rightarrow"</span>,</span>
<span id="cb3-70"><a href="#cb3-70"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb3-71"><a href="#cb3-71"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Mean, <span class="dv">2</span>),</span>
<span id="cb3-72"><a href="#cb3-72"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">sigma^2 = "</span>,</span>
<span id="cb3-73"><a href="#cb3-73"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Variance,<span class="dv">2</span>),</span>
<span id="cb3-74"><a href="#cb3-74"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb3-75"><a href="#cb3-75"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb3-76"><a href="#cb3-76"></a> <span class="st">"$"</span>)</span>
<span id="cb3-77"><a href="#cb3-77"></a> <span class="do">## making the output object</span></span>
<span id="cb3-78"><a href="#cb3-78"></a> output <span class="ot"><-</span> <span class="fu">list</span>(<span class="at">x =</span> x,</span>
<span id="cb3-79"><a href="#cb3-79"></a> <span class="at">Model.Description =</span> Model.Description,</span>
<span id="cb3-80"><a href="#cb3-80"></a> <span class="at">Model.Description.Short =</span> Model.Description.Short,</span>
<span id="cb3-81"><a href="#cb3-81"></a> <span class="at">Model.Parameters =</span> pa,</span>
<span id="cb3-82"><a href="#cb3-82"></a> <span class="at">Empirical.Parameters =</span> Empirical.Parameters)</span>
<span id="cb3-83"><a href="#cb3-83"></a></span>
<span id="cb3-84"><a href="#cb3-84"></a> <span class="fu">return</span>(output)</span>
<span id="cb3-85"><a href="#cb3-85"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
</section>
</div>
<div id="tabset-1-3" class="tab-pane" role="tabpanel" aria-labelledby="tabset-1-3-tab">
<p>The canonical parameters of the BinAR(1) model (which we referred to as <code>BinAR(1)</code> in the simulation) are the survival probability <span class="math inline">\(\alpha\)</span> (default: <code>alpha = 0.5</code>) and the revival probability <span class="math inline">\(\beta\)</span> (default: <code>beta = 0.4</code>). By default, the maximum value on scale <span class="math inline">\(k\)</span> was set to <code>k = 10</code>.</p>
<section id="construction-2" class="level5" data-number="1.1.1.0.5">
<h5 data-number="1.1.1.0.5" class="anchored" data-anchor-id="construction-2"><span class="header-section-number">1.1.1.0.5</span> Construction</h5>
<p>We first calculate the <span class="math inline">\(\theta\)</span> parameter, which characterizes the marginal distribution of the BinaR(1) model:</p>
<p><span class="math display">\[
\theta = \frac{k \beta}{1-(\alpha-\beta)}
\]</span> Then we draw <span class="math inline">\(X_1\)</span> (<code>x[1]</code>) from the marginal distribution of the model:</p>
<p><span class="math display">\[
X_1 \sim Binom(k, \theta)
\]</span></p>
<p>The rest of time series is generated sequentially, for each time point <span class="math inline">\(t\)</span>, by drawing values for the number of survived (<code>S_t[t]</code>) and revived (<code>R_t[t]</code>) elements of the BinAR(1) model based on the previous observations (<span class="math inline">\(X_{t-1}\)</span>), and then adding them:</p>
<p><span class="math display">\[
\begin{aligned}
S_{t} &\sim Binom(X_{t-1}, \alpha) \\
R_t &\sim Binom(k -X_{t-1}, \beta) \\
X_{t} &= S_t + R_t
\end{aligned}
\]</span></p>
</section>
<section id="code-2" class="level5" data-number="1.1.1.0.6">
<h5 data-number="1.1.1.0.6" class="anchored" data-anchor-id="code-2"><span class="header-section-number">1.1.1.0.6</span> Code</h5>
<div class="cell" data-hash="Code-documentations_cache/html/unnamed-chunk-6_ca071e8cb2f50dcca4ea9906233fcc03">
<details>
<summary>Click to expand the code</summary>
<div class="sourceCode cell-code" id="cb4"><pre class="sourceCode numberSource r number-lines code-with-copy"><code class="sourceCode r"><span id="cb4-1"><a href="#cb4-1"></a>dgm_binar <span class="ot"><-</span> <span class="cf">function</span>(...){</span>
<span id="cb4-2"><a href="#cb4-2"></a></span>
<span id="cb4-3"><a href="#cb4-3"></a> pa <span class="ot"><-</span> <span class="fu">list</span>(...)</span>
<span id="cb4-4"><a href="#cb4-4"></a></span>
<span id="cb4-5"><a href="#cb4-5"></a> <span class="cf">if</span>(<span class="fu">is.list</span>(pa<span class="sc">$</span>pa)) pa <span class="ot"><-</span> pa<span class="sc">$</span>pa</span>
<span id="cb4-6"><a href="#cb4-6"></a></span>
<span id="cb4-7"><a href="#cb4-7"></a> <span class="do">## setting default seed if not given</span></span>
<span id="cb4-8"><a href="#cb4-8"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>alpha)) pa<span class="sc">$</span>alpha <span class="ot"><-</span> <span class="fl">0.5</span></span>
<span id="cb4-9"><a href="#cb4-9"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>beta)) pa<span class="sc">$</span>beta <span class="ot"><-</span> <span class="fl">0.4</span></span>
<span id="cb4-10"><a href="#cb4-10"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>k)) pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">10</span></span>
<span id="cb4-11"><a href="#cb4-11"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>T)) pa<span class="sc">$</span>T <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb4-12"><a href="#cb4-12"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>seed)) pa<span class="sc">$</span>seed <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb4-13"><a href="#cb4-13"></a></span>
<span id="cb4-14"><a href="#cb4-14"></a> <span class="do">## making other parameters</span></span>
<span id="cb4-15"><a href="#cb4-15"></a> pa<span class="sc">$</span>rho <span class="ot"><-</span> pa<span class="sc">$</span>alpha <span class="sc">-</span> pa<span class="sc">$</span>beta</span>
<span id="cb4-16"><a href="#cb4-16"></a> pa<span class="sc">$</span>theta <span class="ot"><-</span> pa<span class="sc">$</span>beta<span class="sc">/</span>(<span class="dv">1</span><span class="sc">-</span>pa<span class="sc">$</span>rho)</span>
<span id="cb4-17"><a href="#cb4-17"></a></span>
<span id="cb4-18"><a href="#cb4-18"></a></span>
<span id="cb4-19"><a href="#cb4-19"></a> <span class="do">## drawing the first sample x_1</span></span>
<span id="cb4-20"><a href="#cb4-20"></a> <span class="fu">set.seed</span>(pa<span class="sc">$</span>seed)</span>
<span id="cb4-21"><a href="#cb4-21"></a> x <span class="ot"><-</span> <span class="fu">rep</span>(<span class="cn">NA</span>, pa<span class="sc">$</span>T)</span>
<span id="cb4-22"><a href="#cb4-22"></a> x[<span class="dv">1</span>] <span class="ot"><-</span> <span class="fu">rbinom</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb4-23"><a href="#cb4-23"></a> <span class="at">size =</span> pa<span class="sc">$</span>k ,</span>
<span id="cb4-24"><a href="#cb4-24"></a> <span class="at">prob =</span> pa<span class="sc">$</span>theta)</span>
<span id="cb4-25"><a href="#cb4-25"></a></span>
<span id="cb4-26"><a href="#cb4-26"></a> <span class="do">## making the rest of the time series</span></span>
<span id="cb4-27"><a href="#cb4-27"></a> <span class="cf">for</span> (t <span class="cf">in</span> <span class="dv">2</span><span class="sc">:</span>pa<span class="sc">$</span>T){</span>
<span id="cb4-28"><a href="#cb4-28"></a> S_t <span class="ot"><-</span> <span class="fu">rbinom</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb4-29"><a href="#cb4-29"></a> <span class="at">size =</span> x[t<span class="dv">-1</span>],</span>
<span id="cb4-30"><a href="#cb4-30"></a> <span class="at">prob =</span> pa<span class="sc">$</span>alpha)</span>
<span id="cb4-31"><a href="#cb4-31"></a> R_t <span class="ot"><-</span> <span class="fu">rbinom</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb4-32"><a href="#cb4-32"></a> <span class="at">size =</span> pa<span class="sc">$</span>k <span class="sc">-</span> x[t<span class="dv">-1</span>],</span>
<span id="cb4-33"><a href="#cb4-33"></a> <span class="at">prob =</span> pa<span class="sc">$</span>beta)</span>
<span id="cb4-34"><a href="#cb4-34"></a> x[t] <span class="ot"><-</span> S_t <span class="sc">+</span> R_t</span>
<span id="cb4-35"><a href="#cb4-35"></a> }</span>
<span id="cb4-36"><a href="#cb4-36"></a></span>
<span id="cb4-37"><a href="#cb4-37"></a> <span class="do">## quick output of raw time series without book-keeping variables/parameters</span></span>
<span id="cb4-38"><a href="#cb4-38"></a> <span class="cf">if</span>(<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>only.ts))</span>
<span id="cb4-39"><a href="#cb4-39"></a> <span class="cf">if</span>(pa<span class="sc">$</span>only.ts<span class="sc">==</span><span class="cn">TRUE</span>) <span class="fu">return</span>(x)</span>
<span id="cb4-40"><a href="#cb4-40"></a></span>
<span id="cb4-41"><a href="#cb4-41"></a> Empirical.Parameters <span class="ot">=</span> <span class="fu">list</span>(<span class="at">Mean =</span> <span class="fu">mean</span>(x),</span>
<span id="cb4-42"><a href="#cb4-42"></a> <span class="at">Variance =</span> <span class="fu">var</span>(x),</span>
<span id="cb4-43"><a href="#cb4-43"></a> <span class="at">Skewness =</span> moments<span class="sc">::</span><span class="fu">skewness</span>(x),</span>
<span id="cb4-44"><a href="#cb4-44"></a> <span class="at">AR =</span> <span class="fu">acf</span>(x, <span class="at">lag.max =</span> <span class="dv">1</span>, <span class="at">plot =</span> <span class="cn">FALSE</span>)<span class="sc">$</span>acf[<span class="dv">2</span>]</span>
<span id="cb4-45"><a href="#cb4-45"></a> )</span>
<span id="cb4-46"><a href="#cb4-46"></a></span>
<span id="cb4-47"><a href="#cb4-47"></a> <span class="do">## making a LaTeX-ready list description of the model</span></span>
<span id="cb4-48"><a href="#cb4-48"></a> Model.Description <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"$BinAR(1):"</span>,</span>
<span id="cb4-49"><a href="#cb4-49"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb4-50"><a href="#cb4-50"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb4-51"><a href="#cb4-51"></a> <span class="st">"("</span>,</span>
<span id="cb4-52"><a href="#cb4-52"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb4-53"><a href="#cb4-53"></a> <span class="st">")"</span>,</span>
<span id="cb4-54"><a href="#cb4-54"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb4-55"><a href="#cb4-55"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb4-56"><a href="#cb4-56"></a> <span class="st">"("</span>,</span>
<span id="cb4-57"><a href="#cb4-57"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb4-58"><a href="#cb4-58"></a> <span class="st">")"</span>,</span>
<span id="cb4-59"><a href="#cb4-59"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">rho = "</span>,</span>
<span id="cb4-60"><a href="#cb4-60"></a> <span class="fu">round</span>(pa<span class="sc">$</span>rho,<span class="dv">3</span>),</span>
<span id="cb4-61"><a href="#cb4-61"></a> <span class="st">"("</span>,</span>
<span id="cb4-62"><a href="#cb4-62"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>AR,<span class="dv">2</span>),</span>
<span id="cb4-63"><a href="#cb4-63"></a> <span class="st">")"</span>,</span>
<span id="cb4-64"><a href="#cb4-64"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">alpha = "</span>,</span>
<span id="cb4-65"><a href="#cb4-65"></a> <span class="fu">round</span>(pa<span class="sc">$</span>alpha,<span class="dv">2</span>),</span>
<span id="cb4-66"><a href="#cb4-66"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">beta = "</span>,</span>
<span id="cb4-67"><a href="#cb4-67"></a> <span class="fu">round</span>(pa<span class="sc">$</span>beta,<span class="dv">2</span>),</span>
<span id="cb4-68"><a href="#cb4-68"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">theta = "</span>,</span>
<span id="cb4-69"><a href="#cb4-69"></a> <span class="fu">round</span>(pa<span class="sc">$</span>theta,<span class="dv">3</span>),</span>
<span id="cb4-70"><a href="#cb4-70"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; T = "</span>,</span>
<span id="cb4-71"><a href="#cb4-71"></a> pa<span class="sc">$</span>T,</span>
<span id="cb4-72"><a href="#cb4-72"></a> <span class="st">"$"</span>)</span>
<span id="cb4-73"><a href="#cb4-73"></a></span>
<span id="cb4-74"><a href="#cb4-74"></a> Model.Description.Short <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"$k ="</span>,</span>
<span id="cb4-75"><a href="#cb4-75"></a> <span class="fu">round</span>(pa<span class="sc">$</span>k,<span class="dv">2</span>),</span>
<span id="cb4-76"><a href="#cb4-76"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">alpha = "</span>,</span>
<span id="cb4-77"><a href="#cb4-77"></a> <span class="fu">round</span>(pa<span class="sc">$</span>alpha,<span class="dv">2</span>),</span>
<span id="cb4-78"><a href="#cb4-78"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">beta = "</span>,</span>
<span id="cb4-79"><a href="#cb4-79"></a> <span class="fu">round</span>(pa<span class="sc">$</span>beta,<span class="dv">2</span>),</span>
<span id="cb4-80"><a href="#cb4-80"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">rightarrow"</span>,</span>
<span id="cb4-81"><a href="#cb4-81"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">rho = "</span>,</span>
<span id="cb4-82"><a href="#cb4-82"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>AR, <span class="dv">2</span>),</span>
<span id="cb4-83"><a href="#cb4-83"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb4-84"><a href="#cb4-84"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Mean, <span class="dv">2</span>),</span>
<span id="cb4-85"><a href="#cb4-85"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">sigma^2 = "</span>,</span>
<span id="cb4-86"><a href="#cb4-86"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Variance,<span class="dv">2</span>),</span>
<span id="cb4-87"><a href="#cb4-87"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb4-88"><a href="#cb4-88"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb4-89"><a href="#cb4-89"></a> <span class="st">"$"</span>)</span>
<span id="cb4-90"><a href="#cb4-90"></a> <span class="do">## making the output object</span></span>
<span id="cb4-91"><a href="#cb4-91"></a> output <span class="ot"><-</span> <span class="fu">list</span>(<span class="at">x =</span> x,</span>
<span id="cb4-92"><a href="#cb4-92"></a> <span class="at">Model.Description =</span> Model.Description,</span>
<span id="cb4-93"><a href="#cb4-93"></a> <span class="at">Model.Description.Short =</span> Model.Description.Short,</span>
<span id="cb4-94"><a href="#cb4-94"></a> <span class="at">Model.Parameters =</span> pa,</span>
<span id="cb4-95"><a href="#cb4-95"></a> <span class="at">Empirical.Parameters =</span> Empirical.Parameters)</span>
<span id="cb4-96"><a href="#cb4-96"></a></span>
<span id="cb4-97"><a href="#cb4-97"></a> <span class="fu">return</span>(output)</span>
<span id="cb4-98"><a href="#cb4-98"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
</section>
</div>
<div id="tabset-1-4" class="tab-pane" role="tabpanel" aria-labelledby="tabset-1-4-tab">
<p>The canonical parameters of the PoDAR(1) model (which we referred to as <code>PoDAR(1)</code> in the simulation) are the persistence probability <span class="math inline">\(\tau\)</span> (default: <code>tau = 0.7</code>) and the average rate <span class="math inline">\(\lambda\)</span> (default: <code>lambda = 0.5</code>).</p>
<section id="construction-3" class="level5" data-number="1.1.1.0.7">
<h5 data-number="1.1.1.0.7" class="anchored" data-anchor-id="construction-3"><span class="header-section-number">1.1.1.0.7</span> Construction</h5>
<p>To generate the time series, we first draw the first observation <span class="math inline">\(X_1\)</span> (<code>x[1]</code>) from a Poisson distribution with rate <span class="math inline">\(\lambda\)</span>:</p>
<p><span class="math display">\[
X_1 \sim Poisson(\lambda)
\]</span></p>
<p>And generate the rest of the time series by first drawing <span class="math inline">\(Z_t\)</span> from a Poisson distribution with rate <span class="math inline">\(\lambda\)</span> and <span class="math inline">\(P_t\)</span> from a binomial distribution with size probability of success <span class="math inline">\(\tau\)</span> (that is equivalent to a Bernoulli distribution with probability <span class="math inline">\(\tau\)</span>). Then, we calculate <span class="math inline">\(X_t\)</span> based on the previous observation (<code>x[t-1]</code>) and values of <span class="math inline">\(Z_t\)</span> (<code>Z_t[t]</code>) and <span class="math inline">\(P_t\)</span> (<code>P_t[t]</code>), using the definition of the PoDAR(1) model:</p>
<p><span class="math display">\[
\begin{aligned}
Z_t &\sim Poisson(\lambda) \\
P_t &\sim Binom(1, \tau) \\
X_t &= P_t X_{t-1} + (1-P_t) Z_t
\end{aligned}
\]</span></p>
</section>
<section id="code-3" class="level5" data-number="1.1.1.0.8">
<h5 data-number="1.1.1.0.8" class="anchored" data-anchor-id="code-3"><span class="header-section-number">1.1.1.0.8</span> Code</h5>
<div class="cell" data-hash="Code-documentations_cache/html/unnamed-chunk-7_f5b8e0d8b0a07458adafe9bc78fbd885">
<details>
<summary>Click to expand the code</summary>
<div class="sourceCode cell-code" id="cb5"><pre class="sourceCode numberSource r number-lines code-with-copy"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1"></a>dgm_podar <span class="ot"><-</span> <span class="cf">function</span>(...){</span>
<span id="cb5-2"><a href="#cb5-2"></a></span>
<span id="cb5-3"><a href="#cb5-3"></a> pa <span class="ot"><-</span> <span class="fu">list</span>(...)</span>
<span id="cb5-4"><a href="#cb5-4"></a></span>
<span id="cb5-5"><a href="#cb5-5"></a> <span class="cf">if</span>(<span class="fu">is.list</span>(pa<span class="sc">$</span>pa)) pa <span class="ot"><-</span> pa<span class="sc">$</span>pa</span>
<span id="cb5-6"><a href="#cb5-6"></a></span>
<span id="cb5-7"><a href="#cb5-7"></a> <span class="do">## setting default seed if not given</span></span>
<span id="cb5-8"><a href="#cb5-8"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>tau)) pa<span class="sc">$</span>tau <span class="ot"><-</span> <span class="fl">0.7</span></span>
<span id="cb5-9"><a href="#cb5-9"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>lambda)) pa<span class="sc">$</span>lambda <span class="ot"><-</span> <span class="fl">0.5</span></span>
<span id="cb5-10"><a href="#cb5-10"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>k)) pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">6</span></span>
<span id="cb5-11"><a href="#cb5-11"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>T)) pa<span class="sc">$</span>T <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb5-12"><a href="#cb5-12"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>seed)) pa<span class="sc">$</span>seed <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb5-13"><a href="#cb5-13"></a></span>
<span id="cb5-14"><a href="#cb5-14"></a></span>
<span id="cb5-15"><a href="#cb5-15"></a> <span class="do">## drawing the first sample x_1</span></span>
<span id="cb5-16"><a href="#cb5-16"></a> <span class="fu">set.seed</span>(pa<span class="sc">$</span>seed)</span>
<span id="cb5-17"><a href="#cb5-17"></a> x <span class="ot"><-</span> <span class="fu">rep</span>(<span class="cn">NA</span>, pa<span class="sc">$</span>T)</span>
<span id="cb5-18"><a href="#cb5-18"></a> x[<span class="dv">1</span>] <span class="ot"><-</span> <span class="fu">rpois</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb5-19"><a href="#cb5-19"></a> <span class="at">lambda =</span> pa<span class="sc">$</span>lambda)</span>
<span id="cb5-20"><a href="#cb5-20"></a></span>
<span id="cb5-21"><a href="#cb5-21"></a> <span class="do">## making the rest of the time series</span></span>
<span id="cb5-22"><a href="#cb5-22"></a> <span class="cf">for</span> (t <span class="cf">in</span> <span class="dv">2</span><span class="sc">:</span>pa<span class="sc">$</span>T){</span>
<span id="cb5-23"><a href="#cb5-23"></a> V_t <span class="ot"><-</span> <span class="fu">rbinom</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb5-24"><a href="#cb5-24"></a> <span class="at">size =</span> <span class="dv">1</span>,</span>
<span id="cb5-25"><a href="#cb5-25"></a> <span class="at">prob =</span> pa<span class="sc">$</span>tau)</span>
<span id="cb5-26"><a href="#cb5-26"></a> Z_t <span class="ot"><-</span> <span class="fu">rpois</span>(<span class="at">n =</span> <span class="dv">1</span>,</span>
<span id="cb5-27"><a href="#cb5-27"></a> <span class="at">lambda =</span> pa<span class="sc">$</span>lambda)</span>
<span id="cb5-28"><a href="#cb5-28"></a> x[t] <span class="ot"><-</span> V_t<span class="sc">*</span>x[t<span class="dv">-1</span>] <span class="sc">+</span> (<span class="dv">1</span><span class="sc">-</span>V_t)<span class="sc">*</span>Z_t</span>
<span id="cb5-29"><a href="#cb5-29"></a> }</span>
<span id="cb5-30"><a href="#cb5-30"></a></span>
<span id="cb5-31"><a href="#cb5-31"></a> <span class="do">## quick output of raw time series without book-keeping variables/parameters</span></span>
<span id="cb5-32"><a href="#cb5-32"></a> <span class="cf">if</span>(<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>only.ts))</span>
<span id="cb5-33"><a href="#cb5-33"></a> <span class="cf">if</span>(pa<span class="sc">$</span>only.ts<span class="sc">==</span><span class="cn">TRUE</span>) <span class="fu">return</span>(x)</span>
<span id="cb5-34"><a href="#cb5-34"></a></span>
<span id="cb5-35"><a href="#cb5-35"></a> Empirical.Parameters <span class="ot">=</span> <span class="fu">list</span>(<span class="at">Mean =</span> <span class="fu">mean</span>(x),</span>
<span id="cb5-36"><a href="#cb5-36"></a> <span class="at">Variance =</span> <span class="fu">var</span>(x),</span>
<span id="cb5-37"><a href="#cb5-37"></a> <span class="at">Skewness =</span> moments<span class="sc">::</span><span class="fu">skewness</span>(x),</span>
<span id="cb5-38"><a href="#cb5-38"></a> <span class="at">AR =</span> <span class="fu">acf</span>(x, <span class="at">lag.max =</span> <span class="dv">1</span>, <span class="at">plot =</span> <span class="cn">FALSE</span>)<span class="sc">$</span>acf[<span class="dv">2</span>]</span>
<span id="cb5-39"><a href="#cb5-39"></a> )</span>
<span id="cb5-40"><a href="#cb5-40"></a></span>
<span id="cb5-41"><a href="#cb5-41"></a> <span class="do">## making a LaTeX-ready list description of the model</span></span>
<span id="cb5-42"><a href="#cb5-42"></a> Model.Description <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"$PoDAR(1):"</span>,</span>
<span id="cb5-43"><a href="#cb5-43"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb5-44"><a href="#cb5-44"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb5-45"><a href="#cb5-45"></a> <span class="st">"("</span>,</span>
<span id="cb5-46"><a href="#cb5-46"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Mean,<span class="dv">2</span>),</span>
<span id="cb5-47"><a href="#cb5-47"></a> <span class="st">")"</span>,</span>
<span id="cb5-48"><a href="#cb5-48"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb5-49"><a href="#cb5-49"></a> <span class="fu">round</span>(pa<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb5-50"><a href="#cb5-50"></a> <span class="st">"("</span>,</span>
<span id="cb5-51"><a href="#cb5-51"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb5-52"><a href="#cb5-52"></a> <span class="st">")"</span>,</span>
<span id="cb5-53"><a href="#cb5-53"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">tau = "</span>,</span>
<span id="cb5-54"><a href="#cb5-54"></a> <span class="fu">round</span>(pa<span class="sc">$</span>tau,<span class="dv">3</span>),</span>
<span id="cb5-55"><a href="#cb5-55"></a> <span class="st">"("</span>,</span>
<span id="cb5-56"><a href="#cb5-56"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>AR,<span class="dv">2</span>),</span>
<span id="cb5-57"><a href="#cb5-57"></a> <span class="st">")"</span>,</span>
<span id="cb5-58"><a href="#cb5-58"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">lambda = "</span>,</span>
<span id="cb5-59"><a href="#cb5-59"></a> <span class="fu">round</span>(pa<span class="sc">$</span>lambda,<span class="dv">3</span>),</span>
<span id="cb5-60"><a href="#cb5-60"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; T = "</span>,</span>
<span id="cb5-61"><a href="#cb5-61"></a> pa<span class="sc">$</span>T,</span>
<span id="cb5-62"><a href="#cb5-62"></a> <span class="st">"$"</span>)</span>
<span id="cb5-63"><a href="#cb5-63"></a></span>
<span id="cb5-64"><a href="#cb5-64"></a></span>
<span id="cb5-65"><a href="#cb5-65"></a> Model.Description.Short <span class="ot"><-</span> <span class="fu">paste0</span>(<span class="st">"$</span><span class="sc">\\</span><span class="st">lambda = "</span>,</span>
<span id="cb5-66"><a href="#cb5-66"></a> <span class="fu">round</span>(pa<span class="sc">$</span>lambda,<span class="dv">3</span>),</span>
<span id="cb5-67"><a href="#cb5-67"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">tau = "</span>,</span>
<span id="cb5-68"><a href="#cb5-68"></a> <span class="fu">round</span>(pa<span class="sc">$</span>tau, <span class="dv">2</span>),</span>
<span id="cb5-69"><a href="#cb5-69"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">rightarrow"</span>,</span>
<span id="cb5-70"><a href="#cb5-70"></a> <span class="st">"</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">rho = "</span>,</span>
<span id="cb5-71"><a href="#cb5-71"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>AR, <span class="dv">2</span>),</span>
<span id="cb5-72"><a href="#cb5-72"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">mu = "</span>,</span>
<span id="cb5-73"><a href="#cb5-73"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Mean, <span class="dv">2</span>),</span>
<span id="cb5-74"><a href="#cb5-74"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">sigma^2 = "</span>,</span>
<span id="cb5-75"><a href="#cb5-75"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Variance,<span class="dv">2</span>),</span>
<span id="cb5-76"><a href="#cb5-76"></a> <span class="st">",</span><span class="sc">\\</span><span class="st">; </span><span class="sc">\\</span><span class="st">gamma = "</span>,</span>
<span id="cb5-77"><a href="#cb5-77"></a> <span class="fu">round</span>(Empirical.Parameters<span class="sc">$</span>Skewness,<span class="dv">2</span>),</span>
<span id="cb5-78"><a href="#cb5-78"></a> <span class="st">"$"</span>)</span>
<span id="cb5-79"><a href="#cb5-79"></a> <span class="do">## making the output object</span></span>
<span id="cb5-80"><a href="#cb5-80"></a> output <span class="ot"><-</span> <span class="fu">list</span>(<span class="at">x =</span> x,</span>
<span id="cb5-81"><a href="#cb5-81"></a> <span class="at">Model.Description =</span> Model.Description,</span>
<span id="cb5-82"><a href="#cb5-82"></a> <span class="at">Model.Description.Short =</span> Model.Description.Short,</span>
<span id="cb5-83"><a href="#cb5-83"></a> <span class="at">Model.Parameters =</span> pa,</span>
<span id="cb5-84"><a href="#cb5-84"></a> <span class="at">Empirical.Parameters =</span> Empirical.Parameters)</span>
<span id="cb5-85"><a href="#cb5-85"></a></span>
<span id="cb5-86"><a href="#cb5-86"></a></span>
<span id="cb5-87"><a href="#cb5-87"></a> <span class="fu">return</span>(output)</span>
<span id="cb5-88"><a href="#cb5-88"></a>}</span>
<span id="cb5-89"><a href="#cb5-89"></a></span>
<span id="cb5-90"><a href="#cb5-90"></a><span class="co"># General DGM wrappers ----------------------------------------------------</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
</section>
</div>
</div>
</div>
</section>
<section id="sec-general-dgm-wrappers" class="level3" data-number="1.1.2">
<h3 data-number="1.1.2" class="anchored" data-anchor-id="sec-general-dgm-wrappers"><span class="header-section-number">1.1.2</span> General DGM wrappers</h3>
<p>Given that, in each model, two canonical parameters characterize the dynamic and marginal features of the generated time series, and given that we have analytic formulas that link the canonical parameters to the model-implied <span class="math inline">\(\phi\)</span>, <span class="math inline">\(\mu\)</span>, <span class="math inline">\(\sigma^2\)</span>, and <span class="math inline">\(\gamma\)</span>, we use a function (<code>dgm_parameterizer</code>) to calculate canonical parameters from two given parameters, and make a complete list of parameters (called <code>pa</code>). This list also includes non-parameter variables, importantly, the time series length <span class="math inline">\(T\)</span> (saved in <code>pa$T</code>) and the random seed used in the <code>dgm_*</code> functions (saved in <code>pa$seed</code>). A wrapper function (<code>dgm_generator</code>) is used as an interface to all <code>dgm_*</code> functions, which first makes sure the given parameters are sufficient for data generation, makes a complete parameter list <code>pa</code> with the help of <code>dgm_parameterizer</code>, and passes <code>pa</code> to the respective DGM generating function.</p>
<div class="panel-tabset">
<ul class="nav nav-tabs" role="tablist"><li class="nav-item" role="presentation"><a class="nav-link active" id="tabset-2-1-tab" data-bs-toggle="tab" data-bs-target="#tabset-2-1" role="tab" aria-controls="tabset-2-1" aria-selected="true" href="">Parameter conversions</a></li><li class="nav-item" role="presentation"><a class="nav-link" id="tabset-2-2-tab" data-bs-toggle="tab" data-bs-target="#tabset-2-2" role="tab" aria-controls="tabset-2-2" aria-selected="false" href="">Wrapper around <code>dgm_*</code> functions</a></li></ul>
<div class="tab-content">
<div id="tabset-2-1" class="tab-pane active" role="tabpanel" aria-labelledby="tabset-2-1-tab">
<p>The function <code>dgm_parameterizer</code> calculates canonical/model-implied parameters of a given DGM (specified using the <code>Model</code> argument) based on the parameters given to it as arguments, and saves them in a list of parameters (<code>pa</code>), which s returned by the function. The function makes sure that the set of parameters provided are sufficient to characterize the dynamic parameter of the model (i.e., the autoregression <span class="math inline">\(\phi\)</span>) and at least one of the marginal parameters (importantly, the mean <span class="math inline">\(\mu\)</span>) but giving default values to some parameters.</p>
<div class="cell" data-hash="Code-documentations_cache/html/unnamed-chunk-8_f1167f9a0adebb5a1622c008060fe05c">
<details>
<summary>Click to expand the code</summary>
<div class="sourceCode cell-code" id="cb6"><pre class="sourceCode numberSource r number-lines code-with-copy"><code class="sourceCode r"><span id="cb6-1"><a href="#cb6-1"></a>dgm_parameterizer <span class="ot"><-</span> <span class="cf">function</span>(...){</span>
<span id="cb6-2"><a href="#cb6-2"></a></span>
<span id="cb6-3"><a href="#cb6-3"></a> pa <span class="ot"><-</span> <span class="fu">list</span>(...)</span>
<span id="cb6-4"><a href="#cb6-4"></a></span>
<span id="cb6-5"><a href="#cb6-5"></a> <span class="cf">if</span>(<span class="fu">is.list</span>(pa<span class="sc">$</span>pa)) pa <span class="ot"><-</span> pa<span class="sc">$</span>pa</span>
<span id="cb6-6"><a href="#cb6-6"></a></span>
<span id="cb6-7"><a href="#cb6-7"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>Model)) pa<span class="sc">$</span>Model <span class="ot"><-</span> <span class="st">"ChiAR(1)"</span></span>
<span id="cb6-8"><a href="#cb6-8"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>phi)) pa<span class="sc">$</span>phi <span class="ot"><-</span> <span class="fl">0.2</span></span>
<span id="cb6-9"><a href="#cb6-9"></a></span>
<span id="cb6-10"><a href="#cb6-10"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-11"><a href="#cb6-11"></a> <span class="do">## for NAR(1)</span></span>
<span id="cb6-12"><a href="#cb6-12"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-13"><a href="#cb6-13"></a></span>
<span id="cb6-14"><a href="#cb6-14"></a> <span class="cf">if</span>(<span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"nar(1)"</span> <span class="sc">|</span> <span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"nar"</span>){</span>
<span id="cb6-15"><a href="#cb6-15"></a></span>
<span id="cb6-16"><a href="#cb6-16"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>k)) pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb6-17"><a href="#cb6-17"></a></span>
<span id="cb6-18"><a href="#cb6-18"></a> <span class="do">## Calculating model parameters</span></span>
<span id="cb6-19"><a href="#cb6-19"></a> <span class="do">## if mean is given</span></span>
<span id="cb6-20"><a href="#cb6-20"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>Mean)) {</span>
<span id="cb6-21"><a href="#cb6-21"></a> <span class="co"># then from the mean formula</span></span>
<span id="cb6-22"><a href="#cb6-22"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> pa<span class="sc">$</span>Mean <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-23"><a href="#cb6-23"></a> }</span>
<span id="cb6-24"><a href="#cb6-24"></a> <span class="do">## if mean is not given, we get it from the intercept</span></span>
<span id="cb6-25"><a href="#cb6-25"></a> <span class="cf">else</span>{</span>
<span id="cb6-26"><a href="#cb6-26"></a> <span class="co"># set intercept to zero, if already not defined</span></span>
<span id="cb6-27"><a href="#cb6-27"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>c)) pa<span class="sc">$</span>c <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb6-28"><a href="#cb6-28"></a> <span class="co"># then from the mean formula</span></span>
<span id="cb6-29"><a href="#cb6-29"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> pa<span class="sc">$</span>c <span class="sc">/</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-30"><a href="#cb6-30"></a> }</span>
<span id="cb6-31"><a href="#cb6-31"></a></span>
<span id="cb6-32"><a href="#cb6-32"></a> <span class="do">## Now we certainly have the mean calculated</span></span>
<span id="cb6-33"><a href="#cb6-33"></a></span>
<span id="cb6-34"><a href="#cb6-34"></a> <span class="do">## The within-person (marginal) variance is more interpretable. So</span></span>
<span id="cb6-35"><a href="#cb6-35"></a> <span class="do">## we always calculate the residual variance based on the marginal</span></span>
<span id="cb6-36"><a href="#cb6-36"></a> <span class="do">## variance, unless var.resid is known BUT var.marginal is NOT.</span></span>
<span id="cb6-37"><a href="#cb6-37"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>var.resid) <span class="sc">&</span> <span class="fu">is.null</span>(pa<span class="sc">$</span>var.marginal)) {</span>
<span id="cb6-38"><a href="#cb6-38"></a> <span class="co"># from the marginal variance formula</span></span>
<span id="cb6-39"><a href="#cb6-39"></a> pa<span class="sc">$</span>var.marginal <span class="ot"><-</span> pa<span class="sc">$</span>var.resid <span class="sc">/</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>)</span>
<span id="cb6-40"><a href="#cb6-40"></a> }</span>
<span id="cb6-41"><a href="#cb6-41"></a> <span class="do">## If none of them are known we first give var.marginal a default</span></span>
<span id="cb6-42"><a href="#cb6-42"></a> <span class="do">## value and calculate var.resid based on that.</span></span>
<span id="cb6-43"><a href="#cb6-43"></a> <span class="cf">if</span> (<span class="fu">is.null</span>(pa<span class="sc">$</span>var.resid) <span class="sc">&</span> <span class="fu">is.null</span>(pa<span class="sc">$</span>var.marginal)) {</span>
<span id="cb6-44"><a href="#cb6-44"></a> <span class="co"># we set a default value for the marginal variance</span></span>
<span id="cb6-45"><a href="#cb6-45"></a> pa<span class="sc">$</span>var.marginal <span class="ot"><-</span> <span class="dv">4</span></span>
<span id="cb6-46"><a href="#cb6-46"></a> }</span>
<span id="cb6-47"><a href="#cb6-47"></a> <span class="do">## Then, we calculate var.resid (again) with var.marginal to make</span></span>
<span id="cb6-48"><a href="#cb6-48"></a> <span class="do">## sure that the marginal variance had been given more importance</span></span>
<span id="cb6-49"><a href="#cb6-49"></a> pa<span class="sc">$</span>var.resid <span class="ot"><-</span> pa<span class="sc">$</span>var.marginal <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>)</span>
<span id="cb6-50"><a href="#cb6-50"></a></span>
<span id="cb6-51"><a href="#cb6-51"></a> <span class="do">## Returning the parameter list</span></span>
<span id="cb6-52"><a href="#cb6-52"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-53"><a href="#cb6-53"></a></span>
<span id="cb6-54"><a href="#cb6-54"></a> }</span>
<span id="cb6-55"><a href="#cb6-55"></a></span>
<span id="cb6-56"><a href="#cb6-56"></a></span>
<span id="cb6-57"><a href="#cb6-57"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-58"><a href="#cb6-58"></a> <span class="do">## for ChiAR(1)</span></span>
<span id="cb6-59"><a href="#cb6-59"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-60"><a href="#cb6-60"></a></span>
<span id="cb6-61"><a href="#cb6-61"></a> <span class="cf">if</span>(<span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"chiar(1)"</span> <span class="sc">|</span> <span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"chiar"</span> <span class="sc">|</span></span>
<span id="cb6-62"><a href="#cb6-62"></a> <span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"chi2ar(1)"</span> <span class="sc">|</span> <span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"chi2ar"</span>){</span>
<span id="cb6-63"><a href="#cb6-63"></a></span>
<span id="cb6-64"><a href="#cb6-64"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>k)) pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">100</span></span>
<span id="cb6-65"><a href="#cb6-65"></a></span>
<span id="cb6-66"><a href="#cb6-66"></a> <span class="do">## Calculating model parameters</span></span>
<span id="cb6-67"><a href="#cb6-67"></a> <span class="do">## if mean and skewness are given</span></span>
<span id="cb6-68"><a href="#cb6-68"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>Mean) <span class="sc">&</span> <span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>Skewness)) {</span>
<span id="cb6-69"><a href="#cb6-69"></a> <span class="co"># from the skewness formula</span></span>
<span id="cb6-70"><a href="#cb6-70"></a> pa<span class="sc">$</span>nu <span class="ot"><-</span></span>
<span id="cb6-71"><a href="#cb6-71"></a> <span class="dv">8</span> <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>) <span class="sc">^</span> <span class="dv">3</span> <span class="sc">/</span> ((pa<span class="sc">$</span>Skewness <span class="sc">^</span> <span class="dv">2</span>) <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span></span>
<span id="cb6-72"><a href="#cb6-72"></a> <span class="dv">3</span>) <span class="sc">^</span> <span class="dv">2</span>)</span>
<span id="cb6-73"><a href="#cb6-73"></a> <span class="co"># then from the mean formula</span></span>
<span id="cb6-74"><a href="#cb6-74"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> pa<span class="sc">$</span>Mean <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi) <span class="sc">-</span> pa<span class="sc">$</span>nu</span>
<span id="cb6-75"><a href="#cb6-75"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-76"><a href="#cb6-76"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-77"><a href="#cb6-77"></a> }</span>
<span id="cb6-78"><a href="#cb6-78"></a> <span class="do">## if mean and skewness are not given at the same time, we use c</span></span>
<span id="cb6-79"><a href="#cb6-79"></a> <span class="do">## and one other parameter</span></span>
<span id="cb6-80"><a href="#cb6-80"></a> <span class="cf">else</span>{</span>
<span id="cb6-81"><a href="#cb6-81"></a></span>
<span id="cb6-82"><a href="#cb6-82"></a> <span class="co"># set intercept to zero, if already not defined</span></span>
<span id="cb6-83"><a href="#cb6-83"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>c)) pa<span class="sc">$</span>c <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb6-84"><a href="#cb6-84"></a></span>
<span id="cb6-85"><a href="#cb6-85"></a> <span class="do">## if mean is given</span></span>
<span id="cb6-86"><a href="#cb6-86"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>Mean)) {</span>
<span id="cb6-87"><a href="#cb6-87"></a> <span class="co"># then from the mean formula</span></span>
<span id="cb6-88"><a href="#cb6-88"></a> pa<span class="sc">$</span>nu <span class="ot"><-</span> pa<span class="sc">$</span>Mean <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi) <span class="sc">-</span> pa<span class="sc">$</span>c</span>
<span id="cb6-89"><a href="#cb6-89"></a> <span class="co"># from the skewness formula</span></span>
<span id="cb6-90"><a href="#cb6-90"></a> pa<span class="sc">$</span>Skewness <span class="ot"><-</span></span>
<span id="cb6-91"><a href="#cb6-91"></a> <span class="dv">2</span> <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>) <span class="sc">^</span> <span class="fl">1.5</span> <span class="sc">/</span> (<span class="fu">sqrt</span>(pa<span class="sc">$</span>nu <span class="sc">/</span> <span class="dv">2</span>) <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span></span>
<span id="cb6-92"><a href="#cb6-92"></a> <span class="dv">3</span>))</span>
<span id="cb6-93"><a href="#cb6-93"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-94"><a href="#cb6-94"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-95"><a href="#cb6-95"></a> }</span>
<span id="cb6-96"><a href="#cb6-96"></a></span>
<span id="cb6-97"><a href="#cb6-97"></a> <span class="do">## if skewness is given</span></span>
<span id="cb6-98"><a href="#cb6-98"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>Skewness)) {</span>
<span id="cb6-99"><a href="#cb6-99"></a> <span class="co"># from the skewness formula</span></span>
<span id="cb6-100"><a href="#cb6-100"></a> pa<span class="sc">$</span>nu <span class="ot"><-</span></span>
<span id="cb6-101"><a href="#cb6-101"></a> <span class="dv">8</span> <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>) <span class="sc">^</span> <span class="dv">3</span> <span class="sc">/</span> ((pa<span class="sc">$</span>Skewness <span class="sc">^</span> <span class="dv">2</span>) <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span></span>
<span id="cb6-102"><a href="#cb6-102"></a> <span class="dv">3</span>) <span class="sc">^</span> <span class="dv">2</span>)</span>
<span id="cb6-103"><a href="#cb6-103"></a> <span class="co"># then from the mean formula</span></span>
<span id="cb6-104"><a href="#cb6-104"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> (pa<span class="sc">$</span>c <span class="sc">+</span> pa<span class="sc">$</span>nu) <span class="sc">/</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-105"><a href="#cb6-105"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-106"><a href="#cb6-106"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-107"><a href="#cb6-107"></a> }</span>
<span id="cb6-108"><a href="#cb6-108"></a></span>
<span id="cb6-109"><a href="#cb6-109"></a> <span class="do">## if nu is given</span></span>
<span id="cb6-110"><a href="#cb6-110"></a> <span class="cf">if</span>(<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>nu)){</span>
<span id="cb6-111"><a href="#cb6-111"></a> <span class="co"># from the mean formula</span></span>
<span id="cb6-112"><a href="#cb6-112"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> (pa<span class="sc">$</span>c <span class="sc">+</span> pa<span class="sc">$</span>nu)<span class="sc">/</span>(<span class="dv">1</span><span class="sc">-</span>pa<span class="sc">$</span>phi)</span>
<span id="cb6-113"><a href="#cb6-113"></a> <span class="co"># from the skewness formula</span></span>
<span id="cb6-114"><a href="#cb6-114"></a> pa<span class="sc">$</span>Skewness <span class="ot"><-</span></span>
<span id="cb6-115"><a href="#cb6-115"></a> <span class="dv">2</span> <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span> <span class="dv">2</span>) <span class="sc">^</span> <span class="fl">1.5</span> <span class="sc">/</span> (<span class="fu">sqrt</span>(pa<span class="sc">$</span>nu <span class="sc">/</span> <span class="dv">2</span>) <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi <span class="sc">^</span></span>
<span id="cb6-116"><a href="#cb6-116"></a> <span class="dv">3</span>))</span>
<span id="cb6-117"><a href="#cb6-117"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-118"><a href="#cb6-118"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-119"><a href="#cb6-119"></a> }</span>
<span id="cb6-120"><a href="#cb6-120"></a></span>
<span id="cb6-121"><a href="#cb6-121"></a> }</span>
<span id="cb6-122"><a href="#cb6-122"></a></span>
<span id="cb6-123"><a href="#cb6-123"></a></span>
<span id="cb6-124"><a href="#cb6-124"></a> }</span>
<span id="cb6-125"><a href="#cb6-125"></a></span>
<span id="cb6-126"><a href="#cb6-126"></a></span>
<span id="cb6-127"><a href="#cb6-127"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-128"><a href="#cb6-128"></a> <span class="do">## for BinAR(1)</span></span>
<span id="cb6-129"><a href="#cb6-129"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-130"><a href="#cb6-130"></a></span>
<span id="cb6-131"><a href="#cb6-131"></a> <span class="cf">if</span> (<span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"binar(1)"</span> <span class="sc">|</span> <span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"binar"</span>) {</span>
<span id="cb6-132"><a href="#cb6-132"></a></span>
<span id="cb6-133"><a href="#cb6-133"></a> <span class="cf">if</span> (<span class="fu">is.null</span>(pa<span class="sc">$</span>k))</span>
<span id="cb6-134"><a href="#cb6-134"></a> pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">10</span></span>
<span id="cb6-135"><a href="#cb6-135"></a></span>
<span id="cb6-136"><a href="#cb6-136"></a> <span class="do">## Calculating model parameters</span></span>
<span id="cb6-137"><a href="#cb6-137"></a></span>
<span id="cb6-138"><a href="#cb6-138"></a> <span class="do">## if mean is given</span></span>
<span id="cb6-139"><a href="#cb6-139"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>Mean)) {</span>
<span id="cb6-140"><a href="#cb6-140"></a> <span class="co"># from skewness formula: m = k*theta</span></span>
<span id="cb6-141"><a href="#cb6-141"></a> pa<span class="sc">$</span>theta <span class="ot"><-</span> pa<span class="sc">$</span>Mean <span class="sc">/</span> pa<span class="sc">$</span>k</span>
<span id="cb6-142"><a href="#cb6-142"></a> <span class="do">## we then calculate skewness based on theta</span></span>
<span id="cb6-143"><a href="#cb6-143"></a> pa<span class="sc">$</span>Skewness <span class="ot"><-</span> (<span class="dv">1</span> <span class="sc">-</span> <span class="dv">2</span> <span class="sc">*</span> pa<span class="sc">$</span>theta) <span class="sc">/</span> <span class="fu">sqrt</span>(pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>theta))</span>
<span id="cb6-144"><a href="#cb6-144"></a> <span class="do">## we then calculate beta based on theta and phi</span></span>
<span id="cb6-145"><a href="#cb6-145"></a> pa<span class="sc">$</span>beta <span class="ot"><-</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-146"><a href="#cb6-146"></a> <span class="co"># then we calculate alpha</span></span>
<span id="cb6-147"><a href="#cb6-147"></a> pa<span class="sc">$</span>alpha <span class="ot"><-</span> pa<span class="sc">$</span>phi <span class="sc">+</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-148"><a href="#cb6-148"></a> <span class="co"># finally we calculate k*beta, which is equivalent to c</span></span>
<span id="cb6-149"><a href="#cb6-149"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-150"><a href="#cb6-150"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-151"><a href="#cb6-151"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-152"><a href="#cb6-152"></a> }</span>
<span id="cb6-153"><a href="#cb6-153"></a></span>
<span id="cb6-154"><a href="#cb6-154"></a> <span class="do">## if skewness is given</span></span>
<span id="cb6-155"><a href="#cb6-155"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>Skewness)) {</span>
<span id="cb6-156"><a href="#cb6-156"></a> <span class="co"># from skewness formula: skewness = (1-2*theta)/sqrt(k*theta*(1-theta))</span></span>
<span id="cb6-157"><a href="#cb6-157"></a> <span class="co"># it is easier to write with ks = k*(skewness^2)</span></span>
<span id="cb6-158"><a href="#cb6-158"></a> ks <span class="ot"><-</span> pa<span class="sc">$</span>k<span class="sc">*</span>pa<span class="sc">$</span>Skewness<span class="sc">^</span><span class="dv">2</span></span>
<span id="cb6-159"><a href="#cb6-159"></a> pa<span class="sc">$</span>theta <span class="ot"><-</span> (ks <span class="sc">+</span> <span class="dv">4</span> <span class="sc">-</span></span>
<span id="cb6-160"><a href="#cb6-160"></a> <span class="fu">sqrt</span>(ks <span class="sc">*</span> (ks <span class="sc">+</span> <span class="dv">4</span>))) <span class="sc">/</span> (<span class="dv">2</span> <span class="sc">*</span> ks <span class="sc">+</span> <span class="dv">8</span>)</span>
<span id="cb6-161"><a href="#cb6-161"></a> <span class="co"># The above formula only gives theta < 0.5, thus for negative skewness</span></span>
<span id="cb6-162"><a href="#cb6-162"></a> <span class="co"># we must use 1-theta instead</span></span>
<span id="cb6-163"><a href="#cb6-163"></a> <span class="cf">if</span>(pa<span class="sc">$</span>Skewness <span class="sc"><</span> <span class="dv">0</span>) pa<span class="sc">$</span>theta <span class="ot"><-</span> <span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>theta</span>
<span id="cb6-164"><a href="#cb6-164"></a> <span class="co"># we then calculate mean based on theta</span></span>
<span id="cb6-165"><a href="#cb6-165"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta</span>
<span id="cb6-166"><a href="#cb6-166"></a> <span class="co"># we then calculate beta based on theta and phi</span></span>
<span id="cb6-167"><a href="#cb6-167"></a> pa<span class="sc">$</span>beta <span class="ot"><-</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-168"><a href="#cb6-168"></a> <span class="co"># then we calculate alpha</span></span>
<span id="cb6-169"><a href="#cb6-169"></a> pa<span class="sc">$</span>alpha <span class="ot"><-</span> pa<span class="sc">$</span>phi <span class="sc">+</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-170"><a href="#cb6-170"></a> <span class="co"># finally we calculate k*beta, which is equivalent to c</span></span>
<span id="cb6-171"><a href="#cb6-171"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-172"><a href="#cb6-172"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-173"><a href="#cb6-173"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-174"><a href="#cb6-174"></a> }</span>
<span id="cb6-175"><a href="#cb6-175"></a></span>
<span id="cb6-176"><a href="#cb6-176"></a> <span class="do">## if theta is given</span></span>
<span id="cb6-177"><a href="#cb6-177"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>theta)) {</span>
<span id="cb6-178"><a href="#cb6-178"></a> <span class="co"># from theta formula</span></span>
<span id="cb6-179"><a href="#cb6-179"></a> pa<span class="sc">$</span>beta <span class="ot"><-</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-180"><a href="#cb6-180"></a> <span class="co"># from beta formula</span></span>
<span id="cb6-181"><a href="#cb6-181"></a> pa<span class="sc">$</span>alpha <span class="ot"><-</span> pa<span class="sc">$</span>phi <span class="sc">+</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-182"><a href="#cb6-182"></a> <span class="co"># we then calculate mean based on theta</span></span>
<span id="cb6-183"><a href="#cb6-183"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta</span>
<span id="cb6-184"><a href="#cb6-184"></a> <span class="co"># we then calculate skewness based on theta</span></span>
<span id="cb6-185"><a href="#cb6-185"></a> pa<span class="sc">$</span>Skewness <span class="ot"><-</span> (<span class="dv">1</span> <span class="sc">-</span> <span class="dv">2</span> <span class="sc">*</span> pa<span class="sc">$</span>theta) <span class="sc">/</span> <span class="fu">sqrt</span>(pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>theta))</span>
<span id="cb6-186"><a href="#cb6-186"></a> <span class="co"># finally we calculate k*beta, which is equivalent to c</span></span>
<span id="cb6-187"><a href="#cb6-187"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-188"><a href="#cb6-188"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-189"><a href="#cb6-189"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-190"><a href="#cb6-190"></a> }</span>
<span id="cb6-191"><a href="#cb6-191"></a></span>
<span id="cb6-192"><a href="#cb6-192"></a> <span class="do">## if alpha is given</span></span>
<span id="cb6-193"><a href="#cb6-193"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>alpha)) {</span>
<span id="cb6-194"><a href="#cb6-194"></a> <span class="co"># from beta formula</span></span>
<span id="cb6-195"><a href="#cb6-195"></a> pa<span class="sc">$</span>beta <span class="ot"><-</span> pa<span class="sc">$</span>alpha <span class="sc">-</span> pa<span class="sc">$</span>phi</span>
<span id="cb6-196"><a href="#cb6-196"></a> <span class="co"># from theta formula</span></span>
<span id="cb6-197"><a href="#cb6-197"></a> pa<span class="sc">$</span>theta <span class="ot"><-</span> pa<span class="sc">$</span>beta <span class="sc">/</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-198"><a href="#cb6-198"></a> <span class="co"># we then calculate mean based on theta</span></span>
<span id="cb6-199"><a href="#cb6-199"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta</span>
<span id="cb6-200"><a href="#cb6-200"></a> <span class="co"># we then calculate skewness based on theta</span></span>
<span id="cb6-201"><a href="#cb6-201"></a> pa<span class="sc">$</span>Skewness <span class="ot"><-</span> (<span class="dv">1</span> <span class="sc">-</span> <span class="dv">2</span> <span class="sc">*</span> pa<span class="sc">$</span>theta) <span class="sc">/</span> <span class="fu">sqrt</span>(pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>theta))</span>
<span id="cb6-202"><a href="#cb6-202"></a> <span class="co"># finally we calculate k*beta, which is equivalent to c</span></span>
<span id="cb6-203"><a href="#cb6-203"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-204"><a href="#cb6-204"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-205"><a href="#cb6-205"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-206"><a href="#cb6-206"></a> }</span>
<span id="cb6-207"><a href="#cb6-207"></a></span>
<span id="cb6-208"><a href="#cb6-208"></a> <span class="do">## if beta is given</span></span>
<span id="cb6-209"><a href="#cb6-209"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>beta)) {</span>
<span id="cb6-210"><a href="#cb6-210"></a> <span class="co"># from beta formula</span></span>
<span id="cb6-211"><a href="#cb6-211"></a> pa<span class="sc">$</span>alpha <span class="ot"><-</span> pa<span class="sc">$</span>phi <span class="sc">+</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-212"><a href="#cb6-212"></a> <span class="co"># from theta formula</span></span>
<span id="cb6-213"><a href="#cb6-213"></a> pa<span class="sc">$</span>theta <span class="ot"><-</span> pa<span class="sc">$</span>beta <span class="sc">/</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-214"><a href="#cb6-214"></a> <span class="co"># we then calculate mean based on theta</span></span>
<span id="cb6-215"><a href="#cb6-215"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta</span>
<span id="cb6-216"><a href="#cb6-216"></a> <span class="co"># we then calculate skewness based on theta</span></span>
<span id="cb6-217"><a href="#cb6-217"></a> pa<span class="sc">$</span>Skewness <span class="ot"><-</span> (<span class="dv">1</span> <span class="sc">-</span> <span class="dv">2</span> <span class="sc">*</span> pa<span class="sc">$</span>theta) <span class="sc">/</span> <span class="fu">sqrt</span>(pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>theta))</span>
<span id="cb6-218"><a href="#cb6-218"></a> <span class="co"># finally we calculate k*beta, which is equivalent to c</span></span>
<span id="cb6-219"><a href="#cb6-219"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-220"><a href="#cb6-220"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-221"><a href="#cb6-221"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-222"><a href="#cb6-222"></a> }</span>
<span id="cb6-223"><a href="#cb6-223"></a></span>
<span id="cb6-224"><a href="#cb6-224"></a> <span class="do">## if intercept given</span></span>
<span id="cb6-225"><a href="#cb6-225"></a> <span class="cf">if</span> (<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>c)) {</span>
<span id="cb6-226"><a href="#cb6-226"></a> <span class="co"># from intercept formula c = k*beta</span></span>
<span id="cb6-227"><a href="#cb6-227"></a> pa<span class="sc">$</span>beta <span class="ot"><-</span> pa<span class="sc">$</span>c <span class="sc">/</span> pa<span class="sc">$</span>k</span>
<span id="cb6-228"><a href="#cb6-228"></a> <span class="co"># from beta formula</span></span>
<span id="cb6-229"><a href="#cb6-229"></a> pa<span class="sc">$</span>alpha <span class="ot"><-</span> pa<span class="sc">$</span>phi <span class="sc">+</span> pa<span class="sc">$</span>beta</span>
<span id="cb6-230"><a href="#cb6-230"></a> <span class="co"># from theta formula</span></span>
<span id="cb6-231"><a href="#cb6-231"></a> pa<span class="sc">$</span>theta <span class="ot"><-</span> pa<span class="sc">$</span>beta <span class="sc">/</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>phi)</span>
<span id="cb6-232"><a href="#cb6-232"></a> <span class="co"># we then calculate mean based on theta</span></span>
<span id="cb6-233"><a href="#cb6-233"></a> pa<span class="sc">$</span>Mean <span class="ot"><-</span> pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta</span>
<span id="cb6-234"><a href="#cb6-234"></a> <span class="co"># we then calculate skewness based on theta</span></span>
<span id="cb6-235"><a href="#cb6-235"></a> pa<span class="sc">$</span>Skewness <span class="ot"><-</span> (<span class="dv">1</span> <span class="sc">-</span> <span class="dv">2</span> <span class="sc">*</span> pa<span class="sc">$</span>theta) <span class="sc">/</span> <span class="fu">sqrt</span>(pa<span class="sc">$</span>k <span class="sc">*</span> pa<span class="sc">$</span>theta <span class="sc">*</span> (<span class="dv">1</span> <span class="sc">-</span> pa<span class="sc">$</span>theta))</span>
<span id="cb6-236"><a href="#cb6-236"></a> <span class="co"># returning the parameter list</span></span>
<span id="cb6-237"><a href="#cb6-237"></a> <span class="fu">return</span>(pa)</span>
<span id="cb6-238"><a href="#cb6-238"></a> }</span>
<span id="cb6-239"><a href="#cb6-239"></a></span>
<span id="cb6-240"><a href="#cb6-240"></a> }</span>
<span id="cb6-241"><a href="#cb6-241"></a></span>
<span id="cb6-242"><a href="#cb6-242"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-243"><a href="#cb6-243"></a> <span class="do">## for PoDAR(1)</span></span>
<span id="cb6-244"><a href="#cb6-244"></a> <span class="do">## %%%%%%%%%%%%</span></span>
<span id="cb6-245"><a href="#cb6-245"></a></span>
<span id="cb6-246"><a href="#cb6-246"></a> <span class="cf">if</span> (<span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"podar(1)"</span> <span class="sc">|</span> <span class="fu">tolower</span>(pa<span class="sc">$</span>Model) <span class="sc">==</span> <span class="st">"podar"</span>) {</span>
<span id="cb6-247"><a href="#cb6-247"></a></span>
<span id="cb6-248"><a href="#cb6-248"></a> <span class="cf">if</span> (<span class="fu">is.null</span>(pa<span class="sc">$</span>k))</span>
<span id="cb6-249"><a href="#cb6-249"></a> pa<span class="sc">$</span>k <span class="ot"><-</span> <span class="dv">50</span></span>
<span id="cb6-250"><a href="#cb6-250"></a></span>
<span id="cb6-251"><a href="#cb6-251"></a> <span class="co"># phi and tau are the same, then if tau is defined, it overrules phi</span></span>
<span id="cb6-252"><a href="#cb6-252"></a> <span class="cf">if</span>(<span class="sc">!</span><span class="fu">is.null</span>(pa<span class="sc">$</span>tau)) pa<span class="sc">$</span>phi <span class="ot"><-</span> pa<span class="sc">$</span>tau</span>
<span id="cb6-253"><a href="#cb6-253"></a> <span class="co"># and if tau is not defined, then tau will get the value of phi</span></span>
<span id="cb6-254"><a href="#cb6-254"></a> <span class="cf">if</span>(<span class="fu">is.null</span>(pa<span class="sc">$</span>tau)) pa<span class="sc">$</span>tau <span class="ot"><-</span> pa<span class="sc">$</span>phi</span>
<span id="cb6-255"><a href="#cb6-255"></a> <span class="co"># DAR(1) intercept is zero</span></span>
<span id="cb6-256"><a href="#cb6-256"></a> pa<span class="sc">$</span>c <span class="ot"><-</span> <span class="dv">0</span></span>
<span id="cb6-257"><a href="#cb6-257"></a></span>