forked from joshuacfowler/Endodemog_2019
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathendodemog_figures.R
6774 lines (5730 loc) · 646 KB
/
endodemog_figures.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
## Title: Grass endophyte vital rate variance with a bayesian framework
## Purpose: Visualizes output of endophyte effect on year variance from
## model outputs for Survival, Growth and Flowering models
## Authors: Joshua and Tom
library(tidyverse)
library(rstan)
library(StanHeaders)
library(shinystan)
library(bayesplot)
library(devtools)
library(reshape2)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(RColorBrewer)
invlogit<-function(x){exp(x)/(1+exp(x))}
logit = function(x) { log(x/(1-x)) }
# source in raw output from models
## Read in the survival model output for all species
survPOAL <- read_rds(path = "/Users/joshuacfowler/Dropbox/EndodemogData/Model_Runs/endodemog_surv_POAL_withplot.rds")
survPOSY <- read_rds(path = "/Users/joshuacfowler/Dropbox/EndodemogData/Model_Runs/endodemog_surv_POSY_withplot.rds")
survLOAR <- readRDS(file = "model_run_MAR7/endodemog_surv_LOAR_withplot.rds")
survELVI <- readRDS(file = "model_run_MAR7/endodemog_surv_ELVI_withplot.rds")
survELRI <- readRDS(file = "model_run_MAR7/endodemog_surv_ELRI_withplot.rds")
survFESU <- readRDS(file = "model_run_MAR7/endodemog_surv_FESU_withplot.rds")
survAGPE <- readRDS(file = "model_run_MAR7/endodemog_surv_AGPE_withplot.rds")
## Read in the flower model output for all species
flwPOAL <- readRDS(file = "model_run_MAR7/endodemog_flw_POAL_withplot.rds")
flwPOSY <- readRDS(file = "model_run_MAR7/endodemog_flw_POSY_withplot.rds")
flwLOAR <- readRDS(file = "model_run_MAR7/endodemog_flw_LOAR_withplot.rds")
flwELVI <- readRDS(file = "model_run_MAR7/endodemog_flw_ELVI_withplot.rds")
flwELRI <- readRDS(file = "model_run_MAR7/endodemog_flw_ELRI_withplot.rds")
flwFESU <- readRDS(file = "model_run_MAR7/endodemog_flw_FESU_withplot.rds")
flwAGPE <- readRDS(file = "model_run_MAR7/endodemog_flw_AGPE_withplot.rds")
## Read in the growth model output for all species
growPOAL <- readRDS(file = "model_run_MAR7/endodemog_grow_POAL_withplot.rds")
growPOSY <- readRDS(file = "model_run_MAR7/endodemog_grow_POSY_withplot.rds")
growLOAR <- readRDS(file = "model_run_MAR7/endodemog_grow_LOAR_withplot.rds")
growELVI <- readRDS(file = "model_run_MAR7/endodemog_grow_ELVI_withplot.rds")
growELRI <- readRDS(file = "model_run_MAR7/endodemog_grow_ELRI_withplot.rds")
growFESU <- readRDS(file = "model_run_MAR7/endodemog_grow_FESU_withplot.rds")
growAGPE <- readRDS(file = "model_run_MAR7/endodemog_grow_AGPE_withplot.rds")
# save posteriors and summaries within dataframes
params = c("beta", "tau_year", "sigma_e", "tau_plot")
post_survPOAL <- as.data.frame(survPOAL, pars = params)
post_survPOSY <- as.data.frame(survPOSY, pars = params)
post_survLOAR <- as.data.frame(survLOAR, pars = params)
post_survELVI <- as.data.frame(survELVI, pars = params)
post_survELRI <- as.data.frame(survELRI, pars = params)
post_survFESU <- as.data.frame(survFESU, pars = params)
post_survAGPE <- as.data.frame(survAGPE, pars = params)
sum_survPOAL <- summary(survPOAL, pars = params)
sum_survPOSY <- summary(survPOSY, pars = params)
sum_survLOAR <- summary(survLOAR, pars = params)
sum_survELVI <- summary(survELVI, pars = params)
sum_survELRI <- summary(survELRI, pars = params)
sum_survFESU <- summary(survFESU, pars = params)
sum_survAGPE <- summary(survAGPE, pars = params)
post_flwPOAL <- as.data.frame(flwPOAL, pars = params)
post_flwPOSY <- as.data.frame(flwPOSY, pars = params)
post_flwLOAR <- as.data.frame(flwLOAR, pars = params)
post_flwELVI <- as.data.frame(flwELVI, pars = params)
post_flwELRI <- as.data.frame(flwELRI, pars = params)
post_flwFESU <- as.data.frame(flwFESU, pars = params)
post_flwAGPE <- as.data.frame(flwAGPE, pars = params)
sum_flwPOAL <- summary(flwPOAL, pars = params)
sum_flwPOSY <- summary(flwPOSY, pars = params)
sum_flwLOAR <- summary(flwLOAR, pars = params)
sum_flwELVI <- summary(flwELVI, pars = params)
sum_flwELRI <- summary(flwELRI, pars = params)
sum_flwFESU <- summary(flwFESU, pars = params)
sum_flwAGPE <- summary(flwAGPE, pars = params)
post_growPOAL <- as.data.frame(growPOAL, pars = params)
post_growPOSY <- as.data.frame(growPOSY, pars = params)
post_growLOAR <- as.data.frame(growLOAR, pars = params)
post_growELVI <- as.data.frame(growELVI, pars = params)
post_growELRI <- as.data.frame(growELRI, pars = params)
post_growFESU <- as.data.frame(growFESU, pars = params)
post_growAGPE <- as.data.frame(growAGPE, pars = params)
sum_growPOAL <- summary(growPOAL, pars = params)
sum_growPOSY <- summary(growPOSY, pars = params)
sum_growLOAR <- summary(growLOAR, pars = params)
sum_growELVI <- summary(growELVI, pars = params)
sum_growELRI <- summary(growELRI, pars = params)
sum_growFESU <- summary(growFESU, pars = params)
sum_growAGPE <- summary(growAGPE, pars = params)
#Looking at different color palettes
# Classic palette BuPu, with 4 colors
purps = brewer.pal(9, "Purples")
bupu = brewer.pal(9,"BuPu")
oran = brewer.pal(9, "Oranges")
brewer.pal(9, "Oranges")
orrd = brewer.pal(9, "OrRd")
# I can add more tones to this palette :
purps = colorRampPalette(purps)(11)
# pie(rep(1, length(purps)), col = purps , main="")
bupu = colorRampPalette(bupu)(11)
# pie(rep(1, length(bupu)), col = bupu , main="")
oran = colorRampPalette(oran)(11)
# pie(rep(1, length(oran)), col = oran , main="")
orrd = colorRampPalette(orrd)(11)
# pie(rep(1, length(orrd)), col = orrd , main="")
yearcolors <- c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99",
"#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a",
"#ffff99")
colors2 <- c("#ff7f00","#6a3d9a") #"#ff7f00" = E-, "#e31a1c" = E+
# Variance posterior distributions by species
# survival
AGPE_s <- ggplot() +
geom_histogram(data = post_survAGPE, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_survAGPE, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "Posterior density", title = "Variance Estimate") + theme_classic()
ELRI_s <- ggplot() +
geom_histogram(data = post_survELRI, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_survELRI, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "Posterior density", title = "Variance Estimate") + theme_classic()
ELVI_s <- ggplot() +
geom_histogram(data = post_survELVI, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_survELVI, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "Posterior density", title = "Variance Estimate") + theme_classic()
FESU_s <- ggplot() +
geom_histogram(data = post_survFESU, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_survFESU, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "Posterior density", title = "Variance Estimate") + theme_classic()
LOAR_s <- ggplot() +
geom_histogram(data = post_survLOAR, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_survLOAR, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "Posterior density", title = "Variance Estimate") + theme_classic()
POAL_s <- ggplot() +
geom_histogram(data = post_survPOAL, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_survPOAL, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "Posterior density", title = "Variance Estimate") + theme_classic()
POSY_s <- ggplot() +
geom_histogram(data = post_survPOSY, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_survPOSY, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "Posterior density", title = "Variance Estimate") + theme_classic()
# flowering
AGPE_f <- ggplot() +
geom_histogram(data = post_flwAGPE, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_flwAGPE, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "AGPE") + theme_classic()
ELRI_f <- ggplot() +
geom_histogram(data = post_flwELRI, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_flwELRI, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "ELRI") + theme_classic()
ELVI_f <- ggplot() +
geom_histogram(data = post_flwELVI, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_flwELVI, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "ELVI") + theme_classic()
FESU_f <- ggplot() +
geom_histogram(data = post_flwFESU, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_flwFESU, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "FESU") + theme_classic()
LOAR_f <- ggplot() +
geom_histogram(data = post_flwLOAR, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_flwLOAR, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "LOAR") + theme_classic()
POAL_f <- ggplot() +
geom_histogram(data = post_flwPOAL, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_flwPOAL, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "POAL") + theme_classic()
POSY_f <- ggplot() +
geom_histogram(data = post_flwPOSY, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_flwPOSY, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "POSY") + theme_classic()
# growth
AGPE_g <- ggplot() +
geom_histogram(data = post_growAGPE, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_growAGPE, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "AGPE") + theme_classic()
ELRI_g <- ggplot() +
geom_histogram(data = post_growELRI, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_growELRI, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "ELRI") + theme_classic()
ELVI_g <- ggplot() +
geom_histogram(data = post_growELVI, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_growELVI, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "ELVI") + theme_classic()
FESU_g <- ggplot() +
geom_histogram(data = post_growFESU, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_growFESU, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "FESU") + theme_classic()
LOAR_g <- ggplot() +
geom_histogram(data = post_growLOAR, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .5, bins = 200) +
geom_histogram(data = post_growLOAR, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .5, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "LOAR") + theme_classic()
POAL_g <- ggplot() +
geom_histogram(data = post_growPOAL, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_growPOAL, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "POAL") + theme_classic()
POSY_g <- ggplot() +
geom_histogram(data = post_growPOSY, aes(`sigma_e[1]`), fill = "#ff7f00", alpha = .4, bins = 200) +
geom_histogram(data = post_growPOSY, aes(`sigma_e[2]`), fill = "#6a3d9a", alpha = .4, bins = 200) +
labs(x = expression(σ[year]), y = "", title = "POSY") + theme_classic()
survvar <- grid.arrange(AGPE_s, ELRI_s, ELVI_s, FESU_s, LOAR_s, POAL_s, POSY_s, ncol= 1)
titlesurv <- annotate_figure(survvar, top = "Survival", left = "")
flwrvar <- grid.arrange(AGPE_f, ELRI_f, ELVI_f, FESU_f, LOAR_f, POAL_f, POSY_f, ncol= 1)
titleflw <- annotate_figure(flwrvar, top = "Flowering", left = "")
growvar <- grid.arrange(AGPE_g, ELRI_g, ELVI_g, FESU_g, LOAR_g, POAL_g, POSY_g, ncol= 1)
titlegrow <- annotate_figure(growvar, top = "Growth", left = "")
var <- grid.arrange(titlesurv, titleflw, titlegrow, ncol = 3)
titlevar <- annotate_figure(var, top = "Interannual Variance", left = "Posterior Density")
titlevar
# boxplots of variance estimates
AGPE_melt_surv <- post_survAGPE %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "AGPE") %>%
melt()
sumAGPE_melt_surv <- rownames_to_column(as.data.frame(sum_survAGPE$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "AGPE")
ELVI_melt_surv <- post_survELVI %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "ELVI") %>%
melt()
sumELVI_melt_surv <- rownames_to_column(as.data.frame(sum_survELVI$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "ELVI")
ELRI_melt_surv <- post_survELRI %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "ELRI") %>%
melt()
sumELRI_melt_surv <- rownames_to_column(as.data.frame(sum_survELRI$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "ELRI")
FESU_melt_surv <- post_survFESU %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "FESU") %>%
melt()
sumFESU_melt_surv <- rownames_to_column(as.data.frame(sum_survFESU$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "FESU")
LOAR_melt_surv <- post_survLOAR %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "LOAR") %>%
melt()
sumLOAR_melt_surv <- rownames_to_column(as.data.frame(sum_survLOAR$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "LOAR")
POAL_melt_surv <- post_survPOAL %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "POAL") %>%
melt()
sumPOAL_melt_surv <- rownames_to_column(as.data.frame(sum_survPOAL$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "POAL")
POSY_melt_surv <- post_survPOSY %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "POSY") %>%
melt()
sumPOSY_melt_surv <- rownames_to_column(as.data.frame(sum_survPOSY$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "POSY")
surv_species <- AGPE_melt_surv %>%
merge(ELRI_melt_surv, by = c("species", "value", "variable"), all = TRUE) %>%
merge(ELVI_melt_surv, by = c("species", "value", "variable"), all = TRUE) %>%
merge(FESU_melt_surv, by = c("species", "value", "variable"), all = TRUE) %>%
merge(LOAR_melt_surv, by = c("species", "value", "variable"), all = TRUE) %>%
merge(POAL_melt_surv, by = c("species", "value", "variable"), all = TRUE) %>%
merge(POSY_melt_surv, by = c("species", "value", "variable"), all = TRUE)
sum_surv_species <- sumAGPE_melt_surv %>%
merge(sumELRI_melt_surv, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumELVI_melt_surv, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumFESU_melt_surv, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumLOAR_melt_surv, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumPOAL_melt_surv, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumPOSY_melt_surv, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE)
AGPE_melt_flw <- post_flwAGPE %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "AGPE") %>%
melt()
sumAGPE_melt_flw <- rownames_to_column(as.data.frame(sum_flwAGPE$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "AGPE")
ELVI_melt_flw <- post_flwELVI %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "ELVI") %>%
melt()
sumELVI_melt_flw <- rownames_to_column(as.data.frame(sum_flwELVI$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "ELVI")
ELRI_melt_flw <- post_flwELRI %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "ELRI") %>%
melt()
sumELRI_melt_flw <- rownames_to_column(as.data.frame(sum_flwELRI$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "ELRI")
FESU_melt_flw <- post_flwFESU %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "FESU") %>%
melt()
sumFESU_melt_flw <- rownames_to_column(as.data.frame(sum_flwFESU$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "FESU")
LOAR_melt_flw <- post_flwLOAR %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "LOAR") %>%
melt()
sumLOAR_melt_flw <- rownames_to_column(as.data.frame(sum_flwLOAR$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "LOAR")
POAL_melt_flw <- post_flwPOAL %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "POAL") %>%
melt()
sumPOAL_melt_flw <- rownames_to_column(as.data.frame(sum_flwPOAL$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "POAL")
POSY_melt_flw <- post_flwPOSY %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "POSY") %>%
melt()
sumPOSY_melt_flw <- rownames_to_column(as.data.frame(sum_flwPOSY$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "POSY")
flw_species <- AGPE_melt_flw %>%
merge(ELRI_melt_flw, by = c("species", "value", "variable"), all = TRUE) %>%
merge(ELVI_melt_flw, by = c("species", "value", "variable"), all = TRUE) %>%
merge(FESU_melt_flw, by = c("species", "value", "variable"), all = TRUE) %>%
merge(LOAR_melt_flw, by = c("species", "value", "variable"), all = TRUE) %>%
merge(POAL_melt_flw, by = c("species", "value", "variable"), all = TRUE) %>%
merge(POSY_melt_flw, by = c("species", "value", "variable"), all = TRUE)
sum_flw_species <- sumAGPE_melt_flw %>%
merge(sumELRI_melt_flw, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumELVI_melt_flw, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumFESU_melt_flw, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumLOAR_melt_flw, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumPOAL_melt_flw, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumPOSY_melt_flw, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE)
AGPE_melt_grow <- post_growAGPE %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "AGPE") %>%
melt()
sumAGPE_melt_grow <- rownames_to_column(as.data.frame(sum_growAGPE$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "AGPE")
ELVI_melt_grow <- post_growELVI %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "ELVI") %>%
melt()
sumELVI_melt_grow <- rownames_to_column(as.data.frame(sum_growELVI$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "ELVI")
ELRI_melt_grow <- post_growELRI %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "ELRI") %>%
melt()
sumELRI_melt_grow <- rownames_to_column(as.data.frame(sum_growELRI$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "ELRI")
FESU_melt_grow <- post_growFESU %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "FESU") %>%
melt()
sumFESU_melt_grow <- rownames_to_column(as.data.frame(sum_growFESU$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "FESU")
LOAR_melt_grow <- post_growLOAR %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "LOAR") %>%
melt()
sumLOAR_melt_grow <- rownames_to_column(as.data.frame(sum_growLOAR$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "LOAR")
POAL_melt_grow <- post_growPOAL %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "POAL") %>%
melt()
sumPOAL_melt_grow <- rownames_to_column(as.data.frame(sum_growPOAL$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "POAL")
POSY_melt_grow <- post_growPOSY %>%
select("sigma_e[1]", "sigma_e[2]") %>%
mutate(species = "POSY") %>%
melt()
sumPOSY_melt_grow <- rownames_to_column(as.data.frame(sum_growPOSY$summary)) %>%
select("rowname", "mean", "2.5%", "97.5%") %>%
filter(rowname == "sigma_e[1]" | rowname == "sigma_e[2]") %>%
mutate(species = "POSY")
grow_species <- AGPE_melt_grow %>%
merge(ELRI_melt_grow, by = c("species", "value", "variable"), all = TRUE) %>%
merge(ELVI_melt_grow, by = c("species", "value", "variable"), all = TRUE) %>%
merge(FESU_melt_grow, by = c("species", "value", "variable"), all = TRUE) %>%
merge(LOAR_melt_grow, by = c("species", "value", "variable"), all = TRUE) %>%
merge(POAL_melt_grow, by = c("species", "value", "variable"), all = TRUE) %>%
merge(POSY_melt_grow, by = c("species", "value", "variable"), all = TRUE)
sum_grow_species <- sumAGPE_melt_grow %>%
merge(sumELRI_melt_grow, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumELVI_melt_grow, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumFESU_melt_grow, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumLOAR_melt_grow, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumPOAL_melt_grow, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE) %>%
merge(sumPOSY_melt_grow, by = c("rowname", "species", "mean", "2.5%", "97.5%"), all = TRUE)
# options
ggplot() +
geom_boxplot(data = surv_species, aes(y = value, x = variable, fill = variable)) +
facet_wrap(~species, nrow = 2) +
labs(x = "species", y = expression(σ[year]), title = "POAL Survival Variance") + theme_classic() +
scale_fill_manual(values = colors2)
ggplot() +
geom_boxplot(data = surv_species, aes(y = value, fill = variable)) +
facet_grid(row = vars(species), col = vars(variable)) +
labs(x = "species", y = expression(σ[year]), title = "POAL Survival Variance") + theme_classic() +
scale_fill_manual(values = colors2)
ggplot(data = sum_surv_species) +
geom_point(aes(y = mean, x = species, color = rowname)) +
geom_errorbar(aes(x = species, ymin = `2.5%`, ymax = `97.5%`, color = rowname), width = .5) +
labs(x = "species", y = expression(σ[year]), title = "POAL Survival Variance") + theme_classic() +
scale_color_manual(values = colors2)
# means with 90%CI
surv <- ggplot(data = sum_surv_species) +
geom_point(aes(y = mean, x = species, color = rowname), position = position_dodge( width = .7), size = 2) +
geom_errorbar(aes(x = species, ymin = `2.5%`, ymax = `97.5%`, color = rowname), width = 0, position = position_dodge(width = .7), lwd = 1, alpha = .5) +
labs(x = "", y = expression(σ[year]), title = "Survival") + theme_classic() +
scale_color_manual(values = colors2)
flw <- ggplot(data = sum_flw_species) +
geom_point(aes(y = mean, x = species, color = rowname), position = position_dodge( width = .7), size = 2) +
geom_errorbar(aes(x = species, ymin = `2.5%`, ymax = `97.5%`, color = rowname), width = 0, position = position_dodge(width = .7), lwd = 1, alpha = .5) +
labs(x = "", y = expression(σ[year]), title = "Flowering") + theme_classic() +
scale_color_manual(values = colors2)
grow <- ggplot(data = sum_grow_species) +
geom_point(aes(y = mean, x = species, color = rowname), position = position_dodge( width = .7), size = 2) +
geom_errorbar(aes(x = species, ymin = `2.5%`, ymax = `97.5%`, color = rowname), width = 0, position = position_dodge(width = .7), lwd = 1, alpha = .5) +
labs(x = "", y = expression(σ[year]), title = "Growth") + theme_classic() +
scale_color_manual(values = colors2)
mCIplot <- grid.arrange(surv, grow, nrow= 2)
titlemCIplot <- annotate_figure(mCIplot, top = "", bottom = "Species", left = "")
titlemCIplot
###### plots of Survival model fits#######
######## POAL model fits
# actual data points for probability of survival
POALsurv_bin0 <- POAL_data %>%
select(surv_t1, logsize_t, year_t, endo) %>%
mutate(endo = recode(endo, "0" = 0, "1" = 1, minus = 0, plus = 1)) %>%
filter(endo == 0) %>%
mutate(Year = as.factor(as.integer(as.factor(year_t)))) %>%
mutate(size_bin = cut_interval(logsize_t, n = 15)) %>%
group_by(size_bin, Year) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
POALsurv_bin1 <- POAL_data %>%
select(surv_t1, logsize_t, year_t, endo) %>%
mutate(endo = recode(endo, "0" = 0, "1" = 1, minus = 0, plus = 1)) %>%
filter(endo == 1) %>%
mutate(Year = as.factor(as.integer(as.factor(year_t)))) %>%
mutate(size_bin = cut_interval(logsize_t, n = 15)) %>%
group_by(size_bin, Year) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
POALsurv_bin <- POAL_data %>%
select(surv_t1, logsize_t, year_t, endo) %>%
mutate(endo = recode(endo, "0" = 0, "1" = 1, minus = 0, plus = 1)) %>%
mutate(Year = as.factor(as.integer(as.factor(year_t)))) %>%
mutate(size_bin = cut_interval(logsize_t, n = 15)) %>%
group_by(size_bin, Year, endo) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
POALsurv_binmean <- POAL_data %>%
select(surv_t1, logsize_t, endo) %>%
mutate(endo = recode(endo, "0" = "E-", "1" = "E+", minus = "E-", plus = "E+")) %>%
mutate(size_bin = cut_interval(logsize_t, n = 12)) %>%
group_by(size_bin, endo) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
# fitted model line
nvalues <- length(POAL_data$logsize_t)
xdummy <- seq(min(POAL_data$logsize_t), max(POAL_data$logsize_t), length.out = nvalues)
# E-
POAL_ydummy_eminus <- as.vector(invlogit(sum_survPOAL$summary["beta[1]","mean"] + (sum_survPOAL$summary["beta[2]","mean"])*xdummy + (sum_survPOAL$summary["beta[3]","mean"])*0 + (sum_survPOAL$summary["beta[4]","mean"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","mean"])*xdummy*0))
POAL_ydummy_eminus2.5 <- as.vector(invlogit(sum_survPOAL$summary["beta[1]","2.5%"] + (sum_survPOAL$summary["beta[2]","2.5%"])*xdummy + (sum_survPOAL$summary["beta[3]","2.5%"])*0 + (sum_survPOAL$summary["beta[4]","2.5%"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","2.5%"])*xdummy*0))
POAL_ydummy_eminus97.5 <- as.vector(invlogit(sum_survPOAL$summary["beta[1]","97.5%"] + (sum_survPOAL$summary["beta[2]","97.5%"])*xdummy + (sum_survPOAL$summary["beta[3]","97.5%"])*0 + (sum_survPOAL$summary["beta[4]","97.5%"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","97.5%"])*xdummy*0))
POAL_ydummy_eminus_semin <- as.vector(invlogit((sum_survPOAL$summary["beta[1]","mean"]-sum_survPOAL$summary["beta[1]","se_mean"]) + (sum_survPOAL$summary["beta[2]","mean"]-sum_survPOAL$summary["beta[2]","se_mean"])*xdummy + (sum_survPOAL$summary["beta[3]","mean"]-sum_survPOAL$summary["beta[3]","se_mean"])*0 + (sum_survPOAL$summary["beta[4]","mean"]-sum_survPOAL$summary["beta[4]","se_mean"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","mean"]-sum_survPOAL$summary["beta[5]","se_mean"])*xdummy*0))
POAL_ydummy_eminus_semax <- as.vector(invlogit((sum_survPOAL$summary["beta[1]","mean"]+sum_survPOAL$summary["beta[1]","se_mean"]) + (sum_survPOAL$summary["beta[2]","mean"]+sum_survPOAL$summary["beta[2]","se_mean"])*xdummy + (sum_survPOAL$summary["beta[3]","mean"]+sum_survPOAL$summary["beta[3]","se_mean"])*0 + (sum_survPOAL$summary["beta[4]","mean"]+sum_survPOAL$summary["beta[4]","se_mean"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","mean"]+sum_survPOAL$summary["beta[5]","se_mean"])*xdummy*0))
POAL_ydummy_eminus_y1 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,1]')))
POAL_ydummy_eminus_y2 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,2]')))
POAL_ydummy_eminus_y3 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,3]')))
POAL_ydummy_eminus_y4 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,4]')))
POAL_ydummy_eminus_y5 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,5]')))
POAL_ydummy_eminus_y6 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,6]')))
POAL_ydummy_eminus_y7 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,7]')))
POAL_ydummy_eminus_y8 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,8]')))
POAL_ydummy_eminus_y9 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,9]')))
POAL_ydummy_eminus_y10 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,10]')))
POAL_ydummy_eminus_y11 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*0 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*0
+ mean(post_survPOAL$'tau_year[1,11]')))
POALfitsminus0 <- as.data.frame(cbind(xdummy, POAL_ydummy_eminus, POAL_ydummy_eminus2.5, POAL_ydummy_eminus97.5))
POALfitsminus1 <- as.data.frame(cbind(xdummy,POAL_ydummy_eminus_y1,POAL_ydummy_eminus_y2, POAL_ydummy_eminus_y3, POAL_ydummy_eminus_y4, POAL_ydummy_eminus_y5, POAL_ydummy_eminus_y6, POAL_ydummy_eminus_y7, POAL_ydummy_eminus_y8, POAL_ydummy_eminus_y9, POAL_ydummy_eminus_y10, POAL_ydummy_eminus_y11))
POALfitsminus <- melt(POALfitsminus1, id = "xdummy",
measure.vars = c("POAL_ydummy_eminus_y1","POAL_ydummy_eminus_y2",
"POAL_ydummy_eminus_y3", "POAL_ydummy_eminus_y4",
"POAL_ydummy_eminus_y5", "POAL_ydummy_eminus_y6",
"POAL_ydummy_eminus_y7", "POAL_ydummy_eminus_y8",
"POAL_ydummy_eminus_y9", "POAL_ydummy_eminus_y10",
"POAL_ydummy_eminus_y11")) %>%
mutate(Year = recode(variable, "POAL_ydummy_eminus_y1" = '1',"POAL_ydummy_eminus_y2" = "2",
"POAL_ydummy_eminus_y3" = "3", "POAL_ydummy_eminus_y4" = "4",
"POAL_ydummy_eminus_y5" = "5", "POAL_ydummy_eminus_y6" = "6",
"POAL_ydummy_eminus_y7" = "7", "POAL_ydummy_eminus_y8" = "8",
"POAL_ydummy_eminus_y9" = "9", "POAL_ydummy_eminus_y10" = "10",
"POAL_ydummy_eminus_y11" = "11") )
# E+
POAL_ydummy_eplus <- as.vector(invlogit(sum_survPOAL$summary["beta[1]","mean"] + (sum_survPOAL$summary["beta[2]","mean"])*xdummy + (sum_survPOAL$summary["beta[3]","mean"])*1 + (sum_survPOAL$summary["beta[4]","mean"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","mean"])*xdummy*1))
POAL_ydummy_eplus2.5 <- as.vector(invlogit(sum_survPOAL$summary["beta[1]","2.5%"] + (sum_survPOAL$summary["beta[2]","2.5%"])*xdummy + (sum_survPOAL$summary["beta[3]","2.5%"])*1 + (sum_survPOAL$summary["beta[4]","2.5%"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","2.5%"])*xdummy*1))
POAL_ydummy_eplus97.5 <- as.vector(invlogit(sum_survPOAL$summary["beta[1]","97.5%"] + (sum_survPOAL$summary["beta[2]","97.5%"])*xdummy + (sum_survPOAL$summary["beta[3]","97.5%"])*1 + (sum_survPOAL$summary["beta[4]","97.5%"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","97.5%"])*xdummy*1))
POAL_ydummy_eplus_semin <- as.vector(invlogit((sum_survPOAL$summary["beta[1]","mean"]-sum_survPOAL$summary["beta[1]","se_mean"]) + (sum_survPOAL$summary["beta[2]","mean"]-sum_survPOAL$summary["beta[2]","se_mean"])*xdummy + (sum_survPOAL$summary["beta[3]","mean"]-sum_survPOAL$summary["beta[3]","se_mean"])*1 + (sum_survPOAL$summary["beta[4]","mean"]-sum_survPOAL$summary["beta[4]","se_mean"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","mean"]-sum_survPOAL$summary["beta[5]","se_mean"])*xdummy*1))
POAL_ydummy_eplus_semax <- as.vector(invlogit((sum_survPOAL$summary["beta[1]","mean"]+sum_survPOAL$summary["beta[1]","se_mean"]) + (sum_survPOAL$summary["beta[2]","mean"]+sum_survPOAL$summary["beta[2]","se_mean"])*xdummy + (sum_survPOAL$summary["beta[3]","mean"]+sum_survPOAL$summary["beta[3]","se_mean"])*1 + (sum_survPOAL$summary["beta[4]","mean"]+sum_survPOAL$summary["beta[4]","se_mean"])*mean(POAL_data$origin_01)
+ (sum_survPOAL$summary["beta[5]","mean"]+sum_survPOAL$summary["beta[5]","se_mean"])*xdummy*1))
POAL_ydummy_eplus_y1 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,1]')))
POAL_ydummy_eplus_y2 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,2')))
POAL_ydummy_eplus_y3 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,3]')))
POAL_ydummy_eplus_y4 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,4]')))
POAL_ydummy_eplus_y5 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,5]')))
POAL_ydummy_eplus_y6 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,6]')))
POAL_ydummy_eplus_y7 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,7]')))
POAL_ydummy_eplus_y8 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,8]')))
POAL_ydummy_eplus_y9 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,9]')))
POAL_ydummy_eplus_y10 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,10]')))
POAL_ydummy_eplus_y11 <- as.vector(invlogit(mean(post_survPOAL$'beta[1]') + mean(post_survPOAL$'beta[2]')*xdummy + mean(post_survPOAL$'beta[3]')*1 + mean(post_survPOAL$'beta[4]')*mean(POAL_data$origin_01)
+ mean(post_survPOAL$'beta[5]')*xdummy*1
+ mean(post_survPOAL$'tau_year[2,11]')))
POAL_fitsplus0 <- as.data.frame(cbind(xdummy, POAL_ydummy_eplus, POAL_ydummy_eplus2.5, POAL_ydummy_eplus97.5))
POALfitsplus1 <- as.data.frame(cbind(xdummy,POAL_ydummy_eplus_y1,POAL_ydummy_eplus_y2, POAL_ydummy_eplus_y3, POAL_ydummy_eplus_y4, POAL_ydummy_eplus_y5, POAL_ydummy_eplus_y6, POAL_ydummy_eplus_y7, POAL_ydummy_eplus_y8, POAL_ydummy_eplus_y9, POAL_ydummy_eplus_y10, POAL_ydummy_eplus_y11))
POALfitsplus <- melt(POALfitsplus1, id = "xdummy",
measure.vars = c("POAL_ydummy_eplus_y1","POAL_ydummy_eplus_y2",
"POAL_ydummy_eplus_y3", "POAL_ydummy_eplus_y4",
"POAL_ydummy_eplus_y5", "POAL_ydummy_eplus_y6",
"POAL_ydummy_eplus_y7", "POAL_ydummy_eplus_y8",
"POAL_ydummy_eplus_y9", "POAL_ydummy_eplus_y10",
"POAL_ydummy_eplus_y11")) %>%
mutate(Year = recode(variable, "POAL_ydummy_eplus_y1" = '1',"POAL_ydummy_eplus_y2" = "2",
"POAL_ydummy_eplus_y3" = "3", "POAL_ydummy_eplus_y4" = "4",
"POAL_ydummy_eplus_y5" = "5", "POAL_ydummy_eplus_y6" = "6",
"POAL_ydummy_eplus_y7" = "7", "POAL_ydummy_eplus_y8" = "8",
"POAL_ydummy_eplus_y9" = "9", "POAL_ydummy_eplus_y10" = "10",
"POAL_ydummy_eplus_y11" = "11") )
# E+ by year
POALEplusbyyear <- ggplot(data = POALfitsplus) +
geom_line(aes(x = xdummy, y = value, color = Year), lwd = .8) +
geom_point(data = POALsurv_bin1, aes(x = mean_size, y = mean_surv, color = Year, lwd = samplesize)) +
ggtitle("E+ ") + xlab("log(size_t)") + ylab("Prob. of Survival") +
scale_color_brewer(palette = "Paired") + guides(lwd = "none") +theme_classic()
ggplot(data = POALfitsplus) +
geom_line(aes(x = xdummy, y = value, color = Year), lwd = .8) +
ggtitle("E+") + xlab("log(size_t)") + ylab("Prob. of Survival")
scale_color_manual(values=yearcolors)
# E- by year
POALEminusbyyear <- ggplot(data = POALfitsminus) +
geom_line(aes(x = xdummy, y = value, color = Year), lwd = .8) +
geom_point(data = POALsurv_bin0, aes(x = mean_size, y = mean_surv, color = Year, lwd = samplesize)) +
ggtitle("E-") + xlab("log(size_t)") + ylab("Prob. of Survival") +
scale_color_brewer(palette = "Paired") + guides(lwd = "none")+theme_classic()
ggplot(data = POALfitsminus) +
geom_line(aes(x = xdummy, y = value, color = Year), lwd = .8) +
ggtitle("POAL E- Survival Probability") + xlab("log(size_t)") + ylab("Prob. of Survival")
scale_color_manual(values=yearcolors)
# effect on mean
POALfits <- as.data.frame(cbind(xdummy, POAL_ydummy_eplus, POAL_ydummy_eplus_semin, POAL_ydummy_eplus_semax, POAL_ydummy_eplus2.5, POAL_ydummy_eplus97.5, POAL_ydummy_eminus, POAL_ydummy_eminus_semin, POAL_ydummy_eminus_semax, POAL_ydummy_eminus2.5, POAL_ydummy_eminus97.5))
# without point
ggplot(data = POALfits) +
geom_ribbon(aes(x = xdummy, ymin = POAL_ydummy_eplus2.5, ymax = POAL_ydummy_eplus97.5), fill = "#6a3d9a", alpha = .2) +
geom_ribbon(aes(x = xdummy, ymin = POAL_ydummy_eminus2.5, ymax = POAL_ydummy_eminus97.5), fill = "#ff7f00", alpha = .2)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus), color = "#ff7f00") +
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus), color = "#6a3d9a")+
ggtitle("POAL Survival Probability") + xlab("log(size_t)") + ylab("Prob. of Survival") +theme_classic()
ggplot(data = POALfits) +
geom_ribbon(aes(x = xdummy, ymin = POAL_ydummy_eplus_semin, ymax = POAL_ydummy_eplus_semax), fill = "#6a3d9a", alpha = .2) +
geom_ribbon(aes(x = xdummy, ymin = POAL_ydummy_eminus_semin, ymax = POAL_ydummy_eminus_semax), fill = "#ff7f00", alpha = .2)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus), color = "#6a3d9a")+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus), color = "#ff7f00") +
ggtitle("POAL Survival Probability") + xlab("log(size_t)") + ylab("Prob. of Survival") +theme_classic()
# with points
POALmean <- ggplot(data = POALfits) +
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus), color = "#6a3d9a")+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus), color = "#ff7f00") +
geom_point(data = POALsurv_binmean, aes(x = mean_size, y = mean_surv, color = endo, lwd = samplesize)) +
ggtitle("Mean Survival Probability") + xlab("log(size_t)") + ylab("Prob. of Survival") +
scale_color_manual(values = colors2) + theme_classic() + guides(lwd = "none")
ggplot(data = POALfits) +
geom_ribbon(aes(x = xdummy, ymin = POAL_ydummy_eplus_semin, ymax = POAL_ydummy_eplus_semax), fill = "#6a3d9a", alpha = .2) +
geom_ribbon(aes(x = xdummy, ymin = POAL_ydummy_eminus_semin, ymax = POAL_ydummy_eminus_semax), fill = "#ff7f00", alpha = .2)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus), color = "#6a3d9a")+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus), color = "#ff7f00") +
geom_point(data = POALsurv_binmean, aes(x = mean_size, y = mean_surv, color = endo, lwd = samplesize)) +
ggtitle("POAL Survival Probability") + xlab("log(size_t)") + ylab("Prob. of Survival") +
scale_color_manual(values = colors2) +theme_classic() + guides(lwd = "none")
POALsurv4panel <- grid.arrange(POALmean, POAL_s,POALEplusbyyear,POALEminusbyyear, ncol= 4)
titlePOALsurv4panel <- annotate_figure(POALsurv4panel, top = "Poa alsodes")
titlePOALsurv4panel
# split up plus and minus with datapoints and add mean line
POAL_minus <- ggplot(data = POALfitsminus1)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y1), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y2), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y3), col = "#ff7f00", alpha = .5) +
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y4), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y5), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y6), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y7), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y8), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y9), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y10), col = "#ff7f00", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eminus_y11), col = "#ff7f00", alpha = .5)+
geom_line(data = POALfitsminus0, aes(x = xdummy, y = POAL_ydummy_eminus), col = "#ff7f00", lwd = 1) +
geom_point(data = POALsurv_bin0, aes(x = mean_size, y = mean_surv), color ="#ff7f00") +
theme_classic()+
labs(title = "E-", x = "log(size_t)", y = "")
POAL_plus <- ggplot(data = POALfitsplus1)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y1), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y2), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y3), col = "#6a3d9a", alpha = .5) +
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y4), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y5), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y6), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y7), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y8), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y9), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y10), col = "#6a3d9a", alpha = .5)+
geom_line(aes(x = xdummy, y = POAL_ydummy_eplus_y11), col = "#6a3d9a", alpha = .5)+
geom_line(data = POAL_fitsplus0, aes(x = xdummy, y = POAL_ydummy_eplus), col = "#6a3d9a", lwd = 1) +
geom_point(data = POALsurv_bin1, aes(x = mean_size, y = mean_surv), color ="#6a3d9a") +
theme_classic()+
labs(title = "E+", x = "log(size_t)", y = "")
POALsurv <- grid.arrange(POAL_plus, POAL_minus, ncol= 2)
titlePOALsurv <- annotate_figure(POALsurv, top = "POAL", left = "Probability of Survival")
titlePOALsurv
# with data points
ggplot(data = POALfitsminus) +
geom_line(aes(x = xdummy, y = value, color = Year), lwd = .7) +
geom_point(data = POALsurv_bin0, aes(x = mean_size, y = mean_surv, color = Year)) +
labs(title = "POAL E- Survival Probability", x = "log(size_t)", y = "Prob. of Survival") +
scale_color_manual(values=yearcolors)
# without datapoints
ggplot(data = POALfitsminus) +
geom_line(aes(x = xdummy, y = value, color = Year), lwd = .7) +
labs(title = "POAL E- Survival Probability", x = "log(size_t)", y = "Prob. of Survival") + scale_color_manual(values=yearcolors)
######## POSY model fits
# actual data points for probability of survival
POSYsurv_bin0 <- POSY_data %>%
select(surv_t1, logsize_t, year_t, endo) %>%
mutate(endo = recode(endo, "0" = 0, "1" = 1, minus = 0, plus = 1)) %>%
filter(endo == 0) %>%
mutate(Year = as.factor(as.integer(as.factor(year_t)))) %>%
mutate(size_bin = cut_interval(logsize_t, n = 15)) %>%
group_by(size_bin, Year) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
POSYsurv_bin1 <- POSY_data %>%
select(surv_t1, logsize_t, year_t, endo) %>%
mutate(endo = recode(endo, "0" = 0, "1" = 1, minus = 0, plus = 1)) %>%
filter(endo == 1) %>%
mutate(Year = as.factor(as.integer(as.factor(year_t)))) %>%
mutate(size_bin = cut_interval(logsize_t, n = 15)) %>%
group_by(size_bin, Year) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
POSYsurv_bin <- POSY_data %>%
select(surv_t1, logsize_t, year_t, endo) %>%
mutate(endo = recode(endo, "0" = 0, "1" = 1, minus = 0, plus = 1)) %>%
mutate(Year = as.factor(as.integer(as.factor(year_t)))) %>%
mutate(size_bin = cut_interval(logsize_t, n = 15)) %>%
group_by(size_bin, Year, endo) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
POSYsurv_binmean <- POSY_data %>%
select(surv_t1, logsize_t, endo) %>%
mutate(endo = recode(endo, "0" = "E-", "1" = "E+", minus = "E-", plus = "E+")) %>%
mutate(size_bin = cut_interval(logsize_t, n = 12)) %>%
group_by(size_bin, endo) %>%
summarise(mean_size = mean((logsize_t),na.rm=T),
mean_surv = mean(surv_t1,na.rm=T),
samplesize = n())
# fitted model line
nvalues <- length(POSY_data$logsize_t)
xdummy <- seq(min(POSY_data$logsize_t), max(POSY_data$logsize_t), length.out = nvalues)
# E-
POSY_ydummy_eminus <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0))
POSY_ydummy_eminus2.5 <- as.vector(invlogit(sum_survPOSY$summary["beta[1]","2.5%"] + (sum_survPOSY$summary["beta[2]","2.5%"])*xdummy + (sum_survPOSY$summary["beta[3]","2.5%"])*0 + (sum_survPOSY$summary["beta[4]","2.5%"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","2.5%"])*xdummy*0))
POSY_ydummy_eminus97.5 <- as.vector(invlogit(sum_survPOSY$summary["beta[1]","97.5%"] + (sum_survPOSY$summary["beta[2]","97.5%"])*xdummy + (sum_survPOSY$summary["beta[3]","97.5%"])*0 + (sum_survPOSY$summary["beta[4]","97.5%"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","97.5%"])*xdummy*0))
POSY_ydummy_eminus_semin <- as.vector(invlogit((sum_survPOSY$summary["beta[1]","mean"]-sum_survPOSY$summary["beta[1]","se_mean"]) + (sum_survPOSY$summary["beta[2]","mean"]-sum_survPOSY$summary["beta[2]","se_mean"])*xdummy + (sum_survPOSY$summary["beta[3]","mean"]-sum_survPOSY$summary["beta[3]","se_mean"])*0 + (sum_survPOSY$summary["beta[4]","mean"]-sum_survPOSY$summary["beta[4]","se_mean"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","mean"]-sum_survPOSY$summary["beta[5]","se_mean"])*xdummy*0))
POSY_ydummy_eminus_semax <- as.vector(invlogit((sum_survPOSY$summary["beta[1]","mean"]+sum_survPOSY$summary["beta[1]","se_mean"]) + (sum_survPOSY$summary["beta[2]","mean"]+sum_survPOSY$summary["beta[2]","se_mean"])*xdummy + (sum_survPOSY$summary["beta[3]","mean"]+sum_survPOSY$summary["beta[3]","se_mean"])*0 + (sum_survPOSY$summary["beta[4]","mean"]+sum_survPOSY$summary["beta[4]","se_mean"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","mean"]+sum_survPOSY$summary["beta[5]","se_mean"])*xdummy*0))
POSY_ydummy_eminus_y1 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,1]')))
POSY_ydummy_eminus_y2 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,2]')))
POSY_ydummy_eminus_y3 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,3]')))
POSY_ydummy_eminus_y4 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,4]')))
POSY_ydummy_eminus_y5 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,5]')))
POSY_ydummy_eminus_y6 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,6]')))
POSY_ydummy_eminus_y7 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,7]')))
POSY_ydummy_eminus_y8 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,8]')))
POSY_ydummy_eminus_y9 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,9]')))
POSY_ydummy_eminus_y10 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,10]')))
POSY_ydummy_eminus_y11 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*0 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*0
+ mean(post_survPOSY$'tau_year[1,11]')))
POSYfitsminus0 <- as.data.frame(cbind(xdummy, POSY_ydummy_eminus))
POSYfitsminus1 <- as.data.frame(cbind(xdummy,POSY_ydummy_eminus_y1,POSY_ydummy_eminus_y2, POSY_ydummy_eminus_y3, POSY_ydummy_eminus_y4, POSY_ydummy_eminus_y5, POSY_ydummy_eminus_y6, POSY_ydummy_eminus_y7, POSY_ydummy_eminus_y8, POSY_ydummy_eminus_y9, POSY_ydummy_eminus_y10, POSY_ydummy_eminus_y11))
POSYfitsminus <- melt(POSYfitsminus1, id = "xdummy",
measure.vars = c("POSY_ydummy_eminus_y1","POSY_ydummy_eminus_y2",
"POSY_ydummy_eminus_y3", "POSY_ydummy_eminus_y4",
"POSY_ydummy_eminus_y5", "POSY_ydummy_eminus_y6",
"POSY_ydummy_eminus_y7", "POSY_ydummy_eminus_y8",
"POSY_ydummy_eminus_y9", "POSY_ydummy_eminus_y10",
"POSY_ydummy_eminus_y11")) %>%
mutate(Year = recode(variable, "POSY_ydummy_eminus_y1" = '1',"POSY_ydummy_eminus_y2" = "2",
"POSY_ydummy_eminus_y3" = "3", "POSY_ydummy_eminus_y4" = "4",
"POSY_ydummy_eminus_y5" = "5", "POSY_ydummy_eminus_y6" = "6",
"POSY_ydummy_eminus_y7" = "7", "POSY_ydummy_eminus_y8" = "8",
"POSY_ydummy_eminus_y9" = "9", "POSY_ydummy_eminus_y10" = "10",
"POSY_ydummy_eminus_y11" = "11") )
# E+
POSY_ydummy_eplus <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1))
POSY_ydummy_eplus2.5 <- as.vector(invlogit(sum_survPOSY$summary["beta[1]","2.5%"] + (sum_survPOSY$summary["beta[2]","2.5%"])*xdummy + (sum_survPOSY$summary["beta[3]","2.5%"])*1 + (sum_survPOSY$summary["beta[4]","2.5%"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","2.5%"])*xdummy*1))
POSY_ydummy_eplus97.5 <- as.vector(invlogit(sum_survPOSY$summary["beta[1]","97.5%"] + (sum_survPOSY$summary["beta[2]","97.5%"])*xdummy + (sum_survPOSY$summary["beta[3]","97.5%"])*1 + (sum_survPOSY$summary["beta[4]","97.5%"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","97.5%"])*xdummy*1))
POSY_ydummy_eplus_semin <- as.vector(invlogit((sum_survPOSY$summary["beta[1]","mean"]-sum_survPOSY$summary["beta[1]","se_mean"]) + (sum_survPOSY$summary["beta[2]","mean"]-sum_survPOSY$summary["beta[2]","se_mean"])*xdummy + (sum_survPOSY$summary["beta[3]","mean"]-sum_survPOSY$summary["beta[3]","se_mean"])*1 + (sum_survPOSY$summary["beta[4]","mean"]-sum_survPOSY$summary["beta[4]","se_mean"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","mean"]-sum_survPOSY$summary["beta[5]","se_mean"])*xdummy*1))
POSY_ydummy_eplus_semax <- as.vector(invlogit((sum_survPOSY$summary["beta[1]","mean"]+sum_survPOSY$summary["beta[1]","se_mean"]) + (sum_survPOSY$summary["beta[2]","mean"]+sum_survPOSY$summary["beta[2]","se_mean"])*xdummy + (sum_survPOSY$summary["beta[3]","mean"]+sum_survPOSY$summary["beta[3]","se_mean"])*1 + (sum_survPOSY$summary["beta[4]","mean"]+sum_survPOSY$summary["beta[4]","se_mean"])*mean(POSY_data$origin_01)
+ (sum_survPOSY$summary["beta[5]","mean"]+sum_survPOSY$summary["beta[5]","se_mean"])*xdummy*1))
POSY_ydummy_eplus_y1 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,1]')))
POSY_ydummy_eplus_y2 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,2')))
POSY_ydummy_eplus_y3 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,3]')))
POSY_ydummy_eplus_y4 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,4]')))
POSY_ydummy_eplus_y5 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,5]')))
POSY_ydummy_eplus_y6 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,6]')))
POSY_ydummy_eplus_y7 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,7]')))
POSY_ydummy_eplus_y8 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,8]')))
POSY_ydummy_eplus_y9 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,9]')))
POSY_ydummy_eplus_y10 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,10]')))
POSY_ydummy_eplus_y11 <- as.vector(invlogit(mean(post_survPOSY$'beta[1]') + mean(post_survPOSY$'beta[2]')*xdummy + mean(post_survPOSY$'beta[3]')*1 + mean(post_survPOSY$'beta[4]')*mean(POSY_data$origin_01)
+ mean(post_survPOSY$'beta[5]')*xdummy*1
+ mean(post_survPOSY$'tau_year[2,11]')))