Difference between revisions of "EU-kalat"

# This code is Op_en on page [[EU-kalat]].

library(OpasnetUtils)
library(ggplot2)

eu <- Ovariable("eu", ddata = "Op_en3104", subset = "POPs")
eu <- EvalOutput(eu)
#levels(eu$POP)
eu <- eu[eu$Fish_species == "Baltic herring" , ]

tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values")
tef <- EvalOutput(tef)
colnames(eu@output)[colnames(eu@output) == "POP"] <- "Congener"
levels(eu$Congener) <- gsub("HCDD", "HxCDD", levels(eu$Congener))
levels(eu$Congener) <- gsub("HCDF", "HxCDF", levels(eu$Congener))
levels(eu$Congener) <- gsub("CoPCB", "PCB", levels(eu$Congener))
#levels(eu$Congener)

euteq <- eu * tef
#head(euteq@output)
temp <- aggregate(euteq@output["Result"], by = c(euteq@output["Congener"], euteq@output["Catch_location"]), FUN = mean)

temp2 <- temp
temp2$Result <- temp2$Result / sum(temp2$Result)

cat("Congener TEQ profile of Baltic herring.\n")
oprint(temp2)
 
ggplot(temp, aes(x = Congener, y = Result, fill = Congener)) + geom_bar(stat = "identity") + facet_wrap(~ Catch_location) + 
labs(x = "Congener", y = "TEQ concentration (pg/g fat)") + theme_gray(base_size = 18) + coord_flip()

# This is code Op_en3104/ on page [[EU-kalat]]
library(rjags)
library(OpasnetUtils)
library(ggplot2)
library(reshape2)

eu <- Ovariable("eu", ddata = "Op_en3104", subset = "POPs")
eu <- EvalOutput(eu)
eu@output <- eu@output[c(1:4, 18, 19)] # THL code, Matrix, Congener, Fish species
eu@marginal <- eu@marginal[c(1:4, 18, 19)]
eu@output <- eu@output[result(eu) != 0 , ] # Zeros cannot be used in ratio estimates
colnames(eu@output)[colnames(eu@output) == "POP"] <- "Congener"
colnames(eu@output)[colnames(eu@output) == "Fish_species"] <- "Fish"

#> unique(eu@output$Congener)
#[1] 2378TCDD     12378PeCDD   123478HCDD   123678HCDD   123789HCDD   1234678HpCDD
#[7] OCDD         2378TCDF     12378PeCDF   23478PeCDF   123478HCDF   123678HCDF  
#[13] 123789HCDF   234678HCDF   1234678HpCDF 1234789HpCDF OCDF         CoPCB77 ...     

# Remove the four with too little data (>70% BDL) and all non-PCDDF
# aggregate(eu@data$euResult, by = eu@data["POP"], FUN = function(x) mean(x == 0))
conl <- unique(eu@output$Congener)[c(1:12, 14, 15)] # 7 OCDD should be removed
eu@output <- eu@output[eu@output$Congener %in% conl , ] 

#[1] Baltic herring Sprat          Salmon         Sea trout      Vendace       
#[6] Roach          Perch          Pike           Pike-perch     Burbot        
#[11] Whitefish      Flounder       Bream          River lamprey  Cod           
#[16] Trout          Rainbow trout  Arctic char   

fisl <- unique(eu@output$Fish)[c(1:4, 6:14, 17)] 
eu@output <- eu@output[eu@output$Fish %in% fisl , ] # Remove four with too little data

eut <- eu
eut@output <- eut@output[
  eut@output$Congener == "2378TCDD" & eut@output$Matrix == "Muscle"  , ]
eut@output$Congener <- NULL
eut@marginal[colnames(eut@output) == "Congener"] <- NULL
eut <- log10(eu / eut)@output
conl <- conl[conl != "2378TCDD"]
eut <- eut[eut$Congener %in% conl , ]

# Analysis: a few rows disappear here, as shown by numbers per fish species. Why?
# aggregate(eut@output, by = eut@output["Fish"], FUN = length)[[2]]
# [1] 1181  141  131   55  158   51   96   30   36   40  102   49   61  180 eu
# [1] 1173  141  131   55  158   51   96   27   36   40  102   49   53  180 eu/eut
# There are samples where TCDD concenctration is zero.
#eu@data[eu@data$POP == "2378TCDD" & eu@data$euResult == 0 , c(1,4)][[1]]
#eu@data[eu@data[[1]] %in% c("09K0585", "09K0698", "09K0740", "09K0748") , c(1,3,4,18)]
# Some rows disappear because there is no TCDD measurement to compare with:
#8 Baltic herrings, 8 sprats, 3 rainbow trouts.
# Conclusion: this is ok. Total 2292 rows.

ggplot(eut, aes(x = Result, colour = Fish))+geom_density()+
  facet_wrap(~ Congener) 

modf <- textConnection("
  model{
    for(i in 1:N) { # i = observation
      eutt[i] ~ dnorm(mu[con[i],fis[i]], tau[con[i],fis[i]])
    }
    for(j in 1:J) { # j = congener
      for(k in 1:K) { # k = fish species
        mu[j,k] ~ dunif(-3,3)
        tau[j,k] <- pow(sigma[j,k], -2)
        sigma[j,k] ~ dunif(0, 10)
        pred[j,k] ~ dnorm(mu[j,k], tau[j,k]) # Model prediction
      }
    }
  }
")

# A binomial distribution is assumed for bins of answer choices.
jags <- jags.model(
  modf,
  data = list(
    N = nrow(eut),
    J = length(conl),
    K = length(fisl),
    eutt = eut$Result,
    con = match(eut$Congener, conl),
    fis = match(eut$Fish, fisl)
  ),
  n.chains = 4,
  n.adapt = 100
)

update(jags, 1000)

samps <- jags.samples(jags, c('mu', 'pred'), 1000)
pl <- melt(array(samps$pred, dim = c(length(conl), length(fisl), 1000, 4), 
  dimnames = list(
  Congener = conl, 
  Fish = fisl,
  Iter = 1:1000,
  Seed = 1:4
  )
))

ggplot(pl, aes(x = value, colour = Fish))+geom_density(size = 1) + 
  facet_wrap(~ Congener, scale = "free_y") + coord_cartesian(xlim = c(-1.5,1.5))

ggplot(pl[pl$Fish %in% c("Baltic herring", "Rainbow trout", "Pike-perch", "Pike", "Salmon", "Sprat") , ], aes(x = value, colour = Fish))+geom_density(size = 1) + 
  facet_wrap(~ Congener, scale = "free_y") + coord_cartesian(xlim = c(-1.5,1.5))

samps.coda <- coda.samples(jags, c('mu', 'pred'), 1000)

out <- data.frame(
  Congener = conl,
  Fish = rep(fisl, each = length(conl)),
  Param = rep(c("mu", "pred"), each = length(conl)*length(fisl)),
  Value = (summary(samps.coda)[[1]][,1])
)

oprint(out)

ggplot(out[out$Param == "mu", ], aes(x = Congener, y = Value, colour = Fish, group = Fish)) + 
  geom_line() + labs(title = "Congener amounts in fish compared with TCDD", y = "Log_10 ratio") 

if(FALSE){
  tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values")
  tef <- EvalOutput(tef)
  
  #levels(eu$Congener) <- gsub("HCDD", "HxCDD", levels(eu$Congener))
  #levels(eu$Congener) <- gsub("HCDF", "HxCDF", levels(eu$Congener))
  #levels(eu$Congener) <- gsub("CoPCB", "PCB", levels(eu$Congener))
  euteq <- eu * tef
} # IF FALSE