Dissolved Oxygen Recovery Project

AZMP cruise data was searched manually to find cruises where there were both Winkler and probe measurements. Within Gordana’s files there were three cruises which qualified, HUD2011004, HUD2011043 and HUD2012042

These files were read in as csv’s using the BCD flagged files to provide the most accurate and up to date QC. Once the files were copied from the shared drive to the local drive they were read in as seen below (caution this section will not run unless file paths are exact)

#read in BCD files
H11004_bcd <- read.csv('HUD2011004/HUD2011004_BCD_flagged.csv')
H11043_bcd <- read.csv('HUD2011043/HUD2011043_BCD_flagged.csv')
H12042_bcd <- read.csv('HUD2012042/HUD2012042_BCD_flagged.csv')

Once files were read into R, Winkler and probe data was found and pulled.

#find rows with winkler data
H11004_wink <- grep(H11004_bcd$DATA_TYPE_METHOD, pattern = 'O2_Winkler_Auto')
H11043_wink <- grep(H11043_bcd$DATA_TYPE_METHOD, pattern = 'O2_Winkler_Auto')
H12042_wink <- grep(H12042_bcd$DATA_TYPE_METHOD, pattern = 'O2_Winkler_Auto')
#find rows with probe data
H11004_elec <- grep(H11004_bcd$DATA_TYPE_METHOD, pattern = 'O2_Electrode_mll')
H11043_elec <- grep(H11043_bcd$DATA_TYPE_METHOD, pattern = 'O2_Electrode_mll')
H12042_elec <- grep(H12042_bcd$DATA_TYPE_METHOD, pattern = 'O2_Electrode_mll')
#pull winkler data
H11004_wink_data <- H11004_bcd$DIS_DETAIL_DATA_VALUE[H11004_wink]
H11004_wink_dpth <- H11004_bcd$DIS_HEADER_START_DEPTH[H11004_wink]
H11043_wink_data <- H11043_bcd$DIS_DETAIL_DATA_VALUE[H11043_wink]
H11043_wink_dpth <- H11043_bcd$DIS_HEADER_START_DEPTH[H11043_wink]
H12042_wink_data <- H12042_bcd$DIS_DETAIL_DATA_VALUE[H12042_wink]
H12042_wink_dpth <- H12042_bcd$DIS_HEADER_START_DEPTH[H12042_wink]
#pull probe data
H11004_elec_data <- H11004_bcd$DIS_DETAIL_DATA_VALUE[H11004_elec]
H11004_elec_dpth <- H11004_bcd$DIS_HEADER_START_DEPTH[H11004_elec]
H11043_elec_data <- H11043_bcd$DIS_DETAIL_DATA_VALUE[H11043_elec]
H11043_elec_dpth <- H11043_bcd$DIS_HEADER_START_DEPTH[H11043_elec]
H12042_elec_data <- H12042_bcd$DIS_DETAIL_DATA_VALUE[H12042_elec]
H12042_elec_dpth <- H12042_bcd$DIS_HEADER_START_DEPTH[H12042_elec]

Once probe and Winkler data has been isolated, an initial plot gives a general idea of the data set for each cruise.

plot(H11004_wink_data, H11004_wink_dpth, ylim = c(5000, 0), xlim = c(3, 9), col = 'blue', ylab = 'Depth (m)', xlab = 'O2 concentration (ml/l)',main = 'Oxygen data from HUD2011004', sub = 'Winkler vs Probe(Orion 842)')
par(new = TRUE)

plot(H11004_elec_data, H11004_elec_dpth, ylim = c(5000, 0), xlim = c(3, 9), axes = FALSE, col = 'black', xlab = '', ylab = '')
legend('bottomleft', legend = c('WINKLER', 'PROBE'), col = c('blue', 'black'), lty = 1, cex = 0.6)

plot(H11043_wink_data, H11043_wink_dpth, ylim = c(5000, 0), xlim = c(2, 9), col = 'blue', ylab = 'Depth (m)', xlab = 'O2 concentration (ml/l)',main = 'Oxygen data from HUD2011043', sub = 'Winkler vs Probe(Orion 842)')
par(new = TRUE)

plot(H11043_elec_data, H11043_elec_dpth, ylim = c(5000, 0), xlim = c(2, 9), axes = FALSE, col = 'black', xlab = '', ylab = '')
legend('bottomleft', legend = c('WINKLER', 'PROBE'), col = c('blue', 'black'), lty = 1, cex = 0.6)

plot(H12042_wink_data, H12042_wink_dpth, ylim = c(5000, 0), xlim = c(1, 9), col = 'blue', ylab = 'Depth (m)', xlab = 'O2 concentration (ml/l)',main = 'Oxygen data from HUD2012042', sub = 'Winkler vs Probe(Orion 842)')
par(new = TRUE)

plot(H12042_elec_data, H12042_elec_dpth, ylim = c(5000, 0), xlim = c(1, 9), axes = FALSE, col = 'black', xlab = '', ylab = '')
legend('bottomleft', legend = c('WINKLER', 'PROBE'), col = c('blue', 'black'), lty = 1, cex = 0.6)

