-
Notifications
You must be signed in to change notification settings - Fork 0
/
H2O2_incubation_experiment.Rmd
6189 lines (5589 loc) · 341 KB
/
H2O2_incubation_experiment.Rmd
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: "H2O2 incubation experiment analysis"
output: html_notebook
---
This notebook is for analysis of the data from H2O2 production and decay bottle experiments performed on whole water, 100 um filtered water, and 0.22 um filtered water collected from western Lake Erie.
Analysis ran on smitdere laptop unless otherwise indicated.
Regression analyses with enviornmental parameters and H2O2 production
----
This first section of the document is for dealing solely with regression plots of H2O2 production and decay with respiration, biomass, light vs. dark, and 100 um filtered water vs whole water. The 16S rRNA data will be covered in a following section.
Load the needed libraries:
```{r}
library(ggplot2)
library(patchwork)
library(dplyr)
library(tidyr)
library(vegan)
library(lubridate)
library(plotly)
library(reticulate)
library(modelr)
library(purrr)
library(broom)
library(reshape2)
library(sjPlot)
setwd("E:/Research/2019_Erie_Bloom/Prod_Decay_Experiments")
```
First, load the H2O2 and field data into an R dataframe:
```{r}
#Load the text-delimited data tables into R:
Prod_Decay_df <- read.table("Prod_Decay_Data.txt", header=TRUE, sep="\t")
Environ_df <- read.table("Prod_Decay_Environ_Data.txt", header=TRUE, sep="\t")
```
In the Prod_Decay_df is the modeled gross H2O2 production, the absolute Kloss_H2O2, net production rate of H2O2 in control bottles, and net decay rate of H2O2 in spiked bottles. It also has the summed error of squares between observed vs. modeled H2O2 concentrations. The replicate bottles are listed individually.
Samples from 2017 only had whole water and 0.22 um filtered treatments exposed to light.
Samples from 2018 and 2019 had the same as above, but also included dark bottles and 105 um filtered bottles (on certain dates, usually alternating).
The Environ_df has the other associated data that is paired with each set of H2O2 measurements. This includes Chlorophyll a, pH, CDOM (a305), DIC, DOC, respiration, primary production, and nutrient concentrations.
DIC, DOC, respiration, primary production, and UV data only exist for 2018 and 2019.
The goal is to determine how much of the total H2O2 production in water column is attributed to biological sources and how the biological production changes with growth rates, algal density, nutrient availability, and bacterial community composition.
Gross H2O2 production was estimated using a model described in Vermilyea et al. 2010 and Marsico et al. 2015. Environ. Sci. Technol. and Mar Sci. This model relies on paired 2L incubations of unamended water and water spiked with an H2O2 standard. The decay in the spiked incubation are used to correct the H2O2 production rates measured in the control bottles (the net observed production is lower than the gross due to simultaneous decay, which changes as a function of H2O2 concentration as described in the above references).
The models above assume that gross H2O2 production remains constant, which is likely invalid due to light dependent processes changing with solar zenith angle. This is particularly evident in a few of our bottle incubations, where the model cannot accurately fit the dynamics in H2O2 observed over the 9 hour experiments on a few dates. Dixon et al. 2013 have attempted to allow gross H2O2 production to change nonlinearly with using curve fitting parameters, but their experimental set up is not applicable to our data (they don't model the decreasing portion of diel peak). Their model parameters also hold no real meaning, as the terms behind any changing biological production over time are are unknown. To avoid aribtrary curve fitting, I am going to calculate gross H2O2 production assuming that H2O2 production rates and Kloss do not change, only using those experiments where the model fits the data. Observed net production and decay will also be investigated to include the data from all the experiments.
First, some formating; calculate average H2O2 production and decay rates (both gross and net) along with 95% confidence intervals for each date, then combine the averages and ranges with the other environmental data into one dataframe:
```{r}
Avg_Prod_Decay_df <- Prod_Decay_df %>%
group_by(Experiment_Date, Site, Condition, Model_Fit) %>%
#calculates number of observations, averages, and standard deviation of each column along the grouping specified above
summarise(n=n(), Sum_Error_Squares_avg=mean(Sum_Error_Squares, na.rm = TRUE),
Sum_Error_Squares_sd=sd(Sum_Error_Squares, na.rm = TRUE),
PH2O2_avg=mean(PH2O2, na.rm = TRUE), PH2O2_sd=sd(PH2O2, na.rm = TRUE),
Kloss_avg=mean(Kloss, na.rm = TRUE), Kloss_sd=sd(Kloss, na.rm = TRUE),
Net_production_avg=mean(Net_production, na.rm = TRUE),
Net_production_sd=sd(Net_production, na.rm = TRUE),
Net_decay_avg=mean(Net_decay, na.rm = TRUE), Net_decay_sd=sd(Net_decay, na.rm = TRUE),
Max_H2O2_avg=mean(Max_H2O2, na.rm = TRUE), Max_H2O2_sd=sd(Max_H2O2, na.rm = TRUE),
FC_Net_production_avg=mean(FC_Net_production, na.rm = TRUE),
FC_Net_production_sd=sd(FC_Net_production, na.rm = TRUE),
FC_Net_decay_avg=mean(FC_Net_decay, na.rm = TRUE), FC_Net_decay_sd=sd(FC_Net_decay, na.rm = TRUE),
FC_Max_H2O2_avg=mean(FC_Max_H2O2, na.rm = TRUE),
FC_Max_H2O2_sd=sd(FC_Max_H2O2, na.rm = TRUE)) %>%
#This part calculates 95% confidence intervals from the standard deviation
mutate(Sum_Error_Squares_CI=Sum_Error_Squares_sd/sqrt(n)*1.96) %>%
mutate(PH2O2_CI=PH2O2_sd/sqrt(n)*1.96) %>%
mutate(Kloss_CI=Kloss_sd/sqrt(n)*1.96) %>%
mutate(Net_production_CI=Net_production_sd/sqrt(n)*1.96) %>%
mutate(Net_decay_CI=Net_decay_sd/sqrt(n)*1.96) %>%
mutate(Max_H2O2_CI=Max_H2O2_sd/sqrt(n)*1.96) %>%
mutate(FC_Net_production_CI=FC_Net_production_sd/sqrt(n)*1.96) %>%
mutate(FC_Net_decay_CI=FC_Net_decay_sd/sqrt(n)*1.96)
```
How much of the total H2O2 produced was attributed to biology? Using the model published in Vermilyea et al. and Marisco et al?
I want to calculate the biotic production as the gross PH2O2 in unfiltered water - the net production in 0.22 um filtered water. This assumes that the net production observed in 0.22 um filtered water is gross abiotic production. While some H2O2 could be lost to abiotic processes, decay in spiked 0.22 um filtered bottles was only significantly nonzero on one date -- so it should be negligible here.
Calculate the biotic production in each experiment:
```{r}
#Calculate gross giotic H2O2 production rates:
Avg_Prod_Decay_df$Biotic_PH2O2 <- Avg_Prod_Decay_df$PH2O2_avg - Avg_Prod_Decay_df$FC_Net_production_avg
#Calculate the 95% CI of gross biotic H2O2 production rate using the propagation of error formula:
#I need to first set NANs in the confidence interval for filtered net production to zero, so that the error from PH2O2 is applied to biotic H2O2 when possible:
Avg_Prod_Decay_df$FC_Net_production_CI[is.na(Avg_Prod_Decay_df$FC_Net_production_CI)] <- 0
#Now calculate the 95% CI
Avg_Prod_Decay_df$Biotic_PH2O2_CI <- sqrt((Avg_Prod_Decay_df$PH2O2_CI^2) + (Avg_Prod_Decay_df$FC_Net_production_CI^2))
#What percentage of gross H2O2 production was attributed to biotic production?
Avg_Prod_Decay_df$Perc_Biotic_PH2O2 <- Avg_Prod_Decay_df$Biotic_PH2O2 / Avg_Prod_Decay_df$PH2O2_avg * 100
#Merge the H2O2 dataframe with the environmental dataframe so we can relate this information to the other measured parameters later:
Merged_Prod_Decay_df <- merge(Avg_Prod_Decay_df, Environ_df, by=c("Experiment_Date", "Site", "Condition"), all=TRUE)
#What is the mean and range of biotic PH2O2 for all experiments in which it could be calculated?
#We only want to consider whole water production in the light for now, so create a data frame of just those samples:
Merged_Prod_Decay_WL_only <- Merged_Prod_Decay_df[Merged_Prod_Decay_df$Condition == "WL", ]
#Calculate the number of complete observations
num_obs <- length(Merged_Prod_Decay_WL_only$Biotic_PH2O2[!(is.na(Merged_Prod_Decay_WL_only$Biotic_PH2O2))])
#Calculate the stats:
mean(Merged_Prod_Decay_WL_only$Biotic_PH2O2, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$Biotic_PH2O2, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$Biotic_PH2O2, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$Biotic_PH2O2, na.rm=TRUE)
```
Biotic production ranged from 9 - 244 nM/hr, with a mean of 73 +/- 24 nM/hr.
What does total gross H2O2 look like in comparison?
```{r}
mean(Merged_Prod_Decay_WL_only$PH2O2_avg, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$PH2O2_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$PH2O2_avg, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$PH2O2_avg, na.rm=TRUE)
```
What percentage of the total gross production can be attributed to biotic production?
```{r}
#Get the mean, 95% CI, and range for the percent of gross H2O2 production attributed to biotic production in the experiments:
mean(Merged_Prod_Decay_WL_only$Perc_Biotic_PH2O2, na.rm = TRUE)
(sd(Merged_Prod_Decay_WL_only$Perc_Biotic_PH2O2, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$Perc_Biotic_PH2O2, na.rm = TRUE)
max(Merged_Prod_Decay_WL_only$Perc_Biotic_PH2O2, na.rm = TRUE)
```
The percentage of total production attributed to biotic sources was on average 66 +/- 5%. Biotic production ranged from 44 - 94 %.
What do the numbers for Kloss look like?
```{r}
mean(Merged_Prod_Decay_WL_only$Kloss_avg, na.rm = TRUE)
(sd(Merged_Prod_Decay_WL_only$Kloss_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$Kloss_avg, na.rm = TRUE)
max(Merged_Prod_Decay_WL_only$Kloss_avg, na.rm = TRUE)
```
Let's compare net production in whole water and 0.22 um filtered water to check that this is working.
```{r}
#To make a grouped barplot, we need to rearrange the dataframe a little bit.
#Let's subset the dataframe to only include Experiment date and net production rates:
Net_prod_df <- subset(Merged_Prod_Decay_WL_only, select=c(Experiment_Date, Net_production_avg,
FC_Net_production_avg))
Net_prod_CI <- subset(Merged_Prod_Decay_WL_only, select=c(Experiment_Date, Net_production_CI,
FC_Net_production_CI))
#Need to convert to long format with a key for filtered and whole water, which will make plotting easier:
Net_prod_long <- gather(Net_prod_df, Prod_Type, Net_rate, c(Net_production_avg,FC_Net_production_avg))
#Change the names in the prod type so that they merge better.
Net_prod_long$Prod_Type <- gsub("FC_Net_production_avg", "Filt", Net_prod_long$Prod_Type)
Net_prod_long$Prod_Type <- gsub("Net_production_avg", "WW", Net_prod_long$Prod_Type)
#Make sure that the CI prod_type columns match those in the net_prod_long dataframe so that the two dataframes merge by this column
Net_prod_CI_long <- gather(Net_prod_CI, Prod_Type, CI, c(Net_production_CI,FC_Net_production_CI))
Net_prod_CI_long$Prod_Type <- gsub("FC_Net_production_CI", "Filt", Net_prod_CI_long$Prod_Type)
Net_prod_CI_long$Prod_Type <- gsub("Net_production_CI", "WW", Net_prod_CI_long$Prod_Type)
Net_prod_df <- merge(Net_prod_long, Net_prod_CI_long, by=c("Experiment_Date", "Prod_Type"), all = TRUE)
rm(Net_prod_long)
rm(Net_prod_CI_long)
#Plot:
Net_prod_barplot <- ggplot(Net_prod_df, aes(fill=Prod_Type, y=Net_rate, x=Experiment_Date)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=Net_rate-CI, ymax=Net_rate+CI), width=0.2,
position=position_dodge(0.9), size = 0.1) +
scale_x_discrete(limits=c("31-May-17", "13-Jun-17", "27-Jun-17", "6-Jul-17", "12-Jul-17", "18-Jul-17", "25-Jul-17", "1-Aug-17", "15-Aug-17", "22-Aug-17", "30-Aug-17", "31-Aug-17", "6-Sep-17", "12-Sep-17", "19-Sep-17", "26-Sep-17", "4-Oct-17", "5-Oct-17", "10-Jul-18", "24-Jul-18", "31-Jul-18", "3-Aug-18", "7-Aug-18", "10-Aug-18", "14-Aug-18", "21-Aug-18", "14-Sep-18", "18-Sep-18", "23-Jul-19", "2-Aug-19", "6-Aug-19", "24-Aug-19", "17-Sep-19", "20-Sep-19")) +
scale_fill_manual(values=c("red", "lightblue"), labels=c("0.22 um filtered", "Whole water")) +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 5, b = 0,
l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 5, b = 0,
l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_text(size = 14, color = "black", angle = 45, hjust = 1,
margin = margin(t = 5, r = 5, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "top",
legend.title = element_blank()) +
scale_y_continuous(breaks=seq(-50,300, by=50)) +
coord_cartesian(ylim=c(-50,300)) +
ylab(expression("Net H"[2]*"O"[2]*" production (nM/hr)"))
Net_prod_barplot
```
In many cases, net production in whole water is comparable to or greater than net production in 0.22 um filtered water. This is despite simultaneous decomposition of H2O2 in whole water bottles that is absent from 0.22 um filtered water, suggesting that total production in whole water is higher than that in the filtered water and there is some particle dependent production.
Plot of decay in whole and 0.22 um filtered water:
```{r}
#To make a grouped barplot, we need to rearrange the dataframe a little bit.
#Let's subset the dataframe to only include Experiment date and net decay rates:
Net_decay_df <- subset(Merged_Prod_Decay_WL_only, select=c(Experiment_Date, Net_decay_avg,
FC_Net_decay_avg))
Net_decay_CI <- subset(Merged_Prod_Decay_WL_only, select=c(Experiment_Date, Net_decay_CI,
FC_Net_decay_CI))
#Need to convert to long format with a key for filtered and whole water, which will make plotting easier:
Net_decay_long <- gather(Net_decay_df, Prod_Type, Net_decay, c(Net_decay_avg,FC_Net_decay_avg))
#Change the names in the prod type so that they merge better.
Net_decay_long$Prod_Type <- gsub("FC_Net_decay_avg", "Filt", Net_decay_long$Prod_Type)
Net_decay_long$Prod_Type <- gsub("Net_decay_avg", "WW", Net_decay_long$Prod_Type)
#Make sure that the CI prod_type columns match those in the net_prod_long dataframe so that the two dataframes merge by this column
Net_decay_CI_long <- gather(Net_decay_CI, Prod_Type, CI, c(Net_decay_CI,FC_Net_decay_CI))
Net_decay_CI_long$Prod_Type <- gsub("FC_Net_decay_CI", "Filt", Net_decay_CI_long$Prod_Type)
Net_decay_CI_long$Prod_Type <- gsub("Net_decay_CI", "WW", Net_decay_CI_long$Prod_Type)
Net_decay_df <- merge(Net_decay_long, Net_decay_CI_long, by=c("Experiment_Date", "Prod_Type"), all = TRUE)
rm(Net_decay_long)
rm(Net_decay_CI_long)
#Plot:
#Only plot 2017 for which there is Decay data for both whole water and 0.22 um filtered water:
Net_decay_barplot <- ggplot(Net_decay_df, aes(fill=Prod_Type, y=Net_decay, x=Experiment_Date)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=Net_decay-CI, ymax=Net_decay+CI), width=0.2, position=position_dodge(0.9),
size = 0.1) +
scale_x_discrete(limits=c("31-May-17", "13-Jun-17", "27-Jun-17", "6-Jul-17", "12-Jul-17", "18-Jul-17", "25-Jul-17", "1-Aug-17", "15-Aug-17", "22-Aug-17", "30-Aug-17", "31-Aug-17", "6-Sep-17", "12-Sep-17", "19-Sep-17", "26-Sep-17", "4-Oct-17", "5-Oct-17", "10-Jul-18", "24-Jul-18", "31-Jul-18", "3-Aug-18", "7-Aug-18", "10-Aug-18", "14-Aug-18", "21-Aug-18", "14-Sep-18", "18-Sep-18", "23-Jul-19", "2-Aug-19", "6-Aug-19", "24-Aug-19", "17-Sep-19", "20-Sep-19")) +
scale_fill_manual(values=c("red", "lightblue"), labels=c("0.22 um filtered", "Whole water")) +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 5, b = 0,
l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 5, b = 0,
l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_text(size = 14, color = "black", angle = 45, hjust = 1,
margin = margin(t = 5, r = 5, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "none") +
scale_y_continuous(breaks=seq(-100,500, by=100)) +
coord_cartesian(ylim=c(-100,500)) +
ylab(expression("Net H"[2]*"O"[2]*" decay (nM/hr)"))
Combo_net_plot <- Net_prod_barplot + Net_decay_barplot + plot_layout(ncol = 1)
Combo_net_plot
ggsave("Combo_net_plot.pdf", Combo_net_plot, width = 12, height = 8, units = "in", dpi=300)
ggsave("Net_prod_barplot.pdf", Net_prod_barplot, width = 3.5, height = 3.5, units = "in", dpi=300)
ggsave("Net_decay_barplot.pdf", Net_decay_barplot, width = 3.5, height = 3.5, units = "in", dpi=300)
```
I calculated photochemical production calculated from CDOM absorbance, whole water absorbance, average light intensity over the time window used to calculate net H2O2 production, and pathlength = bottle width. See the excel file named "Photo_H2O2_RATE_EXPERIMENT" for the calculations and input data. The goal was to compare these values to the net production measured in whole water to see if photochemistry could account for the net production (ignoring the decay in spiked bottles).
I copied the numbers in the excel file into a tab-delimited text file to import into R to make plots for the manuscript. The below chunk imports the data and makes a plot:
```{r}
#Import the text file into an R object:
Photo_vs_measured <- read.table("Photo_vs_Net.txt", header=TRUE, sep="\t")
#Make a barplot:
Photo_vs_measured_barplot <- ggplot(Photo_vs_measured, aes(fill=Type, y=H2O2_production_rate, x=Date)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(aes(ymin=H2O2_production_rate-CI, ymax=H2O2_production_rate+CI), width=0.2,
position=position_dodge(0.9), size = 0.1) +
scale_fill_manual(values=c("red", "lightblue")) +
scale_x_discrete(limits=c("23-Jul-19", "2-Aug-19", "6-Aug-19", "24-Aug-19", "17-Sep-19", "20-Sep-19")) +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 5, b = 0,
l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 5, b = 0,
l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_text(size = 14, color = "black", angle = 45, hjust = 1,
margin = margin(t = 5, r = 5, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "top",
legend.title = element_blank()) +
coord_cartesian(ylim=c(0,250)) +
ylab(expression("H"[2]*"O"[2]*" production rate (nM/hr)"))
Photo_vs_measured_barplot
ggsave("Photo_vs_measured_barplot.pdf", Photo_vs_measured_barplot, width = 12, height = 8, units = "in", dpi=300)
```
Only three dates had significantly different net measured and calculated photochemical H2O2 production rates. There is some uncertainty in the light pathlenght, and because higher pathlengths increase photochemical reaction rates, the true photochemical production rates could be higher if the light path was underestimated due to upwelling light and changes in solar zenith angle over the course of the day.
I calculated new photochemical production rates after multiplying the pathlength by a factor of 1.1 - 2 (see the "Photo_H2O2_RATE_EXPERIMENT" excel file). There is a summary tab in the tab labeled as "Path_multiplier_comparison" that I transfered to this markdown file to make it nicer looking for the manuscript.
```{r}
#Import the table from the excel sheet:
Path_Compare_Table <- read.table("Path_Compare_Table.txt", header=TRUE, sep="\t")
tab_df(Path_Compare_Table, alternate.rows = T, file="TableS1.doc") #print to a file.
```
Summarize net production and decay:
```{r}
#Net production in Whole water:
print("Net production in whole water")
mean(Merged_Prod_Decay_WL_only$Net_production_avg, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$Net_production_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$Net_production_avg, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$Net_production_avg, na.rm=TRUE)
#Net production in 0.22 um filtered water:
print("Net production in 0.22 um filtered water")
mean(Merged_Prod_Decay_WL_only$FC_Net_production_avg, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$FC_Net_production_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$FC_Net_production_avg, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$FC_Net_production_avg, na.rm=TRUE)
#Net decay in whole water:
print("Net decay in whole water")
mean(Merged_Prod_Decay_WL_only$Net_decay_avg, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$Net_decay_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$Net_decay_avg, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$Net_decay_avg, na.rm=TRUE)
#Net decay in 0.22 um filtered water:
print("Net decay in 0.22 um filtered water")
mean(Merged_Prod_Decay_WL_only$FC_Net_decay_avg, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$FC_Net_decay_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$FC_Net_decay_avg, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$FC_Net_decay_avg, na.rm=TRUE)
```
How much higher are gross H2O2 production rates than net H2O2 production rates on average?
```{r}
#Divide gross H2O2 production by net production:
Gross_vs_net <- Merged_Prod_Decay_WL_only$PH2O2_avg / Merged_Prod_Decay_WL_only$Net_production_avg
#Remove dates where there was zero net production:
Gross_vs_net <- Gross_vs_net[Gross_vs_net != Inf]
Gross_vs_net <- Gross_vs_net[Gross_vs_net > 0]
min(Gross_vs_net, na.rm = TRUE)
max(Gross_vs_net, na.rm = TRUE)
mean(Gross_vs_net, na.rm = TRUE)
sd(Gross_vs_net, na.rm = TRUE)
#Create a summary dataframe for the paper:
Gross_vs_net_df <- array(numeric(),c(34,6))
Gross_vs_net_df[,1] <- Merged_Prod_Decay_WL_only$Experiment_Date
Gross_vs_net_df[,2] <- Merged_Prod_Decay_WL_only$PH2O2_avg
Gross_vs_net_df[,3] <- Merged_Prod_Decay_WL_only$PH2O2_CI
Gross_vs_net_df[,4] <- Merged_Prod_Decay_WL_only$Net_production_avg
Gross_vs_net_df[,5] <- Merged_Prod_Decay_WL_only$Net_production_CI
Gross_vs_net_df[,6] <- Merged_Prod_Decay_WL_only$PH2O2_avg / Merged_Prod_Decay_WL_only$Net_production_avg
Gross_vs_net_df <- as.data.frame(Gross_vs_net_df)
colnames(Gross_vs_net_df) <- c("Experiment Date", "Total gross H2O2 production", "95 % CI", "Net H2O2 production", "95 % CI", "Fold difference")
Gross_vs_net_df$`Experiment Date` <- dmy(Gross_vs_net_df$`Experiment Date`)
Gross_vs_net_df <- Gross_vs_net_df[ order(Gross_vs_net_df$`Experiment Date`), ] #sort the dataframe by experiment date
#Remove entries where gross H2O2 production could not be calculated:
Gross_vs_net_df <- Gross_vs_net_df[ Gross_vs_net_df$`Total gross H2O2 production` != "NaN", ]
tab_df(Gross_vs_net_df, alternate.rows = T, file="TableS2.doc") #print to a file.
```
The Absolute H2O2 production CI for 8-30-2017 experiment is "NA" because this date only has n=1 because data from one replicate did not fit the model.
The Absolute H2O2 production CI for 9-14-2018 experiment is "NA" because this date only has n=1 because data from one spike replicate did not get H2O2 addition (by mistake). This precluded calculation of Absolute H2O2 production and Kloss.
Note that for the publication, I corrected the reported values for significant digits. I then fixed the fold differences calculated above to match what would be calculated from hand based on the reported values in the table.
Re-calulate mean and range for the fold difference using the values calculated by hand:
```{r}
mean(2.2, 2.9, 2.2, 2.4, 3.9, 18.4, 2.8, 8.3, 9.4, 23.3, 13.3, 5.6, 2.7, 3.0,
3.5, 2.1, 7.6, 3.6, 6.0, 2.7, 2.1, 3.5, 5.8, 1.9, 2.4)
min(2.2, 2.9, 2.2, 2.4, 3.9, 18.4, 2.8, 8.3, 9.4, 23.3, 13.3, 5.6, 2.7, 3.0,
3.5, 2.1, 7.6, 3.6, 6.0, 2.7, 2.1, 3.5, 5.8, 1.9, 2.4)
max(2.2, 2.9, 2.2, 2.4, 3.9, 18.4, 2.8, 8.3, 9.4, 23.3, 13.3, 5.6, 2.7, 3.0,
3.5, 2.1, 7.6, 3.6, 6.0, 2.7, 2.1, 3.5, 5.8, 1.9, 2.4)
```
Are net whole water and 0.22 um filtered production rates significantly different from each other on average?
```{r}
t.test(Merged_Prod_Decay_WL_only$Net_production_avg, Merged_Prod_Decay_WL_only$FC_Net_production_avg, paired = FALSE, alternative = "two.sided")
```
On which dates (if any), was net H2O2 production significantly different from 0.22 um filtered production?
```{r}
#Get a vector of dates to loop through:
dates <- unique(Prod_Decay_df$Experiment_Date)
#Remove the 1-Aug-17 and 24-Aug-19 samples, which only has n=1 for 0.22 um filtered production:
dates <- dates[ dates != "1-Aug-17"]
dates <- dates[ dates != "24-Aug-19"]
#For each date, do a t-test of net H2O2 production in whole water vs 0.22 um filtered water:
for (i in dates){
#Create a dataframe of samples from only that date
t_test_df <- filter(Prod_Decay_df, Experiment_Date == i & Condition == "WL")
#Print the date
print(i)
#run the t-test and print the result:
print(t.test(t_test_df$Net_production, t_test_df$FC_Net_production, paired = FALSE, alternative = "two.sided"))
}
```
Are net decay rates in whole water significantly different from those in 0.22 um filtered water?
```{r}
t.test(Merged_Prod_Decay_WL_only$Net_decay_avg, Merged_Prod_Decay_WL_only$FC_Net_decay_avg, paired = FALSE, alternative = "two.sided", na.rm = TRUE)
```
On which dates (if any) was net decay in 0.22 um filtered water significantly different from zero?
```{r}
#Get a vector of dates to loop through:
dates <- unique(Prod_Decay_df$Experiment_Date)
#Remove the dates from 2018 and 2019, which did not have filtered spike bottles:
drop <- c("10-Jul-18", "24-Jul-18", "31-Jul-18", "3-Aug-18", "7-Aug-18", "10-Aug-18",
"14-Aug-18", "21-Aug-18", "14-Sep-18", "18-Sep-18", "23-Jul-19", "2-Aug-19" , "6-Aug-19",
"24-Aug-19", "17-Sep-19", "20-Sep-19")
dates <- dates[!(dates %in% drop)]
#Create an empty list to store the t-test results:
FC_decay_t_tests <- list()
#For each date, do a t-test of net H2O2 decay in whole water vs 0.22 um filtered water:
for (i in dates){
#Create a dataframe of samples from only that date
t_test_df <- filter(Prod_Decay_df, Experiment_Date == i & Condition == "WL")
#run the t-test and print the result:
FC_decay_t_tests[[i]] <- t.test(t_test_df$FC_Net_decay, mu = 0, paired = FALSE, alternative = "two.sided")
print(FC_decay_t_tests[[i]])
}
```
Are net production and decay rates in whole water correlated?
```{r}
cor(Merged_Prod_Decay_WL_only$Net_production_avg, Merged_Prod_Decay_WL_only$Net_decay_avg, method="pearson", use="complete.obs")
WW_net_prod_vs_net_decay <- lm(Merged_Prod_Decay_WL_only$Net_decay_avg ~ Merged_Prod_Decay_WL_only$Net_production_avg,
na.action = na.omit)
summary(WW_net_prod_vs_net_decay)
```
Net decay is significantly correlated with net production rate:
```{r}
#Plot the relationship:
NetDecay_vs_Net_Prod <- ggplot(Merged_Prod_Decay_WL_only, aes(x=Net_production_avg, y=Net_decay_avg, color=as.factor(Year))) +
geom_smooth(method=lm, color="navy", fill="lightsteelblue", se=TRUE, size=0.25) +
geom_errorbar(aes(ymin=Net_decay_avg-Net_decay_CI, ymax=Net_decay_avg+Net_decay_CI), width=2,
size = 0.1) +
geom_errorbarh(aes(xmin=Net_production_avg-Net_production_CI, xmax=Net_production_avg+Net_production_CI), height=9, size = 0.1) +
geom_point(size = 0.8, alpha = 0.8) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Year") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 12, color = "black", margin = margin(t = 0, r = 5, b = 0, l = 0)),
axis.title.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 5, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 14, color = "black", margin = margin(t = 5, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 12, color = "black", margin = margin(t = 5, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position="none") +
#scale_y_continuous(breaks=c(0,50,100,150,200,250,300,350)) +
coord_cartesian(ylim=c(-50,500), xlim=c(-50,300)) +
xlab(expression("Net H"[2]*"O"[2]*" production (nM/hr)")) +
ylab(expression("Net H"[2]*"O"[2]*" decay (nM/hr)"))
NetDecay_vs_Net_Prod
```
Is absolute Kloss correlated with total gross H2O2 production?
```{r}
cor(Merged_Prod_Decay_WL_only$PH2O2_avg, Merged_Prod_Decay_WL_only$Kloss_avg, method="pearson", use="complete.obs")
WW_gross_prod_vs_kloss <- lm(Merged_Prod_Decay_WL_only$Kloss_avg ~ Merged_Prod_Decay_WL_only$PH2O2_avg,
na.action = na.omit)
summary(WW_gross_prod_vs_kloss)
```
There is a significant correlation between decay constants and total gross H2O2 production:
```{r}
#Plot the relationship:
AbsDecay_vs_GrossProd <- filter(Merged_Prod_Decay_WL_only, Model_Fit == "Yes") %>%
ggplot(aes(x=PH2O2_avg, y=Kloss_avg, color=as.factor(Year))) +
geom_smooth(method=lm, color="navy", fill="lightsteelblue", se=TRUE, size=0.25) +
geom_errorbar(aes(ymin=Kloss_avg-Kloss_CI, ymax=Kloss_avg+Kloss_CI), width=5, size = 0.1) +
geom_errorbarh(aes(xmin=PH2O2_avg-PH2O2_CI, xmax=PH2O2_avg+PH2O2_CI), height=0.02, size= 0.1) +
geom_point(size = 1, alpha = 0.8) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Year") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 12, color = "black", margin = margin(t = 0, r = 5, b = 0, l = 0)),
axis.title.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 5, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 14, color = "black", margin = margin(t = 5, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 12, color = "black", margin = margin(t = 5, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position="none") +
coord_cartesian(ylim=c(0,1), xlim=c(0,500)) +
xlab(expression("Total gross H"[2]*"O"[2]*" production (nM/hr)")) +
ylab(expression("Kloss,H2O2 (hr-1)"))
Combined_Prod_vs_Decay_plot <- NetDecay_vs_Net_Prod + AbsDecay_vs_GrossProd
Combined_Prod_vs_Decay_plot
ggsave("Combined_Prod_vs_Decay_plot.pdf", Combined_Prod_vs_Decay_plot, width = 8, height = 3.5, units = "in", dpi=300)
```
When did maximum H2O2 concentrations usually occur?
```{r}
#Convert exact times into approximate times for binning purposes:
Prod_Decay_df$Time_Max_H2O2 <- gsub('10:.*', '11:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('8:40', '9:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('8:01', '8:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('13:.*', '14:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('14:.*', '14:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('16:.*', '17:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('17:.*', '17:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('7:4.*', '8:00', Prod_Decay_df$Time_Max_H2O2)
Prod_Decay_df$Time_Max_H2O2 <- gsub('7:5.*', '8:00', Prod_Decay_df$Time_Max_H2O2)
#make a histogram of Time Max H2O2:
#Remove the dark bottles, where there was always no net change in H2O2
Time_Max_H2O2_histogram <- filter(Prod_Decay_df, Condition != "WD") %>%
ggplot(aes(x=Time_Max_H2O2)) +
geom_histogram(stat="count", color="black", fill="white") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
strip.background = element_blank(),
panel.spacing = unit(5, "mm"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.x = element_text(size = 16, color = "black", angle = 45, hjust = 1,
margin = margin(t = 6, r = 0, b = 0, l = 0)),
axis.title.x = element_text(size = 14, color = "black",
margin = margin(t = 6, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.y = element_text(size = 16, color = "black",
margin = margin(t = 0, r = 8, b = 0, l = 0)),
axis.text.y = element_text(size = 14, color = "black",
margin = margin(t = 0, r = 6, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.title = element_blank()) +
scale_x_discrete(limits=c("8:00", "9:00", "11:00", "14:00", "17:00", "18:00")) +
coord_cartesian(ylim=c(0,40)) +
xlab(expression("Time of Max H"[2]*"O"[2]*" (EDT)")) +
ylab("Count")
Time_Max_H2O2_histogram
ggsave("Time_Max_H2O2_histogram.pdf", Time_Max_H2O2_histogram, width = 3.5, height = 3.5, units = "in", dpi=300)
```
Most of the time, max H2O2 occurred between 14:00 and 17:00 EDT. In bottles were max H2O2 occurred at 8-9 am, there was net decay.
The time of max H2O2 in each replicate was the same on all but 6 dates.
Next, make a similar histogram of time of maximum H2O2 concentration, but for 0.22 um filtered bottles:
```{r}
#Convert exact times into approximate times for binning purposes:
Prod_Decay_df$FC_Time_Max_H2O2 <- gsub('17:.*', '17:00', Prod_Decay_df$FC_Time_Max_H2O2)
Prod_Decay_df$FC_Time_Max_H2O2 <- gsub('16:.*', '17:00', Prod_Decay_df$FC_Time_Max_H2O2)
Prod_Decay_df$FC_Time_Max_H2O2 <- gsub('7:5.*', '8:00', Prod_Decay_df$FC_Time_Max_H2O2)
Prod_Decay_df$FC_Time_Max_H2O2 <- gsub('13:5.*', '14:00', Prod_Decay_df$FC_Time_Max_H2O2)
#make a histogram of Time Max H2O2:
#Remove the dark bottles, where there was always no net change in H2O2
#Remove the FL bottles, because that would created redundant 0.22 um filtered entries
Time_Max_H2O2_histogram_022um <- filter(Prod_Decay_df, Condition == "WL") %>%
ggplot(aes(x=FC_Time_Max_H2O2)) +
geom_histogram(stat="count", color="black", fill="white") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
strip.background = element_blank(),
panel.spacing = unit(5, "mm"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.x = element_text(size = 16, color = "black", angle = 45, hjust = 1,
margin = margin(t = 6, r = 0, b = 0, l = 0)),
axis.title.x = element_text(size = 14, color = "black",
margin = margin(t = 6, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.y = element_text(size = 16, color = "black",
margin = margin(t = 0, r = 8, b = 0, l = 0)),
axis.text.y = element_text(size = 14, color = "black",
margin = margin(t = 0, r = 6, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.title = element_blank()) +
scale_x_discrete(limits=c("8:00", "9:00", "11:00", "14:00", "17:00", "18:00")) +
coord_cartesian(ylim=c(0,60)) +
xlab(expression("Time of Max H"[2]*"O"[2]*" (EDT)")) +
ylab("Count")
Time_Max_H2O2_histogram_022um
ggsave("Time_Max_H2O2_histogram_022um.pdf", Time_Max_H2O2_histogram_022um, width = 3.5, height = 3.5, units = "in", dpi=300)
```
The two rows removed had NAs in their fields. These were the two replicate bottles with problems in the concentration data.
What was max H2O2 concentration in whole water on average?
```{r}
print("Whole water concentrations")
mean(Merged_Prod_Decay_WL_only$Max_H2O2_avg, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$Max_H2O2_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$Max_H2O2_avg, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$Max_H2O2_avg, na.rm=TRUE)
print("Filtered control concentrations")
mean(Merged_Prod_Decay_WL_only$FC_Max_H2O2_avg, na.rm=TRUE)
(sd(Merged_Prod_Decay_WL_only$FC_Max_H2O2_avg, na.rm=TRUE)/sqrt(num_obs))*1.96
min(Merged_Prod_Decay_WL_only$FC_Max_H2O2_avg, na.rm=TRUE)
max(Merged_Prod_Decay_WL_only$FC_Max_H2O2_avg, na.rm=TRUE)
```
Is the error between observed H2O2 concentrations and model fit related to any environmental parameters?
```{r}
#I want to include 105 um filtered water in the regression analysis, but remove the dark bottles which always fit the data, so I'll make a dataframe without the dark data:
Merged_Prod_Decay_no_WD <- Merged_Prod_Decay_df[Merged_Prod_Decay_df$Condition != "WD", ]
#Calculate the correlations between error sum of squares (SSE) and the environmental variables:
#This is a vector of columns to regress Biotic PH2O2 over:
vars <- c("Kloss_avg", "Chla", "DIC", "H2CO3", "HCO3", "CO3", "DOC", "Resp", "PrimProd", "CDOM", "Day_Integrated_UVA", "Day_Integrated_UVB", "Day_Integrated_UV", "pH", "TP", "TDP", "Nitrate", "NH4", "SRP", "Incubation_Temp", "Incubation_Temp_SD", "peakA", "peakC", "peakT", "C_A_ratio", "T_A_ratio", "IntFlour", "FI", "SlopeRatio")
#Create empty lists to save results of the loop into:
WL_cor_results_SSE <- list()
WL_lm_results_SSE <- list()
WL_lm_results_SSE_table <- list()
#Loop through each item of the vector and find the correlation with SSE in whole water light samples:
for (i in vars){
#Calculate the Pearson's R correlation statistic, ignoring sample pairs which have NAs
WL_cor_results_SSE[[i]] <- cor(Merged_Prod_Decay_no_WD[, colnames(Merged_Prod_Decay_no_WD) %in% i],
Merged_Prod_Decay_no_WD$Sum_Error_Squares_avg, method = "pearson", use = "complete.obs")
#Build a linear model for each correlation:
WL_lm_results_SSE[[i]] <- lm(Merged_Prod_Decay_no_WD$Sum_Error_Squares_avg ~ Merged_Prod_Decay_no_WD[, colnames(Merged_Prod_Decay_no_WD) %in% i], na.action = na.omit)
#Get the statistics into a table format
WL_lm_results_SSE_table[[i]] <- glance(WL_lm_results_SSE[[i]])
}
#print out the stats for the ones that are significant:
for (i in 1:length(WL_lm_results_SSE)){
if (WL_lm_results_SSE_table[[i]]$p.value < 0.05){
print(vars[i])
print(WL_cor_results_SSE[[i]])
print(WL_lm_results_SSE_table[[i]]$p.value)
print(WL_lm_results_SSE_table[[i]]$r.squared)
}
}
```
There is a significant relationship between SSE, chlorophyll concentration, and primary production. I want to compare this to the relationship with CDOM, so get that stats for that one as well.
```{r}
WL_cor_results_SSE[["CDOM"]]
WL_lm_results_SSE_table[["CDOM"]]$p.value
WL_lm_results_SSE_table[["CDOM"]]$r.squared
```
Plot the regression:
```{r}
#Plot Chlorophyll vs SSE
Chla_vs_SSE <- ggplot(Merged_Prod_Decay_no_WD, aes(x=Chla, y=Sum_Error_Squares_avg,
color=as.factor(Year))) +
geom_smooth(method=lm, color="navy", fill="lightsteelblue", se=TRUE, size=0.5) +
geom_errorbar(aes(ymin=Sum_Error_Squares_avg-Sum_Error_Squares_CI,
ymax=Sum_Error_Squares_avg+Sum_Error_Squares_CI), width=1.5, size=0.08) +
geom_errorbarh(aes(xmin=Chla-Chla_CI, xmax=Chla+Chla_CI), height=20000, size=0.08) +
geom_point(size = 0.5, alpha = 0.8) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Year") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position="none") +
coord_cartesian(ylim=c(-200,1200000), xlim=c(0,200)) +
xlab(expression("Chlorophyll a ("*mu*"g/L)")) +
ylab(expression("Sum error of squares in PH2O2 model"))
#Plot Primary Production vs SSE
PrimProd_vs_SSE <- ggplot(Merged_Prod_Decay_no_WD, aes(x=PrimProd, y=Sum_Error_Squares_avg,
color=as.factor(Year))) +
geom_smooth(method=lm, color="navy", fill="lightsteelblue", se=TRUE, size=0.5) +
geom_errorbar(aes(ymin=Sum_Error_Squares_avg-Sum_Error_Squares_CI,
ymax=Sum_Error_Squares_avg+Sum_Error_Squares_CI), width=1.5, size=0.08) +
geom_errorbarh(aes(xmin=PrimProd-PrimProd_CI, xmax=PrimProd+PrimProd_CI), height=11, size=0.08) +
geom_point(size = 0.5, alpha = 0.8) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Year") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "none") +
coord_cartesian(ylim=c(-200,1200000), xlim=c(0,90)) +
scale_x_continuous(breaks=seq(0,90, by=15)) +
xlab(expression("Primary Production ("*mu*"M C/hr)")) +
ylab(expression("Sum error of squares in PH2O2 model"))
#Plot CDOM vs SSE
CDOM_vs_SSE <- ggplot(Merged_Prod_Decay_no_WD, aes(x=CDOM, y=Sum_Error_Squares_avg,
color=as.factor(Year))) +
geom_smooth(method=lm, color="navy", fill="lightsteelblue", se=TRUE, size=0.5) +
geom_errorbar(aes(ymin=Sum_Error_Squares_avg-Sum_Error_Squares_CI,
ymax=Sum_Error_Squares_avg+Sum_Error_Squares_CI), width=0.5, size=0.08) +
geom_errorbarh(aes(xmin=CDOM-CDOM_CI, xmax=CDOM+CDOM_CI), height=11, size=0.08) +
geom_point(size = 0.5, alpha = 0.8) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Year") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_blank(),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "none") +
coord_cartesian(ylim=c(-200,1200000), xlim=c(0,30)) +
scale_x_continuous(breaks=seq(0,30, by=5)) +
xlab(expression("CDOM absorbance (a305)")) +
ylab(expression("Sum error of squares in PH2O2 model"))
Combo_SSE_regression_plot <- Chla_vs_SSE + PrimProd_vs_SSE + CDOM_vs_SSE
Combo_SSE_regression_plot
```
Repeat the above analysis, but without the points where the model had a poor fit.
```{r}
#Exclude the data points where there was poor fit to the model
Merged_Prod_Decay_no_WD_no_poor_fit <- Merged_Prod_Decay_no_WD[Merged_Prod_Decay_no_WD$Model_Fit != "No", ]
#Calculate the correlations between error sum of squares (SSE) and the environmental variables:
#This is a vector of columns to regress Biotic PH2O2 over:
vars <- c("Kloss_avg", "Chla", "DIC", "H2CO3", "HCO3", "CO3", "DOC", "Resp", "PrimProd", "CDOM", "Day_Integrated_UVA", "Day_Integrated_UVB", "Day_Integrated_UV", "pH", "TP", "TDP", "Nitrate", "NH4", "SRP", "Incubation_Temp", "Incubation_Temp_SD", "peakA", "peakC", "peakT", "C_A_ratio", "T_A_ratio", "IntFlour", "FI", "SlopeRatio")
#Create empty lists to save results of the loop into:
WL_cor_results_SSE_no_poor_fit <- list()
WL_lm_results_SSE_no_poor_fit <- list()
WL_lm_results_SSE_no_poor_fit_table <- list()
#Loop through each item of the vector and find the correlation with SSE in whole water light samples:
for (i in vars){
#Calculate the Pearson's R correlation statistic, ignoring sample pairs which have NAs
WL_cor_results_SSE_no_poor_fit[[i]] <- cor(Merged_Prod_Decay_no_WD_no_poor_fit[, colnames(Merged_Prod_Decay_no_WD_no_poor_fit) %in% i],
Merged_Prod_Decay_no_WD_no_poor_fit$Sum_Error_Squares_avg, method = "pearson", use = "complete.obs")
#Build a linear model for each correlation:
WL_lm_results_SSE_no_poor_fit[[i]] <- lm(Merged_Prod_Decay_no_WD_no_poor_fit$Sum_Error_Squares_avg ~ Merged_Prod_Decay_no_WD_no_poor_fit[, colnames(Merged_Prod_Decay_no_WD_no_poor_fit) %in% i], na.action = na.omit)
#Get the statistics into a table format
WL_lm_results_SSE_no_poor_fit_table[[i]] <- glance(WL_lm_results_SSE_no_poor_fit[[i]])
}
#print out the stats for the ones that were cosindered previously:
for (i in c("Chla", "PrimProd", "CDOM")){
print(i)
print(WL_cor_results_SSE_no_poor_fit[[i]])
print(WL_lm_results_SSE_no_poor_fit_table[[i]]$p.value)
print(WL_lm_results_SSE_no_poor_fit_table[[i]]$r.squared)
}
```
Import the dataframe:
```{r}
Poor_fit_curves_df <- read.table("Poor_fit_curves_df.txt", header=TRUE, sep="\t")
```
Plot:
```{r}
#Plot for WL 22-Aug-17:
WL_22Aug17_1 <- filter(Poor_fit_curves_df, Date == "22-Aug-17" & Rep == 1) %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 22-Aug-17 Rep 1") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "top") +
coord_cartesian(ylim=c(0,2400), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,2400, by=400)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_22Aug17_2 <- filter(Poor_fit_curves_df, Date == "22-Aug-17" & Rep == 2) %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 22-Aug-17 Rep 2") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "none") +
coord_cartesian(ylim=c(0,2400), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,2400, by=400)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_22Aug17_plot <- WL_22Aug17_1 + WL_22Aug17_2
WL_22Aug17_plot
ggsave("WL_22Aug17_plot.pdf", WL_22Aug17_plot, width = 8, height = 3.5, units = "in", dpi=300)
```
```{r}
#Plot for WL 30-Aug-17:
WL_30Aug17_1 <- filter(Poor_fit_curves_df, Date == "30-Aug-17" & Rep == 1) %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 30-Aug-17 Rep 1") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "top") +
coord_cartesian(ylim=c(0,1200), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,1200, by=200)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_30Aug17_2 <- filter(Poor_fit_curves_df, Date == "30-Aug-17" & Rep == 2) %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 30-Aug-17 Rep 2") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "none") +
coord_cartesian(ylim=c(0,1200), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,1200, by=200)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_30Aug17_plot <- WL_30Aug17_1 + WL_30Aug17_2
WL_30Aug17_plot
ggsave("WL_30Aug17_plot.pdf", WL_30Aug17_plot, width = 8, height = 3.5, units = "in", dpi=300)
```
```{r}
#Plot for WL 23-Jul-19:
WL_23Jul19_1 <- filter(Poor_fit_curves_df, Date == "23-Jul-19" & Rep == 1) %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 23-Jul-19 Rep 1") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "top") +
coord_cartesian(ylim=c(0,1400), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,1400, by=200)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_23Jul19_2 <- filter(Poor_fit_curves_df, Date == "23-Jul-19" & Rep == 2) %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 23-Jul-19 Rep 2") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "none") +
coord_cartesian(ylim=c(0,1400), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,1400, by=200)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_23Jul19_plot <- WL_23Jul19_1 + WL_23Jul19_2
WL_23Jul19_plot
ggsave("WL_23Jul19_plot.pdf", WL_23Jul19_plot, width = 8, height = 3.5, units = "in", dpi=300)
```
```{r}
#Plot for WL 6-Aug-19:
WL_6Aug19_1 <- filter(Poor_fit_curves_df, Date == "6-Aug-19" & Rep == 1 & Condition == "WL") %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 6-Aug-19 Rep 1") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "top") +
coord_cartesian(ylim=c(0,800), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,800, by=200)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_6Aug19_2 <- filter(Poor_fit_curves_df, Date == "6-Aug-19" & Rep == 2 & Condition == "WL") %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,
ymax=measured_H2O2+measured_H2O2_se), width=0.5, size=0.08) +
scale_color_manual(values=c("red", "blue", "orange"), name = "Bottle:") +
ggtitle("Whole water 6-Aug-19 Rep 2") +
theme_classic() +
theme(plot.background = element_rect(color = "NA"),
axis.line.x = element_line(size=0.1),
axis.line.y = element_line(size=0.1),
axis.text.y = element_text(size = 14, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.y = element_text(size = 16, color = "black", margin = margin(t = 0, r = 10, b = 0, l = 0)),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 16, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10, r = 0, b = 0, l = 0)),
axis.ticks.length = unit(-0.1, "cm"),
axis.ticks = element_line(size=0.1),
legend.position = "none") +
coord_cartesian(ylim=c(0,800), xlim=c(0,9)) +
scale_x_continuous(breaks = seq(0,9, by=3)) +
scale_y_continuous(breaks = seq(0,800, by=200)) +
xlab(expression("Incubation time (hours)")) +
ylab(expression("[H"[2]*"O"[2]*"] (nM)"))
WL_6Aug19_plot <- WL_6Aug19_1 + WL_6Aug19_2
WL_6Aug19_plot
ggsave("WL_6Aug19_plot.pdf", WL_6Aug19_plot, width = 8, height = 3.5, units = "in", dpi=300)
```
```{r}
#Plot for Fl 6-Aug-19:
FL_6Aug19_1 <- filter(Poor_fit_curves_df, Date == "6-Aug-19" & Rep == 1 & Condition == "FL") %>%
ggplot(aes(x=Hours_after_spike, y=measured_H2O2, color=Bottle_type)) +
geom_point() +
geom_line(aes(x=Hours_after_spike, y=model_H2O2, color=Bottle_type), linetype="dashed") +
geom_errorbar(aes(ymin=measured_H2O2-measured_H2O2_se,