3
$\begingroup$

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!

$\endgroup$
7
  • 1
    $\begingroup$ If you don't know how many people recovered (G1 and G2), I do not see how you can compare the effectiveness of the treatment regimens. The way you have it set up, it is equally possible for everyone who took medicine #1 to have died while everyone who took medicine #2 to have recovered as it is for everyone who took medicine #1 to have recovered while everyone who took medicne #2 died. // If you have the usual information in a contingency table like this, my favorite test is called the G-test. Another (more) popular one is called the Chi-squared test. $\endgroup$
    – Dave
    Oct 26, 2021 at 16:20
  • $\begingroup$ @Dave Sry my table was incorrect and I just update it. I know how many ppl recover, but I do not know how many ppl are assigned into each group. $\endgroup$
    – Ryan
    Oct 26, 2021 at 16:34
  • 3
    $\begingroup$ Agree with @Dave in that you don't have enough information. Suppose you see 10 people in G1 recover and 1000 people in G2 recover. If there are 1000 people in both G1 and G2, then G2 is clearly better (1% recovery vs. 100%). If there are 20 people in G1 and 1980 people in G2, there is virtually no difference (50% recovery vs. 50.5%). Note that in both cases, N1=10, N2=1000, and G=2000. Given only this information, you cannot compare recovery rates between groups, since you don't know the denominator for either group. $\endgroup$ Oct 26, 2021 at 16:43
  • $\begingroup$ @NuclearHoagie but it is highly unlikely that there are 20 in G1 and 1980 in G2 (prob < 0.000000001) $\endgroup$
    – Ryan
    Oct 26, 2021 at 16:46
  • $\begingroup$ I mean you can compute the probability of the realization of the table given the marginal sum, and you also know the marginal sums are Binomials. That means the probability of seeing N1 and N2 recoveries in each group is fully available..? $\endgroup$
    – Ryan
    Oct 26, 2021 at 16:53

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ is there a way to formulate a hypothesis test and calculate a p-value? $\endgroup$
    – Ryan
    Oct 26, 2021 at 18:54
  • $\begingroup$ I am not sure about that and I wouldn't trust p-values. There are a lot of criticism about p-values. Bayesian is a much better approach to reasonably measure probabilities. $\endgroup$
    – Amin Shn
    Oct 26, 2021 at 18:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.