Difference between revisions of "EU-kalat"

From Testiwiki
Jump to: navigation, search
(Bayes model for dioxin concentrations: first draft of new model)
(Bayes model for dioxin concentrations)
Line 129: Line 129:
 
* Model run 28.2.2017 with Mu estimates [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=TwY2bAIiWr037zqb]
 
* Model run 28.2.2017 with Mu estimates [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=TwY2bAIiWr037zqb]
 
* Model run 1.3.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=3Xu19vkWK1lyWVg3]
 
* Model run 1.3.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=3Xu19vkWK1lyWVg3]
 +
* Model run 23.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=yblH8XsmSUp0yRPL] produces list conc.param and ovariable concantration
  
<rcode name="bayes" label="Sample Bayes model (for developers only)">
+
<rcode name="bayes" label="Sample Bayes model (for developers only)" graphics=1>
 
# This is code Op_en3014/bayes on page [[EU-kalat]]
 
# This is code Op_en3014/bayes on page [[EU-kalat]]
  
Line 136: Line 137:
 
library(reshape2)
 
library(reshape2)
 
library(rjags)
 
library(rjags)
 +
library(MASS) # mvrnorm
 +
library(car) # scatterplotMatrix
  
 
objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]
 
objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]
Line 150: Line 153:
 
conl <- as.character(unique(eu@output$Congener))
 
conl <- as.character(unique(eu@output$Congener))
 
fisl <- sort(as.character(unique(eu@output$Fish)))
 
fisl <- sort(as.character(unique(eu@output$Fish)))
 +
C <- length(conl)
 +
Fi <- length(fisl)
 +
N <- 1000
 
conl
 
conl
 
fisl
 
fisl
Line 165: Line 171:
 
names(LOQ) <- conl
 
names(LOQ) <- conl
 
cong <- data.matrix(fishsamples[3:ncol(fishsamples)])
 
cong <- data.matrix(fishsamples[3:ncol(fishsamples)])
cong[cong == 0] <- 0.001 # NA # Needed for dinterval
+
cong <- sapply(1:length(LOQ), FUN = function(x) ifelse(cong[,x]==0, 0.5*LOQ[x], cong[,x]))
  
 
mod <- textConnection("
 
mod <- textConnection("
 
   model{
 
   model{
 
     for(i in 1:S) { # s = fish sample
 
     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])
        #        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:F) { # F = fish species
+
     for(i in 1:Fi) { # Fi = fish species
       mu[i,1:C] ~ dmnorm(mu0[1:C], S2[1:C,1:C])
+
       mu[i,1:C] ~ dmnorm(mu0[1:C], Omega0[1:C,1:C])
       Omega[i,1:C,1:C] ~ dwish(S3[1:C,1:C],S)
+
       Omega[i,1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
 
       ans.pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction
 
       ans.pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction
 
     }
 
     }
Line 183: Line 187:
 
")
 
")
  
C <- length(conl)
 
 
jags <- jags.model(
 
jags <- jags.model(
 
   mod,
 
   mod,
Line 189: Line 192:
 
     S = nrow(fishsamples),
 
     S = nrow(fishsamples),
 
     C = C,
 
     C = C,
     F = length(fisl),
+
     Fi = Fi,
 
     cong = log(cong),
 
     cong = log(cong),
 
     #    LOQ = LOQ,
 
     #    LOQ = LOQ,
 
     #    below.LOQ = is.na(cong)*1,
 
     #    below.LOQ = is.na(cong)*1,
 
     fis = match(fishsamples$Fish, fisl),
 
     fis = match(fishsamples$Fish, fisl),
     mu0 = rep(2,C),
+
     mu0 = rep(0,C),
     S2 = diag(C)/100000,
+
     Omega0 = 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.chains = 4,
 
   n.adapt = 100
 
   n.adapt = 100
 
)
 
)
}
 
  
 
update(jags, 100)
 
update(jags, 100)
Line 240: Line 209:
 
   jags,  
 
   jags,  
 
   c('mu', 'Omega', 'ans.pred'),  
 
   c('mu', 'Omega', 'ans.pred'),  
   1000
+
   N
)
 
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")
 
  )
 
 
)
 
)
 +
