Question
What are the levels of persistent organic pollutants (POPs) in Baltic sea salmon?
Answer
POP concentrations in Baltic sea fish have been measured from samples collected in EU-kalat project.
+ Show code- Hide code
library(OpasnetUtils)
library(ggplot2)
library(rjags)
dat <- opbase.data("Op_en3104", subset = "POPs")
dat <- dat[dat$Fish_species == "Salmon" , ]
congeners <- levels(dat$POP) #names of different congeners in data
Y <- length(congeners) #number of congeners and j's in for loop
#resultsall <-c('congener', 'meanorig', 'sdorig', 'meanlog', 'meanpost') #comined results table structure to store results of individual congener runs
compdat <- dat[dat$POP %in% congeners[1:3] , ] #data for current congener
Compound <- log(compdat$Result+1E-2) #+1E-2 because zero concentrations are not allowed
mo <- textConnection("model{
for (j in 1 : 3 ) { # for testing reduced to 3 congeners, real loop amount is Y (=79)
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])
}
}
")
# shap <- shapiro.test(compdat$Result) #test of normality
#filename <- paste("shap", j, sep="") #save shapiro test result for each congener, just for fun
#assign(filename, shap)
# if(shap$p.value < 0.05) { #if skewed, logtransformation, if not continue with original data. Or should the logtransformation be done in any case?
#Compound <- log(compdat$Result+1E-2) #+1E-2 because zero concentrations are not allowed
# } else {
# Compound <- compdat$Result
# }
dataList = list(
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)
#filename <- paste("out", j, sep="") #save posterior for each congener
#assign(filename, out)
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
Meanpost = c(mean(out[[4]][[1]]), #calculate mean of mu for posterior (test 4)
mean(out[[4]][[2]]),
mean(out[[4]][[3]])
),
Sdpost = c(mean(out[[4]][[4]]), #calculate mean of mu for posterior (test 4)
mean(out[[4]][[5]]),
mean(out[[4]][[6]])
)
)
oprint(resultsall)
ggplot(dat, aes(x = Result, colour = Catch_location))+geom_density()+scale_x_log10()
| |
Rationale
Dependencies
Formula
See also
Keywords
References
Related files
<mfanonymousfilelist></mfanonymousfilelist>
- POPs in Baltic salmon. Opasnet . [1]. Accessed 20 May 2024.
