Created
May 16, 2019 12:45
-
-
Save ljcolling/9e7483de3ab64cea7c7759735de7107d to your computer and use it in GitHub Desktop.
R gist for calculating Bayes factors with arbitrary location parameter
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| BayesFactor::ttest.tstat(t = sqrt(N) * mean(x - 1) / sd(x), n1 = 10, simple = T) | |
| #rnorm(10, mean = 1, sd = .5) -> x | |
| N=10 | |
| BayesFactor::ttestBF(x = x, mu = 1) | |
| my.es=mean(x) / sd(x) | |
| nullvalue = 0 | |
| ############################ | |
| #Specify Alternative (up to constant of propertionality) | |
| #Change this section to change alternative | |
| lo=-Inf #lower bound of support | |
| hi=Inf #upper bound of support | |
| altDens=function(delta){dcauchy(delta,location = nullvalue, scale = 0.707)} | |
| ########################### | |
| #Normalize alternative density in case user does not, | |
| K=1/integrate(altDens,lower=lo,upper=hi)$value | |
| f=function(delta) K*altDens(delta) | |
| delta=seq(-10,10,.01) | |
| #Plot Alternative as a density and Null as a point arrow | |
| maxAlt=max(f(delta)) | |
| plot(delta,f(delta),typ='n',xlab="Effect Size Parameter Delta",ylab="Density",ylim=c(0,1.4*maxAlt),main="Models") | |
| arrows(0,0,0,1.3*maxAlt,col='darkblue',lwd=2) | |
| lines(delta,f(delta),col='darkgreen',lwd=2) | |
| legend("topright",legend=c("Null","Alternative"),col=c('darkblue','darkgreen'),lwd=2) | |
| #Prediction Function Under Null | |
| nullPredF=function(obs,N) dt(sqrt(N)*obs,N-1,ncp = nullvalue) | |
| #Prediction Function Under the Alternative | |
| altPredIntegrand=function(delta,obs,N) | |
| dt(sqrt(N)*obs,N-1,sqrt(N)*delta)*f(delta) | |
| altPredF=function(obs,N) | |
| integrate(altPredIntegrand,lower=lo,upper=hi,obs=obs,N=N)$value | |
| obs=delta | |
| I=length(obs) | |
| nullPred=nullPredF(obs,N) | |
| altPred=1:I | |
| for (i in 1:I) altPred[i]=altPredF(obs[i],N) | |
| #Plot The Predictions | |
| top=max(altPred,nullPred) | |
| plot(type='l',obs,nullPred,ylim=c(0,top),xlab="Observed Effect Size",ylab="Density",main="Predictions",col='darkblue',lwd=2) | |
| lines(obs,altPred,col='darkgreen',lwd=2) | |
| legend("topright",legend=c("Null","Alternative"),col=c('darkblue','darkgreen'),lwd=2) | |
| #Evaluate Predicted Density at Observed Value my.es | |
| abline(v=my.es,lty=2,lwd=2,col='red') | |
| valNull=nullPredF(my.es,N) | |
| valAlt=altPredF(my.es,N) | |
| points(pch=19,c(my.es,my.es),c(valNull,valAlt)) | |
| cat("Bayes factor (alt/null) is ",valAlt/valNull) | |
| BayesFactor::ttestBF(x = x, nullInterval = c(.99999999999999999999999,1.0000000000000000000000001)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment