Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- progesterone <- c(1513, 2025)
- placebo <- c(1459, 2013)
- (tab <- as.data.frame(rbind(progesterone, placebo)))
- names(tab) <- c("pass", "total")
- tab$fail <- tab$total - tab$pass
- tab$treatment <- c(1, 0)
- chisq.test(tab[c("pass", "fail")], correct = FALSE)
- # chisq.test(long.tab$treatment, long.tab$pass, correct = FALSE)
- summary(glm(cbind(pass, fail) ~ treatment, binomial, tab))
- long.tab <- data.frame(treatment = c(rep(1, tab$total[tab$treatment == 1]),
- rep(0, tab$total[tab$treatment == 0])),
- pass = c(rep(1, tab$pass[tab$treatment == 1]),
- rep(0, tab$fail[tab$treatment == 1]),
- rep(1, tab$pass[tab$treatment == 0]),
- rep(0, tab$fail[tab$treatment == 0])))
- summary(lm(pass ~ treatment, long.tab))
- summary(glm(pass ~ treatment, quasipoisson, long.tab))
- library(rstan)
- options(mc.cores = parallel::detectCores())
- rstan_options(auto_write = TRUE)
- # Place stan script in stan_scripts folder
- # Uncomment next two lines to compile stan script
- # saveRDS(stan_model(stanc_ret = stanc(file = "stan_scripts/two_group_logistic.stan")),
- # "stan_scripts/two_group_logistic.rds")
- two.grp.stan <- readRDS("stan_scripts/two_group_logistic.rds")
- dat <- list(sd_m = 4 / qnorm(.975), sd_m_diff = 2 / qnorm(.975),
- pass = rev(tab$pass), total = rev(tab$total))
- two.grp.fit <- sampling(two.grp.stan, data = dat, iter = 4e3)
- print(two.grp.fit, digits_summary = 3,
- c("m0", "m_diff", "means_prob", "odds_ratio", "prob_ratio", "prob_diff"))
- eff.df <- as.data.frame(two.grp.fit, c("odds_ratio", "prob_ratio", "prob_diff"))
- apply(eff.df, 2, function (x) quantile(x, .99))
- eff.df$odds_ratio.eff <- seq(1, 1.35, length.out = nrow(eff.df))
- eff.df$odds_ratio.p <- sapply(eff.df$odds_ratio.eff, function (x) mean(eff.df$odds_ratio > x))
- eff.df$prob_ratio.eff <- seq(1, 1.08, length.out = nrow(eff.df))
- eff.df$prob_ratio.p <- sapply(eff.df$prob_ratio.eff, function (x) mean(eff.df$prob_ratio > x))
- eff.df$prob_diff.eff <- seq(0, .06, length.out = nrow(eff.df))
- eff.df$prob_diff.p <- sapply(eff.df$prob_diff.eff, function (x) mean(eff.df$prob_diff > x))
- eff.df$post.ID <- 1:nrow(eff.df)
- names(eff.df)
- eff.df <- reshape(eff.df, direction = "long", varying = list(1:3, c(4, 6, 8), c(5, 7, 9)),
- idvar = "post.ID", times = c("Odds Ratio", "Risk Ratio", "Risk Difference"),
- v.names = c("Obs", "ES", "P"))
- eff.df$time <- factor(eff.df$time, levels = c("Odds Ratio", "Risk Ratio", "Risk Difference"))
- library(ggplot2)
- library(scales)
- theme_set(theme_bw())
- library(dplyr)
- eff.df %>%
- group_by(time) %>% mutate(eff.mean = median(Obs)) %>%
- ggplot(aes(ES, P)) + geom_line() +
- geom_vline(aes(xintercept = eff.mean), linetype = 2, alpha = .5) +
- facet_wrap(~ time, scales = "free_x") +
- scale_y_continuous(labels = percent_format(), breaks = seq(0, 1, .1)) +
- labs(x = "Plausible Effect Size", y = "Prob.(Observed Effect Size > Plausible Effect Size)",
- subtitle = "Probability progesterone had an effect exceeding varying effect thresholds") +
- theme(panel.border = element_blank(), axis.ticks = element_blank(),
- strip.background = element_blank(), panel.spacing = unit(1, "cm"))
- ggsave("progesterone_eff.pdf", height = 5, width = 6.5)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement