Difference between revisions of "EU-kalat"

Revision as of 13:15, 21 April 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]

+ Show code - Hide code

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

Bayes model for dioxin concentrations

Model run 28.2.2017 [6]
Model run 28.2.2017 with corrected survey model [7]
Model run 28.2.2017 with Mu estimates [8]
Model run 1.3.2017 [9]

+ Show code - Hide code

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

library(OpasnetUtils)
library(reshape2)
library(rjags)

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

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

conl <- as.character(unique(eu@output$Congener))
fisl <- sort(as.character(unique(eu@output$Fish)))
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[cong == 0] <- 0.001 # NA # Needed for dinterval

mod <- textConnection("
  model{
    for(i in 1:S) { # s = fish sample
#      for(j in 1:C) { # C = congener
        #        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:F) { # F = fish species
      mu[i,1:C] ~ dmnorm(mu0[1:C], S2[1:C,1:C])
      Omega[i,1:C,1:C] ~ dwish(S3[1:C,1:C],S)
      ans.pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction
    }
  }
")

C <- length(conl)
jags <- jags.model(
  mod,
  data = list(
    S = nrow(fishsamples),
    C = C,
    F = length(fisl),
    cong = log(cong),
    #    LOQ = LOQ,
    #    below.LOQ = is.na(cong)*1,
    fis = match(fishsamples$Fish, fisl),
    mu0 = rep(2,C),
    S2 = diag(C)/100000,
    S3 = diag(C)/10000
  ),
  n.chains = 4,
  n.adapt = 100
)
if(FALSE) {
mod <- textConnection("
  model{
                      for(i in 1:S) {
                      survs[i,1:4] ~ dmnorm(mus[], Omegas[,])
                      }
                      for(j in 1:H) {
                      survh[j,1:4] ~ dmnorm(muh[], Omegah[,])
                      }
                      mus[1:4] ~ dmnorm(mu0[1:4], S2[1:4,1:4])
                      Omegas[1:4,1:4] ~ dwish(S3[1:4,1:4],S)
                      anss.pred ~ dmnorm(mus[], Omegas[,])
                      muh[1:4] ~ dmnorm(mu0[1:4], S2[1:4,1:4])
                      Omegah[1:4,1:4] ~ dwish(S3[1:4,1:4],H)
                      ansh.pred ~ dmnorm(muh[], Omegah[,])
                      }
                      ")

jags <- jags.model(
  mod,
  data = list(
    survs = surv[surv$Eatsalm,c(8:11)],
    S = sum(surv$Eatsalm),
    survh = surv[surv$Eatherr,c(13:16)],
    H = sum(surv$Eatherr),
    mu0 = rep(2,4),
    S2 = diag(4)/100000,
    S3 = diag(4)/10000
  ),
  n.chains = 4,
  n.adapt = 100
)
}

update(jags, 100)

samps.j <- jags.samples(
  jags, 
  c('mu', 'Omega', 'ans.pred'), 
  1000
)
js <- array(
  c(
    samps.j$mu[,,,1],
    samps.j$Omega[,,1,,1],
    samps.j$Omega[,,2,,1],
    samps.j$Omega[,,3,,1],
    samps.j$Omega[,,4,,1],
    samps.j$Omega[,,5,,1],
    samps.j$Omega[,,6,,1],
    samps.j$Omega[,,7,,1],
    samps.j$Omega[,,8,,1],
    samps.j$Omega[,,9,,1],
    samps.j$Omega[,,10,,1],
    samps.j$Omega[,,11,,1],
    samps.j$Omega[,,12,,1],
    samps.j$Omega[,,13,,1],
    samps.j$Omega[,,14,,1],
    samps.j$ans.pred[,,,1]
  ),
  dim = c(14,14,1000,16),
  dimnames = list(
    Compound = 1:14,
    Fish = 1:14,
    Iter = 1:1000,
    Parameter = c("mu", paste("Omega", 1:14, sep=""), "ans.pred")
  )
)

# Mu for all congeners of fish1
scatterplotMatrix(t(js[,1,,1]))
# All parameters for TCDD of fish1
scatterplotMatrix(js[1,1,,])
# Mu for all fishes for TCDD
scatterplotMatrix(t(js[1,,,1]))

jsd <- melt(js)
ggplot(jsd, aes(x=value, colour=Question))+geom_density()+facet_grid(Parameter ~ Fish)
#ggplot(as.data.frame(js), aes(x = anss.pred, y = Sampled))+geom_point()+stat_ellipse()

coda.j <- coda.samples(
  jags, 
  c('mus', 'Omegas', 'anss.pred'), 
  1000
)

plot(coda.j)

######## fish.param contains expected values of the distribution parameters from the model

fish.param <- list(
  mu = apply(js[,,1,], MARGIN = c(1,3), FUN = mean),
  Omega = lapply(
    1:2,
    FUN = function(x) {
      solve(apply(js[,,2:5,], MARGIN = c(1,3,4), FUN = mean)[,,x])
    } # solve matrix: precision->covariace
  )
)

jsp <- Ovariable(
  "jsp",
  dependencies = data.frame(Name = "fish.param"),
  formula = function(...) {
    jsp <- lapply(
      1:2,
      FUN = function(x) {
        mvrnorm(openv$N, fish.param$mu[,x], fish.param$Omega[[x]])
      }
    )
    
    jsp <- rbind(
      cbind(
        Fish = "Salmon",
        Iter = 1:nrow(jsp[[1]]),
        as.data.frame(jsp[[1]])
      ),
      cbind(
        Fish = "Herring",
        Iter = 1:nrow(jsp[[2]]),
        as.data.frame(jsp[[2]])
      )
    )
    jsp <- melt(jsp, id.vars = c("Iter", "Fish"), variable.name = "Question", value.name = "Result")
    jsp <- Ovariable(output=jsp, marginal = colnames(jsp) %in% c("Iter", "Fish", "Question"))
    return(jsp)
  }
)

# Combine modelled survey answers with estimated amounts and frequencies by:
# Rounding the modelled result and merging that with value in ovariable assump

often <- Ovariable(
  "often",
  dependencies = data.frame(Name=c("jsp","assump")),
  formula = function(...) {
    out <- jsp[jsp$Question == "1" , !colnames(jsp@output) %in% c("Question")]
    out$Value <- round(result(out))
    out <- merge(
      assump@output[assump$Variable == "freq",],
      out@output
    )
    colnames(out)[colnames(out) == "assumpResult"] <- "Result"
    return(out[!colnames(out) %in% c("Value", "Variable")])
  }
)

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

update(jags, 100)
samps <- jags.samples(jags, c('mu', 'pcd.pred'), 1000)
#samps.coda <- coda.samples(jags, c('mu', 'pcd.pred', 'ans.pred'), 1000)
#objects.store(samps)

#library(plyr)
#temp <- adply(samps$mu, c(1,2,3,4))

pcd.pred <- array(
  samps$pcd.pred, 
  dim = c(length(fisl), length(conl), 1000, 4), 
  dimnames = list(
    Fish = fisl,
    Congener = conl, 
    Iter = 1:1000,
    Seed = c("S1","S2","S3","S4")
  )
)

mu.pred <- array(
  samps$mu, 
  dim = c(length(fisl), length(conl), 1000, 4), 
  dimnames = list(
    Fish = fisl,
    Congener = conl, 
    Iter = 1:1000,
    Seed = c("S1","S2","S3","S4")
  )
)

#objects.store(pcd.pred, mu.pred)
#cat("Arrays pcd.pred, mu.pred 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]

Difference between revisions of "EU-kalat"

Revision as of 13:15, 21 April 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