/
RMark-package.R
4764 lines (4658 loc) · 255 KB
/
RMark-package.R
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
#' Multistate Jolly Seber
#' @name MSJollySeber
#' @author Jeff Laake (simulated data from Gary White)
#' @examples
#' \donttest{
# three state example provided by Gary White
#'# import the sample data
#'pathtodata=paste(path.package("RMark"),"extdata",sep="/")
#'msjs3=convert.inp(paste(pathtodata,"MSJS3state.inp",sep="/"))
#'# process with MSJollySeber model and strata labels
#'dp=process.data(msjs3,model="MSJollySeber",strata.labels=c("E","F","G"))
#'# create design data and use time pims for Psi and pent because global model used
#'# to simulate data only had time/strata effects
#'ddl=make.design.data(dp,parameters=list(Psi=list(pim.type="time"),pent=list(pim.type="time")))
#'# p for first occasion is confounded because pent since it is time dependent
#'# p for last occasion is confounded because survival is time dependent
#'# Thus times 1 and 2 are combined as 1 and times 3 and 4 are combined as time 2
#'ddl$p$time[ddl$p$time==2]=1
#'ddl$p$time[ddl$p$time==3]=2
#'ddl$p$time[ddl$p$time==4]=2
#'ddl$p=droplevels(ddl$p)
#'# run global model that was used to generate the data. These results match what was produced in
#'# MARK .dbf that Gary provided to me.
#'model=mark(dp,ddl=ddl,model.parameters=list(S=list(formula=~-1+stratum:time,link="sin"),
#' Psi=list(formula=~-1+stratum:tostratum:time),pent=list(formula=~stratum*time),
#' p=list(formula=~-1+stratum:time,link="sin"),pi=list(formula=~-1+stratum)))
#'# Examine derived abundance estimates
#'model$results$derived
#'#$`N*`
#'# estimate se lcl ucl
#'#1 1200 22.88454 1159.6 1249.699
#'#
#'#$`N_TSA(s*t)`
#'#estimate se lcl ucl
#'#1 180.00027 55.22876 126.98046 381.1247
#'#2 249.59448 73.28850 177.32650 511.2946
#'#3 377.20653 92.15176 276.37043 681.2139
#'#4 422.14994 105.74592 307.97920 774.6687
#'#5 119.99866 38.18683 84.16549 261.3716
#'#6 193.39459 57.72862 137.03323 401.0494
#'#7 296.19950 71.38362 217.52522 530.3700
#'#8 350.09106 86.97180 255.76981 639.0182
#'#9 60.00109 21.23089 41.40889 142.5095
#'#10 120.99888 38.29786 84.93901 262.4367
#'#11 194.19013 45.93703 143.06755 343.7371
#'#12 248.60043 61.37312 181.81724 451.9540
#'msjs3$grp=factor(as.numeric(runif(1686)<0.5))
#'dp=process.data(msjs3,model="MSJollySeber",strata.labels=c("E","F","G"),groups="grp")
#'# create design data and use time pims for Psi and pent because global model used
#'# to simulate data only had time/strata effects
#'ddl=make.design.data(dp,parameters=list(Psi=list(pim.type="time"),pent=list(pim.type="time")))
#'# p for first occasion is confounded because pent since it is time dependent
#'# p for last occasion is confounded because survival is time dependent
#'# Thus times 1 and 2 are combined as 1 and times 3 and 4 are combined as time 2
#'ddl$p$time[ddl$p$time==2]=1
#'ddl$p$time[ddl$p$time==3]=2
#'ddl$p$time[ddl$p$time==4]=2
#'ddl$p=droplevels(ddl$p)
#'model=mark(dp,ddl=ddl,model.parameters=list(S=list(formula=~-1+stratum:time,link="sin"),
#' Psi=list(formula=~-1+stratum:tostratum:time),pent=list(formula=~stratum*time),
#' p=list(formula=~-1+stratum:time,link="sin"),pi=list(formula=~-1+stratum)))
#'# Note that I don't have groups separated out; probably need to add these but first 12 are for group 1 and
#'# second 12 are for group 2
#'model$results$derived
#'#$`N*`
#'#estimate se lcl ucl
#'#1 576.5519 10.99515 557.1415 600.4303
#'#2 623.4481 11.88948 602.4589 649.2687
#'#
#'#$`N_TSA(s*t)`
#'#estimate se lcl ucl
#'#1 86.48235 26.53543 59.74870 179.73562
#'#2 119.91930 35.21259 87.08794 250.88028
#'#3 181.23323 44.27568 158.25938 412.81875
#'#4 202.82678 50.80719 138.39582 351.52314
#'#5 57.65490 18.34742 40.00284 124.36399
#'#6 92.91895 27.73659 69.10044 202.16049
#'#7 142.31136 34.29734 108.00098 263.47936
#'#8 168.20397 41.78690 132.40312 333.04237
#'#9 28.82821 10.20072 22.06997 76.00783
#'#10 58.13528 18.40095 45.08338 139.96467
#'#11 93.30069 22.07105 72.55341 174.90726
#'#12 119.44262 29.48746 89.82760 223.31249
#'#13 93.51673 28.69379 67.33555 201.88406
#'#14 129.67341 38.07675 90.39324 261.12212
#'#15 195.97456 47.87702 131.43369 328.98380
#'#16 219.32450 54.93979 175.00065 444.87035
#'#17 62.34450 19.83979 44.17969 137.09149
#'#18 100.47688 29.99266 68.46238 201.34980
#'#19 153.88682 37.08705 110.11103 269.06238
#'#20 181.88551 45.18580 126.26801 317.46649
#'#21 31.17307 11.03043 19.84077 69.37125
#'#22 62.86395 19.89767 40.98810 128.26377
#'#23 100.88968 23.86628 71.50113 172.47586
#'#24 129.15797 31.88594 92.33709 229.95112
#'#
#'#$`N_TSA(t)`
#'#estimate se lcl ucl
#'#1 172.9655 17.57292 146.1440 216.6946
#'#2 270.9735 17.86549 242.2069 313.2976
#'#3 416.8453 21.54149 382.6336 468.5324
#'#4 490.4734 29.73746 440.1191 557.6749
#'#5 187.0343 19.00228 158.0680 234.3750
#'#6 293.0142 19.31866 260.1141 336.4346
#'#7 450.7511 23.29365 409.1688 500.7981
#'#8 530.3680 32.15627 476.2905 603.5088
#'#
#'#$`N(s*t)`
#'#estimate se lcl ucl
#'#1 81.67410 23.937921 58.04354 167.08272
#'#2 126.83858 37.175211 90.14070 259.47682
#'#3 259.54924 63.478782 190.12917 469.06140
#'#4 148.55131 36.331666 108.81919 268.46423
#'#5 56.08549 16.675093 39.76587 115.96055
#'#6 104.09903 30.950270 73.80855 215.23180
#'#7 157.41243 37.901060 115.62023 281.69787
#'#8 204.11663 49.146288 149.92471 365.27751
#'#9 34.90008 10.950259 24.53295 75.18089
#'#10 71.01092 22.280408 49.91700 152.96997
#'#11 109.13438 25.654572 80.49076 192.43787
#'#12 130.02727 30.565932 95.90006 229.27854
#'#13 98.32499 28.818144 69.87688 201.14587
#'#14 122.75425 35.978133 87.23809 251.12141
#'#15 117.65857 28.776129 86.18914 212.63438
#'#16 273.60001 66.915223 200.42187 494.45417
#'#17 63.91392 19.002608 45.31640 132.14636
#'#18 89.29676 26.549323 63.31340 184.62710
#'#19 138.78578 33.416218 101.93886 248.36450
#'#20 145.97290 35.146702 107.21784 261.22623
#'#21 25.10117 7.875752 17.64482 54.07233
#'#22 49.98828 15.684338 35.13917 107.68352
#'#23 85.05591 19.994369 62.73197 149.98004
#'#24 118.57332 27.873415 87.45234 209.08166
#'#
#'#$`N(t)`
#'#estimate se lcl ucl
#'#1 172.6597 12.84355 151.7126 202.7258
#'#2 301.9485 20.11213 268.5098 348.1893
#'#3 526.0961 35.72379 466.8913 608.4782
#'#4 482.6952 31.12914 430.6676 553.9137
#'#5 187.3401 16.99719 160.7135 228.6425
#'#6 262.0393 19.37969 230.3986 307.3613
#'#7 341.5003 20.10434 307.4281 386.8953
#'#8 538.1462 38.90093 474.3711 628.7710
#'#
#'# For better labels and interpretation see MARK output
#'#Total Abundance Estimates
#'#Grp. N*-hat Standard Error Lower Upper
#'#---- -------------- -------------- -------------- --------------
#'#1 576.55192 10.995146 557.14154 600.43028
#'#2 623.44807 11.889480 602.45888 649.26868
#'#
#'#Abundance Estimates
#'#Grp. Str. Occ. Abundance Standard Error Lower Upper
#'#---- ---- ---- -------------- -------------- -------------- --------------
#'#1 E 1 86.482346 26.535427 59.748702 179.73562
#'#1 E 2 119.91930 35.212594 87.087938 250.88028
#'#1 E 3 181.23323 44.275677 158.25938 412.81875
#'#1 E 4 202.82678 50.807187 138.39582 351.52314
#'#1 F 1 57.654903 18.347425 40.002836 124.36399
#'#1 F 2 92.918950 27.736591 69.100442 202.16049
#'#1 F 3 142.31136 34.297342 108.00098 263.47936
#'#1 F 4 168.20397 41.786899 132.40312 333.04237
#'#1 G 1 28.828208 10.200716 22.069972 76.007829
#'#1 G 2 58.135284 18.400953 45.083379 139.96467
#'#1 G 3 93.300695 22.071048 72.553406 174.90726
#'#1 G 4 119.44262 29.487455 89.827598 223.31249
#'#2 E 1 93.516733 28.693792 67.335545 201.88406
#'#2 E 2 129.67341 38.076751 90.393237 261.12212
#'#2 E 3 195.97456 47.877016 131.43369 328.98380
#'#2 E 4 219.32450 54.939793 175.00065 444.87035
#'#2 F 1 62.344496 19.839786 44.179691 137.09149
#'#2 F 2 100.47688 29.992658 68.462381 201.34980
#'#2 F 3 153.88682 37.087053 110.11103 269.06238
#'#2 F 4 181.88551 45.185803 126.26801 317.46649
#'#2 G 1 31.173066 11.030432 19.840768 69.371247
#'#2 G 2 62.863950 19.897668 40.988101 128.26377
#'#2 G 3 100.88968 23.866285 71.501133 172.47586
#'#2 G 4 129.15797 31.885935 92.337089 229.95112
#'#
#'#Summed Abundance Estimates
#'#Grp. Occ. Abundance Standard Error Lower Upper
#'#---- ---- -------------- -------------- -------------- --------------
#'#1 1 172.96546 17.572916 146.14398 216.69465
#'#1 2 270.97353 17.865494 242.20692 313.29762
#'#1 3 416.84528 21.541486 382.63361 468.53243
#'#1 4 490.47337 29.737455 440.11912 557.67492
#'#2 1 187.03430 19.002280 158.06800 234.37501
#'#2 2 293.01425 19.318655 260.11406 336.43459
#'#2 3 450.75106 23.293649 409.16876 500.79814
#'#2 4 530.36798 32.156270 476.29049 603.50885
#'#
#'#Abundance Estimates without TSA
#'#Grp. Str. Occ. Abundance Standard Error Lower Upper
#'#---- ---- ---- -------------- -------------- -------------- --------------
#'#1 E 1 81.674097 23.937921 58.043545 167.08272
#'#1 E 2 126.83858 37.175211 90.140704 259.47682
#'#1 E 3 259.54924 63.478782 190.12917 469.06140
#'#1 E 4 148.55131 36.331666 108.81919 268.46423
#'#1 F 1 56.085487 16.675093 39.765869 115.96055
#'#1 F 2 104.09903 30.950270 73.808547 215.23180
#'#1 F 3 157.41243 37.901060 115.62023 281.69787
#'#1 F 4 204.11663 49.146288 149.92471 365.27751
#'#1 G 1 34.900079 10.950259 24.532949 75.180890
#'#1 G 2 71.010918 22.280408 49.917000 152.96997
#'#1 G 3 109.13438 25.654572 80.490758 192.43787
#'#1 G 4 130.02727 30.565932 95.900063 229.27854
#'#2 E 1 98.324993 28.818144 69.876881 201.14587
#'#2 E 2 122.75425 35.978133 87.238086 251.12141
#'#2 E 3 117.65857 28.776129 86.189140 212.63438
#'#2 E 4 273.60001 66.915223 200.42187 494.45417
#'#2 F 1 63.913917 19.002608 45.316401 132.14636
#'#2 F 2 89.296755 26.549323 63.313404 184.62710
#'#2 F 3 138.78578 33.416218 101.93886 248.36450
#'#2 F 4 145.97290 35.146702 107.21784 261.22623
#'#2 G 1 25.101172 7.8757515 17.644825 54.072326
#'#2 G 2 49.988280 15.684338 35.139173 107.68352
#'#2 G 3 85.055913 19.994369 62.731974 149.98004
#'#2 G 4 118.57332 27.873415 87.452338 209.08166
#'#
#'#Summed Abundance Estimates without TSA
#'#Grp. Occ. Abundance Standard Error Lower Upper
#'#---- ---- -------------- -------------- -------------- --------------
#'#1 1 172.65966 12.843549 151.71260 202.72578
#'#1 2 301.94852 20.112126 268.50983 348.18933
#'#1 3 526.09605 35.723790 466.89129 608.47820
#'#1 4 482.69521 31.129138 430.66757 553.91371
#'#2 1 187.34008 16.997192 160.71347 228.64245
#'#2 2 262.03929 19.379689 230.39863 307.36128
#'#2 3 341.50026 20.104344 307.42812 386.89528
#'#2 4 538.14624 38.900929 474.37108 628.77101
#'#
#'#Expected Ingress, Egress, and Residence Times 1&2
#'#Grp. Estimate Standard Error Lower Upper
#'#---- -------------- -------------- -------------- --------------
#'#1 2.4000144 0.0650048 2.2726051 2.5274237
#'#1 3.7437021 0.0583736 3.6292898 3.8581145
#'#1 2.3436877 0.0811437 2.1846461 2.5027293
#'#1 3.4389990 0.1604949 3.1244290 3.7535689
#'#2 2.4000144 0.0650048 2.2726051 2.5274237
#'#2 3.7437021 0.0583736 3.6292898 3.8581145
#'#2 2.3436877 0.0811437 2.1846461 2.5027293
#'#2 3.4389990 0.1604949 3.1244290 3.7535689
#'}
NULL
#' Hidden Markov Model
#' @name HidMarkov
#' @author Chris Oosthuizen
#' @examples
#' \donttest{
#'# posted by Chris on phidot 15 Sept 2021
#'# ----------------------------------------------------------------------------------
#'# Fit HMM model with 2 unobservable states
#'# 5 states:
#'# 1 = successful breeder
#'# 2 = failed breeder
#'# 3 = post-successful breeder (unobservable)
#'# 4 = post-failed breeder (unobservable)
#'# 5 = dead
#'
#'# 4 field observations:
#'# 0 = not seen
#'# 1 = seen as successful breeder
#'# 2 = seen as failed breeder
#'# 5 = seen with uncertain breeding state (** this is an event in MARK terms)
#'
#'#-----------------------------------------------------------------------------------
#'
#'# Add a few uncertainty events (event = '5') to the above input data
#'df = c('1010001',
#' '1050000', '1010000', '1010000', '1010100', '1000000', '1000001', '1000001',
#' '1000001', '1010105', '1050101', '1010101', '1050101', '1010101', '1010101',
#' '1010101', '1010101', '1010001', '1000001', '1000001', '1000001', '1010101',
#' '1010101', '1010101', '1050101', '1010101', '1010101', '1010201', '1010201',
#' '1010201', '1010201', '1010201', '1010201', '1010201', '1010201', '1010201',
#' '1010201', '1010201', '1010201', '1010201', '1010201', '1050201', '1010205',
#' '2010201', '2010000', '2000000', '2000000', '2000000', '2000000', '2000000',
#' '2050220', '2010221', '2010221', '2050221', '2010221', '2010221', '2010221',
#' '2010221', '2010521', '2010251', '2012221', '2012222', '2012222', '2012222',
#' '2012251', '5012221', '2012020', '2022010', '2025010', '2022010', '2022010',
#' '2052010', '2022010', '2022010', '2022010', '2022010', '2021010', '2021010',
#' '2021010', '2051010', '2051010', '5021010', '2001000', '2025010', '2021010',
#' '2021010')
#'
#'df=as.data.frame(df); names(df) = "ch"; df$freq = 1; df$ch = as.character(df$ch)
#'head(df)
#'
#'# 'event' is only for the uncertain observation - not the 'state' observations!
#'# Error in process.data: unused argument (events) if package marked is loaded in session
#'# states 3 & 4 = unobservable
#'
#'# Bill Kendall: For pi, MARK defaults to deriving the probability of the last state listed by subtraction.
#'# Since we like to fix pi for states 3 and 4 to 0, the solution is to list states 1 or 2
#'# last when you specify the states.
#'
#'dp = process.data(df, model = "HidMarkov", strata.labels = c("3", "4","1", "2"),
#' events = c("5"))
#'
#'## create design data
#'ddl = make.design.data(dp)
#'
#'# Set movement parameters (PSI) to 0
#'ddl$Psi$fix[ddl$Psi$stratum=="1" & ddl$Psi$tostratum=="4"]=0
#'ddl$Psi$fix[ddl$Psi$stratum=="2" & ddl$Psi$tostratum=="3"]=0
#'
#'ddl$Psi$fix[ddl$Psi$stratum=="3" & ddl$Psi$tostratum=="4"]=0
#'ddl$Psi$fix[ddl$Psi$stratum=="4" & ddl$Psi$tostratum=="3"]=0
#'
#'# Set unobservable p to 0
#'ddl$p$fix[ddl$p$stratum=="3"]=0
#'ddl$p$fix[ddl$p$stratum=="4"]=0
#'
#'# animals can only enter as breeders, not as post-breeders (unobservable)
#'# In HMMs, the pi parameter is the probability of the event on the first capture.
#'# Note: P(event)! and not as in ESurge, where it is P(state on first capture)!
#'
#'# Set unobservable pi to 0
#'ddl$pi$fix[ddl$pi$stratum=="3"]=0
#'ddl$pi$fix[ddl$pi$stratum=="4"]=0
#'
#'# The Delta parameters are the probability of determining the state given detection
#'# with probability p (MARK defaults to obtaining the probability of correctly assigning
#'# the state by subtraction)
#'ddl$Delta$fix[ddl$Delta$stratum=="3"]=0
#'ddl$Delta$fix[ddl$Delta$stratum=="4"]=0
#'
#'
#'# define a model to fit
#'# make a survival class (successful = post-successful != failed = post-failed)
#'ddl$S$class = ifelse(ddl$S$`3`=="1" |ddl$S$`1`=="1" , 1, 0)
#'
#'Phi.ct = list(formula=~class)
#'p.ct = list(formula=~stratum)
#'Psi.class =list(formula =~ -1+stratum:tostratum) # see page 918 C.17. Multi-strata example in Mark book
#'pi.ct = list(formula=~stratum)
#'Delta.class = list(formula=~stratum)
#'
#'table(ddl$Psi[,c("stratum","tostratum")])
#'table(ddl$pi[,c("stratum")])
#'table(ddl$Delta[,c("event")])
#'
#'
#'# Aim: run a HMM where uncertain events (=5) could be successful breeders or unsuccesful breeders
#'
#'Model.4 = mark(dp, ddl, model.parameters=list(S = Phi.ct,
#' p = p.ct,
#' Psi = Psi.class,
#' pi = pi.ct,
#' Delta = Delta.class),
#' output = FALSE,delete=T, model ="HidMarkov", mlogit0 = T)
#'
#'summary(Model.4)
#'
#'
#'Psilist=get.real(Model.4,"Psi",vcv=T)
#'Psivalues=Psilist$estimates
#'TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)
#'}
NULL
#' Multi-scale dynamic occupancy models in RMark
#' @name RDMultScalOcc
#' @author Connor M. Wood, University of Wisconsin-Madison <cwood9 at wisc.edu>
#' @examples
#' \donttest{
#'## Multi-scale dynamic occupancy models in RMark: a worked example ##
#'
#'# Study design and data structure:
#'# two sessions (i.e., seasons) and 346 sampling sites
#'# up to three secondary sampling periods per season
#'# up to three survey devices per sampling site
#'
#'# import the sample data, RDMultScalOcc.sampledata.csv
#'pathtodata=paste(path.package("RMark"),"extdata",sep="/")
#'dt=read.csv(paste(pathtodata,"RDMultScalOcc.sampledata.csv",sep="/"))
#'dt[is.na(dt)]=0 # replace NAs with 0
#'dt$ch=as.character(dt$ch) # encounter histories (dt$ch) must be characters
#'
#'# note: habitat variables (amount of open forest and average slope) were collected at
#'# two spatial scales
#'# 'sOpen' and 'sSlope' represent the entire sampling site
#'# 'pOpen##' and 'pSlope##' represent conditions relevant to individual devices at each
#'# sampling period
#'# the p-scale variables are coded [name][session][primary];
#'# dt.ddl$p$primary indicates how this should be entered
#'# in this case the values for p$primary varied among devices and between sessions,
#'# but were constant between secondary sampling periods
#'
#'# create the Process Data MARK object
#'dt.pr=process.data(dt,model="RDMultScalOcc",
#' time.intervals=c(0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0),mixtures=3)
#'# note: time.intervals refers to seasons, not secondary sampling periods
#'# note: mixtures refers to the number of devices
#'
#'# create the design data object
#'dt.ddl=make.design.data(dt.pr)
#'# Examine the built-in covariates for each parameter
#'str(dt.ddl)
#'
#'### Approach 1: Manually create (a somewhat arbitrary) set of models #######################
#'# Name variables for simplicity
#'fit.models=function()
#'{
#'Open=list(formula=~sOpen)
#'pOpen=list(formula=~pOpen)
#'Slope=list(formula=~sSlope)
#'pSlope=list(formula=~pSlope)
#'Time=list(formula=~Time) #use carefully because different covariates use 'Time' in different ways
#'
#'null.model=mark(dt.pr,dt.ddl,delete=TRUE)
#'mod1=mark(dt.pr,dt.ddl,model.parameters=list(p=Open),delete=TRUE)
#'#'mod2=mark(dt.pr,dt.ddl,model.parameters=list(p=pOpen),delete=TRUE)
#'mod3=mark(dt.pr,dt.ddl,model.parameters=list(p=Slope),delete=TRUE)
#'mod4=mark(dt.pr,dt.ddl,model.parameters=list(p=pSlope),delete=TRUE)
#'mod5=mark(dt.pr,dt.ddl,model.parameters=list(Psi=Open),delete=TRUE)
#'mod6=mark(dt.pr,dt.ddl,model.parameters=list(Psi=Slope),delete=TRUE)
#'mod7=mark(dt.pr,dt.ddl,model.parameters=list(Gamma=Open),delete=TRUE)
#'mod8=mark(dt.pr,dt.ddl,model.parameters=list(Epsilon=Slope),delete=TRUE)
#'mod9=mark(dt.pr,dt.ddl,model.parameters=list(Theta=list(formula=~time)),delete=TRUE)
#'
#'return(collect.models()) # View results sorted by AICc
#'}
#'results=fit.models()
#'###########################################################################################
#'
#'### Approach 2: use mark.wrapper to create your models ####################################
#'fit.models=function()
#'{
#' # Models for p
#' p.null=list(formula=~1)
#' p.open=list(formula=~sOpen)
#' p.popen=list(formula=~pOpen)
#' p.slope=list(formula=~sSlope)
#' p.pslope=list(formula=~pSlope)
#' # Models for Psi
#' Psi.null=list(formula=~1)
#' Psi.open=list(formula=~sOpen)
#' Psi.Slope=list(formula=~sSlope)
#' #Model for Gamma
#' Gamma.null=list(formula=~1)
#' Gamma.open=list(formula=~sOpen)
#' # Model for Epsilon
#' Epsilon.null=list(formula=~1)
#' Epsilon.slope=list(formula=~sSlope)
#' # Model for Theta
#' # note: 'time' is defined in dt.ddl$Theta (use str(dt) to see all predefined variables)
#' Theta.null=list(formula=~1)
#' Theta.time=list(formula=~time)
#' # create all combinations of these sub-models
#' cml=create.model.list("RDMultScalOcc")
#' results=mark.wrapper(cml,data=dt.pr,ddl=dt.ddl,delete=TRUE)
#' return(results)
#'}
#'
#'fit.models() # creates all combinations of variables and prints AICc table
#'results
#'}
NULL
#' An example of the Multistrata (multi-state) model in which states are routes taken by migrating fish.
#'
#' @name skagit
#' @author Megan Moore <megan.moore at noaa.gov>
#' @examples
#'# There are just two states which correspond to route A and route B. There are 6 occasions
#'# which are the locations rather than times. After release at 1=A there is no movement
#'# between states for the first segment, fish are migrating downriver together and all pass 2A.
#'# Then after occasion 2, migrants go down the North Fork (3A) or the South Fork (3B),
#'# which both empty into Skagit Bay. Once in saltwater, they can go north to Deception Pass (4A)
#'# or South to a receiver array exiting South Skagit Bay (4B). Fish in route A can then only go
#'# to the Strait of Juan de Fuca, while fish in route B must pass by Admiralty Inlet (5B).
#'# Then both routes end with the array at the Strait of Juan de Fuca.
#'#
#'# 1A
#'# |
#'# 2A
#'# / \
#'# 3A 3B
#'# / \ / \
#'# 4A 4B 4A 4B
#'# | \ / |
#'# 5A 5B 5A 5B
#'# \ \ / /
#'# 6
#'#
#'# from 3A and 3B they can branch to either 4A or 4B; branches merge at 6
#'# 5A does not exist so p=0; only survival from 4A to 6 can be
#'# estimated which is done by setting survival from 4A to 5A to 1 and
#'# estimating survival from 5A to 6 which is then total survival from 4A to 6.
#'\donttest{
#'pathtodata=paste(path.package("RMark"),"extdata",sep="/")
#'skagit=import.chdata(paste(pathtodata,"skagit.txt",sep="/"),field.types=c("f"),header=TRUE)
#'skagit.processed=process.data(skagit,model="Multistrata",groups=c("tag"))
#'skagit.ddl=make.design.data(skagit.processed)
#'#
#'# p
#'#
#'# Can't be seen at 5A or 2B,6B (the latter 2 don't exist)
#'skagit.ddl$p$fix=ifelse((skagit.ddl$p$stratum=="A"&skagit.ddl$p$time==5) |
#' (skagit.ddl$p$stratum=="B"&skagit.ddl$p$time%in%c(2,6)),0,NA)
#'# Estimated externally from current data to allow estimation of survival at last interval
#'skagit.ddl$p$fix[skagit.ddl$p$tag=="v7"&skagit.ddl$p$time==6&skagit.ddl$p$stratum=="A"]=0.687
#'skagit.ddl$p$fix[skagit.ddl$p$tag=="v9"&skagit.ddl$p$time==6&skagit.ddl$p$stratum=="A"]=0.975
#'#
#'# Psi
#'#
#'# only 3 possible transitions are A to B at time interval 2 to 3 and
#'# for time interval 3 to 4 from A to B and from B to A
#'# rest are fixed values
#'skagit.ddl$Psi$fix=NA
#'# stay in A for intervals 1-2, 4-5 and 5-6
#'skagit.ddl$Psi$fix[skagit.ddl$Psi$stratum=="A"&
#' skagit.ddl$Psi$tostratum=="B"&skagit.ddl$Psi$time%in%c(1,4,5)]=0
#'# stay in B for interval 4-5
#'skagit.ddl$Psi$fix[skagit.ddl$Psi$stratum=="B"&skagit.ddl$Psi$tostratum=="A"
#' &skagit.ddl$Psi$time==4]=0
#'# leave B to go to A for interval 5-6
#'skagit.ddl$Psi$fix[skagit.ddl$Psi$stratum=="B"&skagit.ddl$Psi$tostratum=="A"&
#' skagit.ddl$Psi$time==5]=1
#'# "stay" in B for interval 1-2 and 2-3 because none will be in B
#'skagit.ddl$Psi$fix[skagit.ddl$Psi$stratum=="B"&skagit.ddl$Psi$tostratum=="A"&
#' skagit.ddl$Psi$time%in%1:2]=0
#'#
#'# S
#'#
#'# None in B, so fixing S to 1
#'skagit.ddl$S$fix=ifelse(skagit.ddl$S$stratum=="B"&skagit.ddl$S$time%in%c(1,2),1,NA)
#'skagit.ddl$S$fix[skagit.ddl$S$stratum=="A"&skagit.ddl$S$time==4]=1
#'# fit model
#'p.timexstratum.tag=list(formula=~time:stratum+tag,remove.intercept=TRUE)
#'Psi.sxtime=list(formula=~-1+stratum:time)
#'S.stratumxtime=list(formula=~-1+stratum:time)
#'#
#'S.timexstratum.p.timexstratum.Psi.sxtime=mark(skagit.processed,skagit.ddl,
#' model.parameters=list(S=S.stratumxtime,p= p.timexstratum.tag,Psi=Psi.sxtime),delete=TRUE)
#'# calculation of cummulative survival for entire route
#'Sest=plogis(coef(S.timexstratum.p.timexstratum.Psi.sxtime)$estimate)
#'# A
#'prod(Sest[c(1:3,6)])
#'#[1] 0.1644
#'# B
#'prod(Sest[c(1,2,4,5,7)])
#'#[1] 0.1154
#' }
NULL
#' Burnham Live-Dead Model
#'
#' An example of the Burnham live-dead model using simulated data LD1.inp from Chapter 9 of Cooch and White
#'
#' @name Burnham
#' @author Luke Eberhart-Phillips<luke.eberhart at gmail.com>
#' @examples
#' \donttest{
#'###############################################################################
#'#### RMARK script for conducting the Burnham model tutorial in Chapter 9.3 ####
#'#################### the of the Cooch and White MARK book #####################
#'###############################################################################
#'##################### Code by: Luke Eberhart-Phillips #########################
#'###### Dept. Animal Behaviour, Bielefeld University, Bielefeld, Germany #######
#'##################### email: luke.eberhart at gmail.com ##########################
#'###############################################################################
#'
#'# import/convert the simulated "LD1.inp" MARK capture history into an RMARK
#'# dataframe, while defining the two groups as "Y" for individuals marked as
#'# young, and "A" for individuals marked as adults
#'# NOTE: the "LD1.inp" file is found in the zipped folder downloaded when you
#'# click on "Example data files" in the drop-down menu of the MARK book webpage
#'# \url{http://www.phidot.org/software/mark/docs/book/}
#' pathtodata=paste(path.package("RMark"),"extdata",sep="/")
#' LD=convert.inp(paste(pathtodata,"ld1",sep="/"),
#' group.df=data.frame(age_marked=c("Y","A")))
#'# process the data by defining the model type as "Burnham" and the groups in
#'# the data. In this case the only group is the age at which individuals were
#'# marked
#'LD.proc=process.data(data = LD,
#' model = "Burnham",
#' groups=c("age_marked"),
#' age.var=1,
#' initial.age=c(1,0))
#'
#'# make the design data from the process data above
#'LD.ddl=make.design.data(LD.proc)
#'
#'# add the correct binning to the design data so that individuals that were
#'# marked as young are adults in their second year of life, where as those
#'# marked as adults are adults for their entire life.
#'LD.ddl=add.design.data(data = LD.proc,
#' ddl = LD.ddl,
#' parameter="S",
#' type = "age",
#' bins = c(0,1,8),
#' right = FALSE,
#' name = "age",
#' replace = TRUE)
#'
#'# do the same to the F parameter
#'LD.ddl=add.design.data(data = LD.proc,
#' ddl = LD.ddl,
#' parameter="F",
#' type = "age",
#' bins = c(0,1,8),
#' right = FALSE,
#' name = "age",
#' replace = TRUE)
#'
#'# check parameter matrix to see if groups were binned correctly in the S matrix
#'PIMS(mark(data = LD.proc,
#' ddl = LD.ddl,
#' model.parameters=list(S=list(formula=~age)),
#' delete=TRUE,
#' model = "Burnham"),
#' "S")
#'
#'# Create the formulas that describe variation in the parameter we want to test.
#'# In this case we want to test for an age effect on survival and fidelity,
#'# while keeping recapture and recovery probabilities constant.
#'S.age=list(formula=~age) # S(age)
#'p.dot=list(formula=~1) # p(.)
#'F.age=list(formula=~age) # F(age)
#'r.dot=list(formula=~1) # r(.)
#'
#'# Run the model
#'LD.model.age.F.S=mark(data = LD.proc,
#' ddl = LD.ddl,
#' model.parameters = list(S = S.age, p = p.dot,
#' F =F.age, r = r.dot),
#' invisible = FALSE,
#' model = "Burnham",delete=TRUE)
#'
#'# Check the parameter estimates, they should be the same as those generated
#'# when doing the tutorial in chapter 9.3 of the in MARK Book (table on pg 9-8)
#'LD.model.age.F.S$results$real
#'
#'# Clean your working directory
#'cleanup(ask=FALSE)
#'}
NULL
#' Lark Sparrow
#'
#' An example of Multiple Scale Occupancy model for some lark sparrow data that was contributed by David Pavlacky
#' at Rocky Mountain bird observatory. The study design was a GRTS selection of paired "Deferred" and "Grazed" pastures.
#' The point count locations within each pasture were a random selection of systematic point count locations separated by 250 m. Each point count had
#' a radius of 125m. A removal design was used to the set the data to missing after the first detection.
#' The point count data were set to missing when fewer than nine points were surveyed.
#'
#' @name larksparrow
#' @aliases LASP
#' @docType data
#' @format A data frame with 52 observations on the following 20 variables.
#' \describe{\item{ch}{a character vector containing the encounter history}
#' \item{ceap}{a factor with two levels "Deferred" and "Grazed" corresponding to a rest rotation grazing system with pastures either rested (Deferred) or grazed (Grazed) during the spring bird breeding season.}
#' \item{cwx}{a continuous covariate for primary occasion x, representing an ocular estimate of the proportion of area covered by crested wheatgrass in a 50-m radius around the point count location.}
#' \item{tdx}{ a continuous covariate for primary occasion x, representing the starting time (h) of each 6-min point count survey measured on the ratio scale (1.5 h = 1 h 30 min).}
#' }
#' @keywords datasets
#' @examples
#' \donttest{
#'
#' # This example is excluded from testing to reduce package check time
#'# Create dataframe
#'data(LASP)
#' mscale=LASP
#'
#'# Process data with MultScalOcc model and use group variables
#'
#'mscale.proc=process.data(mscale,model="MultScalOcc",groups=c("ceap"),begin.time=1,mixtures=3)
#'
#'# Create design data
#'
#'ddl=make.design.data(mscale.proc)
#'
#'# Create function to build models
#'
#'do.Species=function()
#'{
#' p.1=list(formula=~1)
#' p.2=list(formula=~ceap)
#' p.3=list(formula=~td)
#'
#' Theta.1=list(formula=~1)
#' Theta.2=list(formula=~ceap)
#' Theta.3=list(formula=~cw)
#'
#' Psi.1=list(formula=~1)
#' Psi.2=list(formula=~ceap)
#'
#' cml=create.model.list("MultScalOcc")
#' return(mark.wrapper(cml,data=mscale.proc,ddl=ddl,adjust=FALSE,realvcv=TRUE,delete=TRUE))
#'}
#'
#'# Run function to get results
#'
#'Species.results=do.Species()
#'
#'# Output model table and estimates
#'
#'Species.results$model.table
#'
#'Species.results[[as.numeric(rownames(Species.results$model.table[1,]))]]$results$real
#'Species.results[[as.numeric(rownames(Species.results$model.table[1,]))]]$results$beta
#'
#'#write.csv(Species.results$model.table,file="lasp_model_selection.csv",row.names=FALSE)
#'
#'#write.csv(Species.results[[as.numeric(rownames(Species.results$model.table[1,]))]]$results$real,
#'# file="lasp_m1_real.csv")
#'#write.csv(Species.results[[as.numeric(rownames(Species.results$model.table[1,]))]]$results$beta,
#'# file="lasp_m1_beta.csv")
#'
#'# Covariate prediction and model averaging
#'
#'# p(time of day)
#'
#'mintd <- min(mscale[,12:20])
#'maxtd <- max(mscale[,12:20])
#'td.values <- mintd+(0:100)*(maxtd-mintd)/100
#'
#'PIMS(Species.results[[1]],"p",simplified=FALSE)
#'
#'td <- covariate.predictions(Species.results,data=data.frame(td1=td.values),indices=c(21))
#'
#'#write.table(td$estimates,file="lasp_cov_pred_p_td.csv",sep=",",col.names=TRUE,
#'# row.names=FALSE)
#'
#'# Theta(crested wheatgrass cover)
#'
#'mincw <- min(mscale[,3:11])
#'maxcw <- max(mscale[,3:11])
#'cw.values <- mincw+(0:100)*(maxcw-mincw)/100
#'
#'PIMS(Species.results[[1]],"Theta",simplified=FALSE)
#'
#'cw <- covariate.predictions(Species.results,data=data.frame(cw1=cw.values),indices=c(3))
#'
#'#write.table(cw$estimates,file="lasp_cov_pred_theta_cw.csv",sep=",",col.names=TRUE,
#'# row.names=FALSE)
#'
#'# Psi(ceap grazing for wildlife practice)
#'
#'ceap.values <- as.data.frame(matrix(c(1,2),ncol=1))
#'names(ceap.values) <- c("index")
#'
#'PIMS(Species.results[[1]],"Psi",simplified=FALSE)
#'
#'ceap <- covariate.predictions(Species.results,data=data.frame(ceap=ceap.values))
#'
#'#write.table(ceap$estimates,file="lasp_cov_pred_psi_ceap.csv",sep=",",col.names=TRUE,
#'# row.names=FALSE)
#'
#' }
NULL
#' San Luis Valley mallard data
#'
#' A recovery data set for mallards in San Luis Valley Colorado
#'
#' This is a data set that accompanies program MARK as an example for Brownie and
#' Seber recovery model. In those input files it is in a summarized format. Here
#' it is in the LD encounter history format. The data can be stratified using ReleaseAge (Adult, Young)
#' as a grouping variable.
#'
#' Note that in the MARK example the variable is named Age. In the R code, the
#' fields "age" and "Age" have specific meanings in the design data related to
#' time since release. These will override the use of a field with the same
#' name in the individual covariate data, so the names "time", "Time",
#' "cohort", "Cohort", "age", and "Age" should not be used in the individual
#' covariate data with possibly the exception of "cohort" which is not defined
#' for models with "Square" PIMS such as POPAN and other Jolly-Seber type
#' models.
#'
#' @name brownie
#' @docType data
#' @format A data frame with 108 observations on the following 5 variables.
#' \describe{ \item{ch}{a character vector containing the encounter history of
#' each bird} \item{freq}{frequency of that encounter history} \item{ReleaseAge}{the age of the bird when it was
#' released}}
#' @keywords datasets
#' @examples
#' # brownie=import.chdata("brownie.inp",field.types=c("n","f"))
#' \donttest{
#' # This example is excluded from testing to reduce package check time
#' data(brownie)
#' # default ordering of ReleaseAge is alphabetic so it is
#' # Adult, Young which is why initial.ages=c(1,0)
#' # Seber Recovery
#' br=process.data(brownie,model="Recovery",groups="ReleaseAge",age.var=1,initial.ages=c(1,0))
#' br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
#' r=list(age.bins=c(0,1,10))),right=FALSE)
#' mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~-1+age:time,link="sin"),
#' r=list(formula=~-1+age:time,link="sin")),delete=TRUE)
#' # Brownie Recovery
#' br=process.data(brownie,model="Brownie",groups="ReleaseAge",age.var=1,initial.ages=c(1,0))
#' br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
#' f=list(age.bins=c(0,1,10))),right=FALSE)
#' mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~-1+age:time,link="sin"),
#' f=list(formula=~-1+age:time,link="sin")),delete=TRUE)
#' mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~-1+age,link="sin"),
#' f=list(formula=~-1+age,link="sin")),delete=TRUE)
#' #Random effects Seber recovery (not run to save computation time)
#' #br=process.data(brownie,model="REDead",groups="ReleaseAge",age.var=1,initial.ages=c(1,0))
#' #br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
#' # r=list(age.bins=c(0,1,10))),right=FALSE)
#' #mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~age),r=list(formula=~age)),delete=TRUE)
#' #Pledger Mixture Seber recovery
#' br=process.data(brownie,model="PMDead",groups="ReleaseAge",
#' mixtures=3,age.var=1,initial.ages=c(1,0))
#' br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
#' r=list(age.bins=c(0,1,10))),right=FALSE)
#' mod=mark(br,br.ddl,model.parameters=list(pi=list(formula=~mixture),
#' S=list(formula=~age+mixture),r=list(formula=~age)),delete=TRUE)
#' br=process.data(brownie,model="PMDead",groups="ReleaseAge",
#' mixtures=2,age.var=1,initial.ages=c(1,0))
#' br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
#' r=list(age.bins=c(0,1,10))),right=FALSE)
#' mod=mark(br,br.ddl,model.parameters=list(pi=list(formula=~age),
#' S=list(formula=~age+mixture),r=list(formula=~age)),delete=TRUE)
#' }
NULL
#' Black duck known fate data
#'
#' A known fate data set on Black ducks that accompanies MARK as an example
#' analysis using the Known model.
#'
#' This is a data set that accompanies program MARK as an example for Known
#' fate. The data can be stratified using BirdAge as a grouping variable. The
#' function \code{run.Blackduck} defined below in the examples creates some of
#' the models used in the dbf file that accompanies MARK.
#'
#' Note that in the MARK example the variable is named Age. In the R code, the
#' fields "age" and "Age" have specific meanings in the design data related to
#' time since release. These will override the use of a field with the same
#' name in the individual covariate data, so the names "time", "Time",
#' "cohort", "Cohort", "age", and "Age" should not be used in the individual
#' covariate data with possibly the exception of "cohort" which is not defined
#' for models with "Square" PIMS such as POPAN and other Jolly-Seber type
#' models.
#'
#' @name Blackduck
#' @docType data
#' @format A data frame with 48 observations on the following 5 variables.
#' \describe{ \item{ch}{a character vector containing the encounter history of
#' each bird} \item{BirdAge}{the age of the bird: a factor with levels \code{0}
#' \code{1} for young and adult} \item{Weight}{the weight of the bird at time
#' of marking} \item{Wing_Len}{the wing-length of the bird at time of marking}
#' \item{condix}{the condition index of the bird at time of marking} }
#' @keywords datasets
#' @examples
#'
#' data(Blackduck)
#' # Change BirdAge to numeric; starting with version 1.6.3 factor variables are
#' # no longer allowed. They can work as in this example but they can be misleading
#' # and fail if the levels are non-numeric. The real parameters will remain
#' # unchanged but the betas will be different.
#' Blackduck$BirdAge=as.numeric(Blackduck$BirdAge)-1
#' run.Blackduck=function()
#' {
#' #
#' # Process data
#' #
#' bduck.processed=process.data(Blackduck,model="Known")
#' #
#' # Create default design data
#' #
#' bduck.ddl=make.design.data(bduck.processed)
#' #
#' # Add occasion specific data min < 0; I have no idea what it is
#' #
#' bduck.ddl$S$min=c(4,6,7,7,7,6,5,5)
#' #
#' # Define range of models for S
#' #
#' S.dot=list(formula=~1)
#' S.time=list(formula=~time)
#' S.min=list(formula=~min)
#' S.BirdAge=list(formula=~BirdAge)
#' #
#' # Note that in the following model in the MARK example, the covariates
#' # have been standardized. That means that the beta parameters will be different
#' # for BirdAge, Weight and their interaction but the likelihood and real parameter
#' # estimates are the same.
#' #
#' S.BirdAgexWeight.min=list(formula=~min+BirdAge*Weight)
#' S.BirdAge.Weight=list(formula=~BirdAge+Weight)
#' #
#' # Create model list and run assortment of models
#' #
#' model.list=create.model.list("Known")
#' bduck.results=mark.wrapper(model.list,data=bduck.processed,ddl=bduck.ddl,
#' invisible=FALSE,threads=1,delete=TRUE)
#'
#' #
#' # Return model table and list of models
#' #
#' return(bduck.results)
#' }
#' bduck.results=run.Blackduck()
#' bduck.results
#'
#'
#'
NULL
#' White-tailed deer double observer spotlight capture-recapture analysis
#'
#' This data represents a set of independent double observer road-transect survey data of white-tailed deer
#' on Brosnan Forest, South Carolina surveyed in August, 2005-2009. The primary reason for
#' this package is to provide a completely reproducible example of the analysis from Collier et al. (2012). We used a
#' Huggins closed capture model implemented in MARK \url{http://www.phidot.org/software/mark/} via RMark
#' both of which will need to be installed on the system to use this package. The data have 2 time periods (primary observer (t1) was a thermal imager, secondary observer (t2) was
#' a spotlight observer in the same vehicle on the same side) with the primary objective of the study being to evaluate
#' the detection (recapture) rates of white-tailed deer using spotlights as a survey method.
#'
#' @details In addition to detailing the analysis used by Collier et al. (2012), this example documents the
#' use of the \code{share} argument in the RMark parameter specification because there is presently very little
#' documentation on the use of \code{share}. Parameters in MARK models rarely share columns of the design matrix. For
#' example while you might want to use the same covariate for survival and capture probability, you would never use the
#' same beta (same column of the design matrix) for each parameter. However, there are exceptions when the parameters
#' represent similar quantities and that is when the \code{share} argument is useful. For example, in the closed capture models
#' p is initial capture probability and c is recapture probability. In this case, it would make perfect sense to use the same
#' column of the design matrix for both parameters. The most obvious case is to fit a model in which p=c.
#'
#' In RMark, certain pairs of parameters have been identified as similar and shareable. These can be found in the file parameters.txt
#' which is in the RMark directory in your R library. With each pair that is shareable, the first one listed is the primary parameter.
#' When you want to share columns in the design matrix, share=TRUE is added to the specification of the primary parameter. A parameter
#' specification is not given for the other secondary parameter when they are shared. When RMark, sees that the parameters are to be shared it
#' creates a pooled set of design data and adds a column with the name of the secondary parameter and its value is 0 for the rows
#' for the primary parameter and 1 for the rows for the secondary parameter. For example, with the closed capture model if share=TRUE is
#' added to the parameter specification for p, a model is not specified for c, and the pooled design data set contains a field called c.
#' The added field allows construction of models where there are restricted differences between the parameters. For example, p=list(formula=~time+c,share=TRUE)
#' will fit a model in which capture probability varies by time and recapture probability includes an additive difference on the link scale.
#' Because the design data are pooled when you share parameters, if you modify design data for one of the parameters, the other most be modified as
#' as well, so the columns of the design data for both parameters are the same or RMark will give an error.
#'
#' The argument \code{share} is used in all the candidate models in the below example analysis. As a simplified example of how
#' \code{share} works, look at the candidate models in the \code{bfrun{}} function call named \code{mod.2} and \code{mod.2a} (note that
#' \code{mod.2a} was not included in the supplemental file available from the Journal of Wildlife Management and is only included in this package). Both
#' of these models are conducting the exact same analysis, with the first \code{mod.2}, we used the formula \code{~time} (if you don't
#' know what this means go read the MARKBOOK at \url{http://www.phidot.org/software/mark/}. Notice, however, we used the argument
#' \code{share} in \code{mod.2}, which tells RMark to share columns of the MARK design matrix. For comparison, so you can evaluate how
#' \code{share} works for yourself, \code{mod.2a} recreates the same analysis as \code{mod.2}, but uses the approach more typical to MARK
#' analyses where each parameter is specified independently and uniquely.
#'
#' @name deer
#' @docType data
#' @format The format is a data frame with 4508 observations on the following 7 variables.
#' \describe{ \item{SL (spotlight)}{0/1 whether deer was missed/seen by the spotlight observer}
#' \item{TI (thermal imager)}{0/1 whether deer was missed/seen by the thermal imager observer}
#' \item{Group}{Factor with 79 levels representing each unique paired (TI-SL) survey conducted}
#' \item{Year}{Factor with 5 levels for year of survey}
#' \item{MaxCount}{Count of maximum number of deer seen for each survey, only needed
#' for bootstrapp analysis in MARK, not used in bfdeeR package}
#' \item{Cluster}{Value assigning each deer to a specific observation cluster, only needed
#' for bootstrapp analysis in MARK, not used in bfdeeR package}
#' \item{MgmtUnit}{Management unit identification} }
#' @references Collier, B. A., S. S. Ditchkoff, J. B. Raglin, and C. R. Ruth. 2012. Spotlight surveys
#' for white-tailed deer: monitoring panacea or exercise in futility? Journal of Wildlife Management, In Press.
#' @keywords datasets
#' @author Bret Collier
#' @examples
#'
#' \donttest{
#' # This example is excluded from testing to reduce package check time
#' data(deer)
#' x=data.frame(ch=paste(deer$TI, deer$SL, sep=""), Survey=factor(deer$Group),
#' Year=factor(deer$Year), Cluster=deer$Cluster, MgtUnit=factor(deer$MgmtUnit))
#' x$ch=as.character(x$ch)
#' bfrun=function(){
#' x.proc=process.data(x, model="Huggins", groups=c("Survey", "Year", "MgtUnit"))
#' x.ddl=make.design.data(x.proc)
#'
#' #Silly Null model, constant p & c sharing 1 parameter (one detection estimate)
#' p.shared=list(formula=~1,share=TRUE)
#' mod.1=mark(x.proc, x.ddl, model.parameters=list(p=p.shared), invisible=FALSE,delete=TRUE)
#'
#' #2 Parameter Null Model, constant p, constant c, different p and c (one estimate for each; p ne c)
#' #p(time), c(-), share=TRUE, detection is time dependent, with recapture parameter shared
#' p.sharetime=list(formula=~time, share=TRUE)
#' mod.2=mark(x.proc, x.ddl, model.parameters=list(p=p.sharetime), invisible=FALSE,delete=TRUE)
#'
#' #2a Parameter Null Model, constant p, constant c,
#' # different p and c (one estimate for each; p ne c) not using share
#' mod.2a=mark(x.proc, x.ddl, model.parameters=list(p=list(formula=~1), c=list(formula=~1)),
#' delete=TRUE)
#'
#' #Fully parameterized model, different p and c for each survey transect replicate,
#' # management unit, method (TI or SL) and any observers
#' p.survey=list(formula=~Survey*time, share=TRUE)
#' mod.3=mark(x.proc, x.ddl, model.parameters=list(p=p.survey), invisible=FALSE,delete=TRUE)
#'
#' #p(MU), c(MU), initial detection and recapture differ and are management unit dependent
#' p.mu=list(formula=~MgtUnit*time, share=TRUE)