+ Show code- Hide code
# 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
| |