Difference between revisions of "EU-kalat"

From Testiwiki
Jump to: navigation, search
(→‎Calculations)
(→‎Bayes model for dioxin concentrations)
 
(18 intermediate revisions by the same user not shown)
Line 58: Line 58:
 
=== Calculations ===
 
=== Calculations ===
  
−
Model run 25.1.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=wzisMQHAqcF30zcl]
+
* Preprocess model 22.2.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=p53h4gxtFfDMOtHz]
 +
** 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 22.5.2017 with new ovariables euRaw, euAll, euMain, and euRatio [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=7uTqQeaekwRFwA2J]
 +
* Model run 23.5.2017 with adjusted ovariables euRaw, eu, euRatio [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=qkseWM9rmRysGwKM]
  
−
<rcode graphics=1>
+
<rcode name="preprocess" label="Preprocess (for developers only)">
−
# This is code Op_en3104/ on page [[EU-kalat]]
+
# This is code Op_en3104/preprocess on page [[EU-kalat]]
−
library(rjags)
 
 
library(OpasnetUtils)
 
library(OpasnetUtils)
 
library(ggplot2)
 
library(ggplot2)
 
library(reshape2)
 
library(reshape2)
  
−
eu <- Ovariable("eu", ddata = "Op_en3104", subset = "POPs")
+
euRaw <- Ovariable("euRaw", ddata = "Op_en3104", subset = "POPs")
−
eu <- EvalOutput(eu)
 
−
eu@output <- eu@output[c(1:4, 18, 19)] # THL code, Matrix, Congener, Fish species
 
−
eu@marginal <- eu@marginal[c(1:4, 18, 19)]
 
−
eu@output <- eu@output[result(eu) != 0 , ] # Zeros cannot be used in ratio estimates
 
−
colnames(eu@output)[colnames(eu@output) == "POP"] <- "Congener"
 
−
colnames(eu@output)[colnames(eu@output) == "Fish_species"] <- "Fish"
 
  
−
#> unique(eu@output$Congener)
+
# levels(TEF$Group)
−
#[1] 2378TCDD    12378PeCDD  123478HCDD  123678HCDD  123789HCDD  1234678HpCDD
+
#[1] "Chlorinated dibenzo-p-dioxins" "Chlorinated dibenzofurans"   
−
#[7] OCDD        2378TCDF    12378PeCDF  23478PeCDF  123478HCDF  123678HCDF 
+
#[3] "Mono-ortho–substituted PCBs"   "Non-ortho–substituted PCBs"    
−
#[13] 123789HCDF   234678HCDF   1234678HpCDF 1234789HpCDF OCDF        CoPCB77 ...   
 
  
−
# Remove the four with too little data (>70% BDL) and all non-PCDDF
+
eu <- Ovariable(
−
# aggregate(eu@data$euResult, by = eu@data["POP"], FUN = function(x) mean(x == 0))
+
  "eu",
−
conl <- unique(eu@output$Congener)[c(1:12, 14, 15)] # 7 OCDD should be removed
+
  dependencies = data.frame(
−
eu@output <- eu@output[eu@output$Congener %in% conl , ]
+
    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"
 +
    eu <- combine(eu, temp)
  
−
#[1] Baltic herring Sprat          Salmon        Sea trout      Vendace        
+
    eu$Compound <- factor( # Compound levels are ordered based on the data table on [[TEF]]
−
#[6] Roach          Perch          Pike          Pike-perch     Burbot       
+
      eu$Compound,
−
#[11] Whitefish      Flounder      Bream          River lamprey  Cod         
+
       levels = unique(c(levels(TEF$Compound), levels(eu$Compound)))
−
#[16] Trout          Rainbow trout  Arctic char    
+
    )
 +
    eu$Compound <- eu$Compound[,drop=TRUE]
 +
      
 +
    return(eu)
 +
   }
 +
)
  
−
fisl <- unique(eu@output$Fish)[c(1:4, 6:14, 17)]
+
euRatio <- Ovariable(
−
eu@output <- eu@output[eu@output$Fish %in% fisl , ] # Remove four with too little data
+
  "euRatio",
−
 
+
  dependencies = data.frame(Name=c("eu")),
−
eut <- eu
+
  formula = function(...) {
−
eut@output <- eut@output[
+
    euRatio <- eu[
−
  eut@output$Congener == "2378TCDD" & eut@output$Matrix == "Muscle" , ]
+
      eu$Compound == "2378TCDD" & eu$Matrix == "Muscle" & result(eu) != 0 , ] # Zeros cannot be used in ratio estimates
−
eut@output$Congener <- NULL
+
    euRatio$Compound <- NULL
−
eut@marginal[colnames(eut@output) == "Congener"] <- NULL
+
    euRatio <- log10(eu / euRatio)@output
−
eut <- log10(eu / eut)@output
+
    euRatio <- euRatio[!euRatio$Compound %in% c("2378TCDD", "2378-TCDD", "TCDD") , ]
−
conl <- conl[conl != "2378TCDD"]
+
    return(euRatio)
−
eut <- eut[eut$Congener %in% conl , ]
+
  }
 +
)
  
 
# 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 113: Line 124:
 
# Conclusion: this is ok. Total 2292 rows.
 
# Conclusion: this is ok. Total 2292 rows.
  
−
ggplot(eut, aes(x = Result, colour = Fish))+geom_density()+
+
################## Data for the main congeners and species only
−
   facet_wrap(~ Congener)  
+
 
 +
#> 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("PCDDF", "PCB"),
 +
  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")
 +
</rcode>
 +
 
 +
==== Bayes model for dioxin concentrations ====
 +
 
 +
* Model run 28.2.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=KrjCTGZmB8JkCH75]
 +
* Model run 28.2.2017 with corrected survey model [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=ZeO0SdlshPgOjqdL]
 +
* 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 23.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=8DnCPAKsMxGALkjs] produces list conc.param and ovariable concentration
 +
* Model run 24.4.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=ujtyawudKqJ7mmjn]
 +
* 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)}}
 +
* 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] [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=1P7ZPBbghEfisEcH] [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=BcZDhfjpv3fa4IRU]
 +
* Model run 24.5.2017 TEQdx, TECpcb -> PCDDF, PCB [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=kNNzEMTSD4N2f0Yy]
 +
 
 +
<rcode name="bayes" label="Sample Bayes model (for developers only)" graphics=1>
 +
# 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", "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 <- 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 (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(
 +
  eu3@output,
 +
  v.names = "euResult",
 +
  idvar = "THLcode",
 +
  timevar = "Compound",
 +
  drop = c("Matrix"),
 +
  direction = "wide"
 +
)
 +
 
 +
oprint(head(eu3))
 +
 
 +
#> colnames(eu3)
 +
#[1] "THLcode"        "Fish"            "euResult.PCDDF"  "euResult.PCB"
 +
 
 +
# 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])
 +
)
  
