-
Notifications
You must be signed in to change notification settings - Fork 1
/
2_data_analysis_figures.R
1716 lines (1414 loc) · 91.4 KB
/
2_data_analysis_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
# Fig. 1 maps ---------
# ArcGIS
# Fig 2. schematic ----------
dev.off()
rm(list=ls())
library(ggplot2)
library(cowplot)
library(stringr)
library(reshape2)
library(ggthemes)
library(RColorBrewer)
setwd("...6_soil_maxes")
cvs_full_data<-read.csv("...cvs_full_data_5cm.csv")
names(cvs_full_data)
mydata<-cvs_full_data[,colnames(cvs_full_data)%in%c("all_sr_1000","all_sr_100","all_sr_10","cc","soil_fert1")]
mydata2<-melt(mydata, id.vars=c("cc", "soil_fert1"))
mydata2$variable<-str_split_fixed(mydata2$variable, "_", 3)[,3]
mydata2$variable <- factor(mydata2$variable, levels = c("1000","100","10"))
sz<-10
mycols<-brewer.pal(n = 11, name = "RdYlBu")[c(2,4,10)]
sf<-ggplot(mydata2, aes(soil_fert1,value,group=variable,colour = variable)) +
geom_point(size=.1,alpha=0.5) +
geom_smooth(aes(fill = variable),method="loess", se=T, fullrange=T, level=0.95,span=.9,size=.5) +
scale_color_manual(values=mycols,name = bquote('scale ('*m^2*")")) +
scale_fill_manual(values=mycols,name = bquote('scale ('*m^2*")")) +
theme_bw() +
ggtitle("(a)")+
theme(legend.position = c(0.2, 0.82),
plot.title=element_text(hjust=-0.07,face="bold",vjust=-.5),
text=element_text(size=sz,family = "Times"),
axis.text=element_text(size=sz,family = "Times"),
legend.background = element_blank(),
legend.title=element_text(face="bold"),
legend.key = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank()) +
ylab("species richness (SR)") + xlab("soil fertility (SF)")
cc<-ggplot(mydata2, aes(cc,value,group=variable,colour = variable)) +
geom_point(size=.1, alpha=0.5) +
geom_smooth(aes(fill = variable),method="loess", se=TRUE, fullrange=FALSE, level=0.95,span=1,size=.5) +
scale_color_manual(values=mycols) +
scale_fill_manual(values=mycols) +
theme_bw() +
ggtitle("(b)")+
theme(plot.title=element_text(hjust=-0.07,face="bold",vjust=-.5),
text=element_text(size=sz,family = "Times"),
axis.text.x=element_text(size=sz,family = "Times"),
legend.background = element_blank(),
legend.key = element_blank(),
axis.text.y=element_blank(),
plot.background = element_blank(),
panel.grid = element_blank()) +
ylab("")+xlab("canopy cover (CC)")
legend_b <- get_legend(sf + theme(plot.margin=margin(l=0.1,r=0.5,unit="cm"),legend.text=element_text(family="Times"),legend.position="right"))
prow3<-plot_grid(sf+theme(plot.margin=margin(c(.2,0,.2,0.2),unit="cm")),
cc+theme(legend.position="none",plot.margin=margin(c(.2,0.2,.2,0),unit="cm")),
ncol = 2,
rel_widths = c(1,.9))
prow3
pdf("...Fig_2_loess.pdf",height=4,width=5)
prow3
dev.off()
# Fig. 3. Five-way polynomial - SF1 - all plots across scales -------
rm(list=ls())
library(ggplot2)
library(cowplot)
library(ecodist)
library(quantreg)
library(Hmisc)
library(MASS)
library(dplyr)
library(plyr)
library(rstan)
library(rstanarm)
setwd("...5_interaction_models")
cvs_full_data<-read.csv("...cvs_full_data_5cm.csv") # or the "all" one where SR is also scaled
names(cvs_full_data)
tmp1<-cvs_full_data[,c(2:22,63:68,61,50)]
tmp1<-tmp1[complete.cases(tmp1),]
tmp1[,1:27]<-apply(tmp1[,1:27],2,function(x){(x-min(x)) /(max(x)-min(x))})
tmp1[,28:29]<-scale(tmp1[,28:29])
poly_model_results<-as.data.frame(matrix(NA,27,28))
names(poly_model_results)<-c("plot_subset","scale","species_subset","soil_2.5","soil_25","soil_50","soil_75","soil_97.5",
"cc_2.5","cc_25","cc_50","cc_75","cc_97.5",
"soil_quad_2.5","soil_quad_25","soil_quad_50","soil_quad_75","soil_quad_97.5",
"cc_quad_2.5","cc_quad_25","cc_quad_50","cc_quad_75","cc_quad_97.5",
"soil_cc_2.5","soil_cc_25","soil_cc_50","soil_cc_75","soil_cc_97.5")
g<-1;i<-22
for (i in 1:27){
tmp<-tmp1[,c(i,28,29)]
names(tmp)<-c("SR","SF","CC")
poly_model<-stan_glm(SR~SF*CC+I(SF^2)+I(CC^2), data = tmp, QR = TRUE,chains = 4, iter = 10000)
beta_hat_CIs<-as.data.frame(summary(poly_model))[2:6,4:8]
poly_model_results[g,1]<-"all"
poly_model_results[g,2]<-strsplit(names(tmp1)[i], "_")[[1]][3]
poly_model_results[g,3]<-strsplit(names(tmp1)[i], "_")[[1]][1]
poly_model_results[g,4:28]<-cbind(beta_hat_CIs[1,],beta_hat_CIs[2,],beta_hat_CIs[3,],beta_hat_CIs[4,],beta_hat_CIs[5,])
print(i);g<-g+1
}
poly_model_results$scale <- factor(poly_model_results$scale, levels = rev(c("1000","400","100","10","1","01","001")))
poly_model_results$scale<-revalue(poly_model_results$scale, c("01"="0.1", "001"="0.01"))
poly_model_results$species_subset <- factor(poly_model_results$species_subset, levels = c( "ground","nw","all","w","tree"))
poly_model_results$species_subset<-revalue(poly_model_results$species_subset, c("nw"="herbaceous", "w"="woody","all"="total"))
write.csv(poly_model_results,"...poly_model_results.csv")
poly_model_results<-read.csv("...poly_model_results.csv")
library(RColorBrewer)
mycols<-c("chartreuse4", brewer.pal(n = 11, name = "PiYG")[8],"deepskyblue3",
brewer.pal(n = 11, name = "Spectral")[c(4,2)]) #11 before the 4
mysizes<-c(.25,.5,1,1.25,1.5)
poly_model_results$soil_cc_sig<-poly_model_results$cc_quad_sig<-poly_model_results$soil_quad_sig<-poly_model_results$cc_sig<-poly_model_results$soil_sig<-as.character(poly_model_results$species_subset)
i<-1
for (i in 1:dim(poly_model_results)[1]){
if (poly_model_results[i,4]<0 & poly_model_results[i,8]>0) {poly_model_results[i,29]<-"non"}
if (poly_model_results[i,9]<0 & poly_model_results[i,13]>0) {poly_model_results[i,30]<-"non"}
if (poly_model_results[i,14]<0 & poly_model_results[i,18]>0) {poly_model_results[i,31]<-"non"}
if (poly_model_results[i,19]<0 & poly_model_results[i,23]>0) {poly_model_results[i,32]<-"non"}
if (poly_model_results[i,24]<0 & poly_model_results[i,28]>0) {poly_model_results[i,33]<-"non"}
}
size_tmp<-14
soil_fert1_plot<-ggplot(poly_model_results, aes(x = as.factor(scale),y=soil_50,colour=factor(species_subset),shape=factor(species_subset),fill=soil_sig)) +
theme_bw()+geom_hline(yintercept = 0, linetype = "longdash") + coord_flip()+
geom_errorbar(aes(ymin=soil_2.5, ymax=soil_97.5), width=0,position=position_dodge(width=-0.6)) +
geom_point(position=position_dodge(width=-0.6),size=1.5) +
guides(fill=FALSE,shape=FALSE,colour = guide_legend(override.aes = list(shape =21:25,fill=mycols))) +
scale_shape_manual(values = c(21:25)) +
scale_fill_manual(values = c("non"="white", "ground"=mycols[1], "herbaceous"=mycols[2], "total"=mycols[3], "woody"=mycols[4], "tree"=mycols[5])) +
scale_colour_manual(values=mycols) +
ggtitle("(a)")+
theme(plot.title=element_text(hjust=-0.13,face="bold",vjust=-1.5,family="Times",size=size_tmp),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.text=element_text(family="Times",size=(size_tmp+2)),
legend.title=element_blank(),
axis.title=element_text(family="Times",size=size_tmp),
axis.text.y=element_text(family="Times"),
axis.text.x=element_text(family="Times",angle=25, hjust = 1))+
xlab(bquote('scale ('*m^2*")")) +
ylab(expression(paste("SF (", beta[1],")")))
soil2_plot<-ggplot(poly_model_results, aes(x = as.factor(scale),y=soil_quad_50,colour=factor(species_subset),shape=factor(species_subset),fill=soil_quad_sig)) +
theme_bw()+geom_hline(yintercept = 0, linetype = "longdash") + coord_flip()+
geom_errorbar(aes(ymin=soil_quad_2.5, ymax=soil_quad_97.5), width=0,position=position_dodge(width =-0.6)) +
geom_point(position=position_dodge(width=-0.6),size=1.5) +
guides(fill=FALSE,shape=FALSE,colour = guide_legend(override.aes = list(shape =21:25,fill=mycols))) +
scale_shape_manual(values = c(21:25)) +
scale_fill_manual(values = c("non"="white","ground"=mycols[1],"herbaceous"=mycols[2],"total"=mycols[3],"woody"=mycols[4],"tree"=mycols[5])) +
scale_size_manual(values=mysizes) + scale_colour_manual(values=mycols) +
ggtitle("(b)")+
theme(plot.title=element_text(hjust=-0.13,face="bold",vjust=-1.5,family="Times",size=size_tmp),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title=element_text(family="Times",size=size_tmp),
axis.text.y=element_text(family="Times"),
axis.text.x=element_text(family="Times",angle=25, hjust = 1))+
xlab(" ") + ylab(expression(paste(SF^2," (", beta[2],")")))+labs(colour="Species")
cc_plot<-ggplot(poly_model_results, aes(x = as.factor(scale),y=cc_50,colour=factor(species_subset),shape=factor(species_subset),fill=cc_sig)) +
theme_bw()+geom_hline(yintercept = 0, linetype = "longdash") + coord_flip()+
geom_errorbar(aes(ymin=cc_2.5, ymax=cc_97.5), width=0,position=position_dodge(width =-0.6)) +
geom_point(position=position_dodge(width=-0.6),size=1.5) +
guides(fill=FALSE,shape=FALSE,colour = guide_legend(override.aes = list(shape =21:25,fill=mycols))) +
scale_shape_manual(values = c(21:25)) +
scale_fill_manual(values = c("non"="white","ground"=mycols[1],"herbaceous"=mycols[2],"total"=mycols[3],"woody"=mycols[4],"tree"=mycols[5])) +
scale_size_manual(values=mysizes) + scale_colour_manual(values=mycols) +
ggtitle("(c)")+
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=size_tmp),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title=element_text(family="Times",size=size_tmp),
axis.text.y=element_text(family="Times"),
axis.text.x=element_text(family="Times",angle=25, hjust = 1))+
xlab(" ") + ylab(expression(paste("CC (", beta[3],")")))+labs(colour="Species")
cc2_plot<-ggplot(poly_model_results, aes(x = as.factor(scale),y=cc_quad_50,colour=factor(species_subset),shape=factor(species_subset),fill=cc_quad_sig)) +
theme_bw()+geom_hline(yintercept = 0, linetype = "longdash") + coord_flip()+
geom_errorbar(aes(ymin=cc_quad_2.5, ymax=cc_quad_97.5), width=0,position=position_dodge(width =-0.6)) +
geom_point(position=position_dodge(width=-0.6),size=1.5) +
guides(fill=FALSE,shape=FALSE,colour = guide_legend(override.aes = list(shape =21:25,fill=mycols))) +
scale_shape_manual(values = c(21:25)) +
scale_fill_manual(values = c("non"="white","ground"=mycols[1],"herbaceous"=mycols[2],"total"=mycols[3],"woody"=mycols[4],"tree"=mycols[5])) +
scale_size_manual(values=mysizes) + scale_colour_manual(values=mycols) +
ggtitle("(d)")+
theme(plot.title=element_text(hjust=-0.13,face="bold",vjust=-1.5,family="Times",size=size_tmp),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title=element_text(family="Times",size=size_tmp),
axis.text.y=element_text(family="Times"),
axis.text.x=element_text(family="Times",angle=25, hjust = 1))+
xlab(" ") + ylab(expression(paste(CC^2," (", beta[4],")")))+labs(colour="Species")
soil_cc_plot<-ggplot(poly_model_results, aes(x = as.factor(scale),y=soil_cc_50,colour=factor(species_subset),shape=factor(species_subset),fill=soil_cc_sig)) +
theme_bw()+geom_hline(yintercept = 0, linetype = "longdash") + coord_flip()+
geom_errorbar(aes(ymin=soil_cc_2.5, ymax=soil_cc_97.5), width=0,position=position_dodge(width =-0.6)) +
geom_point(position=position_dodge(width=-0.6),size=1.5) +
guides(fill=FALSE,shape=FALSE,colour = guide_legend(override.aes = list(shape =21:25,fill=mycols))) +
scale_shape_manual(values = c(21:25)) +
scale_fill_manual(values = c("non"="white","ground"=mycols[1],"herbaceous"=mycols[2],"total"=mycols[3],"woody"=mycols[4],"tree"=mycols[5])) +
scale_size_manual(values=mysizes) + scale_colour_manual(values=mycols) +
ggtitle("(e)")+
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=size_tmp),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title=element_text(family="Times",size=size_tmp),
axis.text.y=element_text(family="Times"),
axis.text.x=element_text(family="Times",angle=25, hjust = 1))+
xlab(" ") + ylab(expression(paste("SF:CC (", beta[5],")")))+labs(colour="Species")
setwd("...figs_tabs")
prow <- plot_grid(soil_fert1_plot + theme(legend.position="none",plot.margin=margin(l=0.1,r=0,unit="cm")),
soil2_plot + theme(legend.position="none",axis.text.y=element_blank(),plot.margin=margin(l=0,unit="cm")),
cc_plot + theme(legend.position="none",axis.text.y=element_blank(),plot.margin=margin(l=0,unit="cm")),
cc2_plot + theme(legend.position="none",axis.text.y=element_blank(),plot.margin=margin(l=0,unit="cm")),
soil_cc_plot + theme(legend.position="none",axis.text.y=element_blank(),plot.margin=margin(l=0,r=0.1,unit="cm")),
nrow = 1,rel_widths = c(1,.85,.85,.85,.85),align="h")
legend_b <- get_legend(soil_fert1_plot + theme(legend.title=element_blank(),legend.text=element_text(family="Times"),legend.position="bottom"))
pdf("Fig_3_five_way_poly_scale_class4.pdf",height=4.5,width=11)
plot_grid( prow, legend_b, ncol = 1, rel_heights = c(1, .07))
dev.off()
# Fig. 4 SR~SF+SF2 across cc and scale --------
dev.off()
rm(list=ls())
library(ecodist)
library(quantreg)
library(forcats)
library(Hmisc)
library(dplyr)
library(RColorBrewer)
library(cowplot)
library(plyr)
library(rstan)
library(rstanarm)
setwd("...6_soil_maxes")
cvs_full_data<-read.csv("...cvs_full_data_5cm.csv")
group_num<-3
cvs.dat <- cvs_full_data %>% mutate(cc_grp = cut2(cc, g=group_num))
my_groups<-unique(cvs.dat$cc_grp)
cvs_full_data<-cvs_full_data[,c(2:22,61,50)]
cvs_full_data<-cvs_full_data[complete.cases(cvs_full_data),]
cvs_full_data[,1:21]<-apply(cvs_full_data[,1:21],2,function(x){(x-min(x)) /(max(x)-min(x))})
cvs_full_data[,22:23]<-scale(cvs_full_data[,22:23])
sz=12
names(cvs_full_data)
soil_fert_plots<-list();models_range<-list();poly_models<-list()
scale<-1
for(scale in 1:7){
poly_model_results<-as.data.frame(matrix(NA,group_num*3,6))
names(poly_model_results)<-c("scale","species_subset","cc_level","soil_2.5","soil_50","soil_97.5")
tmp<-cvs_full_data[,c(scale,scale+7,scale+14,22,23)]
groups<-split(tmp, cut2(tmp$cc, g=group_num))
k<-3;i<-2;g<-1;gf<-2
for (k in 1:group_num){
for (gf in 1:3){
tmp_subset<-groups[[k]][,c(gf,4,5)]
names(tmp_subset)<-c("SR","SF","CC")
poly_model<-stan_glm(SR~SF+I(SF^2), data = tmp_subset, QR = TRUE,chains = 4, iter = 10000)
print(plot(poly_model, plotfun = "combo", prob = 0.5))
poly_model_results[g,1]<-strsplit(names(groups[[k]][,c(gf,4,5)])[1], "_")[[1]][3]
poly_model_results[g,2]<-strsplit(names(groups[[k]][,c(gf,4,5)])[1], "_")[[1]][1]
poly_model_results[g,3]<-as.character(sort(my_groups)[k])
poly_model_results[g,4:6]<-as.data.frame(summary(poly_model))[2,c(4,6,8)]
g<-g+1
}
}
poly_model_results$species_subset[poly_model_results$species_subset=="w"]<-"woody"
poly_model_results$species_subset[poly_model_results$species_subset=="nw"]<-"nonwoody"
poly_models[[scale]]<-poly_model_results
models_range[[scale]]<-c(min(poly_model_results$soil_2.5),max(poly_model_results$soil_97.5))
}
soil_betas<-do.call(rbind,poly_models)
names(soil_betas)<-c("scale","species_subset","group","min","mean","max")
write.csv(soil_betas,"...soil_betas.csv")
soil_betas<-read.csv("...oil_betas.csv")
soil_betas<-soil_betas[,-1]
soil_betas$scale[37:63]<-rep(c(1,0.1,0.01),each=9)
soil_betas$scale<-as.factor(soil_betas$scale)
soil_betas$scale <- factor(soil_betas$scale, levels = c(1000,400,100,10,1,0.1,0.01))
soil_betas_all<-subset(soil_betas,species_subset=="all")
mycolors<- brewer.pal(n = 11, name = "RdYlBu")[c(1:4,9:11)]
soil_betas_all_tmp<-soil_betas_all
SF<- ggplot(soil_betas_all, aes(x = as.factor(group),y=mean,col=factor(scale),shape=factor(scale))) +
theme_bw() + coord_flip() +
geom_point(size=1.5,position=position_dodge(width = -0.6)) +
guides(shape=FALSE,colour = guide_legend(override.aes = list(shape =c(15:19,15,18),fill=mycolors))) +
scale_shape_manual(values = c(15:19,15,18)) +
geom_errorbar(aes(ymin=min, ymax=max), size=0.5,position=position_dodge(width = -0.6), width=0) +
theme(plot.title=element_text(hjust=-0.1,vjust=-1.5,size=11,face="bold"),
text=element_text(size=sz,family = "Times"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.key.height = unit(.2, "cm"))+
geom_hline(yintercept = 0, linetype = "longdash") +
ggtitle("(a)") +
xlab("canopy cover levels") + labs(col=bquote('scale ('*m^2*")"),y=expression(SF,"~",SR (beta)))+
ylab(expression(paste("SF (", beta[1],")"))) + scale_color_manual(values = mycolors) +
geom_point(aes(x=1.255, y=soil_betas_all_tmp[1,5]), fill="white",colour=mycolors[1],size=1.5,shape=22) +
geom_point(aes(x=1.17, y=soil_betas_all_tmp[4,5]), fill="white",colour=mycolors[2],size=1.5,shape=21)
# SR~SF2 across cc and scale
for(scale in 1:7){
poly_model_results<-as.data.frame(matrix(NA,group_num*3,6))
names(poly_model_results)<-c("scale","species_subset","cc_level","soil_2.5","soil_50","soil_97.5")
tmp<-cvs_full_data[,c(scale,scale+7,scale+14,22,23)]
groups<-split(tmp, cut2(tmp$cc, g=group_num))
k<-3;i<-2;g<-1;gf<-2
for (k in 1:group_num){
for (gf in 1:3){
tmp_subset<-groups[[k]][,c(gf,4,5)]
names(tmp_subset)<-c("SR","SF","CC")
poly_model<-stan_glm(SR~SF+I(SF^2), data = tmp_subset, QR = TRUE,chains = 4, iter = 10000)
print(plot(poly_model, plotfun = "combo", prob = 0.5))
poly_model_results[g,1]<-strsplit(names(groups[[k]][,c(gf,4,5)])[1], "_")[[1]][3]
poly_model_results[g,2]<-strsplit(names(groups[[k]][,c(gf,4,5)])[1], "_")[[1]][1]
poly_model_results[g,3]<-as.character(sort(my_groups)[k])
poly_model_results[g,4:6]<-as.data.frame(summary(poly_model))[3,c(4,6,8)]
g<-g+1
}
}
poly_model_results$species_subset[poly_model_results$species_subset=="w"]<-"woody"
poly_model_results$species_subset[poly_model_results$species_subset=="nw"]<-"nonwoody"
poly_models[[scale]]<-poly_model_results
models_range[[scale]]<-c(min(poly_model_results$soil_2.5),max(poly_model_results$soil_97.5))
}
soil_betas2<-do.call(rbind,poly_models)
names(soil_betas2)<-c("scale","species_subset","group","min","mean","max")
write.csv(soil_betas2,"..soil_betas2.csv")
soil_betas2<-read.csv("...soil_betas2.csv")
soil_betas2<-soil_betas2[,-1]
soil_betas2$scale[37:63]<-rep(c(1,0.1,0.01),each=9)
soil_betas2$scale<-as.factor(soil_betas2$scale)
soil_betas2$scale <- factor(soil_betas2$scale, levels = c(1000,400,100,10,1,0.1,0.01))
soil_betas_all<-subset(soil_betas2,species_subset=="all")
mycolors<- brewer.pal(n = 11, name = "RdYlBu")[c(1:4,9:11)]
soil_betas_all_tmp2<-soil_betas_all
SF2<- ggplot(soil_betas_all, aes(x = as.factor(group),y=mean,col=factor(scale),shape=factor(scale))) +
theme_bw() + coord_flip() +
geom_point(size=1.5,position=position_dodge(width = -0.6)) +
guides(shape=FALSE,colour = guide_legend(override.aes = list(shape =c(15:19,15,18),fill=mycolors))) +
scale_shape_manual(values = c(15:19,15,18)) +
geom_errorbar(aes(ymin=min, ymax=max), size=0.5,position=position_dodge(width = -0.6), width=0) +
theme(plot.title=element_text(hjust=-0.1,vjust=-1.5,size=11,face="bold"),
text=element_text(size=sz,family = "Times"),
axis.text.y=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.key.height = unit(.2, "cm"))+
geom_hline(yintercept = 0, linetype = "longdash") +
ggtitle("(b)") +
xlab("") +
ylab(expression(paste(SF^2," (", beta[2],")"))) + scale_color_manual(values = mycolors)+
geom_point(aes(x=2.83, y=soil_betas_all_tmp2[18,5]), fill="white",colour=mycolors[6],size=1.5,shape=22) +
geom_point(aes(x=2.745, y=soil_betas_all_tmp2[21,5]), fill="white",colour=mycolors[7],size=1,shape=23) +
geom_point(aes(x=1.745, y=soil_betas_all_tmp2[20,5]), fill="white",colour=mycolors[7],size=1,shape=23)
setwd("...figs_tabs")
prow <- plot_grid(SF + theme(legend.position="none",plot.margin=margin(l=0.1,r=.1,unit="cm")),
SF2 + theme(legend.position="none",plot.margin=margin(l=0.1,r=.1,unit="cm")),rel_widths = c(1,0.8))
legend_b <- get_legend(SF + theme(plot.margin=margin(l=0,r=0,unit="cm"),legend.text=element_text(family="Times"),legend.position="right"))
pdf("Fig_4_five_way_poly_scale_class3.pdf",height=4,width=6)
plot_grid( prow, legend_b, ncol = 2, rel_widths = c(1, .2))
dev.off()
# Fig. 5 1000m2 comp across levels of cc (quartile) -----------------
dev.off()
rm(list=ls())
library(ecodist)
library(ggplot2)
library(quantreg)
library(Hmisc)
library(reshape2)
library(grid)
library(cowplot)
library(MASS)
library(dplyr)
library(plyr)
setwd("...6_soil_maxes")
cvs_full_data<-read.csv("...cvs_full_data_5cm.csv")
group_num<-3
cvs.dat <- cvs_full_data %>% mutate(cc_grp = cut2(cc, g=group_num))
my_groups<-unique(cvs.dat$cc_grp)
cvs_full_data[,c(61,50)]<-scale(cvs_full_data[,c(61,50)])
tmp<-cvs_full_data[,c(2,61,50)]
cvs.dat <- tmp %>% mutate(cc_grp = cut2(cc, g=group_num))
groups <- cvs.dat %>% split(.$cc_grp)
plot_list<-list();my_max1<-list();my_max2<-list()
mycols<- brewer.pal(n = 11, name = "RdYlBu")[c(1,9,10)]
sz<-10
g.num<-2
plots<-list()
for (g.num in 1:group_num){
tmp_data<-groups[[g.num]]
quart_model<-rq(tmp_data[,1]~tmp_data[,2]+I(tmp_data[,2]^2), tau=.9)
quad_model<-lm(all_sr_1000~soil_fert1+I(soil_fert1^2), data=tmp_data)
pframe <- with(tmp_data,expand.grid(soil_fert1=unique(soil_fert1)))
tmp_data$quadratic <- round(predict(quad_model,newdata=pframe,type="response"))
pframe <- with(tmp_data,expand.grid(soil_fert1=unique(soil_fert1)))
tmp_data$quartile <- round(predict(quart_model,newdata=pframe,type="response"))
tmp_data2<-tmp_data[,-c(3,4)]
head(tmp_data)
tmp_data3 <- melt(tmp_data2, id=c("all_sr_1000","soil_fert1"))
names(tmp_data3)[3]<-"Model"
tmp_data3$Model<-revalue(tmp_data3$Model, c("quadratic"="mean", "quartile"="90%"))
my_max1[[g.num]]<-round(c(sort(tmp_data[which.max(fitted(quart_model)),2]),fitted(quart_model)[which(fitted(quart_model)==max(fitted(quart_model)))]),1)
my_max2[[g.num]]<-round(c(sort(tmp_data[which.max(fitted(quad_model)),2]),fitted(quad_model)[which(fitted(quad_model)==max(fitted(quad_model)))]),1)
tmp_data3$Model <- factor(tmp_data3$Model, levels = c("90%", "mean"))
plots[[g.num]]<- ggplot(tmp_data3, aes(soil_fert1,all_sr_1000,group=Model)) +
geom_point(size=.1,color="grey80") +
geom_smooth(aes(y=value,group = Model,color=Model),size=0.5) +
scale_color_manual(values=mycols) +
ylim(min(tmp$all_sr_1000),max(tmp$all_sr_1000)) +
ylab("") + xlab("") +
coord_cartesian(clip = 'off') +
theme(text=element_text(size=sz,family = "Times"),
axis.text=element_text(size=sz,family = "Times"),
legend.position=c(.1,.8),legend.title=element_blank(),
legend.text=element_text(size=sz,family = "Times"),
legend.key.height = unit(.1, "cm"),
plot.title = element_text(size=sz,hjust = 0.5,face="plain",margin = margin(t = 0, b = -10))) +
scale_x_continuous(expand = c(0,0),limits = c(min(tmp$soil_fert1),max(tmp$soil_fert1))) +
guides(color=guide_legend(override.aes=list(fill=NA)))
}
dev.off()
prow <- plot_grid( plots[[3]]+ylab("")+theme(axis.text.x = element_blank(),plot.margin = unit(c(.2, .8, 0, .3), "cm")) + ggtitle(paste("CC: ",sort(my_groups)[3],sep="")) +
geom_segment(x = my_max1[[3]][1], y = my_max1[[3]][2], xend = my_max1[[3]][1], yend = 0,color=mycols[1], linetype = 2) +
geom_segment(x = my_max2[[3]][1], y = my_max2[[3]][2], xend = my_max2[[3]][1], yend = 0,color=mycols[2], linetype = 2) +
geom_point(aes(x=my_max1[[3]][1], y=my_max1[[3]][2]), colour=mycols[1],size=1) +
geom_point(aes(x=my_max2[[3]][1], y=my_max2[[3]][2]), colour=mycols[2],size=1) +
annotate("text", x = my_max1[[3]][1], y = my_max1[[3]][2]+20, label = round(my_max1[[3]][1],1),size=sz/3,color=mycols[1]) +
annotate("text", x = my_max2[[3]][1], y = my_max2[[3]][2]+20, label = round(my_max2[[3]][1],1),size=sz/3,color=mycols[3]),
plots[[2]]+ ylab(bquote('SR (1000'~m^2*")")) + theme(axis.text.x = element_blank(),legend.position = "none",plot.margin = unit(c(-.3, .8, 0, .15), "cm"))+ ggtitle(paste("CC: ",sort(my_groups)[2],sep="")) +
geom_segment(x = my_max1[[2]][1], y = my_max1[[2]][2], xend = my_max1[[2]][1], yend = 0,color=mycols[1], linetype = 2) +
geom_segment(x = my_max2[[2]][1], y = my_max2[[2]][2], xend = my_max2[[2]][1], yend = 0,color=mycols[2], linetype = 2) +
geom_point(aes(x=my_max1[[2]][1], y=my_max1[[2]][2]), colour=mycols[1],size=1) +
geom_point(aes(x=my_max2[[2]][1], y=my_max2[[2]][2]), colour=mycols[2],size=1) +
annotate("text", x = my_max1[[2]][1], y = my_max1[[2]][2]+20, label = round(my_max1[[2]][1],1),size=sz/3,color=mycols[1]) +
annotate("text", x = my_max2[[2]][1], y = my_max2[[2]][2]+20, label = round(my_max2[[2]][1],1),size=sz/3,color=mycols[3]),
plots[[1]]+ ylab("")+xlab("SF")+ theme(legend.position = "none",plot.margin = unit(c(-.3, .8, 0, .3), "cm"))+ ggtitle(paste("CC: ",sort(my_groups)[1],sep="")) +
geom_segment(x = my_max1[[1]][1], y = my_max1[[1]][2], xend = my_max1[[1]][1], yend = 0,color=mycols[1], linetype = 2) +
geom_segment(x = my_max2[[1]][1], y = my_max2[[1]][2], xend = my_max2[[1]][1], yend = 0,color=mycols[2], linetype = 2) +
geom_point(aes(x=my_max1[[1]][1], y=my_max1[[1]][2]), colour=mycols[1],size=1) +
geom_point(aes(x=my_max2[[1]][1], y=my_max2[[1]][2]), colour=mycols[2],size=1) +
annotate("text", x = my_max1[[1]][1], y = my_max1[[1]][2]+20, label = round(my_max1[[1]][1],1),size=sz/3,color=mycols[1]) +
annotate("text", x = my_max2[[1]][1], y = my_max2[[1]][2]+20, label = round(my_max2[[1]][1],1),size=sz/3,color=mycols[3]),
rel_heights=c(1,.8,.8),ncol = 1)
prow
pdf("...Fig_5_opt_3_levels_quart.pdf",height=4,width=5)
prow
dev.off()
#Fig. 6 SF_opt across scales -----------
rm(list = ls(all = TRUE))
library(tidyverse)
library(rstan)
library(rstanarm)
library(Hmisc)
cvs_full_data1<-read_csv("...cvs_full_data_5cm.csv") # or the "all" one where SR is also scaled
cvs_full_data1<-cvs_full_data1[,c(1:8,61,50)]
cvs_full_data<-cvs_full_data1[complete.cases(cvs_full_data1),]
cvs_full_data[,2:8]<-apply(cvs_full_data[,2:8],2,function(x){(x-min(x)) /(max(x)-min(x))})
cvs_full_data[,9:10]<-scale(cvs_full_data[,9:10])
n.grp<-3;scale<-6
my_boxes<-list();results<-list()
for (scale in c(1:7)){
cvs.dat<-cvs_full_data[,c(1,(scale+1),9,10)]
names(cvs.dat)[2]<-"SR"
cvs.dat <- cvs.dat %>% mutate(cc_grp = cut2(cc, g=n.grp))
cvs.dist <- cvs.dat %>% split(.$cc_grp)
fit_stan_lm <- function(dat) {stan_glm(SR ~ soil_fert1 + I(soil_fert1^2), data = dat)}
get_opt_fert_NB <- function(stan) {
stan.mcmc <- as.data.frame(stan) %>%
as_tibble()
colnames(stan.mcmc) <- c('int', 'soil_fert1', 'soil_fert12','reciprocal_dispersion')
stan.mcmc %>%
mutate(opt_fert = -0.5*soil_fert1/soil_fert12)
}
cvs.post <- cvs.dat %>%
split(.$cc_grp) %>%
map(. %>% fit_stan_lm())
cvs.fert_opt<- cvs.post %>%
map(. %>% get_opt_fert_NB())
nonsig_index<-lapply(cvs.post,function(x){
all(as.data.frame(summary(x))[3,4]<0,as.data.frame(summary(x))[3,8]>0)
})
pe<-cvs.fert_opt %>%
lapply(., function(x) {t(apply(x , 2 , quantile , probs = c(0.025, .5, 0.975) , na.rm = TRUE ))} ) %>%
lapply(., `[`,5,)
pe2<-as.data.frame(do.call(rbind,pe))
pe3<-cbind(row.names(pe2),pe2)
pe3[which(nonsig_index==TRUE),2:4]<-c(NA,NA,NA)
names(pe3)<-c("group","min","mean","max")
my_boxes[[scale]]<-pe3
results[[scale]]<-lapply(cvs.post,as.data.frame(summary))
}
pe4<-do.call(rbind,my_boxes)
pe5<-cbind(pe4,rep(c(1000,400,100,10,1,0.1,0.01), each=3))
unique(pe5$group)
names(pe5)[5]<-"scale"
pe5$scale<-as.factor(pe5$scale)
pe5$scale<-fct_rev(pe5$scale)
cvs.dat <- cvs_full_data1 %>% mutate(cc_grp = cut2(cc, g=n.grp))
levels(pe5$group)<-rev(unique(cvs.dat$cc_grp))
pe5$group<- levels(pe5$group)[c(2,3,1)]
write.csv(pe5,"...figs_tabs/Fig_6.csv")
library(RColorBrewer)
mycolors<- brewer.pal(n = 11, name = "RdYlBu")[c(1:4,9:11)]
sz=10
fert_opt_plot<-ggplot(pe5, aes(x = as.factor(group),y=mean,col=as.factor(scale),shape=factor(scale))) +
theme_bw() + coord_flip() +
geom_point(size=1,position=position_dodge(width = -0.6)) +
guides(shape=FALSE,colour = guide_legend(override.aes = list(shape =c(15:19,15,18),fill=mycolors))) +
geom_errorbar(aes(ymin=min, ymax=max), size=0.3,position=position_dodge(width = -0.6), width=0) +
scale_shape_manual(values = c(15:19,15,18)) +
theme(plot.title=element_text(hjust=0.5,size=11),
text=element_text(size=sz,family = "Times"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.key.height = unit(.2, "cm"))+
scale_color_manual(values = mycolors) +
xlab("canopy cover levels") + labs(col=bquote('scale ('*m^2*")"),y=expression(SF[maxSR]))
pdf("...Fig_6_opt.pdf",height=3,width=3)
fert_opt_plot
dev.off()
# Fig 7 Soil and cc by plot segmentation ########
dev.off()
rm(list=ls())
library(ecodist)
library(quantreg)
library(Hmisc)
library(MASS)
library(stringr)
library(reshape)
library(cowplot)
cvs_full_data<-read.csv("...cvs_full_data_5cm.csv")
setwd("...figs_tabs")
names(cvs_full_data)
test<-cvs_full_data[,c(61,57)] #23,51
test2<-test %>% melt(id.vars=c("ecogeo3"))
test2
library(RColorBrewer)
mycols<- brewer.pal(n = 11, name = "Spectral")[c(1,3,5)]
sz1<-11
sz2<-12
means<-aggregate(test2$value, by=list(Category=test2$ecogeo3), FUN=median)
means[,2]<-round(means[,2],1)
v4<-
ggplot(aes(y = value, x = factor(ecogeo3), fill = ecogeo3,colour=ecogeo3), data = subset(test2,variable=="soil_fert1")) +
geom_violin() +
xlab("")+ ylab("soil fertility (index)") +
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
plot.margin = unit(c(0.1, 0.2, 0, 0.2), "cm"),
axis.title=element_text(size = sz2,family="Times"),
axis.text.x = element_blank(),
axis.text.y = element_text(size = sz2,family="Times")) + scale_colour_manual(values=mycols) +
scale_fill_manual(values=mycols) +
ggtitle("(a)") +
annotate("text",x=1,y=means[1,2],label=means[1,2],colour="black",hjust=0.5,family="Times") +
annotate("text",x=2,y=means[2,2],label=means[2,2],colour="black",hjust=0.5,family="Times") +
annotate("text",x=3,y=means[3,2],label=means[3,2],colour="black",hjust=0.5,family="Times")
test<-cvs_full_data[,c(50,57)] #23,51
test2<-test %>% melt(id.vars=c("ecogeo3"))
test2$ecogeo3<-revalue(test2$ecogeo3, c("App Forests"="Appalachian Forests","Piedmont Forests"="Piedmont Forests"))
test2$ecogeo3<-factor(test2$ecogeo3, levels = c("Coastal Forests", "Piedmont Forests", "Appalachian Forests"))
means2<-aggregate(test2$value, by=list(Category=test2$ecogeo3), FUN=median)
means2[,2]<-round(means2[,2],1)
v5<- ggplot(aes(y = value, x = factor(ecogeo3), fill = ecogeo3,colour=ecogeo3), data = subset(test2,variable=="cc")) +
geom_violin() +
xlab("")+ ylab("canopy closure") +
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
plot.margin = unit(c(0, 0.2, 0, 0.2), "cm"),
axis.title=element_text(size = sz2,family="Times"),
axis.text.x = element_text(angle = 15,hjust = 1,size = sz2,family="Times"),
axis.text.y = element_text(size = sz2,family="Times")) + scale_colour_manual(values=mycols) +
scale_fill_manual(values=mycols) +
ggtitle("(b)") +
annotate("text",x=1,y=means2[1,2],label=means2[1,2],colour="black",hjust=0.5,family="Times") +
annotate("text",x=2,y=means2[2,2],label=means2[2,2],colour="black",hjust=0.5,family="Times") +
annotate("text",x=3,y=means2[3,2],label=means2[3,2],colour="black",hjust=0.5,family="Times")
test<-cvs_full_data[,c(61,60)] #23,51
test2<-test %>% melt(id.vars=c("ecogeo6"))
test2$ecogeo6<-as.character(test2$ecogeo6)
test2[test2$ecogeo6!="Longleaf Woodlands",1]<-"Not Longleaf"
mycols<- brewer.pal(n = 11, name = "Spectral")[c(8,11)]
means3<-aggregate(test2$value, by=list(Category=test2$ecogeo6), FUN=median)
means3[,2]<-round(means3[,2],1)
v6<- ggplot(aes(y = value, x = factor(ecogeo6), fill = ecogeo6,colour=ecogeo6), data = subset(test2,variable=="soil_fert1")) +
geom_violin(outlier.colour="grey",outlier.size=.1,position=position_dodge(.9)) +
xlab("")+ ylab("") +
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
plot.margin = unit(c(0.1, 0.2, 0, 0.2), "cm"),
axis.title=element_text(size = sz2,family="Times"),
axis.text.x = element_blank(),
axis.text.y = element_text(size = sz2,family="Times")) + scale_colour_manual(values=mycols) +
scale_fill_manual(values=mycols) +
ggtitle("(c)") +
annotate("text",x=1,y=means3[1,2],label=means3[1,2],colour="black",hjust=0.5,family="Times") +
annotate("text",x=2,y=means3[2,2],label=means3[2,2],colour="black",hjust=0.5,family="Times")
test<-cvs_full_data[,c(50,60)] #23,51
test2<-test %>% melt(id.vars=c("ecogeo6"))
test2$ecogeo6<-as.character(test2$ecogeo6)
test2[test2$ecogeo6!="Longleaf Woodlands",1]<-"Not Longleaf"
test2$ecogeo6<-as.factor(test2$ecogeo6)
means4<-aggregate(test2$value, by=list(Category=test2$ecogeo6), FUN=median)
means4[,2]<-round(means4[,2],1)
v7<-ggplot(aes(y = value, x = ecogeo6, fill = ecogeo6,colour=ecogeo6), data = subset(test2,variable=="cc")) +
geom_violin(outlier.colour="grey",outlier.size=.1,position=position_dodge(.9)) +
xlab("")+ ylab("") +
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
plot.margin = unit(c(0, 0.2, 0, 0.2), "cm"),
axis.title=element_text(size = sz2,family="Times"),
axis.text.x = element_text(angle = 15,hjust = 1,size = sz2,family="Times"),
axis.text.y = element_text(size = sz2,family="Times")) + scale_colour_manual(values=mycols) +
scale_fill_manual(values=mycols) +
ggtitle("(d)") +
annotate("text",x=1,y=means4[1,2],label=means4[1,2],colour="black",hjust=0.5,family="Times") +
annotate("text",x=2,y=means4[2,2],label=means4[2,2],colour="black",hjust=0.5,family="Times")
prow<-plot_grid(v4,v6,v5,v7,align="v",hjust = -.2, ncol = 2,rel_heights=c(.83,1))
pdf("fig_7_plots_2_box_plots_ecogeo3_2way.pdf",height=5,width=6)
prow
dev.off()
# Fig. 8 boxplots and violins of ecoregions ecogeo3 ------
dev.off()
rm(list=ls())
library(ecodist)
library(quantreg)
library(Hmisc)
library(MASS)
library(stringr)
library(reshape)
library(cowplot)
library(plyr)
cvs_full_data<-read.csv("...cvs_full_data_5cm.csv")
setwd("...figs_tabs")
names(cvs_full_data)
cvs_full_data$cc<-cvs_full_data$cc/100
table(cvs_full_data$ecogeo6)
annotate_size<-2.8
levels(cvs_full_data$ecogeo6)[5]<-"Riparian"
op <- par(mfrow = c(4,2), mar=c(1, 4, 1, 1) + 0.1)
test<-cvs_full_data[,c(57,2,4,8,9,11,15,16,18,22,63,65,66,68)] #23,51
test2<-test %>% melt(id.vars=c("ecogeo3"))
test3<-cbind(test2[,1],as.data.frame(str_split_fixed(test2$variable, "_", 3))[,-2],test2[,3])
head(test3)
names(test3)<-c("ecoregion","class","scale","value")
test3$class <- factor(test3$class, levels = c("all","ground","nw","w","tree"))
test3$class<-revalue(test3$class, c("nw"="herbaceous", "w"="woody"))
test3$ecoregion <- factor(test3$ecoregion, levels = c("Longleaf Woodlands","Coastal Forests","Riparian","Piedmont Forests","App Forests"))
test3$ecoregion<-revalue(test3$ecoregion, c("Riparian"="Swamp Floodplains", "App Forests"="Appalachian Forests"))
library(RColorBrewer)
mycols<-c("deepskyblue3","chartreuse4", brewer.pal(n = 11, name = "PiYG")[8],
brewer.pal(n = 11, name = "Spectral")[c(4,2)]) #11 before the 4
mycols<-brewer.pal(n = 11, name = "Spectral")[7:11]
ghj<-subset(test3,scale=="1000")
max(ghj$value)
sz1<-11
sz2<-11
test4<-test3[test3$scale=="1000",]
means<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=median)
maxes<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=max)
v1<-ggplot(aes(y = value, x = ecoregion, fill = class,colour=class), data = subset(test3,scale=="1000")) +
geom_boxplot(outlier.colour="grey",outlier.size=.1,position=position_dodge(.9))+
xlab("")+ ylab(bquote('SR (1000'~m^2*")")) +
#coord_cartesian(ylim = c(0, max(maxes[,3])))+
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
plot.margin = unit(c(0, 0.2, 0, -.1), "cm"),
axis.text.x = element_blank(),
axis.title=element_text(size = sz1,family="Times"),
axis.text.y = element_text(size = sz2,family="Times")) + scale_colour_manual(values=mycols)+
scale_fill_manual(values=mycols) +
annotate("text",x=.64,y=maxes[1,3],label=maxes[1,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.64,y=maxes[2,3],label=maxes[2,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.64,y=maxes[3,3],label=maxes[3,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=.82,y=maxes[4,3],label=maxes[4,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.82,y=maxes[5,3],label=maxes[5,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.82,y=maxes[6,3],label=maxes[6,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1,y=maxes[7,3],label=maxes[7,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2,y=maxes[8,3],label=maxes[8,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3,y=maxes[9,3],label=maxes[9,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.18,y=maxes[10,3],label=maxes[10,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.18,y=maxes[11,3],label=maxes[11,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3.18,y=maxes[12,3],label=maxes[12,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.36,y=maxes[13,3],label=maxes[13,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.36,y=maxes[14,3],label=maxes[14,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3.36,y=maxes[15,3],label=maxes[15,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("segment", x = .565, xend = .715, y = means[1,3], yend = means[1,3],colour = "black",size=.3) +
annotate("segment", x = 1.565, xend = 1.715, y = means[2,3], yend = means[2,3],colour = "black",size=.3) +
annotate("segment", x = 2.565, xend = 2.715, y = means[3,3], yend = means[3,3],colour = "black",size=.3) +
annotate("segment", x = .745, xend = .895, y = means[4,3], yend = means[4,3],colour = "black",size=.3) +
annotate("segment", x = 1.745, xend = 1.895, y = means[5,3], yend = means[5,3],colour = "black",size=.3) +
annotate("segment", x = 2.745, xend = 2.895, y = means[6,3], yend = means[6,3],colour = "black",size=.3) +
annotate("segment", x = .925, xend = 1.075, y = means[7,3], yend = means[7,3],colour = "black",size=.3) +
annotate("segment", x = 1.925, xend = 2.075, y = means[8,3], yend = means[8,3],colour = "black",size=.3) +
annotate("segment", x = 2.925, xend = 3.075, y = means[9,3], yend = means[9,3],colour = "black",size=.3) +
annotate("segment", x = 1.105, xend = 1.255, y = means[10,3], yend = means[10,3],colour = "black",size=.3) +
annotate("segment", x = 2.105, xend = 2.255, y = means[11,3], yend = means[11,3],colour = "black",size=.3) +
annotate("segment", x = 3.105, xend = 3.255, y = means[12,3], yend = means[12,3],colour = "black",size=.3) +
annotate("segment", x = 1.285, xend = 1.435, y = means[13,3], yend = means[13,3],colour = "black",size=.3) +
annotate("segment", x = 2.285, xend = 2.435, y = means[14,3], yend = means[14,3],colour = "black",size=.3) +
annotate("segment", x = 3.285, xend = 3.435, y = means[15,3], yend = means[15,3],colour = "black",size=.3) + ggtitle("(a)")
test4<-test3[test3$scale=="100",]
means2<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=median)
maxes2<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=max)
maxes2[,3]<-round(maxes2[,3])
v2<-ggplot(aes(y = value, x = ecoregion, fill = class,colour=class), data = subset(test3,scale=="100")) +
geom_boxplot(outlier.colour="grey",outlier.size=.1,position=position_dodge(.9))+
#coord_cartesian(ylim = c(0, 100))+
xlab("")+ ylab(bquote('SR (100'~m^2*")")) +
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
axis.title=element_text(size = sz1,family="Times"),
axis.text.x = element_blank(),
plot.margin = unit(c(-.3, 0.2, 0, 0.1), "cm"),
axis.text.y = element_text(size = sz2,family="Times")) +
scale_colour_manual(values=mycols)+ scale_fill_manual(values=mycols) +
annotate("text",x=.64,y=maxes2[1,3],label=maxes2[1,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.64,y=maxes2[2,3],label=maxes2[2,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.64,y=maxes2[3,3],label=maxes2[3,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=.82,y=maxes2[4,3],label=maxes2[4,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.82,y=maxes2[5,3],label=maxes2[5,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.82,y=maxes2[6,3],label=maxes2[6,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1,y=maxes2[7,3],label=maxes2[7,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2,y=maxes2[8,3],label=maxes2[8,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3,y=maxes2[9,3],label=maxes2[9,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.18,y=maxes2[10,3],label=maxes2[10,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.18,y=maxes2[11,3],label=maxes2[11,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3.18,y=maxes2[12,3],label=maxes2[12,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.36,y=maxes2[13,3],label=maxes2[13,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.36,y=maxes2[14,3],label=maxes2[14,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3.36,y=maxes2[15,3],label=maxes2[15,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("segment", x = .565, xend = .715, y = means2[1,3], yend = means2[1,3],colour = "black",size=.3) +
annotate("segment", x = 1.565, xend = 1.715, y = means2[2,3], yend = means2[2,3],colour = "black",size=.3) +
annotate("segment", x = 2.565, xend = 2.715, y = means2[3,3], yend = means2[3,3],colour = "black",size=.3) +
annotate("segment", x = .745, xend = .895, y = means2[4,3], yend = means2[4,3],colour = "black",size=.3) +
annotate("segment", x = 1.745, xend = 1.895, y = means2[5,3], yend = means2[5,3],colour = "black",size=.3) +
annotate("segment", x = 2.745, xend = 2.895, y = means2[6,3], yend = means2[6,3],colour = "black",size=.3) +
annotate("segment", x = .925, xend = 1.075, y = means2[7,3], yend = means2[7,3],colour = "black",size=.3) +
annotate("segment", x = 1.925, xend = 2.075, y = means2[8,3], yend = means2[8,3],colour = "black",size=.3) +
annotate("segment", x = 2.925, xend = 3.075, y = means2[9,3], yend = means2[9,3],colour = "black",size=.3) +
annotate("segment", x = 1.105, xend = 1.255, y = means2[10,3], yend = means2[10,3],colour = "black",size=.3) +
annotate("segment", x = 2.105, xend = 2.255, y = means2[11,3], yend = means2[11,3],colour = "black",size=.3) +
annotate("segment", x = 3.105, xend = 3.255, y = means2[12,3], yend = means2[12,3],colour = "black",size=.3) +
annotate("segment", x = 1.285, xend = 1.435, y = means2[13,3], yend = means2[13,3],colour = "black",size=.3) +
annotate("segment", x = 2.285, xend = 2.435, y = means2[14,3], yend = means2[14,3],colour = "black",size=.3) +
annotate("segment", x = 3.285, xend = 3.435, y = means2[15,3], yend = means2[15,3],colour = "black",size=.3) + ggtitle("(b)")
test4<-test3[test3$scale=="001",]
means3<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=median)
maxes3<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=max)
maxes3[,3]<-round(maxes3[,3],1)
v3<-ggplot(aes(y = value, x = ecoregion, fill = class,colour=class), data = subset(test3,scale=="001")) +
geom_boxplot(outlier.colour="grey",outlier.size=.1,position=position_dodge(.9))+
#coord_cartesian(ylim = c(0, 3))+
xlab("")+ ylab(bquote('SR (0.01'~m^2*")")) +
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
plot.margin = unit(c(-.3, 0.2, 0, 0.2), "cm"),
axis.title=element_text(size = sz1,family="Times"),
axis.text.x = element_text(angle = 10,hjust = 1,size = 10,family="Times"),
axis.text.y = element_text(size = sz2,family="Times")) + scale_colour_manual(values=mycols[c(1,3,4)]) +
scale_fill_manual(values=mycols[c(1,3,4)]) +
annotate("text",x=.7,y=maxes3[1,3],label=maxes3[1,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.7,y=maxes3[2,3],label=maxes3[2,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.7,y=maxes3[3,3],label=maxes3[3,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1,y=maxes3[4,3],label=maxes3[4,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2,y=maxes3[5,3],label=maxes3[5,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3,y=maxes3[6,3],label=maxes3[6,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.3,y=maxes3[7,3],label=maxes3[7,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.3,y=maxes3[8,3],label=maxes3[8,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=3.3,y=maxes3[9,3],label=maxes3[9,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("segment", x = .5725, xend = .8275, y = means3[1,3], yend = means3[1,3],colour = "black",size=.3) +
annotate("segment", x = 1.5725, xend = 1.8275, y = means3[2,3], yend = means3[2,3],colour = "black",size=.3) +
annotate("segment", x = 2.5725, xend = 2.8275, y = means3[3,3], yend = means3[3,3],colour = "black",size=.3) +
annotate("segment", x = .8725, xend = 1.1275, y = means3[4,3], yend = means3[4,3],colour = "black",size=.3) +
annotate("segment", x = 1.875, xend = 2.1275, y = means3[5,3], yend = means3[5,3],colour = "black",size=.3) +
annotate("segment", x = 2.875, xend = 3.1275, y = means3[6,3], yend = means3[6,3],colour = "black",size=.3) +
annotate("segment", x = 1.1725, xend = 1.4275, y = means3[7,3], yend = means3[7,3],colour = "black",size=.3) +
annotate("segment", x = 2.1725, xend = 2.4275, y = means3[8,3], yend = means3[8,3],colour = "black",size=.3) +
annotate("segment", x = 3.1725, xend = 3.4275, y = means3[9,3], yend = means3[9,3],colour = "black",size=.3) + ggtitle("(c)")
prow1 <- plot_grid(v1,v2,v3,hjust = 0,vjust=1,rel_heights=c(.8,.75,.88), nrow= 3,align="v")
test<-cvs_full_data[,c(60,2,4,8,9,11,15,16,18,22,63,65,66,68)] #23,51
test2<-test %>% melt(id.vars=c("ecogeo6"))
test3<-cbind(test2[,1],as.data.frame(str_split_fixed(test2$variable, "_", 3))[,-2],test2[,3])
head(test3)
names(test3)<-c("ecoregion","class","scale","value")
test3$class <- factor(test3$class, levels = c("all","ground","nw","w","tree"))
test3$class<-revalue(test3$class, c("nw"="herbaceous", "w"="woody"))
test3$ecoregion <- factor(test3$ecoregion, levels = c("Longleaf Woodlands","Coastal Forests","Riparian","Piedmont Forests","App Forests"))
test3$ecoregion<-revalue(test3$ecoregion, c("Riparian"="Swamp Floodplains", "App Forests"="Appalachian Forests"))
head(test3)
library(RColorBrewer)
test3$ecoregion<-as.character(test3$ecoregion)
test3[test3$ecoregion!="Longleaf Woodlands",1]<-"Not Longleaf"
test4<-test3[test3$scale=="1000",]
means4<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=median)
maxes4<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=max)
maxes4[,3]<-round(maxes4[,3])
v4<-ggplot(aes(y = value, x = factor(ecoregion), fill = class,colour=class), data = subset(test3,scale=="1000")) +
geom_boxplot(outlier.colour="grey",outlier.size=.1,position=position_dodge(.9))+
xlab("")+ ylab("")+
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",
legend.text=element_text(size=14),
plot.margin = unit(c(0, 0.2, 0, -.1), "cm"),
axis.text.x = element_blank(),
axis.title=element_text(size = sz1,family="Times"),
axis.text.y = element_text(size = sz2,family="Times")) +
scale_colour_manual(values=mycols)+ scale_fill_manual(values=mycols) +
annotate("text",x=.64,y=maxes4[1,3],label=maxes4[1,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.64,y=maxes4[2,3],label=maxes4[2,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=.82,y=maxes4[3,3],label=maxes4[3,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.82,y=maxes4[4,3],label=maxes4[4,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1,y=maxes4[5,3],label=maxes4[5,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2,y=maxes4[6,3],label=maxes4[6,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.18,y=maxes4[7,3],label=maxes4[7,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.18,y=maxes4[8,3],label=maxes4[8,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=1.36,y=maxes4[9,3],label=maxes4[9,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("text",x=2.36,y=maxes4[10,3],label=maxes4[10,3],colour="grey20",hjust=0.5,size=annotate_size,family="Times") +
annotate("segment", x = .565, xend = .715, y = means4[1,3], yend = means4[1,3],colour = "black",size=.3) +
annotate("segment", x = 1.565, xend = 1.715, y = means4[2,3], yend = means4[2,3],colour = "black",size=.3) +
annotate("segment", x = .745, xend = .895, y = means4[3,3], yend = means4[3,3],colour = "black",size=.3) +
annotate("segment", x = 1.745, xend = 1.895, y = means4[4,3], yend = means4[4,3],colour = "black",size=.3) +
annotate("segment", x = .925, xend = 1.075, y = means4[5,3], yend = means4[5,3],colour = "black",size=.3) +
annotate("segment", x = 1.925, xend = 2.075, y = means4[6,3], yend = means4[6,3],colour = "black",size=.3) +
annotate("segment", x = 1.105, xend = 1.255, y = means4[7,3], yend = means4[7,3],colour = "black",size=.3) +
annotate("segment", x = 2.105, xend = 2.255, y = means4[8,3], yend = means4[8,3],colour = "black",size=.3) +
annotate("segment", x = 1.285, xend = 1.435, y = means4[9,3], yend = means4[9,3],colour = "black",size=.3) +
annotate("segment", x = 2.285, xend = 2.435, y = means4[10,3], yend = means4[10,3],colour = "black",size=.3) + ggtitle("(d)")
test4<-test3[test3$scale=="100",]
means5<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=median)
maxes5<-aggregate(test4$value, by=list(test4$ecoregion,test4$class), FUN=max)
maxes5[,3]<-round(maxes5[,3])
v5<-ggplot(aes(y = value, x = ecoregion, fill = class,colour=class), data = subset(test3,scale=="100")) +
geom_boxplot(outlier.colour="grey",outlier.size=.1,position=position_dodge(.9))+
xlab("")+ ylab("")+
theme(plot.title=element_text(hjust=-0.12,face="bold",vjust=-1.5,family="Times",size=sz1),
legend.position = "none",