+ Show code- Hide code # This is code Op_en3104/bayes on page [[EU-kalat]]
library(OpasnetUtils)
library(reshape2)
library(rjags) # JAGS
library(ggplot2)
library(MASS) # mvrnorm
library(car) # scatterplotMatrix
objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]] eu, euRatio, indices
eu2 <- eu <- EvalOutput(eu)
conl <- indices$Compound.TEQ2
fisl <- indices$Fish.Fish14
C <- length(conl)
Fi <- length(fisl)
N <- 1000
conl
fisl
replaces <- list(
c("Chlorinated dibenzo-p-dioxins", "TEQdx"),
c("Chlorinated dibenzofurans", "TEQdx"),
c("Mono-ortho-substituted PCBs", "TEQpcb"),
c("Non-ortho-substituted PCBs", "TEQpcb")
)
for(i in 1:length(replaces)) {
levels(eu2$Compound)[levels(eu2$Compound)==replaces[[i]][1]] <- replaces[[i]][2]
}
eu2 <- unkeep(eu2, prevresults = TRUE, sources = TRUE)
eu2 <- oapply(eu2, cols = "TEFversion", FUN = "sum") # This goes wrong if > 1 TEFversion
# Hierarchical Bayes model.
# PCDD/F concentrations in fish.
# It uses the TEQ sum of PCDD/F (TEQdx) as the total concentration
# of dioxin and TEQpcb respectively for PCB in fish.
# TEQdx depends on age of fish, fish species and catchment area,
# but we only have species now so other variables are omitted.
# cong depends on fish species.
eu3 <- eu2[eu2$Compound %in% conl & eu2$Fish %in% fisl & eu2$Matrix == "Muscle" , ]
eu3 <- reshape(
eu3@output,
v.names = "euResult",
idvar = "THLcode",
timevar = "Compound",
drop = c("Matrix"),
direction = "wide"
)
oprint(head(eu3))
#> colnames(eu3)
#[1] "THLcode" "Fish" "euResult.TEQdx" "euResult.TEQpcb"
# Find the level of quantification for dinterval function
LOQ <- unlist(lapply(eu3[3:ncol(eu3)], FUN = function(x) min(x[x!=0])))
names(LOQ) <- conl
cong <- data.matrix(eu3[3:ncol(eu3)])
cong <- sapply(
1:length(LOQ),
FUN = function(x) ifelse(cong[,x]==0, 0.5*LOQ[x], cong[,x])
)
mod <- textConnection("
model{
for(i in 1:S) { # s = fish sample
# below.LOQ[i,j] ~ dinterval(-cong[i,j], -LOQ[j])
cong[i,1:C] ~ dmnorm(mu[fis[i],], Omega[fis[i],,])
}
for(i in 1:Fi) { # Fi = fish species
for(j in 1:C) {
mu[i,j] ~ dnorm(mu1[j], tau1[j])
}
Omega[i,1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
pred[i,1:C] ~ dmnorm(mu[i,], Omega[i,,]) # Model prediction
}
for(i in 1:C) { # C = Compound
mu1[i] ~ dnorm(0, 0.0001)
tau1[i] ~ dunif(0,10000)
pred1[i] ~ dnorm(mu1[i], tau1[i])
}
Omega1[1:C,1:C] ~ dwish(Omega0[1:C,1:C],S)
}
")
jags <- jags.model(
mod,
data = list(
S = nrow(eu3),
C = C,
Fi = Fi,
cong = log(cong),
fis = match(eu3$Fish, fisl),
Omega0 = diag(C)/100000
),
n.chains = 4,
n.adapt = 100
)
update(jags, 100)
samps.j <- jags.samples(
jags,
c(
'mu', # mean by fish and compound
'Omega', # precision matrix by fish and compound
'pred', # predicted concentration by fish and compound
# 'mu1', # mean prior for mu by compound
'Omega1', # precision matrix by compound
# 'tau1', # precision for prior of all mu
'pred1' # predicted concentration by compound
),
N
)
dimnames(samps.j$mu) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
#dimnames(samps.j$mu1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$pred) <- list(Fish = fisl, Compound = conl, Iter = 1:N, Chain = 1:4)
#dimnames(samps.j$tau1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$pred1) <- list(Compound = conl, Iter = 1:N, Chain = 1:4)
dimnames(samps.j$Omega) <- list(Fish = fisl, Compound = conl, Compound2 = conl, Iter=1:N, Chain=1:4)
dimnames(samps.j$Omega1) <- list(Compound = conl, Compound2 = conl, Iter=1:N, Chain=1:4)
##### conc.param contains expected values of the distribution parameters from the model
conc.param <- list(
mu = apply(samps.j$mu[,,,1], MARGIN = 1:2, FUN = mean),
Omega = apply(samps.j$Omega[,,,,1], MARGIN = 1:3, FUN = mean),
pred.mean = apply(samps.j$pred[,,,1], MARGIN = 1:2, FUN = mean),
pred.sd = apply(samps.j$pred[,,,1], MARGIN = 1:2, FUN = sd),
# mu1 = apply(samps.j$mu1[,,1], MARGIN = 1, FUN = mean),
# tau1 = apply(samps.j$tau1[,,1], MARGIN = 1, FUN = mean),
pred1.mean = apply(samps.j$pred1[,,1], MARGIN = 1, FUN = mean),
pred1.sd = apply(samps.j$pred1[,,1], MARGIN = 1, FUN = sd)
)
objects.store(conc.param, samps.j)
cat("Lists conc.params and samps.j stored.\n")
######################3
cat("Descriptive statistics:\n")
# Leave only the main fish species and congeners and remove others
conl <- indices$Compound.PCDDF14
eu <- eu[eu$Compound %in% conl & eu$Fish %in% fisl , ]
oprint(summary(
eu,
marginals = c("Fish", "Compound"), # Matrix is always 'Muscle'
function_names = c("mean", "sd")
))
euRatio <- EvalOutput(euRatio)
oprint(summary(
euRatio,
marginals = c("Fish", "Compound"), # Matrix is always 'Muscle'
function_names = c("mean", "sd")
))
ggplot(eu@output, aes(x = euResult, colour=Compound))+geom_density()+
facet_wrap( ~ Fish, scales = "free_y")+scale_x_log10()
#stat_ellipse()
ggplot(euRatio@output, aes(x = euRatioResult, colour = Compound))+geom_density()+
facet_wrap(~ Fish, scales = "free_y")
ggplot(melt(exp(samps.j$pred[,,,1])), aes(x=value, colour=Compound))+geom_density()+
facet_wrap( ~ Fish,scales = "free_y")+scale_x_log10()
ggplot(melt(exp(samps.j$pred1[,,1])), aes(x=value, colour=Compound))+geom_density()+
scale_x_log10()
scatterplotMatrix(t(samps.j$pred[1,,,1]), main = "Predictions for all compounds for Baltic herring")
## scatterplotMatrix(t(samps.j$mu1[,,1]), main = "Means for all compounds of the generic fish")
scatterplotMatrix(t(samps.j$pred1[,,1]), main = "Prediction for all compounds of the generic fish")
scatterplotMatrix(t(samps.j$pred[,1,,1]), main = "Predictions for all fish species for TEQdx")
scatterplotMatrix(t(samps.j$Omega[6,2,,,1]), main = "Predictions of Omega for pike and TEQpcb")
coda.j <- coda.samples(
jags,
c('mu', 'pred', 'Omega', 'pred1'),
N
)
plot(coda.j)
| |