EU-kalat: Difference between revisions

From Opasnet
Jump to navigation Jump to search
(→‎Calculations: code updated but not debugged)
Line 61: Line 61:
** Objects used in [[Benefit-risk assessment of Baltic herring and salmon intake]]
** Objects used in [[Benefit-risk assessment of Baltic herring and salmon intake]]
* Model run 25.1.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=wzisMQHAqcF30zcl]
* Model run 25.1.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=wzisMQHAqcF30zcl]
* Model run 22.5.2017 with new ovariables euRaw, euAll, euMain, and euRatio [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=7uTqQeaekwRFwA2J]


<rcode name="preprocess" label="Preprocess" graphics=1>
<rcode name="preprocess" label="Preprocess (for developers only)">
# This is code Op_en3104/preprocess on page [[EU-kalat]]
# This is code Op_en3104/preprocess on page [[EU-kalat]]
library(rjags)
library(OpasnetUtils)
library(OpasnetUtils)
library(ggplot2)
library(ggplot2)
Line 74: Line 74:
#[1] "Chlorinated dibenzo-p-dioxins" "Chlorinated dibenzofurans"     
#[1] "Chlorinated dibenzo-p-dioxins" "Chlorinated dibenzofurans"     
#[3] "Mono-ortho–substituted PCBs"  "Non-ortho–substituted PCBs"   
#[3] "Mono-ortho–substituted PCBs"  "Non-ortho–substituted PCBs"   
TEQgroups <- list(
  In = c(
    "Chlorinated dibenzo-p-dioxins",
    "Chlorinated dibenzofurans",
    "Mono-ortho–substituted PCBs",
    "Non-ortho–substituted PCBs"
  ),
  Out = c("TEQdx", "TEQdx", "TEQpcb", "TEQpcb")
)


euAll <- Ovariable(
euAll <- Ovariable(
Line 94: Line 84:
     euAll <- euRaw[,c(1:4, 18, 19)] # THL code, Matrix, Congener, Fish species
     euAll <- euRaw[,c(1:4, 18, 19)] # THL code, Matrix, Congener, Fish species
     colnames(euAll@output)[1:4] <- c("THLcode", "Matrix", "Compound", "Fish")
     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")
     temp <- oapply(euAll * TEF, cols = "Compound", FUN = "sum")
Line 104: Line 93:
)
)


euAll <- EvalOutput(euAll)
# Leave only the main fish species and congeners and remove others


################## Data for the main congeners and species only
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)
  }
)


