IHW on two covariates

ihw_example.R

library(ggplot2)
library(dplyr)

set.seed(1)
n <- 1e5
cov1 <- runif(n, -1, 1)
cov2 <- runif(n, -1, 1)
prop_alt <- 1/(1 + exp(-1 * (3 * cov1 + 2 * cov2 - 5)))
pvalue <- ifelse(rbinom(n, size=1, prop_alt),
                 rbeta(n, 0.25, 1),
                 runif(n))
dat <- data.frame(cov1=cut(cov1, breaks=quantile(cov1, 0:4/4), include.lowest=TRUE),
                  cov2=cut(cov2, breaks=quantile(cov2, 0:4/4), include.lowest=TRUE),
                  prop_alt,
                  pvalue)

# how many from the alt p-value dist'n
dat %>% group_by(cov1, cov2) %>%
  summarize(mean_prop_alt = mean(prop_alt))

# pvalue histograms
ggplot(dat, aes(pvalue)) +
  geom_histogram() + 
  facet_grid(cov1 ~ cov2)

# more small p-values for larger cov1, cov2
ggplot(dat, aes(cov1, -log10(pvalue))) +
  geom_violin() + facet_wrap(~cov2)

# another visual
dat %>% group_by(cov1, cov2) %>%
  summarize(mean_prop_alt = mean(prop_alt)) %>%
  ggplot(aes(cov1, cov2, fill=mean_prop_alt)) +
  geom_tile()

# how many of each bucket
dat <- dat %>% mutate(group=factor(paste(as.integer(cov1), as.integer(cov2), sep="-")))
table(dat$group)

# run IHW on the grid
library(IHW)
res <- ihw(pvalue ~ group,  data = dat,
           alpha = 0.1,
           covariate_type = "nominal")
rejections(res)
padjBH <- p.adjust(dat$pvalue, method = "BH")
sum(padjBH <= 0.1)

# plot over covariates
plot(res)

# our own plot
ihw_dat <- as.data.frame(res)
ihw_dat$cov1 <- dat$cov1
ihw_dat$cov2 <- dat$cov2
ihw_dat %>% group_by(group) %>%
  summarize(mean_wt = mean(weight), mean_small=mean(pvalue < 1e-4),
            cov1=first(cov1), cov2=first(cov2)) %>%
  ggplot(aes(mean_small, mean_wt, color=cov1, shape=cov2)) + geom_point()