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, 1/sigmaOfLogTCDD^2 )
                    }                   
                     sigmaOfLogTCDD ~ dunif( 0.001*sdOfLogTCDD , 1000*sdOfLogTCDD )
                     muOfLogTCDD ~ 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)

dataList = list(
    TCDD = TCDD,
    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, 10000)

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

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