EU-kalat: Difference between revisions

From Opasnet
Jump to navigation Jump to search
(→‎Bayes model for dioxin concentrations: length and year in model)
Line 230: Line 230:
* 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. [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=k0n2CFnjdGBklm9E]
* 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. [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=k0n2CFnjdGBklm9E]
* Model run 22.3.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=2jX2XxWpiIEZPyzJ] Model does not mix well. Thinning gives little help?
* Model run 22.3.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=2jX2XxWpiIEZPyzJ] Model does not mix well. Thinning gives little help?
* Model run 25.3.2018 with conc.param as ovariable [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=DbcmZJmuZ0h0vaGx]


<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 367: Line 368:
     'lenp',# parameters for length
     'lenp',# parameters for length
     'timep', # parameter for Year
     'timep', # parameter for Year
     'pred' # predicted concentration for year 2009 and length 17 cm  
     'pred' # predicted concentration for year 2009 and length 17 cm
   ),  
   ),  
   #  thin=1000,
   #  thin=1000,
Line 376: Line 377:
dimnames(samps.j$lenp) <- list(Fish = fisl, 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$pred) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$timep) <- list(Param = "time", Iter = 1:N, Chain = 1:4)
dimnames(samps.j$timep) <- list(Dummy = "time", Iter = 1:N, Chain = 1:4)


##### conc.param contains expected values of the distribution parameters from the model
##### conc.param contains expected values of the distribution parameters from the model
conc.param <- list(
 
  Omega = apply(samps.j$Omega, MARGIN = 1:3, FUN = mean),
conc.param <- Ovariable(
  lenp.mean = apply(samps.j$lenp, MARGIN = 1, FUN = mean),
  "conc.param",
  lenp.sd = apply(samps.j$lenp, MARGIN = 1, FUN = sd),
  dependencies = data.frame(Name = "samps.j", Ident=NA),
  mu = apply(samps.j$mu, MARGIN = 1:2, FUN = mean),
  formula = function(...) {
  timep.mean = apply(samps.j$timep, MARGIN = 1, FUN = mean),
    conc.param <- list(
  timep.sd = apply(samps.j$timep, MARGIN = 1, FUN = sd)
      Omega = apply(samps.j$Omega, MARGIN = 1:3, FUN = mean),
      lenp = cbind(
        mean = apply(samps.j$lenp, MARGIN = 1, FUN = mean),
        sd = apply(samps.j$lenp, MARGIN = 1, FUN = sd)
      ),
      mu = apply(samps.j$mu, MARGIN = 1:2, FUN = mean),
      timep = cbind(
        mean = apply(samps.j$timep, MARGIN = 1, FUN = mean),
        sd = apply(samps.j$timep, MARGIN = 1, FUN = sd)
      )
    )
    names(dimnames(conc.param$lenp)) <- c("Fish","Metaparam")
    names(dimnames(conc.param$timep)) <- c("Dummy","Metaparam")
    conc.param <- melt(conc.param)
    colnames(conc.param)[colnames(conc.param)=="value"] <- "Result"
    colnames(conc.param)[colnames(conc.param)=="L1"] <- "Parameter"
    conc.param$Dummy <- NULL
    conc.param$Metaparam <- ifelse(is.na(conc.param$Metaparam), conc.param$Parameter, as.character(conc.param$Metaparam))
    return(Ovariable(output=conc.param, marginal=colnames(conc.param)!="Result"))
  }
)
)


