Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(dplyr)
- library(reshape2)
- library(ggplot2)
- n <- 500
- null <- .75
- obs_n <- seq(1,n-1)
- binom_probs_low <- pbinom(obs_n, n, null)
- binom_probs_upp <- pbinom(obs_n, n, null,lower.tail = FALSE)
- binom_probs <- pmin(binom_probs_low, binom_probs_upp)
- t_probs <- rep(0,n-2)
- for (i in 1:(n-1)) {
- x = c(rep(1,i),rep(0,n-i))
- t_probs[i] <- t.test(x,mu=.75)$p.value
- }
- df <- data.frame(Sample = obs_n, reject = NA, Binomial = binom_probs, t = t_probs)
- df <- melt(df,id.vars = c("Sample","reject"),var = "Distribution")
- df$reject[df$Distribution =='Binomial'] <- ifelse(df$value[df$Distribution =='Binomial'] < .025 |df$value[df$Distribution =='Binomial'] > .975 , 1,0)
- df$reject[df$Distribution =='t'] <- ifelse(df$value[df$Distribution =='t'] < .025, 1,0)
- df <- df %>% group_by(Sample) %>%
- mutate(decision = if (all(reject==1)) { 1} else if (all(reject==0) ) { 0 } else { 2 })
- df$decision <- factor(df$decision)
- ggplot(data = df, aes(x=Sample, y = value, color = decision,shape=Distribution)) +
- geom_point() +
- coord_cartesian(xlim = c(.6*n, .9*n)) +
- ggtitle("Binomial vs. t NHST Results")
- # decision = if (all(reject==1)) { 1} else if (all(reject==0) ) { 0} else { 2}
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement