Nothing
#' Write the JAGS model file
#'
#' \code{write_JAGS_model} creates "MixSIAR_model.txt", which is passed to JAGS
#' by \code{\link{run_model}} when the "RUN MODEL" button is clicked in the GUI.
#' Several model options will have already been specified when loading the mix
#' and source data, but here is where the error term options are selected:
#' \enumerate{
#' \item Residual * Process (resid_err = TRUE, process_err = TRUE)
#' \item Residual only (resid_err = TRUE, process_err = FALSE)
#' \item Process only (resid_err = FALSE, process_err = TRUE)
#' }
#'
#' WARNING messages are displayed if:
#' \itemize{
#' \item resid_err = FALSE and process_err = FALSE are both selected.
#' \item N=1 mix data point and did not choose "Process only" error model (MixSIR)
#' \item Fitting each individual mix data point separately as a Fixed Effect,
#' but did not choose "Process only" error model (MixSIR).
#' }
#'
#' @param filename the JAGS model file is saved in the working directory as 'filename'
#' (default is "MixSIAR_model.txt", but user can specify).
#' @param resid_err T/F: include residual error in the model?
#' @param process_err T/F: include process error in the model?
#' @param mix output from \code{\link{load_mix_data}}
#' @param source output from \code{\link{load_source_data}}
#' @export
write_JAGS_model <- function(filename = "MixSIAR_model.txt", resid_err = TRUE, process_err = TRUE, mix, source){
if(!process_err && !resid_err){
stop(paste("Invalid error structure, must choose one of:
1. Residual * Process (resid_err=TRUE, process_err=TRUE)
2. Residual only (resid_err=TRUE, process_err=FALSE)
3. Process only (resid_err=FALSE, process_err=TRUE)",sep=""))
}
if(resid_err && !process_err){
err_model <- "Residual only"
err <- "resid"
}
if(process_err && !resid_err){
err_model <- "Process only (MixSIR, for N = 1)"
err <- "process"
}
if(resid_err && process_err){
err_model <- "Residual * Process"
err <- "mult"
}
if(mix$N==1 && err!="process"){
stop(paste("Invalid error structure. If N=1 mix data point,
must choose Process only error model (MixSIR).
Set resid_err=FALSE and process_err=TRUE.",sep=""))}
if(mix$n.fe==1 && mix$N==mix$FAC[[1]]$levels && err!="process"){
stop(paste("Invalid error structure. If fitting each individual
mix data point separately, must choose Process only error model (MixSIR).
Set resid_err=FALSE and process_err=TRUE.",sep=""))}
cat(paste("# source$data_type: ",source$data_type,sep=""), file=filename)
cat("
", file=filename, append=T)
cat(paste("# source$by_factor: ",source$by_factor,sep=""), file=filename, append=T)
cat("
", file=filename, append=T)
cat(paste("# random effects: ",mix$n.re,sep=""), file=filename, append=T)
cat("
", file=filename, append=T)
cat(paste("# fixed effects: ",mix$n.fe,sep=""), file=filename, append=T)
cat("
", file=filename, append=T)
cat(paste("# nested factors: ",paste(mix$fac_nested,collapse=" "),sep=""), file=filename, append=T)
cat("
", file=filename, append=T)
cat(paste("# factors: ",paste(mix$factors,collapse=" "),sep=""), file=filename, append=T)
cat("
", file=filename, append=T)
cat(paste("# continuous effects: ",mix$n.ce,sep=""), file=filename, append=T)
cat("
", file=filename, append=T)
cat(paste("# error structure: ",err_model,sep=""), file=filename, append=T)
cat("
", file=filename, append=T)
cat(paste("# source$conc_dep: ",source$conc_dep,sep=""), file=filename, append=T)
if(source$data_type=="raw" && is.na(source$by_factor)){
cat("
var rho[n.sources,n.iso,n.iso], src_cov[n.sources,n.iso,n.iso], src_var[n.sources,n.iso,n.iso], src_Sigma[n.sources,n.iso,n.iso], Sigma.ind[N,n.iso,n.iso], mix.cov[N,n.iso,n.iso];", file=filename, append=T)
}
if(source$data_type=="raw" && !is.na(source$by_factor)){
cat("
var rho[n.sources,source_factor_levels,n.iso,n.iso], src_cov[n.sources,source_factor_levels,n.iso,n.iso], src_var[n.sources,source_factor_levels,n.iso,n.iso], src_Sigma[n.sources,source_factor_levels,n.iso,n.iso], Sigma.ind[N,n.iso,n.iso], mix.cov[N,n.iso,n.iso];", file=filename, append=T)
}
cat("
model{", file=filename, append=T)
################################################################################
# Setup src_mu and src_tau arrays
# Fit source means (if source$data_type=="raw")
################################################################################
# if(source$data_type=="raw" && source$by_factor==TRUE){ # fit the source means and precisions (by factor 1)
# cat("
# # uninformed priors on source means and precisions
# for(src in 1:n.sources){
# for(f1 in 1:source_factor_levels){
# for(iso in 1:n.iso){
# src_mu[src,iso,f1] ~ dnorm(0,.001);
# }
# }
# }
# # each source data point is distributed normally according to the source means and precisions
# for(src in 1:n.sources){
# for(f1 in 1:source_factor_levels){
# for(r in 1:n.rep[src,f1]){
# SOURCE_array[src,,f1,r] ~ dmnorm(src_mu[src,,f1],src_Sigma.inv[src,f1,,]);
# }
# }
# }
# ", file=filename, append=T)
# }
# fit the source means and precisions (by factor)
if(source$data_type=="raw" && !is.na(source$by_factor) && mix$n.iso > 1){
cat("
# fit source data (big for loop over sources)
for(src in 1:n.sources){
for(f1 in 1:source_factor_levels){ # fit source data from each source level separately
# uninformative priors on source means (src_mu vector)
for(iso in 1:n.iso){
src_mu[src,iso,f1] ~ dnorm(0,.001);
}
# uninformative priors on source variances (src_tau matrix)
for(i in 2:n.iso){
for(j in 1:(i-1)){
src_tau[src,f1,i,j] <- 0;
src_tau[src,f1,j,i] <- 0;
}
}
for(i in 1:n.iso){
src_tau[src,f1,i,i] ~ dgamma(.001,.001);
}
# uninformative priors on source correlations (rho matrix)
for(i in 2:n.iso){
for(j in 1:(i-1)){
rho[src,f1,i,j] ~ dunif(-1,1);
rho[src,f1,j,i] <- rho[src,f1,i,j];
}
}
for(i in 1:n.iso){
rho[src,f1,i,i] <- 1;
}
# Construct source precision matrix (src_Sigma)
src_var[src,f1,,] <- inverse(src_tau[src,f1,,]);
src_cov[src,f1,,] <- src_var[src,f1,,] %*% rho[src,f1,,] %*% src_var[src,f1,,];
src_Sigma[src,f1,,] <- inverse(src_cov[src,f1,,]);
# each source data point is distributed normally according to the source means and precisions
for(r in 1:n.rep[src,f1]){
SOURCE_array[src,,f1,r] ~ dmnorm(src_mu[src,,f1],src_Sigma[src,f1,,]);
}
} # end loop over f1
} # end source data fitting loop
", file=filename, append=T)
}
# fit the source means and precisions (not by factor 1)
if(source$data_type=="raw" && is.na(source$by_factor) && mix$n.iso > 1){
cat("
# fit source data (big for loop over sources)
for(src in 1:n.sources){
# uninformative priors on source means (src_mu vector)
for(iso in 1:n.iso){
src_mu[src,iso] ~ dnorm(0,.001);
}
# uninformative priors on source variances (src_tau matrix)
for(i in 2:n.iso){
for(j in 1:(i-1)){
src_tau[src,i,j] <- 0;
src_tau[src,j,i] <- 0;
}
}
for(i in 1:n.iso){
src_tau[src,i,i] ~ dgamma(.001,.001);
}
# uninformative priors on source correlations (rho matrix)
for(i in 2:n.iso){
for(j in 1:(i-1)){
rho[src,i,j] ~ dunif(-1,1);
rho[src,j,i] <- rho[src,i,j];
}
}
for(i in 1:n.iso){
rho[src,i,i] <- 1;
}
# Construct source precision matrix (src_Sigma)
src_var[src,,] <- inverse(src_tau[src,,]);
src_cov[src,,] <- src_var[src,,] %*% rho[src,,] %*% src_var[src,,];
src_Sigma[src,,] <- inverse(src_cov[src,,]);
# each source data point is distributed normally according to the source means and precisions
for(r in 1:n.rep[src]){
SOURCE_array[src,,r] ~ dmnorm(src_mu[src,],src_Sigma[src,,]);
}
} # end source data fitting loop
", file=filename, append=T)
}
if(source$data_type=="raw" && mix$n.iso == 1){ # if one iso, can't call "i in 2:n.iso"
cat("
# fit source data (big for loop over sources)
for(src in 1:n.sources){
# uninformative priors on source means and precisions, rho=1
for(iso in 1:n.iso){
src_mu[src,iso] ~ dnorm(0,.001);
src_tau[src,iso,iso] ~ dgamma(.001,.001);
rho[src,iso,iso] <- 1;
}
# Construct source precision matrix (src_Sigma)
src_var[src,,] <- 1/src_tau[src,,];
src_cov[src,,] <- src_var[src,,] %*% rho[src,,] %*% src_var[src,,];
src_Sigma[src,,] <- 1/src_cov[src,,];
# each source data point is distributed normally according to the source means and precisions
for(r in 1:n.rep[src]){
SOURCE_array[src,,r] ~ dnorm(src_mu[src,],src_Sigma[src,,]);
}
} # end source data fitting loop
", file=filename, append=T)
}
# Here we fit the source means and variances according to Gelman, p.79-80:
# u | sig2,y ~ Normal(m, sig2/n), (Eqn 3.8)
# sig2 | y ~ Inv-X2(n, n-1/n * s2), (Eqn 3.9)
# where u and sig2 are the true source mean and variance,
# m and s2 are the sample mean and variance, and n is the sample size
# We are using an uninformative prior, so k_0 = v_0 = 0
# To simulate sig2, draw tmp.X from chisqr(n), then src_tau = tmp.X/(n-1)*s2 (Gelman p.580)
# If we have only residual error (SIAR model), then we only care about src_mu (src_tau and frac_sig2 go unused)
# If we have process error (MixSIR model), then we use src_mu and src_tau (and frac_sig2, which isn't fit to data yet)
# src_mu and src_tau are used in the X[i,iso] ~ dnorm call at the very bottom of this function
# dim(MU/SIG2_array) is either (n.src,n.iso) or (n.src,n.iso,n.fac), if source$by_factor is FALSE or TRUE, respectively.
if(source$data_type=="means" && is.na(source$by_factor)){
cat("
for(src in 1:n.sources){
for(iso in 1:n.iso){
src_mu[src,iso] ~ dnorm(MU_array[src,iso], n_array[src]/SIG2_array[src,iso]); # Eqn 3.8 but with precision instead of variance
tmp.X[src,iso] ~ dchisqr(n_array[src]);
src_tau[src,iso] <- tmp.X[src,iso]/(SIG2_array[src,iso]*(n_array[src] - 1)); # Eqn 3.9, following the simulation on p.580
}
}
", file=filename, append=T)
}
if(source$data_type=="means" && !is.na(source$by_factor)){
cat("
for(src in 1:n.sources){
for(f1 in 1:source_factor_levels){
for(iso in 1:n.iso){
src_mu[src,iso,f1] ~ dnorm(MU_array[src,iso,f1], n_array[src,f1]/SIG2_array[src,iso,f1]); # Eqn 3.8 but with precision instead of variance
tmp.X[src,iso,f1] ~ dchisqr(n_array[src,f1]);
src_tau[src,iso,f1] <- tmp.X[src,iso,f1]/(SIG2_array[src,iso,f1]*(n_array[src,f1] - 1)); # Eqn 3.9, following the simulation on p.580
}
}
}
", file=filename, append=T)
}
###########################################################################
# ILR transform and factors
###########################################################################
cat("
# draw p.global (global proportion means) from an uninformative Dirichlet,
# then ilr.global is the ILR-transform of p.global
p.global[1:n.sources] ~ ddirch(alpha[1:n.sources]);
for(src in 1:(n.sources-1)){
gmean[src] <- prod(p.global[1:src])^(1/src);
ilr.global[src] <- sqrt(src/(src+1))*log(gmean[src]/p.global[src+1]); # page 296, Egozcue 2003
}
", file=filename, append=T)
if(mix$n.effects > 0 && mix$FAC[[1]]$re){ # factor 1 is random effect
cat("
fac1.sig ~ dunif(0,20);
fac1.invSig2 <- 1/(fac1.sig*fac1.sig);
# draw the fac1 (region) specific ILR terms (random effect)
for(f1 in 1:factor1_levels) {
for(src in 1:(n.sources-1)) {
ilr.fac1[f1,src] ~ dnorm(0,fac1.invSig2);
}
}
", file=filename, append=T)}
if(mix$n.effects > 0 && !mix$FAC[[1]]$re){ # factor 1 is fixed effect
cat("
# draw the fac1 specific ILR terms (fixed effect)
for(src in 1:(n.sources-1)){
ilr.fac1[1,src] <- 0;
for(f1 in 2:factor1_levels){
ilr.fac1[f1,src] ~ dnorm(0,1);
}
}
", file=filename, append=T)}
if(mix$n.effects > 1 && mix$FAC[[2]]$re){ # factor 2 is random effect
cat("
fac2.sig ~ dunif(0,20);
fac2.invSig2 <- 1/(fac2.sig*fac2.sig);
# draw the fac2-specific ILR terms (random effect)
for(f2 in 1:factor2_levels){
for(src in 1:(n.sources-1)){
ilr.fac2[f2,src] ~ dnorm(0,fac2.invSig2);
}
}
", file=filename, append=T)}
if(mix$n.effects > 1 && !mix$FAC[[2]]$re){ # factor 2 is fixed effect
cat("
# draw the fac2 specific ILR terms (fixed effect)
for(src in 1:(n.sources-1)){
ilr.fac2[1,src] <- 0;
for(f2 in 2:factor2_levels){
ilr.fac2[f2,src] ~ dnorm(0,1);
}
}
", file=filename, append=T)}
# Continuous Effects section
ilr.cont.string <- ""
if(mix$n.ce > 0){
cat("
# Priors on all n.ce continuous effects (slopes for a linear regression in ilr-space)", file=filename, append=T)
for(ce in 1:mix$n.ce){
cat("
for(src in 1:(n.sources-1)){
ilr.cont",ce,"[src] ~ dnorm(0,.001)
}
", file=filename, append=T, sep="")
ilr.cont.string <- paste(ilr.cont.string," + ilr.cont",ce,"[src]*Cont.",ce,"[i]",sep="")}
}
# if(indiv_effect){ # if user has chosen to include Individual as a random effect
# cat("
# ind.sig ~ dunif(0,20);
# ind.invSig2 <- 1/(ind.sig*ind.sig);
# # generate individual deviates from the global/region/pack mean
# for(i in 1:N) {
# for(src in 1:(n.sources-1)) {
# ilr.ind[i,src] ~ dnorm(0, ind.invSig2);", file=filename, append=T)
# } else { # user has chosen to NOT include Individual as a random effect - remove ind.sig and set ilr.ind=0
cat("
# DON'T generate individual deviates from the global/region/pack mean (but keep same model structure)
for(i in 1:N) {
for(src in 1:(n.sources-1)) {
ilr.ind[i,src] <- 0;", file=filename, append=T)
# } # end Individual effect block
# Add ilr.tot line, adding ilr.global, ilr.fac's, ilr.ind, and ilr.cont's
# This will be different according to n.effects, but *MUST* be directly after the individual effect section!!
# Can condense in the future by creating a "ilr.rand.string" object in a loop so this line is the same for all n.re
if(mix$n.effects==2){
cat("
ilr.tot[i,src] <- ilr.global[src] + ilr.fac1[Factor.1[i],src] + ilr.fac2[Factor.2[i],src]",ilr.cont.string," + ilr.ind[i,src]; # add all effects together for each individual (in ilr-space)
}
}
", file=filename, append=T, sep="")}
if(mix$n.effects==1){
cat("
ilr.tot[i,src] <- ilr.global[src] + ilr.fac1[Factor.1[i],src]",ilr.cont.string," + ilr.ind[i,src]; # add all effects together for each individual (in ilr-space)
}
}
", file=filename, append=T, sep="")}
if(mix$n.effects==0){
cat("
ilr.tot[i,src] <- ilr.global[src]",ilr.cont.string," + ilr.ind[i,src]; # add all effects together for each individual (in ilr-space)
}
}
", file=filename, append=T, sep="")}
#####################################################################
# Inverse ILR section
#####################################################################
cat("
# Inverse ILR math (equation 24, page 294, Egozcue 2003)
for(i in 1:N){
for(j in 1:(n.sources-1)){
cross[i,,j] <- (e[,j]^ilr.tot[i,j])/sum(e[,j]^ilr.tot[i,j]);
}
for(src in 1:n.sources){
tmp.p[i,src] <- prod(cross[i,src,]);
}
for(src in 1:n.sources){
p.ind[i,src] <- tmp.p[i,src]/sum(tmp.p[i,]);
}
}
for(src in 1:n.sources) {
for(i in 1:N){
# these are weights for variances
p2[i,src] <- p.ind[i,src]*p.ind[i,src];
}
}
", file=filename, append=T)
if(mix$n.effects > 0 & !mix$fere){
if(mix$fac_nested[1]){
cat("
# Transform ilr.fac1 into p.fac1 (fac1 nested within fac2)
for(f1 in 1:factor1_levels) {
for(src in 1:(n.sources-1)) {
ilr.fac1.tot[f1,src] <- ilr.global[src] + ilr.fac2[factor2_lookup[f1],src] + ilr.fac1[f1,src];
cross.fac1[f1,,src] <- (e[,src]^ilr.fac1.tot[f1,src])/sum(e[,src]^ilr.fac1.tot[f1,src]);
}
for(src in 1:n.sources) {
tmp.p.fac1[f1,src] <- prod(cross.fac1[f1,src,]);
}
for(src in 1:n.sources){
p.fac1[f1,src] <- tmp.p.fac1[f1,src]/sum(tmp.p.fac1[f1,]);
}
}
", file=filename, append=T)
} else {
cat("
# Transform ilr.fac1 into p.fac1 (fac1 not nested within fac2)
for(f1 in 1:factor1_levels) {
for(src in 1:(n.sources-1)) {
ilr.fac1.tot[f1,src] <- ilr.global[src] + ilr.fac1[f1,src];
cross.fac1[f1,,src] <- (e[,src]^ilr.fac1.tot[f1,src])/sum(e[,src]^ilr.fac1.tot[f1,src]);
}
for(src in 1:n.sources) {
tmp.p.fac1[f1,src] <- prod(cross.fac1[f1,src,]);
}
for(src in 1:n.sources){
p.fac1[f1,src] <- tmp.p.fac1[f1,src]/sum(tmp.p.fac1[f1,]);
}
}
", file=filename, append=T)
} # end nested ilr.fac1 section
} # end fac1 section
if(mix$n.effects > 1 & !mix$fere){ # i.e. n.re=2
if(mix$fac_nested[2]){
cat("
# Transform ilr.fac2 into p.fac2 (fac2 nested within fac1)
for(f2 in 1:factor2_levels){
for(src in 1:(n.sources-1)){
ilr.fac2.tot[f2,src] <- ilr.global[src] + ilr.fac1[factor1_lookup[f2],src] + ilr.fac2[f2,src];
cross.fac2[f2,,src] <- (e[,src]^ilr.fac2.tot[f2,src])/sum(e[,src]^ilr.fac2.tot[f2,src]);
}
for(src in 1:n.sources) {
tmp.p.fac2[f2,src] <- prod(cross.fac2[f2,src,]);
}
for(src in 1:n.sources){
p.fac2[f2,src] <- tmp.p.fac2[f2,src]/sum(tmp.p.fac2[f2,]);
}
}
", file=filename, append=T)
} else {
cat("
# Transform ilr.fac2 into p.fac2 (fac2 not nested within fac1)
for(f2 in 1:factor2_levels){
for(src in 1:(n.sources-1)){
ilr.fac2.tot[f2,src] <- ilr.global[src] + ilr.fac2[f2,src];
cross.fac2[f2,,src] <- (e[,src]^ilr.fac2.tot[f2,src])/sum(e[,src]^ilr.fac2.tot[f2,src]);
}
for(src in 1:n.sources) {
tmp.p.fac2[f2,src] <- prod(cross.fac2[f2,src,]);
}
for(src in 1:n.sources){
p.fac2[f2,src] <- tmp.p.fac2[f2,src]/sum(tmp.p.fac2[f2,]);
}
}
", file=filename, append=T)
} # end nested ilr.fac2 section
} # end n.re=2 section
# 2 FE or 1 FE + 1 RE section
if(mix$fere){
if(mix$n.re==1){ # if 1 FE and 1 RE, get p.fac1 (fixed)
cat("
# Transform ilr.fac1 into p.fac1 (fac1 fixed, not nested within fac2)
for(f1 in 1:factor1_levels) {
for(src in 1:(n.sources-1)) {
ilr.fac1.tot[f1,src] <- ilr.global[src] + ilr.fac1[f1,src];
cross.fac1[f1,,src] <- (e[,src]^ilr.fac1.tot[f1,src])/sum(e[,src]^ilr.fac1.tot[f1,src]);
}
for(src in 1:n.sources) {
tmp.p.fac1[f1,src] <- prod(cross.fac1[f1,src,]);
}
for(src in 1:n.sources){
p.fac1[f1,src] <- tmp.p.fac1[f1,src]/sum(tmp.p.fac1[f1,]);
}
}
", file=filename, append=T)
}
# for both 2 FE and 1FE + 1RE, don't get p.fac2. Instead, get p.both[f1,f2]
# cat("
# # Transform ilr.fac2 into p.both (fac1 fixed, so don't get p.fac2)
# for(f1 in 1:factor1_levels) {
# for(f2 in factor2_lookup[[f1]]){
# for(src in 1:(n.sources-1)) {
# ilr.both[f1,f2,src] <- ilr.global[src] + ilr.fac1[f1,src] + ilr.fac2[f2,src];
# cross.both[f1,f2,,src] <- (e[,src]^ilr.both[f1,f2,src])/sum(e[,src]^ilr.both[f1,f2,src]);
# }
# for(src in 1:n.sources) {
# tmp.p.both[f1,f2,src] <- prod(cross.both[f1,f2,src,]);
# }
# for(src in 1:n.sources){
# p.both[f1,f2,src] <- tmp.p.both[f1,f2,src]/sum(tmp.p.both[f1,f2,]);
# }
# }
# }
# ", file=filename, append=T)
}
###############################################################################
# mix.mu section
###############################################################################
cat("
# for each isotope and population, calculate the predicted mixtures
for(iso in 1:n.iso) {
for(i in 1:N) {
", file=filename, append=T)
if(!is.na(source$by_factor) && source$conc_dep==T){
if(source$by_factor == 1){
cat("
mix.mu[iso,i] <- (inprod(src_mu[,iso,Factor.1[i]],(p.ind[i,]*conc[,iso])) + inprod(frac_mu[,iso],(p.ind[i,]*conc[,iso]))) / inprod(p.ind[i,],conc[,iso]);", file=filename, append=T)
} else { # by_factor == 2
cat("
mix.mu[iso,i] <- (inprod(src_mu[,iso,Factor.2[i]],(p.ind[i,]*conc[,iso])) + inprod(frac_mu[,iso],(p.ind[i,]*conc[,iso]))) / inprod(p.ind[i,],conc[,iso]);", file=filename, append=T)
}
} else if(!is.na(source$by_factor) && source$conc_dep==F){
if(source$by_factor == 1){
cat("
mix.mu[iso,i] <- inprod(src_mu[,iso,Factor.1[i]],p.ind[i,]) + inprod(frac_mu[,iso],p.ind[i,]);", file=filename, append=T)
} else { # by_factor == 2
cat("
mix.mu[iso,i] <- inprod(src_mu[,iso,Factor.2[i]],p.ind[i,]) + inprod(frac_mu[,iso],p.ind[i,]);", file=filename, append=T)
}
} else if(is.na(source$by_factor) && source$conc_dep==T){
cat("
mix.mu[iso,i] <- (inprod(src_mu[,iso],(p.ind[i,]*conc[,iso])) + inprod(frac_mu[,iso],(p.ind[i,]*conc[,iso]))) / inprod(p.ind[i,],conc[,iso]);", file=filename, append=T)
} else if(is.na(source$by_factor) && source$conc_dep==F){
cat("
mix.mu[iso,i] <- inprod(src_mu[,iso],p.ind[i,]) + inprod(frac_mu[,iso],p.ind[i,]);", file=filename, append=T)
}
cat("
}
}
", file=filename, append=T)
################################################################################
# mix.prcsn section (error terms)
################################################################################
# # Add process/MixSIR error (define mix.var or set equal to 0)
# if(process_err){ # if process_err = TRUE, there are two varieties of mix.var (source$by_factor = T/F)
# if(source$by_factor){ # source$by_factor = TRUE, include process error as:
# cat("
# process.var[iso,i] <- inprod(1/src_tau[,iso,Factor.1[i]],p2[i,]) + inprod(frac_sig2[,iso],p2[i,]);", file=filename, append=T)
# } else { # source$by_factor = FALSE, include process error as:
# cat("
# process.var[iso,i] <- inprod(1/src_tau[,iso],p2[i,]) + inprod(frac_sig2[,iso],p2[i,]);", file=filename, append=T)
# } # end if/else(source$by_factor)
# } else { # process_err = FALSE, set process.var[iso,i] = 0:
# cat("
# process.var[iso,i] <- 0;", file=filename, append=T)
# } # end if/else(process_err)
# # Calculate total variance (mix.var)
# if(resid_err){ # resid_err==TRUE
# cat("
# mix.var[iso,i] <- process.var[iso,i] + resid.var[iso];
# mix.prcsn[iso,i] <- 1/mix.var[iso,i];
# }
# }
# ", file=filename, append=T)
# } else { # resid_err==FALSE
# if(resid_err_mult){
# cat("
# mix.var[iso,i] <- process.var[iso,i] * resid.prop[iso];
# mix.prcsn[iso,i] <- 1/mix.var[iso,i];
# }
# }
# ", file=filename, append=T)
# } else { # resid_err==FALSE AND resid_err_mult==FALSE (process_err only = MixSIR)
# cat("
# mix.var[iso,i] <- process.var[iso,i];
# mix.prcsn[iso,i] <- 1/mix.var[iso,i];
# }
# }
# ", file=filename, append=T)
# }
# } # end resid_err==FALSE
###############################################################################
# Error structure section
###############################################################################
if(err=="mult"){
cat("
# Multiplicative residual error
for(iso in 1:n.iso){
resid.prop[iso] ~ dunif(0,20);
}
", file=filename, append=T)
if(source$data_type=="means"){
cat("
# Calculate process variance for each isotope and population
for(iso in 1:n.iso) {
for(i in 1:N) {
", file=filename, append=T)
if(!is.na(source$by_factor)){ # source$by_factor = TRUE, include process error as:
if(source$by_factor == 1){
cat("
process.var[iso,i] <- inprod(1/src_tau[,iso,Factor.1[i]],p2[i,]) + inprod(frac_sig2[,iso],p2[i,]);", file=filename, append=T)
} else { # by factor 2
cat("
process.var[iso,i] <- inprod(1/src_tau[,iso,Factor.2[i]],p2[i,]) + inprod(frac_sig2[,iso],p2[i,]);", file=filename, append=T)
}
} else { # source$by_factor = NA, include process error as:
cat("
process.var[iso,i] <- inprod(1/src_tau[,iso],p2[i,]) + inprod(frac_sig2[,iso],p2[i,]);", file=filename, append=T)
}
cat("
}
}
# Construct Sigma, the mixture precision matrix
for(ind in 1:N){
for(i in 1:n.iso){
for(j in 1:n.iso){
Sigma.ind[ind,i,j] <- equals(i,j)/(process.var[i,ind]*resid.prop[i]);
}
}
}
# Likelihood
for(i in 1:N) {
", file=filename, append=T)
if(mix$n.iso > 1){
cat("
X_iso[i,] ~ dmnorm(mix.mu[,i], Sigma.ind[i,,]);
loglik[i] <- logdensity.mnorm(X_iso[i,], mix.mu[,i], Sigma.ind[i,,]);", file=filename, append=T)
} else { # n.iso == 1, can't use dmnorm
cat("
X_iso[i,] ~ dnorm(mix.mu[,i], Sigma.ind[i,,]);
loglik[i] <- logdensity.norm(X_iso[i,], mix.mu[,i], Sigma.ind[i,,]);", file=filename, append=T)
}
cat("
}
} # end model
", file=filename, append=T)
} # end source$data_type=="means"
if(source$data_type=="raw"){
cat("
# resid.prop matrix
for(i in 1:n.iso){
for(j in 1:n.iso){
resid.prop.mat[i,j] <- sqrt(resid.prop[i]*resid.prop[j]);
}
}
# Construct mix covariance
for(ind in 1:N){
for(i in 1:n.iso){
for(j in 1:n.iso){
", file=filename, append=T)
if(!is.na(source$by_factor)){ # source$by_factor = 1 or 2, include process error as:
if(source$by_factor == 1){
cat("
mix.cov[ind,i,j] <- equals(i,j)*resid.prop[i]*(inprod(src_cov[,Factor.1[ind],i,j],p2[ind,]) + inprod(frac_sig2[,i],p2[ind,])) + (1-equals(i,j))*inprod(src_cov[,Factor.1[ind],i,j],p2[ind,])*resid.prop.mat[i,j];", file=filename, append=T)
} else { # by factor = 2
cat("
mix.cov[ind,i,j] <- equals(i,j)*resid.prop[i]*(inprod(src_cov[,Factor.2[ind],i,j],p2[ind,]) + inprod(frac_sig2[,i],p2[ind,])) + (1-equals(i,j))*inprod(src_cov[,Factor.2[ind],i,j],p2[ind,])*resid.prop.mat[i,j];", file=filename, append=T)
}
} else { # source$by_factor = NA, include process error as:
cat("
mix.cov[ind,i,j] <- equals(i,j)*resid.prop[i]*(inprod(src_cov[,i,j],p2[ind,]) + inprod(frac_sig2[,i],p2[ind,])) + (1-equals(i,j))*inprod(src_cov[,i,j],p2[ind,])*resid.prop.mat[i,j];", file=filename, append=T)
}
cat("
}
}
Sigma.ind[ind,,] <- inverse(mix.cov[ind,,]);
}
# Likelihood
for(i in 1:N){
", file=filename, append=T)
if(mix$n.iso > 1){
cat("
X_iso[i,] ~ dmnorm(mix.mu[,i], Sigma.ind[i,,]);
loglik[i] <- logdensity.mnorm(X_iso[i,], mix.mu[,i], Sigma.ind[i,,]);", file=filename, append=T)
} else { # n.iso == 1, can't use dmnorm
cat("
X_iso[i,] ~ dnorm(mix.mu[,i], Sigma.ind[i,,]);
loglik[i] <- logdensity.norm(X_iso[i,], mix.mu[,i], Sigma.ind[i,,]);", file=filename, append=T)
}
cat("
}
} # end model
", file=filename, append=T)
} # end source$data_type=="raw"
} # end multiplicative residual error section
# Residual error only section
if(err=="resid" && mix$n.iso>1){ # > 1 iso, have covariance
cat("
# Mixture covariance prior (residual error only model)
Sigma ~ dwish(I,n.iso+1);
# Likelihood
for(i in 1:N) {
X_iso[i,] ~ dmnorm(mix.mu[,i], Sigma);
loglik[i] <- logdensity.mnorm(X_iso[i,], mix.mu[,i], Sigma);
}
} # end model
", file=filename, append=T)
}
if(err=="resid" && mix$n.iso==1){ # if only one iso can't use dwish or dmnorm
cat("
# Mixture precision prior (residual error only model)
Sigma ~ dgamma(.001,.001);
# Likelihood
for(i in 1:N) {
X_iso[i,] ~ dnorm(mix.mu[,i], Sigma);
loglik[i] <- logdensity.norm(X_iso[i,], mix.mu[,i], Sigma);
}
} # end model
", file=filename, append=T)
}
# MixSIR/process error section (for N = 1)
if(err=="process"){
cat("
# calculate mix variance and likelihood
for(i in 1:N){
for(iso in 1:n.iso){
", file=filename, append=T)
if(source$data_type=="raw"){ # source data = raw, src_tau dim = src,iso,iso
cat("
process.var[iso,i] <- inprod(1/src_tau[,iso,iso],p2[i,]) + inprod(frac_sig2[,iso],p2[i,]);", file=filename, append=T)
} else { # source data = means+SD, src_tau dim = src, iso
cat("
process.var[iso,i] <- inprod(1/src_tau[,iso],p2[i,]) + inprod(frac_sig2[,iso],p2[i,]);", file=filename, append=T)
}
cat("
mix.prcsn[iso,i] <- 1/process.var[iso,i];
X_iso[i,iso] ~ dnorm(mix.mu[iso,i], mix.prcsn[iso,i]);
loglik_mat[i,iso] <- logdensity.norm(X_iso[i,iso], mix.mu[iso,i], mix.prcsn[iso,i]);
}
loglik[i] <- sum(loglik_mat[i,])
}
} # end model
", file=filename, append=T)
} # end MixSIR error section (N = 1)
} # end function write_JAGS_model
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.