Difference between revisions of "EU-kalat"

# This code is Op_en on page [[EU-kalat]].

library(OpasnetUtils)
library(ggplot2)

eu <- Ovariable("eu", ddata = "Op_en3104", subset = "POPs")
eu <- EvalOutput(eu)
#levels(eu$POP)
eu <- eu[eu$Fish_species == "Baltic herring" , ]

tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values")
tef <- EvalOutput(tef)
colnames(eu@output)[colnames(eu@output) == "POP"] <- "Congener"
levels(eu$Congener) <- gsub("HCDD", "HxCDD", levels(eu$Congener))
levels(eu$Congener) <- gsub("HCDF", "HxCDF", levels(eu$Congener))
levels(eu$Congener) <- gsub("CoPCB", "PCB", levels(eu$Congener))
#levels(eu$Congener)

euteq <- eu * tef
#head(euteq@output)
temp <- aggregate(euteq@output["Result"], by = c(euteq@output["Congener"], euteq@output["Catch_location"]), FUN = mean)

temp2 <- temp
temp2$Result <- temp2$Result / sum(temp2$Result)

cat("Congener TEQ profile of Baltic herring.\n")
oprint(temp2)
 
ggplot(temp, aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") + facet_wrap(~ Catch_location) + 
labs(x = "Congener", y = "TEQ concentration (pg/g fat)") + theme_gray(base_size = 18) + coord_flip()

# This is code Op_en3104/preprocess on page [[EU-kalat]]
library(rjags)
library(OpasnetUtils)
library(ggplot2)
library(reshape2)

eu <- Ovariable("eu", ddata = "Op_en3104", subset = "POPs")
eu <- EvalOutput(eu)
eu@output <- eu@output[c(1:4, 18, 19)] # THL code, Matrix, Congener, Fish species
eu@marginal <- eu@marginal[c(1:4, 18, 19)]
colnames(eu@output) <- c("THLcode", "Matrix", "Congener", "Fish", "euResult", "euSource")

#> unique(eu@output$Congener)
#[1] 2378TCDD     12378PeCDD   123478HCDD   123678HCDD   123789HCDD   1234678HpCDD
#[7] OCDD         2378TCDF     12378PeCDF   23478PeCDF   123478HCDF   123678HCDF  
#[13] 123789HCDF   234678HCDF   1234678HpCDF 1234789HpCDF OCDF         CoPCB77 ...     

# Remove the four with too little data (>70% BDL) and all non-PCDDF
# aggregate(eu@data$euResult, by = eu@data["POP"], FUN = function(x) mean(x == 0))
conl <- unique(eu@output$Congener)[c(1:12, 14, 15)] # 7 OCDD should be removed
eu@output <- eu@output[eu@output$Congener %in% conl , ] 

#[1] Baltic herring Sprat          Salmon         Sea trout      Vendace       
#[6] Roach          Perch          Pike           Pike-perch     Burbot        
#[11] Whitefish      Flounder       Bream          River lamprey  Cod           
#[16] Trout          Rainbow trout  Arctic char   

fisl <- unique(eu@output$Fish)[c(1:4, 6:14, 17)] 
eu@output <- eu@output[eu@output$Fish %in% fisl , ] # Remove four with too little data

eut <- eu
eut@output <- eut@output[
  eut@output$Congener == "2378TCDD" & eut@output$Matrix == "Muscle" & result(eut) != 0 , ] # Zeros cannot be used in ratio estimates
eut@output$Congener <- NULL
eut@marginal[colnames(eut@output) == "Congener"] <- NULL
eut <- log10(eu / eut)@output
conl <- conl[conl != "2378TCDD"]
eut <- eut[eut$Congener %in% conl , ]

# Analysis: a few rows disappear here, as shown by numbers per fish species. Why?
# aggregate(eut@output, by = eut@output["Fish"], FUN = length)[[2]]
# [1] 1181  141  131   55  158   51   96   30   36   40  102   49   61  180 eu
# [1] 1173  141  131   55  158   51   96   27   36   40  102   49   53  180 eu/eut
# There are samples where TCDD concenctration is zero.
#eu@data[eu@data$POP == "2378TCDD" & eu@data$euResult == 0 , c(1,4)][[1]]
#eu@data[eu@data[[1]] %in% c("09K0585", "09K0698", "09K0740", "09K0748") , c(1,3,4,18)]
# Some rows disappear because there is no TCDD measurement to compare with:
#8 Baltic herrings, 8 sprats, 3 rainbow trouts.
# Conclusion: this is ok. Total 2292 rows.

ggplot(eu@output, aes(x = euResult, colour = Fish))+geom_density()+
  facet_wrap(~ Congener) + scale_x_log10()

ggplot(eut, aes(x = Result, colour = Fish))+geom_density()+
  facet_wrap(~ Congener) 

objects.store(eu,eut)
cat("Ovariable eu and data.frame eut stored.\n")