Difference between revisions of "EU-kalat"

Latest revision as of 07:26, 9 June 2017

This page is a study. The page identifier is Op_en3104
Moderator:Arja (see all)
Give your opinion to the peer rating of the content of this page. {{ #opasnet_rater: }}
Upload data Show results

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]
Model run 23.5.2017 with adjusted ovariables euRaw, eu, euRatio [7]

+ 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"   

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

    eu$Compound <- factor( # Compound levels are ordered based on the data table on [[TEF]]
      eu$Compound, 
      levels = unique(c(levels(TEF$Compound), levels(eu$Compound)))
    )
    eu$Compound <- eu$Compound[,drop=TRUE]
    
    return(eu)
  }
)

euRatio <- Ovariable(
  "euRatio",
  dependencies = data.frame(Name=c("eu")),
  formula = function(...) {
    euRatio <- eu[
      eu$Compound == "2378TCDD" & eu$Matrix == "Muscle" & result(eu) != 0 , ] # Zeros cannot be used in ratio estimates
    euRatio$Compound <- NULL
    euRatio <- log10(eu / euRatio)@output
    euRatio <- euRatio[!euRatio$Compound %in% c("2378TCDD", "2378-TCDD", "TCDD") , ]
    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.

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

#> unique(eu$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 PCDDFs 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   

indices <- list(
  Compound.TEQ2 = c("PCDDF", "PCB"),
  Compound.PCDDF14 = as.character(unique(euRaw@data$POP)[c(1:12, 14, 15)]), # 7 OCDD should be removed
  Fish.Fish14 = 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" 

objects.store(euRaw, eu, euRatio, indices)
cat("Ovariables euRaw, eu, and euRatio and list indices stored.\n")

Bayes model for dioxin concentrations

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

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

Model run 22.5.2017 with TEQdx and TEQpcb as the only Compounds [15]
Model run 23.5.2017 debugged [16] [17] [18]
Model run 24.5.2017 TEQdx, TECpcb -> PCDDF, PCB [19]

+ 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]] eu, euRatio, indices

eu2 <- eu <- EvalOutput(eu)

conl <- indices$Compound.TEQ2
fisl <- indices$Fish.Fish14
C <- length(conl)
Fi <- length(fisl)
N <- 1000
conl
fisl

replaces <- list(
  c("Chlorinated dibenzo-p-dioxins", "PCDDF"),
  c("Chlorinated dibenzofurans", "PCDDF"),
  c("Mono-ortho-substituted PCBs", "PCB"),
  c("Non-ortho-substituted PCBs", "PCB")
)

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

eu2 <- unkeep(eu2, prevresults = TRUE, sources = TRUE)
eu2 <- oapply(eu2, cols = "TEFversion", FUN = "sum") # This goes wrong if > 1 TEFversion

# Hierarchical Bayes model.

# PCDD/F concentrations in fish.
# It uses the TEQ sum of PCDD/F (PCDDF) as the total concentration
# of dioxin and PCB respectively for PCB in fish.
# PCDDF depends on age of fish, fish species and catchment area,
# but we only have species now so other variables are omitted.
# cong depends on fish species.

eu3 <- eu2[eu2$Compound %in% conl & eu2$Fish %in% fisl & eu2$Matrix == "Muscle" , ]
eu3 <- reshape(
  eu3@output, 
  v.names = "euResult", 
  idvar = "THLcode", 
  timevar = "Compound", 
  drop = c("Matrix"), 
  direction = "wide"
)

oprint(head(eu3))

#> colnames(eu3)
#[1] "THLcode"         "Fish"            "euResult.PCDDF"  "euResult.PCB"

# Find the level of quantification for dinterval function
LOQ <- unlist(lapply(eu3[3:ncol(eu3)], FUN = function(x) min(x[x!=0])))
names(LOQ) <- conl
cong <- data.matrix(eu3[3:ncol(eu3)])
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(eu3),
    C = C,
    Fi = Fi,
    cong = log(cong),
    fis = match(eu3$Fish, fisl),
    Omega0 = diag(C)/100000
  ),
  n.chains = 4,
  n.adapt = 100
)

update(jags, 100)

