Two group logistic model

progesterone.R

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)

two_group_logistic.stan

data {
  real<lower = 0> sd_m; // 3 / qnorm(.975) is reasonable default
  real<lower = 0> sd_m_diff; // (2, 3, 10) / qnorm(.975) for different levels of belief
  int pass[2];
  int total[2];
}
parameters {
  real m0;
  real m_diff;
}
transformed parameters {
  vector[2] means_logit;
  means_logit[1] = m0;
  means_logit[2] = m0 + m_diff;
}
model {
  m0 ~ normal(0, sd_m);
  m_diff ~ normal(0, sd_m_diff);

  pass ~ binomial_logit(total, means_logit);
}
generated quantities {
  vector[2] log_lik;
  vector[2] yhat;
  real odds_ratio = exp(m_diff);
  vector[2] means_prob = inv_logit(means_logit);
  real prob_ratio = means_prob[2] / means_prob[1];
  real prob_diff = means_prob[2] - means_prob[1];

  for (i in 1:2) {
    log_lik[i] = binomial_logit_lpmf(pass[i] | total[i], means_logit[i]);
    yhat[i] = binomial_rng(total[i], means_prob[i]);
  }
}