inst/doc/episensr.R

In episensr: Basic Sensitivity Analysis of Epidemiological Results

## ----selection----------------------------------------------------------------
library(episensr)

stang <- selection(matrix(c(136, 107, 297, 165),
						  dimnames = list(c("UM+", "UM-"), c("Mobile+", "Mobile-")),
						  nrow = 2, byrow = TRUE),
				   bias_parms = c(.94, .85, .64, .25))
stang

## ----misclass-----------------------------------------------------------------
misclass(matrix(c(126, 92, 71, 224),
				dimnames = list(c("Case", "Control"),
								c("Smoking +", "Smoking - ")),
				nrow = 2, byrow = TRUE),
		 type = "exposure",
		 bias_parms = c(0.94, 0.94, 0.97, 0.97))

## ----warning=FALSE,fig.cap='Contour plot of point estimates of corrected odds ratio (OR)',fig.width=6,fig.height=4----
dat <- expand.grid(Se = seq(0.582, 1, 0.02),
				   Sp = seq(0.762, 1, 0.02))

dat$OR_c <- apply(dat, 1,
				  function(x) misclass(matrix(c(126, 92, 71, 224),
											  nrow = 2, byrow = TRUE),
									   type = "exposure",
									   bias_parms = c(x[1], x[1], x[2], x[2]))$adj_measures[2, 1])
dat$OR_c <- round(dat$OR_c, 2)
dat$lab <- cut(dat$OR_c, breaks = c(4.32, 6.96, 8.56, 12.79, 23.41, 432.1, Inf), right = FALSE,
				 labels = c("", "6.96", "8.56", "12.79", "23.41", "432.1"))

library(ggplot2)
library(directlabels)
ggplot(dat, aes(x = Se, y = Sp, z = OR_c)) +
	geom_contour(aes(colour = after_stat(level)), breaks = c(4.32, 6.96, 8.56, 12.79, 23.41, 432.1)) +
	annotate("text", x = 1, y = 1, label = "4.32", size = 3) +
	scale_x_continuous(breaks = seq(0.5, 1, .1), expand = c(0, 0)) +
	scale_y_continuous(breaks = seq(0.5, 1, .1), expand = c(0, 0)) +
	coord_fixed(ylim = c(0.5, 1.025), xlim = c(0.5, 1.025)) +
	scale_colour_gradient(guide = "none") +
	xlab("Sensitivity") +
	ylab("Specificity") +
	ggtitle("Estimates of Corrected OR") +
	geom_dl(aes(label = lab), method = list("last.points", hjust = 1, vjust = -.1, cex = .7))
dat$lab <- NULL

## ----warning=FALSE,fig.cap='Contour plot of 95% lower confidence limit of corrected OR',warning=FALSE,fig.width=6,fig.height=4----
dat$ORc_lci <- apply(dat, 1,
					 function(x) misclass(matrix(c(126, 92, 71, 224),
												 nrow = 2, byrow = TRUE),
										  type = "exposure",
										  bias_parms = c(x[1], x[1], x[2], x[2]))$adj_measures[2, 2])
dat$ORc_lci <- round(dat$ORc_lci, 2)

p3 <- ggplot(dat, aes(x = Se, y = Sp, z = ORc_lci)) +
	geom_contour(aes(colour = ..level..),
				 breaks = c(4.05, 4.64, 5.68, 7.00, 9.60)) +
	annotate("text", x = 1, y = 1, label = "2.96", size = 3) +
	scale_fill_gradient(limits = range(dat$ORc_lci), high = 'red', low = 'green') +
	scale_x_continuous(breaks = seq(0.5, 1, .1), expand = c(0,0)) +
	scale_y_continuous(breaks = seq(0.5, 1, .1), expand = c(0,0)) +
	coord_fixed(ylim = c(0.5, 1.025), xlim = c(0.5, 1.025)) +
	scale_colour_gradient(guide = 'none') +
	xlab("Sensitivity") +
	ylab("Specificity") +
	ggtitle("95% LCI of Corrected OR")
direct.label(p3, list("far.from.others.borders", "calc.boxes", "enlarge.box",
						   hjust = 1, vjust = -.5, box.color = NA, cex = .6,
						   fill = "transparent", "draw.rects"))

## ----misclass_pv--------------------------------------------------------------
misclass(matrix(c(599, 4978, 31175, 391851),
				dimnames = list(c("Preterm", "Term"), c("Underweight", "Normal weight")),
				nrow = 2, byrow = TRUE),
		 type = "exposure_pv",
		 bias_parms = c(0.65, 0.74, 1, 0.98))

## ----cov_misclass-------------------------------------------------------------
misclass_cov(array(c(1319, 38054, 5641, 405546, 565, 3583,
					 781, 21958, 754, 34471, 4860, 383588),
				   dimnames = list(c("Twins+", "Twins-"),
								   c("Folic acid+", "Folic acid-"),
								   c("Total", "IVF+", "IVF-")),
				   dim = c(2, 2, 3)),
			 bias_parms = c(.6, .6, .95, .95))

