https://github.com/cran/bayestestR
Raw File
Tip revision: 9985109256c08d654edb46adb9cb20c913fb1888 authored by Dominique Makowski on 29 May 2019, 14:10:02 UTC
version 0.2.0
Tip revision: 9985109
indicesExistenceComparison.R
## ----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")

back to top