2x2 contingency table: what test should be used?

Question

Assume the following experiment: $G$ people with certain disease are randomly assigned into 2 groups with equal probability. People in group $i$ will received medicine $m_i, i=1,2$. After a period of time, it is observed that in group $i$, $N_i$($i=1,2$) people have recovered. Note that we do not know how many people were assigned to each group.

The numbers can be summarized into the following 2x2 table

	recover	not recover	total
group1	N1	G1 - N1	G1
group2	N2	G2 - N2	G2
total	N	G - N	G

Here we assume that only the values N1, N2 and G are known. We don't have the full table. We are interested in which medicine has a higher recovery rate. How can we test it based on the information? I am aware of the exact test by Fisher and Barnard, but I am not sure if they could be directly applied in this simple case.

As an example, assume N1=40, N2=20, G=100. Since the patients are assigned into each group with equal probability, we are likely to see a distribution of ~50/50:

	recover	not recover	total
group1	40	10	50
group2	20	30	50
total	60	40	100

In such scenario it's intuitively clear(Barnard's test has p-value of 0.000039) that m1 is more effective that m2.

PS: If I simply assume that G1( and thus G2) are given, perform Barnard's test to get a p-value $p=p(G_1)$, and aggregate these p-values through $$ p = \sum p(G_1) P(G_1=g) $$ where the summation is over all possible value of $g$, and $G_1 \sim Bin(G,0.5)$, would I get something that makes sense?

Sorry if this is a noob question as I am not familiar in this area. Thanks in advance!

Answer 1

Let's summarise what we know and what we don't know in your question. We know the total number of participants $G$ and the probability of assigning people into two groups being 0.5. But we don't know the exact number of people assigned to each group, we only know how many people recovered within each group. So we should estimate both the probability of recovery and number of people in each group. I would do it this way:

$y_{i} \sim Binomial(N_{i}, p_{i})$

where $y_{i}$ is the number of recovered patients in group i and $p_{i}$ is the probability of recovery in that group and $N_{i}$ is the number of people in group i. So lets estimate both N and P this way:

$N_{1} \sim Binomial(G, 0.5)$ & $N_{2} = G - N_{1}$

We can easily fit this model in R using JAGS and Bayesian framework to estimate the probability of recovery for each group that we are interested in. Below I generated some data according to your model assumption and fit the model I just described:

library(R2jags) # Package needed for Bayesian inference
library(tidyverse)
Some R code to simulate data
G <- 200 #total sample size
T1 <- rbinom(1, G,0.5) # group 1 sample size randomly assigned
T2 = G - T1  # group 2 sample size randomly assigned
set.seed(123)
p1 <- 0.3 # Probability of recovery by medicine 1
p2 <- 0.2 # Probability of recovery by medicine 2
y1 <- rbinom(T1, 1, p1)
y2 <- rbinom(T2, 1, p2)
df1 <- data.frame(y = y1, medicine = 1)
df2 <- data.frame(y = y2, medicine = 2)
df_merged <- rbind(df1, df2) %>% group_by(medicine) %>% summarise(y = sum(y))
Jags code to fit the model to the simulated data
model_code <- "
model
{
Likelihood
for (i in 1:I) {
    y[i] ~ dbin(p[med[i]], N[med[i]])
  }
N[1] ~ dbin(0.5,G)
  N[2] <- G - N[1]
Priors
for(i in 1:I){
  p[i] ~ dunif(0,1)
  }
}
"
Set up the data
model_data <- list(I = nrow(df_merged), 
                   y = df_merged$y, 
               med = df_merged$medicine,
                   G = G)
Choose the parameters to watch
model_parameters <- c("p")
Run the model
model_run <- jags(
  data = model_data,
  parameters.to.save = model_parameters,
  model.file = textConnection(model_code)
)
print(model_run)

And here are the results:

> print(model_run)
Inference for Bugs model at "7", fit using jags,
 3 chains, each with 2000 iterations (first 1000 discarded)
 n.sims = 3000 iterations saved
         mu.vect sd.vect  2.5%   25%    50%    75%  97.5%  Rhat n.eff
p[1]       0.305   0.050 0.212 0.270  0.303  0.338  0.409 1.002  1800
p[2]       0.186   0.041 0.115 0.157  0.184  0.212  0.271 1.001  3000

As you can see the model correctly estimated the values of $p_{1}$ = 0.3 (0.21 - 0.4) and $p_{2}$ = 0.19 (0.11 - 0.27) suggesting that patients in group one had higher chance of recovery. Note that the confidence interval for each probability is also given.