The next step in analysis was to match the sample ID’s for each measurement method so that DO could be directly compared. There was a minor bug with the length of the vectors but this was easily solved.

#isolate individual bottles
#hud2011004
#pull sample ids
H11004_wink_id <- H11004_bcd$DIS_DETAIL_COLLECTOR_SAMP_ID[H11004_wink] #length = 137
H11004_elec_id <- H11004_bcd$DIS_DETAIL_COLLECTOR_SAMP_ID[H11004_elec] #length = 176
#These are different lengths meaning not all sample ids can match
length(H11004_wink_id) <- 176 #fix length issue
#combine lists
H11004_id <- cbind(H11004_wink_id, H11004_elec_id)
#compare lists
H11004_id_match <- H11004_wink_id == H11004_elec_id #useless
#most sample id's match but they are ordered differently, must compare 1 by 1
#create matching list wink:elec
match <- list()
for (i in 1:176){
  h <- H11004_wink_id[i] == H11004_elec_id
  m <- grep(h, pattern = TRUE)
  #m <- grep(pattern = H11004_wink_id[i], H11004_elec_id)
  if(length(m >0)){
  match[i] <- m
  }else{
    match[i] <- NULL
  }
}
#name list wink:elec
    #why is list only 175 items, should be 176, last is NULL?
names(match) <- c(1:175)
H11004_wink_elec <- list()
for (i in 1:175){
  if(!is.null(match[[i]])){
H11004_wink_elec[i] <- H11004_wink_data[[i]] - H11004_elec_data[[match[[i]]]]
  }else{
  H11004_wink_elec[i] <- NULL
}
}
#keeps subtracting one from length?
length(H11004_wink_elec) <- 175

Once individual sample ID’s could be compared it was possible to plot the difference between methods.

#plot differences (wink-elec)
plot(unlist(match), unlist(H11004_wink_elec), type = 'h', xlab = 'sample ID', ylab = 'Winkler - electrode value')
mtext('Difference between Winkler and Electrode values 
      (HUD2011004)')

Some summary statistics were also calculated at this point and presented in the report.

This process was then repeated for each cruise.

HUD2011043

#hud2011043
#pull sample ids
H11043_wink_id <- H11043_bcd$DIS_DETAIL_COLLECTOR_SAMP_ID[H11043_wink] #length = 310
H11043_elec_id <- H11043_bcd$DIS_DETAIL_COLLECTOR_SAMP_ID[H11043_elec] #length = 71
#These are different lengths meaning not all sample ids can match
length(H11043_elec_id) <- 310 #fix length issue
#combine lists
H11043_id <- cbind(H11043_wink_id, H11043_elec_id)
#compare lists
H11043_id_match <- H11043_wink_id == H11043_elec_id #useless
#most sample id's match but they are ordered differently, must compare 1 by 1
#create matching list wink:elec
match2 <- list()
for (i in 1:310){
  h <- H11043_wink_id[i] == H11043_elec_id
  m <- grep(h, pattern = TRUE)
  #m <- grep(pattern = H11004_wink_id[i], H11004_elec_id)
  if(length(m >0)){
    match2[i] <- m
  }else{
    match2[i] <- NULL
  }
}
#name list wink:elec
#why is list one short of full length?
names(match2) <- c(1:309)
H11043_wink_elec <- list()
for (i in 1:309){
  if(!is.null(match2[[i]])){
    H11043_wink_elec[i] <- H11043_wink_data[[i]] - H11043_elec_data[[match2[[i]]]]
  }else{
    H11043_wink_elec[i] <- NULL
  }
}
#keeps subtracting one from length?
length(H11043_wink_elec) <- 310
#plot differences (wink-elec)
plot(unlist(match2), unlist(H11043_wink_elec), type = 'h', xlab = 'sample ID', ylab = 'Winkler - electrode value')
mtext('Difference between Winkler and Electrode values 
      (HUD2011043)')

HUD2012042

#Hud2012042
#pull sample ids
H12042_wink_id <- H12042_bcd$DIS_DETAIL_COLLECTOR_SAMP_ID[H12042_wink] #length = 486
H12042_elec_id <- H12042_bcd$DIS_DETAIL_COLLECTOR_SAMP_ID[H12042_elec] #length = 207
#These are different lengths meaning not all sample ids can match
length(H12042_elec_id) <- 486 #fix length issue
#combine lists
H12042_id <- cbind(H12042_wink_id, H12042_elec_id)
#compare lists
H12042_id_match <- H12042_wink_id == H12042_elec_id #useless
#most sample id's match but they are ordered differently, must compare 1 by 1
#create matching list wink:elec
match3 <- list()
for (i in 1:486){
  h <- H12042_wink_id[i] == H12042_elec_id
  m <- grep(h, pattern = TRUE)
  #m <- grep(pattern = H11004_wink_id[i], H11004_elec_id)
  if(length(m >0)){
    match3[i] <- m
  }else{
    match3[i] <- NULL
  }
}
#name list wink:elec
#why is list only 175 items, should be 176, last is NULL?
names(match3) <- c(1:486)
H12042_wink_elec <- list()
for (i in 1:486){
  if(!is.null(match3[[i]])){
    H12042_wink_elec[i] <- H12042_wink_data[[i]] - H12042_elec_data[[match3[[i]]]]
  }else{
    H12042_wink_elec[i] <- NULL
  }
}
#keeps subtracting one from length?
length(H12042_wink_elec) <- 486
#plot differences (wink-elec)
plot(unlist(match3), unlist(H12042_wink_elec), type = 'h', xlab = 'sample ID', ylab = 'Winkler - electrode value')
mtext('Difference between Winkler and Electrode values 
      (HUD2012042)')

Further statistics were calculated including a t test and confidence intervals

##is the mean difference significantly diffrent from zero?
  #t test
a <- 0.01
H12042_t <- t.test(x = unlist(H12042_wink_elec))
H11043_t <- t.test(x = unlist(H11043_wink_elec))
H11004_t <- t.test(x = unlist(H11004_wink_elec))
ci1 <- H11004_t$conf.int
ci2 <- H11043_t$conf.int
ci3 <- H12042_t$conf.int
x <- c(1: 200)
a <- unlist(H11004_wink_elec)
b <- unlist(H11043_wink_elec)
c <- unlist(H12042_wink_elec)
length(a) <- 200
length(b) <- 200
length(c) <- 200
plot(x, a, col = 'black', ylim = c(-2.5, 2.5), xlab = 'Sample ID', ylab = "Difference Winkler - Electrode (ml/l)" )
mtext('Difference between 
      Winkler and Electrode Methods')

mtext("Black: HUD2011004,
  Blue: HUD2011043, 
  Red: HUD2012042", side = 1, cex = 0.8, line = -2)
par(new = TRUE)

plot(x, b, col = 'blue', ylim = c(-2.5, 2.5), axes = FALSE, xlab = '', ylab = '')
par(new = TRUE)

plot(x, c, col = 'red', ylim = c(-2.5, 2.5), axes = FALSE, xlab = '', ylab = '')
#confidence interval lines
abline(h = ci1[[1]], col = 'black')

abline(h = ci1[[2]], col = 'black')
abline(h = ci2[[1]], col = 'blue')

abline(h = ci2[[2]], col = 'blue')
abline(h = ci3[[1]], col = 'red')

abline(h = ci3[[2]], col = 'red')

The next step of analysis was to remove outliers in the Winkler data using Andrew Cogswell’s method of eliminating points outside 1.5*IQR

#eliminate outliers (past 1.5 IQT, and flagged values)
#calculate iqt
H11004_iqr<- IQR(H11004_wink_data, na.rm = TRUE)
H11043_iqr <- IQR(H11043_wink_data, na.rm = TRUE)
H12042_iqr <- IQR(H12042_wink_data, na.rm = TRUE)
#calculate 1.5 iqt
H11004_iqr_1.5 <- 1.5*(H11004_iqr)
H11043_iqr_1.5 <- 1.5*(H11043_iqr)
H12042_iqr_1.5 <- 1.5*(H12042_iqr)
#calculate mean
H11004_wink_mean <- mean(H11004_wink_data)
H11043_wink_mean <- mean(H11043_wink_data)
H12042_wink_mean <- mean(H12042_wink_data)

This was then plotted so data could be visualized and ensure correct removal of outliers

#hud2011004
plot( c(1:length(H11004_wink_data)), unlist(H11004_wink_data), xlab = 'Sample ID', ylab = 'Oxygen COncentration (ml/L)')
mtext('Dissolved Oxygen Winkler Measurements
      HUD2011004')

par(new = TRUE)
abline(h = H11004_wink_mean, col = 'blue')

uiqr <- 0.5*(H11004_iqr_1.5) + H11004_wink_mean
liqr <- H11004_wink_mean - 0.5*(H11004_iqr_1.5) 
abline(h = uiqr, col = 'red')
abline(h = liqr, col = 'red')

#null values outside 1.5 iqt
H11004_wink_data[H11004_wink_data > uiqr] <- NA
H11004_wink_data[H11004_wink_data < liqr] <- NA
plot(c(1:length(H11004_wink_data)), unlist(H11004_wink_data), ylim = c(3,9), xlab = 'Sample ID', ylab = 'Dissolved OXygen (ml/L)')
mtext('Winkler Oxygen Measurements after QC
      HUD2011004')

mtext('Red lines represent 1.5 IQT', side = 1, cex = 0.8)
abline(h = uiqr, col = 'red')

abline(h = liqr, col = 'red')

#hud2011043
plot( c(1:length(H11043_wink_data)), unlist(H11043_wink_data), xlab = 'Sample ID', ylab = 'Oxygen COncentration (ml/L)')
mtext('Dissolved Oxygen Winkler Measurements
      HUD2011043')

par(new = TRUE)
abline(h = H11043_wink_mean, col = 'blue')

uiqr2 <- 0.5*(H11043_iqr_1.5) + H11043_wink_mean
liqr2 <- H11043_wink_mean - 0.5*(H11043_iqr_1.5) 
abline(h = uiqr2, col = 'red')
abline(h = liqr2, col = 'red')

#null values outside 1.5 iqt
H11043_wink_data[H11043_wink_data > uiqr2] <- NA
H11043_wink_data[H11043_wink_data < liqr2] <- NA
plot(c(1:length(H11043_wink_data)), unlist(H11043_wink_data), ylim = c(3, 7), xlab = 'Sample ID', ylab = 'Dissolved OXygen (ml/L)')
mtext('Winkler Oxygen Measurements after QC
      HUD2011043')

mtext('Red lines represent 1.5 IQT', side = 1, cex = 0.8)
abline(h = uiqr2, col = 'red')

abline(h = liqr2, col = 'red')

#hud2012042
plot( c(1:length(H12042_wink_data)), unlist(H12042_wink_data), xlab = 'Sample ID', ylab = 'Oxygen Concentration (ml/L)')
mtext('Dissolved Oxygen Winkler Measurements
      HUD2012042')

par(new = TRUE)
abline(h = H12042_wink_mean, col = 'blue')

uiqr3 <- 0.5*(H12042_iqr_1.5) + H12042_wink_mean
liqr3 <- H12042_wink_mean - 0.5*(H12042_iqr_1.5) 
abline(h = uiqr3, col = 'red')
abline(h = liqr3, col = 'red')

#null values outside 1.5 iqt
H12042_wink_data[H12042_wink_data > uiqr3] <- NA
H12042_wink_data[H12042_wink_data < liqr3] <- NA
plot(c(1:length(H12042_wink_data)), unlist(H12042_wink_data), ylim = c(2,7), xlab = 'Sample ID', ylab = 'Dissolved Oxygen (ml/L)')
mtext('Winkler Oxygen Measurements after QC
      HUD2012042')

mtext('Red lines represent 1.5 IQT', side = 1, cex = 0.8)
abline(h = uiqr3, col = 'red')

abline(h = liqr3, col = 'red')

Differences within individual sample ID’s were then recalculated, removing any samples where the Winkler value was null due to being an outlier.

#plot differences after removal
#hud2011004
H11004_wink_elec_1 <- list()
for( i in 1:175){
  if(!is.null(match[[i]])){
    H11004_wink_elec_1[i] <- H11004_wink_data[[i]] - H11004_elec_data[[match[[i]]]]
  }else{
    H11004_wink_elec_1[i] <- NULL
  }
}
plot(unlist(match), unlist(H11004_wink_elec_1), type = 'h', xlab = 'Sample ID', ylab = 'Winkler - Electrode (DO , ml/L)')
mtext('Difference between Winkler and Electrode values (post QC)
      HUD2011004')

#hud2011043
H11043_wink_elec_1 <- list()
for( i in 1:309){
  if(!is.null(match2[[i]])){
    H11043_wink_elec_1[i] <- H11043_wink_data[[i]] - H11043_elec_data[[match2[[i]]]]
  }else{
    H11043_wink_elec_1[i] <- NULL
  }
}
plot(unlist(match2), unlist(H11043_wink_elec_1), type = 'h', xlab = 'Sample ID', ylab = 'Winkler - Electrode (DO , ml/L)')
mtext('Difference between Winkler and Electrode values (post QC)
      HUD2011043')

#hud2012042
H12042_wink_elec_1 <- list()
for( i in 1:486){
  if(!is.null(match3[[i]])){
    H12042_wink_elec_1[i] <- H12042_wink_data[[i]] - H12042_elec_data[[match3[[i]]]]
  }else{
    H12042_wink_elec_1[i] <- NULL
  }
}
plot(unlist(match3), unlist(H12042_wink_elec_1), type = 'h', xlab = 'Sample ID', ylab = 'Winkler - Electrode (DO , ml/L)')
mtext('Difference between Winkler and Electrode values (post QC)
      HUD2012042')

Depth profiles were then visualized again to see if there were any significant visual patterns or relationship between depth and offset between Winkler and probe data.

##post QC depth profiles
plot(H11004_wink_data, H11004_wink_dpth, ylim = (c(5000, 0)), xlim = c(4,9), xlab = 'DO ml/L', ylab = 'Depth (m)', main = 'HUD2011004')
par(new = TRUE)

plot(H11004_elec_data, H11004_elec_dpth, ylim = (c(5000, 0)),xlim = c(4,9), axes = FALSE, col = 'blue', xlab = '', ylab = '')
legend('bottomright',  c('Winkler', 'Electrode'), col = c('black', 'blue'), lty = 1)

plot(H11043_wink_data, H11043_wink_dpth, ylim = (c(5000, 0)), xlim = c(4,9), xlab = 'DO ml/L', ylab = 'Depth (m)', main = 'HUD2011043')
par(new = TRUE)

plot(H11043_elec_data, H11043_elec_dpth, ylim = (c(5000, 0)),xlim = c(4,9), axes = FALSE, col = 'blue', xlab = '', ylab = '')
legend('bottomright',  c('Winkler', 'Electrode'), col = c('black', 'blue'), lty = 1)

plot(H12042_wink_data, H12042_wink_dpth, ylim = (c(5000, 0)), xlim = c(4,9), xlab = 'DO ml/L', ylab = 'Depth (m)', main = 'HUD2012042')
par(new = TRUE)

plot(H12042_elec_data, H12042_elec_dpth, ylim = (c(5000, 0)),xlim = c(4,9), axes = FALSE, col = 'blue', xlab = '', ylab = '')
legend('bottomright',  c('Winkler', 'Electrode'), col = c('black', 'blue'), lty = 1, cex = 0.8)

The next step in Analysis was to plot Winkler and Electrode data against each other for each cruise to attempt a linear regression. In an ideal situtation the data would have presented with a 1:1 relationship, high correlation and minimal offset (intercept).

For cruise HUD2011004

H11004_winkelec <- matrix(NA, nrow = 175, ncol = 2)
for( i in 1:175){
  if(!is.null(match[[i]])){
    H11004_winkelec[i,1] <- H11004_wink_data[[i]] 
    H11004_winkelec[i,2] <- H11004_elec_data[[match[[i]]]]
  # }else{
  #   H11004_winkelec[i,] <- NULL
   }
}
plot(H11004_winkelec[,1], H11004_winkelec[,2], xlim = c(0,7), ylim = c(0,7),xlab = 'Winkler Data', ylab = 'Electrode Data')
par(new = TRUE)

lm <- lm(H11004_winkelec[,1] ~ H11004_winkelec[,2])
abline(lm)
abline(coef = c(0, 1), col = 'red')

title(main = 'HUD2011004')
title(sub = 'R^2 = 0.9844')

Where black line represents simple linear regression on the data and red line represents perfect 1:1 relationship with 0 intercept.

Statistic summary

summary(lm)

Call:
lm(formula = H11004_winkelec[, 1] ~ H11004_winkelec[, 2])

Residuals:
     Min       1Q   Median       3Q      Max 
-0.75153 -0.08761  0.01495  0.10951  0.43717 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          -1.30733    0.10351  -12.63   <2e-16 ***
H11004_winkelec[, 2]  1.27511    0.01592   80.10   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1696 on 102 degrees of freedom
  (71 observations deleted due to missingness)
Multiple R-squared:  0.9844,    Adjusted R-squared:  0.9842 
F-statistic:  6417 on 1 and 102 DF,  p-value: < 2.2e-16

plot(lm)

For cruise HUD2011043

H11043_winkelec <- matrix(NA, nrow = 310, ncol = 2)
for( i in 1:309){
  if(!is.null(match2[[i]])){
    H11043_winkelec[i,1] <- H11043_wink_data[[i]] 
    H11043_winkelec[i,2] <- H11043_elec_data[[match2[[i]]]]
  # }else{
  #   H11043_winkelec[i,] <- NULL
   }
}
plot(H11043_winkelec[,1], H11043_winkelec[,2], xlim = c(0, 7), ylim = c(0,7), xlab = 'Winkler Data', ylab = 'Electrode Data')
par(new = TRUE)

lm2 <- lm(H11043_winkelec[,1] ~ H11043_winkelec[,2])
abline(lm2)
abline(coef = c(0, 1), col = 'red')

title(main = 'HUD2011043')
title(sub = 'R^2 = 0.3697')

Where black line represents simple linear regression on the data and red line represents perfect 1:1 relationship with 0 intercept.

Statistic summary

summary(lm2)

Call:
lm(formula = H11043_winkelec[, 1] ~ H11043_winkelec[, 2])

Residuals:
     Min       1Q   Median       3Q      Max 
-0.77188 -0.15051  0.04768  0.19181  0.49000 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)            2.8913     0.6757   4.279 0.000159 ***
H11043_winkelec[, 2]   0.5370     0.1240   4.332 0.000136 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2698 on 32 degrees of freedom
  (276 observations deleted due to missingness)
Multiple R-squared:  0.3697,    Adjusted R-squared:   0.35 
F-statistic: 18.77 on 1 and 32 DF,  p-value: 0.0001365

plot(lm2)

For cruise HUD2012042

length(H12042_elec_data) <- 486
H12042_winkelec <- matrix(NA, nrow = 486, ncol = 2)
for( i in 1:486){
  if(!is.null(match3[[i]])){
    H12042_winkelec[i,1] <- H12042_wink_data[[i]] 
    H12042_winkelec[i,2] <- H12042_elec_data[[match3[[i]]]]
   # }else{
   #   H12042_winkelec[i,] <- NULL
   }
}
plot(H12042_winkelec[,1], H12042_winkelec[,2],xlim = c(0,7), ylim = c(0,7), xlab = 'Winkler Data', ylab = 'Electrode Data')
par(new = TRUE)

lm3 <- lm(H12042_winkelec[,1] ~ H12042_winkelec[,2])
abline(lm3)
abline(coef = c(0, 1), col = 'red')

title(main = 'HUD2012042')
title(sub = 'R^2 = 0.4968')

Where black line represents simple linear regression on the data and red line represents perfect 1:1 relationship with 0 intercept.

Statistic summary

summary(lm3)

Call:
lm(formula = H12042_winkelec[, 1] ~ H12042_winkelec[, 2])

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3420 -0.1528  0.1133  0.2371  0.5601 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           0.89119    0.46932   1.899   0.0607 .  
H12042_winkelec[, 2]  0.82843    0.08693   9.530 2.22e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3978 on 92 degrees of freedom
  (392 observations deleted due to missingness)
Multiple R-squared:  0.4968,    Adjusted R-squared:  0.4913 
F-statistic: 90.83 on 1 and 92 DF,  p-value: 2.221e-15

plot(lm3)

In order to determine if the relationship between Winkler and electrode data was pressure dependent the differences by cruise were plotted over depth along with simple linear regression and summary statistics

For HUD2011004

length(H11004_wink_dpth) <- 131
plot( H11004_wink_dpth, unlist(H11004_wink_elec_1), ylim = c(-2, 2), xlab = 'Depth (m)', ylab = 'Difference (Winkler - electrode)')
abline(h = 0, col = 'blue')

ll <- lm( unlist(H11004_wink_elec_1) ~ H11004_wink_dpth)
abline(ll, col = 'red')
title(main = 'HUD2011004', sub = 'R^2 = 0.01128')

summary(ll)

Call:
lm(formula = unlist(H11004_wink_elec_1) ~ H11004_wink_dpth)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.88809 -0.21947  0.09733  0.25251  0.70152 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.405e-01  3.675e-02  11.988   <2e-16 ***
H11004_wink_dpth 3.770e-05  3.494e-05   1.079    0.283    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3342 on 102 degrees of freedom
  (27 observations deleted due to missingness)
Multiple R-squared:  0.01128,   Adjusted R-squared:  0.00159 
F-statistic: 1.164 on 1 and 102 DF,  p-value: 0.2832

plot(ll)

For HUD2011043

dp2 <- list()
for (i in 1:length(H11043_wink_elec_1)){
  if(is.null(H11043_wink_elec_1[[i]])){
    dp2[i] <- NULL
  } else{
    dp2[i] <- H11043_wink_dpth[[i]]
  }
}
plot(unlist(dp2), unlist(H11043_wink_elec_1), ylim = c(-2, 2), xlab = 'Depth (m)', ylab = 'Difference (Winkler - Electrode)')
abline(h = 0, col = 'blue')

ll2 <- lm(unlist(H11043_wink_elec_1) ~ unlist(dp2))
abline(ll2, col = 'red')
title(main = 'HUD2011043', sub = 'R^2 = 0.00944')

summary(ll2)

Call:
lm(formula = unlist(H11043_wink_elec_1) ~ unlist(dp2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.84757 -0.12671  0.06725  0.16245  0.56277 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.893e-01  6.233e-02   6.245 5.35e-07 ***
unlist(dp2) -2.088e-05  3.780e-05  -0.552    0.585    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3218 on 32 degrees of freedom
  (22 observations deleted due to missingness)
Multiple R-squared:  0.009444,  Adjusted R-squared:  -0.02151 
F-statistic: 0.3051 on 1 and 32 DF,  p-value: 0.5845

plot(ll2)

For HUD2012042

dp3 <- list()
for (i in 1:length(H12042_wink_elec_1)){
  if(is.null(H12042_wink_elec_1[[i]])){
    dp3[i] <- NULL
  } else{
    dp3[i] <- H12042_wink_dpth[[i]]
  }
}
plot(unlist(dp3), unlist(H12042_wink_elec_1), ylim = c(-2, 2), xlab = 'Depth (m)', ylab = 'Difference (Winkler - Electrode)')
abline(h = 0, col = 'blue')

ll3 <- lm(unlist(H12042_wink_elec_1) ~ unlist(dp3))
abline(ll3, col = 'red')
title(main = 'HUD2012042', sub = 'R^2 = 0.7767')

summary(ll3)

Call:
lm(formula = unlist(H12042_wink_elec_1) ~ unlist(dp3))

Residuals:
     Min       1Q   Median       3Q      Max 
-1.33694 -0.07307  0.13191  0.23863  0.49426 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -7.746e-02  4.348e-02  -1.782  0.07812 . 
unlist(dp3)  9.929e-05  3.567e-05   2.783  0.00653 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3901 on 92 degrees of freedom
  (72 observations deleted due to missingness)
Multiple R-squared:  0.07767,   Adjusted R-squared:  0.06765 
F-statistic: 7.748 on 1 and 92 DF,  p-value: 0.006527

plot(ll3)

This procedure was done to mimic CTD QC procedure but may not have signifcant bearing on the statistical results. In theory there should be no relationship between pressure and Winkler - electrode measurements since both measurements are made a t the surface, although there may be some relationship between accuracy/ offset between methods and dissolved oxygen concentration which will be tested next.

Difference (winkler - electrode measurement per sample ID) was plotted over DO concentration for each cruise, blue line are at y = 0 and red lines represent simple linear regression of the data.

For HUD2011004

wd <- list()
for (i in 1:length(H11004_wink_data)){
  if (is.null(H11004_wink_elec_1[[i]])){
    wd[[i]] <- NULL
  } else{
    wd[[i]] <- H11004_wink_data[[i]]
  }
}
length(H11004_wink_elec_1) <- length(wd)
plot(unlist(wd), unlist(H11004_wink_elec_1), ylim = c(-3,3), xlab = 'DO Concentrations, ml/l', ylab = 'Difference (Winkler - electrode)')
abline(h = 0, col = 'blue')

lh <- lm(unlist(H11004_wink_elec_1) ~ unlist(wd))
abline(lh, col = 'red')
title(main = 'HUD2011004')

summary(lh)

Call:
lm(formula = unlist(H11004_wink_elec_1) ~ unlist(wd))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.57249 -0.07691  0.01129  0.09296  0.32293 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.109656   0.067524  -16.43   <2e-16 ***
unlist(wd)   0.228027   0.009637   23.66   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1319 on 102 degrees of freedom
  (27 observations deleted due to missingness)
Multiple R-squared:  0.8459,    Adjusted R-squared:  0.8444 
F-statistic: 559.9 on 1 and 102 DF,  p-value: < 2.2e-16

plot(lh)

For HUD2011043

wd2 <- list()
for (i in 1:length(H11043_wink_elec_1)){
  if (is.null(H11043_wink_elec_1[[i]])){
    wd2[[i]] <- NULL
  } else{
    wd2[[i]] <- H11043_wink_data[[i]]
  }
}
wd2 <- unlist(wd2)
length(wd2) <- 56
plot(wd2, unlist(H11043_wink_elec_1), ylim = c(-3,3), xlab = 'DO Concentrations, ml/l', ylab = 'Difference (Winkler - electrode)')
abline(h = 0, col = 'blue')

lh2 <- lm(unlist(H11043_wink_elec_1) ~ wd2)
abline(lh2, col = 'red')
title(main = 'HUD2011043')

summary(lh2)

Call:
lm(formula = unlist(H11043_wink_elec_1) ~ wd2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.95656 -0.06824  0.02838  0.13253  0.53755 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -1.4374     0.9250  -1.554   0.1301  
wd2           0.3115     0.1589   1.961   0.0587 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3055 on 32 degrees of freedom
  (22 observations deleted due to missingness)
Multiple R-squared:  0.1072,    Adjusted R-squared:  0.07933 
F-statistic: 3.844 on 1 and 32 DF,  p-value: 0.05869

plot(lh2)

For HUD2012042

wd3 <- list()
for (i in 1:length(H12042_wink_elec_1)){
  if (is.null(H12042_wink_elec_1[[i]])){
    wd3[[i]] <- NULL
    }
  else{
    wd3[[i]] <- H12042_wink_data[[i]]
  }
}
plot(unlist(wd3), unlist(H12042_wink_elec_1), ylim = c(-3,3), xlab = 'DO Concentrations, ml/l', ylab = 'Difference (Winkler - electrode)')
abline(h = 0, col = 'blue')

lh3 <- lm(unlist(H12042_wink_elec_1) ~ unlist(wd3))
abline(lh3, col = 'red')
title(main = 'HUD2012042')

summary(lh3)

Call:
lm(formula = unlist(H12042_wink_elec_1) ~ unlist(wd3))

Residuals:
     Min       1Q   Median       3Q      Max 
-1.29387 -0.08654  0.05764  0.21150  0.52617 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.17202    0.33825  -6.421 5.83e-09 ***
unlist(wd3)  0.40032    0.06292   6.362 7.63e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3385 on 92 degrees of freedom
  (72 observations deleted due to missingness)
Multiple R-squared:  0.3055,    Adjusted R-squared:  0.298 
F-statistic: 40.48 on 1 and 92 DF,  p-value: 7.628e-09

plot(lh3)

These plots show some interesting trends of accuracy over DO concentrations.

The trend seems to be a slightly increased under estimation by the electrode at higher DO concentrations. Keeping in mind the small range of DO concentrations that are visible in this data this may be a significant trend.

Groundfish surveys were also analyzed in order to bulk up the data set. A total of 55 new data points were discovered where winkler and lectrode measurments were performed for the same sample ID.

#read in groundfish data
data <- readxl::read_xlsx('C:/Users/ChisholmE/Desktop/BIOCHEM/GROUNDFISH OXYGEN DATA.xlsx' )
#isolate by cruise
NED2013022 <- (5:935)
NED2014018 <- (937:1669)
NED2015017 <- (1674:2161) #CAUTIOUS OF DUPLICATE VALUE, x__4 IS ELECTRODE VALUE AVERAGED BETWEEN REPS
#read data column containing difference (Winkler- electrode)
DIFF <- as.numeric(as.character(data$X__6[-1]))
diff <- na.omit(DIFF)
#plot differences over both cruises due to minimal data points
plot(diff, type = 'h', xlab = 'Sample ID', ylab = 'Difference between methods (Winkler - Electrode)', main = 'Difference in O2 measurments from Groundfish surveys')
m <- mean(diff)
abline(h = m, col = 'red')

Plotting these difference over dissolved oxygen concetrations revealed trends consistent with AZMP data

#isolate concentrations from Winkler measurments for relevant sample IDs
conc <- list()
for ( i in 1:2648){
  if(is.na(data$X__6[-1][i]) ){
    conc[i] <- NA
  }else{
    conc[i] <- as.numeric(data$X__5[-1][i]) 
  }
}
  
c <- na.omit(unlist(conc))
#conc <- data$X__5[!is.na(DIFF) ]
#plot differences over concentration
plot(conc, DIFF, xlab = 'Dissolved Oxygen Concentration (ml/L)', ylab = 'Difference between measurements (Winkler - electrode)')
abline(h = 0, col = 'blue')

#linear regression
cm <- lm(unlist(DIFF) ~unlist(conc))
abline(cm, col = 'red')
#plot potential correction factors 
#of 2, 4, 6 %
twoperc <- list()
for (i in 1:length(conc)){
  twoperc[i] <- as.numeric(conc[i]) *0.02
}
pm <- lm(unlist(twoperc)~unlist(conc))
abline(pm, col = 'green')

fourperc <- list()
for (i in 1:length(conc)){
  fourperc[i] <- as.numeric(conc[i]) *0.04
}
fm <- lm(unlist(fourperc)~unlist(conc))
abline(fm, col = 'orange')
sixperc <- list()
for (i in 1:length(conc)){
  sixperc[i] <- as.numeric(conc[i]) *0.06
}
sm <- lm(unlist(sixperc)~unlist(conc))
abline(sm, col = 'purple')

legend('bottomright', legend = c('2%', '4%', '6%'), col = c('green', 'orange', 'purple'), lty = 2, cex = 0.5 )

Potential correction factors of 2, 4,6 % of dissolved oxygen concentrations were also plotted. Through visual inspection it appears the 4% correction factors agrees best with the data although there is still signifcant unexplained variability particularly at higher DO concentration values which may indicate a non linear trend.

Regression was performed over groundfish surveys and AZMP cruies to find a common trend in the total data set.

#plot original data
plot(wt, et, xlim = c(0, 10), ylim = c(0,10), col = 'red', xlab = 'Winkler (ml/L)', ylab = 'Electrode (ml/L)')
#regression of original data
ttlm <- lm(wt ~as.numeric(et))
abline(ttlm, col = 'red', lty = 2)

#plot 1:1 line reference
abline(coef = c(0,1), col = 'blue')
#apply correction factor following ttlm regression equation
eet <- list()
for (i in 1:length(et)){
  eet[i] <- as.numeric(et[i])*1.1424 - 0.6171
}
#plot corrected data for comparision
par(new = TRUE)

plot(wt, eet, xlim = c(0,10), ylim = c(0,10), xlab = '', ylab = '', axes = FALSE, col = 'blue')
mtext(side = 3, 'Original data (red), Corrected data (blue)')

#regress corrected data to ensure high correlation and 1.00 slope
ttcor <- lm(wt ~ unlist(eet))

Differences were analyzed before and after correction

#calculate original and corrected differences
wt_et <- list()
for (i in 1:length(wt)){
  wt_et[i] <- as.numeric(wt)[i] - as.numeric(et)[i]
}
wt_eet <- list()
for (i in 1:length(wt)){
  wt_eet[i] <- as.numeric(wt[i]) - as.numeric(eet[i])
}
#plot differences
plot(unlist(wt_et), type = 'h', ylim = c(-1,1), col = 'red', xlab = 'Sample ID', ylab = 'Difference (Winkler - electrode)')
par(new = TRUE)

plot(unlist(wt_eet), type = 'h', ylim = c(-1,1), axes = FALSE, xlab = '', ylab = '')
etm <- mean(unlist(wt_et), na.rm = TRUE)
eetm <- mean(unlist(wt_eet), na.rm = TRUE)
abline(h = etm, col = 'red', lty = 2)

abline(h = eetm, col = 'black', lty = 2)
mtext(side = 3, 'Original data differences (red), corrected data differences (black)')

Then differences over dissolved oxygen concentration were analyzed to ensure that there was random scatter around zero rather than increasing trend observed in original data

#plot new corrected data over DO concentration
plot(eet, wt_eet, xlim = c(2,9), ylim = c(-1,1), ylab = 'Differences (Winkler - Electrode**corrected)', xlab = "Dissolved oxygen concentration (ml/L)")

Gordana suggested plotting differences between WInkler and Electrode measurements over lab - in situ temparture differences to see if the temperature difference had a significant effect on the dissolved oxygen readings.

x <- readxl::read_xlsx('C:/Users/ChisholmE/Desktop/BIOCHEM/hud2011004_tdiff_exercise.xlsx')
w_e_diff <- x$`W - E difference`
t_diff <- x$`lab - in situ temp (mean)`
plot(t_diff, w_e_diff, xlab = 'Temperature difference (lab - in situ) degC', ylab = 'Oxygen difference (Winkler - Electrode) ml/l')
abline(h = 0, lty = 2)

This plot shows no significant trends in dissolved oxygen concentration over temperature differences.The points are mostly scattered around zero with some outliers at lower temperature differences, as well as some low values between 3 and 5 deg C temperature differences.

Plotting differences in oxygen over lab temperature was inconclusive. Samples are from HUD2011004

x <- readxl::read_xlsx('C:/Users/ChisholmE/Desktop/BIOCHEM/hud2011004_tdiff_exercise.xlsx' )
w_e_diff <- x$`W - E difference`
plot(x$X__11, w_e_diff, xlab = 'Temperature (degC)', ylab = 'Oxygen Difference (Winkler - electrode) ml/l')

NAs introduced by coercion