Difference between revisions of "POPs in Baltic herring"

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

 

for (j in 1 : 3 ) { # for testing reduced to 3 congeners, real loop amount is Y (=79)

  compdat <- dat[dat$POP == congeners[j] , ] #data for current congener

  mo <- textConnection("model{

                    for( i in 1 : N ) {

                     Compound[i] ~ dnorm(muOfCompound, tau2[i])

                     tau2[i] <- tau1*sqrt(n[i])

                    }                  

                     tau1 ~ dunif(0.001, 1000)

                     muOfCompound ~ dnorm(0, 0.001)

                            }

                     ")

  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+1) #+1 because of negative log values

  }
  else Compound <- compdat$result
  dataList = list(
    Compound = Compound,
    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('tau2', 'muOfCompound'), 500) # Stores a posterior sample
  plot(out)

  #filename <- paste("out", j, sep="") #save posterior for each congener
  #assign(filename, out)

  meanorig <- mean(compdat$result) #calculate mean for original data
  sdorig <- sd(compdat$result) #calcaulte sd of original data
  meanlog <- mean(Compound) #calculate mean for logdata
  sdlog <- sd(Compound) #calculate sd for logdata
  meanpost <- mean(data.frame(out[[4]])$muOfCompound) #calculate mean of mu for posterior (test 4)
  #sdpost <- mean(data.frame(out[[4]])$tau2) #if needed, to be figured out how to do this

  results <- data.frame(congeners[j], meanorig, sdorig, meanlog, meanpost) #results for congener

  resultsall <- rbind(resultsall, results[, 1:5]) #save results of congener to combined results table
}

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