demo/evaluate_identification.R

In djhshih/orient-bias-filter: Orientation Bias Filter

library(ggplot2);
library(io);
library(devtools);

load_all("..");


n <- 200;
phi <-  anti_phred(20);
e.mean <- 30;
e.sd <- 2;

out.fname <- filename("orient-bias-model");
pdf.fname <- insert(out.fname, ext="pdf");

B <- 100;
ps <- unlist(lapply(1:B,
	function(b) {
		# theta = 0
		s <- simulate_orient_bias_read_stats(1000, 0, phi, e.mean, e.sd)
		orient_bias_adjusted_variant_test(s, phi)$p.value
	}
));

#alphas <- quantile(ps, seq(0, 1, 0.1));
alphas <- seq(0, 1, 0.1);
reject.prop <- unlist(lapply(alphas, function(p.cut) mean(ps < p.cut)));
p.df <- data.frame(alpha = alphas, empirical = reject.prop);

qdraw(
	ggplot(p.df, aes(x=alpha, y=empirical)) +
		geom_point() + theme_bw() + xlim(0, 1) + ylim(0, 1) +
		geom_abline(intercept=0, slope=1, colour="grey30") +
		ylab("Pr[p < alpha]") + xlab("alpha")
	,
	file = insert(pdf.fname, "pval-calibration")
);

B <- 100;
Ns <- seq(10, 200, by=10);
alpha <- 0.05;
type1.errors <- unlist(lapply(Ns,
	function(n) {
		ps <- unlist(lapply(1:B, function(b) {
			# theta = 0
			s <- simulate_orient_bias_read_stats(n, 0, phi, e.mean, e.sd)
			orient_bias_adjusted_variant_test(s, phi)$p.value
		}));
		mean(ps < alpha)
	}
));

type1.df <- data.frame(n = Ns, type1_error = type1.errors);

qdraw(
	ggplot(type1.df, aes(x=n, y=type1_error)) +
		geom_point() + theme_bw() + ylim(0, 0.25) +
		geom_hline(yintercept=alpha, colour="grey30") +
		xlab("Read depth") + ylab("Type I error")
	,
	file = insert(pdf.fname, "type1-error")
);


theta <- 0.01;

B <- 100;
Ns <- seq(10, 200, by=10);
alpha <- 0.05;
type2.errors <- unlist(lapply(Ns,
	function(n) {
		ps <- unlist(lapply(1:B, function(b) {
			s <- simulate_orient_bias_read_stats(n, theta, phi, e.mean, e.sd)
			orient_bias_adjusted_variant_test(s, phi)$p.value
		}));
		mean(ps > alpha)
	}
));

type2.low.signal.df <- data.frame(n = Ns, type2_error = type2.errors);

qdraw(
	ggplot(type2.low.signal.df, aes(x=n, y=type2_error)) +
		geom_point() + theme_bw() + ylim(0, 1) +
		xlab("Read depth") + ylab("Type II error")
	,
	file = insert(pdf.fname, c("type2-error", "low-signal"))
);


theta <- 0.05;

type2.errors <- unlist(lapply(Ns,
	function(n) {
		ps <- unlist(lapply(1:B, function(b) {
			s <- simulate_orient_bias_read_stats(n, theta, phi, e.mean, e.sd)
			orient_bias_adjusted_variant_test(s, phi)$p.value
		}));
		mean(ps > alpha)
	}
));

type2.df <- data.frame(n = Ns, type2_error = type2.errors);

qdraw(
	ggplot(type2.df, aes(x=n, y=type2_error)) +
		geom_point() + theme_bw() + ylim(0, 1) +
		xlab("Read depth") + ylab("Type II error")
	,
	file = insert(pdf.fname, c("type2-error", "med-signal"))
);


theta <- 0.1;

type2.errors <- unlist(lapply(Ns,
	function(n) {
		ps <- unlist(lapply(1:B, function(b) {
			s <- simulate_orient_bias_read_stats(n, theta, phi, e.mean, e.sd)
			orient_bias_adjusted_variant_test(s, phi)$p.value
		}));
		mean(ps > alpha)
	}
));

type2.high.signal.df <- data.frame(n = Ns, type2_error = type2.errors);

qdraw(
	ggplot(type2.high.signal.df, aes(x=n, y=type2_error)) +
		geom_point() + theme_bw() + ylim(0, 1) +
			xlab("Read depth") + ylab("Type II error")
	,
	file = insert(pdf.fname, c("type2-error", "high-signal"))
);