R/conditional_inferences.R

In lilymedina/gbiqqtools: Tools for generalized BIQQ

#' Conditional inferences
#'
#' Calculate estimands conditional on observed data (currently, for single-case process tracing) together with data realization probabilities
#' Realization probabilities are the probability of the observed data given data is sought on observed variables
#' @param model A model generated by \code{make_model}.
#' @param parameters A numeric vector. Values of parameters may be specified. By default, it is drawn from priors.
#' @param query A character vector of length 1L. Query for example 'Y[X=1]>Y[X=0]'
#' @param given A conditioning set as a character string that evaluates to a logical, for example 'Y==1'
#'
#' @export
#' @examples
#' model <- make_model("X->Y")
#' model <- set_parameters(model, type = "flat")
#' conditional_inferences(model, query = "Y[X=1]>Y[X=0]")
#'
#' # Example of posteriors given monotonic X -> M -> Y model
#' library(dplyr)
#' model <- make_model("X-> M -> Y")  %>%
#'   set_restrictions(node_restrict = list(M = "10", Y = "10")) %>%
#'   set_parameters(type = "flat")
#' conditional_inferences(model, query = "Y[X=1]>Y[X=0]", given = "Y==1")
#'
#' # Running example
#' model <- make_model("S -> C -> Y <- R <- X; X -> C -> R") %>%
#'    set_restrictions(node_restrict =
#'    list(C = "C1110", R = "R0001", Y = "Y0001"), keep = TRUE)
#' conditional_inferences(model, query = list(COE = "(Y[S=0] > Y[S=1])"),
#' given = "Y==1 & S==0")


conditional_inferences <- function(model, query, parameters=NULL,  given = NULL){

	if(is.null(parameters)) {
		if(is.null(model$parameters)) stop("parameters not provided")
		parameters <- model$parameters }

	vars <- model$variables

	# Possible data
	vals <- data.frame(perm(rep(2,length(model$variables)))) - 1
	vals[vals ==-1] <- NA
	names(vals) <- vars
	if(!is.null(given)) vals <- dplyr::filter(vals, eval(parse(text = given)))

	# Conditions
	conds <- t(apply(vals, 1, function(j) paste(vars, j, sep = "==")))
	conds[is.na(vals)] <- NA
	subsets <- apply(conds, 1, function(j) paste(j[!is.na(j)], collapse = " & "))
	subsets <- as.list(subsets)
	subsets[subsets==""] <- TRUE
	estimands <- query_model(
		model   = model,
		parameters  = parameters,
		using = "parameters",
		queries = query,
		subsets = subsets)$mean

	probs <- unlist(get_data_probs(model, vals))

	# hack to deal with fact that get_data_probs returns missing if all NAs
	p <- allNAs <- apply(vals, 1, function(j) all(is.na(j)))
	p[p] <- 1
	p[!p] <- probs

	out <- data.frame(cbind(vals, estimands, p))

	names(out) <- c(vars, "posterior", "prob")
	rownames(out) <- NULL
	data.frame(out)
}