Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(discreteRV)
- placebo.rate <- 0.5
- mmm.rate <- 0.3
- mmm.power <- power.prop.test(p1 = placebo.rate, p2 = mmm.rate, power = 0.8, alternative = "one.sided")
- n <- as.integer(ceiling(mmm.power$n))
- patients <- seq(from = 0, to = n, by = 1)
- placebo_distribution <- dbinom(patients, size = n, prob = placebo.rate)
- mmm_distribution <- dbinom(patients, size = n, prob = mmm.rate)
- get_pmf <- function(p1, p2) {
- X1 <- RV(patients,p1, fractions = F)
- X2 <- RV(patients,p2, fractions = F)
- pmf <- joint(X1, X2, fractions = F)
- return(pmf)
- }
- extract <- function(string) {
- ints <- unlist(strsplit(string,","))
- x1 <- as.integer(ints[1])
- x2 <- as.integer(ints[2])
- return(x1-x2)
- }
- diff_prob <- function(pmf) {
- diff <- unname(sapply(outcomes(pmf),FUN = extract)/n)
- probabilities <- unname(probs(pmf))
- df <- data.frame(diff,probabilities)
- df <- aggregate(. ~ diff, data = df, FUN = sum)
- return(df)
- }
- most_likely_rate <- function(x) {
- x[which(x$probabilities == max(x$probabilities)),]$diff
- }
- mmm_rate_diffs <- diff_prob(get_pmf(mmm_distribution,placebo_distribution))
- placebo_rate_diffs <- diff_prob(get_pmf(placebo_distribution,placebo_distribution))
- plot(mmm_rate_diffs$diff,mmm_rate_diffs$probabilities * 100, type = "l", lty = 2, xlab = "Rate difference", ylab = "# of trials per 100", main = paste("Trials with",n,"patients per treatment arm",sep = " "))
- lines(placebo_rate_diffs$diff, placebo_rate_diffs$probabilities * 100, lty = 1, xaxs = "i")
- abline(v = c(most_likely_rate(placebo_rate_diffs), most_likely_rate(mmm_rate_diffs)), lty = c(1,2))
- legend("topleft", legend = c("Alternative hypothesis", "Null hypothesis"), lty = c(2,1))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement