EU-kalat: Difference between revisions

From Opasnet
Jump to navigation Jump to search
(→‎Calculations: updated but contains bugs)
Line 77: Line 77:
* Model run 11.10.2017: Small herring and Large herring added as new species [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=cVl6fTatRsE0DwiQ]
* Model run 11.10.2017: Small herring and Large herring added as new species [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=cVl6fTatRsE0DwiQ]
* Model rerun 15.11.2017 because the previous stored run was lost in update [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=njuQREEjuMtY4MlQ]
* Model rerun 15.11.2017 because the previous stored run was lost in update [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=njuQREEjuMtY4MlQ]
* Model run 21.3.2018: Small and large herring replaced by actual fish length [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=rzvSZKM7ekjDqEhu]


<rcode name="preprocess" label="Preprocess (for developers only)" 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(OpasnetUtils)
library(OpasnetUtils)
Line 86: Line 87:
# Divide herring into small and large (as new "species")
# Divide herring into small and large (as new "species")
euRaw <- Ovariable("euRaw", ddata = "Op_en3104", subset = "POPs")
euRaw <- Ovariable("euRaw", ddata = "Op_en3104", subset = "POPs")
euRaw <- EvalOutput(euRaw)
#euRaw <- EvalOutput(euRaw)
temp <- euRaw[euRaw$Fish_species == "Baltic herring" , ]
 
temp$Fish_species <- as.factor(
#temp <- euRaw[euRaw$Fish_species == "Baltic herring" , ]
  ifelse(
#temp$Fish_species <- as.factor(
    as.numeric(as.character(temp$Length_mean_mm)) >= 170,
ifelse(
    "Large herring",
#    as.numeric(as.character(temp$Length_mean_mm)) >= 170,
    "Small herring"
#    "Large herring",
  )
#    "Small herring"
)
)
euRaw <- combine(euRaw, temp, name = "euRaw")
#)
#euRaw <- combine(euRaw, temp, name = "euRaw")


