EU-kalat

This page is a study. The page identifier is Op_en3104
Moderator:Arja (see all)
Give your opinion to the peer rating of the content of this page. {{ #opasnet_rater: }}
Upload data Show results

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]

+ Show code- Hide code # This is code Op_en3104/preprocess on page [[EU-kalat]] library(rjags) library(OpasnetUtils) library(ggplot2) library(reshape2) eu <- Ovariable("eu", 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)] colnames(eu@output) <- c("THLcode", "Matrix", "Congener", "Fish", "euResult", "euSource") #> unique(eu@output$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 with too little data (>70% BDL) and all non-PCDDF # aggregate(eu@data$euResult, by = eu@data["POP"], FUN = function(x) mean(x == 0)) conl <- unique(eu@output$Congener)[c(1:12, 14, 15)] # 7 OCDD should be removed eu@output <- eu@output[eu@output$Congener %in% conl , ] #[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 fisl <- unique(eu@output$Fish)[c(1:4, 6:14, 17)] eu@output <- eu@output[eu@output$Fish %in% fisl , ] # Remove four with too little data eut <- eu eut@output <- eut@output[ eut@output$Congener == "2378TCDD" & eut@output$Matrix == "Muscle" & result(eut) != 0 , ] # Zeros cannot be used in ratio estimates eut@output$Congener <- NULL eut@marginal[colnames(eut@output) == "Congener"] <- NULL eut <- log10(eu / eut)@output conl <- conl[conl != "2378TCDD"] eut <- eut[eut$Congener %in% conl , ] # 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. ggplot(eu@output, aes(x = euResult, colour = Fish))+geom_density()+ facet_wrap(~ Congener) + scale_x_log10() ggplot(eut, aes(x = Result, colour = Fish))+geom_density()+ facet_wrap(~ Congener) objects.store(eu,eut) cat("Ovariable eu and data.frame eut stored.\n")

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]

+ Show code- Hide code # This is code Op_en3014/bayes on page [[EU-kalat]] library(OpasnetUtils) library(reshape2) library(rjags) objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]] # 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. conl <- as.character(unique(eu@output$Congener)) fisl <- sort(as.character(unique(eu@output$Fish))) 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[cong == 0] <- 0.001 # NA # Needed for dinterval mod <- textConnection(" model{ 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]) cong[i,1:C] ~ dmnorm(mu[fis[i],], Omega[fis[i],,]) # } } for(i in 1:F) { # F = fish species mu[i,1:C] ~ dmnorm(mu0[1:C], S2[1:C,1:C]) Omega[i,1:C,1:C] ~ dwish(S3[1:C,1:C],S) ans.pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction } } ") C <- length(conl) jags <- jags.model( mod, data = list( S = nrow(fishsamples), C = C, F = length(fisl), cong = log(cong), # LOQ = LOQ, # below.LOQ = is.na(cong)*1, fis = match(fishsamples$Fish, fisl), mu0 = rep(2,C), S2 = 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.adapt = 100 ) } update(jags, 100) samps.j <- jags.samples( jags, c('mu', 'Omega', 'ans.pred'), 1000 ) 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") ) ) # Mu for all congeners of fish1 scatterplotMatrix(t(js[,1,,1])) # All parameters for TCDD of fish1 scatterplotMatrix(js[1,1,,]) # 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( "jsp", dependencies = data.frame(Name = "fish.param"), formula = function(...) { jsp <- lapply( 1:2, FUN = function(x) { mvrnorm(openv$N, fish.param$mu[,x], fish.param$Omega[[x]]) } ) jsp <- rbind( cbind( Fish = "Salmon", 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) } ) # Combine modelled survey answers with estimated amounts and frequencies by: # Rounding the modelled result and merging that with value in ovariable assump often <- Ovariable( "often", dependencies = data.frame(Name=c("jsp","assump")), formula = function(...) { 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) 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) #temp <- adply(samps$mu, c(1,2,3,4)) pcd.pred <- array( samps$pcd.pred, dim = c(length(fisl), length(conl), 1000, 4), dimnames = list( 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) #cat("Arrays pcd.pred, mu.pred stored.\n")