−
modf <- textConnection("
+
mod <- textConnection("
 
   model{
 
   model{
−
    for(i in 1:N) { # i = observation
+
  for(i in 1:S) { # s = fish sample
−
      eutt[i] ~ dnorm(mu[con[i],fis[i]], tau[con[i],fis[i]])
+
    #        below.LOQ[i,j] ~ dinterval(-cong[i,j], -LOQ[j])
−
    }
+
    cong[i,1:C] ~ dmnorm(mu[fis[i],], Omega[fis[i],,])
−
    for(j in 1:J) { # j = congener
+
  }
−
      for(k in 1:K) { # k = fish species
+
  for(i in 1:Fi) { # Fi = fish species
−
        mu[j,k] ~ dunif(-3,3)
+
    for(j in 1:C) {  
−
        tau[j,k] <- pow(sigma[j,k], -2)
+
      mu[i,j] ~ dnorm(mu1[j], tau1[j])
−
        sigma[j,k] ~ dunif(0, 10)
 
−
        pred[j,k] ~ dnorm(mu[j,k], tau[j,k]) # 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)
 +
    tau1[i] ~ dunif(0,10000)
 +
    pred1[i] ~ dnorm(mu1[i], tau1[i])
 +
  }
 +
  Omega1[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
 
   }
 
   }
 
")
 
")
  
−
# A binomial distribution is assumed for bins of answer choices.
 
 
jags <- jags.model(
 
jags <- jags.model(
−
   modf,
+
   mod,
 
   data = list(
 
   data = list(
−
     N = nrow(eut),
+
     S = nrow(eu3),
−
     J = length(conl),
+
     C = C,
−
     K = length(fisl),
+
     Fi = Fi,
−
     eutt = eut$Result,
+
     cong = log(cong),
−
     con = match(eut$Congener, conl),
+
     fis = match(eu3$Fish, fisl),
−
     fis = match(eut$Fish, fisl)
+
     Omega0 = diag(C)/100000
 
   ),
 
   ),
 
   n.chains = 4,
 
   n.chains = 4,
Line 147: Line 278:
 
)
 
)
  
−
update(jags, 1000)
+
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 , ]
  
−
samps <- jags.samples(jags, c('mu', 'pred'), 1000)
+
oprint(summary(
−
pl <- melt(array(samps$pred, dim = c(length(conl), length(fisl), 1000, 4),  
+
  eu,
−
   dimnames = list(
+
  marginals = c("Fish", "Compound"), # Matrix is always 'Muscle'
−
  Congener = conl,  
+
   function_names = c("mean", "sd")
−
  Fish = fisl,
 
−
  Iter = 1:1000,
 
−
  Seed = 1:4
 
−
  )
 
 
))
 
))
  
−
ggplot(pl, aes(x = value, colour = Fish))+geom_density(size = 1) +  
+
euRatio <- EvalOutput(euRatio)
−
   facet_wrap(~ Congener, scale = "free_y") + coord_cartesian(xlim = c(-1.5,1.5))
+
 
 +
oprint(summary(
 +
  euRatio,
 +
  marginals = c("Fish", "Compound"), # Matrix is always 'Muscle'
 +
  function_names = c("mean", "sd")
 +
))
 +
 
 +
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 = Compound))+geom_density()+
 +
  facet_wrap(~ Fish, 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()
  
−
ggplot(pl[pl$Fish %in% c("Baltic herring", "Rainbow trout", "Pike-perch", "Pike", "Salmon", "Sprat") , ], aes(x = value, colour = Fish))+geom_density(size = 1) +
 
−
  facet_wrap(~ Congener, scale = "free_y") + coord_cartesian(xlim = c(-1.5,1.5))
 
  
−
samps.coda <- coda.samples(jags, c('mu', 'pred'), 1000)
+
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 PCDDF")
 +
scatterplotMatrix(t(samps.j$Omega[6,2,,,1]), main = "Predictions of Omega for pike and PCB")
  
−
out <- data.frame(
+
coda.j <- coda.samples(
−
   Congener = conl,
+
   jags,  
−
   Fish = rep(fisl, each = length(conl)),
+
   c('mu', 'pred', 'Omega', 'pred1'),  
−
  Param = rep(c("mu", "pred"), each = length(conl)*length(fisl)),
+
   N
−
   Value = (summary(samps.coda)[[1]][,1])
 
 
)
 
)
  
−
oprint(out)
+
plot(coda.j)
 +
</rcode>
 +
 
 +
'''Initiate conc_pcddf
 +
 
 +
* Model run 19.5.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=ystfGN6yfNwWNfnq]
 +
* Model run 23.5.2017 with bugs fixed [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=8iYF4GXFO9bUnld4]
 +
 
 +
<rcode name="initiate" label="Initiate conc_pcddf (for developers only)">
 +
# This is code Op_en3104/initiate on page [[EU-kalat]]
  
−
ggplot(out[out$Param == "mu", ], aes(x = Congener, y = Value, colour = Fish, group = Fish)) +
+
library(OpasnetUtils)
−
  geom_line() + labs(title = "Congener amounts in fish compared with TCDD", y = "Log_10 ratio")  
 
  
−
if(FALSE){
+
conc_pcddf <- Ovariable(
−
  tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values")
+
  "conc_pcddf",
−
   tef <- EvalOutput(tef)
+
  dependencies = data.frame(Name = "conc.param", Ident = "Op_en3104/bayes"),
−
 
+
   formula = function(...) {
−
  #levels(eu$Congener) <- gsub("HCDD", "HxCDD", levels(eu$Congener))
+
    require(MASS)
−
  #levels(eu$Congener) <- gsub("HCDF", "HxCDF", levels(eu$Congener))
+
    require(reshape2)
−
  #levels(eu$Congener) <- gsub("CoPCB", "PCB", levels(eu$Congener))
+
    jsp <- lapply(1:length(conc.param$mu[, 1]), FUN = function(x) {
−
   euteq <- eu * tef
+
      temp <- exp(mvrnorm(
−
} # IF FALSE
+
        openv$N,
 +
        conc.param$mu[x, ],
 +
        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"
 +
    jsp <- Ovariable(
 +
      output = jsp,
 +
      marginal = colnames(jsp) != "Result"
 +
    )
 +
    return(jsp)
 +
   }
 +
)
 +
 
 +
objects.store(conc_pcddf)
 +
cat("Ovariable conc_pcddf stored.\n")
 
</rcode>
 
</rcode>
 +
 +
{{attack|# |These codes should be coherent with [[POPs in Baltic herring]].|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 12:14, 7 June 2017 (UTC)}}
  
 
==See also==
 
==See also==

Latest revision as of 07:26, 9 June 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 [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] [18]
  • Model run 24.5.2017 TEQdx, TECpcb -> PCDDF, PCB [19]

+ Show code

Initiate conc_pcddf

  • Model run 19.5.2017 [20]
  • Model run 23.5.2017 with bugs fixed [21]

+ Show code

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

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 "http://www.opasnet.org/testiwiki/index.php?title=EU-kalat&oldid=40926"