samps.j <- jags.samples(
  jags, 
  c(
    'mu', # mean by fish and compound
    'Omega', # precision matrix by fish and compound
    'pred', # predicted concentration by fish and compound
    #    'mu1', # mean prior for mu by compound
    'Omega1', # precision matrix by compound
    #    'tau1', # precision for prior of all mu 
    'pred1' # predicted concentration by compound
  ), 
  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$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")

# Leave only the main fish species and congeners and remove others
conl <- indices$Compound.PCDDF14
eu <- eu[eu$Compound %in% conl & eu$Fish %in% fisl , ] 

oprint(summary(
  eu,
  marginals = c("Fish", "Compound"), # Matrix is always 'Muscle'
  function_names = c("mean", "sd")
))

euRatio <- EvalOutput(euRatio)

oprint(summary(
  euRatio,
  marginals = c("Fish", "Compound"), # Matrix is always 'Muscle'
  function_names = c("mean", "sd")
))

ggplot(eu@output, aes(x = euResult, colour=Compound))+geom_density()+
  facet_wrap( ~ Fish, scales = "free_y")+scale_x_log10()
#stat_ellipse()

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

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

ggplot(melt(exp(samps.j$pred1[,,1])), aes(x=value, colour=Compound))+geom_density()+
  scale_x_log10()


scatterplotMatrix(t(samps.j$pred[1,,,1]), main = "Predictions for all compounds for Baltic herring")
## scatterplotMatrix(t(samps.j$mu1[,,1]), main = "Means for all compounds of the generic fish")
scatterplotMatrix(t(samps.j$pred1[,,1]), main = "Prediction for all compounds of the generic fish")
scatterplotMatrix(t(samps.j$pred[,1,,1]), main = "Predictions for all fish species for PCDDF")
scatterplotMatrix(t(samps.j$Omega[6,2,,,1]), main = "Predictions of Omega for pike and PCB")

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

plot(coda.j)

Initiate conc_pcddf

Model run 19.5.2017 [20]
Model run 23.5.2017 with bugs fixed [21]

+ Show code - Hide code

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

library(OpasnetUtils)

conc_pcddf <- Ovariable(
  "conc_pcddf",
  dependencies = data.frame(Name = "conc.param", Ident = "Op_en3104/bayes"),
  formula = function(...) {
    require(MASS)
    require(reshape2)
    jsp <- lapply(1:length(conc.param$mu[, 1]), FUN = function(x) {
      temp <- exp(mvrnorm(
        openv$N, 
        conc.param$mu[x, ],
        solve(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"
    jsp <- Ovariable(
      output = jsp,
      marginal = colnames(jsp) != "Result"
    )
    return(jsp)
  }
)

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

⇤# : These codes should be coherent with POPs in Baltic herring. --Jouni (talk) 12:14, 7 June 2017 (UTC)

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]

@@ Line 368: / Line 368: @@
 </rcode>
-'''Initiate concentration
+'''Initiate conc_pcddf
 * Model run 19.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=ystfGN6yfNwWNfnq]
 * Model run 23.5.2017 with bugs fixed [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=8iYF4GXFO9bUnld4]
-<rcode name="initiate" label="Initiate concentration (for developers only)">
+<rcode name="initiate" label="Initiate conc_pcddf (for developers only)">
 # This is code Op_en3104/initiate on page [[EU-kalat]]
 library(OpasnetUtils)
-concentration <- Ovariable(
+conc_pcddf <- Ovariable(
-   "concentration",
+   "conc_pcddf",
    dependencies = data.frame(Name = "conc.param", Ident = "Op_en3104/bayes"),
    formula = function(...) {
@@ Line 404: / Line 404: @@
 )
-objects.store(concentration)
+objects.store(conc_pcddf)
-cat("Ovariable concentration stored.\n")
+cat("Ovariable conc_pcddf stored.\n")
 </rcode>
+{{attack|# |These codes should be coherent with [[POPs in Baltic herring]].|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 12:14, 7 June 2017 (UTC)}}
 ==See also==

Difference between revisions of "EU-kalat"

Latest revision as of 07:26, 9 June 2017

Contents

Question

Answer

Rationale

Data

Calculations

Bayes model for dioxin concentrations

See also

References

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Page Tools

Tools

In other websites