|
|
Line 95: |
Line 95: |
| | | |
| == See also == | | == See also == |
| + | |
| + | * [http://www.openbugs.net/Manuals/Tricks.html#HandlingUnbalancedDatasets Tricks for OpenBUGS] |
| | | |
| == References == | | == References == |
| | | |
| <references /> | | <references /> |
Revision as of 13:40, 30 March 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
library(OpasnetUtils)
library(ggplot2)
library(rjags)
dat <- opbase.data("Op_en3104", subset = "POPs")
dat <- dat[dat$Fish_species == "Baltic herring" , ]
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()
| |
See also
References