Difference between revisions of "EU-kalat"

Revision as of 16:02, 23 May 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"
    eu2 <- combine(eu, temp)
    
    return(eu2)
  }
)

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("TEQdx", "TEQpcb"),
  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]

+ 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", "TEQdx"),
  c("Chlorinated dibenzofurans", "TEQdx"),
  c("Mono-ortho-substituted PCBs", "TEQpcb"),
  c("Non-ortho-substituted PCBs", "TEQpcb")
)

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 (TEQdx) as the total concentration
# of dioxin and TEQpcb respectively for PCB in fish.
# TEQdx 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.TEQdx"  "euResult.TEQpcb"

# 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 , ] 

euRatio <- EvalOutput(euRatio)

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 = Fish))+geom_density()+
  facet_wrap(~ Compound, 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 TEQdx")
scatterplotMatrix(t(samps.j$Omega[6,2,,,1]), main = "Predictions of Omega for pike and TEQpcb")

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

#plot(coda.j)

Initiate concentration

Model run 19.5.2017 [18]

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

Difference between revisions of "EU-kalat"

Revision as of 16:02, 23 May 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

@@ Line 164: / Line 164: @@
 ::{{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]
-* Model run 23.5.2017 debugged [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=rMSAZy6PSKzKhHwp]
+* Model run 23.5.2017 debugged [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=rMSAZy6PSKzKhHwp] [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=1P7ZPBbghEfisEcH]
 <rcode name="bayes" label="Sample Bayes model (for developers only)" graphics=1>
@@ Line 176: / Line 176: @@
 library(car) # scatterplotMatrix
-objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]
+objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]] eu, euRatio, indices
+eu2 <- eu <- EvalOutput(eu)
 conl <- indices$Compound.TEQ2
@@ Line 186: / Line 188: @@
 fisl
-eu2 <- EvalOutput(eu)
-eu <- eu2
 replaces <- list(
    c("Chlorinated dibenzo-p-dioxins", "TEQdx"),
@@ Line 196: / Line 196: @@
 for(i in 1:length(replaces)) {
-#  print(replaces[[i]][1])
+   levels(eu2$Compound)[levels(eu2$Compound)==replaces[[i]][1]] <- replaces[[i]][2]
-#  print(levels(eu$Compound)[levels(eu$Compound)==replaces[[i]][1]])
-   levels(eu$Compound)[levels(eu$Compound)==replaces[[i]][1]] <- replaces[[i]][2]
 }
-eu <- oapply(eu, cols = "TEFversion", FUN = "sum") # This goes wrong if > 1 TEFversion
+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 sum of PCDD/F (Pcdsum) as the total concentration of dioxin in fish.
+# It uses the TEQ sum of PCDD/F (TEQdx) as the total concentration
-# Cong_j is the fraction of a congener from pcdsum.
+# of dioxin and TEQpcb respectively for PCB in fish.
-# pcdsum is log-normally distributed. cong_j follows Dirichlet distribution.
+# TEQdx depends on age of fish, fish species and catchment area,
-# pcdsum depends on age of fish, fish species and catchment area, but we only have species now so other variables are omitted.
+# but we only have species now so other variables are omitted.
-# cong_j depends on fish species.
+# cong depends on fish species.
-if(FALSE) { # Multicongener model NOT updated to reflect new names: eu -> euMain
+eu3 <- eu2[eu2$Compound %in% conl & eu2$Fish %in% fisl & eu2$Matrix == "Muscle" , ]
-  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
-  eu3 <- 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(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),
-      #    LOQ = LOQ,
-      #    below.LOQ = is.na(cong)*1,
-      fis = match(eu3$Fish, fisl),
-      Omega0 = diag(C)/100000
-    ),
-    n.chains = 4,
-    n.adapt = 100
-  )
-} # if(FALSE)
-eu3 <- unkeep(eu, prevresults = TRUE, sources = TRUE)@output
-eu3 <- eu3[eu3$Compound %in% conl & eu$Fish %in% fisl & eu3$Matrix == "Muscle" , ]
 eu3 <- reshape(
-   eu3,
+   eu3@output,
    v.names = "euResult",
    idvar = "THLcode",
@@ Line 301: / Line 237: @@
 mod <- textConnection("
    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)
-      pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction
      }
+    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)
@@ Line 317: / Line 253: @@
      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)
    }
 ")
@@ Line 328: / Line 264: @@
      Fi = Fi,
      cong = log(cong),
-    #    LOQ = LOQ,
-    #    below.LOQ = is.na(cong)*1,
      fis = match(eu3$Fish, fisl),
      Omega0 = diag(C)/100000
@@ Line 341: / Line 275: @@
 samps.j <- jags.samples(
    jags,
-   c('mu', 'Omega', 'pred', 'mu1', 'Omega1', 'tau1', 'pred1'),
+   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$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$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)
@@ Line 359: / Line 300: @@
    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),
+#  mu1 = apply(samps.j$mu1[,,1], MARGIN = 1, FUN = mean),
-  tau1 = apply(samps.j$tau1[,,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)
@@ Line 372: / Line 314: @@
 cat("Descriptive statistics:\n")
-eu <- eu2
+# Leave only the main fish species and congeners and remove others
-euRatio <- EvalOutput(euRatio)
 conl <- indices$Compound.PCDDF14
-# Leave only the main fish species and congeners and remove others
 eu <- eu[eu$Compound %in% conl & eu$Fish %in% fisl , ]
-ggplot(eu@output, aes(x = euResult, colour = Fish))+geom_density()+
+euRatio <- EvalOutput(euRatio)
-  facet_wrap(~ Compound) + scale_x_log10()
-euRatio <- euRatio[euRatio$Compound %in% conl & euRatio$Fish %in% fisl , ]
+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 = Fish))+geom_density()+
    facet_wrap(~ Compound, scales = "free_y")
-# 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=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()+
-#stat_ellipse()
+   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 TEQdx")
+scatterplotMatrix(t(samps.j$Omega[6,2,,,1]), main = "Predictions of Omega for pike and TEQpcb")
 coda.j <- coda.samples(