|
|
| Line 117: |
Line 117: |
| | resultsall <- EvalOutput(resultsall) | | resultsall <- EvalOutput(resultsall) |
| | teq = resultsall * tef | | teq = resultsall * tef |
| − | teqpart <- teq@output
| |
| − | teqpart <- teqpart[teqpart$variable == "Meanpost", ]
| |
| | | | |
| | objects.store(resultsall, teqpart) | | objects.store(resultsall, teqpart) |
| − | cat("resultsall and teqpart stored for later use:\n", paste(ls(), collapse = ", "), "\n") | + | cat("resultsall and teq stored for later use:\n", paste(ls(), collapse = ", "), "\n") |
| | | | |
| | </rcode> | | </rcode> |
Revision as of 10:39, 17 June 2016
Question
What are the concentrations of persistent organic pollutants (POPs) in Baltic herring.
Answer
Answer is under work.
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 not shown yet).
Based on the mean concentrations of individual congeners, TEQs are calculated for each congener and plotted below.
Error creating thumbnail: Unable to save thumbnail to destination
Rationale
+ 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
objects.store(resultsall, teqpart)
cat("resultsall and teq stored for later use:\n", paste(ls(), collapse = ", "), "\n")
| |
See also
References
