EU-kalat

[show] This page is a study. The page identifier is Op_en3104
Moderator:Arja (see all)
Upload data {{#opasnet_base_link:Op_en3104}}

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. .

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) (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 [15]

+ 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]] conl <- indices$Compound.TEQ2 fisl <- indices$Fish.Fish14 C <- length(conl) Fi <- length(fisl) N <- 1000 conl fisl eu <- EvalOutput(eu) 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(eu$Compound)[replaces[[i]][1]] <- replaces[[i]][2] } eu <- oapply(eu, INDEX = "THLcode", FUN = "sum") # Hierarchical Bayes model. # PCDD/F concentrations in fish. # It uses the sum of PCDD/F (Pcdsum) as the total concentration of dioxin in fish. # Cong_j is the fraction of a congener from pcdsum. # pcdsum is log-normally distributed. cong_j follows Dirichlet distribution. # pcdsum depends on age of fish, fish species and catchment area, but we only have species now so other variables are omitted. # cong_j depends on fish species. if(FALSE) { # Multicongener model NOT updated to reflect new names: eu -> euMain 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 fishsamples <- 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(fishsamples[3:ncol(fishsamples)], FUN = function(x) min(x[x!=0]))) names(LOQ) <- conl cong <- data.matrix(fishsamples[3:ncol(fishsamples)]) 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(fishsamples), C = C, Fi = Fi, cong = log(cong), # LOQ = LOQ, # below.LOQ = is.na(cong)*1, fis = match(fishsamples$Fish, fisl), Omega0 = diag(C)/100000 ), n.chains = 4, n.adapt = 100 ) } # if(FALSE) fishsamples <- unkeep(eu, prevresults = TRUE, sources = TRUE)@output fishsamples <- fishsamples[fishsamples$Compound %in% conl & fishsamples$Matrix == "Muscle" , ] fishsamples <- reshape( fishsamples, v.names = "euResult", idvar = "THLcode", timevar = "Compound", drop = c("Matrix"), direction = "wide" ) # colnames(eu@output) # "THLcode" "Matrix" "Compound" "Fish" "euRawSource" # [6] "TEFversion" "TEFrawSource" "TEFSource" "Source" "euResult" # [11] "euSource" oprint(head(fishsamples)) # Find the level of quantification for dinterval function LOQ <- unlist(lapply(fishsamples[3:ncol(fishsamples)], FUN = function(x) min(x[x!=0]))) names(LOQ) <- conl cong <- data.matrix(fishsamples[3:ncol(fishsamples)]) 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(fishsamples), C = C, Fi = Fi, cong = log(cong), # LOQ = LOQ, # below.LOQ = is.na(cong)*1, fis = match(fishsamples$Fish, fisl), Omega0 = diag(C)/100000 ), n.chains = 4, n.adapt = 100 ) update(jags, 100) samps.j <- jags.samples( jags, c('mu', 'Omega', 'pred', 'mu1', 'Omega1', 'tau1', 'pred1'), 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$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$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") euRatio <- EvalOutput(euRatio) # Leave only the main fish species and congeners and remove others eu[eu$Compound %in% conl & eu$Fish %in% fisl , ] ggplot(eu@output, aes(x = euResult, colour = Fish))+geom_density()+ facet_wrap(~ Compound) + scale_x_log10() ggplot(euRatio@output, aes(x = euRatioResult, colour = Fish))+geom_density()+ facet_wrap(~ Compound) # 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=Congener))+geom_density()+ facet_wrap( ~ Fish, scales = "free_y")+scale_x_log10() #stat_ellipse() coda.j <- coda.samples( jags, c('mu', 'pred', 'omega1', 'pred1'), N ) plot(coda.j)

Initiate concentration

• Model run 19.5.2017 [16]

+ 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")