Two proportions, p value

X2.R

N = c(16, 16)
y = c(3, 9)
actual_result = prop.test(y, N, correct = FALSE)


compute_p_value = function(p0){
  ## Probability of all outcomes,
  ## assuming independence
  pr_X2s = 
    outer( dbinom(0:N[1], N[1], p0), 
           dbinom(0:N[2], N[2], p0), "*")
  
  ## Compute X2 stats for all outcomes
  X2s = outer( 0:N[1],
               0:N[2],
               Vectorize(function(y1,y2){
                 X2 = suppressWarnings(
                   prop.test(c(y1,y2), N, correct = FALSE)$statistic
                 )
                 ifelse(is.nan(X2),0,X2) ## (0,0) and (N,N) count toward null
               }, c("y1","y2")))
  
  ## Probability of as or more extreme result
  more_extreme = X2s >= actual_result$statistic
  sum(more_extreme * pr_X2s)
}

## Guess worst case null
p0 = sum(y)/sum(N)
compute_p_value(p0)

## Find worst case through optimization
optimise(compute_p_value, interval = c(0,1),
         maximum = TRUE)

## Compare to approximation
prop.test(y, N, correct = FALSE)
prop.test(y, N, correct = TRUE)