dimnames(samps.j$mu) <- list(Fish = fisl, 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$ans.pred) <- list(Fish = fisl, Compound=conl,Iter=1:N, Chain=1:4)
  
# Mu for all congeners of fish1
+
##### cong.param contains expected values of the distribution parameters from the model
scatterplotMatrix(t(js[,1,,1]))
+
cong.param <- list(
# All parameters for TCDD of fish1
+
   mu = apply(samps.j$mu[,,,1], MARGIN = c(1,2), FUN = mean),
scatterplotMatrix(js[1,1,,])
+
   Omega = apply(samps.j$Omega[,,,,1], MARGIN = c(1,2,3), FUN = mean)
# 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(
+
concentration <- Ovariable(
   "jsp",
+
   "concentration",
   dependencies = data.frame(Name = "fish.param"),
+
   dependencies = data.frame(Name = "cong.param"),
 
   formula = function(...) {
 
   formula = function(...) {
 
     jsp <- lapply(
 
     jsp <- lapply(
       1:2,
+
       1:length(cong.param$mu[,1]),
 
       FUN = function(x) {
 
       FUN = function(x) {
         mvrnorm(openv$N, fish.param$mu[,x], fish.param$Omega[[x]])
+
         temp <- exp(mvrnorm(openv$N, cong.param$mu[x,], cong.param$Omega[x,,]))
 +
        dimnames(temp) <- c(list(Iter = 1:openv$N), dimnames(cong.param$mu)[2])
 +
        return(temp)
 
       }
 
       }
 
     )
 
     )
      
+
     names(jsp) <- dimnames(cong.param$mu)[[1]]
    jsp <- rbind(
+
    jsp <- melt(jsp, value.name = "Result")
      cbind(
+
    colnames(jsp)[colnames(jsp)=="L1"] <- "Fish" # Convert automatic name to meaningful
        Fish = "Salmon",
+
     jsp <- Ovariable(output=jsp, marginal = colnames(jsp) != "Result")
        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)
 
     return(jsp)
 
   }
 
   }
 
)
 
)
  
# Combine modelled survey answers with estimated amounts and frequencies by:
+
#objects.store(concentration, conc.param)
# Rounding the modelled result and merging that with value in ovariable assump
+
#cat("Ovariable concentration and list conc.params stored.\n")
  
often <- Ovariable(
+
# Predictions for all congeners of fish1 (Baltic herring)
  "often",
+
scatterplotMatrix(t(samps.j$ans.pred[1,,,1]))
  dependencies = data.frame(Name=c("jsp","assump")),
+
# Predictions for all fish species for TCDD
  formula = function(...) {
+
scatterplotMatrix(t(samps.j$ans.pred[,1,,1]))
    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)
+
concentration <- EvalOutput(concentration)
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)
+
ggplot(melt(exp(samps.j$ans.pred[,,,1])), aes(x=value, colour=Compound))+geom_density()+
#temp <- adply(samps$mu, c(1,2,3,4))
+
  facet_wrap( ~ Fish,scales = "free_y")+scale_x_log10()
 +
ggplot(eu@output, aes(x = euResult, colour=Congener))+geom_density()+
 +
  facet_wrap( ~ Fish, scales = "free_y")+scale_x_log10()
 +
#stat_ellipse()
  
pcd.pred <- array(
+
#coda.j <- coda.samples(
  samps$pcd.pred,
+
#  jags,  
  dim = c(length(fisl), length(conl), 1000, 4),
+
c('mu', 'Omega', 'ans.pred'),  
  dimnames = list(
+
#  N
    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)
+
#plot(coda.j)
#cat("Arrays pcd.pred, mu.pred stored.\n")
 
 
</rcode>
 
</rcode>
  

Revision as of 13:12, 23 April 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 [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]
  • Model run 23.4.2017 [10] produces list conc.param and ovariable concantration

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