inst/doc/bpbounds.R

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
if (!requireNamespace("dplyr", quietly = TRUE)) {
  install.packages("dplyr", repos = "https://cloud.r-project.org/")
}
if (!requireNamespace("tidyr", quietly = TRUE)) {
  install.packages("tidyr", repos = "https://cloud.r-project.org/")
}

## -----------------------------------------------------------------------------
library(bpbounds)
library(dplyr)
library(tidyr)

## -----------------------------------------------------------------------------
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)
xt

## -----------------------------------------------------------------------------
p <- prop.table(xt, margin = 3)
p

## -----------------------------------------------------------------------------
bpres <- bpbounds(p)
sbp = summary(bpres)
print(sbp)

## -----------------------------------------------------------------------------
covyz = cov(tab1inddat$y, tab1inddat$z)

covxz = cov(tab1inddat$x, tab1inddat$z)

ace = covyz / covxz
ace

## -----------------------------------------------------------------------------
condprob = c(.0064, 0, .9936, 0, .0028, .001, .1972, .799)
tabp = array(condprob,
             dim = c(2, 2, 2),
             dimnames = list(
               x = c(0, 1),
               y = c(0, 1),
               z = c(0, 1)
             )) %>%
  as.table()
bpbounds(tabp)

## -----------------------------------------------------------------------------
gtab  = xtabs( ~ y + z, data = tab1inddat)
gp = prop.table(gtab, margin = 2)
gp

ttab  = xtabs( ~ x + z, data = tab1inddat)
tp = prop.table(ttab, margin = 2)
tp

bpres2 = bpbounds(p = gp, t = tp, fmt = "bivariate")
sbp2 = summary(bpres2)
print(sbp2)

## -----------------------------------------------------------------------------
mt3 <- c(.83, .05, .11, .01, .88, .06, .05, .01, .72, .05, .20, .03)
p3 = array(mt3,
           dim = c(2, 2, 3),
           dimnames = list(
             x = c(0, 1),
             y = c(0, 1),
             z = c(0, 1, 2)
           ))
p3 = as.table(p3)
bpres3 = bpbounds(p3)
sbp3 = summary(bpres3)
print(sbp3)

## -----------------------------------------------------------------------------
set.seed(2232011)
n = 10000
z = rbinom(n, 1, .5)
u = rbinom(n, 1, .5)
px = .05 + .1 * z + .1 * u
x = rbinom(n, 1, px)
p1 = .1 + .2 * z + .05 * x + .1 * u
y1 = rbinom(n, 1, p1)
p2 = .1 + .05 * z + .05 * x + .1 * u
y2 = rbinom(n, 1, p2)

tab1 = xtabs(~ x + y1 + z)
p1 <- prop.table(tab1, margin = 3)

bpres1 = bpbounds(p1)
sbp1 = summary(bpres1)
print(sbp1)

## -----------------------------------------------------------------------------
tab2 = xtabs(~ x + y2 + z)
p2 <- prop.table(tab2, margin = 3)

bpres2 = bpbounds(p2)
sbp2 = summary(bpres2)
print(sbp2)

Try the bpbounds package in your browser

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

bpbounds documentation built on May 4, 2023, 1:08 a.m.