EU-kalat: Difference between revisions

Revision as of 05:11, 22 March 2018

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 this link or by running the codel 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()

temp2 <- temp[temp$Congener %in% levels(temp$Congener)[1:17],]

ggplot(temp2, aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") + facet_wrap(~ Catch_location) + 
  labs(
    title="Dioxin concentrations in Baltic herring",
    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

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]
Model run 11.10.2017: Small herring and Large herring added as new species [8]
Model rerun 15.11.2017 because the previous stored run was lost in update [9]
Model run 21.3.2018: Small and large herring replaced by actual fish length [10]

+ Show code - Hide code

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

# Divide herring into small and large (as new "species")
euRaw <- Ovariable("euRaw", ddata = "Op_en3104", subset = "POPs")
#euRaw <- EvalOutput(euRaw)

#temp <- euRaw[euRaw$Fish_species == "Baltic herring" , ]
#temp$Fish_species <- as.factor(
#  ifelse(
#    as.numeric(as.character(temp$Length_mean_mm)) >= 170,
#    "Large herring",
#    "Small herring"
#  )
#)
#euRaw <- combine(euRaw, temp, name = "euRaw")

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

indexguide <- Ovariable(data=data.frame( # This is needed to make these indices marginals in eu (which is in long format).
  Length = 1,
  Year = 1,
  Weight = 1,
  Large = 1,
  Result = 1
))

eu <- Ovariable(
  "eu",
  dependencies = data.frame(
    Name=c("euRaw", "TEF", "indexguide"),
    Ident=c(NA,"Op_en4017/initiate", NA)
  ),
  formula = function(...) {
    euRaw$Length<-as.numeric(as.character(euRaw$Length_mean_mm))
    euRaw$Year <- as.numeric(substr(euRaw$Catch_date, nchar(as.character(euRaw$Catch_date))-3,100))
    euRaw@marginal[colnames(euRaw@output) %in% c("Length","Year")] <- TRUE
#    euRaw$Weight<-as.numeric(as.character(euRaw$Weight_mean_g))
#    euRaw$Large <- euRaw$Length >= 170
    
    eu <- euRaw[,c(1:4, 10, 20, 21, 18)] # See below + Length, Year, Result
    colnames(eu@output)[1:5] <- c("THLcode", "Matrix", "Compound", "Fish", "N")
    
    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   Small herring   Large herring

indices <- list(
  Compound.TEQ2 = c("PCDDF", "PCB"),
  Compound.PCDDF14 = as.character(unique(euRaw$POP)[c(1:12, 14, 15)]), # 7 OCDD should be removed
  Fish.Fish16 = as.character(unique(euRaw$Fish_species)[c(1:4, 6:14, 17, 19, 20)])
)

#> indices
#$Compound.TEQ2
#[1] "PCDDF" "PCB"  
#
#$Compound.PCDDF14
# [1] "2378TCDD"     "12378PeCDD"   "123478HCDD"   "123678HCDD"   "123789HCDD"   "1234678HpCDD"
# [7] "OCDD"         "2378TCDF"     "12378PeCDF"   "23478PeCDF"   "123478HCDF"   "123678HCDF"  
#[13] "234678HCDF"   "1234678HpCDF"
#
#$Fish.Fish16
# [1] "Baltic herring" "Sprat"          "Salmon"         "Sea trout"      "Roach"         
# [6] "Perch"          "Pike"           "Pike-perch"     "Burbot"         "Whitefish"     
#[11] "Flounder"       "Bream"          "River lamprey"  "Rainbow trout"  "Small herring" 
#[16] "Large herring" 

#ggplot(euRaw@output, aes(x=Length, y=Weight))+geom_point()+
#  facet_wrap(~Fish_species, scales="free")

#oprint(summary(
#  euRaw[euRaw$Fish_species %in% c("Baltic herring","Salmon") & euRaw$POP == "2378TCDD" , ],
#  marginals=c("Catch_location","Large","Fish_species"),
#  fun = "length"
#))

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

Bayes model for dioxin concentrations

Model run 28.2.2017 [11]
Model run 28.2.2017 with corrected survey model [12]
Model run 28.2.2017 with Mu estimates [13]
Model run 1.3.2017 [14]
Model run 23.4.2017 [15] produces list conc.param and ovariable concentration
Model run 24.4.2017 [16]
Model run 19.5.2017 without ovariable concentration [17] ⇤--#: . 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 [18]
Model run 23.5.2017 debugged [19] [20] [21]
Model run 24.5.2017 TEQdx, TECpcb -> PCDDF, PCB [22]
Model run 11.10.2017 with small and large herring [23] (removed in update)
Model run 12.3.2018: bugs fixed with data used in Bayes. In addition, redundant fish species removed and Omega assumed to be the same for herring and salmon. [24]
Model run 22.3.2018 [25] Model does not mix well. Use thinning?

Warning messages: 1: In reshapeWide(data, idvar = idvar, timevar = timevar, varying = varying, : multiple rows match for Compound=PCDDF: first taken 2: In reshapeWide(data, idvar = idvar, timevar = timevar, varying = varying, : multiple rows match for Compound=PCB: first taken Error in jags.model(mod, data = list(S = nrow(eu3), C = C, Fi = Fi, cong = log(cong), : RUNTIME ERROR: Compilation error on line 5. Index out of range taking subset of cong

+ 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.Fish16
# fisl <- c("Baltic herring","Large herring","Salmon","Small herring")
fisl <- c("Baltic herring","Salmon")
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@marginal[colnames(eu2@output) %in% c("Length","Year")] <- TRUE # Indexguide should take care of this but it doesn't!
eu2 <- unkeep(eu2, prevresults = TRUE, sources = TRUE)
eu2 <- oapply(eu2, cols = "TEFversion", FUN = "sum") # Sums up dioxin+furan and non+monoortho. 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( # Gives warning that more than one value and first taken. Why?
  eu3@output, 
  v.names = "euResult", 
  idvar = c("THLcode", "Fish"),
  timevar = "Compound", 
  drop = c("Matrix"), 
  direction = "wide"
)

oprint(head(eu3))

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

cong <- data.matrix(eu3[6:ncol(eu3)])

if(FALSE){
  # Find the level of quantification for dinterval function
  LOQ <- unlist(lapply(eu3[6:ncol(eu3)], FUN = function(x) min(x[x!=0]))) 
  # With TEQ, there are no zeros. So this is useful only if there are congener-specific results.
  names(LOQ) <- conl


  cong <- sapply(
    1:length(LOQ),
    FUN = function(x) ifelse(cong[,x]==0, 0.5*LOQ[x], cong[,x])
  )
}
# This version of the model looks only at Baltic herring, Large herring, small herring and salmon.
# It assumes that all fish groups have the same Omega but mu varies.

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(muind[i,], Omega[,])
  muind[i,1:C] <- mu[fis[i],1:C] + lenp[fis[i]]*length[i] + timep*year[i]
  }
  
  # Priors for parameters
  timep ~ dnorm(-0.1, 1) # Time trend

  for(i in 1:Fi) { # Fi = fish species
  lenp[i] ~ dnorm(0,0.0001) # length parameter 
  pred[i,1:C] ~ dmnorm(mu[i,1:C]+lenp[i]*lenpred[i]+timep*timepred, Omega[,]) # Model prediction
  for(j in 1:C) { 
    mu[i,j] ~ dnorm(0, 0.0001) # mu1[j], tau1[j]) # Congener-specific mean for fishes
  }
  }
  Omega[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),
    length = eu3$Length,
    year = eu3$Year,
    fis = match(eu3$Fish, fisl),
    lenpred = c(170, 860),
    timepred = 2009,
    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 compound
    'lenp',# parameters for length
    'timep', # parameter for Year
    'pred' # predicted concentration for year 2009 and 17 cm herring, 80 cm salmon
  ), 
  N
)
dimnames(samps.j$mu) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$lenp) <- list(Fish = fisl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$pred) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$Omega) <- 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(
  Omega = apply(samps.j$Omega[,,,1], MARGIN = 1:2, FUN = mean),
  lenp.mean = apply(samps.j$lenp[,,1], MARGIN = 1, FUN = mean),
  lenp.sd = apply(samps.j$lenp[,,1], MARGIN = 1, FUN = sd),
  mu = apply(samps.j$mu[,,,1], MARGIN = 1:2, FUN = mean),
  timep.mean = apply(samps.j$lenp[,,1], MARGIN = 1, FUN = mean),
  timep.sd = apply(samps.j$lenp[,,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

#oprint(summary(
#  eu2[eu2$Compound %in% indices$Compound.PCDDF14 & eu$Fish %in% fisl , ],
#  marginals = c("Fish", "Compound"), # Matrix is always 'Muscle'
#  function_names = c("mean", "sd")
#))

ggplot(eu2@output[eu2$Fish %in% fisl & eu2$Compound %in% conl,], aes(x = euResult, colour=Compound))+stat_ecdf()+
  facet_wrap( ~ Fish)+scale_x_log10()

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

scatterplotMatrix(t(samps.j$pred[1,,,1]), main = "Predictions for all compounds for Baltic herring")
scatterplotMatrix(t(samps.j$pred[,1,,1]), main = "Predictions for all fish species for PCDDF")

coda.j <- coda.samples(
  jags, 
  c('mu', 'pred', 'timep','lenp'), 
  N
)

plot(coda.j)

Initiate conc_pcddf

Model run 19.5.2017 [26]
Model run 23.5.2017 with bugs fixed [27]
Model run 12.10.2017: TEQ calculation added [28]
Model rerun 15.11.2017 because the previous stored run was lost in update [29]
12.3.2018 adjusted to match the same Omega for all fish species [30]

+ 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[ , ])
      ))
      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"
    )
    
    d <- oapply(jsp, cols = "Compound", FUN = sum)
    d$Compound <- "TEQ"
    
    return(combine(jsp,d))
  }
)

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) (type: truth; paradigms: science: attack)

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: Difference between revisions

Revision as of 05:11, 22 March 2018

Contents

Question

Answer

Rationale

Data

Calculations

Preprocess

Bayes model for dioxin concentrations

Initiate conc_pcddf

See also

References

Navigation menu