Frequentist inference
Fixed-effects ANOVAGeneral linear models: mixing continuous and categorical covariatesOutput plot as PDFSimple linear regressionSkewness & Kurtosis
data(ToothGrowth)
## Example plot from ?ToothGrowth
coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
xlab = "ToothGrowth data: length vs dose, given type of supplement")
## Treat dose as a factor
ToothGrowth$dose = factor(ToothGrowth$dose)
levels(ToothGrowth$dose) = c("Low", "Medium", "High")
summary(aov(len ~ supp*dose, data=ToothGrowth))
#install.packages("xtable")
library(xtable)
xtable(x, caption = NULL, label = NULL, align = NULL, digits = NULL,
display = NULL, auto = FALSE, ...)
print(xtable(d), type="html")
print(xtable(d), type="latex") # anova table to latex
#https://cran.r-project.org/web/packages/xtable/index.html
#https://rmarkdown.rstudio.com/
data(ToothGrowth) # model log2 of dose instead of dose directly ToothGrowth$dose = log2(ToothGrowth$dose) # Classical analysis for comparison lmToothGrowth <- lm(len ~ supp + dose + supp:dose, data=ToothGrowth) summary(lmToothGrowth)
x<- (1:5) y<- (1:111) pdf(file=file.choose()) hist(x) plot(x, type='o') dev.off()
Notes
file.choose() is a very handy command which saves the work associated with defining absolute and relative paths which can be quite cumbersome.
list.files for non-interactive selection.
choose.files for selecting multiple files interactively.
More
https://www.rdocumentation.org/packages/base/versions/3.5.2/topics/file.choose
# http://wiki.math.yorku.ca/index.php/VR4:_Chapter_1_summary
x <- seq( 1, 20, 0.5) # sequence from 1 to 20 by 0.5
x # print it
w <- 1 + x/2 # a weight vector
y <- x + w*rnorm(x) # generate y with standard deviations
# equal to w
dum <- data.frame( x, y, w ) # make a data frame of 3 columns
dum # typing it's name shows the object
rm( x, y, w ) # remove the original variables
fm <- lm ( y ~ x, data = dum) # fit a simple linear regression (between x and y)
summary( fm ) # and look at the analysis
fm1 <- lm( y ~ x, data = dum, weight = 1/w^2 ) # a weighted regression
summary(fm1)
lrf <- loess( y ~ x, dum) # a smooth regression curve using a
# modern local regression function
attach( dum ) # make the columns in the data frame visible
# as variables (Note: before this command, typing x and
# y would yield errors)
plot( x, y) # a standard scatterplot to which we will
# fit regression curves
lines( spline( x, fitted (lrf)), col = 2) # add in the local regression using spline
# interpolation between calculated
# points
abline(0, 1, lty = 3, col = 3) # add in the true regression line (lty=3: read line type=dotted;
# col=3: read colour=green; type "?par" for details)
abline(fm, col = 4) # add in the unweighted regression line
abline(fm1, lty = 4, col = 5) # add in the weighted regression line
plot( fitted(fm), resid(fm),
xlab = "Fitted Values",
ylab = "Residuals") # a standard regression diagnostic plot
# to check for heteroscedasticity.
# The data were generated from a
# heteroscedastic process. Can you see
# it from this plot?
qqnorm( resid(fm) ) # Normal scores to check for skewness,
# kurtosis and outliers. How would you
# expect heteroscedasticity to show up?
search() # Have a look at the search path
detach() # Remove the data frame from the search path
search() # Have another look. How did it change?
rm(fm, fm1, lrf, dum) # Clean up (this is not that important in R
# for reasons we will understand later)
library(rgl)
rgl.open()
rgl.bg(col='white')
x <- sort(rnorm(1000))
y <- rnorm(1000)
z1 <- atan2(x,y)
z <- rnorm(1000) + atan2(x, y)
plot3d(x, y, z1, col = rainbow(1000), size = 2)
bbox3d()
skew(x, na.rm = TRUE)
kurtosi(x, na.rm = TRUE)
#x A data.frame or matrix
#na.rm how to treat missing data
## The function is currently defined as
function (x, na.rm = TRUE)
{
if (length(dim(x)) == 0) {
if (na.rm) {
x <- x[!is.na(x)]
}
mx <- mean(x)
sdx <- sd(x,na.rm=na.rm)
kurt <- sum((x - mx)^4)/(length(x) * sd(x)^4) -3
} else {
kurt <- rep(NA,dim(x)[2])
mx <- colMeans(x,na.rm=na.rm)
sdx <- sd(x,na.rm=na.rm)
for (i in 1:dim(x)[2]) {
kurt[i] <- sum((x[,i] - mx[i])^4, na.rm = na.rm)/((length(x[,i]) - sum(is.na(x[,i]))) * sdx[i]^4) -3
}
}
return(kurt)
}
Bayes Factor analysis
ttestBF function (one sample) BFmcmc functionttestBF function (two sample)Bayesian meta analysisBayes factor robustness analysis, one-sidedAnimating Robustness-Check of Bayes Factor PriorsSkewness & KurtosisanovaBF function
#https://richarddmorey.github.io/BayesFactor/ bf = ttestBF(x = diffScores) ## Equivalently: ## bf = ttestBF(x = sleep$extra[1:10],y=sleep$extra[11:20], paired=TRUE) bf
Comment
ttestBF function, which performs the “JZS” t test described by Rouder, Speckman, Sun, Morey, and Iverson (2009).
#https://richarddmorey.github.io/BayesFactor/ chains2 = recompute(chains, iterations = 10000) plot(chains2[,1:2])
Comment
The posterior function returns a object of type BFmcmc, which inherits the methods of the mcmc class from the coda package.
#https://richarddmorey.github.io/BayesFactor/ ## Compute Bayes factor bf = ttestBF(formula = weight ~ feed, data = chickwts) bf ### chains = posterior(bf, iterations = 10000) plot(chains[,2])
#https://richarddmorey.github.io/BayesFactor/ ## Bem's t statistics from four selected experiments t = c(-.15, 2.39, 2.42, 2.43) N = c(100, 150, 97, 99) ### ## Do analysis again, without nullInterval restriction bf = meta.ttestBF(t=t, n1=N, rscale=1) ## Obtain posterior samples chains = posterior(bf, iterations = 10000) plot(chains)
Comment
ttestBF function, which performs the “JZS” t test described by Rouder, Speckman, Sun, Morey, and Iverson (2009).
## This source code is licensed under the FreeBSD license
## (c) 2013 Felix Schönbrodt
#' @title Plots a comparison of a sequence of priors for t test Bayes factors
#'
#' @details
#'
#'
#' @param ts A vector of t values
#' @param ns A vector of corresponding sample sizes
#' @param rs The sequence of rs that should be tested. r should run up to 2 (higher values are implausible; E.-J. Wagenmakers, personal communication, Aug 22, 2013)
#' @param labels Names for the studies (displayed in the facet headings)
#' @param dots Values of r's which should be marked with a red dot
#' @param plot If TRUE, a ggplot is returned. If false, a data frame with the computed Bayes factors is returned
#' @param sides If set to "two" (default), a two-sided Bayes factor is computed. If set to "one", a one-sided Bayes factor is computed. In this case, it is assumed that positive t values correspond to results in the predicted direction and negative t values to results in the unpredicted direction. For details, see Wagenmakers, E. J., & Morey, R. D. (2013). Simple relation between one-sided and two-sided Bayesian point-null hypothesis tests.
#' @param nrow Number of rows of the faceted plot.
#' @param forH1 Defines the direction of the BF. If forH1 is TRUE, BF > 1 speak in favor of H1 (i.e., the quotient is defined as H1/H0). If forH1 is FALSE, it's the reverse direction.
#'
#' @references
#'
#' Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson, G. (2009). Bayesian t-tests for accepting and rejecting the null hypothesis. Psychonomic Bulletin and Review, 16, 225-237.
#' Wagenmakers, E.-J., & Morey, R. D. (2013). Simple relation between one-sided and two-sided Bayesian point-null hypothesis tests. Manuscript submitted for publication
#' Wagenmakers, E.-J., Wetzels, R., Borsboom, D., Kievit, R. & van der Maas, H. L. J. (2011). Yes, psychologists must change the way they analyze their data: Clarifications for Bem, Utts, & Johnson (2011)
BFrobustplot <- function(
ts, ns, rs=seq(0, 2, length.out=200), dots=1, plot=TRUE,
labels=c(), sides="two", nrow=2, xticks=3, forH1=TRUE)
{
library(BayesFactor)
# compute one-sided p-values from ts and ns
ps <- pt(ts, df=ns-1, lower.tail = FALSE) # one-sided test
# add the dots location to the sequences of r's
rs <- c(rs, dots)
res <- data.frame()
for (r in rs) {
# first: calculate two-sided BF
B_e0 <- c()
for (i in 1:length(ts))
B_e0 <- c(B_e0, exp(ttest.tstat(t = ts[i], n1 = ns[i], rscale=r)$bf))
# second: calculate one-sided BF
B_r0 <- c() for (i in 1:length(ts)) { if (ts[i] > 0) {
# correct direction
B_r0 <- c(B_r0, (2 - 2*ps[i])*B_e0[i])
} else {
# wrong direction
B_r0 <- c(B_r0, (1 - ps[i])*2*B_e0[i])
}
}
res0 <- data.frame(t=ts, n=ns, BF_two=B_e0, BF_one=B_r0, r=r) if (length(labels) > 0) {
res0$labels <- labels
res0$heading <- factor(1:length(labels), labels=paste0(labels, "\n(t = ", ts, ", df = ", ns-1, ")"), ordered=TRUE)
} else {
res0$heading <- factor(1:length(ts), labels=paste0("t = ", ts, ", df = ", ns-1), ordered=TRUE)
}
res <- rbind(res, res0)
}
# define the measure to be plotted: one- or two-sided?
res$BF <- res[, paste0("BF_", sides)]
# Flip BF if requested
if (forH1 == FALSE) {
res$BF <- 1/res$BF
}
if (plot==TRUE) {
library(ggplot2)
p1 <- ggplot(res, aes(x=r, y=log(BF))) + geom_line() + facet_wrap(~heading, nrow=nrow) + theme_bw() + ylab("log(BF)")
p1 <- p1 + geom_hline(yintercept=c(c(-log(c(30, 10, 3)), log(c(3, 10, 30)))), linetype="dotted", color="darkgrey")
p1 <- p1 + geom_hline(yintercept=log(1), linetype="dashed", color="darkgreen")
# add the dots
p1 <- p1 + geom_point(data=res[res$r %in% dots,], aes(x=r, y=log(BF)), color="red", size=2)
# add annotation
p1 <- p1 + annotate("text", x=max(rs)*1.8, y=-2.85, label=paste0("Strong~H[", ifelse(forH1==TRUE,0,1), "]"), hjust=1, vjust=.5, size=3, color="black", parse=TRUE)
p1 <- p1 + annotate("text", x=max(rs)*1.8, y=-1.7 , label=paste0("Moderate~H[", ifelse(forH1==TRUE,0,1), "]"), hjust=1, vjust=.5, size=3, color="black", parse=TRUE)
p1 <- p1 + annotate("text", x=max(rs)*1.8, y=-.55 , label=paste0("Anectodal~H[", ifelse(forH1==TRUE,0,1), "]"), hjust=1, vjust=.5, size=3, color="black", parse=TRUE)
p1 <- p1 + annotate("text", x=max(rs)*1.8, y=2.86 , label=paste0("Strong~H[", ifelse(forH1==TRUE,1,0), "]"), hjust=1, vjust=.5, size=3, color="black", parse=TRUE)
p1 <- p1 + annotate("text", x=max(rs)*1.8, y=1.7 , label=paste0("Moderate~H[", ifelse(forH1==TRUE,1,0), "]"), hjust=1, vjust=.5, size=3, color="black", parse=TRUE)
p1 <- p1 + annotate("text", x=max(rs)*1.8, y=.55 , label=paste0("Anectodal~H[", ifelse(forH1==TRUE,1,0), "]"), hjust=1, vjust=.5, vjust=.5, size=3, color="black", parse=TRUE)
# set scale ticks
p1 <- p1 + scale_y_continuous(breaks=c(c(-log(c(30, 10, 3)), 0, log(c(3, 10, 30)))), labels=c("-log(30)", "-log(10)", "-log(3)", "log(1)", "log(3)", "log(10)", "log(30)"))
p1 <- p1 + scale_x_continuous(breaks=seq(min(rs), max(rs), length.out=xticks))
return(p1)
} else {
return(res)
}
}
References
Rouder, J. N., Speckman, P. L., Sun, D., Morey, R. D., & Iverson, G. (2009). Bayesian t-tests for accepting and rejecting the null hypothesis. Psychonomic Bulletin and Review, 16, 225-237
#https://jimgrange.wordpress.com/2015/11/27/animating-robustness-check-of-bayes-factor-priors/
### plots of prior robustness check
# gif created from PNG plots using http://gifmaker.me/
# clear R's memory
rm(list = ls())
# load the Bayes factor package
library(BayesFactor)
# set working directory (will be used to save files here, so make sure
# this is where you want to save your plots!)
setwd <- "D:/Work//Blog_YouTube code//Blog/Prior Robust Visualisation/plots"
#------------------------------------------------------------------------------
### declare some variables for the analysis
# what is the t-value for the data?
tVal <- 3.098
# how many points in the prior should be explored?
nPoints <- 1000
# what Cauchy rates should be explored?
cauchyRates <- seq(from = 0.01, to = 1.5, length.out = nPoints)
# what effect sizes should be plotted?
effSize <- seq(from = -2, to = 2, length.out = nPoints)
# get the Bayes factor for each prior value
bayesFactors <- sapply(cauchyRates, function(x)
exp(ttest.tstat(t = tVal, n1 = 76, rscale = x)[['bf']]))
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
### do the plotting
# how many plots do we want to produce?
nPlots <- 50
plotWidth <- round(seq(from = 1, to = nPoints, length.out = nPlots), 0)
# loop over each plot
for(i in plotWidth){
# set up the file
currFile <- paste(getwd(), "/plot_", i, ".png", sep = "")
# initiate the png file
png(file = currFile, width = 1200, height = 1000, res = 200)
# change the plotting window so plots appear side-by-side
par(mfrow = c(1, 2))
#----
# do the prior density plot
d <- dcauchy(effSize, scale = cauchyRates[i])
plot(effSize, d, type = "l", ylim = c(0, 5), xlim = c(-2, 2),
ylab = "Density", xlab = "Effect Size (d)", lwd = 2, col = "gray48",
main = paste("Rate (r) = ", round(cauchyRates[i], 3), sep = ""))
#-----
# do the Bayes factor plot
plot(cauchyRates, bayesFactors, type = "l", lwd = 2, col = "gray48",
ylim = c(0, max(bayesFactors)), xaxt = "n",
xlab = "Cauchy Prior Width (r)", ylab = "Bayes Factor (10)")
abline(h = 0, lwd = 1)
abline(h = 6, col = "black", lty = 2, lwd = 2)
axis(1, at = seq(0, 1.5, 0.25))
# add the BF at the default Cauchy point
points(0.707, 9.97, col = "black", cex = 1.5, pch = 21, bg = "red")
# add the BF for the Cauchy prior currently being plotted
points(cauchyRates[i], bayesFactors[i], col = "black", pch = 21, cex = 1.3,
bg = "cyan")
# add legend
legend(x = 0.25, y = 3, legend = c("r = 0.707", paste("r = ",
round(cauchyRates[i], 3),
sep = "")),
pch = c(21, 21), lty = c(NA, NA), lwd = c(NA, NA), pt.cex = c(1, 1),
col = c("black", "black"), pt.bg = c("red", "cyan"), bty = "n")
# save the current plot
dev.off()
}
#------------------------------------------------------------------------------
skew(x, na.rm = TRUE)
kurtosi(x, na.rm = TRUE)
#x A data.frame or matrix
#na.rm how to treat missing data
## The function is currently defined as
function (x, na.rm = TRUE)
{
if (length(dim(x)) == 0) {
if (na.rm) {
x <- x[!is.na(x)]
}
mx <- mean(x)
sdx <- sd(x,na.rm=na.rm)
kurt <- sum((x - mx)^4)/(length(x) * sd(x)^4) -3
} else {
kurt <- rep(NA,dim(x)[2])
mx <- colMeans(x,na.rm=na.rm)
sdx <- sd(x,na.rm=na.rm)
for (i in 1:dim(x)[2]) {
kurt[i] <- sum((x[,i] - mx[i])^4, na.rm = na.rm)/((length(x[,i]) - sum(is.na(x[,i]))) * sdx[i]^4) -3
}
}
return(kurt)
}
#https://richarddmorey.github.io/BayesFactor/
data(ToothGrowth)
## Example plot from ?ToothGrowth
coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
xlab = "ToothGrowth data: length vs dose, given type of supplement")
bf = anovaBF(len ~ supp*dose, data=ToothGrowth)
bf
plot(bf[3:4] / bf[2])
#reduce number of comparisons (alpha inflation)
bf = anovaBF(len ~ supp*dose, data=ToothGrowth, whichModels="top")
bf
bfMainEffects = lmBF(len ~ supp + dose, data = ToothGrowth) bfInteraction = lmBF(len ~ supp + dose + supp:dose, data = ToothGrowth) ## Compare the two models bf = bfInteraction / bfMainEffects bf ## newbf = recompute(bf, iterations = 500000) newbf # posterior function ## Sample from the posterior of the full model chains = posterior(bfInteraction, iterations = 10000) ## 1:13 are the only "interesting" parameters summary(chains[,1:13]) #lot posterior of selected effects plot(chains[,4:6])
Comment
Reduce number of comparisons by specifyig the xact model of interest a priori.
Markov chain Monte Carlo methods
MCMC 1GBEST RBEST - Bayesian MCMC power analysis
# MODIFIED FROM BEST.R FOR ONE GROUP INSTEAD OF TWO.
# Version of Dec 02, 2015.
# John K. Kruschke
# johnkruschke@gmail.com
# http://www.indiana.edu/~kruschke/BEST/
#
# This program is believed to be free of errors, but it comes with no guarantee!
# The user bears all responsibility for interpreting the results.
# Please check the webpage above for updates or corrections.
#
### ***************************************************************
### ******** SEE FILE BEST1Gexample.R FOR INSTRUCTIONS ************
### ***************************************************************
source("openGraphSaveGraph.R") # graphics functions for Windows, Mac OS, Linux
BEST1Gmcmc = function( y , numSavedSteps=100000, thinSteps=1, showMCMC=FALSE) {
# This function generates an MCMC sample from the posterior distribution.
# Description of arguments:
# showMCMC is a flag for displaying diagnostic graphs of the chains.
# If F (the default), no chain graphs are displayed. If T, they are.
require(rjags)
#------------------------------------------------------------------------------
# THE MODEL.
modelString = "
model {
for ( i in 1:Ntotal ) {
y[i] ~ dt( mu , tau , nu )
}
mu ~ dnorm( muM , muP )
tau <- 1/pow( sigma , 2 )
sigma ~ dunif( sigmaLow , sigmaHigh )
nu ~ dexp(1/30)
}
" # close quote for modelString
# Write out modelString to a text file
writeLines( modelString , con="BESTmodel.txt" )
#------------------------------------------------------------------------------
# THE DATA.
# Load the data:
Ntotal = length(y)
# Specify the data in a list, for later shipment to JAGS:
dataList = list(
y = y ,
Ntotal = Ntotal ,
muM = mean(y) ,
muP = 0.000001 * 1/sd(y)^2 ,
sigmaLow = sd(y) / 1000 ,
sigmaHigh = sd(y) * 1000
)
#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Initial values of MCMC chains based on data:
mu = mean(y)
sigma = sd(y)
# Regarding initial values in next line: (1) sigma will tend to be too big if
# the data have outliers, and (2) nu starts at 5 as a moderate value. These
# initial values keep the burn-in period moderate.
initsList = list( mu = mu , sigma = sigma , nu = 5 )
#------------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "mu" , "sigma" , "nu" ) # The parameters to be monitored
adaptSteps = 500 # Number of steps to "tune" the samplers
burnInSteps = 1000
nChains = 3
nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains )
# Create, initialize, and adapt the model:
jagsModel = jags.model( "BESTmodel.txt" , data=dataList , inits=initsList ,
n.chains=nChains , n.adapt=adaptSteps )
# Burn-in:
cat( "Burning in the MCMC chain...\n" )
update( jagsModel , n.iter=burnInSteps )
# The saved MCMC chain:
cat( "Sampling final MCMC chain...\n" )
codaSamples = coda.samples( jagsModel , variable.names=parameters ,
n.iter=nIter , thin=thinSteps )
# resulting codaSamples object has these indices:
# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
#------------------------------------------------------------------------------
# EXAMINE THE RESULTS
if ( showMCMC ) {
openGraph(width=7,height=7)
autocorr.plot( codaSamples[[1]] , ask=FALSE )
show( gelman.diag( codaSamples ) )
effectiveChainLength = effectiveSize( codaSamples )
show( effectiveChainLength )
}
# Convert coda-object codaSamples to matrix object for easier handling.
# But note that this concatenates the different chains into one long chain.
# Result is mcmcChain[ stepIdx , paramIdx ]
mcmcChain = as.matrix( codaSamples )
return( mcmcChain )
} # end function BESTmcmc
#==============================================================================
BEST1Gplot = function( y , mcmcChain , compValm=0.0 , ROPEm=NULL ,
compValsd=NULL , ROPEsd=NULL ,
ROPEeff=NULL , showCurve=FALSE , pairsPlot=FALSE ) {
# This function plots the posterior distribution (and data).
# Description of arguments:
# y is the data vector.
# mcmcChain is a list of the type returned by function BEST1G.
# compValm is hypothesized mean for comparing with estimated mean.
# ROPEm is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the mean.
# ROPEsd is a two element vector, such as c(14,16), specifying the limit
# of the ROPE on the difference of standard deviations.
# ROPEeff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the effect size.
# showCurve is TRUE or FALSE and indicates whether the posterior should
# be displayed as a histogram (by default) or by an approximate curve.
# pairsPlot is TRUE or FALSE and indicates whether scatterplots of pairs
# of parameters should be displayed.
mu = mcmcChain[,"mu"]
sigma = mcmcChain[,"sigma"]
nu = mcmcChain[,"nu"]
if ( pairsPlot ) {
# Plot the parameters pairwise, to see correlations:
openGraph(width=7,height=7)
nPtToPlot = 1000
plotIdx = floor(seq(1,length(mu),by=length(mu)/nPtToPlot))
panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
usr = par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r = (cor(x, y))
txt = format(c(r, 0.123456789), digits=digits)[1]
txt = paste(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex=1.25 ) # was cex=cex.cor*r
}
pairs( cbind( mu , sigma , log10(nu) )[plotIdx,] ,
labels=c( expression(mu) ,
expression(sigma) ,
expression(log10(nu)) ) ,
lower.panel=panel.cor , col="skyblue" )
}
# Define matrix for storing summary info:
summaryInfo = NULL
source("plotPost.R")
# Set up window and layout:
openGraph(width=6.0,height=5.0)
layout( matrix( c(3,3,4,4,5,5, 1,1,1,1,2,2) , nrow=6, ncol=2 , byrow=FALSE ) )
par( mar=c(3.5,3.5,2.5,0.5) , mgp=c(2.25,0.7,0) )
# Select thinned steps in chain for plotting of posterior predictive curves:
chainLength = NROW( mcmcChain )
nCurvesToPlot = 30
stepIdxVec = seq( 1 , chainLength , floor(chainLength/nCurvesToPlot) )
xRange = range( y )
xLim = c( xRange[1]-0.1*(xRange[2]-xRange[1]) ,
xRange[2]+0.1*(xRange[2]-xRange[1]) )
xVec = seq( xLim[1] , xLim[2] , length=200 )
maxY = max( dt( 0 , df=max(nu[stepIdxVec]) ) / min(sigma[stepIdxVec]) )
# Plot data y and smattering of posterior predictive curves:
stepIdx = 1
plot( xVec , dt( (xVec-mu[stepIdxVec[stepIdx]])/sigma[stepIdxVec[stepIdx]] ,
df=nu[stepIdxVec[stepIdx]] )/sigma[stepIdxVec[stepIdx]] ,
ylim=c(0,maxY) , cex.lab=1.75 ,
type="l" , col="skyblue" , lwd=1 , xlab="y" , ylab="p(y)" ,
main="Data w. Post. Pred." )
for ( stepIdx in 2:length(stepIdxVec) ) {
lines(xVec, dt( (xVec-mu[stepIdxVec[stepIdx]])/sigma[stepIdxVec[stepIdx]] ,
df=nu[stepIdxVec[stepIdx]] )/sigma[stepIdxVec[stepIdx]] ,
type="l" , col="skyblue" , lwd=1 )
}
histBinWd = median(sigma)/2
histCenter = mean(mu)
histBreaks = sort( c( seq( histCenter-histBinWd/2 , min(xVec)-histBinWd/2 ,
-histBinWd ),
seq( histCenter+histBinWd/2 , max(xVec)+histBinWd/2 ,
histBinWd ) , xLim ) )
histInfo = hist( y , plot=FALSE , breaks=histBreaks )
yPlotVec = histInfo$density
yPlotVec[ yPlotVec==0.0 ] = NA
xPlotVec = histInfo$mids
xPlotVec[ yPlotVec==0.0 ] = NA
points( xPlotVec , yPlotVec , type="h" , lwd=3 , col="red" )
text( max(xVec) , maxY , bquote(N==.(length(y))) , adj=c(1.1,1.1) )
# Plot posterior distribution of parameter nu:
paramInfo = plotPost( log10(nu) , col="skyblue" , # breaks=30 ,
showCurve=showCurve ,
xlab=bquote("log10("*nu*")") , cex.lab = 1.75 , showMode=TRUE ,
main="Normality" ) # (<0.7 suggests kurtosis)
summaryInfo = rbind( summaryInfo , paramInfo )
rownames(summaryInfo)[NROW(summaryInfo)] = "log10(nu)"
# Plot posterior distribution of parameters mu1, mu2, and their difference:
xlim = range( c( mu , compValm ) )
paramInfo = plotPost( mu , xlim=xlim , cex.lab = 1.75 , ROPE=ROPEm ,
showCurve=showCurve , compVal = compValm ,
xlab=bquote(mu) , main=paste("Mean") ,
col="skyblue" )
summaryInfo = rbind( summaryInfo , paramInfo )
rownames(summaryInfo)[NROW(summaryInfo)] = "mu"
# Plot posterior distribution of sigma:
xlim=range( sigma )
paramInfo = plotPost( sigma , xlim=xlim , cex.lab = 1.75 , ROPE=ROPEsd ,
showCurve=showCurve , compVal=compValsd ,
xlab=bquote(sigma) , main=paste("Std. Dev.") ,
col="skyblue" , showMode=TRUE )
summaryInfo = rbind( summaryInfo , paramInfo )
rownames(summaryInfo)[NROW(summaryInfo)] = "sigma"
# Plot of estimated effect size:
effectSize = ( mu - compValm ) / sigma
paramInfo = plotPost( effectSize , compVal=0 , ROPE=ROPEeff ,
showCurve=showCurve ,
xlab=bquote( (mu-.(compValm)) / sigma ),
showMode=TRUE , cex.lab=1.75 , main="Effect Size" , col="skyblue" )
summaryInfo = rbind( summaryInfo , paramInfo )
rownames(summaryInfo)[NROW(summaryInfo)] = "effSz"
return( summaryInfo )
} # end of function BEST1Gplot
#==============================================================================
# MODIFIED 2015 DEC 02.
# John K. Kruschke
# johnkruschke@gmail.com
# http://www.indiana.edu/~kruschke/BEST/
#
# This program is believed to be free of errors, but it comes with no guarantee!
# The user bears all responsibility for interpreting the results.
# Please check the webpage above for updates or corrections.
#
### ***************************************************************
### ******** SEE FILE BESTexample.R FOR INSTRUCTIONS **************
### ***************************************************************
# Load various essential functions:
source("DBDA2E-utilities.R")
BESTmcmc = function( y1 , y2 ,
priorOnly=FALSE , showMCMC=FALSE ,
numSavedSteps=20000 , thinSteps=1 ,
mu1PriorMean = mean(c(y1,y2)) ,
mu1PriorSD = sd(c(y1,y2))*5 ,
mu2PriorMean = mean(c(y1,y2)) ,
mu2PriorSD = sd(c(y1,y2))*5 ,
sigma1PriorMode = sd(c(y1,y2)) ,
sigma1PriorSD = sd(c(y1,y2))*5 ,
sigma2PriorMode = sd(c(y1,y2)) ,
sigma2PriorSD = sd(c(y1,y2))*5 ,
nuPriorMean = 30 ,
nuPriorSD = 30 ,
runjagsMethod=runjagsMethodDefault ,
nChains=nChainsDefault ) {
# This function generates an MCMC sample from the posterior distribution.
# Description of arguments:
# showMCMC is a flag for displaying diagnostic graphs of the chains.
# If F (the default), no chain graphs are displayed. If T, they are.
#------------------------------------------------------------------------------
# THE DATA.
# Load the data:
y = c( y1 , y2 ) # combine data into one vector
x = c( rep(1,length(y1)) , rep(2,length(y2)) ) # create group membership code
Ntotal = length(y)
# Specify the data and prior constants in a list, for later shipment to JAGS:
if ( priorOnly ) {
dataList = list(
# y = y ,
x = x ,
Ntotal = Ntotal ,
mu1PriorMean = mu1PriorMean ,
mu1PriorSD = mu1PriorSD ,
mu2PriorMean = mu2PriorMean ,
mu2PriorSD = mu2PriorSD ,
Sh1 = gammaShRaFromModeSD( mode=sigma1PriorMode ,
sd=sigma1PriorSD )$shape ,
Ra1 = gammaShRaFromModeSD( mode=sigma1PriorMode ,
sd=sigma1PriorSD )$rate ,
Sh2 = gammaShRaFromModeSD( mode=sigma2PriorMode ,
sd=sigma2PriorSD )$shape ,
Ra2 = gammaShRaFromModeSD( mode=sigma2PriorMode ,
sd=sigma2PriorSD )$rate ,
ShNu = gammaShRaFromMeanSD( mean=nuPriorMean , sd=nuPriorSD )$shape ,
RaNu = gammaShRaFromMeanSD( mean=nuPriorMean , sd=nuPriorSD )$rate
)
} else {
dataList = list(
y = y ,
x = x ,
Ntotal = Ntotal ,
mu1PriorMean = mu1PriorMean ,
mu1PriorSD = mu1PriorSD ,
mu2PriorMean = mu2PriorMean ,
mu2PriorSD = mu2PriorSD ,
Sh1 = gammaShRaFromModeSD( mode=sigma1PriorMode ,
sd=sigma1PriorSD )$shape ,
Ra1 = gammaShRaFromModeSD( mode=sigma1PriorMode ,
sd=sigma1PriorSD )$rate ,
Sh2 = gammaShRaFromModeSD( mode=sigma2PriorMode ,
sd=sigma2PriorSD )$shape ,
Ra2 = gammaShRaFromModeSD( mode=sigma2PriorMode ,
sd=sigma2PriorSD )$rate ,
ShNu = gammaShRaFromMeanSD( mean=nuPriorMean , sd=nuPriorSD )$shape ,
RaNu = gammaShRaFromMeanSD( mean=nuPriorMean , sd=nuPriorSD )$rate
)
}
#----------------------------------------------------------------------------
# THE MODEL.
modelString = "
model {
for ( i in 1:Ntotal ) {
y[i] ~ dt( mu[x[i]] , 1/sigma[x[i]]^2 , nu )
}
mu[1] ~ dnorm( mu1PriorMean , 1/mu1PriorSD^2 ) # prior for mu[1]
sigma[1] ~ dgamma( Sh1 , Ra1 ) # prior for sigma[1]
mu[2] ~ dnorm( mu2PriorMean , 1/mu2PriorSD^2 ) # prior for mu[2]
sigma[2] ~ dgamma( Sh2 , Ra2 ) # prior for sigma[2]
nu ~ dgamma( ShNu , RaNu ) # prior for nu
}
" # close quote for modelString
# Write out modelString to a text file
writeLines( modelString , con="BESTmodel.txt" )
#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Initial values of MCMC chains based on data:
mu = c( mean(y1) , mean(y2) )
sigma = c( sd(y1) , sd(y2) )
# Regarding initial values in next line: (1) sigma will tend to be too big if
# the data have outliers, and (2) nu starts at 5 as a moderate value. These
# initial values keep the burn-in period moderate.
initsList = list( mu = mu , sigma = sigma , nu = 5 )
#------------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "mu" , "sigma" , "nu" ) # The parameters to be monitored
adaptSteps = 500 # Number of steps to "tune" the samplers
burnInSteps = 1000
runJagsOut <- run.jags( method=runjagsMethod , model="BESTmodel.txt" , monitor=parameters , data=dataList , inits=initsList , n.chains=nChains , adapt=adaptSteps , burnin=burnInSteps , sample=ceiling(numSavedSteps/nChains) , thin=thinSteps , summarise=FALSE , plots=FALSE ) codaSamples = as.mcmc.list( runJagsOut ) # resulting codaSamples object has these indices: # codaSamples[[ chainIdx ]][ stepIdx , paramIdx ] ## Original version, using rjags: # nChains = 3 # nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # # Create, initialize, and adapt the model: # jagsModel = jags.model( "BESTmodel.txt" , data=dataList , inits=initsList , # n.chains=nChains , n.adapt=adaptSteps ) # # Burn-in: # cat( "Burning in the MCMC chain...\n" ) # update( jagsModel , n.iter=burnInSteps ) # # The saved MCMC chain: # cat( "Sampling final MCMC chain...\n" ) # codaSamples = coda.samples( jagsModel , variable.names=parameters , # n.iter=nIter , thin=thinSteps ) #------------------------------------------------------------------------------ # EXAMINE THE RESULTS if ( showMCMC ) { for ( paramName in varnames(codaSamples) ) { diagMCMC( codaSamples , parName=paramName ) } } # Convert coda-object codaSamples to matrix object for easier handling. # But note that this concatenates the different chains into one long chain. # Result is mcmcChain[ stepIdx , paramIdx ] mcmcChain = as.matrix( codaSamples ) return( mcmcChain ) } # end function BESTmcmc #============================================================================== BESTsummary = function( y1 , y2 , mcmcChain ) { #source("HDIofMCMC.R") # in DBDA2E-utilities.R mcmcSummary = function( paramSampleVec , compVal=NULL ) { meanParam = mean( paramSampleVec ) medianParam = median( paramSampleVec ) dres = density( paramSampleVec ) modeParam = dres$x[which.max(dres$y)] hdiLim = HDIofMCMC( paramSampleVec ) if ( !is.null(compVal) ) { pcgtCompVal = ( 100 * sum( paramSampleVec > compVal )
/ length( paramSampleVec ) )
} else {
pcgtCompVal=NA
}
return( c( meanParam , medianParam , modeParam , hdiLim , pcgtCompVal ) )
}
# Define matrix for storing summary info:
summaryInfo = matrix( 0 , nrow=9 , ncol=6 , dimnames=list(
PARAMETER=c( "mu1" , "mu2" , "muDiff" , "sigma1" , "sigma2" , "sigmaDiff" ,
"nu" , "nuLog10" , "effSz" ),
SUMMARY.INFO=c( "mean" , "median" , "mode" , "HDIlow" , "HDIhigh" ,
"pcgtZero" )
) )
summaryInfo[ "mu1" , ] = mcmcSummary( mcmcChain[,"mu[1]"] )
summaryInfo[ "mu2" , ] = mcmcSummary( mcmcChain[,"mu[2]"] )
summaryInfo[ "muDiff" , ] = mcmcSummary( mcmcChain[,"mu[1]"]
- mcmcChain[,"mu[2]"] ,
compVal=0 )
summaryInfo[ "sigma1" , ] = mcmcSummary( mcmcChain[,"sigma[1]"] )
summaryInfo[ "sigma2" , ] = mcmcSummary( mcmcChain[,"sigma[2]"] )
summaryInfo[ "sigmaDiff" , ] = mcmcSummary( mcmcChain[,"sigma[1]"]
- mcmcChain[,"sigma[2]"] ,
compVal=0 )
summaryInfo[ "nu" , ] = mcmcSummary( mcmcChain[,"nu"] )
summaryInfo[ "nuLog10" , ] = mcmcSummary( log10(mcmcChain[,"nu"]) )
N1 = length(y1)
N2 = length(y2)
effSzChain = ( ( mcmcChain[,"mu[1]"] - mcmcChain[,"mu[2]"] )
/ sqrt( ( mcmcChain[,"sigma[1]"]^2
+ mcmcChain[,"sigma[2]"]^2 ) / 2 ) )
summaryInfo[ "effSz" , ] = mcmcSummary( effSzChain , compVal=0 )
# Or, use sample-size weighted version:
# effSz = ( mu1 - mu2 ) / sqrt( ( sigma1^2 *(N1-1) + sigma2^2 *(N2-1) )
# / (N1+N2-2) )
# Be sure also to change plot label in BESTplot function, below.
return( summaryInfo )
}
#==============================================================================
BESTplot = function( y1 , y2 , mcmcChain , ROPEm=NULL , ROPEsd=NULL ,
ROPEeff=NULL , showCurve=FALSE , pairsPlot=FALSE ) {
# This function plots the posterior distribution (and data).
# Description of arguments:
# y1 and y2 are the data vectors.
# mcmcChain is a list of the type returned by function BTT.
# ROPEm is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of means.
# ROPEsd is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of standard deviations.
# ROPEeff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the effect size.
# showCurve is TRUE or FALSE and indicates whether the posterior should
# be displayed as a histogram (by default) or by an approximate curve.
# pairsPlot is TRUE or FALSE and indicates whether scatterplots of pairs
# of parameters should be displayed.
mu1 = mcmcChain[,"mu[1]"]
mu2 = mcmcChain[,"mu[2]"]
sigma1 = mcmcChain[,"sigma[1]"]
sigma2 = mcmcChain[,"sigma[2]"]
nu = mcmcChain[,"nu"]
if ( pairsPlot ) {
# Plot the parameters pairwise, to see correlations:
openGraph(width=7,height=7)
nPtToPlot = 1000
plotIdx = floor(seq(1,length(mu1),by=length(mu1)/nPtToPlot))
panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
usr = par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r = (cor(x, y))
txt = format(c(r, 0.123456789), digits=digits)[1]
txt = paste(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex=1.25 ) # was cex=cex.cor*r
}
pairs( cbind( mu1 , mu2 , sigma1 , sigma2 , log10(nu) )[plotIdx,] ,
labels=c( expression(mu[1]) , expression(mu[2]) ,
expression(sigma[1]) , expression(sigma[2]) ,
expression(log10(nu)) ) ,
lower.panel=panel.cor , col="skyblue" )
}
# source("plotPost.R") # in DBDA2E-utilities.R
# Set up window and layout:
openGraph(width=6.0,height=8.0)
layout( matrix( c(4,5,7,8,3,1,2,6,9,10) , nrow=5, byrow=FALSE ) )
par( mar=c(3.5,3.5,2.5,0.5) , mgp=c(2.25,0.7,0) )
# Select thinned steps in chain for plotting of posterior predictive curves:
chainLength = NROW( mcmcChain )
nCurvesToPlot = 30
stepIdxVec = seq( 1 , chainLength , floor(chainLength/nCurvesToPlot) )
xRange = range( c(y1,y2) )
if ( isTRUE( all.equal( xRange[2] , xRange[1] ) ) ) {
meanSigma = mean( c(sigma1,sigma2) )
xRange = xRange + c( -meanSigma , meanSigma )
}
xLim = c( xRange[1]-0.1*(xRange[2]-xRange[1]) ,
xRange[2]+0.1*(xRange[2]-xRange[1]) )
xVec = seq( xLim[1] , xLim[2] , length=200 )
maxY = max( dt( 0 , df=max(nu[stepIdxVec]) ) /
min(c(sigma1[stepIdxVec],sigma2[stepIdxVec])) )
# Plot data y1 and smattering of posterior predictive curves:
stepIdx = 1
plot( xVec , dt( (xVec-mu1[stepIdxVec[stepIdx]])/sigma1[stepIdxVec[stepIdx]] ,
df=nu[stepIdxVec[stepIdx]] )/sigma1[stepIdxVec[stepIdx]] ,
ylim=c(0,maxY) , cex.lab=1.75 ,
type="l" , col="skyblue" , lwd=1 , xlab="y" , ylab="p(y)" ,
main="Data Group 1 w. Post. Pred." )
for ( stepIdx in 2:length(stepIdxVec) ) {
lines(xVec, dt( (xVec-mu1[stepIdxVec[stepIdx]])/sigma1[stepIdxVec[stepIdx]] ,
df=nu[stepIdxVec[stepIdx]] )/sigma1[stepIdxVec[stepIdx]] ,
type="l" , col="skyblue" , lwd=1 )
}
histBinWd = median(sigma1)/2
histCenter = mean(mu1)
histBreaks = sort( c( seq( histCenter-histBinWd/2 , min(xVec)-histBinWd/2 ,
-histBinWd ),
seq( histCenter+histBinWd/2 , max(xVec)+histBinWd/2 ,
histBinWd ) , xLim ) )
histInfo = hist( y1 , plot=FALSE , breaks=histBreaks )
yPlotVec = histInfo$density
yPlotVec[ yPlotVec==0.0 ] = NA
xPlotVec = histInfo$mids
xPlotVec[ yPlotVec==0.0 ] = NA
points( xPlotVec , yPlotVec , type="h" , lwd=3 , col="red" )
text( max(xVec) , maxY , bquote(N[1]==.(length(y1))) , adj=c(1.1,1.1) )
# Plot data y2 and smattering of posterior predictive curves:
stepIdx = 1
plot( xVec , dt( (xVec-mu2[stepIdxVec[stepIdx]])/sigma2[stepIdxVec[stepIdx]] ,
df=nu[stepIdxVec[stepIdx]] )/sigma2[stepIdxVec[stepIdx]] ,
ylim=c(0,maxY) , cex.lab=1.75 ,
type="l" , col="skyblue" , lwd=1 , xlab="y" , ylab="p(y)" ,
main="Data Group 2 w. Post. Pred." )
for ( stepIdx in 2:length(stepIdxVec) ) {
lines(xVec, dt( (xVec-mu2[stepIdxVec[stepIdx]])/sigma2[stepIdxVec[stepIdx]] ,
df=nu[stepIdxVec[stepIdx]] )/sigma2[stepIdxVec[stepIdx]] ,
type="l" , col="skyblue" , lwd=1 )
}
histBinWd = median(sigma2)/2
histCenter = mean(mu2)
histBreaks = sort( c( seq( histCenter-histBinWd/2 , min(xVec)-histBinWd/2 ,
-histBinWd ),
seq( histCenter+histBinWd/2 , max(xVec)+histBinWd/2 ,
histBinWd ) , xLim ) )
histInfo = hist( y2 , plot=FALSE , breaks=histBreaks )
yPlotVec = histInfo$density
yPlotVec[ yPlotVec==0.0 ] = NA
xPlotVec = histInfo$mids
xPlotVec[ yPlotVec==0.0 ] = NA
points( xPlotVec , yPlotVec , type="h" , lwd=3 , col="red" )
text( max(xVec) , maxY , bquote(N[2]==.(length(y2))) , adj=c(1.1,1.1) )
# Plot posterior distribution of parameter nu:
histInfo = plotPost( log10(nu) , col="skyblue" , # breaks=30 ,
showCurve=showCurve ,
xlab=bquote("log10("*nu*")") , cex.lab = 1.75 ,
cenTend=c("mode","median","mean")[1] ,
main="Normality" ) # (<0.7 suggests kurtosis)
# Plot posterior distribution of parameters mu1, mu2, and their difference:
xlim = range( c( mu1 , mu2 ) )
histInfo = plotPost( mu1 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(mu[1]) , main=paste("Group",1,"Mean") ,
col="skyblue" )
histInfo = plotPost( mu2 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(mu[2]) , main=paste("Group",2,"Mean") ,
col="skyblue" )
histInfo = plotPost( mu1-mu2 , compVal=0 , showCurve=showCurve ,
xlab=bquote(mu[1] - mu[2]) , cex.lab = 1.75 , ROPE=ROPEm ,
main="Difference of Means" , col="skyblue" )
# Plot posterior distribution of param's sigma1, sigma2, and their difference:
xlim=range( c( sigma1 , sigma2 ) )
histInfo = plotPost( sigma1 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(sigma[1]) , main=paste("Group",1,"Std. Dev.") ,
col="skyblue" , cenTend=c("mode","median","mean")[1] )
histInfo = plotPost( sigma2 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(sigma[2]) , main=paste("Group",2,"Std. Dev.") ,
col="skyblue" , cenTend=c("mode","median","mean")[1] )
histInfo = plotPost( sigma1-sigma2 ,
compVal=0 , showCurve=showCurve ,
xlab=bquote(sigma[1] - sigma[2]) , cex.lab = 1.75 ,
ROPE=ROPEsd ,
main="Difference of Std. Dev.s" , col="skyblue" ,
cenTend=c("mode","median","mean")[1] )
# Plot of estimated effect size. Effect size is d-sub-a from
# Macmillan & Creelman, 1991; Simpson & Fitter, 1973; Swets, 1986a, 1986b.
effectSize = ( mu1 - mu2 ) / sqrt( ( sigma1^2 + sigma2^2 ) / 2 )
histInfo = plotPost( effectSize , compVal=0 , ROPE=ROPEeff ,
showCurve=showCurve ,
xlab=bquote( (mu[1]-mu[2])
/sqrt((sigma[1]^2 +sigma[2]^2 )/2 ) ),
cenTend=c("mode","median","mean")[1] , cex.lab=1.0 , main="Effect Size" ,
col="skyblue" )
# Or use sample-size weighted version:
# Hedges 1981; Wetzels, Raaijmakers, Jakab & Wagenmakers 2009.
# N1 = length(y1)
# N2 = length(y2)
# effectSize = ( mu1 - mu2 ) / sqrt( ( sigma1^2 *(N1-1) + sigma2^2 *(N2-1) )
# / (N1+N2-2) )
# Be sure also to change BESTsummary function, above.
# histInfo = plotPost( effectSize , compVal=0 , ROPE=ROPEeff ,
# showCurve=showCurve ,
# xlab=bquote( (mu[1]-mu[2])
# /sqrt((sigma[1]^2 *(N[1]-1)+sigma[2]^2 *(N[2]-1))/(N[1]+N[2]-2)) ),
# cenTend=c("mode","median","mean")[1] , cex.lab=1.0 , main="Effect Size" , col="skyblue" )
return( BESTsummary( y1 , y2 , mcmcChain ) )
} # end of function BESTplot
#==============================================================================
BESTpower = function( mcmcChain , N1 , N2 , ROPEm , ROPEsd , ROPEeff ,
maxHDIWm , maxHDIWsd , maxHDIWeff , nRep=200 ,
mcmcLength=10000 , saveName="BESTpower.Rdata" ,
showFirstNrep=0 ) {
# This function estimates power.
# Description of arguments:
# mcmcChain is a matrix of the type returned by function BESTmcmc.
# N1 and N2 are sample sizes for the two groups.
# ROPEm is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of means.
# ROPEsd is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of standard deviations.
# ROPEeff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the effect size.
# maxHDIWm is the maximum desired width of the 95% HDI on the difference
# of means.
# maxHDIWsd is the maximum desired width of the 95% HDI on the difference
# of standard deviations.
# maxHDIWeff is the maximum desired width of the 95% HDI on the effect size.
# nRep is the number of simulated experiments used to estimate the power.
#source("HDIofMCMC.R") # in DBDA2E-utilities.R
#source("HDIofICDF.R") # in DBDA2E-utilities.R
chainLength = NROW( mcmcChain )
# Select thinned steps in chain for posterior predictions:
stepIdxVec = seq( 1 , chainLength , floor(chainLength/nRep) )
goalTally = list( HDIm.gt.ROPE = c(0) ,
HDIm.lt.ROPE = c(0) ,
HDIm.in.ROPE = c(0) ,
HDIm.wd.max = c(0) ,
HDIsd.gt.ROPE = c(0) ,
HDIsd.lt.ROPE = c(0) ,
HDIsd.in.ROPE = c(0) ,
HDIsd.wd.max = c(0) ,
HDIeff.gt.ROPE = c(0) ,
HDIeff.lt.ROPE = c(0) ,
HDIeff.in.ROPE = c(0) ,
HDIeff.wd.max = c(0) )
power = list( HDIm.gt.ROPE = c(0,0,0) ,
HDIm.lt.ROPE = c(0,0,0) ,
HDIm.in.ROPE = c(0,0,0) ,
HDIm.wd.max = c(0,0,0) ,
HDIsd.gt.ROPE = c(0,0,0) ,
HDIsd.lt.ROPE = c(0,0,0) ,
HDIsd.in.ROPE = c(0,0,0) ,
HDIsd.wd.max = c(0,0,0) ,
HDIeff.gt.ROPE = c(0,0,0) ,
HDIeff.lt.ROPE = c(0,0,0) ,
HDIeff.in.ROPE = c(0,0,0) ,
HDIeff.wd.max = c(0,0,0) )
nSim = 0
for ( stepIdx in stepIdxVec ) {
nSim = nSim + 1
cat( "\n:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\n" )
cat( paste( "Power computation: Simulated Experiment" , nSim , "of" ,
length(stepIdxVec) , ":\n\n" ) )
# Get parameter values for this simulation:
mu1Val = mcmcChain[stepIdx,"mu[1]"]
mu2Val = mcmcChain[stepIdx,"mu[2]"]
sigma1Val = mcmcChain[stepIdx,"sigma[1]"]
sigma2Val = mcmcChain[stepIdx,"sigma[2]"]
nuVal = mcmcChain[stepIdx,"nu"]
# Generate simulated data:
y1 = rt( N1 , df=nuVal ) * sigma1Val + mu1Val
y2 = rt( N2 , df=nuVal ) * sigma2Val + mu2Val
# Get posterior for simulated data:
simChain = BESTmcmc( y1, y2, numSavedSteps=mcmcLength, thinSteps=1,
showMCMC=FALSE )
if (nSim <= showFirstNrep ) { BESTplot( y1, y2, simChain, ROPEm=ROPEm, ROPEsd=ROPEsd, ROPEeff=ROPEeff ) } # Assess which goals were achieved and tally them: HDIm = HDIofMCMC( simChain[,"mu[1]"] - simChain[,"mu[2]"] ) if ( HDIm[1] > ROPEm[2] ) {
goalTally$HDIm.gt.ROPE = goalTally$HDIm.gt.ROPE + 1 }
if ( HDIm[2] < ROPEm[1] ) { goalTally$HDIm.lt.ROPE = goalTally$HDIm.lt.ROPE + 1 } if ( HDIm[1] > ROPEm[1] & HDIm[2] < ROPEm[2] ) {
goalTally$HDIm.in.ROPE = goalTally$HDIm.in.ROPE + 1 }
if ( HDIm[2] - HDIm[1] < maxHDIWm ) { goalTally$HDIm.wd.max = goalTally$HDIm.wd.max + 1 } HDIsd = HDIofMCMC( simChain[,"sigma[1]"] - simChain[,"sigma[2]"] ) if ( HDIsd[1] > ROPEsd[2] ) {
goalTally$HDIsd.gt.ROPE = goalTally$HDIsd.gt.ROPE + 1 }
if ( HDIsd[2] < ROPEsd[1] ) { goalTally$HDIsd.lt.ROPE = goalTally$HDIsd.lt.ROPE + 1 } if ( HDIsd[1] > ROPEsd[1] & HDIsd[2] < ROPEsd[2] ) {
goalTally$HDIsd.in.ROPE = goalTally$HDIsd.in.ROPE + 1 }
if ( HDIsd[2] - HDIsd[1] < maxHDIWsd ) { goalTally$HDIsd.wd.max = goalTally$HDIsd.wd.max + 1 } HDIeff = HDIofMCMC( ( simChain[,"mu[1]"] - simChain[,"mu[2]"] ) / sqrt( ( simChain[,"sigma[1]"]^2 + simChain[,"sigma[2]"]^2 ) / 2 ) ) if ( HDIeff[1] > ROPEeff[2] ) {
goalTally$HDIeff.gt.ROPE = goalTally$HDIeff.gt.ROPE + 1 }
if ( HDIeff[2] < ROPEeff[1] ) { goalTally$HDIeff.lt.ROPE = goalTally$HDIeff.lt.ROPE + 1 } if ( HDIeff[1] > ROPEeff[1] & HDIeff[2] < ROPEeff[2] ) {
goalTally$HDIeff.in.ROPE = goalTally$HDIeff.in.ROPE + 1 }
if ( HDIeff[2] - HDIeff[1] < maxHDIWeff ) { goalTally$HDIeff.wd.max = goalTally$HDIeff.wd.max + 1 } for ( i in 1:length(goalTally) ) { s1 = 1 + goalTally[[i]] s2 = 1 + ( nSim - goalTally[[i]] ) power[[i]][1] = s1/(s1+s2) power[[i]][2:3] = HDIofICDF( qbeta , shape1=s1 , shape2=s2 ) } cat( "\nAfter ", nSim , " Simulated Experiments, Power Is: \n( mean, HDI_low, HDI_high )\n" ) for ( i in 1:length(power) ) { cat( names(power[i]) , round(power[[i]],3) , "\n" ) } save( nSim , power , file=saveName ) } return( power ) } # end of function BESTpower #============================================================================== makeData = function( mu1 , sd1 , mu2 , sd2 , nPerGrp , pcntOut=0 , sdOutMult=2.0 , rnd.seed=NULL ) { # Auxilliary function for generating random values from a # mixture of normal distibutions. if(!is.null(rnd.seed)){set.seed(rnd.seed)} # Set seed for random values. nOut = ceiling(nPerGrp*pcntOut/100) # Number of outliers. nIn = nPerGrp - nOut # Number from main distribution. if ( pcntOut > 100 | pcntOut < 0 ) { stop("pcntOut must be between 0 and 100.") } if ( pcntOut > 0 & pcntOut < 1 ) { warning("pcntOut is specified as percentage 0-100, not proportion 0-1.") } if ( pcntOut > 50 ) {
warning("pcntOut indicates more than 50% outliers; did you intend this?")
}
if ( nOut < 2 & pcntOut > 0 ) {
stop("Combination of nPerGrp and pcntOut yields too few outliers.")
}
if ( nIn < 2 ) { stop("Too few non-outliers.") } sdN = function( x ) { sqrt( mean( (x-mean(x))^2 ) ) } # genDat = function( mu , sigma , Nin , Nout , sigmaOutMult=sdOutMult ) { y = rnorm(n=Nin) # Random values in main distribution. y = ((y-mean(y))/sdN(y))*sigma + mu # Translate to exactly realize desired. yOut = NULL if ( Nout > 0 ) {
yOut = rnorm(n=Nout) # Random values for outliers
yOut=((yOut-mean(yOut))/sdN(yOut))*(sigma*sdOutMult)+mu # Realize exactly.
}
y = c(y,yOut) # Concatenate main with outliers.
return(y)
}
y1 = genDat( mu=mu1 , sigma=sd1 , Nin=nIn , Nout=nOut )
y2 = genDat( mu=mu2 , sigma=sd2 , Nin=nIn , Nout=nOut )
#
# Set up window and layout:
openGraph(width=7,height=7) # Plot the data.
layout(matrix(1:2,nrow=2))
histInfo = hist( y1 , main="Simulated Data" , col="pink2" , border="white" ,
xlim=range(c(y1,y2)) , breaks=30 , prob=TRUE )
text( max(c(y1,y2)) , max(histInfo$density) ,
bquote(N==.(nPerGrp)) , adj=c(1,1) )
xlow=min(histInfo$breaks)
xhi=max(histInfo$breaks)
xcomb=seq(xlow,xhi,length=1001)
lines( xcomb , dnorm(xcomb,mean=mu1,sd=sd1)*nIn/(nIn+nOut) +
dnorm(xcomb,mean=mu1,sd=sd1*sdOutMult)*nOut/(nIn+nOut) , lwd=3 )
lines( xcomb , dnorm(xcomb,mean=mu1,sd=sd1)*nIn/(nIn+nOut) ,
lty="dashed" , col="green", lwd=3)
lines( xcomb , dnorm(xcomb,mean=mu1,sd=sd1*sdOutMult)*nOut/(nIn+nOut) ,
lty="dashed" , col="red", lwd=3)
histInfo = hist( y2 , main="" , col="pink2" , border="white" ,
xlim=range(c(y1,y2)) , breaks=30 , prob=TRUE )
text( max(c(y1,y2)) , max(histInfo$density) ,
bquote(N==.(nPerGrp)) , adj=c(1,1) )
xlow=min(histInfo$breaks)
xhi=max(histInfo$breaks)
xcomb=seq(xlow,xhi,length=1001)
lines( xcomb , dnorm(xcomb,mean=mu2,sd=sd2)*nIn/(nIn+nOut) +
dnorm(xcomb,mean=mu2,sd=sd2*sdOutMult)*nOut/(nIn+nOut) , lwd=3)
lines( xcomb , dnorm(xcomb,mean=mu2,sd=sd2)*nIn/(nIn+nOut) ,
lty="dashed" , col="green", lwd=3)
lines( xcomb , dnorm(xcomb,mean=mu2,sd=sd2*sdOutMult)*nOut/(nIn+nOut) ,
lty="dashed" , col="red", lwd=3)
#
return( list( y1=y1 , y2=y2 ) )
}
#==============================================================================
# Version of May 26, 2012. Re-checked on 2015 May 08.
# John K. Kruschke
# johnkruschke@gmail.com
# http://www.indiana.edu/~kruschke/BEST/
#
# This program is believed to be free of errors, but it comes with no guarantee!
# The user bears all responsibility for interpreting the results.
# Please check the webpage above for updates or corrections.
#
### ***************************************************************
### ******** SEE FILE BESTexample.R FOR INSTRUCTIONS **************
### ***************************************************************
# OPTIONAL: Clear R's memory and graphics:
rm(list=ls()) # Careful! This clears all of R's memory!
graphics.off() # This closes all of R's graphics windows.
# Get the functions loaded into R's working memory:
source("BEST.R")
#-------------------------------------------------------------------------------
# RETROSPECTIVE POWER ANALYSIS.
# !! This section assumes you have already run BESTexample.R !!
# Re-load the saved data and MCMC chain from the previously conducted
# Bayesian analysis. This re-loads the variables y1, y2, mcmcChain, etc.
load( "BESTexampleMCMC.Rdata" )
power = BESTpower( mcmcChain , N1=length(y1) , N2=length(y2) ,
ROPEm=c(-0.1,0.1) , ROPEsd=c(-0.1,0.1) , ROPEeff=c(-0.1,0.1) ,
maxHDIWm=2.0 , maxHDIWsd=2.0 , maxHDIWeff=0.2 , nRep=1000 ,
mcmcLength=10000 , saveName = "BESTexampleRetroPower.Rdata" )
#-------------------------------------------------------------------------------
# PROSPECTIVE POWER ANALYSIS, using fictitious strong data.
# Generate large fictitious data set that expresses hypothesis:
prospectData = makeData( mu1=108, sd1=17, mu2=100, sd2=15, nPerGrp=1000,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL )
y1pro = prospectData$y1 # Merely renames simulated data for convenience below.
y2pro = prospectData$y2 # Merely renames simulated data for convenience below.
# Generate Bayesian posterior distribution from fictitious data:
# (uses fewer than usual MCMC steps because it only needs nRep credible
# parameter combinations, not a high-resolution representation)
mcmcChainPro = BESTmcmc( y1pro , y2pro , numSavedSteps=2000 )
postInfoPro = BESTplot( y1pro , y2pro , mcmcChainPro , pairsPlot=TRUE )
save( y1pro, y2pro, mcmcChainPro, postInfoPro,
file="BESTexampleProPowerMCMC.Rdata" )
# Now compute the prospective power for planned sample sizes:
N1plan = N2plan = 50 # specify planned sample size
powerPro = BESTpower( mcmcChainPro , N1=N1plan , N2=N2plan , showFirstNrep=5 ,
ROPEm=c(-1.5,1.5) , ROPEsd=c(-0.0,0.0) , ROPEeff=c(-0.0,0.0) ,
maxHDIWm=15.0 , maxHDIWsd=10.0 , maxHDIWeff=1.0 , nRep=1000 ,
mcmcLength=10000 , saveName = "BESTexampleProPower.Rdata" )
#-------------------------------------------------------------------------------
Various plots
BeanplotsDensity and Rug plotDirac function - Interdisciplinarity versus OverspecialisationGoogle TrendsJASP prior & posterior plot
#https://cran.r-project.org/web/packages/beanplot/beanplot.pdf
par( mfrow = c( 1, 4 ) )
with(dataexp2, beanplot(v00, ylim = c(0,10), col="lightgray", main = "v00", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
with(dataexp2, beanplot(v01, ylim = c(0,10), col="lightgray", main = "v01", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
with(dataexp2, beanplot(v10, ylim = c(0,10), col="darkgray", main = "v10", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
with(dataexp2, beanplot(v11, ylim = c(0,10), col="darkgray", main = "v11", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
#with ylim and titles
with(dataexp2, beanplot(v00, ylab = "VAS Brightness rating", xlab = "beanplot showing distribution and mean",
ylim = c(0,10), col="darkgray", main = "v00", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf,
overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
install.packages ("beanplot")
library("beanplot")
par( mfrow = c( 1, 2 ) )
with(dataexp2, beanplot(v00, col="darkgray", main = "v00", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
with(dataexp2, beanplot(v01, col="darkgray", main = "v01", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
par( mfrow = c( 1, 2 ) )
with(dataexp2, beanplot(v10, col="darkgray", main = "v10", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
with(dataexp2, beanplot(v11, col="darkgray", main = "v11", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
stripchart(variable ~ factor, vertical=TRUE, method="stack",
ylab="variable", data=StackedData)
#two sided
with(dataexp1, beanplot(v00, v10, col="darkgray", main = "v00 vs. v10", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "both", jitter = NULL, beanlinewd = 2))
#verticalpar( mfrow = c( 2, 1 ) )
install.packages ("beanplot")
library("beanplot")
par( mfrow = c( 2, 1 ) ) # layout matrix
with(dataexp1, beanplot(v00, v10, main = "v00 vs. v10", kernel = "gaussian", cut = 3, col = list("black", c("grey", "white")), xlab="VAS",
cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = T, side = "both", jitter = NULL, beanlinewd = 2))
with(dataexp1, beanplot(v01, v11, main = "v01 vs. v11", kernel = "gaussian", cut = 3, col = list("black", c("grey", "white")), xlab="VAS",
cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = T, side = "both", jitter = NULL, beanlinewd = 2))
#exp1
dataexp1 <-
read.table("http://www.irrational-decisions.com/phd-thesis/dataexp1.csv",
header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)
par( mfrow = c( 1, 2) )
with(dataexp1, beanplot(v00, ylim = c(0,7), col="lightgray", main = "v00", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
with(dataexp1, beanplot(v10, ylim = c(0,7), col="lightgray", main = "v01", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
par( mfrow = c( 1, 2) )
with(dataexp1, beanplot(v01, ylim = c(3,10), col="darkgray", main = "v10", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
with(dataexp1, beanplot(v11, ylim = c(3,10), col="darkgray", main = "v11", kernel = "gaussian", cut = 3, cutmin = -Inf, cutmax = Inf, overallline = "mean", horizontal = FALSE, side = "no", jitter = NULL, beanlinewd = 2))
#histogram and rug plot plot(density(dataexp2$v00)) rug(dataexp2$v00, col=2, lwd=3.5) ############## library(MASS) plot(density(dataexp2$v00, bw=5)) rug(dataexp2$v00+ rnorm(length(dataexp2$v00), sd=5), col=2, lwd=3.5)
x = c(5) y = c(9) plot(x,y,type="h",xlim=c(0,10),ylim=c(0,10),lwd=2,col="blue",ylab="knowledge", xlab="diversity")
install.packages("gtrendsR")
library(gtrendsR)
gtrends(c("bayesian inference", "p-value"), time = "all") # since 2004
res <- gtrends(c("Bayesian inference", "Markov Chain Monte Carlo", "ANOVA"), geo = c("US", "GB", "DE", "CA",))
res <- gtrends("Markov Chain Monte Carlo", geo = c("US", "GB", "DE"))
plot(res)
#get vars
length(res)
str(res)
res$interest_over_time
ff<-res$interest_over_time
scatterplot(ff$hits~ff$date| geo, reg.line=FALSE, smooth=FALSE, spread=FALSE, boxplots=FALSE, span=0.5, ellipse=FALSE, levels=c(.5, .9), by.groups=TRUE, data=ff)
res2 <- gtrends(c("p-value", "bayes factor"))
library(plotrix)
.likelihoodShiftedT <- function(par, data) {
-sum(log(dt((data - par[1])/par[2], par[3])/par[2]))
}
.dposteriorShiftedT <- function(x, parameters, oneSided) { # function that returns the posterior density if (oneSided == FALSE) { dt((x - parameters[1])/parameters[2], parameters[3])/parameters[2] } else if (oneSided == "right") { ifelse(x >= 0, (dt((x - parameters[1])/parameters[2], parameters[3])/parameters[2])/pt((0 -
parameters[1])/parameters[2], parameters[3], lower.tail = FALSE),
0)
} else if (oneSided == "left") {
ifelse(x <= 0, (dt((x - parameters[1])/parameters[2], parameters[3])/parameters[2])/pt((0 -
parameters[1])/parameters[2], parameters[3], lower.tail = TRUE),
0)
}
}
.dprior <- function(x, r, oneSided = oneSided) {
# function that returns the prior density
if (oneSided == "right") {
y <- ifelse(x < 0, 0, 2/(pi * r * (1 + (x/r)^2)))
return(y)
}
if (oneSided == "left") {
y <- ifelse(x > 0, 0, 2/(pi * r * (1 + (x/r)^2)))
return(y)
} else {
return(1/(pi * r * (1 + (x/r)^2)))
}
}
.plotPosterior.ttest <- function(x = NULL, y = NULL, paired = FALSE, oneSided = FALSE,
BF, BFH1H0 = TRUE, iterations = 10000, rscale = "medium", lwd = 2,
cexPoints = 1.5, cexAxis = 1.2, cexYlab = 1.5, cexXlab = 1.5, cexTextBF = 1.4,
cexCI = 1.1, cexLegend = 1.2, lwdAxis = 1.2, addInformation = TRUE,
dontPlotData = FALSE) {
if (addInformation) {
par(mar = c(5.6, 5, 7, 4) + 0.1, las = 1)
} else {
par(mar = c(5.6, 5, 4, 4) + 0.1, las = 1)
}
if (dontPlotData) {
plot(1, type = "n", xlim = 0:1, ylim = 0:1, bty = "n", axes = FALSE,
xlab = "", ylab = "")
axis(1, at = 0:1, labels = FALSE, cex.axis = cexAxis, lwd = lwdAxis,
xlab = "")
axis(2, at = 0:1, labels = FALSE, cex.axis = cexAxis, lwd = lwdAxis,
ylab = "")
mtext(text = "Density", side = 2, las = 0, cex = cexYlab, line = 3.25)
mtext(expression(paste("Effect size", ~delta)), side = 1, cex = cexXlab,
line = 2.5)
return()
}
if (rscale == "medium") {
r <- sqrt(2)/2
}
if (rscale == "wide") {
r <- 1
}
if (rscale == "ultrawide") {
r <- sqrt(2)
}
if (mode(rscale) == "numeric") {
r <- rscale
}
if (oneSided == FALSE) {
nullInterval <- NULL
}
if (oneSided == "right") {
nullInterval <- c(0, Inf)
}
if (oneSided == "left") {
nullInterval <- c(-Inf, 0)
}
# sample from delta posterior
samples <- BayesFactor::ttestBF(x = x, y = y, paired = paired, posterior = TRUE,
iterations = iterations, rscale = r)
delta <- samples[, "delta"]
# fit shifted t distribution
if (is.null(y)) {
deltaHat <- mean(x)/sd(x)
N <- length(x)
df <- N - 1
sigmaStart <- 1/N
} else if (paired) {
deltaHat <- mean(x - y)/sd(x - y)
N <- length(x)
df <- N - 1
sigmaStart <- 1/N
} else if (!is.null(y) && !paired) {
N1 <- length(x)
N2 <- length(y)
sdPooled <- sqrt(((N1 - 1) * var(x) + (N2 - 1) * var(y))/(N1 +
N2))
deltaHat <- (mean(x) - mean(y))/sdPooled
df <- N1 + N2 - 2
sigmaStart <- sqrt(N1 * N2/(N1 + N2))
}
if (sigmaStart < 0.01)
sigmaStart <- 0.01
parameters <- try(silent = TRUE, expr = optim(par = c(deltaHat, sigmaStart,
df), fn = .likelihoodShiftedT, data = delta, method = "BFGS")$par)
if (class(parameters) == "try-error") {
parameters <- try(silent = TRUE, expr = optim(par = c(deltaHat,
sigmaStart, df), fn = .likelihoodShiftedT, data = delta, method = "Nelder-Mead")$par)
}
if (BFH1H0) {
BF10 <- BF
BF01 <- 1/BF10
} else {
BF01 <- BF
BF10 <- 1/BF01
}
# set limits plot
xlim <- vector("numeric", 2)
if (oneSided == FALSE) {
xlim[1] <- min(-2, quantile(delta, probs = 0.01)[[1]])
xlim[2] <- max(2, quantile(delta, probs = 0.99)[[1]])
if (length(x) < 10) {
if (addInformation) {
stretch <- 1.52
} else {
stretch <- 1.4
}
} else {
stretch <- 1.2 } } if (oneSided == "right") { # if (length(delta[delta >= 0]) < 10) return('Plotting is not
# possible: To few posterior samples in tested interval')
xlim[1] <- min(-2, quantile(delta[delta >= 0], probs = 0.01)[[1]])
xlim[2] <- max(2, quantile(delta[delta >= 0], probs = 0.99)[[1]])
if (any(is.na(xlim))) {
xlim[1] <- min(-2, .qShiftedT(0.01, parameters, oneSided = "right"))
xlim[2] <- max(2, .qShiftedT(0.99, parameters, oneSided = "right"))
}
stretch <- 1.32
}
if (oneSided == "left") {
# if (length(delta[delta <= 0]) < 10) return('Plotting is not
# possible: To few posterior samples in tested interval')
xlim[1] <- min(-2, quantile(delta[delta <= 0], probs = 0.01)[[1]])
xlim[2] <- max(2, quantile(delta[delta <= 0], probs = 0.99)[[1]])
if (any(is.na(xlim))) {
xlim[1] <- min(-2, .qShiftedT(0.01, parameters, oneSided = "left"))
xlim[2] <- max(2, .qShiftedT(0.99, parameters, oneSided = "left"))
}
stretch <- 1.32
}
xticks <- pretty(xlim)
ylim <- vector("numeric", 2)
ylim[1] <- 0
dmax <- optimize(function(x) .dposteriorShiftedT(x, parameters = parameters,
oneSided = oneSided), interval = range(xticks), maximum = TRUE)$objective
ylim[2] <- max(stretch * .dprior(0, r, oneSided = oneSided), stretch *
dmax) # get maximum density
# calculate position of 'nice' tick marks and create labels
yticks <- pretty(ylim)
xlabels <- formatC(xticks, 1, format = "f")
ylabels <- formatC(yticks, 1, format = "f")
# compute 95% credible interval & median:
if (oneSided == FALSE) {
CIlow <- quantile(delta, probs = 0.025)[[1]]
CIhigh <- quantile(delta, probs = 0.975)[[1]]
medianPosterior <- median(delta)
if (any(is.na(c(CIlow, CIhigh, medianPosterior)))) {
CIlow <- .qShiftedT(0.025, parameters, oneSided = FALSE)
CIhigh <- .qShiftedT(0.975, parameters, oneSided = FALSE)
medianPosterior <- .qShiftedT(0.5, parameters, oneSided = FALSE)
}
}
if (oneSided == "right") {
CIlow <- quantile(delta[delta >= 0], probs = 0.025)[[1]]
CIhigh <- quantile(delta[delta >= 0], probs = 0.975)[[1]]
medianPosterior <- median(delta[delta >= 0])
if (any(is.na(c(CIlow, CIhigh, medianPosterior)))) {
CIlow <- .qShiftedT(0.025, parameters, oneSided = "right")
CIhigh <- .qShiftedT(0.975, parameters, oneSided = "right")
medianPosterior <- .qShiftedT(0.5, parameters, oneSided = "right")
}
}
if (oneSided == "left") {
CIlow <- quantile(delta[delta <= 0], probs = 0.025)[[1]]
CIhigh <- quantile(delta[delta <= 0], probs = 0.975)[[1]]
medianPosterior <- median(delta[delta <= 0])
if (any(is.na(c(CIlow, CIhigh, medianPosterior)))) {
CIlow <- .qShiftedT(0.025, parameters, oneSided = "left")
CIhigh <- .qShiftedT(0.975, parameters, oneSided = "left")
medianPosterior <- .qShiftedT(0.5, parameters, oneSided = "left")
}
}
posteriorLine <- .dposteriorShiftedT(x = seq(min(xticks), max(xticks),
length.out = 1000), parameters = parameters, oneSided = oneSided)
xlim <- c(min(CIlow, range(xticks)[1]), max(range(xticks)[2], CIhigh)) plot(1, 1, xlim = xlim, ylim = range(yticks), ylab = "", xlab = "", type = "n", axes = FALSE) lines(seq(min(xticks), max(xticks), length.out = 1000), posteriorLine, lwd = lwd) lines(seq(min(xticks), max(xticks), length.out = 1000), .dprior(seq(min(xticks), max(xticks), length.out = 1000), r = r, oneSided = oneSided), lwd = lwd, lty = 3) axis(1, at = xticks, labels = xlabels, cex.axis = cexAxis, lwd = lwdAxis) axis(2, at = yticks, labels = ylabels, , cex.axis = cexAxis, lwd = lwdAxis) if (nchar(ylabels[length(ylabels)]) > 4) {
mtext(text = "Density", side = 2, las = 0, cex = cexYlab, line = 4)
} else if (nchar(ylabels[length(ylabels)]) == 4) {
mtext(text = "Density", side = 2, las = 0, cex = cexYlab, line = 3.25)
} else if (nchar(ylabels[length(ylabels)]) < 4) {
mtext(text = "Density", side = 2, las = 0, cex = cexYlab, line = 2.85)
}
mtext(expression(paste("Effect size", ~delta)), side = 1, cex = cexXlab,
line = 2.5)
points(0, .dprior(0, r, oneSided = oneSided), col = "black", pch = 21,
bg = "grey", cex = cexPoints)
if (oneSided == FALSE) {
heightPosteriorAtZero <- .dposteriorShiftedT(0, parameters = parameters,
oneSided = oneSided)
} else if (oneSided == "right") {
posteriorLineLargerZero <- posteriorLine[posteriorLine > 0]
heightPosteriorAtZero <- posteriorLineLargerZero[1]
} else if (oneSided == "left") {
posteriorLineLargerZero <- posteriorLine[posteriorLine > 0]
heightPosteriorAtZero <- posteriorLineLargerZero[length(posteriorLineLargerZero)]
}
points(0, heightPosteriorAtZero, col = "black", pch = 21, bg = "grey",
cex = cexPoints)
### 95% credible interval
# enable plotting in margin
par(xpd = TRUE)
yCI <- grconvertY(dmax, "user", "ndc") + 0.04
yCI <- grconvertY(yCI, "ndc", "user")
arrows(CIlow, yCI, CIhigh, yCI, angle = 90, code = 3, length = 0.1,
lwd = lwd)
medianText <- formatC(medianPosterior, digits = 3, format = "f")
if (addInformation) {
# display BF10 value
offsetTopPart <- 0.06
yy <- grconvertY(0.75 + offsetTopPart, "ndc", "user")
yy2 <- grconvertY(0.806 + offsetTopPart, "ndc", "user")
xx <- min(xticks) if (BF10 >= 1000000 | BF01 >= 1000000) {
BF10t <- formatC(BF10, 3, format = "e")
BF01t <- formatC(BF01, 3, format = "e")
}
if (BF10 < 1000000 & BF01 < 1000000) {
BF10t <- formatC(BF10, 3, format = "f")
BF01t <- formatC(BF01, 3, format = "f")
}
if (oneSided == FALSE) {
text(xx, yy2, bquote(BF[10] == .(BF10t)), cex = cexTextBF,
pos = 4)
text(xx, yy, bquote(BF[0][1] == .(BF01t)), cex = cexTextBF,
pos = 4)
}
if (oneSided == "right") {
text(xx, yy2, bquote(BF["+"][0] == .(BF10t)), cex = cexTextBF,
pos = 4)
text(xx, yy, bquote(BF[0]["+"] == .(BF01t)), cex = cexTextBF,
pos = 4)
}
if (oneSided == "left") {
text(xx, yy2, bquote(BF["-"][0] == .(BF10t)), cex = cexTextBF,
pos = 4)
text(xx, yy, bquote(BF[0]["-"] == .(BF01t)), cex = cexTextBF,
pos = 4)
}
yy <- grconvertY(0.756 + offsetTopPart, "ndc", "user")
yy2 <- grconvertY(0.812 + offsetTopPart, "ndc", "user")
CIText <- paste("95% CI: [", bquote(.(formatC(CIlow, 3, format = "f"))),
", ", bquote(.(formatC(CIhigh, 3, format = "f"))), "]", sep = "")
medianLegendText <- paste("median =", medianText)
text(max(xticks), yy2, medianLegendText, cex = 1.1, pos = 2)
text(max(xticks), yy, CIText, cex = 1.1, pos = 2)
# probability wheel
if (max(nchar(BF10t), nchar(BF01t)) <= 4) {
xx <- grconvertX(0.44, "ndc", "user")
}
if (max(nchar(BF10t), nchar(BF01t)) == 5) {
xx <- grconvertX(0.44 + 0.001 * 5, "ndc", "user")
}
if (max(nchar(BF10t), nchar(BF01t)) == 6) {
xx <- grconvertX(0.44 + 0.001 * 6, "ndc", "user")
}
if (max(nchar(BF10t), nchar(BF01t)) == 7) {
xx <- grconvertX(0.44 + 0.002 * max(nchar(BF10t), nchar(BF01t)),
"ndc", "user")
}
if (max(nchar(BF10t), nchar(BF01t)) == 8) {
xx <- grconvertX(0.44 + 0.003 * max(nchar(BF10t), nchar(BF01t)), "ndc", "user") } if (max(nchar(BF10t), nchar(BF01t)) > 8) {
xx <- grconvertX(0.44 + 0.005 * max(nchar(BF10t), nchar(BF01t)),
"ndc", "user")
}
yy <- grconvertY(0.788 + offsetTopPart, "ndc", "user")
# make sure that colored area is centered
radius <- 0.06 * diff(range(xticks))
A <- radius^2 * pi
alpha <- 2/(BF01 + 1) * A/radius^2
startpos <- pi/2 - alpha/2
# draw probability wheel
plotrix::floating.pie(xx, yy, c(BF10, 1), radius = radius, col = c("darkred",
"white"), lwd = 2, startpos = startpos)
yy <- grconvertY(0.865 + offsetTopPart, "ndc", "user")
yy2 <- grconvertY(0.708 + offsetTopPart, "ndc", "user")
if (oneSided == FALSE) {
text(xx, yy, "data|H1", cex = cexCI)
text(xx, yy2, "data|H0", cex = cexCI)
}
if (oneSided == "right") {
text(xx, yy, "data|H+", cex = cexCI)
text(xx, yy2, "data|H0", cex = cexCI)
}
if (oneSided == "left") {
text(xx, yy, "data|H-", cex = cexCI)
text(xx, yy2, "data|H0", cex = cexCI)
}
# add legend
CIText <- paste("95% CI: [", bquote(.(formatC(CIlow, 3, format = "f"))),
" ; ", bquote(.(formatC(CIhigh, 3, format = "f"))), "]", sep = "")
medianLegendText <- paste("median =", medianText)
}
mostPosterior <- mean(delta > mean(range(xticks)))
if (mostPosterior >= 0.5) {
legendPosition <- min(xticks)
legend(legendPosition, max(yticks), legend = c("Posterior", "Prior"),
lty = c(1, 3), bty = "n", lwd = c(lwd, lwd), cex = cexLegend,
xjust = 0, yjust = 1, x.intersp = 0.6, seg.len = 1.2)
} else {
legendPosition <- max(xticks)
legend(legendPosition, max(yticks), legend = c("Posterior", "Prior"),
lty = c(1, 3), bty = "n", lwd = c(lwd, lwd), cex = cexLegend,
xjust = 1, yjust = 1, x.intersp = 0.6, seg.len = 1.2)
}
}
### generate data ###
set.seed(1)
x <- rnorm(30, 0.15)
### calculate Bayes factor ###
library(BayesFactor)
BF <- extractBF(ttestBF(x, rscale = "medium"), onlybf = TRUE)
### plot ###
.plotPosterior.ttest(x = x, rscale = "medium", BF = BF)
Simulations
#monte carlo t test simulation
#install.packages("MonteCarlo")
library(MonteCarlo)
# Define function that generates data and applies the method of interest
ttest<-function(n,loc,scale){
# generate sample:
sample<-rnorm(n, loc, scale)
# calculate test statistic:
stat<-sqrt(n)*mean(sample)/sd(sample)
# get test decision:
decision<-abs(stat)>1.96
# return result:
return(list("decision"=decision))
}
# define parameter grid:
n_grid<-c(50,100,250,500)
loc_grid<-seq(0,1,0.2)
scale_grid<-c(1,2)
# collect parameter grids in list:
param_list=list("n"=n_grid, "loc"=loc_grid, "scale"=scale_grid)
# run simulation:
MC_result<-MonteCarlo(func=ttest, nrep=1000, param_list=param_list)
summary(MC_result)
R to WordPress
#install.packages("RWordPress", repos="http://www.omegahat.org/R", build=TRUE)
library(RWordPress)
options(WordPressLogin=c(user="password"),
WordPressURL="http://your_wp_installation.org/xmlrpc.php")
#[code lang='r']
#...
#[/code]
knit_hooks$set(output=function(x, options) paste("\\[code\\]\n", x, "\\[/code\\]\n", sep=""))
knit_hooks$set(source=function(x, options) paste("\\[code lang='r'\\]\n", x, "\\[/code\\]\n", sep=""))
knit2wp <- function(file) {
require(XML)
content <- readLines(file)
content <- htmlTreeParse(content, trim=FALSE)
## WP will add the h1 header later based on the title, so delete here
content$children$html$children$body$children$h1 <- NULL
content <- paste(capture.output(print(content$children$html$children$body,
indent=FALSE, tagSeparator="")),
collapse="\n")
content <- gsub("<?.body>", "", content) # remove body tag
## enclose code snippets in SyntaxHighlighter format
content <- gsub("<?pre><code class=\"r\">", "\\[code lang='r'\\]\\\n",
content)
content <- gsub("<?pre><code class=\"no-highlight\">", "\\[code\\]\\\n",
content)
content <- gsub("<?pre><code>", "\\[code\\]\\\n", content)
content <- gsub("<?/code></pre>", "\\[/code\\]\\\n", content)
return(content)
}
newPost(content=list(description=knit2wp('rerWorkflow.html'),
title='Workflow: Post R markdown to WordPress',
categories=c('R')),
publish=FALSE)
postID <- 99 # post id returned by newPost()
editPost(postID,
content=list(description=knit2wp('rerWorkflow.html'),
title='Workflow: Post R markdown to WordPress',
categories=c('R')),
publish=FALSE)
Shiny Apps
Exploring the diagnosticity of the p-valueP-value distribution and power curves for an independent two-tailed t-testBIC approximation for ANOVA designsWhen does a significant p-value indicate a true effect?
RGL
Scatterplot3d
Rcmdr
Video lectures
Packages/Libraries
Ask/answer questions on StackOverflow
Trending on GitHub
External URL: https://github.com/trending/r