# levels(TEF$Group)
# levels(TEF$Group)
Line 108: Line 110:
   ),
   ),
   formula = function(...) {
   formula = function(...) {
     eu <- euRaw[,c(1:4, 18, 19)] # THL code, Matrix, Congener, Fish species
    euRaw$Length<-as.numeric(as.character(euRaw$Length_mean_mm))
     colnames(eu@output)[1:4] <- c("THLcode", "Matrix", "Compound", "Fish")
    euRaw$Year <- as.numeric(substr(euRaw$Catch_date, nchar(as.character(euRaw$Catch_date))-3,100))
    euRaw@marginal[colnames(euRaw@output) %in% c("Length","Year")] <- TRUE
#    euRaw$Weight<-as.numeric(as.character(euRaw$Weight_mean_g))
#    euRaw$Large <- euRaw$Length >= 170
   
     eu <- euRaw[,c(1:4, 10, 20, 21, 18)] # See below + Length, Year, Result
     colnames(eu@output)[1:5] <- c("THLcode", "Matrix", "Compound", "Fish", "N")
      
      
     temp <- oapply(eu * TEF, cols = "Compound", FUN = "sum")
     temp <- oapply(eu * TEF, cols = "Compound", FUN = "sum")
     colnames(temp@output)[colnames(temp@output)=="Group"] <- "Compound"
     colnames(temp@output)[colnames(temp@output)=="Group"] <- "Compound"
     eu <- combine(eu, temp)
     eu <- combine(eu, temp)
 
   
     eu$Compound <- factor( # Compound levels are ordered based on the data table on [[TEF]]
     eu$Compound <- factor( # Compound levels are ordered based on the data table on [[TEF]]
       eu$Compound,  
       eu$Compound,  
Line 185: Line 193:
#[16] "Large herring"  
#[16] "Large herring"  


euRaw$Length<-as.numeric(as.character(euRaw$Length_mean_mm))
#ggplot(euRaw@output, aes(x=Length, y=Weight))+geom_point()+
euRaw$Weight<-as.numeric(as.character(euRaw$Weight_mean_g))
facet_wrap(~Fish_species, scales="free")
euRaw$Large <- euRaw$Length >= 170
 
ggplot(euRaw@output, aes(x=Length, y=Weight))+geom_point()+
  facet_wrap(~Fish_species, scales="free")


#oprint(summary(
#oprint(summary(
Line 217: Line 221:
* Model run 11.10.2017 with small and large herring [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=ICIWZTUZR6rlNwuD] (removed in update)
* Model run 11.10.2017 with small and large herring [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=ICIWZTUZR6rlNwuD] (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. [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]
Warning messages:
1: In reshapeWide(data, idvar = idvar, timevar = timevar, varying = varying, :
multiple rows match for Compound=PCDDF: first taken
2: In reshapeWide(data, idvar = idvar, timevar = timevar, varying = varying, :
multiple rows match for Compound=PCB: first taken
Error in jags.model(mod, data = list(S = nrow(eu3), C = C, Fi = Fi, cong = log(cong), :
RUNTIME ERROR:
Compilation error on line 5.
Index out of range taking subset of cong


<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 234: Line 248:
conl <- indices$Compound.TEQ2
conl <- indices$Compound.TEQ2
# fisl <- indices$Fish.Fish16
# fisl <- indices$Fish.Fish16
fisl <- c("Baltic herring","Large herring","Salmon","Small herring")
# fisl <- c("Baltic herring","Large herring","Salmon","Small herring")
fisl <- c("Baltic herring","Salmon")
C <- length(conl)
C <- length(conl)
Fi <- length(fisl)
Fi <- length(fisl)
Line 253: Line 268:


eu2 <- unkeep(eu2, prevresults = TRUE, sources = TRUE)
eu2 <- unkeep(eu2, prevresults = TRUE, sources = TRUE)
eu2 <- oapply(eu2, cols = "TEFversion", FUN = "sum") # This goes wrong if > 1 TEFversion
#eu2 <- oapply(eu2, cols = "TEFversion", FUN = "sum") # This goes wrong if > 1 TEFversion. Why do we need this, it drops ca. 400 rows?


# Hierarchical Bayes model.
# Hierarchical Bayes model.
Line 277: Line 292:


#> colnames(eu3)
#> colnames(eu3)
#[1] "THLcode"        "Fish"           "euResult.PCDDF" "euResult.PCB"
#[1] "THLcode"        "Fish"          "N"              "Length"        "Year"        
#[6] "euResult.PCDDF" "euResult.PCB"


if(FALSE){
# Find the level of quantification for dinterval function
# Find the level of quantification for dinterval function
LOQ <- unlist(lapply(eu3[3:ncol(eu3)], FUN = function(x) min(x[x!=0])))
LOQ <- unlist(lapply(eu3[6:ncol(eu3)], FUN = function(x) min(x[x!=0])))  
# With TEQ, there are no zeros. So this is useful only if there are congener-specific results.
names(LOQ) <- conl
names(LOQ) <- conl
cong <- data.matrix(eu3[3:ncol(eu3)])
}
cong <- sapply(
 
  1:length(LOQ),
cong <- data.matrix(eu3[6:ncol(eu3)])
  FUN = function(x) ifelse(cong[,x]==0, 0.5*LOQ[x], cong[,x])
#cong <- sapply(
)
1:length(LOQ),
FUN = function(x) ifelse(cong[,x]==0, 0.5*LOQ[x], cong[,x])
#)


# This version of the model looks only at Baltic herring, Large herring, small herring and salmon.
# This version of the model looks only at Baltic herring, Large herring, small herring and salmon.
Line 292: Line 312:


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[,])
  cong[i,1:C] ~ dmnorm(muind[fis[i],], Omega[,])
    }
  muind[i,1:C] <- mu[fis[i],1:C] + b[fis[i]]*length[i] + b[3]*year[i]
    for(i in 1:Fi) { # Fi = fish species
  }
      for(j in 1:C) {  
 
        mu[i,j] ~ dnorm(mu1[j], tau1[j])
  # Prior for betas
      }
  b[1] ~ dnorm(0,0.0001) # length parameter herring
      pred[i,1:C] ~ dmnorm(mu[i,], Omega[,]) # Model prediction
  b[2] ~ dnorm(0,0.0001) # length parameter salmon
    }
  b[3] ~ dnorm(0,0.0001) # Time trend
    for(i in 1:C) { # C = Compound
 
      mu1[i] ~ dnorm(0, 0.0001)
  for(i in 1:Fi) { # Fi = fish species
      tau1[i] ~ dunif(0,10000)
  for(j in 1:C) {  
      pred1[i] ~ dnorm(mu1[i], tau1[i])
    mu[i,j] ~ dnorm(mu1[j], tau1[j]) # Congener-specific mean for fishes
    }
  }
    Omega[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
pred[i,1:C] ~ dmnorm(mu[i,1:C], Omega[,]) # Model prediction
    Omega1[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
  }
  for(i in 1:C) { # C = Compound
  mu1[i] ~ dnorm(0, 0.0001)
  tau1[i] ~ dunif(0,10000)
pred1[i] ~ dnorm(mu[i], tau1[i])
  }
  Omega[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
  Omega1[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
   }
   }
")
")
Line 321: Line 348:
     Fi = Fi,
     Fi = Fi,
     cong = log(cong),
     cong = log(cong),
    length = eu3$Length,
    year = eu3$Year,
     fis = match(eu3$Fish, fisl),
     fis = match(eu3$Fish, fisl),
     Omega0 = diag(C)/100000
     Omega0 = diag(C)/100000
Line 335: Line 364:
     'mu', # mean by fish and compound
     'mu', # mean by fish and compound
     'Omega', # precision matrix by fish and compound
     'Omega', # precision matrix by fish and compound
    'pred', # predicted concentration by fish and compound
'pred', # predicted concentration by fish and compound
     #    'mu1', # mean prior for mu by compound
     #    'mu1', # mean prior for mu by compound
     'Omega1', # precision matrix by compound
     'Omega1', # precision matrix by compound
     #    'tau1', # precision for prior of all mu  
     #    'tau1', # precision for prior of all mu  
    'pred1' # predicted concentration by compound
#    'pred1', # predicted concentration by compound
    'b' # parameters for length and year
   ),  
   ),  
   N
   N
)
)
dimnames(samps.j$mu) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$mu) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$b) <- list(Param = c("Lenherr","Lensalm","Time"), 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$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$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$pred1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$Omega) <- list(Compound = conl, Compound2 = conl, Iter=1:N, Chain=1:4)
dimnames(samps.j$Omega) <- list(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)
dimnames(samps.j$Omega1) <- list(Compound = conl, Compound2 = conl, Iter=1:N, Chain=1:4)
Line 355: Line 386:
   mu = apply(samps.j$mu[,,,1], MARGIN = 1:2, FUN = mean),
   mu = apply(samps.j$mu[,,,1], MARGIN = 1:2, FUN = mean),
   Omega = apply(samps.j$Omega[,,,1], MARGIN = 1:2, FUN = mean),
   Omega = apply(samps.j$Omega[,,,1], MARGIN = 1:2, FUN = mean),
   pred.mean = apply(samps.j$pred[,,,1], MARGIN = 1:2, FUN = mean),
   b = apply(samps.j$b[,,1], MARGIN = 1, FUN = mean)
  pred.sd = apply(samps.j$pred[,,,1], MARGIN = 1:2, FUN = sd),
  #  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.mean = apply(samps.j$pred1[,,1], MARGIN = 1, FUN = mean),
  pred1.sd = apply(samps.j$pred1[,,1], MARGIN = 1, FUN = sd)
pred1.sd = apply(samps.j$pred1[,,1], MARGIN = 1, FUN = sd)
)
)



Revision as of 20:17, 21 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]

Warning messages: 1: In reshapeWide(data, idvar = idvar, timevar = timevar, varying = varying, : multiple rows match for Compound=PCDDF: first taken 2: In reshapeWide(data, idvar = idvar, timevar = timevar, varying = varying, : multiple rows match for Compound=PCB: first taken Error in jags.model(mod, data = list(S = nrow(eu3), C = C, Fi = Fi, cong = log(cong), : RUNTIME ERROR: Compilation error on line 5. Index out of range taking subset of cong

+ Show code

Initiate conc_pcddf

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

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