EtzAlex

set.seed(2)
options(scipen=20) #disable scientific notation for numbers

nSim<-10 #numbber of simulated studies

library(pwr)
library(MBESS)
library(gsDesign) # The group sequential design package
library(BayesFactor)

EtzAlex

## Plots the likelihood function for the data obtained
## h = number of successes (heads), n = number of trials (flips), 
## p1 = prob of success (head) on H1, p2 = prob of success (head) on H2
## Returns the likelihood ratio for p1 over p2. The default values are the ones used in the blog post
LR <- function(h,n,p1=.5,p2=.75){
        L1 <- dbinom(h,n,p1)/dbinom(h,n,h/n) ## Likelihood for p1, standardized vs the MLE
        L2 <- dbinom(h,n,p2)/dbinom(h,n,h/n) ## Likelihood for p2, standardized vs the MLE
        Ratio <- dbinom(h,n,p1)/dbinom(h,n,p2) ## Likelihood ratio for p1 vs p2
        curve((dbinom(h,n,x)/max(dbinom(h,n,x))), xlim = c(0,1), ylab = "Likelihood",xlab = "Probability of heads",las=1,
                main = "Likelihood function for coin flips", lwd = 3)

EtzAlex

library('shiny')
# read this: http://alexanderetz.com/2015/04/15/understanding-bayes-a-look-at-the-likelihood/


shinyApp(
  ui = shinyUI(fluidPage(
    sidebarLayout(
      sidebarPanel(
        sliderInput("p1", label = "p1",
                    min = 0, max = 1, value = 0.7, step = 0.01),

EtzAlex

#The first plot with the null value and the proposed true value
x <- seq(-35,35,.001) #set up for plotting the curve
y <- dnorm(x,0,8.1) #y values for plotting curve
plot(x=x,y=y, main="Type-S and Type-M error example", xlab="Estimated effect size", 
      ylab="Relative frequency", type="l", cex.lab=1.5, axes=F, col="white")
axis(1,at=seq(-30,30,10),pos=0) #make the axis nice
axis(2,at=seq(0,.05,.01),pos=-35,las=1) #make the axis nice
lines(c(0,0),c(0,.05),col="red",lwd=3) ##Add line at null value
lines(c(2,2),c(0,.05),col="blue",lwd=3) ##Add line at population mean
points(17, .001, pch=23, bg="grey",col="black",cex=1.5) ##Add sample mean

EtzAlex

shotData<- c(1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0,
             1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1,
             0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
             1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0)

#figure 1 from blog, likelihood curve for 58/100 shots

x = seq(.001, .999, .001) ##Set up for creating the distributions
y2 = dbeta(x, 1 + 58, 1 + 42) # data for likelihood curve, plotted as the posterior from a beta(1,1)

EtzAlex

## Plots the likelihood function for the data obtained
## h = number of successes (heads), n = number of trials (flips), 
## p1 = prob of success (head) on H1, p2 = prob of success (head) on H0
#the auto plot loop is taken from http://www.r-bloggers.com/automatically-save-your-plots-to-a-folder/
#and then the pngs are combined into a gif online 


LR <- function(h,n,p1=seq(0,1,.01),p2=rep(.5,101)){
        L1 <- dbinom(h,n,p1)/dbinom(h,n,h/n) ## Likelihood for p1, standardized vs the MLE
        L2 <- dbinom(h,n,p2)/dbinom(h,n,h/n) ## Likelihood for p2, standardized vs the MLE

EtzAlex

Study#	N_orig	N_rep	r_orig	r_rep	bfRep       rep_pval	bin #	code
110	278	142	0.55	0.09	3.84E-06	0.277	1	Very strong
97	73	1486	0.38	-0.04	1.35E-03	0.154	2	Strong
8	37	31	0.56	-0.11	1.63E-02	0.540	2	Strong
4	190	268	0.23	-0.01	2.97E-02	0.920	2	Strong
65	41	131	0.43	0.01	3.06E-02	0.893	2	Strong
93	83	68	0.32	-0.14	3.12E-02	0.265	2	Strong
81	90	137	0.27	-0.10	3.24E-02	0.234	2	Strong
151	41	124	0.40	0.00	4.52E-02	0.975	2	Strong
7	99	14	0.72	0.13	5.04E-02	0.314	2	Strong

EtzAlex

## from the reproducibility project code here https://osf.io/vdnrb/
#make sure this file is in your working directory
info <- GET('https://osf.io/fgjvw/?action=download', write_disk('rpp_data.csv', overwrite = TRUE)) #downloads data file from the OSF
MASTER <- read.csv("rpp_data.csv")[1:167, ]
colnames(MASTER)[1] <- "ID" # Change first column name to ID to be able to load .csv file

studies<-MASTER$ID[!is.na(MASTER$T_r..O.) & !is.na(MASTER$T_r..R.)] ##to keep track of which studies are which
studies<-studies[-31]##remove the problem studies (46 and 139)
studies<-studies[-80]

EtzAlex

#This is the code for Alexander Etz's blog at http://wp.me/p4sgtg-SQ

#Code to make the figure and find max oracle BF for normal data
maxL <- function(mean=1.96,se=1,h0=0){
        L1 <-dnorm(mean,mean,se)
        L2 <-dnorm(mean,h0,se)
        Ratio <- L1/L2
        curve(dnorm(x,mean,se), xlim = c(-2*mean,2.5*mean), ylab = "Likelihood", xlab = "Population mean", las=1,
              main = "Likelihood function for the mean", lwd = 3)
        points(mean, L1, cex = 2, pch = 21, bg = "cyan")

EtzAlex

library(MASS)
library(ellipse)

#Helper function to make colors more transluscent
add.alpha <- function(col, alpha=1){#from https://www.r-bloggers.com/how-to-change-the-alpha-value-of-colours-in-r/
  if(missing(col))
    
    stop("Please provide a vector of colours.")
  
  apply(sapply(col, col2rgb)/255, 2,