## ----confounders--------------------------------------------------------------
tyndall <- confounders(matrix(c(105, 85, 527, 93),
							  dimnames = list(c("HIV+", "HIV-"), c("Circ+", "Circ-")),
							  nrow = 2, byrow = TRUE),
					   type = "RR",
					   bias_parms = c(.63, .8, .05))
tyndall

## ----confounders_add----------------------------------------------------------
## Confounder +
tyndall$cfder_data

## Confounder -
tyndall$nocfder_data

## ----miettinen1---------------------------------------------------------------
rash <- matrix(c(15, 94, 52, 1163),
			   dimnames = list(c("Rash +", "Rash -"),
							   c("Allopurinol +", "Allopurinol -")),
			   nrow = 2, byrow = TRUE)

## ----miettinen1-tab, echo=FALSE-----------------------------------------------
knitr::kable(rash)

## ----miettinen1-conf----------------------------------------------------------
## Outcome by confounders
rash_conf <- matrix(c(36, 58, 645, 518),
					dimnames = list(c("Rash +", "Rash -"),
									c("Males", "Females")),
					nrow = 2, byrow = TRUE)
## By confounders: among males
rash_males <- matrix(c(5, 36, 33, 645),
					 dimnames = list(c("Rash +", "Rash -"),
									 c("Allopurinol +", "Allopurinol -")),
					 nrow = 2, byrow = TRUE)
## By confounders: among females
rash_females <- matrix(c(10, 58, 19, 518),
					   dimnames = list(c("Rash +", "Rash -"),
									   c("Allopurinol +", "Allopurinol -")),
					   nrow = 2, byrow = TRUE)

## ----miettinen1-conftab, echo=FALSE-------------------------------------------
knitr::kable(rash_conf, caption = "Outcome by confounders")
knitr::kable(rash_males, caption = "Among males")
knitr::kable(rash_females, caption = "Among females")

## ----miettinen1_parms---------------------------------------------------------
# RR between confounder and outcome among non-exposed
(RR_no <- (36/(36+645))/(58/(58+518)))
# prevalence of confounder among exposed
(p1 <- (5+33)/(15+52))
# prevalence of confounder among unexposed
(p2 <- (36+645)/(94+1163))

## ----miettinen1_res-----------------------------------------------------------
rash %>% confounders(., type = "RR", bias_parms = c(RR_no, p1, p2))

## ----miettinen2---------------------------------------------------------------
thromb <- matrix(c(12, 30, 53, 347),
				 dimnames = list(c("Thrombosis +", "Thrombosis -"),
								 c("OC +", "OC -")),
				 nrow = 2, byrow = TRUE)

## ----miettinen2-tab, echo=FALSE-----------------------------------------------
knitr::kable(thromb)

## ----miettinen2-conf----------------------------------------------------------
## Outcome by confounders
thromb_conf <- matrix(c(18, 12, 158, 189),
					  dimnames = list(c("Thrombosis +", "Thrombosis -"),
									  c("20-29 years old", "30-39 years old")),
					  nrow = 2, byrow = TRUE)
## By confounders: among 20-29 years old
thromb_younger <- matrix(c(10, 18, 39, 158),
					   dimnames = list(c("Thrombosis +", "Thrombosis -"),
									   c("OC +", "OC -")),
					   nrow = 2, byrow = TRUE)
## By confounders: among 30-39 years old
thromb_older <- matrix(c(2, 12, 14, 189),
					 dimnames = list(c("Thrombosis +", "Thrombosis -"),
									 c("OC +", "OC -")),
					 nrow = 2, byrow = TRUE)

## ----miettinen2-conftab, echo=FALSE-------------------------------------------
knitr::kable(thromb_conf, caption = "Outcome by confounders")
knitr::kable(thromb_younger, caption = "Among 20-29 years old")
knitr::kable(thromb_older, caption = "Among 30-39 years old")

## ----miettinen2_parms---------------------------------------------------------
# OR between confounder and outcome among non-exposed
(OR_no <- (18/12)/(158/189))
# prevalence of confounder among exposed
(p1 <- (10+39)/(12+53))
# prevalence of confounder among unexposed
(p2 <- (18+158)/(30+347))

## ----miettinen2_res-----------------------------------------------------------
thromb %>% confounders(., type = "OR", bias_parms = c(OR_no, p1, p2))

## ----boot---------------------------------------------------------------------
library(aplore3)  # to get ICU data
data(icu)

## First run the model
misclass_eval <- misclass(icu$sta, icu$inf,
						  type = "exposure",
						  bias_parms = c(.75, .85, .9, .95))
misclass_eval

## Then bootstrap it
set.seed(123)
misclass_boot <- boot.bias(misclass_eval, R = 10000)
misclass_boot

## ----boot_fig, fig.cap = "Bootstrap replicates and confidence interval.", warning=F----
plot(misclass_boot, "rr")