#> unique(euAll$Congener)
euRatio <- Ovariable(
#[1] 2378TCDD    12378PeCDD   123478HCDD  123678HCDD  123789HCDD  1234678HpCDD
   "euRatio",
#[7] OCDD        2378TCDF    12378PeCDF   23478PeCDF  123478HCDF  123678HCDF 
   dependencies = data.frame(Name="euMain"),
#[13] 123789HCDF  234678HCDF  1234678HpCDF 1234789HpCDF OCDF        CoPCB77 ...   
  formula = function(...) {
 
     euRatio <- euMain[
# Remove the four with too little data (>70% BDL) and all non-PCDDF
      euMain$Compound == "2378TCDD" & euMain$Matrix == "Muscle" & result(euMain) != 0 , ] # Zeros cannot be used in ratio estimates
# aggregate(eu@data$euResult, by = eu@data["POP"], FUN = function(x) mean(x == 0))
    euRatio$Compound <- NULL
conl <- unique(euAll@output$Congener)[c(1:12, 14, 15)] # 7 OCDD should be removed
    euRatio <- log10(euMain / euRatio)@output
euMain <- euAll[euAll$Congener %in% conl , ]
    euRatio <- euRatio[euRatio$Compound %in% setdiff(conl, "2378TCDD") , ]
 
    return(euRatio)
#[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(euAll@output$Fish)[c(1:4, 6:14, 17)]
euMain <- euMain[euMain$Fish %in% fisl , ] # Remove four with too little data
 
euRatio <- euMain[
  euMain$Congener == "2378TCDD" & euMain$Matrix == "Muscle" & result(euMain) != 0 , ] # Zeros cannot be used in ratio estimates
euRatio$Congener <- NULL
euRatio <- log10(euMain / euRatio)@output
euRatio <- euRatio[euRatio$Congener %in% setdiff(conl, "2378TCDD") , ]


# Analysis: a few rows disappear here, as shown by numbers per fish species. Why?
# Analysis: a few rows disappear here, as shown by numbers per fish species. Why?
Line 143: Line 130:
# Conclusion: this is ok. Total 2292 rows.
# Conclusion: this is ok. Total 2292 rows.


ggplot(euAll@output, aes(x = euAllResult, colour = Fish))+geom_density()+
objects.store(euRaw, euAll, euMain, euRatio)
  facet_wrap(~ Compound) + scale_x_log10()
cat("Ovariables euRaw, euAll, euMain, and euRatio stored.\n")
 
ggplot(euRatio, aes(x = Result, colour = Fish))+geom_density()+
  facet_wrap(~ Compound)
 
objects.store(euAll, euMain, euRatio)
cat("Ovariables euAll and euMain, and data.frame eutRatio stored.\n")
</rcode>
</rcode>


Line 163: Line 144:
* Model run 19.5.2017 without ovariable concentration [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=bXhdwkBaQQi1LTcu] {{attack|# |The model does not mix well, so the results should not be used for final results.|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 19:37, 19 May 2017 (UTC)}}
* Model run 19.5.2017 without ovariable concentration [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=bXhdwkBaQQi1LTcu] {{attack|# |The model does not mix well, so the results should not be used for final results.|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 19:37, 19 May 2017 (UTC)}}
::{{comment|# |Maybe we should just estimate TEQs until the problem is fixed.|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 19:37, 19 May 2017 (UTC)}}
::{{comment|# |Maybe we should just estimate TEQs until the problem is fixed.|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 19:37, 19 May 2017 (UTC)}}
* Model run 22.5.2017 with TEQdx and TEQpcb as the only Compounds [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=2vTgALXXTzLgd4l1]


<rcode name="bayes" label="Sample Bayes model (for developers only)" graphics=1>
<rcode name="bayes" label="Sample Bayes model (for developers only)" graphics=1>
Line 176: Line 158:
objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]
objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]


# Hierarchical Bayes model.
################## 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)])


# PCDD/F concentrations in fish.
# conl
# It uses the sum of PCDD/F (Pcdsum) as the total concentration of dioxin in fish.
#[1] "2378TCDD"    "12378PeCDD"  "123478HCDD"  "123678HCDD"  "123789HCDD" 
# Cong_j is the fraction of a congener from pcdsum.
#[6] "1234678HpCDD" "OCDD"        "2378TCDF"    "12378PeCDF"  "23478PeCDF" 
# pcdsum is log-normally distributed. cong_j follows Dirichlet distribution.
#[11] "123478HCDF"  "123678HCDF"  "234678HCDF"  "1234678HpCDF"
# pcdsum depends on age of fish, fish species and catchment area, but we only have species now so other variables are omitted.
#> fisl
# cong_j depends on fish species.
#[1] "Baltic herring" "Sprat"          "Salmon"        "Sea trout"   
#[5] "Roach"          "Perch"          "Pike"          "Pike-perch"   
#[9] "Burbot"        "Whitefish"      "Flounder"      "Bream"       
#[13] "River lamprey"  "Rainbow trout"


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)
C <- length(conl)
Fi <- length(fisl)
Fi <- length(fisl)
Line 193: Line 196:
conl
conl
fisl
fisl
fishsamples <- reshape(
 
  eu@output,
euAll <- EvalOutput(euAll)
   v.names = "euResult",  
 
   idvar = "THLcode",  
replaces <- list(
  timevar = "Congener",  
   c("Chlorinated dibenzo-p-dioxins", "TEQdx"),
   drop = c("Matrix", "euSource"),  
   c("Chlorinated dibenzofurans", "TEQdx"),
   direction = "wide"
   c("Mono-ortho–substituted PCBs", "TEQpcb"),
   c("Non-ortho–substituted PCBs", "TEQpcb")
)
)


# Find the level of quantification for dinterval function
for(i in 1:length(replaces)) {
LOQ <- unlist(lapply(fishsamples[3:ncol(fishsamples)], FUN = function(x) min(x[x!=0])))
  levels(euAll$Compound)[replaces[[i]][1]] <- replaces[[i]][2]
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]))
euRatio <- EvalOutput(euRatio)
 
# Hierarchical Bayes model.


mod <- textConnection("
# PCDD/F concentrations in fish.
  model{
# It uses the sum of PCDD/F (Pcdsum) as the total concentration of dioxin in fish.
    for(i in 1:S) { # s = fish sample
# Cong_j is the fraction of a congener from pcdsum.
      #       below.LOQ[i,j] ~ dinterval(-cong[i,j], -LOQ[j])
# pcdsum is log-normally distributed. cong_j follows Dirichlet distribution.
      cong[i,1:C] ~ dmnorm(mu[fis[i],], Omega[fis[i],,])
# 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.
    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(
if(FALSE) { # Multicongener model NOT updated to reflect new names: eu -> euMain
  mod,
  conl <- as.character(unique(eu@output$Congener))
  data = list(
  fisl <- sort(as.character(unique(eu@output$Fish)))
    S = nrow(fishsamples),
  C <- length(conl)
    C = C,
  Fi <- length(fisl)
    Fi = Fi,
  N <- 1000
    cong = log(cong),
  conl
    #    LOQ = LOQ,
  fisl
    #    below.LOQ = is.na(cong)*1,
  fishsamples <- reshape(
    fis = match(fishsamples$Fish, fisl),
    eu@output,
    Omega0 = diag(C)/100000
    v.names = "euResult",
  ),
    idvar = "THLcode",
  n.chains = 4,
    timevar = "Congener",
  n.adapt = 100
    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)
} # if(FALSE)


conl <- c("TEQdx", "TEQpcb")
fishsamples <- unkeep(euAll, prevresults = TRUE, sources = TRUE)@output
fisl <- sort(as.character(unique(euMain$Fish))) # Take only the main fish species
fishsamples <- fishsamples[fishsamples$Compound %in% conl & fishsamples$Matrix == "Muscle" , ]
C <- length(conl)
Fi <- length(fisl)
N <- 1000
conl
fisl
fishsamples <- euAll@output[euAll$Compound %in% conl , ]
fishsamples <- reshape(
fishsamples <- reshape(
   fishsamples,  
   fishsamples,  
Line 260: Line 290:
   idvar = "THLcode",  
   idvar = "THLcode",  
   timevar = "Compound",  
   timevar = "Compound",  
   drop = c("Matrix", "euRawSource"),  
   drop = c("Matrix"),  
   direction = "wide"
   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
# Find the level of quantification for dinterval function
Line 271: Line 308:


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


jags <- jags.model(
jags <- jags.model(
Line 338: Line 375:
objects.store(conc.param, samps.j)
objects.store(conc.param, samps.j)
cat("Lists conc.params and samps.j stored.\n")
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)
# Predictions for all congeners of fish1 (Baltic herring)
Line 363: Line 412:


plot(coda.j)
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)
</rcode>
</rcode>


Revision as of 21:36, 22 May 2017


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

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

+ Show code

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

Initiate concentration

  • Model run 19.5.2017 [15]

+ Show code

See also

References

  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]
  2. ↑ E-R.Venäläinen, A. Hallikainen, R. Parmanne, P.J. Vuorinen: Kotimaisen järvi- ja merikalan raskasmetallipitoisuudet. Elintarvikeviraston julkaisuja 3/2004. [2]
  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]
Retrieved from "https://en.opasnet.org/index.php?title=EU-kalat&oldid=40760"