|
|
Line 12: |
Line 12: |
| POP concentrations in Baltic sea fish have been measured from samples collected in [[EU-kalat]] project. | | POP concentrations in Baltic sea fish have been measured from samples collected in [[EU-kalat]] project. |
| | | |
− | <rcode graphics=1> | + | <rcode name="pop_bayes" label="Calculate (for developers only)" graphics=1 store=1> |
| + | ## This code is ???/pop_bayes on page [[POPs_in_Baltic_herring]] |
| + | |
| library(OpasnetUtils) | | library(OpasnetUtils) |
| library(ggplot2) | | library(ggplot2) |
Line 59: |
Line 61: |
| #update(jags, 1000) | | #update(jags, 1000) |
| out <- coda.samples(jags, c('tau1', 'muOfCompound'), 500) # Stores a posterior sample | | out <- coda.samples(jags, c('tau1', 'muOfCompound'), 500) # Stores a posterior sample |
− | #plot(out)
| + | plot(out) |
| | | |
| Meanlogpost = c() | | Meanlogpost = c() |
Line 84: |
Line 86: |
| ) | | ) |
| | | |
− | #oprint(resultsall)
| + | oprint(resultsall) |
| | | |
| tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values") | | tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values") |
Line 98: |
Line 100: |
| teqpart <- teqpart[teqpart$variable == "Meanpost", ] | | teqpart <- teqpart[teqpart$variable == "Meanpost", ] |
| | | |
− | ggplot(teqpart, aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") +
| + | objects.store(resultsall, teqpart) |
− | labs(x = "Congener", y = "TEQ (pg/g fat)") + coord_flip()
| + | cat("resultsall and teqpart stored for later use:\n", paste(ls(), collapse = ", "), "\n") |
| | | |
| </rcode> | | </rcode> |
Revision as of 10:21, 16 June 2016
Question
What are the concentrations of persistent organic pollutants (POPs) in Baltic sea fish.
Answer
Rationale
POP concentrations in Baltic sea fish have been measured from samples collected in EU-kalat project.
+ Show code- Hide code
## This code is ???/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