Difference between revisions of "POPs in Baltic herring"

library(OpasnetUtils)
library(ggplot2)
library(rjags)

dat <- opbase.data("Op_en3104", subset = "POPs")
dat <- dat[dat$POP == "2378TCDD" , ]
dat <- dat[dat$Fish_species == "Baltic herring" , ]

#Bayesian test based on http://doingbayesiandataanalysis.blogspot.fi/2013/04/bayesian-estimation-of-log-normal.html

mo <- textConnection("model{
                     for( i in 1 : N ) {
                     TCDD[i] ~ dlnorm(muOfLogTCDD, tau)#1/sigmaOfLogTCDD^2 )
                     }                   
                      tau ~ dunif(0.001, 1000)
                      #sigmaOfLogTCDD ~ dunif( 0.001, 1000 ) #*sdOfLogTCDD , 1000*sdOfLogTCDD )
                     muOfLogTCDD ~ dnorm(0, 0.001) # dnorm(meanOfLogTCDD , (0.001*(1/sdOfLogTCDD^2)))
                     #muOfTCDD <- (exp(muOfLogTCDD +sigmaOfLogTCDD^2/2))
                     #modeOfTCDD <- (exp(muOfLogTCDD-sigmaOfLogTCDD^2))
                     #sigmaOfTCDD <- (sqrt(exp(2*muOfLogTCDD+sigmaOfLogTCDD^2)*(exp(sigmaOfLogTCDD^2)-1)))
                     }
                     ")
#N <- nrow(dat)
#LogTCDD <- log(dat$result+1) #+1 because of negative log values
#TCDD <- 10^(LogTCDD-1)
  TCDD <- dat$Result + 0.0001
  N <- length(TCDD)

dataList = list(
  TCDD = dat$Result + 0.0001,
  N = length(TCDD)
#  meanOfLogTCDD = mean(LogTCDD),
#  sdOfLogTCDD = sd(LogTCDD)
)

jags <- jags.model(mo, data = dataList, n.chains = 4, n.adapt = 1000)
close(mo)
update(jags, 1000)

out <- coda.samples(jags, c('tau', 'muOfLogTCDD'), nrow(dat)) # Stores a posterior sample
plot(out)

ggplot(dat, aes(x = Result, colour = Catch_location))+geom_density()+scale_x_log10()