EU-kalat

This page is a study. The page identifier is Op_en3104
Moderator:Arja (see all)

Upload data {{#opasnet_base_link:Op_en3104}}

EU-kalat is a study, where concentrations of PCDD/Fs, PCBs, PBDEs and heavy metals have been measured from fish

Question

The scope of EU-kalat study was to measure concentrations of persistent organic pollutants (POPs) including dioxin (PCDD/F), PCB and BDE in fish from Baltic sea and Finnish inland lakes and rivers. ^[1] ^[2] ^[3].

Answer

The original sample results can be acquired from Opasnet base. The study showed that levels of PCDD/Fs and PCBs depends especially on the fish species. Highest levels were on salmon and large sized herring. Levels of PCDD/Fs exceeded maximum level of 4 pg TEQ/g fw multiple times. Levels of PCDD/Fs were correlated positively with age of the fish.

Mean congener concentrations as WHO2005-TEQ in Baltic herring can be printed out with the Run code below.

+ Show code - Hide code

# 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()

Rationale

Data

Data was collected between 2009-2010. The study contains years, tissue type, fish species, and fat content for each concentration measurement. Number of observations is 285.

There is a new study EU-kalat 3, which will produce results in 2016.

Calculations

Preprocess model 22.2.2017 [4]
- Objects used in Benefit-risk assessment of Baltic herring and salmon intake
Model run 25.1.2017 [5]
Model run 22.5.2017 with new ovariables euRaw, euAll, euMain, and euRatio [6]

+ Show code - Hide code

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

euRaw <- Ovariable("euRaw", ddata = "Op_en3104", subset = "POPs")

# levels(TEF$Group)
#[1] "Chlorinated dibenzo-p-dioxins" "Chlorinated dibenzofurans"    
#[3] "Mono-ortho–substituted PCBs"   "Non-ortho–substituted PCBs"   

euAll <- Ovariable(
  "euAll",
  dependencies = data.frame(
    Name=c("euRaw", "TEF"),
    Ident=c(NA,"Op_en4017/initiate")
  ),
  formula = function(...) {
    euAll <- euRaw[,c(1:4, 18, 19)] # THL code, Matrix, Congener, Fish species
    colnames(euAll@output)[1:4] <- c("THLcode", "Matrix", "Compound", "Fish")
    
    temp <- oapply(euAll * TEF, cols = "Compound", FUN = "sum")
    colnames(temp@output)[colnames(temp@output)=="Group"] <- "Compound"
    eu2 <- combine(euAll, temp)
    
    return(eu2)
  }
)

# Leave only the main fish species and congeners and remove others

euMain <- Ovariable(
  "euMain",
  dependencies = data.frame(
    Name = c("euAll", "conl", "fisl")
  ),
  formula = function(...) {
    euMain <- euAll[euAll$Compound %in% conl & euAll$Fish %in% fisl , ] 
    return(euMain)
  }
)

euRatio <- Ovariable(
  "euRatio",
  dependencies = data.frame(Name="euMain"),
  formula = function(...) {
    euRatio <- euMain[
      euMain$Compound == "2378TCDD" & euMain$Matrix == "Muscle" & result(euMain) != 0 , ] # Zeros cannot be used in ratio estimates
    euRatio$Compound <- NULL
    euRatio <- log10(euMain / euRatio)@output
    euRatio <- euRatio[euRatio$Compound %in% setdiff(conl, "2378TCDD") , ]
    return(euRatio)
  }
)

# 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.

objects.store(euRaw, euAll, euMain, euRatio)
cat("Ovariables euRaw, euAll, euMain, and euRatio stored.\n")

Bayes model for dioxin concentrations

Model run 28.2.2017 [7]
Model run 28.2.2017 with corrected survey model [8]
Model run 28.2.2017 with Mu estimates [9]
Model run 1.3.2017 [10]
Model run 23.4.2017 [11] produces list conc.param and ovariable concentration
Model run 24.4.2017 [12]
Model run 19.5.2017 without ovariable concentration [13] ⇤--#: . The model does not mix well, so the results should not be used for final results. --Jouni (talk) 19:37, 19 May 2017 (UTC) (type: truth; paradigms: science: attack)

----#: . Maybe we should just estimate TEQs until the problem is fixed. --Jouni (talk) 19:37, 19 May 2017 (UTC) (type: truth; paradigms: science: comment)

Model run 22.5.2017 with TEQdx and TEQpcb as the only Compounds [14]

+ Show code - Hide code

# This is code Op_en3104/bayes on page [[EU-kalat]]

library(OpasnetUtils)
library(reshape2)
library(rjags) # JAGS
library(ggplot2)
library(MASS) # mvrnorm
library(car) # scatterplotMatrix

objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]

################## Data for the main congeners and species only

#> unique(euAll$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))

#[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   

if(TRUE) {
  conl <- c("TEQdx", "TEQpcb")
} else {
  conl <- as.character(unique(euRaw@data$POP)[c(1:12, 14, 15)]) # 7 OCDD should be removed
}

fisl <- as.character(unique(euRaw@data$Fish_species)[c(1:4, 6:14, 17)])

# conl
#[1] "2378TCDD"     "12378PeCDD"   "123478HCDD"   "123678HCDD"   "123789HCDD"  
#[6] "1234678HpCDD" "OCDD"         "2378TCDF"     "12378PeCDF"   "23478PeCDF"  
#[11] "123478HCDF"   "123678HCDF"   "234678HCDF"   "1234678HpCDF"
#> fisl
#[1] "Baltic herring" "Sprat"          "Salmon"         "Sea trout"     
#[5] "Roach"          "Perch"          "Pike"           "Pike-perch"    
#[9] "Burbot"         "Whitefish"      "Flounder"       "Bream"         
#[13] "River lamprey"  "Rainbow trout" 

C <- length(conl)
Fi <- length(fisl)
N <- 1000
conl
fisl

euAll <- EvalOutput(euAll)

replaces <- list(
  c("Chlorinated dibenzo-p-dioxins", "TEQdx"),
  c("Chlorinated dibenzofurans", "TEQdx"),
  c("Mono-ortho–substituted PCBs", "TEQpcb"),
  c("Non-ortho–substituted PCBs", "TEQpcb")
)

for(i in 1:length(replaces)) {
  levels(euAll$Compound)[replaces[[i]][1]] <- replaces[[i]][2]
}

euRatio <- EvalOutput(euRatio)

# Hierarchical Bayes model.

# PCDD/F concentrations in fish.
# It uses the sum of PCDD/F (Pcdsum) as the total concentration of dioxin in fish.
# Cong_j is the fraction of a congener from pcdsum.
# pcdsum is log-normally distributed. cong_j follows Dirichlet distribution.
# pcdsum depends on age of fish, fish species and catchment area, but we only have species now so other variables are omitted.
# cong_j depends on fish species.

if(FALSE) { # Multicongener model NOT updated to reflect new names: eu -> euMain
  conl <- as.character(unique(eu@output$Congener))
  fisl <- sort(as.character(unique(eu@output$Fish)))
  C <- length(conl)
  Fi <- length(fisl)
  N <- 1000
  conl
  fisl
  fishsamples <- reshape(
    eu@output, 
    v.names = "euResult", 
    idvar = "THLcode", 
    timevar = "Congener", 
    drop = c("Matrix", "euSource"), 
    direction = "wide"
  )
  
  # Find the level of quantification for dinterval function
  LOQ <- unlist(lapply(fishsamples[3:ncol(fishsamples)], FUN = function(x) min(x[x!=0])))
  names(LOQ) <- conl
  cong <- data.matrix(fishsamples[3:ncol(fishsamples)])
  cong <- sapply(1:length(LOQ), FUN = function(x) ifelse(cong[,x]==0, 0.5*LOQ[x], cong[,x]))
  
  mod <- textConnection("
                        model{
                        for(i in 1:S) { # s = fish sample
                        #        below.LOQ[i,j] ~ dinterval(-cong[i,j], -LOQ[j])
                        cong[i,1:C] ~ dmnorm(mu[fis[i],], Omega[fis[i],,])
                        }
                        for(i in 1:Fi) { # Fi = fish species
                        for(j in 1:C) { 
                        mu[i,j] ~ dnorm(mu1[j], tau1[j])
                        }
                        Omega[i,1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
                        pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction
                        }
                        for(i in 1:C) { # C = Compound
                        mu1[i] ~ dnorm(0, 0.0001)
                        tau1[i] ~ dunif(0,10000)
                        pred1[i] ~ dnorm(mu1[i], tau1[i])
                        }
                        Omega1[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
                        }
                        ")
  
  jags <- jags.model(
    mod,
    data = list(
      S = nrow(fishsamples),
      C = C,
      Fi = Fi,
      cong = log(cong),
      #    LOQ = LOQ,
      #    below.LOQ = is.na(cong)*1,
      fis = match(fishsamples$Fish, fisl),
      Omega0 = diag(C)/100000
    ),
    n.chains = 4,
    n.adapt = 100
  )
} # if(FALSE)

fishsamples <- unkeep(euAll, prevresults = TRUE, sources = TRUE)@output
fishsamples <- fishsamples[fishsamples$Compound %in% conl & fishsamples$Matrix == "Muscle" , ]
fishsamples <- reshape(
  fishsamples, 
  v.names = "euAllResult", 
  idvar = "THLcode", 
  timevar = "Compound", 
  drop = c("Matrix"), 
  direction = "wide"
)

# colnames(euAll@output)
# "THLcode"      "Matrix"       "Compound"     "Fish"         "euRawSource" 
# [6] "TEFversion"   "TEFrawSource" "TEFSource"    "Source"       "euAllResult" 
# [11] "euAllSource"

oprint(head(fishsamples))

# Find the level of quantification for dinterval function
LOQ <- unlist(lapply(fishsamples[3:ncol(fishsamples)], FUN = function(x) min(x[x!=0])))
names(LOQ) <- conl
cong <- data.matrix(fishsamples[3:ncol(fishsamples)])
cong <- sapply(1:length(LOQ), FUN = function(x) ifelse(cong[,x]==0, 0.5*LOQ[x], cong[,x]))

mod <- textConnection("
                      model{
                      for(i in 1:S) { # s = fish sample
                      #        below.LOQ[i,j] ~ dinterval(-cong[i,j], -LOQ[j])
                      cong[i,1:C] ~ dmnorm(mu[fis[i],], Omega[fis[i],,])
                      }
                      for(i in 1:Fi) { # Fi = fish species
                      for(j in 1:C) { 
                      mu[i,j] ~ dnorm(mu1[j], tau1[j])
                      }
                      Omega[i,1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
                      pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction
                      }
                      for(i in 1:C) { # C = Compound
                      mu1[i] ~ dnorm(0, 0.0001)
                      tau1[i] ~ dunif(0,10000)
                      pred1[i] ~ dnorm(mu1[i], tau1[i])
                      }
                      Omega1[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
                      }
                      ")

jags <- jags.model(
  mod,
  data = list(
    S = nrow(fishsamples),
    C = C,
    Fi = Fi,
    cong = log(cong),
    #    LOQ = LOQ,
    #    below.LOQ = is.na(cong)*1,
    fis = match(fishsamples$Fish, fisl),
    Omega0 = diag(C)/100000
  ),
  n.chains = 4,
  n.adapt = 100
)

update(jags, 100)

samps.j <- jags.samples(
  jags, 
  c('mu', 'Omega', 'pred', 'mu1', 'Omega1', 'tau1', 'pred1'), 
  N
)
dimnames(samps.j$mu) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$mu1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$pred) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$mu1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$tau1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$pred1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$Omega) <- list(Fish = fisl, Compound = conl, Compound2 = conl, Iter=1:N, Chain=1:4)
dimnames(samps.j$Omega1) <- list(Compound = conl, Compound2 = conl, Iter=1:N, Chain=1:4)

##### conc.param contains expected values of the distribution parameters from the model
conc.param <- list(
  mu = apply(samps.j$mu[,,,1], MARGIN = 1:2, FUN = mean),
  Omega = apply(samps.j$Omega[,,,,1], MARGIN = 1:3, FUN = mean),
  pred.mean = apply(samps.j$pred[,,,1], MARGIN = 1:2, FUN = mean),
  pred.sd = apply(samps.j$pred[,,,1], MARGIN = 1:2, FUN = sd),
  mu1 = apply(samps.j$mu1[,,1], MARGIN = 1, FUN = mean),
  tau1 = apply(samps.j$tau1[,,1], MARGIN = 1, FUN = mean),
  pred1.mean = apply(samps.j$pred1[,,1], MARGIN = 1, FUN = mean),
  pred1.sd = apply(samps.j$pred1[,,1], MARGIN = 1, FUN = sd)
)

objects.store(conc.param, samps.j)
cat("Lists conc.params and samps.j stored.\n")

######################3

cat("Descriptive statistics:\n")

euRatio <- EvalOutput(euRatio)

ggplot(euAll@output, aes(x = euAllResult, colour = Fish))+geom_density()+
  facet_wrap(~ Compound) + scale_x_log10()

ggplot(euRatio@output, aes(x = euRatioResult, colour = Fish))+geom_density()+
  facet_wrap(~ Compound) 

# Predictions for all congeners of fish1 (Baltic herring)
scatterplotMatrix(t(samps.j$pred[1,,,1]))
# Means for all congeners of the generic fish
scatterplotMatrix(t(samps.j$mu1[,,1]))
# Prediction for all congeners of the generic fish
scatterplotMatrix(t(samps.j$pred1[,,1]))
# Predictions for all fish species for TCDD
scatterplotMatrix(t(samps.j$pred[,1,,1]))
# Predictions for pike for omegas and PeCDD
scatterplotMatrix(t(samps.j$Omega[6,2,,,1]))

ggplot(melt(exp(samps.j$pred[,,,1])), aes(x=value, colour=Compound))+geom_density()+
  facet_wrap( ~ Fish,scales = "free_y")+scale_x_log10()
ggplot(eu@output, aes(x = euResult, colour=Congener))+geom_density()+
  facet_wrap( ~ Fish, scales = "free_y")+scale_x_log10()
#stat_ellipse()

coda.j <- coda.samples(
  jags, 
  c('mu', 'pred', 'omega1', 'pred1'), 
  N
)

plot(coda.j)

#################

euAll <- euRaw[, c(1:4, 18, 19)]
colnames(euAll@output)[1:4] <- c("THLcode", "Matrix", "Compound", 
                                 "Fish")
if (exists("TEQgroups")) 
  levels(TEF$Group) <- TEQgroups$Out
temp <- oapply(euAll * TEF, cols = "Compound", FUN = "sum")
colnames(temp@output)[colnames(temp@output) == "Group"] <- "Compound"
eu2 <- combine(euAll, temp)
return(eu2)

Initiate concentration

Model run 19.5.2017 [15]

+ Show code - Hide code

# This is code Op_en3104/initiate on page [[EU-kalat]]

library(OpasnetUtils)

concentration <- Ovariable(
  "concentration",
  dependencies = data.frame(Name = "conc.param", Ident = "Op_en3104/bayes"),
  formula = function(...) {
    jsp <- lapply(
      1:length(conc.param$mu[,1]),
      FUN = function(x) {
        temp <- exp(mvrnorm(openv$N, conc.param$mu[x,], conc.param$Omega[x,,]))
        dimnames(temp) <- c(list(Iter = 1:openv$N), dimnames(conc.param$mu)[2])
        return(temp)
      }
    )
    names(jsp) <- dimnames(conc.param$mu)[[1]]
    jsp <- melt(jsp, value.name = "Result")
    colnames(jsp)[colnames(jsp)=="L1"] <- "Fish" # Convert automatic name to meaningful
    jsp <- Ovariable(output=jsp, marginal = colnames(jsp) != "Result")
    return(jsp)
  }
)

objects.store(concentration)
cat("Ovariable concentration stored.\n")

References

↑ A. Hallikainen, H. Kiviranta, P. Isosaari, T. Vartiainen, R. Parmanne, P.J. Vuorinen: Kotimaisen järvi- ja merikalan dioksiinien, furaanien, dioksiinien kaltaisten PCB-yhdisteiden ja polybromattujen difenyylieettereiden pitoisuudet. Elintarvikeviraston julkaisuja 1/2004. [1]
↑ E-R.Venäläinen, A. Hallikainen, R. Parmanne, P.J. Vuorinen: Kotimaisen järvi- ja merikalan raskasmetallipitoisuudet. Elintarvikeviraston julkaisuja 3/2004. [2]
↑ Anja Hallikainen, Riikka Airaksinen, Panu Rantakokko, Jani Koponen, Jaakko Mannio, Pekka J. Vuorinen, Timo Jääskeläinen, Hannu Kiviranta. Itämeren kalan ja muun kotimaisen kalan ympäristömyrkyt: PCDD/F-, PCB-, PBDE-, PFC- ja OT-yhdisteet. Eviran tutkimuksia 2/2011. ISSN 1797-2981 ISBN 978-952-225-083-4 [3]

[eukalatpop-1] A. Hallikainen, H. Kiviranta, P. Isosaari, T. Vartiainen, R. Parmanne, P.J. Vuorinen: Kotimaisen järvi- ja merikalan dioksiinien, furaanien, dioksiinien kaltaisten PCB-yhdisteiden ja polybromattujen difenyylieettereiden pitoisuudet. Elintarvikeviraston julkaisuja 1/2004. [1]

[eukalatmetallit-2] E-R.Venäläinen, A. Hallikainen, R. Parmanne, P.J. Vuorinen: Kotimaisen järvi- ja merikalan raskasmetallipitoisuudet. Elintarvikeviraston julkaisuja 3/2004. [2]

[eukalat2-3] Anja Hallikainen, Riikka Airaksinen, Panu Rantakokko, Jani Koponen, Jaakko Mannio, Pekka J. Vuorinen, Timo Jääskeläinen, Hannu Kiviranta. Itämeren kalan ja muun kotimaisen kalan ympäristömyrkyt: PCDD/F-, PCB-, PBDE-, PFC- ja OT-yhdisteet. Eviran tutkimuksia 2/2011. ISSN 1797-2981 ISBN 978-952-225-083-4 [3]

[1]

[2]

[3]

EU-kalat

Contents

Question

Answer

Rationale

Data

Calculations

Bayes model for dioxin concentrations

See also

References

Navigation menu