# R/bpbounds.R In bpbounds: Nonparametric Bounds for the Average Causal Effect Due to Balke and Pearl and Extensions

#### Documented in bpbounds

#' Nonparametric Bounds for the Average Causal Effect due to Balke and Pearl.
#'
#' @param p Object of class "table" containing either cell counts or conditional
#' probabilities. For trivariate data these are for the
#' phenotype/treatment-outcome association given Z, i.e. P(X, Y | Z).
#'
#' Cell counts could be generated from xtabs(~ x + y + z, data = data).
#' And then conditional probabilities obatained by calling
#' prop.table(..., margins = 3) on your object from xtabs().
#'
#' If you only know the conditional probabilities you can enter these, e.g.
#' for the Balke and Pearl Vitamin A example:
#' 
#' cp   = c(.0064, 0, .9936, 0, .0028, .001, .1972, .799)
#' tabp = as.table(array(
#'   cp,
#'   dim = c(2, 2, 2),
#'   dimnames = list(
#'     x = c(0, 1),
#'     y = c(0, 1),
#'     z = c(0, 1)
#'   )
#' ))
#' 
#' And then call bpbounds() using this object.
#'
#' For bivariate data this object contains cell conditional probabilities for
#' the outcome-instrument (Y|Z) association.
#' @param t Specified for bivariate data. Object with
#' treatment/phenotype-instrument cell counts or conditional probabilities,
#' i.e. (X|Z).
#' @param fmt A character string which sould be either "bivariate"
#' (i.e. X, Z in one dataset and Y, Z in another dataset) or
#' "trivariate" (X, Y, Z in the same dataset).
#'
#' @importFrom methods new
#' @export
#' @return List with the following elements:
#' \describe{
#' \item{fmt}{whether the data is bivariate or trivariate}
#' \item{nzcats}{2 or 3, the no. instrument categories}
#' \item{inequality}{Logical, indicating whether the IV inquality is satisfied}
#' \item{bplb}{Lower bound of ACE}
#' \item{bpub}{Upper bound of ACE}
#' \item{bplower}{Vector of lower bound probabilities}
#' \item{bpupper}{Vector of upper bound probabilities}
#' \item{p11low}{Lower bound of P(Y=1|do(X=1))}
#' \item{p11upp}{Upper bound of P(Y=1|do(X=1))}
#' \item{p10low}{Lower bound of P(Y=1|do(X=0))}
#' \item{p10upp}{Upper bound of P(Y=1|do(X=0))}
#' \item{p11lower}{Vector of probabilities for lower bound of P(Y=1|do(X=1))}
#' \item{p11upper}{Vector of probabilities for upper bound of P(Y=1|do(X=1))}
#' \item{p10lower}{Vector of probabilities for lower bound of P(Y=1|do(X=0))}
#' \item{p10upper}{Vector of probabilities for upper bound of P(Y=1|do(X=0))}
#' \item{crrlb}{Lower bound of CRR}
#' \item{crrub}{Upper bound of CRR}
#' \item{monoinequality}{Logical, indicating whether the monoticity inequality
#' is satisfied}
#' \item{monobplb}{Lower bound of ACE assuming monotonicity}
#' \item{monobpub}{Upper bound of ACE assuming monotonicity}
#' \item{monobplower}{Vector of probabilities for lower bound of ACE assuming
#' monotonicity}
#' \item{monobpupper}{Vector of probabilities for upper bound of ACE assuming
#' monotonicity}
#' \item{monop11low}{Lower bound of P(Y=1|do(X=1)) assuming monotonicity}
#' \item{monop11upp}{Upper bound of P(Y=1|do(X=1)) assuming monotonicity}
#' \item{monop10low}{Lower bound of P(Y=1|do(X=0)) assuming monotonicity}
#' \item{monop10upp}{Upper bound of P(Y=1|do(X=0)) assuming monotonicity}
#' \item{monop11lower}{Vector for corresponding bound above}
#' \item{monop11upper}{Vector for corresponding bound above}
#' \item{monop10lower}{Vector for corresponding bound above}
#' \item{monop10upper}{Vector for corresponding bound above}
#' \item{monocrrlb}{Lower bound of CRR assuming monotonicity}
#' \item{monocrrub}{Upper bound of CRR assuming monotonicity}
#' }
#'
#' @examples
#' \donttest{
#' # Vitamin A example, using cell counts
#'
#' require(tidyr)
#' require(bpbounds)
#'
#' tab1dat <- data.frame(
#'   z = c(0, 0, 1, 1, 1, 1, 0, 0),
#'   x = c(0, 0, 0, 0, 1, 1, 1, 1),
#'   y = c(0, 1, 0, 1, 0, 1, 0, 1),
#'   freq = c(74, 11514, 34, 2385, 12, 9663, 0, 0)
#' )
#'
#' tab1inddat = uncount(tab1dat, freq)
#' xt = xtabs(~ x + y + z, data = tab1inddat)
#' p = prop.table(xt, margin = 3)
#' bpres = bpbounds(p)
#' sbpres = summary(bpres)
#' print(sbpres)
#' }
#'
#' \donttest{
#' # Vitamin A example, using conditional probabilities
#'
#' require(bpbounds)
#' cp = c(.0064, 0, .9936, 0, .0028, .001, .1972, .799)
#' tabp = as.table(array(
#'   cp,
#'   dim = c(2, 2, 2),
#'   dimnames = list(
#'     x = c(0, 1),
#'     y = c(0, 1),
#'     z = c(0, 1)
#'   )
#' ))
#' bpbounds(tabp)
#' }
#'
bpbounds <- function(p, t = NULL, fmt = "trivariate") {
# Check arguments

## check fmt specified correctly
if (fmt != "bivariate" & fmt != "trivariate") {
stop('fmt argument must be either "bivariate" or "trivariate"')
}

## check t has been specified along with if fmt is bivariate
if (fmt == "bivariate" & is.null(t)) {
stop(
't, a matrix of trestment/phenotype-instrument conditional probabilities, must be specified for bivariate data'
)
}

## check p of class "table"
if (!is.table(p)) {
stop('p must be of class "table"')
}

## check length of p
if (fmt == "trivariate") {
if (length(p) != 8 & length(p) != 12) {
stop("The length of p must be 8 or 12, i.e. for a 2 or 3 category instrument.")
}
} else {
if (length(p) != 4 & length(p) != 6) {
stop("p seems to have the incorrect number of dimensions")
}
}

## check length of t
if (fmt == "bivariate") {
if (length(t) != 4 & length(t) != 6) {
stop("t seems to have the incorrect number of dimensions")
}
}

## check if p conditional probabilities
#  if (sum(p <= 1) != length(p)) {
#	  stop('p seems to be a mix of cell counts and probabilities')
#  }
if (sum(p <= 1) == length(p)) {
} else if (all(p == floor(p))) {
if (fmt == "trivariate") {
p <- prop.table(p, margin = 3)
}	else {
p = prop.table(p, margin = 2)
}
} else {
stop("All elements of p must either be conditional probabilities or cell counts.")
}

## p should now be conditional probabilities
## check they sum to approx. no. instrument categories
if (fmt == "trivariate" & length(p) == 8) {
nzcats <- 2
if (abs(sum(p) - 2) > 0.1) {
stop("The conditional probabilities add up to ",
sum(p),
}
} else if (fmt == "trivariate" & length(p) == 12) {
nzcats <- 3
if (abs(sum(p) - 3) > 0.1) {
stop("The conditional probabilities add up to ",
sum(p),
}
}

# check that t is either cell counts or conditional probabilities
if (sum(t <= 1) == length(t)) {
# conditional probabilities
} else if (all(t == floor(t))) {
# convert to conditional probabilities if cell counts
if (fmt == "bivariate") {
t = prop.table(t, margin = 2)
}
} else {
stop("All elements of t must either be conditional probabilities or cell counts.")
}

# Control flow for bivariate or trivariate data
if (fmt == "bivariate") {
if (length(p) == 4) {
cp = numeric(8)
cp[1] = p[1] # g00
cp[2] = p[2] # g10
cp[3] = p[3] # g01
cp[4] = p[4] # g11
cp[5] = t[1] # t00
cp[6] = t[2] # t10
cp[7] = t[3] # t01
cp[8] = t[4] # t11
nzcats <- 2
bp <- bpbounds_biv_x2y2z2(p = cp)
bpres = bpbounds_calc_biv_z2(g = p, t = t)
} else if (length(p) == 6) {
cp = numeric(12)
cp[1]  = p[1] # g00
cp[2]  = p[2] # g10
cp[3]  = p[3] # g01
cp[4]  = p[4] # g11
cp[5]  = p[5] # g02
cp[6]  = p[6] # g12
cp[7]  = t[1] # t00
cp[8]  = t[2] # t10
cp[9]  = t[3] # t01
cp[10] = t[4] # t11
cp[11] = t[5] # t02
cp[12] = t[6] # t12
nzcats <- 3
bp <- bpbounds_biv_x2y2z3(p = cp)
bpres = bpbounds_calc_biv_z3(g = p, t = t)
}
} else if (fmt == "trivariate") {
if (length(p) == 8) {
bp <- bpbounds_tri_x2y2z2(p = p)
bpres = bpbounds_calc_tri_z2(p)
} else if (length(p) == 12) {
bp <- bpbounds_tri_x2y2z3(p = p)
bpres = bpbounds_calc_tri_z3(p)
}
}

# return objects
retlist <- list(
"fmt" = fmt,
"nzcats" = nzcats,
"inequality" = bp$inequality, "bplb" = bp$bplow,
"bpub" = bp$bpupp, "bplower" = bp$bplower,
"bpupper" = bp$bpupper, "p11low" = bpres$p11low,
"p11upp" = bpres$p11upp, "p10low" = bpres$p10low,
"p10upp" = bpres$p10upp, "p11lower" = bpres$p11lower,
"p11upper" = bpres$p11upper, "p10lower" = bpres$p10lower,
"p10upper" = bpres$p10upper, "crrlb" = bpres$crrlb,
"crrub" = bpres$crrub, "monoinequality" = bpres$monoinequality,
"monobplb" = bpres$monobplb, "monobpub" = bpres$monobpub,
"monobplower" = bpres$monolower, "monobpupper" = bpres$monoupper,
"monop11low" = bpres$monop11lb, "monop11upp" = bpres$monop11ub,
"monop10low" = bpres$monop10lb, "monop10upp" = bpres$monop10ub,
"monop11lower" = bpres$monop11lower, "monop11upper" = bpres$monop11upper,
"monop10lower" = bpres$monop10lower, "monop10upper" = bpres$monop10upper,
"monocrrlb" = bpres$monocrrlb, "monocrrub" = bpres$monocrrub
)
class(retlist) <- "bpbounds"
return(retlist)
}

#' @export
summary.bpbounds <- function(object, ...) {
ans = list()
bp = object
ans$fmt = bp$fmt
ans$nzcats = bp$nzcats
ans$inequality = bp$inequality
ans$monoinequality = bp$monoinequality

boundscolnames = c("Causal parameter", "Lower bound", "Upper bound")
causalparameternames = c("ACE", "P(Y|do(X=0))", "P(Y|do(X=1))", "CRR")

if (bp$inequality) { lb = c(bp$bplb, bp$p10low, bp$p11low, bp$crrlb) ub = c(bp$bpub, bp$p10upp, bp$p11upp, bp$crrub) ans$bounds = data.frame(causalparameternames, lb, ub)
colnames(ans$bounds) = boundscolnames if (bp$monoinequality) {
monolb = c(bp$monobplb, bp$monop10low, bp$monop11low, bp$monocrrlb)
monoub = c(bp$monobpub, bp$monop10upp, bp$monop11upp, bp$monocrrub)
ans$monobounds = data.frame(causalparameternames, monolb, monoub) colnames(ans$monobounds) = boundscolnames
}
}

class(ans) <- "summary.bpbounds"
ans
}

#' @export
print.summary.bpbounds <-
function(x, digits = getOption("digits"), ...) {
cat("\n")
cat("Data:                    ", x$fmt, "\n", sep = "") cat("Instrument categories: ", x$nzcats, "\n\n", sep = "")
cat("Instrumental inequality:", x$inequality, "\n") if (x$inequality) {
print(x$bounds, digits = digits, row.names = FALSE, ...) } cat("\nMonotonicity inequality:", x$monoinequality, "\n")
if (x$monoinequality) { print(x$monobounds,
digits = digits,
row.names = FALSE,
...)
}
cat("\n")
invisible(x)
}

#' @export
print.bpbounds <- function(x, digits = getOption("digits"), ...) {
print(summary(x), digits = digits, ...)
}


## Try the bpbounds package in your browser

Any scripts or data that you put into this service are public.

bpbounds documentation built on Jan. 21, 2020, 5:07 p.m.