Difference between revisions of "POPs in Baltic salmon"

## This code is Op_en????/pop_bayes on page [[POPs_in_Baltic_salmon]]

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

dat <- opbase.data("Op_en3104", subset = "POPs")
dat <- dat[dat$Fish_species == "Salmon" , ]

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
  Medianorig = aggregate(compdat$Result, compdat["POP"], median)$x, #calculate median 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,
  Medianpost = 10^Meanlogpost-1E-02,
  Sdpost = 10^Sdlogpost-1E-02 #this might be incorrect
  )

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

objects.store(resultsall, teq)
cat("resultsall and teq stored for later use:\n", paste(ls(), collapse = ", "), "\n")