|
|
Line 29: |
Line 29: |
| objects.latest("Op_en2583", "pop_bayes") | | objects.latest("Op_en2583", "pop_bayes") |
| | | |
− | ggplot(subset(resultsall@output, grepl("PCB", Congener) & grepl("Meanpost", variable)), aes(x = Congener, y = resultsallResult, fill = Congener)) + geom_bar(stat = "identity") + | + | ggplot(subset(resultsall@output, grepl("PCB", Congener) & grepl("Medianpost", variable)), aes(x = Congener, y = resultsallResult, fill = Congener)) + geom_bar(stat = "identity") + |
− | labs(x = "Congener", y = "Concentration (pg/g fat)") + coord_flip() + ggtitle("Mean concentrations of PCB's in Baltic herring") | + | labs(x = "Congener", y = "Concentration (pg/g fat)") + coord_flip() + ggtitle("Median posterior concentrations of PCB's in Baltic herring") |
− | | |
− | ggplot(subset(resultsall@output, grepl("CD", Congener) & grepl("Meanpost", variable)), aes(x = Congener, y = resultsallResult, fill = Congener)) + geom_bar(stat = "identity") +
| |
− | labs(x = "Congener", y = "Concentration (pg/g fat)") + coord_flip() + ggtitle("Mean concentrations of PCDD/F's in Baltic herring")
| |
− | | |
− | ggplot(subset(resultsall@output, grepl("BDE", Congener) & grepl("Meanpost", variable)), aes(x = Congener, y = resultsallResult, fill = Congener)) + geom_bar(stat = "identity") +
| |
− | labs(x = "Congener", y = "Concentration (pg/g fat)") + coord_flip() + ggtitle("Mean concentrations of BDE's in Baltic herring")
| |
− | | |
− | ggplot(subset(teq@output, grepl("Meanpost", variable)), aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") +
| |
− | labs(x = "Congener", y = "TEQ (pg/g fat)") + coord_flip() + ggtitle("Mean TEQ of POPs in Baltic herring")
| |
− | | |
| | | |
| + | ggplot(subset(resultsall@output, grepl("CD", Congener) & grepl("Medianpost", variable)), aes(x = Congener, y = resultsallResult, fill = Congener)) + geom_bar(stat = "identity") + |
| + | labs(x = "Congener", y = "Concentration (pg/g fat)") + coord_flip() + ggtitle("Median posterior concentrations of PCDD/F's in Baltic herring") |
| | | |
| + | ggplot(subset(resultsall@output, grepl("BDE", Congener) & grepl("Medianpost", variable)), aes(x = Congener, y = resultsallResult, fill = Congener)) + geom_bar(stat = "identity") + |
| + | labs(x = "Congener", y = "Concentration (pg/g fat)") + coord_flip() + ggtitle("Median posterior concentrations of BDE's in Baltic herring") |
| | | |
| + | ggplot(subset(teq@output, grepl("Medianpost", variable)), aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") + |
| + | labs(x = "Congener", y = "TEQ (pg/g fat)") + coord_flip() + ggtitle("Median posterior TEQ of POPs in Baltic herring") |
| | | |
| </rcode> | | </rcode> |
Revision as of 12:11, 21 June 2016
Question
What are the concentrations of persistent organic pollutants (POPs) in Baltic herring.
Answer
Answer is under work and results are preliminary.
POP concentrations in Baltic sea fish have been measured from samples collected in EU-kalat project. The original data of individual fish samples si accessible through Opasnet base. This data is used here for a Bayesian model to calculate posterior distributions (mean and SD) for each congener. This data is then translated into TEQ to be used for health impact assessment of Baltic herring.
Posterior congener mean concentrations are presented below for each compound group (PCDD/F, PCB, BDE) taken under analysis in EU-kalat.
Error creating thumbnail: Unable to save thumbnail to destination
Error creating thumbnail: Unable to save thumbnail to destination
Error creating thumbnail: Unable to save thumbnail to destination
Based on the mean posterior concentrations of individual congeners, TEQs are calculated for each congener by using TEF values by WHO and plotted below.
Error creating thumbnail: Unable to save thumbnail to destination
Rationale
This model takes in measured congener concentrations of POPs in Baltic herring in northern bart of Baltic sea. Measured data is used for Bayesian model to produce posterior medians and sds for each congener and also to calculate TEQ values. Numerical results are saved as variables to Opasnetbase and result figures are presented above in the Answer section.
+ 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
Medianorig = aggregate(compdat$Result, compdat["POP"], median)$x, #calculate median 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,
Medianpost = 10^Meanlogpost-1E-02,
Sdpost = 10^Sdlogpost-1E-02 #this might be incorrect
)
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
objects.store(resultsall, teq)
cat("resultsall and teq stored for later use:\n", paste(ls(), collapse = ", "), "\n")
| |
See also
References