Question
What are the concentrations of persistent organic pollutants (POPs) in Baltic sea fish.
Answer
Distributions (mean and SD) of congener concentrations are calculated with a Bayes model and presented in a table below (not shown yet). Based on the mean concentrations of individual congeners, TEQs are calculated for each congener and plotted below.
Rationale
POP concentrations in Baltic sea fish have been measured from samples collected in EU-kalat project.
+ Show code- Hide code
## This code is Op_en2583/pop_bayes on page [[POPs_in_Baltic_herring]]
library(OpasnetUtils)
library(ggplot2)
library(rjags)
library(reshape2)
dat <- opbase.data("Op_en3104", subset = "POPs")
dat <- dat[dat$Fish_species == "Baltic herring" , ]
dat <- subset(dat,!(is.na(dat["Result"])))
dat <- dropall(dat)
levels(dat$POP) <- gsub("HCDD", "HxCDD", levels(dat$POP))
levels(dat$POP) <- gsub("HCDF", "HxCDF", levels(dat$POP))
levels(dat$POP) <- gsub("CoPCB", "PCB", levels(dat$POP))
congeners <- levels(dat$POP) #names of different congeners in data
Y <- length(congeners) #number of congeners and j's in for loop
compdat <- dat[dat$POP %in% congeners[1:Y] , ] #data for current congener
Compound <- log10(compdat$Result+1E-2) #+1E-2 because zero concentrations are not allowed
mo <- textConnection("model{
for (j in 1 : Y ) {
tau1[j] ~ dunif(0.001, 1000)
muOfCompound[j] ~ dnorm(0, 0.001)
}
for( i in 1 : N ) {
Compound[i] ~ dnorm(muOfCompound[POP[i]], tau2[i])
tau2[i] <- tau1[POP[i]]*sqrt(n[i])
}
}
")
dataList = list(
Y = Y,
Compound = Compound,
POP = as.numeric(compdat$POP),
N = length(Compound),
n = compdat$N_individuals #number of fishes in sample
)
jags <- jags.model(mo, data = dataList, n.chains = 4, n.adapt = 1000)
close(mo)
#update(jags, 1000)
out <- coda.samples(jags, c('tau1', 'muOfCompound'), 500) # Stores a posterior sample
plot(out)
Meanlogpost = c()
for (j in 1 : Y) {
logmean <- mean(out[[4]][,j]) #calculate mean of logmu for posterior (test 4)
Meanlogpost = c(Meanlogpost, logmean)
}
Sdlogpost = c()
for (j in 1 : Y) {
logsd <- sqrt(1/(mean(out[[4]][,j+Y])))
Sdlogpost = c(Sdlogpost, logsd)
}
resultsall <- data.frame(
Meanorig = aggregate(compdat$Result, compdat["POP"], mean), #calculate mean for original data
Sdorig = aggregate(compdat$Result, compdat["POP"], sd)$x, #calcaulte sd of original data
Meanlog = aggregate(Compound, compdat["POP"], mean)$x, #calculate mean for logdata
Sdlog = aggregate(Compound, compdat["POP"], sd)$x, #calculate sd for logdata
Meanlogpost = Meanlogpost,
Sdlogpost = Sdlogpost,
Meanpost = 10^Meanlogpost-1E-02,
Sdpost = 10^Sdlogpost-1E-02
)
oprint(resultsall)
tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values")
tef <- EvalOutput(tef)
colnames(resultsall)[1] <- "Congener"
resultsall <- melt(resultsall)
colnames(resultsall)[3] <- "Result"
resultsall <- Ovariable("resultsall", data = resultsall)
resultsall <- EvalOutput(resultsall)
teq = resultsall * tef
teqpart <- teq@output
teqpart <- teqpart[teqpart$variable == "Meanpost", ]
objects.store(resultsall, teqpart)
cat("resultsall and teqpart stored for later use:\n", paste(ls(), collapse = ", "), "\n")
| |
See also
References
