|
|
Line 7: |
Line 7: |
| | | |
| == Answer == | | == 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. | + | 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. |
| + | |
| + | [[File:Mean_TEQ_herring.jpg|800px]] |
| + | |
| | | |
| <rcode label="Show results" embed=1 graphics=1 > | | <rcode label="Show results" embed=1 graphics=1 > |
Line 15: |
Line 20: |
| | | |
| objects.latest("Op_en2583", "pop_bayes") | | objects.latest("Op_en2583", "pop_bayes") |
| + | |
| ggplot(teqpart, aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") + | | ggplot(teqpart, aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") + |
− | labs(x = "Congener", y = "TEQ (pg/g fat)") + coord_flip() | + | labs(x = "Congener", y = "TEQ (pg/g fat)") + coord_flip() + ggtitle("Mean TEQ in Baltic herring") |
| | | |
| </rcode> | | </rcode> |
Revision as of 11:01, 16 June 2016
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.
Error creating thumbnail: Unable to save thumbnail to destination
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