if(FALSE){
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")
Line 406: Line 425:
ggplot(eu2@output, aes(x = euResult, colour=Compound))+stat_ecdf()+
ggplot(eu2@output, aes(x = euResult, colour=Compound))+stat_ecdf()+
   facet_wrap( ~ Fish)+scale_x_log10()
   facet_wrap( ~ Fish)+scale_x_log10()
}
fislen <- c(233, 170)
jsp <- lapply(1:length(conc.param$mu[, 1]), FUN = function(x) {
  temp <- exp(mvrnorm(
    openv$N,
    conc.param$mu[x, ]+conc.param$lenp.mean*(fislen[x]-170)+conc.param$timep.mean*(2009-2009),
    solve(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"
   
ggplot()+
  stat_ecdf(data=eu2@output, aes(x=euResult, colour="Data"))+
  stat_ecdf(data=melt(exp(samps.j$pred[,,,1])), aes(x=value, colour="Bayes estimate"))+
  stat_ecdf(data=jsp, aes(x=Result, colour="MC estimate"))+
  facet_grid(Compound ~ Fish)+scale_x_log10()


ggplot(eu2@output,  
ggplot(eu2@output,  
Line 436: Line 432:
scatterplotMatrix(t(exp(samps.j$pred[1,,,1])), main = "Predictions for all compounds for Baltic herring")
scatterplotMatrix(t(exp(samps.j$pred[1,,,1])), main = "Predictions for all compounds for Baltic herring")
scatterplotMatrix(t(exp(samps.j$pred[,1,,1])), main = "Predictions for all fish species for PCDDF")
scatterplotMatrix(t(exp(samps.j$pred[,1,,1])), main = "Predictions for all fish species for PCDDF")
scatterplotMatrix(t(samps.j$Omega[,1,1,,1]))
#scatterplotMatrix(t(cbind(samps.j$Omega[1,1,1,,1],samps.j$mu[1,1,,1])))


plot(coda.samples(jags, 'Omega', N))
plot(coda.samples(jags, 'Omega', N))
Line 459: Line 457:
conc_pcddf <- Ovariable(
conc_pcddf <- Ovariable(
   "conc_pcddf",
   "conc_pcddf",
   dependencies = data.frame(Name = "conc.param", Ident = "Op_en3104/bayes"),
   dependencies = data.frame(
   formula = function(...) {
    Name=c("conc.param","length","time"),
     require(MASS)
    Ident=c("Op_en3104/bayes",NA,NA)
    require(reshape2)
  ),
    jsp <- lapply(1:length(conc.param$mu[, 1]), FUN = function(x) {
   formula=function(...) {
      temp <- exp(mvrnorm(
     tmp1 <- length + time + conc.param + Ovariable(data=data.frame(Result="0-1")) # Ensures Iter
        openv$N,
    tmp2 <- unique(tmp1@output[setdiff(
        conc.param$mu[x, ],
      colnames(tmp1@output)[tmp1@marginal],
        solve(conc.param$Omega[ , ])
      c("Compound","Compound2","Metaparam","Parameter")
      ))
    )])
      dimnames(temp) <- c(list(Iter = 1:openv$N), dimnames(conc.param$mu)[2])
    tmp2$Row <- 1:nrow(tmp2)
      return(temp)
    tmp3 <- merge(tmp2,tmp1@output)
     })
    out <- data.frame()
     names(jsp) <- dimnames(conc.param$mu)[[1]]
     con <- levels(tmp1$Compound)
    jsp <- melt(jsp, value.name = "Result")
     for(i in 1:nrow(tmp2)) {
    colnames(jsp)[colnames(jsp) == "L1"] <- "Fish"
      tmp <- tmp3[tmp3$Row == i , ]
      mu <- tmp$conc.paramResult[tmp$Parameter=="mu"][match(con,tmp$Compound[tmp$Parameter=="mu"])] + # baseline
        rnorm(1,
          tmp$conc.paramResult[tmp$Parameter=="lenp" & tmp$Metaparam=="mean"][1],
          tmp$conc.paramResult[tmp$Parameter=="lenp" & tmp$Metaparam=="sd"][1]
        )* (tmp$lengthResult[1]-170) + # length
        rnorm(1,
          tmp$conc.paramResult[tmp$Parameter=="timep" & tmp$Metaparam=="mean"][1],
          tmp$conc.paramResult[tmp$Parameter=="timep" & tmp$Metaparam=="sd"][1]
        )* (tmp$timeResult[1]-2009) # time
      
      
    jsp <- Ovariable(
      Omega <- solve(tapply( # Is it sure that PCDDF and PCB are not mixed to wrong order?
      output = jsp,
        tmp$conc.paramResult[tmp$Parameter=="Omega"],
      marginal = colnames(jsp) != "Result"
        tmp[tmp$Parameter=="Omega", c("Compound","Compound2")],
    )
        sum # Equal to identity because only 1 row per cell.
      )) # Precision matrix
      
      
    d <- oapply(jsp, cols = "Compound", FUN = sum)
      rnd <- exp(mvrnorm(1, mu, Omega))
     d$Compound <- "TEQ"
      out <- rbind(out, merge(tmp2[i,], data.frame(Compound=names(rnd),Result=rnd)))
   
    }
     return(combine(jsp,d))
     out$Row <- NULL
  }
     return(out)
    }
)
)



Revision as of 20:19, 25 March 2018


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

Dioxin concentrations in Baltic herring.

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

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

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. Thinning gives little help?
  • Model run 25.3.2018 with conc.param as ovariable [26]

+ Show code

Initiate conc_pcddf

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

+ Show code

⇤--#: . These codes should be coherent with POPs in Baltic herring. --Jouni (talk) 12:14, 7 June 2017 (UTC) (type: truth; paradigms: science: attack)

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]