##### these are internal functions.
#' Write .bug model file for model without regression
#'
#' `write_model_NoReg` automatically generates model file according to
#' `model_options`
#'
#' @inheritParams insert_bugfile_chunk_noreg_meas
#'
#' @return a long character string to be written into .bug file.
#'
#' @seealso
#' [insert_bugfile_chunk_noreg_meas] for inserting .bug file
#' chunk for measurements (plug-and-play); [insert_bugfile_chunk_noreg_etiology]
#' for inserting .bug file chunk for distribution of latent status (etiology).
#'
#' @family model generation functions
write_model_NoReg <- function(k_subclass,Mobs,prior,cause_list,
use_measurements,ppd=NULL,use_jags=FALSE){
## 1) accommodate singletons, combos, and NoA;
## 2) If use_jags=FALSE, check-bit to prevent case data informing FPR (see the use of cut() functions in "plug-and-play.R");
chunk1 <- insert_bugfile_chunk_noreg_meas(k_subclass,Mobs,
prior,cause_list,use_measurements,ppd,use_jags)
chunk2 <- insert_bugfile_chunk_noreg_etiology(ppd)
paste0("model{#BEGIN OF MODEL:\n",
chunk1,"\n",
chunk2,"\n",
"}#END OF MODEL.")
}
#' Write .bug model file for regression model without nested subclasses
#'
#' `write_model_Reg_NoNest` automatically generates model file according to
#' `model_options`
#'
#' @inheritParams insert_bugfile_chunk_reg_nonest_meas
#' @param Eti_formula Etiology regression formula; Check `model_options$likelihood$Eti_formula`.
#' @param FPR_formula A list of FPR regression formula; check
#' `model_options$likelihood$FPR_formula`
#' @return a long character string to be written into .bug file.
#'
#' @seealso
#' [insert_bugfile_chunk_noreg_meas] for inserting .bug file
#' chunk for measurements (plug-and-play); [insert_bugfile_chunk_noreg_etiology]
#' for inserting .bug file chunk for distribution of latent status (etiology).
#'
#' @family model generation functions
write_model_Reg_NoNest <- function(Mobs,prior,cause_list,Eti_formula,FPR_formula,
use_measurements,ppd=NULL,use_jags=FALSE){
chunk1 <- insert_bugfile_chunk_reg_nonest_meas(Mobs,
prior,cause_list,FPR_formula,
use_measurements,ppd,use_jags)
chunk2 <- insert_bugfile_chunk_reg_etiology(Eti_formula,length(cause_list),ppd) #DONE.
paste0("model{#BEGIN OF MODEL:\n",
chunk1,"\n",
chunk2,"\n",
"}#END OF MODEL.")
}
#' Write .bug model file for regression model without nested subclasses
#'
#' `write_model_Reg_discrete_predictor_NoNest` automatically generates model file according to
#' `model_options`
#'
#' @inheritParams insert_bugfile_chunk_reg_nonest_meas
#'
#' @return a long character string to be written into .bug file.
#'
#' @seealso
#' [insert_bugfile_chunk_noreg_meas] for inserting .bug file
#' chunk for measurements (plug-and-play); [insert_bugfile_chunk_noreg_etiology]
#' for inserting .bug file chunk for distribution of latent status (etiology).
#'
#' @family model generation functions
write_model_Reg_discrete_predictor_NoNest <- function(Mobs,prior,cause_list,
use_measurements,ppd=NULL,use_jags=FALSE){
chunk1 <- insert_bugfile_chunk_reg_discrete_predictor_nonest_meas(Mobs,
prior,cause_list,
use_measurements,ppd,use_jags)
chunk2 <- insert_bugfile_chunk_reg_discrete_predictor_etiology(length(cause_list),ppd) #DONE.
paste0("model{#BEGIN OF MODEL:\n",
chunk1,"\n",
chunk2,"\n",
"}#END OF MODEL.")
}
#' Write `.bug` model file for regression model WITH nested subclasses
#'
#' `write_model_Reg_Nest` automatically generates model file according to
#' `model_options`; This is called within [nplcm_fit_Reg_Nest].
#'
#' @inheritParams insert_bugfile_chunk_reg_nonest_meas
#' @param Eti_formula Etiology regression formula; Check `model_options$likelihood$Eti_formula`.
#' @param FPR_formula A list of FPR regression formula; check
#' `model_options$likelihood$FPR_formula`
#' @return a long character string to be written into .bug file.
#'
#' @seealso
#' [insert_bugfile_chunk_noreg_meas] for inserting .bug file
#' chunk for measurements (plug-and-play.R); [insert_bugfile_chunk_noreg_etiology]
#' for inserting .bug file chunk for distribution of latent status (etiology).
#'
#' @family model generation functions
write_model_Reg_Nest <- function(Mobs,prior,cause_list,Eti_formula,FPR_formula,
use_measurements,ppd=NULL,use_jags=FALSE){
chunk1 <- insert_bugfile_chunk_reg_nest_meas(Mobs,
prior,cause_list,FPR_formula,
use_measurements,ppd,use_jags)
chunk2 <- insert_bugfile_chunk_reg_etiology(Eti_formula,length(cause_list),ppd)
paste0("model{#BEGIN OF MODEL:\n",
chunk1,"\n",
chunk2,"\n",
"}#END OF MODEL.")
}
##############################################################################
############### insert-bugfile-chunk.R #################
##############################################################################
#' Insert measurement likelihood (without regression) code chunks into .bug model file
#'
#' @param k_subclass the number of subclasses for the slices that require
#' conditional dependence modeling (only applicable to BrS data); its length is
#' of the same value as the number of BrS slices.
#' @param Mobs measurement data in the form of `data_nplcm`
#' @param prior prior specification from `model_options`
#' @param cause_list a list of latent status names (crucial for building templates;
#' see [make_template()])
#' @param use_measurements "BrS", or "SS"
#' @param ppd Default is NULL; set to TRUE for posterior predictive checking
#' @param use_jags Default is FALSE; set to TRUE if want to use JAGS for model fitting.
#'
#' @return a long character string to be inserted into .bug model file as measurement
#' likelihood
#'
#' @seealso It is used in [write_model_NoReg] for constructing a .bug file along with
#' specification of latent status distribution ([insert_bugfile_chunk_noreg_etiology])
insert_bugfile_chunk_noreg_meas <-
function(k_subclass,Mobs,prior,cause_list,use_measurements = "BrS",ppd=NULL,use_jags=FALSE) {
if (!("BrS" %in% use_measurements) && !("SS" %in% use_measurements)){
stop("==[baker] No BrS or SS measurements specified in the model! ==")
}
for (s in seq_along(Mobs$MBS)){
if (k_subclass[s]>1 && ncol(Mobs$MBS[[s]])==1){
stop("==[baker] Cannot do nested modeling for BrS measurements with only one column! ==")
}
}
# generate file:
if ("BrS" %in% use_measurements){
chunk_BrS_case <- ""
chunk_BrS_ctrl <- ""
chunk_BrS_subclass <- ""
chunk_BrS_param <- ""
for (s in seq_along(Mobs$MBS)){# begin iterate over slices:
k_curr <- k_subclass[s]
if (k_curr == 1 ){
if (!use_jags){ # use winbugs:
chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_NoNest_Slice(s,Mobs,cause_list,ppd)$plug)
chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_NoNest_Slice(s,Mobs,cause_list)$plug)
} else{# use jags:
chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_NoNest_Slice_jags(s,Mobs,prior,cause_list,ppd)$plug)
chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_NoNest_Slice_jags(s,Mobs,prior,cause_list)$plug)
}
chunk_BrS_ctrl <- paste0(chunk_BrS_ctrl, add_meas_BrS_ctrl_NoNest_Slice(s,Mobs,cause_list,ppd)$plug)
}
if ((k_curr > 1 )){
if (!use_jags){# use winbugs:
chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_Nest_Slice(s,Mobs,cause_list)$plug)
}else{#use jags:
chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_Nest_Slice_jags(s,Mobs,cause_list,ppd)$plug)
chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_Nest_Slice_jags(s,Mobs,cause_list)$plug)
}
chunk_BrS_subclass <- paste0(chunk_BrS_subclass, add_meas_BrS_subclass_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
chunk_BrS_ctrl <- paste0(chunk_BrS_ctrl, add_meas_BrS_ctrl_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
}
}# end iterate over slices.
}
if ("BrS" %in% use_measurements & !("SS" %in% use_measurements)) {
chunk <- paste0(
"
# BrS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# bronze-standard measurement characteristics:
",chunk_BrS_param
)
}
if ("BrS" %in% use_measurements & ("SS" %in% use_measurements)) {
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# BrS and SS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# BrS and SS measurement characteristics:
",chunk_BrS_param,
add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
if (!("BrS" %in% use_measurements) & ("SS" %in% use_measurements)) {
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# SS measurements:
for (i in 1:Nd){
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
# SS measurement characteristics:
",add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
paste0(chunk,"\n")
}
#' insert distribution for latent status code chunk into .bug file
#'
#' @param ppd Default is NULL; set to TRUE for posterior predictive checking
#'
#' @return a long character string to be inserted into .bug model file
#' as distribution specification for latent status
insert_bugfile_chunk_noreg_etiology <- function(ppd = NULL){
ppd_seg <- ""
if (!is.null(ppd) && ppd){
ppd_seg <- paste0("Icat.new[i] ~ dcat(pEti[1:Jcause])")
}
chunk_etiology <- paste0(
"
# etiology priors
for (i in 1:Nd){
Icat[i] ~ dcat(pEti[1:Jcause])
",
ppd_seg,"
}
pEti[1:Jcause]~ddirch(alphaEti[])"
)
paste0(chunk_etiology,"\n")
}
##############################################################################
############### REGRESSION MODELS #################
##############################################################################
#' Insert measurement likelihood (with regression) code chunks into .bug model file
#'
#' @param Mobs Measurement data in the form of `data_nplcm`
#' @param prior Prior specification from `model_options`
#' @param cause_list A list of latent status names (crucial for building templates;
#' see [make_template()])
#' @param FPR_formula A list of FPR regression formula; check `model_options$likelihood$FPR_formula`
#' @param use_measurements "BrS", or "SS"
#' @param ppd Default is NULL; set to TRUE for posterior predictive checking
#' @param use_jags Default is FALSE; set to TRUE if want to use JAGS for model fitting.
#'
#' @return A long character string to be inserted into .bug model file as measurement
#' likelihood
#'
#' @seealso It is used in [write_model_Reg_NoNest] for constructing a .bug file along with
#' specification of latent status regression ([insert_bugfile_chunk_reg_etiology])
#'
insert_bugfile_chunk_reg_nonest_meas <-
function(Mobs,prior,cause_list,FPR_formula,use_measurements = "BrS",ppd=NULL,use_jags=FALSE) {
if (!("BrS" %in% use_measurements) && !("SS" %in% use_measurements)){
stop("==[baker] No BrS or SS measurements specified in the model! ==")
}
# for (s in seq_along(Mobs$MBS)){
# if (k_subclass[s]>1 && ncol(Mobs$MBS[[s]])==1){
# stop("==[baker] Cannot do nested modeling for BrS measurements with only one column! ==")
# }
# }
# generate file:
if ("BrS" %in% use_measurements){
chunk_BrS_case <- ""
chunk_BrS_ctrl <- ""
chunk_BrS_subclass <- ""
chunk_BrS_param <- ""
for (s in seq_along(Mobs$MBS)){# begin iterate over slices:
# k_curr <- k_subclass[s]
# if (k_curr == 1 ){
if (!use_jags){ # use winbugs:
stop("==[baker] WinBUGS code not developed. Set 'use_jags=TRUE'. ==\n")
#chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_NoNest_reg_Slice_jags(s,Mobs,cause_list,ppd)$plug)
#chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_NoNest_reg_Slice_jags(s,Mobs,cause_list)$plug)
} else{# use jags:
chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_NoNest_reg_Slice_jags(s,Mobs,prior,cause_list,ppd)$plug)
chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_NoNest_reg_Slice_jags(s,Mobs,prior,cause_list,FPR_formula)$plug)
}
chunk_BrS_ctrl <- paste0(chunk_BrS_ctrl, add_meas_BrS_ctrl_NoNest_reg_Slice_jags(s,Mobs,cause_list,ppd)$plug)
# }
# if ((k_curr > 1 )){
# if (!use_jags){# use winbugs:
# chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
# chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_Nest_Slice(s,Mobs,cause_list)$plug)
# }else{#use jags:
# chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_Nest_Slice_jags(s,Mobs,cause_list,ppd)$plug)
# chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_Nest_Slice_jags(s,Mobs,cause_list)$plug)
# }
# chunk_BrS_subclass <- paste0(chunk_BrS_subclass, add_meas_BrS_subclass_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
# chunk_BrS_ctrl <- paste0(chunk_BrS_ctrl, add_meas_BrS_ctrl_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
# }
}# end iterate over slices.
}
if ("BrS" %in% use_measurements & !("SS" %in% use_measurements)) {
chunk <- paste0(
"
# BrS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# BrS measurement characteristics:
",chunk_BrS_param
)
}
if ("BrS" %in% use_measurements & ("SS" %in% use_measurements)) {
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# BrS and SS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# BrS and SS measurement characteristics:
",chunk_BrS_param,
add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
if (!("BrS" %in% use_measurements) & ("SS" %in% use_measurements)) {
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# SS measurements:
for (i in 1:Nd){
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
# SS measurement characteristics:
",add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
paste0(chunk,"\n")
}
#' Insert measurement likelihood (nested model+regression) code chunks into .bug model file
#'
#'
#' @param Mobs Measurement data in the form of `data_nplcm`
#' @param prior Prior specification from `model_options`
#' @param cause_list A list of latent status names (crucial for building templates;
#' see [make_template()])
#' @param FPR_formula A list of FPR regression formula; check `model_options$likelihood$FPR_formula`
#' @param use_measurements "BrS", or "SS"
#' @param ppd Default is NULL; set to TRUE for posterior predictive checking
#' @param use_jags Default is FALSE; set to TRUE if want to use JAGS for model fitting.
#'
#' @return A long character string to be inserted into .bug model file as measurement
#' likelihood
#'
#' @seealso Called by [write_model_Reg_NoNest] for constructing a `.bug` file.
#' This is usually called along with specification of latent status regression
#' ([insert_bugfile_chunk_reg_etiology]).
#'
insert_bugfile_chunk_reg_nest_meas <-
function(Mobs,prior,cause_list,FPR_formula,use_measurements = "BrS",ppd=NULL,use_jags=FALSE) {
if (!("BrS" %in% use_measurements) && !("SS" %in% use_measurements)){
stop("==[baker] No BrS or SS measurements specified in the model! ==")
}
# generate file:
if ("BrS" %in% use_measurements){
chunk_BrS_case <- ""
chunk_BrS_ctrl <- ""
chunk_BrS_subclass <- ""
chunk_BrS_param <- ""
for (s in seq_along(Mobs$MBS)){# begin iterate over slices:
# k_curr <- k_subclass[s]
# if (k_curr == 1 ){
if (!use_jags){ # use winbugs:
stop("==[baker] WinBUGS code not developed. Set 'use_jags=TRUE'. ==\n")
#chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_NoNest_reg_Slice_jags(s,Mobs,cause_list,ppd)$plug)
#chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_NoNest_reg_Slice_jags(s,Mobs,cause_list)$plug)
} else{# use jags:
# chunk_BrS_case is using case_Nest_Slice - no separate case_Nest_Slice_Reg needed.
chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_Nest_reg_Slice_jags(s,Mobs,prior,cause_list,FPR_formula)$plug)
}
# the following are using Nest_Slice - no separate Nest_Slice_Reg needed:
chunk_BrS_subclass <- paste0(chunk_BrS_subclass, add_meas_BrS_subclass_Nest_Slice(s,Mobs,cause_list,ppd,reg=TRUE)$plug)
chunk_BrS_ctrl <- paste0(chunk_BrS_ctrl, add_meas_BrS_ctrl_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
}# end iterate over slices.
}
if ("BrS" %in% use_measurements & !("SS" %in% use_measurements)){
chunk <- paste0(
"
# BrS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# BrS measurement characteristics:
",chunk_BrS_param
)
}
if ("BrS" %in% use_measurements & ("SS" %in% use_measurements)){
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# BrS and SS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# BrS and SS measurement characteristics:
",chunk_BrS_param,
add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
if (!("BrS" %in% use_measurements) & ("SS" %in% use_measurements)) {
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# SS measurements:
for (i in 1:Nd){
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
# SS measurement characteristics:
",add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
paste0(chunk,"\n")
}
#' insert etiology regression for latent status code chunk into .bug file
#'
#' @param Eti_formula Etiology regression formula; Check `model_options$likelihood$Eti_formula`.
#' @param ppd Default is NULL; set to TRUE for posterior predictive checking
#' @param Jcause The number of distinct causes, i.e., categories of latent health
#' status; equals `length(model_options$likelihood$cause_list)`.
#'
#' @return a long character string to be inserted into .bug model file
#' as distribution specification for latent status
#'
insert_bugfile_chunk_reg_etiology <- function(Eti_formula, Jcause, ppd = NULL){
constant_Eti <- is_intercept_only(Eti_formula)
# ER: etiology regression: #<---- add basis related operation here for PS.
ER_has_basis <- ifelse(length(grep("s_date",Eti_formula))==0,FALSE,TRUE)
ER_has_non_basis <- has_non_basis(Eti_formula)
ppd_seg <- ""
if (!is.null(ppd) && ppd){ppd_seg <-
"
Icat.new[i] ~ dcat(pEti[i,1:Jcause])
"}
if (!constant_Eti){
chunk_etiology <- paste0(
"
# etiology priors:
mu_Eti_mat <- Z_Eti%*%betaEti # <--- Z_Eti with rows for cases, columns for covariates; betaEti: rows for covariates, columns for 1:Jcause.
")
} else{
chunk_etiology <- paste0(
"
# etiology priors:
for (j in 1:Jcause){
mu_Eti_mat[1:Nd,j] <- Z_Eti*betaEti[1,j]
}
")
}
ER_basis_seg <- ""
if (ER_has_basis){
ER_basis_seg <- paste0(ER_basis_seg,
"
# B-spline basis coefficients for etiology regression:
betaEti[ER_basis_id[1],j] ~ dnorm(0,ER_prec_first)
for (l in 2:ER_n_basis){# iterate over the vector of B-spline basis:
betaEti[ER_basis_id[l],j] ~ dnorm(betaEti[ER_basis_id[l-1],j], ER_taubeta[j])
}
# select flexible semiparametric regression:
ER_taubeta0[j,1] ~ dgamma(3,2) ## <-- flexible fit.
ER_taubeta_inv[j] ~ dpar(1.5,0.0025) ## <--- constant fit.
ER_taubeta0[j,2] <- pow(ER_taubeta_inv[j],-1)
ER_flexible_select[j] ~ dbern(ER_p_flexible_select)
ER_ind_flex_select[j] <- 2- ER_flexible_select[j]
ER_taubeta[j] <- ER_taubeta0[j,ER_ind_flex_select[j]]
")
}
ER_non_basis_seg <- ""
if (ER_has_non_basis){
ER_non_basis_seg <- paste0(ER_non_basis_seg,
"
#non_basis_coefficients (e.g., intercept):
for (l in ER_non_basis_id){
betaEti[l,j] ~ dnorm(0,ER_prec_nonbasis) # <- shared with the ps ones.
}
")
}
if (Jcause > 2) {
chunk_etiology <- paste0(chunk_etiology,
"
for (i in 1:Nd){
Icat[i] ~ dcat(pEti[i,1:Jcause])
",
ppd_seg,
"for (j in 1:Jcause){
pEti[i,j] <- phi[i,j]/sum(phi[i,])
log(phi[i,j]) <- mu_Eti_mat[i,j]
}
}
",
c("","for (j in 1:Jcause){")[any(c(ER_has_basis,ER_has_non_basis))+1],
ER_basis_seg,
ER_non_basis_seg,
c("","}")[any(c(ER_has_basis,ER_has_non_basis))+1],
c("",
"#hyperprior of smoothness:
ER_p_flexible_select ~ dbeta(ER_alpha,ER_beta)
ER_prec_first <- pow(sd_betaEti_basis,-2) #1/4 #precision for ps coefficients - redundant if no ps bases.
")[ER_has_basis+1],
c("",
"ER_prec_nonbasis <- pow(sd_betaEti_nonbasis,-2) #1/4 #precision nonbasis coefficients
")[ER_has_non_basis+1]
#"for (p in 1:(ncol_dm_Eti)){",
# # betaEti[p,1:(Jcause-1)] ~ dmnorm(zero_Jcause_1,0.1*I_Jcause_1)
# "
# for (j in 1:Jcause){
# betaEti[p,j] ~ dnorm(0,1/sd_betaEti^2)
# }
#for (j in 1:(Jcause-1)){
# betaEti[p,j] ~ dnorm(0,1/sd_betaEti^2)
# }
# betaEti[p,Jcause] <- 0
#}"
)
} else{
chunk_etiology <- paste0(chunk_etiology,
"
for (i in 1:Nd){
Icat[i] ~ dcat(pEti[i,1:Jcause])
",
ppd_seg,
"for (j in 1:Jcause){
pEti[i,j] <- phi[i,j]/sum(phi[i,])
log(phi[i,j]) <- mu_Eti_mat[i,j]
}
}
for (p in 1:(ncol_dm_Eti)){
for (j in 1:Jcause){
betaEti[p,j] ~ dnorm(0,1/sd_betaEti^2)
}
#betaEti[p,1] ~ dnorm(0,1/2.25)
#betaEti[p,Jcause] <- 0
}"
)
}
paste0(chunk_etiology,"\n")
}
##############################################################################
############### REGRESSION MODELS: DISCRETE #################
##############################################################################
#' Insert measurement likelihood (with regression; discrete) code chunks into .bug model file
#'
#' @param Mobs Measurement data in the form of `data_nplcm`
#' @param prior Prior specification from `model_options`
#' @param cause_list A list of latent status names (crucial for building templates;
#' see [make_template()])
#' @param use_measurements "BrS", or "SS"
#' @param ppd Default is NULL; set to TRUE for posterior predictive checking
#' @param use_jags Default is FALSE; set to TRUE if want to use JAGS for model fitting.
#'
#' @return A long character string to be inserted into .bug model file as measurement
#' likelihood
#'
#' @seealso It is used in [write_model_Reg_NoNest] for constructing a .bug file along with
#' specification of latent status regression ([insert_bugfile_chunk_reg_etiology])
#'
#'
insert_bugfile_chunk_reg_discrete_predictor_nonest_meas <-
function(Mobs,prior,cause_list,use_measurements = "BrS",ppd=NULL,use_jags=FALSE) {
if (!("BrS" %in% use_measurements) && !("SS" %in% use_measurements)){
stop("==[baker] No BrS or SS measurements specified in the model! ==")
}
# for (s in seq_along(Mobs$MBS)){
# if (k_subclass[s]>1 && ncol(Mobs$MBS[[s]])==1){
# stop("==[baker] Cannot do nested modeling for BrS measurements with only one column! ==")
# }
# }
# generate file:
if ("BrS" %in% use_measurements){
chunk_BrS_case <- ""
chunk_BrS_ctrl <- ""
chunk_BrS_subclass <- ""
chunk_BrS_param <- ""
for (s in seq_along(Mobs$MBS)){# begin iterate over slices:
# k_curr <- k_subclass[s]
# if (k_curr == 1 ){
if (!use_jags){ # use winbugs:
stop("==[baker] WinBUGS code coming soon. ==\n")
#chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_NoNest_reg_Slice_jags(s,Mobs,cause_list,ppd)$plug)
#chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_NoNest_reg_Slice_jags(s,Mobs,cause_list)$plug)
} else{# use jags:
chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_NoNest_reg_discrete_predictor_Slice_jags(s,Mobs,prior,cause_list,ppd)$plug)
chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_NoNest_reg_discrete_predictor_Slice_jags(s,Mobs,prior,cause_list)$plug)
}
chunk_BrS_ctrl <- paste0(chunk_BrS_ctrl, add_meas_BrS_ctrl_NoNest_reg_discrete_predictor_Slice_jags(s,Mobs,cause_list,ppd)$plug)
# }
# if ((k_curr > 1 )){
# if (!use_jags){# use winbugs:
# chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
# chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_Nest_Slice(s,Mobs,cause_list)$plug)
# }else{#use jags:
# chunk_BrS_case <- paste0(chunk_BrS_case, add_meas_BrS_case_Nest_Slice_jags(s,Mobs,cause_list,ppd)$plug)
# chunk_BrS_param <- paste0(chunk_BrS_param, add_meas_BrS_param_Nest_Slice_jags(s,Mobs,cause_list)$plug)
# }
# chunk_BrS_subclass <- paste0(chunk_BrS_subclass, add_meas_BrS_subclass_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
# chunk_BrS_ctrl <- paste0(chunk_BrS_ctrl, add_meas_BrS_ctrl_Nest_Slice(s,Mobs,cause_list,ppd)$plug)
# }
}# end iterate over slices.
}
if ("BrS" %in% use_measurements & !("SS" %in% use_measurements)) {
chunk <- paste0(
"
# BrS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# BrS measurement characteristics:
",chunk_BrS_param
)
}
if ("BrS" %in% use_measurements & ("SS" %in% use_measurements)) {
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# BrS and SS measurements:
for (i in 1:Nd){
",chunk_BrS_case,"
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
for (i in (Nd+1):(Nd+Nu)){
",chunk_BrS_ctrl,"
}
",chunk_BrS_subclass,"
# BrS and SS measurement characteristics:
",chunk_BrS_param,
add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
if (!("BrS" %in% use_measurements) & ("SS" %in% use_measurements)) {
nslice_SS <- length(Mobs$MSS)
chunk <- paste0(
"
# SS measurements:
for (i in 1:Nd){
",add_meas_SS_case(nslice_SS,Mobs,prior,cause_list)$plug,"
}
# SS measurement characteristics:
",add_meas_SS_param(nslice_SS,Mobs,prior,cause_list)$plug
)
}
paste0(chunk,"\n")
}
#' insert etiology regression for latent status code chunk into .bug file; discrete predictors
#'
#' @param ppd Default is NULL; set to TRUE for posterior predictive checking
#' @param Jcause The number of distinct causes, i.e., categories of latent health
#' status; equals `length(model_options$likelihood$cause_list)`.
#'
#' @return a long character string to be inserted into .bug model file
#' as distribution specification for latent status
#'
#'
insert_bugfile_chunk_reg_discrete_predictor_etiology <- function(Jcause, ppd = NULL){
ppd_seg <- ""
if (!is.null(ppd) && ppd){ppd_seg <-
"
Icat.new[i] ~ dcat(pEti[Eti_stratum_id[i],1:Jcause])
"}
chunk_etiology <- paste0(
"
")
chunk_etiology <- paste0(chunk_etiology,
"
for (i in 1:Nd){
Icat[i] ~ dcat(pEti[Eti_stratum_id[i],1:Jcause])
",
ppd_seg,
"
}
for (s in 1:n_unique_Eti_level){
pEti[s, 1:Jcause]~ddirch(alphaEti[s,])
}
")
paste0(chunk_etiology,"\n")
}
##############################################################################
############### plug-and-play.R #################
##############################################################################
#
# 1. Bronze-standard data: conditional independence model (no nested structure).
#
#' add a likelihood component for a BrS measurement slice among cases (conditional independence)
#'
#' @param s the slice
#' @param Mobs See `data_nplcm` described in [nplcm()]
#' @param cause_list the list of causes in `data_nplcm` described in [nplcm()]
#' @param ppd Default is NULL; Set to TRUE for enabling posterior predictive checking.
#' @return a list of two elements: the first is `plug`, the .bug code;
#' the second is `parameters` that stores model parameters introduced by this
#' plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_case_NoNest_Slice <- function(s,Mobs,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_") #
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
templateBS_nm <- paste("templateBS",seq_along(BrS_nm),sep = "_")#
indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")#
Icat_nm <- "Icat"
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<-", indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-", indBS_nm[s],"[i,j])*",psiBS.cut_nm[s],"[j]
}
"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<-", indBS_nm[s],"[i]*",thetaBS_nm[s],"+(1-", indBS_nm[s],"[i])*",psiBS.cut_nm[s],"
"
)
}
parameters <- c(Icat_nm, thetaBS_nm[s],psiBS.cut_nm[s])
# if posterior predictive distribution is requested:
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
indBS_nm.new <- paste("indBS.new",seq_along(BrS_nm),sep = "_")#
Icat_nm.new <- "Icat.new"
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<-", indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-", indBS_nm[s],"[i,j])*",psiBS.cut_nm[s],"[j]
# posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i],j]
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j]<-", indBS_nm.new[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-", indBS_nm.new[s],"[i,j])*",psiBS.cut_nm[s],"[j]
}
"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<-", indBS_nm[s],"[i]*",thetaBS_nm[s],"+(1-", indBS_nm[s],"[i])*",psiBS.cut_nm[s],"
# posterior predictive distribution:
",indBS_nm.new[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i]]
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i]<-", indBS_nm.new[s],"[i]*",thetaBS_nm[s],"+(1-", indBS_nm.new[s],"[i])*",psiBS.cut_nm[s],"
"
)
}
parameters <- c(Icat_nm,Icat_nm.new, thetaBS_nm[s],psiBS.cut_nm[s])
}
make_list(plug,parameters)
}
#' add a likelihood component for a BrS measurement slice among controls (conditional independence)
#'
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_ctrl_NoNest_Slice <- function(s, Mobs,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")#
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
## control BrS measurements; no subclass:
for (j in 1:",JBrS_nm[s],"){
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<- ",psiBS_nm[s],"[j]
}
"
)
} else{
plug <- paste0(
"
## control BrS measurements; no subclass (only one column):
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<- ",psiBS_nm[s],"\n"
)
}
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
## control BrS measurements; no subclass:
for (j in 1:",JBrS_nm[s],"){
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<- ",psiBS_nm[s],"[j]
## posterior predictive distribution
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j]<- ",psiBS_nm[s],"[j]
}
"
)
} else{
plug <- paste0(
"
## control BrS measurements; no subclass (only one column):
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<- ",psiBS_nm[s],"
## posterior predictive distribution:
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i]<- ",psiBS_nm[s],"
"
)
}
}
parameters <- c(psiBS_nm[s])
make_list(plug,parameters)
}
#' add parameters for a BrS measurement slice among cases and controls (conditional independence)
#'
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_param_NoNest_Slice <- function(s,Mobs,cause_list) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")#
alphaB_nm <- paste("alphaB",seq_along(BrS_nm),sep = "_")#
betaB_nm <- paste("betaB",seq_along(BrS_nm),sep = "_")#
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
# BrS measurement characteristics - non-nested:
for (j in 1:",JBrS_nm[s],"){
",thetaBS_nm[s],"[j]~ dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
",psiBS_nm[s],"[j] ~ dbeta(1,1)
",psiBS.cut_nm[s],"[j]<-cut(",psiBS_nm[s],"[j])
}"
)
} else{
plug <-
paste0(
"
# BrS measurement characteristics - non-nested (only one column):
",thetaBS_nm[s],"~ dbeta(",alphaB_nm[s],",",betaB_nm[s],")
",psiBS_nm[s]," ~ dbeta(1,1)
",psiBS.cut_nm[s],"<-cut(",psiBS_nm[s],")\n"
)
}
parameters <- c(thetaBS_nm[s],psiBS_nm[s],alphaB_nm[s],betaB_nm[s])
make_list(plug,parameters)
}
#
# 2. Bronze-standard data: conditional dependence model (nested structure).
#
#' add likelihood for a BrS measurement slice among cases (conditional dependence)
#'
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
#'
add_meas_BrS_case_Nest_Slice <- function(s,Mobs,cause_list,ppd=NULL){
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs.bound_nm <- paste("mu_bs.bound",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
PR_BS_nm <- paste("PR_BS",seq_along(BrS_nm),sep = "_")#
PsiBS_nm <- paste("PsiBS",seq_along(BrS_nm),sep = "_")#
#PsiBS.cut_nm <- paste("PsiBS.cut",seq_along(BrS_nm),sep = "_")#
ThetaBS_nm <- paste("ThetaBS",seq_along(BrS_nm),sep="_")#
Z_nm <- paste("Z",seq_along(BrS_nm),sep="_")#
K_nm <- paste("K",seq_along(BrS_nm),sep="_")#
templateBS_nm <-
paste("templateBS",seq_along(BrS_nm),sep = "_")#
indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")#
Icat_nm <- "Icat"
if (length(patho_BrS_list[[s]]) == 1){stop("==cannot do nested modeling for BrS measurement with 1 dimension!==")}
plug <-
paste0(
"
## case BrS measurements; with subclasses:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j]~dbern(",mu_bs.bound_nm[s],"[i,j])
",mu_bs.bound_nm[s],"[i,j]<-max(0.000001,min(0.999999,",mu_bs_nm[s],"[i,j]))
",mu_bs_nm[s],"[i,j]<-",PR_BS_nm[s],"[i,j,",Z_nm[s],"[i]]
for (s in 1:",K_nm[s],"){
",PR_BS_nm[s],"[i,j,s]<-",PsiBS_nm[s],"[j,s]*(1-",indBS_nm[s],"[i,j])+",ThetaBS_nm[s],"[j,s]*",indBS_nm[s],"[i,j]
}
}
"
)
parameters <- c(Icat_nm,ThetaBS_nm[s])
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs.bound_nm.new <- paste("mu_bs.bound.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
PR_BS_nm.new <- paste("PR_BS.new",seq_along(BrS_nm),sep = "_")#
Z_nm.new <- paste("Z.new",seq_along(BrS_nm),sep="_")#
Icat_nm.new <- "Icat.new"
indBS_nm.new <- paste("indBS.new",seq_along(BrS_nm),sep = "_")#
plug <-
paste0(
"
## case BrS measurements; with subclasses:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j]~dbern(",mu_bs.bound_nm[s],"[i,j])
",mu_bs.bound_nm[s],"[i,j]<-max(0.000001,min(0.999999,",mu_bs_nm[s],"[i,j]))
",mu_bs_nm[s],"[i,j]<-",PR_BS_nm[s],"[i,j,",Z_nm[s],"[i]]
for (s in 1:",K_nm[s],"){
",PR_BS_nm[s],"[i,j,s]<-",PsiBS_nm[s],"[j,s]*(1-",indBS_nm[s],"[i,j])+",ThetaBS_nm[s],"[j,s]*",indBS_nm[s],"[i,j]
}
## posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i],j]
",MBS_nm.new[s],"[i,j]~dbern(",mu_bs.bound_nm.new[s],"[i,j])
",mu_bs.bound_nm.new[s],"[i,j]<-max(0.000001,min(0.999999,",mu_bs_nm.new[s],"[i,j]))
",mu_bs_nm.new[s],"[i,j]<-",PR_BS_nm.new[s],"[i,j,",Z_nm.new[s],"[i]]
for (s in 1:",K_nm[s],"){
",PR_BS_nm.new[s],"[i,j,s]<-",PsiBS_nm[s],"[j,s]*(1-",indBS_nm.new[s],"[i,j])+",ThetaBS_nm[s],"[j,s]*",indBS_nm.new[s],"[i,j]
}
}
"
)
parameters <- c(Icat_nm, Icat_nm.new, ThetaBS_nm[s])
}
make_list(plug,parameters)
}
#' add likelihood for a BrS measurement slice among controls (conditional independence)
#'
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
#'
add_meas_BrS_ctrl_Nest_Slice <- function(s, Mobs,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs.bound_nm <- paste("mu_bs.bound",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
Z_nm <- paste("Z",seq_along(BrS_nm),sep="_")#
PsiBS_nm <- paste("PsiBS",seq_along(BrS_nm),sep="_")#
#templateBS_nm <-
# paste("templateBS",seq_along(BrS_nm),sep = "_")
if (length(patho_BrS_list[[s]]) == 1){stop("==[baker]cannot do nested modeling for BrS measurement with 1 dimension!==")}
plug <- paste0(
"
## control BrS measurements; with subclasses:
for (j in 1:",JBrS_nm[s],"){
",MBS_nm[s],"[i,j]~dbern(",mu_bs.bound_nm[s],"[i,j])
",mu_bs.bound_nm[s],"[i,j] <-max(0.000001,min(0.999999,",mu_bs_nm[s],"[i,j]))
",mu_bs_nm[s],"[i,j]<-",PsiBS_nm[s],"[j,",Z_nm[s],"[i]]
}
"
)
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs.bound_nm.new <- paste("mu_bs.bound.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
Z_nm.new <- paste("Z.new",seq_along(BrS_nm),sep="_")#
plug <- paste0(
"
## control BrS measurements; with subclasses:
for (j in 1:",JBrS_nm[s],"){
",MBS_nm[s],"[i,j]~dbern(",mu_bs.bound_nm[s],"[i,j])
",mu_bs.bound_nm[s],"[i,j] <-max(0.000001,min(0.999999,",mu_bs_nm[s],"[i,j]))
",mu_bs_nm[s],"[i,j]<-",PsiBS_nm[s],"[j,",Z_nm[s],"[i]]
## posterior predictive distribution:
",MBS_nm.new[s],"[i,j]~dbern(",mu_bs.bound_nm.new[s],"[i,j])
",mu_bs.bound_nm.new[s],"[i,j] <-max(0.000001,min(0.999999,",mu_bs_nm.new[s],"[i,j]))
",mu_bs_nm.new[s],"[i,j]<-",PsiBS_nm[s],"[j,",Z_nm.new[s],"[i]]
}
"
)
}
parameters <- c(PsiBS_nm[s])
make_list(plug,parameters)
}
#' add parameters for a BrS measurement slice among cases and controls (conditional dependence)
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
#'
add_meas_BrS_param_Nest_Slice <- function(s,Mobs,cause_list) { #note: has separated case and controls subclass weights.
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
PsiBS.cut_nm <- paste("PsiBS.cut",seq_along(BrS_nm),sep = "_")#
ThetaBS_nm <- paste("ThetaBS",seq_along(BrS_nm),sep="_")#
PsiBS_nm <- paste("PsiBS",seq_along(BrS_nm),sep="_")#
Lambda_nm <- paste("Lambda",seq_along(BrS_nm),sep="_")#
Lambda0_nm <- paste("Lambda0",seq_along(BrS_nm),sep="_")#
r0_nm <- paste("r0",seq_along(BrS_nm),sep="_")#
alphadp0_nm <- paste("alphadp0",seq_along(BrS_nm),sep="_")#
alphadp0_case_nm <- paste("alphadp0_case",seq_along(BrS_nm),sep="_")# <--- cases' subclass weights.
Eta0_nm <- paste("Eta0",seq_along(BrS_nm),sep="_")#
r1_nm <- paste("r1",seq_along(BrS_nm),sep="_")#
Eta_nm <- paste("Eta",seq_along(BrS_nm),sep="_")#
alphaB_nm <- paste("alphaB",seq_along(BrS_nm),sep = "_")#
betaB_nm <- paste("betaB",seq_along(BrS_nm),sep = "_")#
K_nm <- paste("K",seq_along(BrS_nm),sep="_")#
templateBS_nm <-
paste("templateBS",seq_along(BrS_nm),sep = "_")
indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")
if (length(patho_BrS_list[[s]]) == 1){stop("==cannot do nested modeling for BrS measurement with 1 dimension!==")}
plug <- paste0(
"
## cut the feedback from case model to FPR:
for (j in 1:",JBrS_nm[s],"){
for (s in 1:",K_nm[s],"){
",PsiBS.cut_nm[s],"[j,s]<-cut(",PsiBS_nm[s],"[j,s])
}
}
####################################
### stick-breaking prior specification
####################################
# control subclass mixing weights:
",Lambda0_nm[s],"[1]<-",r0_nm[s],"[1]
",r0_nm[s],"[",K_nm[s],"]<-1
for(j in 2:",K_nm[s],") {",Lambda0_nm[s],"[j]<-",r0_nm[s],"[j]*(1-",r0_nm[s],"[j-1])*",Lambda0_nm[s],"[j-1]/",r0_nm[s],"[j-1]}
for(k in 1:(",K_nm[s],"-1)){
",r0_nm[s],"[k]~dbeta(1,",alphadp0_nm[s],")T(0.000001,0.999999)
}
for (k in 1:(",K_nm[s],"-1)){",Lambda_nm[s],"[k]<-max(0.000001,min(0.999999,",Lambda0_nm[s],"[k]))}
",Lambda_nm[s],"[",K_nm[s],"]<-1-sum(",Lambda_nm[s],"[1:(",K_nm[s],"-1)])
# case subclass mixing weights:
",Eta0_nm[s],"[1]<-",r1_nm[s],"[1]
",r1_nm[s],"[",K_nm[s],"]<-1
for(j in 2:",K_nm[s],") {",Eta0_nm[s],"[j]<-",r1_nm[s],"[j]*(1-",r1_nm[s],"[j-1])*",Eta0_nm[s],"[j-1]/",r1_nm[s],"[j-1]}
for(k in 1:(",K_nm[s],"-1)){
",r1_nm[s],"[k]~dbeta(1,",alphadp0_case_nm[s],")T(0.000001,0.999999)
}
for (k in 1:(",K_nm[s],"-1)){",Eta_nm[s],"[k]<-max(0.000001,min(0.999999,",Eta0_nm[s],"[k]))}
",Eta_nm[s],"[",K_nm[s],"]<-1-sum(",Eta_nm[s],"[1:(",K_nm[s],"-1)])
",alphadp0_nm[s],"~dgamma(.25,.25)I(0.001,20)
",alphadp0_case_nm[s],"~dgamma(.25,.25)I(0.001,20)
#########################
## priors on TPR and FPR:
#########################
for (j in 1:",JBrS_nm[s],"){
for (s in 1:",K_nm[s],"){
",PsiBS_nm[s],"[j,s]~dbeta(1,1)
#ThetaBS[j,s]~dbeta(1,1)
",ThetaBS_nm[s],"[j,s]~dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
}
}
"
)
parameters <- c(PsiBS_nm[s], ThetaBS_nm[s], Lambda_nm[s],Eta_nm[s],alphaB_nm[s],betaB_nm[s])
make_list(plug,parameters)
}
#' add subclass indicators for a BrS measurement slice among cases and controls (conditional independence)
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#' @param reg Default is NULL; set to TRUE if doing regression (double index of subclass weights: subject and subclass)
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_subclass_Nest_Slice <- function(s,Mobs,cause_list,ppd=NULL,reg=NULL){
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
Z_nm <- paste("Z",seq_along(BrS_nm),sep="_")#
Eta_nm <- paste("Eta",seq_along(BrS_nm),sep="_")#
Lambda_nm <- paste("Lambda",seq_along(BrS_nm),sep="_")#
K_nm <- paste("K",seq_along(BrS_nm),sep="_")#
if (is.null(reg)){
plug <- paste0(
"
# cases' subclass indicators:
for (i in 1:Nd){
",Z_nm[s],"[i] ~ dcat(",Eta_nm[s],"[1:",K_nm[s],"])
}
# controls' subclass indicators:
for (i in (Nd+1):(Nd+Nu)){
",Z_nm[s],"[i] ~ dcat(",Lambda_nm[s],"[1:",K_nm[s],"])
}
"
)
parameters <- c(Eta_nm[s],Lambda_nm[s],Z_nm[s])
if (!is.null(ppd) && ppd){
Z_nm.new <- paste("Z.new",seq_along(BrS_nm),sep="_")#
plug <- paste0(
"
# cases' subclass indicators:
for (i in 1:Nd){
",Z_nm[s],"[i] ~ dcat(",Eta_nm[s],"[1:",K_nm[s],"])
",Z_nm.new[s],"[i] ~ dcat(",Eta_nm[s],"[1:",K_nm[s],"])
}
# controls' subclass indicators:
for (i in (Nd+1):(Nd+Nu)){
",Z_nm[s],"[i] ~ dcat(",Lambda_nm[s],"[1:",K_nm[s],"])
",Z_nm.new[s],"[i] ~ dcat(",Lambda_nm[s],"[1:",K_nm[s],"])
}
"
)
parameters <- c(parameters,Z_nm.new[s])
}
return(make_list(plug,parameters))
}
if (!is.null(reg) && reg){ #if do regression:
plug <- paste0(
"
# cases' subclass indicators:
for (i in 1:Nd){
",Z_nm[s],"[i] ~ dcat(",Eta_nm[s],"[i,1:",K_nm[s],"])
}
# controls' subclass indicators:
for (i in (Nd+1):(Nd+Nu)){
",Z_nm[s],"[i] ~ dcat(",Lambda_nm[s],"[i,1:",K_nm[s],"])
}
"
)
parameters <- c(Eta_nm[s],Lambda_nm[s],Z_nm[s])
if (!is.null(ppd) && ppd){
Z_nm.new <- paste("Z.new",seq_along(BrS_nm),sep="_")#
plug <- paste0(
"
# cases' subclass indicators:
for (i in 1:Nd){
",Z_nm[s],"[i] ~ dcat(",Eta_nm[s],"[i,1:",K_nm[s],"])
",Z_nm.new[s],"[i] ~ dcat(",Eta_nm[s],"[i,1:",K_nm[s],"])
}
# controls' subclass indicators:
for (i in (Nd+1):(Nd+Nu)){
",Z_nm[s],"[i] ~ dcat(",Lambda_nm[s],"[i,1:",K_nm[s],"])
",Z_nm.new[s],"[i] ~ dcat(",Lambda_nm[s],"[i,1:",K_nm[s],"])
}
"
)
parameters <- c(parameters,Z_nm.new[s])
}
return(make_list(plug,parameters))
}
}
#
# 3. Silver standard data:
#
#' add likelihood for a SS measurement slice among cases (conditional independence)
#'
#'
#' @param nslice the total number of SS measurement slices
#' @param Mobs see `data_nplcm` described in [nplcm()]
#' @param prior see `model_options` described in [nplcm()]
#' @param cause_list the list of causes in `model_options` described in [nplcm()]
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_SS_case <- function(nslice,Mobs,prior,cause_list) {
# mapping template (by `make_template` function):
patho_SS_list <- lapply(Mobs$MSS,colnames)
template_SS_list <-
lapply(patho_SS_list,make_template,cause_list) # key.
# create variable names:
SS_nm <- names(Mobs$MSS)
# index measurement slices by numbers:
JSS_nm <- paste("JSS",seq_along(SS_nm),sep = "_")#
MSS_nm <- paste("MSS",seq_along(SS_nm),sep = "_")#
mu_ss_nm <- paste("mu_ss",seq_along(SS_nm),sep = "_")#
thetaSS_nm <- paste("thetaSS",seq_along(SS_nm),sep = "_")#
psiSS_nm <- paste("psiSS",seq_along(SS_nm),sep = "_")#
templateSS_nm <-
paste("templateSS",seq_along(SS_nm),sep = "_")#
indSS_nm <- paste("indSS",seq_along(SS_nm),sep = "_")#
# for SS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
SS_TPR_strat <- FALSE
SS_TPR_grp_nm <- paste("SS_TPR_grp",seq_along(SS_nm),sep = "_")
GSS_TPR_nm <- paste("GSS_TPR",seq_along(SS_nm),sep = "_") # level of groups within each slice.
prior_SS <- prior$TPR_prior$SS
if (!is.null(prior_SS$grp) && length(unique(prior_SS$grp)) >1 ){
SS_TPR_strat <- TRUE
}
res <- rep(NA,nslice)
for (s in 1:nslice) {
if (!SS_TPR_strat){# without stratified TPR in SS:
if (length(patho_SS_list[[s]]) > 1) { # if more dim>1.
res[s] <-
paste0(
"
# cases' SS measurements:
for (j in 1:",JSS_nm[s],"){
",indSS_nm[s],"[i,j] <- ",templateSS_nm[s],"[Icat[i],j]
",MSS_nm[s],"[i,j] ~ dbern(",mu_ss_nm[s],"[i,j])
",mu_ss_nm[s],"[i,j]<-", indSS_nm[s],"[i,j]*",thetaSS_nm[s],"[j]+(1-", indSS_nm[s],"[i,j])*",psiSS_nm[s],"[j]
}
")
} else{
res[s] <-
paste0(
"
# cases' SS measurements (only one column):
",indSS_nm[s],"[i] <- ",templateSS_nm[s],"[Icat[i]]
",MSS_nm[s],"[i] ~ dbern(",mu_ss_nm[s],"[i])
",mu_ss_nm[s],"[i]<-", indSS_nm[s],"[i]*",thetaSS_nm[s],"+(1-",indSS_nm[s],"[i])*",psiSS_nm[s],"\n"
)
}
} else{# WITH stratified TPR in SS:
if (length(patho_SS_list[[s]]) > 1) {
res[s] <-
paste0(
"
# cases' SS measurements:
for (j in 1:",JSS_nm[s],"){
",indSS_nm[s],"[i,j] <- ",templateSS_nm[s],"[Icat[i],j]
",MSS_nm[s],"[i,j] ~ dbern(",mu_ss_nm[s],"[i,j])
",mu_ss_nm[s],"[i,j]<-", indSS_nm[s],"[i,j]*",thetaSS_nm[s],"[",SS_TPR_grp_nm[s],"[i], j]+(1-", indSS_nm[s],"[i,j])*",psiSS_nm[s],"[j]
}
")
} else{
res[s] <-
paste0(
"
# cases' SS measurements (only one column):
",indSS_nm[s],"[i] <- ",templateSS_nm[s],"[Icat[i]]
",MSS_nm[s],"[i] ~ dbern(",mu_ss_nm[s],"[i])
",mu_ss_nm[s],"[i]<-", indSS_nm[s],"[i]*",thetaSS_nm[s],"[",SS_TPR_grp_nm[s],"[i]]+(1-", indSS_nm[s],"[i])*",psiSS_nm[s],"
"
)
}
}
}
plug <- paste0(res,collapse = "")
parameters <- c("Icat",thetaSS_nm[s],psiSS_nm[s])
make_list(plug,parameters)
}
#' add parameters for a SS measurement slice among cases (conditional independence)
#'
#'
#' @inheritParams add_meas_SS_case
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_SS_param <- function(nslice,Mobs,prior,cause_list) {
# mapping template (by `make_template` function):
patho_SS_list <- lapply(Mobs$MSS,colnames)
template_SS_list <-
lapply(patho_SS_list,make_template,cause_list) # key.
# create variable names:
SS_nm <- names(Mobs$MSS)
# index measurement slices by numbers:
JSS_nm <- paste("JSS",seq_along(SS_nm),sep = "_")#
thetaSS_nm <- paste("thetaSS",seq_along(SS_nm),sep = "_")#
psiSS_nm <- paste("psiSS",seq_along(SS_nm),sep = "_")#
alphaS_nm <- paste("alphaS",seq_along(SS_nm),sep = "_")#
betaS_nm <- paste("betaS",seq_along(SS_nm),sep = "_")#
# for SS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
SS_TPR_strat <- FALSE
SS_TPR_grp_nm <- paste("SS_TPR_grp",seq_along(SS_nm),sep = "_")
GSS_TPR_nm <- paste("GSS_TPR",seq_along(SS_nm),sep = "_") # level of groups within each slice.
prior_SS <- prior$TPR_prior$SS
if (!is.null(prior_SS$grp) && length(unique(prior_SS$grp)) >1 ){
SS_TPR_strat <- TRUE
}
res <- rep(NA,nslice)
for (s in 1:nslice) {
if (!SS_TPR_strat){# without stratified TPR in SS:
if (length(patho_SS_list[[s]]) > 1) {
res[s] <- paste0(
"
for (j in 1:",JSS_nm[s],"){
",thetaSS_nm[s],"[j]~dbeta(",alphaS_nm[s],"[j],",betaS_nm[s],"[j])
",psiSS_nm[s],"[j]<-0
","
}"
)
} else{
res[s] <-
paste0(
"
",thetaSS_nm[s],"~dbeta(",alphaS_nm[s],",",betaS_nm[s],")
",psiSS_nm[s],"<- 0
"
)
}
}else{# WITH stratified TPR in SS:
if (length(patho_SS_list[[s]]) > 1) {
res[s] <- paste0(
"
for (j in 1:",JSS_nm[s],"){
for (g in 1:",GSS_TPR_nm[s],"){
",thetaSS_nm[s],"[g,j]~dbeta(",alphaS_nm[s],"[g,j],",betaS_nm[s],"[g,j])
","
}
",psiSS_nm[s],"[j]<-0
}"
)
} else{
res[s] <-
paste0(
"
for (g in 1:",GSS_TPR_nm[s],"){
",thetaSS_nm[s],"[g]~dbeta(",alphaS_nm[s],"[g,1],",betaS_nm[s],"[g,1])
","
}
",psiSS_nm[s],"<- 0
"
)
}
}
}
plug <- paste0(res,collapse = "")
parameters <- c(thetaSS_nm[s],psiSS_nm[s], alphaS_nm[s],betaS_nm[s])
make_list(plug,parameters)
}
##
##
##
##
##
##
##
##
##
##
## FOR JAGS:
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
#
# 1. Bronze-standard data: conditional independence model (no nested structure).
#
#' add a likelihood component for a BrS measurement slice among cases (conditional independence)
#'
#' @param s the slice
#' @param Mobs See `data_nplcm` described in [nplcm()]
#' @param prior Prior specifications.
#' @param cause_list the list of causes in `data_nplcm` described in [nplcm()]
#' @param ppd Default is NULL; Set to TRUE for enabling posterior predictive checking.
#' @return a list of two elements: the first is `plug`, the .bug code;
#' the second is `parameters` that stores model parameters introduced by this
#' plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_case_NoNest_Slice_jags <- function(s,Mobs,prior,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_") #
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")#
#psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
templateBS_nm <- paste("templateBS",seq_along(BrS_nm),sep = "_")#
indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")#
Icat_nm <- "Icat"
# for BrS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
BrS_TPR_strat <- FALSE
BrS_TPR_grp_nm <- paste("BrS_TPR_grp",seq_along(BrS_nm),sep = "_")
GBrS_TPR_nm <- paste("GBrS_TPR",seq_along(BrS_nm),sep = "_") # level of groups within each slice.
prior_BrS <- prior$TPR_prior$BrS
if (!is.null(prior_BrS$grp) && length(unique(prior_BrS$grp)) >1 ){
BrS_TPR_strat <- TRUE
}
if (!BrS_TPR_strat){ # no TPR stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<-", indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-", indBS_nm[s],"[i,j])*",psiBS_nm[s],"[j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<-", indBS_nm[s],"[i]*",thetaBS_nm[s],"+(1-", indBS_nm[s],"[i])*",psiBS_nm[s],"\n"
)
}
} else{ # TPR stratified.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<-", indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j]+(1-", indBS_nm[s],"[i,j])*",psiBS_nm[s],"[j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<-", indBS_nm[s],"[i]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]]+(1-", indBS_nm[s],"[i])*",psiBS_nm[s],"\n"
)
}
}
parameters <- c(Icat_nm, thetaBS_nm[s],psiBS_nm[s])
# if posterior predictive distribution is requested:
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
indBS_nm.new <- paste("indBS.new",seq_along(BrS_nm),sep = "_")#
Icat_nm.new <- "Icat.new"
if (BrS_TPR_strat){ # TPR stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<-", indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j]+(1-", indBS_nm[s],"[i,j])*",psiBS_nm[s],"[j]
# posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i],j]
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j]<-", indBS_nm.new[s],"[i,j]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j]+(1-", indBS_nm.new[s],"[i,j])*",psiBS_nm[s],"[j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<-", indBS_nm[s],"[i]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]]+(1-", indBS_nm[s],"[i])*",psiBS_nm[s],"
# posterior predictive distribution:
",indBS_nm.new[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i]]
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i]<-", indBS_nm.new[s],"[i]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]]+(1-", indBS_nm.new[s],"[i])*",psiBS_nm[s],"
\n"
)
}
} else{ # without TPR stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j]<-", indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-", indBS_nm[s],"[i,j])*",psiBS_nm[s],"[j]
# posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i],j]
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j]<-", indBS_nm.new[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-", indBS_nm.new[s],"[i,j])*",psiBS_nm[s],"[j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i]<-", indBS_nm[s],"[i]*",thetaBS_nm[s],"+(1-", indBS_nm[s],"[i])*",psiBS_nm[s],"
# posterior predictive distribution:
",indBS_nm.new[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i]]
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i]<-", indBS_nm.new[s],"[i]*",thetaBS_nm[s],"+(1-", indBS_nm.new[s],"[i])*",psiBS_nm[s],"
\n"
)
}
}
parameters <- c(Icat_nm,Icat_nm.new, thetaBS_nm[s],psiBS_nm[s])
}
make_list(plug,parameters)
}
#' add parameters for a BrS measurement slice among cases and controls (conditional independence)
#'
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#' @param prior Prior specifications.
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_param_NoNest_Slice_jags <- function(s,Mobs,prior,cause_list) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
#psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")#
alphaB_nm <- paste("alphaB",seq_along(BrS_nm),sep = "_")#
betaB_nm <- paste("betaB",seq_along(BrS_nm),sep = "_")#
# for BrS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
BrS_TPR_strat <- FALSE
BrS_TPR_grp_nm <- paste("BrS_TPR_grp",seq_along(BrS_nm),sep = "_")
GBrS_TPR_nm <- paste("GBrS_TPR",seq_along(BrS_nm),sep = "_") # level of groups within each slice.
prior_BrS <- prior$TPR_prior$BrS
if (!is.null(prior_BrS$grp) && length(unique(prior_BrS$grp)) >1 ){
BrS_TPR_strat <- TRUE
}
if (!BrS_TPR_strat){ # no TPR stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
# BrS measurement characteristics - non-nested:
for (j in 1:",JBrS_nm[s],"){
",thetaBS_nm[s],"[j]~ dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
",psiBS_nm[s],"[j] ~ dbeta(1,1)
}"
)
} else{
plug <-
paste0(
"
# BrS measurement characteristics - non-nested (only one column):
",thetaBS_nm[s],"~ dbeta(",alphaB_nm[s],",",betaB_nm[s],")
",psiBS_nm[s]," ~ dbeta(1,1)
\n"
)
}
}else{ # with TPR stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
# BrS measurement characteristics - non-nested:
for (j in 1:",JBrS_nm[s],"){
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g,j]~ dbeta(",alphaB_nm[s],"[g,j],",betaB_nm[s],"[g,j])
}
",psiBS_nm[s],"[j] ~ dbeta(1,1)
}"
)
} else{
plug <-
paste0(
"
# BrS measurement characteristics - non-nested (only one column):
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g]~ dbeta(",alphaB_nm[s],"[g,1],",betaB_nm[s],"[g,1])
}
",psiBS_nm[s]," ~ dbeta(1,1)
\n"
)
}
}
parameters <- c(thetaBS_nm[s],psiBS_nm[s],alphaB_nm[s],betaB_nm[s])
make_list(plug,parameters)
}
#
# 2. Bronze-standard data: conditional dependence model (nested structure).
#
#' add likelihood for a BrS measurement slice among cases (conditional dependence)
#'
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#'
#' @return a list of two elements: the first is `plug`, the `.bug` code;
#' the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
#'
add_meas_BrS_case_Nest_Slice_jags <- function(s,Mobs,cause_list,ppd=NULL){
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs.bound_nm <- paste("mu_bs.bound",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
PR_BS_nm <- paste("PR_BS",seq_along(BrS_nm),sep = "_")#
#PsiBS.cut_nm <- paste("PsiBS.cut",seq_along(BrS_nm),sep = "_")#
PsiBS_nm <- paste("PsiBS",seq_along(BrS_nm),sep = "_")#
ThetaBS_nm <- paste("ThetaBS",seq_along(BrS_nm),sep="_")#
Z_nm <- paste("Z",seq_along(BrS_nm),sep="_")#
K_nm <- paste("K",seq_along(BrS_nm),sep="_")#
templateBS_nm <-
paste("templateBS",seq_along(BrS_nm),sep = "_")#
indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")#
Icat_nm <- "Icat"
if (length(patho_BrS_list[[s]]) == 1){stop("==cannot do nested modeling for BrS measurement with 1 dimension!==")}
plug <-
paste0(
"
## case BrS measurements; with subclasses:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j]~dbern(",mu_bs.bound_nm[s],"[i,j])
",mu_bs.bound_nm[s],"[i,j]<-max(0.000001,min(0.999999,",mu_bs_nm[s],"[i,j]))
",mu_bs_nm[s],"[i,j]<-",PR_BS_nm[s],"[i,j,",Z_nm[s],"[i]]
for (s in 1:",K_nm[s],"){
",PR_BS_nm[s],"[i,j,s]<-",PsiBS_nm[s],"[j,s]*(1-",indBS_nm[s],"[i,j])+",ThetaBS_nm[s],"[j,s]*",indBS_nm[s],"[i,j]
}
}
"
)
parameters <- c(Icat_nm,ThetaBS_nm[s])
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs.bound_nm.new <- paste("mu_bs.bound.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
PR_BS_nm.new <- paste("PR_BS.new",seq_along(BrS_nm),sep = "_")#
Z_nm.new <- paste("Z.new",seq_along(BrS_nm),sep="_")#
Icat_nm.new <- "Icat.new"
indBS_nm.new <- paste("indBS.new",seq_along(BrS_nm),sep = "_")#
plug <-
paste0(
"
## case BrS measurements; with subclasses:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j]~dbern(",mu_bs.bound_nm[s],"[i,j])
",mu_bs.bound_nm[s],"[i,j]<-max(0.000001,min(0.999999,",mu_bs_nm[s],"[i,j]))
",mu_bs_nm[s],"[i,j]<-",PR_BS_nm[s],"[i,j,",Z_nm[s],"[i]]
for (s in 1:",K_nm[s],"){
",PR_BS_nm[s],"[i,j,s]<-",PsiBS_nm[s],"[j,s]*(1-",indBS_nm[s],"[i,j])+",ThetaBS_nm[s],"[j,s]*",indBS_nm[s],"[i,j]
}
## posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new,"[i],j]
",MBS_nm.new[s],"[i,j]~dbern(",mu_bs.bound_nm.new[s],"[i,j])
",mu_bs.bound_nm.new[s],"[i,j]<-max(0.000001,min(0.999999,",mu_bs_nm.new[s],"[i,j]))
",mu_bs_nm.new[s],"[i,j]<-",PR_BS_nm.new[s],"[i,j,",Z_nm.new[s],"[i]]
for (s in 1:",K_nm[s],"){
",PR_BS_nm.new[s],"[i,j,s]<-",PsiBS_nm[s],"[j,s]*(1-",indBS_nm.new[s],"[i,j])+",ThetaBS_nm[s],"[j,s]*",indBS_nm.new[s],"[i,j]
}
}
"
)
parameters <- c(Icat_nm, Icat_nm.new, ThetaBS_nm[s])
}
make_list(plug,parameters)
}
#' add parameters for a BrS measurement slice among cases and controls (conditional dependence)
#'
#'
#' @inheritParams add_meas_BrS_case_NoNest_Slice
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_param_Nest_Slice_jags <- function(s,Mobs,cause_list) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
#PsiBS.cut_nm <- paste("PsiBS.cut",seq_along(BrS_nm),sep = "_")#
ThetaBS_nm <- paste("ThetaBS",seq_along(BrS_nm),sep="_")#
PsiBS_nm <- paste("PsiBS",seq_along(BrS_nm),sep="_")#
Lambda_nm <- paste("Lambda",seq_along(BrS_nm),sep="_")#
Lambda0_nm <- paste("Lambda0",seq_along(BrS_nm),sep="_")#
r0_nm <- paste("r0",seq_along(BrS_nm),sep="_")#
alphadp0_nm <- paste("alphadp0",seq_along(BrS_nm),sep="_")#
alphadp0_case_nm <- paste("alphadp0_case",seq_along(BrS_nm),sep="_")# <--- added case subclass weights.
Eta0_nm <- paste("Eta0",seq_along(BrS_nm),sep="_")#
r1_nm <- paste("r1",seq_along(BrS_nm),sep="_")#
Eta_nm <- paste("Eta",seq_along(BrS_nm),sep="_")#
alphaB_nm <- paste("alphaB",seq_along(BrS_nm),sep = "_")#
betaB_nm <- paste("betaB",seq_along(BrS_nm),sep = "_")#
K_nm <- paste("K",seq_along(BrS_nm),sep="_")#
#templateBS_nm <-
# paste("templateBS",seq_along(BrS_nm),sep = "_")
#indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")
if (length(patho_BrS_list[[s]]) == 1){stop("==cannot do nested modeling for BrS measurement with 1 dimension!==")}
plug <- paste0(
"
####################################
### stick-breaking prior specification
####################################
# control subclass mixing weights:
",Lambda0_nm[s],"[1]<-",r0_nm[s],"[1]
",r0_nm[s],"[",K_nm[s],"]<-1
for(j in 2:",K_nm[s],") {",Lambda0_nm[s],"[j]<-",r0_nm[s],"[j]*(1-",r0_nm[s],"[j-1])*",Lambda0_nm[s],"[j-1]/",r0_nm[s],"[j-1]}
for(k in 1:(",K_nm[s],"-1)){
",r0_nm[s],"[k]~dbeta(1,",alphadp0_nm[s],")T(0.000001,0.999999)
}
for (k in 1:(",K_nm[s],"-1)){",Lambda_nm[s],"[k]<-max(0.000001,min(0.999999,",Lambda0_nm[s],"[k]))}
",Lambda_nm[s],"[",K_nm[s],"]<-1-sum(",Lambda_nm[s],"[1:(",K_nm[s],"-1)])
# case subclass mixing weights:
",Eta0_nm[s],"[1]<-",r1_nm[s],"[1]
",r1_nm[s],"[",K_nm[s],"]<-1
for(j in 2:",K_nm[s],") {",Eta0_nm[s],"[j]<-",r1_nm[s],"[j]*(1-",r1_nm[s],"[j-1])*",Eta0_nm[s],"[j-1]/",r1_nm[s],"[j-1]}
for(k in 1:(",K_nm[s],"-1)){
",r1_nm[s],"[k]~dbeta(1,",alphadp0_case_nm[s],")T(0.000001,0.999999)
}
for (k in 1:(",K_nm[s],"-1)){",Eta_nm[s],"[k]<-max(0.000001,min(0.999999,",Eta0_nm[s],"[k]))}
",Eta_nm[s],"[",K_nm[s],"]<-1-sum(",Eta_nm[s],"[1:(",K_nm[s],"-1)])
",alphadp0_nm[s],"~dgamma(.25,.25)T(0.001,20)
",alphadp0_case_nm[s],"~dgamma(.25,.25)T(0.001,20)
#########################
## priors on TPR and FPR:
#########################
for (j in 1:",JBrS_nm[s],"){
for (s in 1:",K_nm[s],"){
",PsiBS_nm[s],"[j,s]~dbeta(1,1)
#ThetaBS[j,s]~dbeta(1,1)
",ThetaBS_nm[s],"[j,s]~dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
}
}
"
)
parameters <- c(PsiBS_nm[s], ThetaBS_nm[s], Lambda_nm[s],Eta_nm[s],alphaB_nm[s],betaB_nm[s])
make_list(plug,parameters)
}
##############################################################################
############### PLUG-AND-PLAY for REGRESSION MODELS: NON-DISCRETE#################
##############################################################################
#' add likelihood component for a BrS measurement slice among cases
#'
#' regression model with no nested subclasses
#'
#' @param s the slice
#' @param Mobs See `data_nplcm` described in [nplcm()]
#' @param prior Prior specifications.
#' @param cause_list the list of causes in `data_nplcm` described in [nplcm()]
#' @param ppd Default is NULL; Set to TRUE for enabling posterior predictive checking.
#' @return a list of two elements: the first is `plug`, the .bug code;
#' the second is `parameters` that stores model parameters introduced by this
#' plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_case_NoNest_reg_Slice_jags <- function(s,Mobs,prior,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_") #
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
#psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
templateBS_nm <- paste("templateBS",seq_along(BrS_nm),sep = "_")#
indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")#
Icat_nm <- "Icat"
linpred_psiBS_nm <- paste("linpred_psiBS",seq_along(BrS_nm),sep = "_")
# for BrS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
BrS_TPR_strat <- FALSE
BrS_TPR_grp_nm <- paste("BrS_TPR_grp",seq_along(BrS_nm),sep = "_")
GBrS_TPR_nm <- paste("GBrS_TPR",seq_along(BrS_nm),sep = "_") # level of groups within each slice.
prior_BrS <- prior$TPR_prior$BrS
if (!is.null(prior_BrS$grp) && length(unique(prior_BrS$grp)) >1 ){
BrS_TPR_strat <- TRUE
}
if (!BrS_TPR_strat){ # no stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurements; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ilogit(",indBS_nm[s],"[i,j]*logit(",thetaBS_nm[s],"[j])+(1-",indBS_nm[s],"[i,j])*",linpred_psiBS_nm[s],"[i,j])
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <- ilogit(",indBS_nm[s],"[i]*logit(",thetaBS_nm[s],")+(1-",indBS_nm[s],"[i])*",linpred_psiBS_nm[s],"[i,1])
\n"
)
}
} else{ # with stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurements; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ilogit(",indBS_nm[s],"[i,j]*logit(",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j])+(1-",indBS_nm[s],"[i,j])*",linpred_psiBS_nm[s],"[i,j])
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <- ilogit(",indBS_nm[s],"[i]*logit(",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]])+(1-",indBS_nm[s],"[i])*",linpred_psiBS_nm[s],"[i,1])
\n"
)
}
}
parameters <- c(Icat_nm, thetaBS_nm[s],linpred_psiBS_nm[s])
# if posterior predictive distribution is requested:
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
indBS_nm.new <- paste("indBS.new",seq_along(BrS_nm),sep = "_")#
Icat_nm.new <- "Icat.new"
if (!BrS_TPR_strat){ # no stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ilogit(",indBS_nm[s],"[i,j]*logit(",thetaBS_nm[s],"[j])+(1-",indBS_nm[s],"[i,j])*",linpred_psiBS_nm[s],"[i,j])
# posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i],j]
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j] <- ilogit(",indBS_nm.new[s],"[i,j]*logit(",thetaBS_nm[s],"[j])+(1-",indBS_nm.new[s],"[i,j])*",linpred_psiBS_nm[s],"[i,j])
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <-ilogit( ",indBS_nm[s],"[i]*logit(",thetaBS_nm[s],")+(1-",indBS_nm[s],"[i])*",linpred_psiBS_nm[s],"[i,1])
# posterior predictive distribution:
",indBS_nm.new[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i]]
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i] <- ilogit(",indBS_nm.new[s],"[i]*logit(",thetaBS_nm[s],")+(1-",indBS_nm.new[s],"[i])*",linpred_psiBS_nm[s],"[i,1])
"
)
}
} else{# with stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ilogit(",indBS_nm[s],"[i,j]*logit(",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j])+(1-",indBS_nm[s],"[i,j])*",linpred_psiBS_nm[s],"[i,j])
# posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i],j]
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j] <- ilogit(",indBS_nm.new[s],"[i,j]*logit(",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j])+(1-",indBS_nm.new[s],"[i,j])*",linpred_psiBS_nm[s],"[i,j])
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <- ilogit(",indBS_nm[s],"[i]*logit(",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]])+(1-",indBS_nm[s],"[i])*",linpred_psiBS_nm[s],"[i,1])
# posterior predictive distribution:
",indBS_nm.new[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i]]
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i] <- ilogit(",indBS_nm.new[s],"[i]*logit(",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]])+(1-",indBS_nm.new[s],"[i])*",linpred_psiBS_nm[s],"[i,1])
"
)
}
}
parameters <- c(Icat_nm,Icat_nm.new, thetaBS_nm[s],linpred_psiBS_nm[s])
}
make_list(plug,parameters)
}
#' add parameters for a BrS measurement slice among cases and controls
#'
#' regression model with no nested subclasses
#'
#' @inheritParams add_meas_BrS_case_NoNest_reg_Slice_jags
#' @param FPR_formula False positive regression formula for slice s of BrS data.
#' Check `model_options$likelihood$FPR_formula[[s]]`.
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_param_NoNest_reg_Slice_jags <- function(s,Mobs,prior,cause_list,FPR_formula) {
constant_FPR <- is_intercept_only(FPR_formula[[s]])
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
has_basis <- ifelse(length(grep("s_date",FPR_formula[[s]]))==0,FALSE,TRUE)
exists_non_basis <- has_non_basis(FPR_formula[[s]])
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
#psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
alphaB_nm <- paste("alphaB",seq_along(BrS_nm),sep = "_")#
betaB_nm <- paste("betaB",seq_along(BrS_nm),sep = "_")#
betaFPR_nm <- paste("betaFPR",seq_along(BrS_nm),sep = "_")#
basis_id_nm <- paste("basis_id",seq_along(BrS_nm),sep = "_")#
non_basis_id_nm <- paste("non_basis_id",seq_along(BrS_nm),sep = "_")#
n_basis_nm <- paste("n_basis",seq_along(BrS_nm),sep = "_")#
prec_first_nm <- paste("prec_first",seq_along(BrS_nm),sep = "_")#
prec_non_basis_nm <- paste("prec_non_basis_nm",seq_along(BrS_nm),sep = "_")#
taubeta_nm <- paste("taubeta",seq_along(BrS_nm),sep = "_")#
taubeta0_nm <- paste("taubeta0",seq_along(BrS_nm),sep = "_")#
taubeta_inv_nm <- paste("taubeta_inv",seq_along(BrS_nm),sep = "_")#
flexible_select_nm <- paste("flexible_select",seq_along(BrS_nm),sep = "_")#
ind_flex_select_nm <- paste("ind_flex_select",seq_along(BrS_nm),sep = "_")#
p_flexible_nm <- paste("p_flexible",seq_along(BrS_nm),sep = "_")#
linpred_psiBS_nm <- paste("linpred_psiBS",seq_along(BrS_nm),sep = "_")#
Z_FPR_nm <- paste("Z_FPR",seq_along(BrS_nm),sep = "_")#
I_JBrS_nm <- paste("I_JBrS",seq_along(BrS_nm),sep = "_")#
zero_JBrS_nm <- paste("zero_JBrS",seq_along(BrS_nm),sep = "_")#
# for BrS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
BrS_TPR_strat <- FALSE
BrS_TPR_grp_nm <- paste("BrS_TPR_grp",seq_along(BrS_nm),sep = "_")
GBrS_TPR_nm <- paste("GBrS_TPR",seq_along(BrS_nm),sep = "_") # level of groups within each slice.
prior_BrS <- prior$TPR_prior$BrS
if (!is.null(prior_BrS$grp) && length(unique(prior_BrS$grp)) >1 ){
BrS_TPR_strat <- TRUE
}
if (BrS_TPR_strat){ # with stratification.
if (length(patho_BrS_list[[s]]) > 1) { # <--- if the dimension is higher than 2:
if (!constant_FPR){
plug <- paste0(
"
# BrS measurement characteristics - non-nested:
",linpred_psiBS_nm[s]," <- ",Z_FPR_nm[s],"%*%",betaFPR_nm[s]," # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
for (j in 1:",JBrS_nm[s],"){
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g,j]~ dbeta(",alphaB_nm[s],"[g,j],",betaB_nm[s],"[g,j])
}
")
}else{
plug <- paste0(
"
for (j in 1:",JBrS_nm[s],"){
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g,j]~ dbeta(",alphaB_nm[s],"[g,j],",betaB_nm[s],"[g,j])
}
# BrS measurement characteristics - non-nested:
",linpred_psiBS_nm[s],"[1:(Nd+Nu),j] <- ",Z_FPR_nm[s],"*",betaFPR_nm[s],"[1,j] # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
")
}
if(has_basis){
plug <- paste0(plug,
"
# B-spline basis coefficients:
",#betaFPR_nm[s],"[",basis_id_nm[s],"[1],1:",JBrS_nm[s],"] ~ dmnorm(",zero_JBrS_nm[s],",",prec_first_nm[s],")
betaFPR_nm[s],"[",basis_id_nm[s],"[1]",",j] ~ dnorm(0,",prec_first_nm[s],")
for (l in 2:",n_basis_nm[s],"){# iterate over the vector of B-spline basis.
",betaFPR_nm[s],"[",basis_id_nm[s],"[l],j] ~ dnorm(",betaFPR_nm[s],"[",basis_id_nm[s],"[l-1],j],",taubeta_nm[s],"[j])
}
# select flexible semiparametric regression:
",taubeta0_nm[s],"[j,1] ~ dgamma(3,2) # <-------- flexible fit.
",taubeta_inv_nm[s],"[j] ~ dpar(1.5,0.0025) # <--------constant fit.
",taubeta0_nm[s],"[j,2] <- pow(",taubeta_inv_nm[s],"[j],-1)
",flexible_select_nm[s],"[j] ~ dbern(",p_flexible_nm[s],")
",ind_flex_select_nm[s],"[j] <- 2-",flexible_select_nm[s],"[j]
",taubeta_nm[s],"[j] <- ",taubeta0_nm[s],"[j,",ind_flex_select_nm[s],"[j]]
}
# hyperprior of smoothness:
",p_flexible_nm[s]," ~ dbeta(3,3)#flexible prob
",#prec_first_nm[s]," <- 1/4*",I_JBrS_nm[s],"
prec_first_nm[s]," <- pow(sd_betaFPR_basis,-2) #1/4 #precision for spline coefficients
")
} else{
plug <- paste0(plug,
"
}
")
}
if (exists_non_basis){
plug <- paste0(plug,
"
for (l in ",non_basis_id_nm[s],"){
",#betaFPR_nm[s],"[l,1:",JBrS_nm[s],"] ~ dmnorm(",zero_JBrS_nm[s],",",prec_first_nm[s],")
"for (j in 1:",JBrS_nm[s],"){
",betaFPR_nm[s],"[l,j] ~ dnorm(0,",prec_non_basis_nm[s],")
}
}
",#prec_first_nm[s]," <- 1/4*",I_JBrS_nm[s],"
prec_non_basis_nm[s]," <- pow(sd_betaFPR_nonbasis,-2) #1/4 #precision for spline coefficients
"
)
}
} else{ # <-- if the dimension equals 1:
if (!constant_FPR){
plug <-
paste0(
"
",linpred_psiBS_nm[s],"[1:(Nd+Nu),1] <- ",Z_FPR_nm[s],"%*%",betaFPR_nm[s],"[,1] # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
# BrS measurement characteristics - non-nested (only one column):
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g]~ dbeta(",alphaB_nm[s],"[g,1],",betaB_nm[s],"[g,1])
}
")
}else{
plug <-
paste0(
"
# BrS measurement characteristics - non-nested (only one column):
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g]~ dbeta(",alphaB_nm[s],"[g,1],",betaB_nm[s],"[g,1])
}
",linpred_psiBS_nm[s],"[1:(Nd+Nu),1] <- ",Z_FPR_nm[s],"*",betaFPR_nm[s],"[,1] # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
")
}
if(has_basis){
plug <-
paste0(plug,
"
# B-spline basis coefficients:
",betaFPR_nm[s],"[",basis_id_nm[s],"[1],1] ~ dnorm(0,",prec_first_nm[s],")
for (l in 2:",n_basis_nm[s],"){# iterate over the vector of B-spline basis.
",betaFPR_nm[s],"[",basis_id_nm[s],"[l],1] ~ dnorm(",betaFPR_nm[s],"[",basis_id_nm[s],"[l-1],1],",taubeta_nm[s],")
}
# select flexible semiparametric regression:
",taubeta0_nm[s],"[1] ~ dgamma(3,2) # <-------- flexible fit.
",taubeta_inv_nm[s]," ~ dpar(1.5,0.0025) # <--------constant fit.
",taubeta0_nm[s],"[2] <- pow(",taubeta_inv_nm[s],",-1)
",flexible_select_nm[s]," ~ dbern(",p_flexible_nm[s],")
",ind_flex_select_nm[s]," <- 2-",flexible_select_nm[s],"
",taubeta_nm[s]," <- ",taubeta0_nm[s],"[",ind_flex_select_nm[s],"]
# hyperprior of smoothness:
",p_flexible_nm[s]," ~ dbeta(3,3)#flexible prob
"
,prec_first_nm[s]," <- pow(sd_betaFPR_basis,-2) #1/4 #precision for spline coefficients
"
)
}
if (exists_non_basis){
plug <- # issue: can we test if we have non_basis coefficients, so then we include the following segment?
paste0(plug,
"
# non-basis coefficients (e.g., intercept):
for (l in ",non_basis_id_nm[s],"){
",betaFPR_nm[s],"[l,1] ~ dnorm(0,",prec_non_basis_nm[s],")
}
",prec_non_basis_nm[s]," <- pow(sd_betaFPR_nonbasis,-2) #1/4 #precision for spline coefficients
"
)
}
}
} else{ # no stratification.
if (length(patho_BrS_list[[s]]) > 1) { # <--- if the dimension is higher than 2:
if (!constant_FPR){
plug <- paste0(
"
# BrS measurement characteristics - non-nested:
",linpred_psiBS_nm[s]," <- ",Z_FPR_nm[s],"%*%",betaFPR_nm[s]," # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
for (j in 1:",JBrS_nm[s],"){
",thetaBS_nm[s],"[j]~ dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
")
}else{
plug <- paste0(
"
for (j in 1:",JBrS_nm[s],"){
",thetaBS_nm[s],"[j]~ dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
# BrS measurement characteristics - non-nested:
",linpred_psiBS_nm[s],"[1:(Nd+Nu),j] <- ",Z_FPR_nm[s],"*",betaFPR_nm[s],"[1,j] # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
")
}
if(has_basis){
plug <- paste0(plug,
"
# B-spline basis coefficients:
",#betaFPR_nm[s],"[",basis_id_nm[s],"[1],1:",JBrS_nm[s],"] ~ dmnorm(",zero_JBrS_nm[s],",",prec_first_nm[s],")
betaFPR_nm[s],"[",basis_id_nm[s],"[1]",",j] ~ dnorm(0,",prec_first_nm[s],")
for (l in 2:",n_basis_nm[s],"){# iterate over the vector of B-spline basis.
",betaFPR_nm[s],"[",basis_id_nm[s],"[l],j] ~ dnorm(",betaFPR_nm[s],"[",basis_id_nm[s],"[l-1],j],",taubeta_nm[s],"[j])
}
# select flexible semiparametric regression:
",taubeta0_nm[s],"[j,1] ~ dgamma(3,2) # <-------- flexible fit.
",taubeta_inv_nm[s],"[j] ~ dpar(1.5,0.0025) # <--------constant fit.
",taubeta0_nm[s],"[j,2] <- pow(",taubeta_inv_nm[s],"[j],-1)
",flexible_select_nm[s],"[j] ~ dbern(",p_flexible_nm[s],")
",ind_flex_select_nm[s],"[j] <- 2-",flexible_select_nm[s],"[j]
",taubeta_nm[s],"[j] <- ",taubeta0_nm[s],"[j,",ind_flex_select_nm[s],"[j]]
}
# hyperprior of smoothness:
",p_flexible_nm[s]," ~ dbeta(3,3)#flexible prob
",#prec_first_nm[s]," <- 1/4*",I_JBrS_nm[s],"
prec_first_nm[s]," <- pow(sd_betaFPR_basis,-2) #1/4 #precision for spline coefficients
")
} else{
plug <- paste0(plug,
"
}
")
}
if (exists_non_basis){
plug <- paste0(plug,
"
for (l in ",non_basis_id_nm[s],"){
",#betaFPR_nm[s],"[l,1:",JBrS_nm[s],"] ~ dmnorm(",zero_JBrS_nm[s],",",prec_first_nm[s],")
"for (j in 1:",JBrS_nm[s],"){
",betaFPR_nm[s],"[l,j] ~ dnorm(0,",prec_non_basis_nm[s],")
}
}
",#prec_first_nm[s]," <- 1/4*",I_JBrS_nm[s],"
prec_non_basis_nm[s]," <- pow(sd_betaFPR_nonbasis,-2) #1/4 #precision for spline coefficients
"
)
}
} else{ # <-- if the dimension equals 1:
if (!constant_FPR){
plug <-
paste0(
"
",linpred_psiBS_nm[s],"[1:(Nd+Nu),1] <- ",Z_FPR_nm[s],"%*%",betaFPR_nm[s],"[,1] # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
# BrS measurement characteristics - non-nested (only one column):
",thetaBS_nm[s],"~ dbeta(",alphaB_nm[s],",",betaB_nm[s],")
")
}else{
plug <-
paste0(
"
# BrS measurement characteristics - non-nested (only one column):
",thetaBS_nm[s],"~ dbeta(",alphaB_nm[s],",",betaB_nm[s],")
",linpred_psiBS_nm[s],"[1:(Nd+Nu),1] <- ",Z_FPR_nm[s],"*",betaFPR_nm[s],"[,1] # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
")
}
if(has_basis){
plug <-
paste0(plug,
"
# B-spline basis coefficients:
",betaFPR_nm[s],"[",basis_id_nm[s],"[1],1] ~ dnorm(0,",prec_first_nm[s],")
for (l in 2:",n_basis_nm[s],"){# iterate over the vector of B-spline basis.
",betaFPR_nm[s],"[",basis_id_nm[s],"[l],1] ~ dnorm(",betaFPR_nm[s],"[",basis_id_nm[s],"[l-1],1],",taubeta_nm[s],")
}
# select flexible semiparametric regression:
",taubeta0_nm[s],"[1] ~ dgamma(3,2) # <-------- flexible fit.
",taubeta_inv_nm[s]," ~ dpar(1.5,0.0025) # <--------constant fit.
",taubeta0_nm[s],"[2] <- pow(",taubeta_inv_nm[s],",-1)
",flexible_select_nm[s]," ~ dbern(",p_flexible_nm[s],")
",ind_flex_select_nm[s]," <- 2-",flexible_select_nm[s],"
",taubeta_nm[s]," <- ",taubeta0_nm[s],"[",ind_flex_select_nm[s],"]
# hyperprior of smoothness:
",p_flexible_nm[s]," ~ dbeta(3,3)#flexible prob
",prec_first_nm[s]," <- pow(sd_betaFPR_basis,-2) #1/4 #precision for spline coefficients
")
}
if (exists_non_basis){
plug <-
paste0(plug,
"
# non-basis coefficients:
for (l in ",non_basis_id_nm[s],"){
",betaFPR_nm[s],"[l,1] ~ dnorm(0,",prec_non_basis_nm[s],")
}
",prec_non_basis_nm[s]," <- pow(sd_betaFPR_nonbasis,-2) #1/4 #precision for spline coefficients
"
)
}
}
}
parameters <- c(thetaBS_nm[s],betaFPR_nm[s],alphaB_nm[s],betaB_nm[s],taubeta_nm[s],p_flexible_nm[s])
make_list(plug,parameters)
}
#' add a likelihood component for a BrS measurement slice among controls
#'
#' regression model without nested subclasses
#'
#' @inheritParams add_meas_BrS_case_NoNest_reg_Slice_jags
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_ctrl_NoNest_reg_Slice_jags <- function(s, Mobs,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")#
linpred_psiBS_nm <- paste("linpred_psiBS",seq_along(BrS_nm),sep = "_")#
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
## control BrS measurements; no subclass:
for (j in 1:",JBrS_nm[s],"){
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
logit(",mu_bs_nm[s],"[i,j]) <- ",linpred_psiBS_nm[s],"[i,j]
}
"
)
} else{
plug <- paste0(
"
## control BrS measurements; no subclass (only one column):
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
logit(",mu_bs_nm[s],"[i]) <- ",linpred_psiBS_nm[s],"[i,1]
"
)
}
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
## control BrS measurements; no subclass:
for (j in 1:",JBrS_nm[s],"){
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
logit(",mu_bs_nm[s],"[i,j]) <- ",linpred_psiBS_nm[s],"[i,j]
## posterior predictive distribution
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
}
"
)
} else{
plug <- paste0(
"
## control BrS measurements; no subclass (only one column):
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
logit(",mu_bs_nm[s],"[i]) <- ",linpred_psiBS_nm[s],"[i,1]
## posterior predictive distribution:
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
"
)
}
}
parameters <- c(psiBS_nm[s])
make_list(plug,parameters)
}
##############################################################################
############### PLUG-AND-PLAY for REGRESSION MODELS: NESTED ################
############### only needs to add parameters.
##############################################################################
#' add parameters for a BrS measurement slice among cases and controls
#'
#' regression model with nested subclasses; called by [insert_bugfile_chunk_reg_nest_meas]
#'
#' @inheritParams add_meas_BrS_case_NoNest_reg_Slice_jags
#' @param FPR_formula False positive regression formula for slice s of BrS data.
#' Check `model_options$likelihood$FPR_formula[[s]]`.
#'
#' @return a list of two elements: the first is `plug`, the .bug code;
#' the second is `parameters` that stores model parameters introduced by
#' this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_param_Nest_reg_Slice_jags <- function(s,Mobs,prior,cause_list,FPR_formula=NULL) {
#constant_FPR <- is_intercept_only(FPR_formula[[s]]) # formula for the subclass weights.
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
has_basis <- ifelse(length(grep("s_date",FPR_formula[[s]]))==0,FALSE,TRUE)
exists_non_basis <- has_non_basis(FPR_formula[[s]])
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
ThetaBS_nm <- paste("ThetaBS",seq_along(BrS_nm),sep = "_")#
PsiBS_nm <- paste("PsiBS",seq_along(BrS_nm),sep="_")#
#psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
alphaB_nm <- paste("alphaB",seq_along(BrS_nm),sep = "_")# for TPR.
betaB_nm <- paste("betaB",seq_along(BrS_nm),sep = "_")# for TPR.
# the number of subclasses:
K_nm <- paste("K",seq_along(BrS_nm),sep="_")#
# control - subclass weight regression:
Lambda_nm <- paste("Lambda",seq_along(BrS_nm),sep="_")# this is Lambda0 truncated to 0.000001 ~ 0.999999.
#Lambda0_nm <- paste("Lambda0",seq_along(BrS_nm),sep="_")# every control has an Eta0.
r0_nm <- paste("r0",seq_along(BrS_nm),sep="_")#
betaFPR_nm <- paste("betaFPR",seq_along(BrS_nm),sep = "_")# coeffients
basis_id_nm <- paste("basis_id",seq_along(BrS_nm),sep = "_")# which ones are penalized splines basis for one variable.
non_basis_id_nm <- paste("non_basis_id",seq_along(BrS_nm),sep = "_")#which ones are not.
n_basis_nm <- paste("n_basis",seq_along(BrS_nm),sep = "_")# the number of basis expansions in p-spline for one variable.
prec_first_nm <- paste("prec_first",seq_along(BrS_nm),sep = "_")# inv variance for the first coef in p-spline.
prec_non_basis_nm <- paste("prec_non_basis",seq_along(BrS_nm),sep = "_")# inv var for coefs not in p-spline.
taubeta0_nm <- paste("taubeta0",seq_along(BrS_nm),sep = "_")#inv vars for constant and flexible curves (to be selected below).
taubeta_nm <- paste("taubeta",seq_along(BrS_nm),sep = "_")# the inv var for adjacent p-spline basis coefs.
taubeta_inv_nm <- paste("taubeta_inv",seq_along(BrS_nm),sep = "_")# inv of pareto.
flexible_select_nm <- paste("flexible_select",seq_along(BrS_nm),sep = "_")#
ind_flex_select_nm <- paste("ind_flex_select",seq_along(BrS_nm),sep = "_")#
p_flexible_nm <- paste("p_flexible",seq_along(BrS_nm),sep = "_")#
#linpred_wt_nm <- paste("linpred_wt",seq_along(BrS_nm),sep = "_")#
Z_FPR_nm <- paste("Z_FPR",seq_along(BrS_nm),sep = "_")#
mu_ctrl_nm <- paste("mu_ctrl",seq_along(BrS_nm),sep = "_")#
mu0_ctrl_nm <- paste("mu0_ctrl",seq_along(BrS_nm),sep = "_")#
inv_scale_mu0_ctrl_nm <- paste("inv_scale_mu0_ctrl",seq_along(BrS_nm),sep = "_")#
half_s2_ctrl_nm <- paste("half_s2_ctrl",seq_along(BrS_nm),sep = "_")#
half_nu_ctrl_nm <- paste("half_nu_ctrl",seq_along(BrS_nm),sep = "_")#
# case - subclass weight regression:
Eta_nm <- paste("Eta",seq_along(BrS_nm),sep="_")# this is Eta0 truncated to 0.000001 ~ 0.999999.
#Eta0_nm <- paste("Eta0",seq_along(BrS_nm),sep="_")# every case has an Eta0.
r1_nm <- paste("r1",seq_along(BrS_nm),sep="_")#
case_betaFPR_nm <- paste("case_betaFPR",seq_along(BrS_nm),sep = "_")# coeffients
#prec_first_nm <- paste("case_prec_first",seq_along(BrS_nm),sep = "_")# inv variance for the first coef in p-spline.
#prec_non_basis_nm <- paste("case_prec_non_basis_nm",seq_along(BrS_nm),sep = "_")# inv var for coefs not in p-spline.
case_taubeta0_nm <- paste("case_taubeta0",seq_along(BrS_nm),sep = "_")#inv vars for constant and flexible curves (to be selected below).
case_taubeta_nm <- paste("case_taubeta",seq_along(BrS_nm),sep = "_")# the inv var for adjacent p-spline basis coefs.
case_taubeta_inv_nm <- paste("case_taubeta_inv",seq_along(BrS_nm),sep = "_")# inv of pareto.
case_flexible_select_nm <- paste("case_flexible_select",seq_along(BrS_nm),sep = "_")#
case_ind_flex_select_nm <- paste("case_ind_flex_select",seq_along(BrS_nm),sep = "_")#
case_p_flexible_nm <- paste("case_p_flexible",seq_along(BrS_nm),sep = "_")#
#case_linpred_wt_nm <- paste("case_linpred_wt",seq_along(BrS_nm),sep = "_")#
mu_case_nm <- paste("mu_case",seq_along(BrS_nm),sep = "_")#
#mu0_case_nm <- paste("mu0_case",seq_along(BrS_nm),sep = "_")#
#inv_scale_mu0_case_nm <- paste("inv_scale_mu0_case",seq_along(BrS_nm),sep = "_")#
#half_s2_case_nm <- paste("half_s2_case",seq_along(BrS_nm),sep = "_")#
#Z_FPR_nm <- paste("Z_FPR",seq_along(BrS_nm),sep = "_")# we let case and control have the same design matrices.
d_FPR_nm <- paste("d_FPR",seq_along(BrS_nm),sep = "_")#
# for BrS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
BrS_TPR_strat <- FALSE
BrS_TPR_grp_nm <- paste("BrS_TPR_grp",seq_along(BrS_nm),sep = "_")
GBrS_TPR_nm <- paste("GBrS_TPR",seq_along(BrS_nm),sep = "_") # level of groups within each slice.
prior_BrS <- prior$TPR_prior$BrS
if (!is.null(prior_BrS$grp) && length(unique(prior_BrS$grp)) >1 ){
BrS_TPR_strat <- TRUE
}
if (BrS_TPR_strat){stop("==[baker] cannot do true positive rate stratification for nested regression.
Use non nested regression now. And check back later. ==\n")}
# no stratification:
if (length(patho_BrS_list[[s]]) == 1){stop("==[baker] cannot do nested regression for
a slice of bronze-standard data (only 1 dimensional measure). ==\n")}
plug <- paste0(
"
# BrS measurement characteristics - nested:
",mu_ctrl_nm[s]," <- ",Z_FPR_nm[s],"%*%",betaFPR_nm[s]," # <--- Z_FPR_1: rows for cases and controls, columns for covariates; betaFPR_1: rows for covariates, columns for 1:JBrS_1, i.e., pathogens.
",mu_case_nm[s]," <- ",Z_FPR_nm[s],"%*%",case_betaFPR_nm[s],"
for (j in 1:",JBrS_nm[s],"){ # begin iterations over bronze-standard measures:
#########################
## priors on TPR and FPR:
#########################
for (s in 1:",K_nm[s],"){
",PsiBS_nm[s],"[j,s] ~ dbeta(1,1)
",ThetaBS_nm[s],"[j,s]~ dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
}
}# end iterations over bronze-standard measures.
###########################
## stick breaking prior for each person's probability of falling into K subclasses
##########################
for (i in 1:Nd){
",Eta_nm[s],"[i,1] <- ",r1_nm[s],"[i,1]
",r1_nm[s],"[i,",K_nm[s],"] <- 1
for (k in 2:",K_nm[s],"){",Eta_nm[s],"[i,k] <- ",r1_nm[s],"[i,k]*(1-",r1_nm[s],"[i,k-1])*",Eta_nm[s],"[i,k-1]/",r1_nm[s],"[i,k-1]}
for (j in 1:(",K_nm[s],"-1)){",r1_nm[s],"[i,j] <- max(0.000001,min(0.999999,ilogit(sum(",mu0_ctrl_nm,"[1:j])+",mu_case_nm,"[i,j])))} # <--- prevent extreme values that makes the division above undefined.
}
for (i in (Nd+1):(Nd+Nu)){
",Lambda_nm[s],"[i,1] <- ",r0_nm[s],"[i,1]
",r0_nm[s],"[i,",K_nm[s],"] <- 1
for (k in 2:",K_nm[s],"){",Lambda_nm[s],"[i,k] <- ",r0_nm[s],"[i,k]*(1-",r0_nm[s],"[i,k-1])*",Lambda_nm[s],"[i,k-1]/",r0_nm[s],"[i,k-1]}
for (j in 1:(",K_nm[s],"-1)){",r0_nm[s],"[i,j] <- max(0.000001,min(0.999999,ilogit(sum(",mu0_ctrl_nm,"[1:j])+",mu_ctrl_nm,"[i,j])))} # <--- prevent extreme values that makes the division above undefined.
}
for (j in 1:(",K_nm[s],"-1)){
",mu0_ctrl_nm[s],"[j] ~ dnorm(0,",inv_scale_mu0_ctrl_nm[s],"[j])T(0,) # can have penalty in the intercepts whether or not we have basis
",inv_scale_mu0_ctrl_nm[s],"[j] ~ dgamma(",half_nu_ctrl_nm[s],",",half_s2_ctrl_nm[s],")
}
#",half_s2_ctrl_nm[s]," ~ dgamma(2,0.04)
for (l in 1:",d_FPR_nm[s],"){
",betaFPR_nm[s],"[l,",K_nm[s],"] <- 0
",case_betaFPR_nm[s],"[l,",K_nm[s],"] <- 0
}
")
if(has_basis){
plug <- paste0(plug,
"
# BrS measurement characteristics - nested:
for (j in 1:(",K_nm[s],"-1)){
## control: B-spline basis coefficients:
",#betaFPR_nm[s],"[",basis_id_nm[s],"[1],1:",JBrS_nm[s],"] ~ dmnorm(",zero_JBrS_nm[s],",",prec_first_nm[s],")
betaFPR_nm[s],"[",basis_id_nm[s],"[1]",",j] ~ dnorm(0,",prec_first_nm[s],")
for (l in 2:",n_basis_nm[s],"){# iterate over the vector of B-spline basis.
",betaFPR_nm[s],"[",basis_id_nm[s],"[l],j] ~ dnorm(",betaFPR_nm[s],"[",basis_id_nm[s],"[l-1],j],",taubeta_nm[s],"[j])
}
# select flexible semiparametric regression:
",taubeta0_nm[s],"[j,1] ~ dgamma(3,2) # <-------- flexible fit.
",taubeta_inv_nm[s],"[j] ~ dpar(1.5,0.0025) # <--------constant fit.
",taubeta0_nm[s],"[j,2] <- pow(",taubeta_inv_nm[s],"[j],-1)
",flexible_select_nm[s],"[j] ~ dbern(",p_flexible_nm[s],")
",ind_flex_select_nm[s],"[j] <- 2-",flexible_select_nm[s],"[j]
",taubeta_nm[s],"[j] <- ",taubeta0_nm[s],"[j,",ind_flex_select_nm[s],"[j]]
## case: B-spline basis coefficients:
",#betaFPR_nm[s],"[",basis_id_nm[s],"[1],1:",JBrS_nm[s],"] ~ dmnorm(",zero_JBrS_nm[s],",",prec_first_nm[s],")
case_betaFPR_nm[s],"[",basis_id_nm[s],"[1]",",j] ~ dnorm(0,",prec_first_nm[s],")
for (l in 2:",n_basis_nm[s],"){# iterate over the vector of B-spline basis.
",case_betaFPR_nm[s],"[",basis_id_nm[s],"[l],j] ~ dnorm(",case_betaFPR_nm[s],"[",basis_id_nm[s],"[l-1],j],",case_taubeta_nm[s],"[j])
}
# select flexible semiparametric regression:
",case_taubeta0_nm[s],"[j,1] ~ dgamma(3,2) # <-------- flexible fit.
",case_taubeta_inv_nm[s],"[j] ~ dpar(1.5,0.0025) # <--------constant fit.
",case_taubeta0_nm[s],"[j,2] <- pow(",case_taubeta_inv_nm[s],"[j],-1)
",case_flexible_select_nm[s],"[j] ~ dbern(",case_p_flexible_nm[s],")
",case_ind_flex_select_nm[s],"[j] <- 2-",case_flexible_select_nm[s],"[j]
",case_taubeta_nm[s],"[j] <- ",case_taubeta0_nm[s],"[j,",case_ind_flex_select_nm[s],"[j]]
}
# control: hyperprior of smoothness:
",p_flexible_nm[s]," ~ dbeta(ctrl_flex_alpha,ctrl_flex_beta)#flexible prob
",#prec_first_nm[s]," <- 1/4*",I_JBrS_nm[s],"
prec_first_nm[s]," <- pow(sd_betaFPR_basis,-2) #1/4 #precision for spline coefficients
# case: hyperprior of smoothness:
",case_p_flexible_nm[s]," ~ dbeta(case_flex_alpha,case_flex_beta)#flexible prob
"
)
}
if (exists_non_basis){
plug <- paste0(plug,
"
for (l in ",non_basis_id_nm[s],"){
",#betaFPR_nm[s],"[l,1:",JBrS_nm[s],"] ~ dmnorm(",zero_JBrS_nm[s],",",prec_first_nm[s],")
"for (j in 1:(",K_nm[s],"-1)){
",betaFPR_nm[s],"[l,j] ~ dnorm(0,",prec_non_basis_nm[s],") # control.
",case_betaFPR_nm[s],"[l,j] ~ dnorm(0,",prec_non_basis_nm[s],") # case.
}
}
",#prec_first_nm[s]," <- 1/4*",I_JBrS_nm[s],"
prec_non_basis_nm[s]," <- pow(sd_betaFPR_nonbasis,-2) #1/4 #precision for spline coefficients
"
)
}
# NB: need to add parameters:
if (has_basis){
parameters <- c(Lambda_nm[s],taubeta_nm[s],p_flexible_nm[s],flexible_select_nm[s],mu_ctrl_nm[s],betaFPR_nm[s],mu0_ctrl_nm[s],half_s2_ctrl_nm[s],#mu0_case_nm[s],
Eta_nm[s],case_taubeta_nm[s],case_p_flexible_nm[s],case_flexible_select_nm[s],mu_case_nm[s],case_betaFPR_nm[s])
}else{
parameters <- c(Lambda_nm[s],mu_ctrl_nm[s],betaFPR_nm[s],mu0_ctrl_nm[s],half_s2_ctrl_nm[s],#mu0_case_nm[s],
Eta_nm[s],mu_case_nm[s],case_betaFPR_nm[s])
}
make_list(plug,parameters)
}
##############################################################################
############### PLUG-AND-PLAY for REGRESSION MODELS: DISCRETE################
##############################################################################
#' add likelihood component for a BrS measurement slice among cases
#'
#' regression model with no nested subclasses; discrete predictors
#'
#' @param s the slice
#' @param Mobs See `data_nplcm` described in [nplcm()]
#' @param prior Prior specifications.
#' @param cause_list the list of causes in `data_nplcm` described in [nplcm()]
#' @param ppd Default is NULL; Set to TRUE for enabling posterior predictive checking.
#' @return a list of two elements: the first is `plug`, the .bug code;
#' the second is `parameters` that stores model parameters introduced by this
#' plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_case_NoNest_reg_discrete_predictor_Slice_jags <- function(s,Mobs,prior,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_") #
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
#psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
templateBS_nm <- paste("templateBS",seq_along(BrS_nm),sep = "_")#
indBS_nm <- paste("indBS",seq_along(BrS_nm),sep = "_")#
Icat_nm <- "Icat"
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")
FPR_stratum_id_nm <- paste("FPR_stratum_id",seq_along(BrS_nm),sep = "_")
# for BrS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
BrS_TPR_strat <- FALSE
BrS_TPR_grp_nm <- paste("BrS_TPR_grp",seq_along(BrS_nm),sep = "_")
GBrS_TPR_nm <- paste("GBrS_TPR",seq_along(BrS_nm),sep = "_") # level of groups within each slice.
prior_BrS <- prior$TPR_prior$BrS
if (!is.null(prior_BrS$grp) && length(unique(prior_BrS$grp)) >1 ){
BrS_TPR_strat <- TRUE
}
if (!BrS_TPR_strat){ # no stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurements; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ",indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-",indBS_nm[s],"[i,j])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <- ",indBS_nm[s],"[i]*",thetaBS_nm[s],"+(1-",indBS_nm[s],"[i])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
\n"
)
}
} else{ # with stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurements; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ",indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j]+(1-",indBS_nm[s],"[i,j])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <- ",indBS_nm[s],"[i]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]]+(1-",indBS_nm[s],"[i])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
\n"
)
}
}
parameters <- c(Icat_nm, thetaBS_nm[s],psiBS_nm[s])
# if posterior predictive distribution is requested:
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
indBS_nm.new <- paste("indBS.new",seq_along(BrS_nm),sep = "_")#
Icat_nm.new <- "Icat.new"
if (!BrS_TPR_strat){ # no stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ",indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-",indBS_nm[s],"[i,j])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
# posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i],j]
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j] <- ",indBS_nm.new[s],"[i,j]*",thetaBS_nm[s],"[j]+(1-",indBS_nm.new[s],"[i,j])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <- ",indBS_nm[s],"[i]*",thetaBS_nm[s],"+(1-",indBS_nm[s],"[i])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
# posterior predictive distribution:
",indBS_nm.new[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i]]
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i] <- ",indBS_nm.new[s],"[i]*",thetaBS_nm[s],"+(1-",indBS_nm.new[s],"[i])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
"
)
}
} else{# with stratification.
if (length(patho_BrS_list[[s]]) > 1) {
plug <-
paste0(
"
# case BrS measurement; non-nested:
for (j in 1:",JBrS_nm[s],"){
",indBS_nm[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm,"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
",mu_bs_nm[s],"[i,j] <- ",indBS_nm[s],"[i,j]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j]+(1-",indBS_nm[s],"[i,j])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
# posterior predictive distribution:
",indBS_nm.new[s],"[i,j] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i],j]
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm.new[s],"[i,j])
",mu_bs_nm.new[s],"[i,j] <- ",indBS_nm.new[s],"[i,j]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i],j]+(1-",indBS_nm.new[s],"[i,j])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
}","\n"
)
} else{
plug <-
paste0(
"
# case BrS measurement; non-nested (with only one column):
",indBS_nm[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm,"[i]]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
",mu_bs_nm[s],"[i] <- ",indBS_nm[s],"[i]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]]+(1-",indBS_nm[s],"[i])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
# posterior predictive distribution:
",indBS_nm.new[s],"[i] <- ",templateBS_nm[s],"[",Icat_nm.new[s],"[i]]
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm.new[s],"[i])
",mu_bs_nm.new[s],"[i] <- ",indBS_nm.new[s],"[i]*",thetaBS_nm[s],"[",BrS_TPR_grp_nm[s],"[i]]+(1-",indBS_nm.new[s],"[i])*",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
"
)
}
}
parameters <- c(Icat_nm,Icat_nm.new, thetaBS_nm[s],psiBS_nm[s])
}
make_list(plug,parameters)
}
#' add parameters for a BrS measurement slice among cases and controls
#'
#' regression model with no nested subclasses; discrete
#'
#' @inheritParams add_meas_BrS_case_NoNest_reg_Slice_jags
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_param_NoNest_reg_discrete_predictor_Slice_jags <- function(s,Mobs,prior,cause_list) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
thetaBS_nm <- paste("thetaBS",seq_along(BrS_nm),sep = "_")#
#psiBS.cut_nm <- paste("psiBS.cut",seq_along(BrS_nm),sep = "_")#
alphaB_nm <- paste("alphaB",seq_along(BrS_nm),sep = "_")#
betaB_nm <- paste("betaB",seq_along(BrS_nm),sep = "_")#
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")#
n_unique_FPR_level_nm <- paste("n_unique_FPR_level",seq_along(BrS_nm),sep = "_")#
# for BrS TPR across groups (currently only allows uniform grouping across dimensions, i.e.,
# only allow the same way of splitting cases for every pathogen):
BrS_TPR_strat <- FALSE
BrS_TPR_grp_nm <- paste("BrS_TPR_grp",seq_along(BrS_nm),sep = "_")
GBrS_TPR_nm <- paste("GBrS_TPR",seq_along(BrS_nm),sep = "_") # level of groups within each slice.
prior_BrS <- prior$TPR_prior$BrS
if (!is.null(prior_BrS$grp) && length(unique(prior_BrS$grp)) >1 ){
BrS_TPR_strat <- TRUE
}
if (BrS_TPR_strat){ # with stratification.
if (length(patho_BrS_list[[s]]) > 1) { # <--- if the dimension is higher than 2:
plug <- paste0(
"
# BrS measurement characteristics - non-nested:
for (j in 1:",JBrS_nm[s],"){
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g,j]~ dbeta(",alphaB_nm[s],"[g,j],",betaB_nm[s],"[g,j])
}
for (s in 1:", n_unique_FPR_level_nm[s],"){
",psiBS_nm[s],"[s,j] ~ dbeta(1,1)
}
}
")
} else{ # <-- if the dimension equals 1:
plug <-
paste0(
"
# BrS measurement characteristics - non-nested (only one column):
for (g in 1:",GBrS_TPR_nm[s],"){
",thetaBS_nm[s],"[g]~ dbeta(",alphaB_nm[s],"[g,1],",betaB_nm[s],"[g,1])
}
for (s in 1:", n_unique_FPR_level_nm[s],"){
",psiBS_nm[s],"[s,1] ~ dbeta(1,1)
}
")
}
} else{ # no stratification.
if (length(patho_BrS_list[[s]]) > 1) { # <--- if the dimension is higher than 2:
plug <- paste0(
"
# BrS measurement characteristics - non-nested:
for (j in 1:",JBrS_nm[s],"){
",thetaBS_nm[s],"[j]~ dbeta(",alphaB_nm[s],"[j],",betaB_nm[s],"[j])
for (s in 1:", n_unique_FPR_level_nm[s],"){
",psiBS_nm[s],"[s,j] ~ dbeta(1,1)
}
}
")
} else{ # <-- if the dimension equals 1:
plug <-
paste0(
"
# BrS measurement characteristics - non-nested (only one column):
",thetaBS_nm[s],"~ dbeta(",alphaB_nm[s],",",betaB_nm[s],")
for (s in 1:", n_unique_FPR_level_nm[s],"){
",psiBS_nm[s],"[s,1] ~ dbeta(1,1)
}
")
}
}
parameters <- c(thetaBS_nm[s],psiBS_nm[s],alphaB_nm[s],betaB_nm[s])
make_list(plug,parameters)
}
#' add a likelihood component for a BrS measurement slice among controls
#'
#' regression model without nested subclasses; discrete
#'
#' @inheritParams add_meas_BrS_case_NoNest_reg_Slice_jags
#'
#' @return a list of two elements: the first is `plug`, the .bug code; the second is `parameters`
#' that stores model parameters introduced by this plugged measurement slice
#' @family likelihood specification functions
#' @family plug-and-play functions
#'
add_meas_BrS_ctrl_NoNest_reg_discrete_predictor_Slice_jags <- function(s, Mobs,cause_list,ppd=NULL) {
# mapping template (by `make_template` function):
patho_BrS_list <- lapply(Mobs$MBS,colnames)
template_BrS_list <-
lapply(patho_BrS_list,make_template,cause_list) # key.
# create variable names:
BrS_nm <- names(Mobs$MBS)
# index measurement slices by numbers:
JBrS_nm <- paste("JBrS",seq_along(BrS_nm),sep = "_")#
MBS_nm <- paste("MBS",seq_along(BrS_nm),sep = "_")#
mu_bs_nm <- paste("mu_bs",seq_along(BrS_nm),sep = "_")#
psiBS_nm <- paste("psiBS",seq_along(BrS_nm),sep = "_")#
FPR_stratum_id_nm <- paste("FPR_stratum_id",seq_along(BrS_nm),sep = "_")#
n_unique_FPR_level_nm <- paste("n_unique_FPR_level",seq_along(BrS_nm),sep = "_")#
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
## control BrS measurements; no subclass:
for (j in 1:",JBrS_nm[s],"){
",mu_bs_nm[s],"[i,j] <- ",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
}
"
)
} else{
plug <- paste0(
"
## control BrS measurements; no subclass (only one column):
",mu_bs_nm[s],"[i] <- ",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
"
)
}
if (!is.null(ppd) && ppd){
MBS_nm.new <- paste("MBS.new",seq_along(BrS_nm),sep = "_")#
mu_bs_nm.new <- paste("mu_bs.new",seq_along(BrS_nm),sep = "_")#
if (length(patho_BrS_list[[s]]) > 1) {
plug <- paste0(
"
## control BrS measurements; no subclass:
for (j in 1:",JBrS_nm[s],"){
",mu_bs_nm[s],"[i,j] <- ",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],j]
",MBS_nm[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
## posterior predictive distribution
",MBS_nm.new[s],"[i,j] ~ dbern(",mu_bs_nm[s],"[i,j])
}
"
)
} else{
plug <- paste0(
"
## control BrS measurements; no subclass (only one column):
",mu_bs_nm[s],"[i] <- ",psiBS_nm[s],"[",FPR_stratum_id_nm[s],"[i],1]
",MBS_nm[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
## posterior predictive distribution:
",MBS_nm.new[s],"[i] ~ dbern(",mu_bs_nm[s],"[i])
"
)
}
}
parameters <- c(psiBS_nm[s])
make_list(plug,parameters)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.