Skip to content

Instantly share code, notes, and snippets.

@richarddmorey

richarddmorey/X2.R

Created Jul 29, 2020
Embed
What would you like to do?
Two proportions, p value
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.