## ----message=FALSE, warning=FALSE, include=FALSE-------------------------
library(knitr)
options(knitr.kable.NA = '')
knitr::opts_chunk$set(comment=">")
options(digits=2)
## ----message=FALSE, warning=FALSE----------------------------------------
library(ggplot2)
library(dplyr)
library(tidyr)
df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayesSim_study1.csv") %>%
mutate(bayesfactor_0707 = log(bayesfactor_0707),
bayesfactor_1 = log(bayesfactor_1))
## ---- message=FALSE, warning=FALSE, fig.height=25, fig.width=15----------
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_90, ROPE_95, ROPE_full, bayesfactor_0707, bayesfactor_1) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_value", "p_direction", "p_MAP", "p_ROPE", "ROPE_90", "ROPE_95", "ROPE_full", "bayesfactor_0707", "bayesfactor_1"))) %>%
mutate(temp = as.factor(cut(error, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(error_group = round(mean(error), 1)) %>%
ungroup() %>%
ggplot(aes(x = error_group, y = value, fill = index, colour=true_effect, group = interaction(index, true_effect, error_group))) +
# geom_jitter(shape=16, alpha=0.02) +
geom_boxplot(outlier.shape = NA) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_text(face="bold")) +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A", name="Effect")) +
scale_fill_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_90" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "bayesfactor_0707"="#9C27B0", "bayesfactor_1"="#673AB7"), guide=FALSE) +
ylab("Index Value\n") +
xlab("\nNoise")
## ---- message=FALSE, warning=FALSE, fig.height=25, fig.width=15----------
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_90, ROPE_95, ROPE_full, bayesfactor_0707, bayesfactor_1) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_value", "p_direction", "p_MAP", "p_ROPE", "ROPE_90", "ROPE_95", "ROPE_full", "bayesfactor_0707", "bayesfactor_1"))) %>%
mutate(temp = as.factor(cut(sample_size, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = round(mean(sample_size))) %>%
ungroup() %>%
ggplot(aes(x = size_group, y = value, fill = index, colour=true_effect, group = interaction(index, true_effect, size_group))) +
# geom_jitter(shape=16, alpha=0.02) +
geom_boxplot(outlier.shape = NA) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_text(face="bold")) +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A", name="Effect")) +
scale_fill_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_90" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "bayesfactor_0707"="#9C27B0", "bayesfactor_1"="#673AB7"), guide=FALSE) +
ylab("Index Value\n") +
xlab("\nSample Size")
## ---- message=FALSE, warning=FALSE, fig.height=25, fig.width=15----------
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_90, ROPE_95, ROPE_full, bayesfactor_0707, bayesfactor_1) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type, -p_value) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_direction", "p_MAP", "p_ROPE", "ROPE_90", "ROPE_95", "ROPE_full", "bayesfactor_0707", "bayesfactor_1"))) %>%
mutate(temp = as.factor(cut(sample_size, 3, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = as.factor(round(mean(sample_size)))) %>%
ungroup() %>%
ggplot(aes(x = p_value, y = value, color = true_effect, shape=size_group)) +
geom_point(alpha=0.025, stroke = 0, shape=16) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_text(face="bold")) +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A"), name="Effect") +
guides(colour = guide_legend(override.aes = list(alpha = 1)),
shape = guide_legend(override.aes = list(alpha = 1), title="Sample Size"))
## ---- message=FALSE, warning=FALSE, fig.height=15, fig.width=10----------
df$sig_1 <- factor(ifelse(df$p_value >= .1, "n.s.", "-"), levels=c("n.s.", "-"))
df$sig_05 <- factor(ifelse(df$p_value >= .05, "n.s.", "*"), levels=c("n.s.", "*"))
df$sig_01 <- factor(ifelse(df$p_value >= .01, "n.s.", "**"), levels=c("n.s.", "**"))
df$sig_001 <- factor(ifelse(df$p_value >= .001, "n.s.", "***"), levels=c("n.s.", "***"))
get_data <- function(predictor, outcome, lbound=0, ubound=0.3){
fit <- glm(paste(outcome, "~ outcome_type * ", predictor), data=df, family = "binomial")
# data <- data.frame(x=rep(1:100, 2))
data <- data.frame(outcome_type=rep(c("linear", "binary"), each=100))
data[predictor] <- rep(seq(lbound, ubound, length.out = 100), 2)
data$index <- predictor
predict_fit <- predict(fit, newdata=data, type="response", se.fit = TRUE)
data[outcome] <- predict_fit$fit
data$CI_lower <- predict_fit$fit - (qnorm(0.99) * predict_fit$se.fit)
data$CI_upper <- predict_fit$fit + (qnorm(0.99) * predict_fit$se.fit)
data <- select_(data, "value"=predictor, outcome, "outcome_type", "index", "CI_lower", "CI_upper")
return(data)
}
rbind(
rbind(
get_data(predictor="p_direction", outcome="sig_001", lbound=99.5, ubound=100),
get_data(predictor="p_MAP", outcome="sig_001", lbound=0, ubound=0.01),
get_data(predictor="p_ROPE", outcome="sig_001", lbound=97, ubound=100),
get_data(predictor="ROPE_90", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="ROPE_95", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="ROPE_full", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="bayesfactor_0707", outcome="sig_001", lbound=0, ubound=10),
get_data(predictor="bayesfactor_1", outcome="sig_001", lbound=0, ubound=10)
) %>%
rename("sig"=sig_001) %>%
mutate(threshold="p < .001"),
rbind(
get_data(predictor="p_direction", outcome="sig_01", lbound=98, ubound=100),
get_data(predictor="p_MAP", outcome="sig_01", lbound=0, ubound=0.1),
get_data(predictor="p_ROPE", outcome="sig_01", lbound=85, ubound=100),
get_data(predictor="ROPE_90", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="ROPE_95", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="ROPE_full", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="bayesfactor_0707", outcome="sig_01", lbound=0, ubound=5),
get_data(predictor="bayesfactor_1", outcome="sig_01", lbound=0, ubound=5)
) %>%
rename("sig"=sig_01) %>%
mutate(threshold="p < .01"),
rbind(
get_data(predictor="p_direction", outcome="sig_05", lbound=95, ubound=100),
get_data(predictor="p_MAP", outcome="sig_05", lbound=0, ubound=0.3),
get_data(predictor="p_ROPE", outcome="sig_05", lbound=50, ubound=100),
get_data(predictor="ROPE_90", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="ROPE_95", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="ROPE_full", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="bayesfactor_0707", outcome="sig_05", lbound=0, ubound=2),
get_data(predictor="bayesfactor_1", outcome="sig_05", lbound=0, ubound=2)
) %>%
rename("sig"=sig_05) %>%
mutate(threshold="p < .05"),
rbind(
get_data(predictor="p_direction", outcome="sig_1", lbound=90, ubound=100),
get_data(predictor="p_MAP", outcome="sig_1", lbound=0, ubound=0.5),
get_data(predictor="p_ROPE", outcome="sig_1", lbound=25, ubound=100),
get_data(predictor="ROPE_90", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="ROPE_95", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="ROPE_full", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="bayesfactor_0707", outcome="sig_1", lbound=0, ubound=1),
get_data(predictor="bayesfactor_1", outcome="sig_1", lbound=0, ubound=1)
) %>%
rename("sig"=sig_1) %>%
mutate(threshold="p < .1")
) %>%
mutate(index = as.factor(index)) %>%
ggplot(aes(x=value, y=sig)) +
geom_ribbon(aes(ymin=CI_lower, ymax=CI_upper), alpha=0.1) +
geom_line(aes(color=index), size=1) +
facet_wrap(~ index * threshold * outcome_type, scales = "free", ncol=8) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_text(face="bold")) +
scale_color_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_90" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "bayesfactor_0707"="#9C27B0", "bayesfactor_1"="#673AB7"), guide=FALSE) +
ylab("Probability of being significant\n") +
xlab("\nIndex Value")
## ---- message=FALSE, warning=FALSE, fig.height=15, fig.width=10----------
df$equivalence_95 <- factor(ifelse(df$ROPE_95 == 0, "significant", "n.s."), levels=c("n.s.", "significant"))
df$equivalence_90 <- factor(ifelse(df$ROPE_90 == 0, "significant", "n.s."), levels=c("n.s.", "significant"))
rbind(
rbind(
get_data(predictor="p_direction", outcome="equivalence_95", lbound=99.5, ubound=100),
get_data(predictor="p_MAP", outcome="equivalence_95", lbound=0, ubound=0.02),
get_data(predictor="p_ROPE", outcome="equivalence_95", lbound=0, ubound=100),
get_data(predictor="ROPE_90", outcome="equivalence_95", lbound=0, ubound=0.5),
get_data(predictor="ROPE_full", outcome="equivalence_95", lbound=0, ubound=0.5),
get_data(predictor="bayesfactor_0707", outcome="equivalence_95", lbound=0, ubound=3),
get_data(predictor="bayesfactor_1", outcome="equivalence_95", lbound=0, ubound=3)
) %>%
rename("equivalence"=equivalence_95) %>%
mutate(level="95 HDI"),
rbind(
get_data(predictor="p_direction", outcome="equivalence_90", lbound=99.5, ubound=100),
get_data(predictor="p_MAP", outcome="equivalence_90", lbound=0, ubound=0.02),
get_data(predictor="p_ROPE", outcome="equivalence_90", lbound=0, ubound=100),
get_data(predictor="ROPE_95", outcome="equivalence_90", lbound=0, ubound=0.5),
get_data(predictor="ROPE_full", outcome="equivalence_90", lbound=0, ubound=0.5),
get_data(predictor="bayesfactor_0707", outcome="equivalence_90", lbound=0, ubound=3),
get_data(predictor="bayesfactor_1", outcome="equivalence_90", lbound=0, ubound=3)
) %>%
rename("equivalence"=equivalence_90) %>%
mutate(level="90 HDI")
) %>%
ggplot(aes(x=value, y=equivalence)) +
geom_ribbon(aes(ymin=CI_lower, ymax=CI_upper), alpha=0.1) +
geom_line(aes(color=index), size=1) +
facet_wrap(~ index * level * outcome_type, scales = "free", nrow=7) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_text(face="bold")) +
scale_color_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_90" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "bayesfactor_0707"="#9C27B0", "bayesfactor_1"="#673AB7"), guide=FALSE) +
ylab("Probability of rejecting H0 with the equivalence test\n") +
xlab("\nIndex Value")
## ---- message=FALSE, warning=FALSE, fig.height=6, fig.width=10.08--------
df %>%
mutate(true_effect = as.factor(true_effect)) %>%
ggplot(aes(x=p_direction, y=ROPE_full, color=true_effect)) +
geom_point(alpha=0.025, stroke = 0, shape=16) +
facet_wrap(~ outcome_type, scales = "free", ncol=2) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_text(face="bold")) +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A"), name="Effect") +
ylab("ROPE (full)\n") +
xlab("\nProbability of Direction")