Skip to content

Instantly share code, notes, and snippets.

View OmnibusParadox
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,
@EtzAlex
EtzAlex / MaxBF.R
Last active Jun 19, 2016
This is the code for Alexander Etz's blog at http://wp.me/p4sgtg-SQ
View MaxBF.R
#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")
View Bayes_reproducibility.R
## 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]
View Rep_BF_table.R
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
View BF_visuals.R
## 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
EtzAlex / updating.R
Created Jul 25, 2015
Updating priors via likelihood
View updating.R
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
EtzAlex / gist:7979b0e467f7593b9731
Created May 22, 2015
TypeS and TypeM blog code
View gist:7979b0e467f7593b9731
#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
View shiny_likelihoods.R
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
EtzAlex / LikelihoodFunctions
Created Apr 15, 2015
Looking at Likelihoods function
View LikelihoodFunctions
## 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)
View SampleSizesforStudies.R
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)